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