Kstars

fitssepdetector.cpp
1 /*
2  SPDX-FileCopyrightText: 2004 Jasem Mutlaq
3  SPDX-FileCopyrightText: 2020 Eric Dejouhanet <[email protected]>
4 
5  SPDX-License-Identifier: GPL-2.0-or-later
6 
7  Some code fragments were adapted from Peter Kirchgessner's FITS plugin
8  SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
9 */
10 
11 #include "config-kstars.h"
12 
13 #include "fits_debug.h"
14 #include "fitssepdetector.h"
15 #include "fitsdata.h"
16 #include "Options.h"
17 #include "kspaths.h"
18 
19 #include <memory>
20 #include <math.h>
21 #include <QPointer>
22 #include <QtConcurrent>
23 
24 #ifdef HAVE_STELLARSOLVER
25 #include "ekos/auxiliary/stellarsolverprofileeditor.h"
26 #include <stellarsolver.h>
27 #else
28 #include <cstring>
29 #include "sep/sep.h"
30 #endif
31 
32 //void FITSSEPDetector::configure(const QString &param, const QVariant &value)
33 //{
34 // if (param == "numStars")
35 // numStars = value.toInt();
36 // else if (param == "fractionRemoved")
37 // fractionRemoved = value.toDouble();
38 // else if (param == "deblendNThresh")
39 // deblendNThresh = value.toInt();
40 // else if (param == "deblendMincont")
41 // deblendMincont = value.toDouble();
42 // else if (param == "radiusIsBoundary")
43 // radiusIsBoundary = value.toBool();
44 // else
45 // qCDebug(KSTARS_FITS) << "Bad SEP Parameter!!!!! " << param;
46 //}
47 
48 // TODO: (hy 4/11/2020)
49 // The api into these star detection methods should be generalized so that various parameters
50 // could be controlled by the caller. For instance, saturation removal, removal of N% of the largest stars
51 // number of stars desired, some parameters to the SEP algorithms, as all of these may have different
52 // optimal values for different uses. Also, there may be other desired outputs
53 // (e.g. unfiltered number of stars detected,background sky level). Waiting on rlancaste's
54 // investigations into SEP before doing this.
55 
56 QFuture<bool> FITSSEPDetector::findSources(QRect const &boundary)
57 {
58  return QtConcurrent::run(this, &FITSSEPDetector::findSourcesAndBackground, boundary);
59 }
60 
61 bool FITSSEPDetector::findSourcesAndBackground(QRect const &boundary)
62 {
63 #ifndef HAVE_STELLARSOLVER
64  Q_UNUSED(boundary)
65  return false;
66 #else
67  QList<Edge*> starCenters;
68  SkyBackground skyBG;
69  int maxStarsCount = getValue("maxStarsCount", 100000).toInt();
70 
71  int optionsProfileIndex = getValue("optionsProfileIndex", -1).toInt();
72  Ekos::ProfileGroup group = static_cast<Ekos::ProfileGroup>(getValue("optionsProfileGroup", 1).toInt());
73  QScopedPointer<StellarSolver, QScopedPointerDeleteLater> solver(new StellarSolver(m_ImageData->getStatistics(),
74  m_ImageData->getImageBuffer()));
75  QString filename = "";
76  QPointer<FITSData> image(m_ImageData);
77  switch(group)
78  {
79  case Ekos::AlignProfiles:
80  //So it should not be here if it is Align.
81  break;
82  case Ekos::GuideProfiles:
83  filename = "SavedGuideProfiles.ini";
84  break;
85  case Ekos::FocusProfiles:
86  filename = "SavedFocusProfiles.ini";
87  break;
88  case Ekos::HFRProfiles:
89  filename = "SavedHFRProfiles.ini";
90  break;
91  }
92 
93  QString savedOptionsProfiles = QDir(KSPaths::writableLocation(QStandardPaths::AppLocalDataLocation)).filePath(filename);
94  QList<SSolver::Parameters> optionsList;
95  if(QFile(savedOptionsProfiles).exists())
96  optionsList = StellarSolver::loadSavedOptionsProfiles(savedOptionsProfiles);
97  else
98  {
99  switch(group)
100  {
101  case Ekos::AlignProfiles:
102  optionsList = Ekos::getDefaultAlignOptionsProfiles();
103  break;
104  case Ekos::GuideProfiles:
105  optionsList = Ekos::getDefaultGuideOptionsProfiles();
106  break;
107  case Ekos::FocusProfiles:
108  optionsList = Ekos::getDefaultFocusOptionsProfiles();
109  break;
110  case Ekos::HFRProfiles:
111  optionsList = Ekos::getDefaultHFROptionsProfiles();
112  break;
113  }
114  }
115  if (optionsProfileIndex >= 0 && optionsList.count() > optionsProfileIndex)
116  {
117  auto params = optionsList[optionsProfileIndex];
118  params.partition = Options::stellarSolverPartition();
119  solver->setParameters(params);
120  qCDebug(KSTARS_FITS) << "Sextract with: " << optionsList[optionsProfileIndex].listName;
121  }
122  else
123  {
124  auto params = SSolver::Parameters(); // This is default
125  params.partition = Options::stellarSolverPartition();
126  solver->setParameters(params);
127  }
128 
130  const bool runHFR = group != Ekos::AlignProfiles;
131 
132  solver->setLogLevel(SSolver::LOG_NONE);
133  solver->setSSLogLevel(SSolver::LOG_OFF);
134 
135  if (boundary.isValid())
136  solver->extract(runHFR, boundary);
137  else
138  solver->extract(runHFR);
139 
140  stars = solver->getStarList();
141 
142  // If m_ImageData goes out of scope, also return.
143  if (stars.empty() || image.isNull())
144  return false;
145 
146  auto bg = solver->getBackground();
147 
148  skyBG.mean = bg.global;
149  skyBG.sigma = bg.globalrms;
150  skyBG.numPixelsInSkyEstimate = bg.bw * bg.bh;
151  skyBG.setStarsDetected(bg.num_stars_detected);
152  m_ImageData->setSkyBackground(skyBG);
153 
154  //There is more information that can be obtained by the Stellarsolver->
155  //Background info, Star positions(if a plate solve was done before), etc
156  //The information is available as long as the StellarSolver exists.
157 
158  // Let's sort edges, starting with widest
159  if (runHFR)
160  std::sort(stars.begin(), stars.end(), [](const FITSImage::Star & star1, const FITSImage::Star & star2) -> bool { return star1.HFR > star2.HFR;});
161  else
162  std::sort(stars.begin(), stars.end(), [](const FITSImage::Star & star1, const FITSImage::Star & star2) -> bool { return star1.flux > star2.flux;});
163 
164 
165  // Take only the first maxNumCenters stars
166  int starCount = qMin(maxStarsCount, stars.count());
167  starCenters.reserve(starCount);
168  for (int i = 0; i < starCount; i++)
169  {
170  Edge *oneEdge = new Edge();
171  oneEdge->x = stars[i].x;
172  oneEdge->y = stars[i].y;
173  oneEdge->val = stars[i].peak;
174  oneEdge->sum = stars[i].flux;
175  oneEdge->HFR = stars[i].HFR;
176  oneEdge->width = stars[i].a;
177  if (stars[i].a > 0)
178  // See page 63 to find the ellipticity equation for SEP.
179  // http://astroa.physics.metu.edu.tr/MANUALS/sextractor/Guide2source_extractor.pdf
180  oneEdge->ellipticity = 1 - stars[i].b / stars[i].a;
181  else
182  oneEdge->ellipticity = 0;
183 
184  starCenters.append(oneEdge);
185  }
186  m_ImageData->setStarCenters(starCenters);
187  return true;
188 #endif
189 }
190 
191 template <typename T>
192 void FITSSEPDetector::getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const *data) const
193 {
194  auto * rawBuffer = reinterpret_cast<T const *>(data->getImageBuffer());
195 
196  if (buffer == nullptr)
197  return;
198 
199  float * floatPtr = buffer;
200 
201  int x2 = x + w;
202  int y2 = y + h;
203 
204  FITSImage::Statistic const &stats = data->getStatistics();
205 
206  for (int y1 = y; y1 < y2; y1++)
207  {
208  int offset = y1 * stats.width;
209  for (int x1 = x; x1 < x2; x1++)
210  {
211  *floatPtr++ = rawBuffer[offset + x1];
212  }
213  }
214 }
215 
216 SkyBackground::SkyBackground(double mean_, double sigma_, double numPixels_)
217 {
218  initialize(mean_, sigma_, numPixels_);
219 }
220 
221 void SkyBackground::initialize(double mean_, double sigma_,
222  double numPixelsInSkyEstimate_, int numStars_)
223 {
224  mean = mean_;
225  sigma = sigma_;
226  numPixelsInSkyEstimate = numPixelsInSkyEstimate_;
227  varSky = sigma_ * sigma_;
228  starsDetected = numStars_;
229 }
230 
231 // Taken from: http://www1.phys.vt.edu/~jhs/phys3154/snr20040108.pdf
232 double SkyBackground::SNR(double flux, double numPixels, double gain) const
233 {
234  if (numPixelsInSkyEstimate <= 0 || gain <= 0)
235  return 0;
236  // const double numer = flux - numPixels * mean;
237  const double numer = flux; // It seems SEP flux subtracts out the background.
238  const double denom = sqrt( numer / gain + numPixels * varSky * (1 + 1 / numPixelsInSkyEstimate));
239  if (denom <= 0)
240  return 0;
241  return numer / denom;
242 
243 }
void append(const T &value)
QFuture< T > run(Function function,...)
int count(const T &value) const const
void reserve(int alloc)
void initialize(StandardShortcut id)
bool empty() const const
bool isValid() const const
QString filePath(const QString &fileName) const const
QList::iterator begin()
QList::iterator end()
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Thu Aug 11 2022 03:59:58 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.