Kstars

fitscentroiddetector.cpp
1 /*
2  SPDX-FileCopyrightText: 2004 Jasem Mutlaq
3  SPDX-FileCopyrightText: 2020 Eric Dejouhanet <[email protected]>
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 
34 bool 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 */
53 QFuture<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 
86 template <typename T>
87 bool 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 
96  int initStdDev = MINIMUM_STDVAR;
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 
106  QList<Edge *> edges;
107 
108  if (JMIndex < DIFFUSE_THRESHOLD)
109  {
110  minEdgeWidth = JMIndex * 35 + 1;
111  minimumEdgeCount = minEdgeWidth - 1;
112  }
113  else
114  {
115  minEdgeWidth = 6;
116  minimumEdgeCount = 4;
117  }
118 
119  while (initStdDev >= 1)
120  {
121  minEdgeWidth--;
122  minimumEdgeCount--;
123 
124  minEdgeWidth = qMax(3, minEdgeWidth);
125  minimumEdgeCount = qMax(3, minimumEdgeCount);
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
146  dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2;
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
204  if (starDiameter >= minEdgeWidth)
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) >=
214  dispersion_ratio) &&
215  ((buffer[i_center + (i * stats.width)] - min) /
216  (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >=
217  dispersion_ratio))
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();
260  qDeleteAll(edges);
261  return -1;
262  }
263 
264  if (edges.count() >= minimumEdgeCount)
265  break;
266 
267  qDeleteAll(edges);
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
446  qDeleteAll(edges);
447 
448  return true;
449 }
450 
451 
void append(const T &value)
QFuture< T > run(Function function,...)
int count(const T &value) const const
int width() const const
int x() const const
int y() const const
bool isNull() const const
int height() const const
void clear()
QList::iterator begin()
QTextStream & center(QTextStream &stream)
QList::iterator end()
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Sun Aug 14 2022 04:13:55 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.