Kstars

fitsthresholddetector.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 "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 
26 QFuture<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 
59 template <typename T>
60 bool 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 
106  while (running_threshold >= critical_threshold)
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;
118  if (testX < subX || testX > subW || testY < subY || testY > subH)
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%
138  running_threshold -= running_threshold * 0.1;
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 
159  QVector<double> subPixels;
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
177  TF = subPixels[subPixelCenter];
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  {
183  TF += subPixels[subPixelCenter + k];
184  TF += subPixels[subPixelCenter - k];
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 }
void append(const T &value)
QFuture< T > run(Function function,...)
void append(const T &value)
int width() const const
int x() const const
int y() const const
void reserve(int size)
int height() const const
QTextStream & center(QTextStream &stream)
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 19 2022 03:57:50 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.