Kstars

fitsgradientdetector.cpp
1 /*
2  SPDX-FileCopyrightText: 2004 Jasem Mutlaq
3  SPDX-FileCopyrightText: 2020 Eric Dejouhanet <[email protected]>
4  SPDX-FileCopyrightText: 2015 Gonzalo Exequiel Pedone <[email protected]>
5 
6  SPDX-License-Identifier: GPL-2.0-or-later AND GPL-3.0-or-later
7 
8  Some code fragments were adapted from Peter Kirchgessner's FITS plugin
9  SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
10 */
11 
12 #include <math.h>
13 #include <cmath>
14 #include <QtConcurrent>
15 
16 #include "fits_debug.h"
17 #include "fitsgradientdetector.h"
18 #include "fitsdata.h"
19 
20 QFuture<bool> FITSGradientDetector::findSources(const QRect &boundary)
21 {
22  FITSImage::Statistic const &stats = m_ImageData->getStatistics();
23  switch (stats.dataType)
24  {
25 
26  case TBYTE:
27  default:
28  return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint8_t>, boundary);
29 
30  case TSHORT:
31  return QtConcurrent::run(this, &FITSGradientDetector::findSources<int16_t>, boundary);
32 
33  case TUSHORT:
34  return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint16_t>, boundary);
35 
36  case TLONG:
37  return QtConcurrent::run(this, &FITSGradientDetector::findSources<int32_t>, boundary);
38 
39  case TULONG:
40  return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint16_t>, boundary);
41 
42  case TFLOAT:
43  return QtConcurrent::run(this, &FITSGradientDetector::findSources<float>, boundary);
44 
45  case TLONGLONG:
46  return QtConcurrent::run(this, &FITSGradientDetector::findSources<int64_t>, boundary);
47 
48  case TDOUBLE:
49  return QtConcurrent::run(this, &FITSGradientDetector::findSources<double>, boundary);
50  }
51 }
52 
53 template <typename T>
54 bool FITSGradientDetector::findSources(const QRect &boundary)
55 {
56  int subX = qMax(0, boundary.isNull() ? 0 : boundary.x());
57  int subY = qMax(0, boundary.isNull() ? 0 : boundary.y());
58  int subW = (boundary.isNull() ? m_ImageData->width() : boundary.width());
59  int subH = (boundary.isNull() ? m_ImageData->height() : boundary.height());
60 
61  int BBP = m_ImageData->getBytesPerPixel();
62 
63  uint16_t dataWidth = m_ImageData->width();
64 
65  // #1 Find offsets
66  uint32_t size = subW * subH;
67  uint32_t offset = subX + subY * dataWidth;
68 
69  // #2 Create new buffer
70  auto * buffer = new uint8_t[size * BBP];
71  // If there is no offset, copy whole buffer in one go
72  if (offset == 0)
73  memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP);
74  else
75  {
76  uint8_t * dataPtr = buffer;
77  uint8_t const * origDataPtr = m_ImageData->getImageBuffer();
78  uint32_t lineOffset = 0;
79  // Copy data line by line
80  for (int height = subY; height < (subY + subH); height++)
81  {
82  lineOffset = (subX + height * dataWidth) * BBP;
83  memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
84  dataPtr += (subW * BBP);
85  }
86  }
87 
88  // #3 Create new FITSData to hold it
89  auto * boundedImage = new FITSData();
90  FITSImage::Statistic stats;
91  stats.width = subW;
92  stats.height = subH;
93  stats.dataType = m_ImageData->getStatistics().dataType;
94  stats.bytesPerPixel = m_ImageData->getStatistics().bytesPerPixel;
95  stats.samples_per_channel = size;
96  stats.ndim = 2;
97  boundedImage->restoreStatistics(stats);
98 
99  // #4 Set image buffer and calculate stats.
100  boundedImage->setImageBuffer(buffer);
101 
102  boundedImage->calculateStats(true);
103 
104  // #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain
105  boundedImage->applyFilter(FITS_MEDIAN);
106  boundedImage->applyFilter(FITS_HIGH_CONTRAST);
107 
108  // #6 Perform Sobel to find gradients and their directions
109  QVector<float> gradients;
110  QVector<float> directions;
111 
112  // TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed
113  // and discarded whenever necessary. It won't work on noisy images unless this is done.
114  sobel<T>(boundedImage, gradients, directions);
115 
116  QVector<int> ids(gradients.size());
117 
118  int maxID = partition(subW, subH, gradients, ids);
119 
120  // Not needed anymore
121  delete boundedImage;
122 
123  if (maxID == 0)
124  return 0;
125 
126  typedef struct
127  {
128  float massX = 0;
129  float massY = 0;
130  float totalMass = 0;
131  } massInfo;
132 
133  QMap<int, massInfo> masses;
134 
135  // #7 Calculate center of mass for all detected regions
136  for (int y = 0; y < subH; y++)
137  {
138  for (int x = 0; x < subW; x++)
139  {
140  int index = x + y * subW;
141 
142  int regionID = ids[index];
143  if (regionID > 0)
144  {
145  float pixel = gradients[index];
146 
147  masses[regionID].totalMass += pixel;
148  masses[regionID].massX += x * pixel;
149  masses[regionID].massY += y * pixel;
150  }
151  }
152  }
153 
154  // Compare multiple masses, and only select the highest total mass one as the desired star
155  int maxRegionID = 1;
156  int maxTotalMass = masses[1].totalMass;
157  double totalMassRatio = 1e6;
158  for (auto key : masses.keys())
159  {
160  massInfo oneMass = masses.value(key);
161  if (oneMass.totalMass > maxTotalMass)
162  {
163  totalMassRatio = oneMass.totalMass / maxTotalMass;
164  maxTotalMass = oneMass.totalMass;
165  maxRegionID = key;
166  }
167  }
168 
169  // If image has many regions and there is no significant relative center of mass then it's just noise and no stars
170  // are probably there above a useful threshold.
171  if (maxID > 10 && totalMassRatio < 1.5)
172  return false;
173 
174  auto * center = new Edge;
175  center->width = -1;
176  center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
177  center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
178  center->HFR = 1;
179 
180  // Maximum Radius
181  int maxR = qMin(subW - 1, subH - 1) / 2;
182 
183  for (int r = maxR; r > 1; r--)
184  {
185  int pass = 0;
186 
187  for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
188  {
189  int testX = center->x + std::cos(theta) * r;
190  int testY = center->y + std::sin(theta) * r;
191 
192  // if out of bound, break;
193  if (testX < 0 || testX >= subW || testY < 0 || testY >= subH)
194  break;
195 
196  if (gradients[testX + testY * subW] > 0)
197  //if (thresholded[testX + testY * subW] > 0)
198  {
199  if (++pass >= 24)
200  {
201  center->width = r * 2;
202  // Break of outer loop
203  r = 0;
204  break;
205  }
206  }
207  }
208  }
209 
210  qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width;
211 
212  // If no stars were detected
213  if (center->width == -1)
214  {
215  delete center;
216  return false;
217  }
218 
219  // 30% fuzzy
220  //center->width += center->width*0.3 * (running_threshold / threshold);
221 
222  double FSum = 0, HF = 0, TF = 0;
223  const double resolution = 1.0 / 20.0;
224 
225  int cen_y = qRound(center->y);
226 
227  double rightEdge = center->x + center->width / 2.0;
228  double leftEdge = center->x - center->width / 2.0;
229 
230  QVector<double> subPixels;
231  subPixels.reserve(center->width / resolution);
232 
233  const T * origBuffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer()) + offset;
234 
235  for (double x = leftEdge; x <= rightEdge; x += resolution)
236  {
237  double slice = resolution * (origBuffer[static_cast<int>(floor(x)) + cen_y * dataWidth]);
238  FSum += slice;
239  subPixels.append(slice);
240  }
241 
242  // Half flux
243  HF = FSum / 2.0;
244 
245  int subPixelCenter = (center->width / resolution) / 2;
246 
247  // Start from center
248  TF = subPixels[subPixelCenter];
249  double lastTF = TF;
250  // Integrate flux along radius axis until we reach half flux
251  //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
252  for (int k = 1; k < subPixelCenter; k++)
253  {
254  TF += subPixels[subPixelCenter + k];
255  TF += subPixels[subPixelCenter - k];
256 
257  if (TF >= HF)
258  {
259  // We overpassed HF, let's calculate from last TF how much until we reach HF
260 
261  // #1 Accurate calculation, but very sensitive to small variations of flux
262  center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
263 
264  // #2 Less accurate calculation, but stable against small variations of flux
265  //center->HFR = (k - 1) * resolution;
266  break;
267  }
268 
269  lastTF = TF;
270  }
271 
272  // Correct center for subX and subY
273  center->x += subX;
274  center->y += subY;
275 
276  QList<Edge*> starCenters;
277  starCenters.append(center);
278  m_ImageData->setStarCenters(starCenters);
279 
280  qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR;
281 
282  return true;
283 }
284 
285 /* CannyDetector, Implementation of Canny edge detector in Qt/C++.
286  * Web-Site: https://github.com/hipersayanX/CannyDetector
287  */
288 
289 template <typename T>
290 void FITSGradientDetector::sobel(FITSData const *data, QVector<float> &gradient, QVector<float> &direction) const
291 {
292  if (data == nullptr)
293  return;
294 
295  FITSImage::Statistic const &stats = data->getStatistics();
296 
297  //int size = image.width() * image.height();
298  gradient.resize(stats.samples_per_channel);
299  direction.resize(stats.samples_per_channel);
300 
301  for (int y = 0; y < stats.height; y++)
302  {
303  size_t yOffset = y * stats.width;
304  const T * grayLine = reinterpret_cast<T const *>(data->getImageBuffer()) + yOffset;
305 
306  const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width;
307  const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width;
308 
309  float * gradientLine = gradient.data() + yOffset;
310  float * directionLine = direction.data() + yOffset;
311 
312  for (int x = 0; x < stats.width; x++)
313  {
314  int x_m1 = x < 1 ? x : x - 1;
315  int x_p1 = x >= stats.width - 1 ? x : x + 1;
316 
317  int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] -
318  2 * grayLine[x_m1] - grayLine_p1[x_m1];
319 
320  int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] -
321  2 * grayLine_p1[x] - grayLine_p1[x_p1];
322 
323  gradientLine[x] = qAbs(gradX) + qAbs(gradY);
324 
325  /* Gradient directions are classified in 4 possible cases
326  *
327  * dir 0
328  *
329  * x x x
330  * - - -
331  * x x x
332  *
333  * dir 1
334  *
335  * x x /
336  * x / x
337  * / x x
338  *
339  * dir 2
340  *
341  * \ x x
342  * x \ x
343  * x x \
344  *
345  * dir 3
346  *
347  * x | x
348  * x | x
349  * x | x
350  */
351  if (gradX == 0 && gradY == 0)
352  directionLine[x] = 0;
353  else if (gradX == 0)
354  directionLine[x] = 3;
355  else
356  {
357  qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
358 
359  if (a >= -22.5 && a < 22.5)
360  directionLine[x] = 0;
361  else if (a >= 22.5 && a < 67.5)
362  directionLine[x] = 2;
363  else if (a >= -67.5 && a < -22.5)
364  directionLine[x] = 1;
365  else
366  directionLine[x] = 3;
367  }
368  }
369  }
370 }
371 
372 int FITSGradientDetector::partition(int width, int height, QVector<float> &gradient, QVector<int> &ids) const
373 {
374  int id = 0;
375 
376  for (int y = 1; y < height - 1; y++)
377  {
378  for (int x = 1; x < width - 1; x++)
379  {
380  int index = x + y * width;
381  float val = gradient[index];
382  if (val > 0 && ids[index] == 0)
383  {
384  trace(width, height, ++id, gradient, ids, x, y);
385  }
386  }
387  }
388 
389  // Return max id
390  return id;
391 }
392 
393 void FITSGradientDetector::trace(int width, int height, int id, QVector<float> &image, QVector<int> &ids, int x,
394  int y) const
395 {
396  int yOffset = y * width;
397  float * cannyLine = image.data() + yOffset;
398  int * idLine = ids.data() + yOffset;
399 
400  if (idLine[x] != 0)
401  return;
402 
403  idLine[x] = id;
404 
405  for (int j = -1; j < 2; j++)
406  {
407  int nextY = y + j;
408 
409  if (nextY < 0 || nextY >= height)
410  continue;
411 
412  float * cannyLineNext = cannyLine + j * width;
413 
414  for (int i = -1; i < 2; i++)
415  {
416  int nextX = x + i;
417 
418  if (i == j || nextX < 0 || nextX >= width)
419  continue;
420 
421  if (cannyLineNext[nextX] > 0)
422  {
423  // Trace neighbors.
424  trace(width, height, id, image, ids, nextX, nextY);
425  }
426  }
427  }
428 }
void append(const T &value)
QFuture< T > run(Function function,...)
const T value(const Key &key, const T &defaultValue) const const
void append(const T &value)
int width() const const
int x() const const
int y() const const
T * data()
bool isNull() const const
void resize(int size)
void reserve(int size)
int height() const const
QList< Key > keys() const const
QTextStream & center(QTextStream &stream)
int size() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:53 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.