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 return QtConcurrent::run(this, &FITSSEPDetector::findSourcesAndBackground, boundary);
59}
60
61bool 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);
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 oneEdge->numPixels = stars[i].numPixels;
178 if (stars[i].a > 0)
179 // See page 63 to find the ellipticity equation for SEP.
180 // http://astroa.physics.metu.edu.tr/MANUALS/sextractor/Guide2source_extractor.pdf
181 oneEdge->ellipticity = 1 - stars[i].b / stars[i].a;
182 else
183 oneEdge->ellipticity = 0;
184
185 starCenters.append(oneEdge);
186 }
187 m_ImageData->setStarCenters(starCenters);
188 return true;
189#endif
190}
191
192template <typename T>
193void FITSSEPDetector::getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const *data) const
194{
195 auto * rawBuffer = reinterpret_cast<T const *>(data->getImageBuffer());
196
197 if (buffer == nullptr)
198 return;
199
200 float * floatPtr = buffer;
201
202 int x2 = x + w;
203 int y2 = y + h;
204
205 FITSImage::Statistic const &stats = data->getStatistics();
206
207 for (int y1 = y; y1 < y2; y1++)
208 {
209 int offset = y1 * stats.width;
210 for (int x1 = x; x1 < x2; x1++)
211 {
212 *floatPtr++ = rawBuffer[offset + x1];
213 }
214 }
215}
216
217SkyBackground::SkyBackground(double mean_, double sigma_, double numPixels_)
218{
219 initialize(mean_, sigma_, numPixels_);
220}
221
222void SkyBackground::initialize(double mean_, double sigma_,
223 double numPixelsInSkyEstimate_, int numStars_)
224{
225 mean = mean_;
226 sigma = sigma_;
227 numPixelsInSkyEstimate = numPixelsInSkyEstimate_;
228 varSky = sigma_ * sigma_;
229 starsDetected = numStars_;
230}
231
232// Taken from: http://www1.phys.vt.edu/~jhs/phys3154/snr20040108.pdf
233double SkyBackground::SNR(double flux, double numPixels, double gain) const
234{
235 if (numPixelsInSkyEstimate <= 0 || gain <= 0)
236 return 0;
237 // const double numer = flux - numPixels * mean;
238 const double numer = flux; // It seems SEP flux subtracts out the background.
239 const double denom = sqrt( numer / gain + numPixels * varSky * (1 + 1 / numPixelsInSkyEstimate));
240 if (denom <= 0)
241 return 0;
242 return numer / denom;
243
244}
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 Oct 4 2024 12:01:03 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.