13#include "fits_debug.h"
14#include "fitsbahtinovdetector.h"
15#include "hough/houghline.h"
18#include <QElapsedTimer>
19#include <QtConcurrent>
49 FITSImage::Statistic
const &stats = m_ImageData->getStatistics();
50 switch (stats.dataType)
53 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<int16_t>, boundary);
56 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<uint16_t>, boundary);
59 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<int32_t>, boundary);
62 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<uint32_t>, boundary);
65 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<float>, boundary);
68 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<int64_t>, boundary);
71 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<double>, boundary);
75 return QtConcurrent::run(
this, &FITSBahtinovDetector::findBahtinovStar<uint8_t>, boundary);
81bool FITSBahtinovDetector::findBahtinovStar(
const QRect &boundary)
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());
92 int BBP = m_ImageData->getBytesPerPixel();
93 uint16_t dataWidth = m_ImageData->width();
96 uint32_t size = subW * subH;
97 uint32_t offset = subX + subY * dataWidth;
100 auto * buffer =
new uint8_t[size * BBP];
104 memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP);
108 uint8_t * dataPtr = buffer;
109 const uint8_t * origDataPtr = m_ImageData->getImageBuffer();
111 for (
int height = subY; height < (subY + subH); height++)
113 uint32_t lineOffset = (subX + height * dataWidth) * BBP;
114 memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
115 dataPtr += (subW * BBP);
120 FITSImage::Statistic stats;
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);
131 boundedImage->setImageBuffer(buffer);
133 boundedImage->calculateStats(
true);
139 const int steps = 180;
140 double radPerStep = M_PI / steps;
144 for (
int angle = 0; angle < steps; angle++)
147 BahtinovLineAverage lineAverage = calculateMaxAverage<T>(boundedImage, angle);
149 lineAveragesPerAngle.
insert(angle, lineAverage);
152 qCDebug(KSTARS_FITS) <<
"Getting max average for all 180 rotations took" << timer1.
elapsed() <<
"milliseconds";
161 for (
int index1 = 0; index1 < 3; index1++)
163 double maxAverage = 0.0;
164 double maxAngle = 0.0;
165 int maxAverageOffset = 0;
166 for (
int angle = 0; angle < steps; angle++)
168 BahtinovLineAverage lineAverage = lineAveragesPerAngle[angle];
169 if (lineAverage.average > maxAverage)
171 maxAverage = lineAverage.average;
172 maxAverageOffset = lineAverage.offset;
176 HoughLine* pHoughLine =
new HoughLine(maxAngle * radPerStep, maxAverageOffset, subW, subH, maxAverage);
177 if (pHoughLine !=
nullptr)
179 bahtinov_angles.
append(pHoughLine);
183 int minBahtinovAngleOffset = 18;
184 for (
int subAngle = maxAngle - minBahtinovAngleOffset; subAngle < maxAngle + minBahtinovAngleOffset; subAngle++)
186 int angleInRange = subAngle;
187 if (angleInRange < 0)
191 if (angleInRange >= 180)
195 lineAveragesPerAngle.
remove(angleInRange);
201 if (bahtinov_angles.
size() >= 3)
203 HoughLine::getSortedTopThreeLines(bahtinov_angles, top3Lines);
206 qCDebug(KSTARS_FITS) <<
"Sorted bahtinov angles:";
209 ln->printHoughLine();
216 HoughLine::IntersectResult result = oneLine->
Intersect(*otherLine, intersection);
217 if (result == HoughLine::INTERESECTING)
220 qCDebug(KSTARS_FITS) <<
"Intersection: " << intersection.
x() <<
", " << intersection.
y();
226 if (midLine->DistancePointLine(intersection, intersectionOnMidLine, distance))
228 qCDebug(KSTARS_FITS) <<
"Distance between intersection and midline is " <<
distance
229 <<
" at mid line point " << intersectionOnMidLine.
x() <<
", "
230 << intersectionOnMidLine.
y();
234 int maxR = qMin(subW - 1, subH - 1) / 2;
235 BahtinovEdge*
center =
new BahtinovEdge();
237 center->x = subX + intersection.
x();
238 center->y = subY + intersection.
y();
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);
254 qCWarning(KSTARS_FITS) <<
"Closest point does not fall within the line segment.";
259 qCWarning(KSTARS_FITS) <<
"Lines are not intersecting (result: " << result <<
")";
264 for (
int index = 0; index < bahtinov_angles.
size(); index++)
266 HoughLine* pLineAverage = bahtinov_angles[index];
267 if (pLineAverage !=
nullptr)
272 bahtinov_angles.
clear();
275 lineAveragesPerAngle.
clear();
279 m_ImageData->setStarCenters(starCenters);
285BahtinovLineAverage FITSBahtinovDetector::calculateMaxAverage(
const FITSData *data,
int angle)
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();
293 BahtinovLineAverage lineAverage;
294 auto * rotimage =
new T[size * BBP];
296 rotateImage(data, angle, rotimage);
299 auto * rotBuffer =
reinterpret_cast<T *
>(rotimage);
303 int NUMBER_OF_AVERAGE_ROWS = getValue(
"NUMBER_OF_AVERAGE_ROWS", 1).
toInt();
304 if (NUMBER_OF_AVERAGE_ROWS % 2 == 0)
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;
311 if (NUMBER_OF_AVERAGE_ROWS < 1)
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;
318 for (
int y = 0; y < height; y++)
320 int yMin = y - ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
321 int yMax = y + ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
323 unsigned long multiRowSum = 0;
325 for (
int y1 = yMin; y1 <= yMax; y1++)
336 if (y2 < 0 || y2 >= height)
338 qCWarning(KSTARS_FITS) <<
"Y still out of bounds: 0 <=" << y2 <<
"<" << height;
341 for (
int x1 = 0; x1 < width; x1++)
343 int index = y2 * width + x1;
344 unsigned long channelAverage = 0;
345 for (
int i = 0; i < numChannels; i++)
347 int offset = size * i;
348 channelAverage += rotBuffer[index + offset];
350 multiRowSum += qRound(channelAverage /
static_cast<double>(numChannels));
355 double average = multiRowSum /
static_cast<double>(width * NUMBER_OF_AVERAGE_ROWS);
356 if (average > lineAverage.average)
358 lineAverage.average = average;
359 lineAverage.offset = y;
379bool FITSBahtinovDetector::rotateImage(
const FITSData *data,
int angle, T * rotImage)
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();
391 if (rotImage ==
nullptr)
393 qWarning() <<
"No memory allocated for rotated image buffer!";
406 hx = qFloor((width + 1) / 2.0);
407 hy = qFloor((height + 1) / 2.0);
409 bufferSize = size * numChannels * BBP;
410 memset(rotImage, 0, bufferSize);
412 auto * rotBuffer =
reinterpret_cast<T *
>(rotImage);
413 auto * buffer =
reinterpret_cast<const T *
>(data->getImageBuffer());
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);
424 for (
int i = 0; i < numChannels; i++)
426 int offset = size * i;
427 for (
int x1 = leftEdge; x1 < rightEdge; x1++)
429 for (
int y1 = topEdge; y1 < bottomEdge; y1++)
436 double xnew = x2 * cosAngle - y2 * sinAngle;
437 double ynew = x2 * sinAngle + y2 * cosAngle;
443 int orgIndex = y1 * height + x1;
444 int newIndex = qRound(y2) * height + qRound(x2);
446 if (newIndex >= 0 && newIndex <
static_cast<int>(bufferSize))
448 rotBuffer[newIndex + offset] = buffer[orgIndex + offset];
Line representation for HoughTransform Based on the java implementation found on http://vase....
IntersectResult Intersect(const HoughLine &other_line, QPointF &intersection)
Sources for intersection and distance calculations came from http://paulbourke.net/geometry/pointline...
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
qint64 elapsed() const const
void append(QList< T > &&value)
qsizetype size() const const
iterator insert(const Key &key, const T &value)
size_type remove(const Key &key)
bool setProperty(const char *name, QVariant &&value)
bool isEmpty() const const
bool isNull() const const
QTextStream & center(QTextStream &stream)
QFuture< T > run(Function function,...)
int toInt(bool *ok) const const