Kstars

fitsthresholddetector.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 <math.h>
12#include <cmath>
13#include <QtConcurrent>
14
15#include "fits_debug.h"
16#include "fitsthresholddetector.h"
17#include "fitsdata.h"
18
19//void FITSThresholdDetector::configure(const QString &setting, const QVariant &value)
20//{
21// if (!setting.compare("THRESHOLD_PERCENTAGE", Qt::CaseInsensitive))
22// if (value.canConvert<int>())
23// THRESHOLD_PERCENTAGE = value.value<int>();
24//}
25
26QFuture<bool> FITSThresholdDetector::findSources(QRect const &boundary)
27{
28 FITSImage::Statistic const &stats = m_ImageData->getStatistics();
29 switch (stats.dataType)
30 {
31 case TSHORT:
32 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int16_t>, boundary);
33
34 case TUSHORT:
35 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint16_t>, boundary);
36
37 case TLONG:
38 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int32_t>, boundary);
39
40 case TULONG:
41 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint32_t>, boundary);
42
43 case TFLOAT:
44 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<float>, boundary);
45
46 case TLONGLONG:
47 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int64_t>, boundary);
48
49 case TDOUBLE:
50 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<double>, boundary);
51
52 case TBYTE:
53 default:
54 return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint8_t>, boundary);
55 }
56
57}
58
59template <typename T>
60bool FITSThresholdDetector::findOneStar(const QRect &boundary) const
61{
62 FITSImage::Statistic const &stats = m_ImageData->getStatistics();
63
64 int subX = boundary.x();
65 int subY = boundary.y();
66 int subW = subX + boundary.width();
67 int subH = subY + boundary.height();
68
69 float massX = 0, massY = 0, totalMass = 0;
70
71 auto * buffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer());
72
73 double THRESHOLD_PERCENTAGE = getValue("THRESHOLD_PERCENTAGE", 120.0).toDouble();
74 // TODO replace magic number with something more useful to understand
75 double threshold = stats.mean[0] * THRESHOLD_PERCENTAGE / 100.0;
76
77 for (int y = subY; y < subH; y++)
78 {
79 for (int x = subX; x < subW; x++)
80 {
81 T pixel = buffer[x + y * stats.width];
82 if (pixel > threshold)
83 {
84 totalMass += pixel;
85 massX += x * pixel;
86 massY += y * pixel;
87 }
88 }
89 }
90
91 qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << massX / totalMass << " Y: " << massY / totalMass;
92
93 auto * center = new Edge;
94 center->width = -1;
95 center->x = massX / totalMass + 0.5;
96 center->y = massY / totalMass + 0.5;
97 center->HFR = 1;
98
99 // Maximum Radius
100 int maxR = qMin(subW - 1, subH - 1) / 2;
101
102 // Critical threshold
103 double critical_threshold = threshold * 0.7;
104 double running_threshold = threshold;
105
107 {
108 for (int r = maxR; r > 1; r--)
109 {
110 int pass = 0;
111
112 for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0)
113 {
114 int testX = center->x + std::cos(theta) * r;
115 int testY = center->y + std::sin(theta) * r;
116
117 // if out of bound, break;
119 break;
120
121 if (buffer[testX + testY * stats.width] > running_threshold)
122 pass++;
123 }
124
125 //qDebug() << Q_FUNC_INFO << "Testing for radius " << r << " passes # " << pass << " @ threshold " << running_threshold;
126 //if (pass >= 6)
127 if (pass >= 5)
128 {
129 center->width = r * 2;
130 break;
131 }
132 }
133
134 if (center->width > 0)
135 break;
136
137 // Increase threshold fuzziness by 10%
139 }
140
141 // If no stars were detected
142 if (center->width == -1)
143 {
144 delete center;
145 return false;
146 }
147
148 // 30% fuzzy
149 //center->width += center->width*0.3 * (running_threshold / threshold);
150
151 double FSum = 0, HF = 0, TF = 0, min = stats.min[0];
152 const double resolution = 1.0 / 20.0;
153
154 int cen_y = qRound(center->y);
155
156 double rightEdge = center->x + center->width / 2.0;
157 double leftEdge = center->x - center->width / 2.0;
158
160 subPixels.reserve(center->width / resolution);
161
162 for (double x = leftEdge; x <= rightEdge; x += resolution)
163 {
164 //subPixels[x] = resolution * (image_buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
165 double slice = resolution * (buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
166 FSum += slice;
167 subPixels.append(slice);
168 }
169
170 // Half flux
171 HF = FSum / 2.0;
172
173 //double subPixelCenter = center->x - fmod(center->x,resolution);
174 int subPixelCenter = (center->width / resolution) / 2;
175
176 // Start from center
178 double lastTF = TF;
179 // Integrate flux along radius axis until we reach half flux
180 //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
181 for (int k = 1; k < subPixelCenter; k++)
182 {
185
186 if (TF >= HF)
187 {
188 // 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.
189 // The second method is not truly HFR but is much more resistant to noise.
190
191 // #1 Approximate HFR, accurate and reliable but quite variable to small changes in star flux
192 center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
193
194 // #2 Not exactly HFR, but much more stable
195 //center->HFR = (k*resolution) * (HF/TF);
196 break;
197 }
198
199 lastTF = TF;
200 }
201
202 QList<Edge*> starCenters;
203 starCenters.append(center);
204 m_ImageData->setStarCenters(starCenters);
205
206 return true;
207}
int height() const const
int width() const const
int x() const const
int y() const const
QTextStream & center(QTextStream &stream)
QFuture< T > run(Function function,...)
double toDouble(bool *ok) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri May 17 2024 11:48:26 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.