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 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
80template <typename T>
81bool 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;
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
136
137 // Rotate image 180 degrees in steps of 1 degree
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
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;
173 maxAngle = angle;
174 }
175 }
177 if (pHoughLine != nullptr)
178 {
180 }
181
182 // Remove data around peak to prevent it from being detected again
183 int minBahtinovAngleOffset = 18;
185 {
187 if (angleInRange < 0)
188 {
189 angleInRange += 180;
190 }
191 if (angleInRange >= 180)
192 {
193 angleInRange -= 180;
194 }
196 }
197 }
198
199 // Proceed with focus offset calculation, but only when at least 3 lines have been detected
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
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
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() << ", "
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 {
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
284template <typename T>
285BahtinovLineAverage 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 {
307 qCWarning(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to "
309 }
310 // Rows must be a positive number!
312 {
314 qCWarning(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to "
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 */
378template <typename T>
379bool 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);
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}
Line representation for HoughTransform Based on the java implementation found on http://vase....
Definition houghline.h:23
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
qreal x() const const
qreal y() const const
int height() const const
bool isEmpty() const const
bool isNull() const const
int width() const const
QTextStream & center(QTextStream &stream)
QFuture< T > run(Function function,...)
int toInt(bool *ok) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri May 24 2024 11:49:22 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.