Differential D28292 Diff 79624 kstars/fitsviewer/fitsthresholddetector.cpp

Changeset View

Standalone View

View Options

kstars/fitsviewer/fitsthresholddetector.cpp

  • This file was added.
1​/***************************************************************************
2​ fitsthresholddetector.cpp - FITS Image
3​ -------------------
4​ begin : Sat March 28 2020
5​ copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet
6​ email : eric.dejouhanet@gmail.com
7​ ***************************************************************************/
8​
9​/***************************************************************************
10​ * *
11​ * This program is free software; you can redistribute it and/or modify *
12​ * it under the terms of the GNU General Public License as published by *
13​ * the Free Software Foundation; either version 2 of the License, or *
14​ * (at your option) any later version. *
15​ * *
16​ * Some code fragments were adapted from Peter Kirchgessner's FITS plugin*
17​ * See http://members.aol.com/pkirchg for more details. *
18​ ***************************************************************************/
19​
20​#include <math.h>
21​
22​#include "fits_debug.h"
23​#include "fitsthresholddetector.h"
24​
25​FITSStarDetector& FITSThresholdDetector::configure(const QString &setting, const QVariant &value)
26​{
27​ bool ok = false;
28​
29​ if (!setting.compare("threshold", Qt::CaseInsensitive))
30​ {
31​ double const _focusThreshold = value.toDouble(&ok);
32​ if (ok)
33​ focusThreshold = _focusThreshold;
34​ }
35​
36​ return *this;
37​}
38​
39​int FITSThresholdDetector::findSources(QList<Edge*> &starCenters, QRect const &boundary)
40​{
41​ switch (parent()->property("dataType").toInt())
42​ {
43​ case TBYTE:
44​ return findOneStar<uint8_t>(starCenters, boundary);
45​
46​ case TSHORT:
47​ return findOneStar<int16_t>(starCenters, boundary);
48​
49​ case TUSHORT:
50​ return findOneStar<uint16_t>(starCenters, boundary);
51​
52​ case TLONG:
53​ return findOneStar<int32_t>(starCenters, boundary);
54​
55​ case TULONG:
56​ return findOneStar<uint32_t>(starCenters, boundary);
57​
58​ case TFLOAT:
59​ return findOneStar<float>(starCenters, boundary);
60​
61​ case TLONGLONG:
62​ return findOneStar<int64_t>(starCenters, boundary);
63​
64​ case TDOUBLE:
65​ return findOneStar<double>(starCenters, boundary);
66​
67​ default:
68​ break;
69​ }
70​
71​ return 0;
72​}
73​
74​template <typename T>
75​int FITSThresholdDetector::findOneStar(QList<Edge*> &starCenters, const QRect &boundary) const
76​{
77​ FITSData const * const image_data = reinterpret_cast<FITSData const *>(parent());
78​
79​ if (image_data == nullptr)
80​ return 0;
81​
82​ FITSData::Statistic const &stats = image_data->getStatistics();
83​
84​ if (boundary.isEmpty())
85​ return -1;
86​
87​ int subX = boundary.x();
88​ int subY = boundary.y();
89​ int subW = subX + boundary.width();
90​ int subH = subY + boundary.height();
91​
92​ float massX = 0, massY = 0, totalMass = 0;
93​
94​ auto * buffer = reinterpret_cast<T const *>(image_data->getImageBuffer());
95​
96​ // TODO replace magic number with something more useful to understand
97​ double threshold = stats.mean[0] * focusThreshold / 100.0;
98​
99​ for (int y = subY; y < subH; y++)
100​ {
101​ for (int x = subX; x < subW; x++)
102​ {
103​ T pixel = buffer[x + y * stats.width];
104​ if (pixel > threshold)
105​ {
106​ totalMass += pixel;
107​ massX += x * pixel;
108​ massY += y * pixel;
109​ }
110​ }
111​ }
112​
113​ qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << massX / totalMass << " Y: " << massY / totalMass;
114​
115​ auto * center = new Edge;
116​ center->width = -1;
117​ center->x = massX / totalMass + 0.5;
118​ center->y = massY / totalMass + 0.5;
119​ center->HFR = 1;
120​
121​ // Maximum Radius
122​ int maxR = qMin(subW - 1, subH - 1) / 2;
123​
124​ // Critical threshold
125​ double critical_threshold = threshold * 0.7;
126​ double running_threshold = threshold;
127​
128​ while (running_threshold >= critical_threshold)
129​ {
130​ for (int r = maxR; r > 1; r--)
131​ {
132​ int pass = 0;
133​
134​ for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0)
135​ {
136​ int testX = center->x + std::cos(theta) * r;
137​ int testY = center->y + std::sin(theta) * r;
138​
139​ // if out of bound, break;
140​ if (testX < subX || testX > subW || testY < subY || testY > subH)
141​ break;
142​
143​ if (buffer[testX + testY * stats.width] > running_threshold)
144​ pass++;
145​ }
146​
147​ //qDebug() << "Testing for radius " << r << " passes # " << pass << " @ threshold " << running_threshold;
148​ //if (pass >= 6)
149​ if (pass >= 5)
150​ {
151​ center->width = r * 2;
152​ break;
153​ }
154​ }
155​
156​ if (center->width > 0)
157​ break;
158​
159​ // Increase threshold fuzziness by 10%
160​ running_threshold -= running_threshold * 0.1;
161​ }
162​
163​ // If no stars were detected
164​ if (center->width == -1)
165​ {
166​ delete center;
167​ return 0;
168​ }
169​
170​ // 30% fuzzy
171​ //center->width += center->width*0.3 * (running_threshold / threshold);
172​
173​ double FSum = 0, HF = 0, TF = 0, min = stats.min[0];
174​ const double resolution = 1.0 / 20.0;
175​
176​ int cen_y = qRound(center->y);
177​
178​ double rightEdge = center->x + center->width / 2.0;
179​ double leftEdge = center->x - center->width / 2.0;
180​
181​ QVector<double> subPixels;
182​ subPixels.reserve(center->width / resolution);
183​
184​ for (double x = leftEdge; x <= rightEdge; x += resolution)
185​ {
186​ //subPixels[x] = resolution * (image_buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
187​ double slice = resolution * (buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
188​ FSum += slice;
189​ subPixels.append(slice);
190​ }
191​
192​ // Half flux
193​ HF = FSum / 2.0;
194​
195​ //double subPixelCenter = center->x - fmod(center->x,resolution);
196​ int subPixelCenter = (center->width / resolution) / 2;
197​
198​ // Start from center
199​ TF = subPixels[subPixelCenter];
200​ double lastTF = TF;
201​ // Integrate flux along radius axis until we reach half flux
202​ //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
203​ for (int k = 1; k < subPixelCenter; k++)
204​ {
205​ TF += subPixels[subPixelCenter + k];
206​ TF += subPixels[subPixelCenter - k];
207​
208​ if (TF >= HF)
209​ {
210​ // We have two ways to calculate HFR. The first is the correct method but it can get quite variable within 10% due to random fluctuations of the measured star.
211​ // The second method is not truly HFR but is much more resistant to noise.
212​
213​ // #1 Approximate HFR, accurate and reliable but quite variable to small changes in star flux
214​ center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
215​
216​ // #2 Not exactly HFR, but much more stable
217​ //center->HFR = (k*resolution) * (HF/TF);
218​ break;
219​ }
220​
221​ lastTF = TF;
222​ }
223​
224​ starCenters.append(center);
225​
226​ return 1;
227​}