Kstars

fitsbahtinovdetector.cpp
1 /*
2  SPDX-FileCopyrightText: 2020 Patrick Molenaar <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 
6  Some code fragments were adapted from Peter Kirchgessner's FITS plugin
7  SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
8 */
9 
10 #include <math.h>
11 #include <cmath>
12 
13 #include "fits_debug.h"
14 #include "fitsbahtinovdetector.h"
15 #include "hough/houghline.h"
16 #include "fitsdata.h"
17 
18 #include <QElapsedTimer>
19 #include <QtConcurrent>
20 
21 //void FITSBahtinovDetector::configure(const QString &setting, const QVariant &value)
22 //{
23 // if (!setting.compare("NUMBER_OF_AVERAGE_ROWS", Qt::CaseInsensitive))
24 // {
25 // if (value.canConvert <int> ())
26 // {
27 // NUMBER_OF_AVERAGE_ROWS = value.value <int> ();
28 
29 // // Validate number of average rows value
30 // if (NUMBER_OF_AVERAGE_ROWS % 2 == 0)
31 // {
32 // NUMBER_OF_AVERAGE_ROWS--;
33 // qCInfo(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to "
34 // << NUMBER_OF_AVERAGE_ROWS;
35 // }
36 // // Rows must be a positive number!
37 // if (NUMBER_OF_AVERAGE_ROWS < 1)
38 // {
39 // NUMBER_OF_AVERAGE_ROWS = 1;
40 // qCInfo(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to "
41 // << NUMBER_OF_AVERAGE_ROWS;
42 // }
43 // }
44 // }
45 //}
46 
47 QFuture<bool> FITSBahtinovDetector::findSources(QRect const &boundary)
48 {
49  FITSImage::Statistic const &stats = m_ImageData->getStatistics();
50  switch (stats.dataType)
51  {
52  case TSHORT:
53  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int16_t>, boundary);
54 
55  case TUSHORT:
56  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint16_t>, boundary);
57 
58  case TLONG:
59  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int32_t>, boundary);
60 
61  case TULONG:
62  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint32_t>, boundary);
63 
64  case TFLOAT:
65  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<float>, boundary);
66 
67  case TLONGLONG:
68  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int64_t>, boundary);
69 
70  case TDOUBLE:
71  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<double>, boundary);
72 
73  default:
74  case TBYTE:
75  return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint8_t>, boundary);
76 
77  }
78 }
79 
80 template <typename T>
81 bool FITSBahtinovDetector::findBahtinovStar(const QRect &boundary)
82 {
83  if (boundary.isEmpty())
84  return false;
85 
86  QList<Edge*> starCenters;
87  int subX = qMax(0, boundary.isNull() ? 0 : boundary.x());
88  int subY = qMax(0, boundary.isNull() ? 0 : boundary.y());
89  int subW = (boundary.isNull() ? m_ImageData->width() : boundary.width());
90  int subH = (boundary.isNull() ? m_ImageData->height() : boundary.height());
91 
92  int BBP = m_ImageData->getBytesPerPixel();
93  uint16_t dataWidth = m_ImageData->width();
94 
95  // #1 Find offsets
96  uint32_t size = subW * subH;
97  uint32_t offset = subX + subY * dataWidth;
98 
99  // #2 Create new buffer
100  auto * buffer = new uint8_t[size * BBP];
101  // If there is no offset, copy whole buffer in one go
102  if (offset == 0)
103  {
104  memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP);
105  }
106  else
107  {
108  uint8_t * dataPtr = buffer;
109  const uint8_t * origDataPtr = m_ImageData->getImageBuffer();
110  // Copy data line by line
111  for (int height = subY; height < (subY + subH); height++)
112  {
113  uint32_t lineOffset = (subX + height * dataWidth) * BBP;
114  memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
115  dataPtr += (subW * BBP);
116  }
117  }
118 
119  // #3 Create new FITSData to hold it
120  FITSImage::Statistic stats;
121  stats.width = subW;
122  stats.height = subH;
123  stats.dataType = m_ImageData->getStatistics().dataType;
124  stats.bytesPerPixel = BBP;
125  stats.samples_per_channel = size;
126  FITSData* boundedImage = new FITSData();
127  boundedImage->restoreStatistics(stats);
128  boundedImage->setProperty("dataType", stats.dataType);
129 
130  // #4 Set image buffer
131  boundedImage->setImageBuffer(buffer);
132 
133  boundedImage->calculateStats(true);
134 
135  QElapsedTimer timer1;
136 
137  // Rotate image 180 degrees in steps of 1 degree
138  QMap<int, BahtinovLineAverage> lineAveragesPerAngle;
139  const int steps = 180;
140  double radPerStep = M_PI / steps;
141 
142  timer1.start();
143 
144  for (int angle = 0; angle < steps; angle++)
145  {
146  // TODO Apply multi threading to speed up calculation
147  BahtinovLineAverage lineAverage = calculateMaxAverage<T>(boundedImage, angle);
148  // Store line average in map
149  lineAveragesPerAngle.insert(angle, lineAverage);
150  }
151 
152  qCDebug(KSTARS_FITS) << "Getting max average for all 180 rotations took" << timer1.elapsed() << "milliseconds";
153 
154  // Not needed anymore
155  delete boundedImage;
156 
157  // Calculate Bahtinov angles
158  QVector<HoughLine*> bahtinov_angles;
159 
160  // For all three Bahtinov angles
161  for (int index1 = 0; index1 < 3; index1++)
162  {
163  double maxAverage = 0.0;
164  double maxAngle = 0.0;
165  int maxAverageOffset = 0;
166  for (int angle = 0; angle < steps; angle++)
167  {
168  BahtinovLineAverage lineAverage = lineAveragesPerAngle[angle];
169  if (lineAverage.average > maxAverage)
170  {
171  maxAverage = lineAverage.average;
172  maxAverageOffset = lineAverage.offset;
173  maxAngle = angle;
174  }
175  }
176  HoughLine* pHoughLine = new HoughLine(maxAngle * radPerStep, maxAverageOffset, subW, subH, maxAverage);
177  if (pHoughLine != nullptr)
178  {
179  bahtinov_angles.append(pHoughLine);
180  }
181 
182  // Remove data around peak to prevent it from being detected again
183  int minBahtinovAngleOffset = 18;
184  for (int subAngle = maxAngle - minBahtinovAngleOffset; subAngle < maxAngle + minBahtinovAngleOffset; subAngle++)
185  {
186  int angleInRange = subAngle;
187  if (angleInRange < 0)
188  {
189  angleInRange += 180;
190  }
191  if (angleInRange >= 180)
192  {
193  angleInRange -= 180;
194  }
195  lineAveragesPerAngle.remove(angleInRange);
196  }
197  }
198 
199  // Proceed with focus offset calculation, but only when at least 3 lines have been detected
200  QVector<HoughLine*> top3Lines;
201  if (bahtinov_angles.size() >= 3)
202  {
203  HoughLine::getSortedTopThreeLines(bahtinov_angles, top3Lines);
204 
205  // Debug output
206  qCDebug(KSTARS_FITS) << "Sorted bahtinov angles:";
207  foreach (HoughLine* ln, top3Lines)
208  {
209  ln->printHoughLine();
210  }
211 
212  // Determine intersection between outer lines
213  HoughLine* oneLine = top3Lines[0];
214  HoughLine* otherLine = top3Lines[2];
215  QPointF intersection;
216  HoughLine::IntersectResult result = oneLine->Intersect(*otherLine, intersection);
217  if (result == HoughLine::INTERESECTING)
218  {
219 
220  qCDebug(KSTARS_FITS) << "Intersection: " << intersection.x() << ", " << intersection.y();
221 
222  // Determine offset between intersection and middle line
223  HoughLine* midLine = top3Lines[1];
224  QPointF intersectionOnMidLine;
225  double distance;
226  if (midLine->DistancePointLine(intersection, intersectionOnMidLine, distance))
227  {
228  qCDebug(KSTARS_FITS) << "Distance between intersection and midline is " << distance
229  << " at mid line point " << intersectionOnMidLine.x() << ", "
230  << intersectionOnMidLine.y();
231 
232  // Add star center to selected stars
233  // Maximum Radius
234  int maxR = qMin(subW - 1, subH - 1) / 2;
235  BahtinovEdge* center = new BahtinovEdge();
236  center->width = maxR / 3;
237  center->x = subX + intersection.x();
238  center->y = subY + intersection.y();
239  // Set distance value in HFR
240  center->HFR = distance;
241 
242  center->offset.setX(subX + intersectionOnMidLine.x());
243  center->offset.setY(subY + intersectionOnMidLine.y());
244  oneLine->Offset(subX, subY);
245  midLine->Offset(subX, subY);
246  otherLine->Offset(subX, subY);
247  center->line.append(*oneLine);
248  center->line.append(*midLine);
249  center->line.append(*otherLine);
250  starCenters.append(center);
251  }
252  else
253  {
254  qCWarning(KSTARS_FITS) << "Closest point does not fall within the line segment.";
255  }
256  }
257  else
258  {
259  qCWarning(KSTARS_FITS) << "Lines are not intersecting (result: " << result << ")";
260  }
261  }
262 
263  // Clean up Bahtinov line array (of pointers) as they are no longer needed
264  for (int index = 0; index < bahtinov_angles.size(); index++)
265  {
266  HoughLine* pLineAverage = bahtinov_angles[index];
267  if (pLineAverage != nullptr)
268  {
269  delete pLineAverage;
270  }
271  }
272  bahtinov_angles.clear();
273 
274  // Clean up line averages array as they are no longer needed
275  lineAveragesPerAngle.clear();
276 
277  top3Lines.clear();
278 
279  m_ImageData->setStarCenters(starCenters);
280 
281  return true;
282 }
283 
284 template <typename T>
285 BahtinovLineAverage FITSBahtinovDetector::calculateMaxAverage(const FITSData *data, int angle)
286 {
287  int BBP = data->getBytesPerPixel();
288  int size = data->getStatistics().samples_per_channel;
289  int width = data->width();
290  int height = data->height();
291  int numChannels = data->channels();
292 
293  BahtinovLineAverage lineAverage;
294  auto * rotimage = new T[size * BBP];
295 
296  rotateImage(data, angle, rotimage);
297 
298  // Calculate average pixel value for each row
299  auto * rotBuffer = reinterpret_cast<T *>(rotimage);
300 
301  // printf("Angle;%d;Width;%d;Height;%d;Rows;%d;;RowSum;", angle, width, height, NUMBER_OF_AVERAGE_ROWS);
302 
303  int NUMBER_OF_AVERAGE_ROWS = getValue("NUMBER_OF_AVERAGE_ROWS", 1).toInt();
304  if (NUMBER_OF_AVERAGE_ROWS % 2 == 0)
305  {
306  NUMBER_OF_AVERAGE_ROWS--;
307  qCWarning(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to "
308  << NUMBER_OF_AVERAGE_ROWS;
309  }
310  // Rows must be a positive number!
311  if (NUMBER_OF_AVERAGE_ROWS < 1)
312  {
313  NUMBER_OF_AVERAGE_ROWS = 1;
314  qCWarning(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to "
315  << NUMBER_OF_AVERAGE_ROWS;
316  }
317 
318  for (int y = 0; y < height; y++)
319  {
320  int yMin = y - ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
321  int yMax = y + ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
322 
323  unsigned long multiRowSum = 0;
324  // Calculate average over multiple rows
325  for (int y1 = yMin; y1 <= yMax; y1++)
326  {
327  int y2 = y1;
328  if (y2 < 0)
329  {
330  y2 += height;
331  }
332  if (y2 >= height)
333  {
334  y2 -= height;
335  }
336  if (y2 < 0 || y2 >= height)
337  {
338  qCWarning(KSTARS_FITS) << "Y still out of bounds: 0 <=" << y2 << "<" << height;
339  }
340 
341  for (int x1 = 0; x1 < width; x1++)
342  {
343  int index = y2 * width + x1;
344  unsigned long channelAverage = 0;
345  for (int i = 0; i < numChannels; i++)
346  {
347  int offset = size * i;
348  channelAverage += rotBuffer[index + offset];
349  }
350  multiRowSum += qRound(channelAverage / static_cast<double>(numChannels));
351  }
352  }
353  // printf("%lu;", multiRowSum);
354 
355  double average = multiRowSum / static_cast<double>(width * NUMBER_OF_AVERAGE_ROWS);
356  if (average > lineAverage.average)
357  {
358  lineAverage.average = average;
359  lineAverage.offset = y;
360  }
361  }
362  // printf(";;MaxAverage;%.3f;MaxAverageIndex;%lu\r\n", lineAverage.average, lineAverage.offset);
363  // fflush(stdout);
364 
365  rotBuffer = nullptr;
366  delete[] rotimage;
367  rotimage = nullptr;
368 
369  return lineAverage;
370 }
371 
372 /** Rotate an image by angle degrees.
373  * verbose generates extra info on stdout.
374  * return true if successful and rotated image.
375  * @param angle The angle over which the image needs to be rotated
376  * @param rotImage The image that needs to be rotated, also the rotated image return value
377  */
378 template <typename T>
379 bool FITSBahtinovDetector::rotateImage(const FITSData *data, int angle, T * rotImage)
380 {
381  int BBP = data->getBytesPerPixel();
382  int size = data->getStatistics().samples_per_channel;
383  int width = data->width();
384  int height = data->height();
385  int numChannels = data->channels();
386 
387  int hx, hy;
388  size_t bufferSize;
389 
390  /* Check allocation buffer for rotated image */
391  if (rotImage == nullptr)
392  {
393  qWarning() << "No memory allocated for rotated image buffer!";
394  return false;
395  }
396 
397  while (angle < 0)
398  {
399  angle = angle + 360;
400  }
401  while (angle >= 360)
402  {
403  angle = angle - 360;
404  }
405 
406  hx = qFloor((width + 1) / 2.0);
407  hy = qFloor((height + 1) / 2.0);
408 
409  bufferSize = size * numChannels * BBP;
410  memset(rotImage, 0, bufferSize);
411 
412  auto * rotBuffer = reinterpret_cast<T *>(rotImage);
413  auto * buffer = reinterpret_cast<const T *>(data->getImageBuffer());
414 
415  double innerCircleRadius = (0.5 * qSqrt(2.0) * qMin(hx, hy));
416  double angleInRad = angle * M_PI / 180.0;
417  double sinAngle = qSin(angleInRad);
418  double cosAngle = qCos(angleInRad);
419  int leftEdge = qCeil(hx - innerCircleRadius);
420  int rightEdge = qFloor(hx + innerCircleRadius);
421  int topEdge = qCeil(hy - innerCircleRadius);
422  int bottomEdge = qFloor(hy + innerCircleRadius);
423 
424  for (int i = 0; i < numChannels; i++)
425  {
426  int offset = size * i;
427  for (int x1 = leftEdge; x1 < rightEdge; x1++)
428  {
429  for (int y1 = topEdge; y1 < bottomEdge; y1++)
430  {
431  // translate point back to origin:
432  double x2 = x1 - hx;
433  double y2 = y1 - hy;
434 
435  // rotate point
436  double xnew = x2 * cosAngle - y2 * sinAngle;
437  double ynew = x2 * sinAngle + y2 * cosAngle;
438 
439  // translate point back:
440  x2 = xnew + hx;
441  y2 = ynew + hy;
442 
443  int orgIndex = y1 * height + x1;
444  int newIndex = qRound(y2) * height + qRound(x2);
445 
446  if (newIndex >= 0 && newIndex < static_cast<int>(bufferSize))
447  {
448  rotBuffer[newIndex + offset] = buffer[orgIndex + offset];
449  } // else index out of bounds, do not update pixel
450  }
451  }
452  }
453 
454  return true;
455 }
void append(const T &value)
QFuture< T > run(Function function,...)
void clear()
void append(const T &value)
int width() const const
int x() const const
int y() const const
int remove(const Key &key)
QMap::iterator insert(const Key &key, const T &value)
void clear()
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
bool isNull() const const
qint64 elapsed() const const
int height() const const
qreal x() const const
qreal y() const const
bool isEmpty() const const
QTextStream & center(QTextStream &stream)
int size() const const
IntersectResult Intersect(const HoughLine &other_line, QPointF &intersection)
Sources for intersection and distance calculations came from http://paulbourke.net/geometry/pointline...
Definition: houghline.cpp:94
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.