Kstars

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

KDE's Doxygen guidelines are available online.