Kstars

fitscentroiddetector.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 "fitscentroiddetector.h"
16#include "fits_debug.h"
17#include "fitsdata.h"
18
19//void FITSCentroidDetector::configure(const QString &setting, const QVariant &value)
20//{
21// if (!setting.compare("MINIMUM_STDVAR", Qt::CaseInsensitive))
22// if (value.canConvert <int> ())
23// MINIMUM_STDVAR = value.value <int> ();
24
25// if (!setting.compare("MINIMUM_PIXEL_RANGE", Qt::CaseInsensitive))
26// if (value.canConvert <int> ())
27// MINIMUM_PIXEL_RANGE = value.value <int> ();
28
29// if (!setting.compare("JMINDEX", Qt::CaseInsensitive))
30// if (value.canConvert <double> ())
31// JMINDEX = value.value <double> ();
32//}
33
34bool FITSCentroidDetector::checkCollision(Edge * s1, Edge * s2) const
35{
36 int dis; //distance
37
38 int diff_x = s1->x - s2->x;
39 int diff_y = s1->y - s2->y;
40
41 dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y));
42 dis -= s1->width / 2;
43 dis -= s2->width / 2;
44
45 if (dis <= 0) //collision
46 return true;
47
48 //no collision
49 return false;
50}
51
52/*** Find center of stars and calculate Half Flux Radius */
53QFuture<bool> FITSCentroidDetector::findSources(const QRect &boundary)
54{
55 FITSImage::Statistic const &stats = m_ImageData->getStatistics();
56 switch (stats.dataType)
57 {
58#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
59 case TBYTE:
60 default:
61 return QtConcurrent::run(&FITSCentroidDetector::findSources<uint8_t const>, this, boundary);
62
63 case TSHORT:
64 return QtConcurrent::run(&FITSCentroidDetector::findSources<int16_t const>, this, boundary);
65
66 case TUSHORT:
67 return QtConcurrent::run(&FITSCentroidDetector::findSources<uint16_t const>, this, boundary);
68
69 case TLONG:
70 return QtConcurrent::run(&FITSCentroidDetector::findSources<int32_t const>, this, boundary);
71
72 case TULONG:
73 return QtConcurrent::run(&FITSCentroidDetector::findSources<uint32_t const>, this, boundary);
74
75 case TFLOAT:
76 return QtConcurrent::run(&FITSCentroidDetector::findSources<float const>, this, boundary);
77
78 case TLONGLONG:
79 return QtConcurrent::run(&FITSCentroidDetector::findSources<int64_t const>, this, boundary);
80
81 case TDOUBLE:
82 return QtConcurrent::run(&FITSCentroidDetector::findSources<double const>, this, boundary);
83#else
84 case TBYTE:
85 default:
86 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint8_t const>, boundary);
87
88 case TSHORT:
89 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int16_t const>, boundary);
90
91 case TUSHORT:
92 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint16_t const>, boundary);
93
94 case TLONG:
95 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int32_t const>, boundary);
96
97 case TULONG:
98 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint32_t const>, boundary);
99
100 case TFLOAT:
101 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<float const>, boundary);
102
103 case TLONGLONG:
104 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int64_t const>, boundary);
105
106 case TDOUBLE:
107 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<double const>, boundary);
108#endif
109 }
110}
111
112template <typename T>
113bool FITSCentroidDetector::findSources(const QRect &boundary)
114{
115 FITSImage::Statistic const &stats = m_ImageData->getStatistics();
116 FITSMode const m_Mode = static_cast<FITSMode>(m_ImageData->property("mode").toInt());
117
118 int MINIMUM_STDVAR = getValue("MINIMUM_STDVAR", 5).toInt();
119 int minEdgeWidth = getValue("MINIMUM_PIXEL_RANGE", 5).toInt();
120 double JMIndex = getValue("JMINDEX", 100.0).toDouble();
121
122 int initStdDev = MINIMUM_STDVAR;
123 double threshold = 0, sum = 0, avg = 0, min = 0;
124 int starDiameter = 0;
125 int pixVal = 0;
126 int minimumEdgeCount = MINIMUM_EDGE_LIMIT;
127
128 auto * buffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer());
129
130 float dispersion_ratio = 1.5;
131
132 QList<Edge *> edges;
133
134 if (JMIndex < DIFFUSE_THRESHOLD)
135 {
136 minEdgeWidth = JMIndex * 35 + 1;
137 minimumEdgeCount = minEdgeWidth - 1;
138 }
139 else
140 {
141 minEdgeWidth = 6;
142 minimumEdgeCount = 4;
143 }
144
145 while (initStdDev >= 1)
146 {
147 minEdgeWidth--;
148 minimumEdgeCount--;
149
150 minEdgeWidth = qMax(3, minEdgeWidth);
151 minimumEdgeCount = qMax(3, minimumEdgeCount);
152
153 if (JMIndex < DIFFUSE_THRESHOLD)
154 {
155 // Taking the average out seems to have better result for noisy images
156 threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1);
157
158 min = stats.min[0];
159 if (threshold - min < 0)
160 {
161 threshold = stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1);
162 min = 0;
163 }
164
165 dispersion_ratio = 1.4 - (MINIMUM_STDVAR - initStdDev) * 0.08;
166 }
167 else
168 {
169 threshold = stats.mean[0] + stats.stddev[0] * initStdDev * (0.3 - (MINIMUM_STDVAR - initStdDev) * 0.05);
170 min = stats.min[0];
171 // Ratio between centeroid center and edge
172 dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2;
173 }
174
175 qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR;
176 qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min
177 << ") minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount;
178
179 threshold -= min;
180
181 int subX, subY, subW, subH;
182
183 if (boundary.isNull())
184 {
185 if (m_Mode == FITS_GUIDE || m_Mode == FITS_FOCUS)
186 {
187 // Only consider the central 70%
188 subX = round(stats.width * 0.15);
189 subY = round(stats.height * 0.15);
190 subW = stats.width - subX;
191 subH = stats.height - subY;
192 }
193 else
194 {
195 // Consider the complete area 100%
196 subX = 0;
197 subY = 0;
198 subW = stats.width;
199 subH = stats.height;
200 }
201 }
202 else
203 {
204 subX = boundary.x();
205 subY = boundary.y();
206 subW = subX + boundary.width();
207 subH = subY + boundary.height();
208 }
209
210 // Detect "edges" that are above threshold
211 for (int i = subY; i < subH; i++)
212 {
213 starDiameter = 0;
214
215 for (int j = subX; j < subW; j++)
216 {
217 pixVal = buffer[j + (i * stats.width)] - min;
218
219 // If pixel value > threshold, let's get its weighted average
220 if (pixVal >= threshold)
221 {
222 avg += j * pixVal;
223 sum += pixVal;
224 starDiameter++;
225 }
226 // Value < threshold but avg exists
227 else if (sum > 0)
228 {
229 // We found a potential centroid edge
230 if (starDiameter >= minEdgeWidth)
231 {
232 float center = avg / sum + 0.5;
233 if (center > 0)
234 {
235 int i_center = std::floor(center);
236
237 // Check if center is 10% or more brighter than edge, if not skip
238 if (((buffer[i_center + (i * stats.width)] - min) /
239 (buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >=
240 dispersion_ratio) &&
241 ((buffer[i_center + (i * stats.width)] - min) /
242 (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >=
243 dispersion_ratio))
244 {
245 qCDebug(KSTARS_FITS)
246 << "Edge center is " << buffer[i_center + (i * stats.width)] - min
247 << " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min
248 << " and ratio is "
249 << ((buffer[i_center + (i * stats.width)] - min) /
250 (buffer[i_center + (i * stats.width) - starDiameter / 2] - min))
251 << " located at X: " << center << " Y: " << i + 0.5;
252
253 auto * newEdge = new Edge();
254
255 newEdge->x = center;
256 newEdge->y = i + 0.5;
257 newEdge->scanned = 0;
258 newEdge->val = buffer[i_center + (i * stats.width)] - min;
259 newEdge->width = starDiameter;
260 newEdge->HFR = 0;
261 newEdge->sum = sum;
262
263 edges.append(newEdge);
264 }
265 }
266 }
267
268 // Reset
269 avg = sum = starDiameter = 0;
270 }
271 }
272 }
273
274 qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count();
275
276 // In case of hot pixels
277 if (edges.count() == 1 && initStdDev > 1)
278 {
279 initStdDev--;
280 continue;
281 }
282
283 if (edges.count() >= MAX_EDGE_LIMIT)
284 {
285 qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count();
286 qDeleteAll(edges);
287 return -1;
288 }
289
290 if (edges.count() >= minimumEdgeCount)
291 break;
292
293 qDeleteAll(edges);
294 edges.clear();
295 initStdDev--;
296 }
297
298 int cen_count = 0;
299 int cen_x = 0;
300 int cen_y = 0;
301 int cen_v = 0;
302 int cen_w = 0;
303 int width_sum = 0;
304
305 // Let's sort edges, starting with widest
306 auto const greaterThan = [](Edge const * a, Edge const * b)
307 {
308 return a->sum > b->sum;
309 };
310 std::sort(edges.begin(), edges.end(), greaterThan);
311
312 QList<Edge*> starCenters;
313 // Now, let's scan the edges and find the maximum centroid vertically
314 for (int i = 0; i < edges.count(); i++)
315 {
316 qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of "
317 << edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum;
318
319 // If edge scanned already, skip
320 if (edges[i]->scanned == 1)
321 {
322 qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted";
323 continue;
324 }
325
326 qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ...";
327
328 // Get X, Y, and Val of edge
329 cen_x = edges[i]->x;
330 cen_y = edges[i]->y;
331 cen_v = edges[i]->sum;
332 cen_w = edges[i]->width;
333
334 float avg_x = 0;
335 float avg_y = 0;
336
337 sum = 0;
338 cen_count = 0;
339
340 // Now let's compare to other edges until we hit a maxima
341 for (int j = 0; j < edges.count(); j++)
342 {
343 if (edges[j]->scanned)
344 continue;
345
346 if (checkCollision(edges[j], edges[i]))
347 {
348 if (edges[j]->sum >= cen_v)
349 {
350 cen_v = edges[j]->sum;
351 cen_w = edges[j]->width;
352 }
353
354 edges[j]->scanned = 1;
355 cen_count++;
356
357 avg_x += edges[j]->x * edges[j]->val;
358 avg_y += edges[j]->y * edges[j]->val;
359 sum += edges[j]->val;
360
361 continue;
362 }
363 }
364
365 int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev));
366
367 if (edges.count() < LOW_EDGE_CUTOFF_1)
368 {
369 if (edges.count() < LOW_EDGE_CUTOFF_2)
370 cen_limit = 1;
371 else
372 cen_limit = 2;
373 }
374
375 qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is "
376 << cen_limit;
377
378 if (cen_limit < 1)
379 continue;
380
381 // If centroid count is within acceptable range
382 //if (cen_limit >= 2 && cen_count >= cen_limit)
383 if (cen_count >= cen_limit)
384 {
385 // We detected a centroid, let's init it
386 auto * rCenter = new Edge();
387
388 rCenter->x = avg_x / sum;
389 rCenter->y = avg_y / sum;
390 width_sum += rCenter->width;
391 rCenter->width = cen_w;
392
393 qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")";
394
395 // Calculate Total Flux From Center, Half Flux, Full Summation
396 double TF = 0;
397 double HF = 0;
398 double FSum = 0;
399
400 cen_x = (int)std::floor(rCenter->x);
401 cen_y = (int)std::floor(rCenter->y);
402
403 if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height)
404 {
405 delete rCenter;
406 continue;
407 }
408
409 // Complete sum along the radius
410 //for (int k=0; k < rCenter->width; k++)
411 for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--)
412 {
413 FSum += buffer[cen_x - k + (cen_y * stats.width)] - min;
414 //qDebug() << Q_FUNC_INFO << image_buffer[cen_x-k+(cen_y*stats.width)] - min;
415 }
416
417 // Half flux
418 HF = FSum / 2.0;
419
420 // Total flux starting from center
421 TF = buffer[cen_y * stats.width + cen_x] - min;
422
423 int pixelCounter = 1;
424
425 // Integrate flux along radius axis until we reach half flux
426 for (int k = 1; k < rCenter->width / 2; k++)
427 {
428 if (TF >= HF)
429 {
430 qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels.";
431 break;
432 }
433
434 TF += buffer[cen_y * stats.width + cen_x + k] - min;
435 TF += buffer[cen_y * stats.width + cen_x - k] - min;
436
437 pixelCounter++;
438 }
439
440 // Calculate weighted Half Flux Radius
441 rCenter->HFR = pixelCounter * (HF / TF);
442 // Store full flux
443 rCenter->val = FSum;
444
445 qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum;
446
447 starCenters.append(rCenter);
448 }
449 }
450
451 if (starCenters.count() > 1 && m_Mode != FITS_FOCUS)
452 {
453 float width_avg = (float)width_sum / starCenters.count();
454 float lsum = 0, sdev = 0;
455
456 for (auto &center : starCenters)
457 lsum += (center->width - width_avg) * (center->width - width_avg);
458
459 sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4;
460
461 // Reject stars > 4 * stddev
462 foreach (Edge * center, starCenters)
463 if (center->width > sdev)
464 starCenters.removeOne(center);
465
466 //foreach(Edge *center, starCenters)
467 //qDebug() << Q_FUNC_INFO << center->x << "," << center->y << "," << center->width << "," << center->val << Qt::endl;
468 }
469
470 m_ImageData->setStarCenters(starCenters);
471 // Release memory
472 qDeleteAll(edges);
473
474 return true;
475}
476
477
QList::append
void append(QList< T > &&value)
QList::begin
iterator begin()
QList::clear
void clear()
QList::count
qsizetype count() const const
QList::end
iterator end()
QObject::property
QVariant property(const char *name) const const
QRect::height
int height() const const
QRect::isNull
bool isNull() const const
QRect::width
int width() const const
QRect::x
int x() const const
QRect::y
int y() const const
Qt::center
QTextStream & center(QTextStream &stream)
QtConcurrent::run
QFuture< T > run(Function function,...)
QVariant::toDouble
double toDouble(bool *ok) const const
QVariant::toInt
int toInt(bool *ok) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Dec 13 2024 11:47:11 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.