Kstars

fitssepdetector.cpp
1/*
2 SPDX-FileCopyrightText: 2004 Jasem Mutlaq
3 SPDX-FileCopyrightText: 2020 Eric Dejouhanet <eric.dejouhanet@gmail.com>
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
56QFuture<bool> FITSSEPDetector::findSources(QRect const &boundary)
57{
58#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
59 return QtConcurrent::run(&FITSSEPDetector::findSourcesAndBackground, this, boundary);
60#else
61 return QtConcurrent::run(this, &FITSSEPDetector::findSourcesAndBackground, boundary);
62#endif
63}
64
65bool FITSSEPDetector::findSourcesAndBackground(QRect const &boundary)
66{
67#ifndef HAVE_STELLARSOLVER
68 Q_UNUSED(boundary)
69 return false;
70#else
71 QList<Edge*> starCenters;
72 SkyBackground skyBG;
73 int maxStarsCount = getValue("maxStarsCount", 100000).toInt();
74
75 int optionsProfileIndex = getValue("optionsProfileIndex", -1).toInt();
76 Ekos::ProfileGroup group = static_cast<Ekos::ProfileGroup>(getValue("optionsProfileGroup", 1).toInt());
77 QScopedPointer<StellarSolver, QScopedPointerDeleteLater> solver(new StellarSolver(m_ImageData->getStatistics(),
78 m_ImageData->getImageBuffer()));
79 QString filename = "";
80 QPointer<FITSData> image(m_ImageData);
81 switch(group)
82 {
83 case Ekos::AlignProfiles:
84 //So it should not be here if it is Align.
85 break;
86 case Ekos::GuideProfiles:
87 filename = "SavedGuideProfiles.ini";
88 break;
89 case Ekos::FocusProfiles:
90 filename = "SavedFocusProfiles.ini";
91 break;
92 case Ekos::HFRProfiles:
93 filename = "SavedHFRProfiles.ini";
94 break;
95 }
96
97 QString savedOptionsProfiles = QDir(KSPaths::writableLocation(QStandardPaths::AppLocalDataLocation)).filePath(filename);
99 if(QFile(savedOptionsProfiles).exists())
100 optionsList = StellarSolver::loadSavedOptionsProfiles(savedOptionsProfiles);
101 else
102 {
103 switch(group)
104 {
105 case Ekos::AlignProfiles:
106 optionsList = Ekos::getDefaultAlignOptionsProfiles();
107 break;
108 case Ekos::GuideProfiles:
109 optionsList = Ekos::getDefaultGuideOptionsProfiles();
110 break;
111 case Ekos::FocusProfiles:
112 optionsList = Ekos::getDefaultFocusOptionsProfiles();
113 break;
114 case Ekos::HFRProfiles:
115 optionsList = Ekos::getDefaultHFROptionsProfiles();
116 break;
117 }
118 }
119 if (optionsProfileIndex >= 0 && optionsList.count() > optionsProfileIndex)
120 {
121 auto params = optionsList[optionsProfileIndex];
122 params.partition = Options::stellarSolverPartition();
123 solver->setParameters(params);
124 qCDebug(KSTARS_FITS) << "Sextract with: " << optionsList[optionsProfileIndex].listName;
125 }
126 else
127 {
128 auto params = SSolver::Parameters(); // This is default
129 params.partition = Options::stellarSolverPartition();
130 solver->setParameters(params);
131 }
132
134 const bool runHFR = group != Ekos::AlignProfiles;
135
136 solver->setLogLevel(SSolver::LOG_NONE);
137 solver->setSSLogLevel(SSolver::LOG_OFF);
138
139 if (boundary.isValid())
140 solver->extract(runHFR, boundary);
141 else
142 solver->extract(runHFR);
143
144 stars = solver->getStarList();
145
146 // If m_ImageData goes out of scope, also return.
147 if (stars.empty() || image.isNull())
148 return false;
149
150 auto bg = solver->getBackground();
151
152 skyBG.mean = bg.global;
153 skyBG.sigma = bg.globalrms;
154 skyBG.numPixelsInSkyEstimate = bg.bw * bg.bh;
155 skyBG.setStarsDetected(bg.num_stars_detected);
156 m_ImageData->setSkyBackground(skyBG);
157
158 //There is more information that can be obtained by the Stellarsolver->
159 //Background info, Star positions(if a plate solve was done before), etc
160 //The information is available as long as the StellarSolver exists.
161
162 // Let's sort edges, starting with widest
163 if (runHFR)
164 std::sort(stars.begin(), stars.end(), [](const FITSImage::Star & star1, const FITSImage::Star & star2) -> bool { return star1.HFR > star2.HFR;});
165 else
166 std::sort(stars.begin(), stars.end(), [](const FITSImage::Star & star1, const FITSImage::Star & star2) -> bool { return star1.flux > star2.flux;});
167
168
169 // Take only the first maxNumCenters stars
170 int starCount = qMin(maxStarsCount, stars.count());
171 starCenters.reserve(starCount);
172 for (int i = 0; i < starCount; i++)
173 {
174 Edge *oneEdge = new Edge();
175 oneEdge->x = stars[i].x;
176 oneEdge->y = stars[i].y;
177 oneEdge->val = stars[i].peak;
178 oneEdge->sum = stars[i].flux;
179 oneEdge->HFR = stars[i].HFR;
180 oneEdge->width = stars[i].a;
181 oneEdge->numPixels = stars[i].numPixels;
182 if (stars[i].a > 0)
183 // See page 63 to find the ellipticity equation for SEP.
184 // http://astroa.physics.metu.edu.tr/MANUALS/sextractor/Guide2source_extractor.pdf
185 oneEdge->ellipticity = 1 - stars[i].b / stars[i].a;
186 else
187 oneEdge->ellipticity = 0;
188
189 starCenters.append(oneEdge);
190 }
191 m_ImageData->setStarCenters(starCenters);
192 return true;
193#endif
194}
195
196template <typename T>
197void FITSSEPDetector::getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const *data) const
198{
199 auto * rawBuffer = reinterpret_cast<T const *>(data->getImageBuffer());
200
201 if (buffer == nullptr)
202 return;
203
204 float * floatPtr = buffer;
205
206 int x2 = x + w;
207 int y2 = y + h;
208
209 FITSImage::Statistic const &stats = data->getStatistics();
210
211 for (int y1 = y; y1 < y2; y1++)
212 {
213 int offset = y1 * stats.width;
214 for (int x1 = x; x1 < x2; x1++)
215 {
216 *floatPtr++ = rawBuffer[offset + x1];
217 }
218 }
219}
220
221SkyBackground::SkyBackground(double mean_, double sigma_, double numPixels_)
222{
223 initialize(mean_, sigma_, numPixels_);
224}
225
226void SkyBackground::initialize(double mean_, double sigma_,
227 double numPixelsInSkyEstimate_, int numStars_)
228{
229 mean = mean_;
230 sigma = sigma_;
231 numPixelsInSkyEstimate = numPixelsInSkyEstimate_;
232 varSky = sigma_ * sigma_;
233 starsDetected = numStars_;
234}
235
236// Taken from: http://www1.phys.vt.edu/~jhs/phys3154/snr20040108.pdf
237double SkyBackground::SNR(double flux, double numPixels, double gain) const
238{
239 if (numPixelsInSkyEstimate <= 0 || gain <= 0)
240 return 0;
241 // const double numer = flux - numPixels * mean;
242 const double numer = flux; // It seems SEP flux subtracts out the background.
243 const double denom = sqrt( numer / gain + numPixels * varSky * (1 + 1 / numPixelsInSkyEstimate));
244 if (denom <= 0)
245 return 0;
246 return numer / denom;
247
248}
QString filePath(const QString &fileName) const const
void append(QList< T > &&value)
iterator begin()
qsizetype count() const const
bool empty() const const
iterator end()
void reserve(qsizetype size)
bool isValid() const const
QFuture< T > run(Function function,...)
int toInt(bool *ok) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Nov 29 2024 11:57:48 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.