14#include <QtConcurrent>
16#include "fits_debug.h"
17#include "fitsgradientdetector.h"
22 FITSImage::Statistic
const &stats = m_ImageData->getStatistics();
23 switch (stats.dataType)
28 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<uint8_t>, boundary);
31 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<int16_t>, boundary);
34 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<uint16_t>, boundary);
37 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<int32_t>, boundary);
40 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<uint16_t>, boundary);
43 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<float>, boundary);
46 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<int64_t>, boundary);
49 return QtConcurrent::run(
this, &FITSGradientDetector::findSources<double>, boundary);
54bool FITSGradientDetector::findSources(
const QRect &boundary)
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());
61 int BBP = m_ImageData->getBytesPerPixel();
63 uint16_t dataWidth = m_ImageData->width();
66 uint32_t size = subW * subH;
67 uint32_t offset = subX + subY * dataWidth;
70 auto * buffer =
new uint8_t[size * BBP];
73 memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP);
76 uint8_t * dataPtr = buffer;
77 uint8_t
const * origDataPtr = m_ImageData->getImageBuffer();
78 uint32_t lineOffset = 0;
80 for (
int height = subY; height < (subY + subH); height++)
82 lineOffset = (subX + height * dataWidth) * BBP;
83 memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
84 dataPtr += (subW * BBP);
89 auto * boundedImage =
new FITSData();
90 FITSImage::Statistic stats;
93 stats.dataType = m_ImageData->getStatistics().dataType;
94 stats.bytesPerPixel = m_ImageData->getStatistics().bytesPerPixel;
95 stats.samples_per_channel = size;
97 boundedImage->restoreStatistics(stats);
100 boundedImage->setImageBuffer(buffer);
102 boundedImage->calculateStats(
true);
105 boundedImage->applyFilter(FITS_MEDIAN);
106 boundedImage->applyFilter(FITS_HIGH_CONTRAST);
114 sobel<T>(boundedImage, gradients, directions);
118 int maxID = partition(subW, subH, gradients, ids);
126 typedef struct massInfo
136 for (
int y = 0; y < subH; y++)
138 for (
int x = 0; x < subW; x++)
140 int index = x + y * subW;
142 int regionID = ids[index];
145 float pixel = gradients[index];
147 masses[regionID].totalMass += pixel;
148 masses[regionID].massX += x * pixel;
149 masses[regionID].massY += y * pixel;
156 int maxTotalMass = masses[1].totalMass;
157 double totalMassRatio = 1e6;
158 for (
auto key : masses.keys())
160 massInfo oneMass = masses.
value(key);
161 if (oneMass.totalMass > maxTotalMass)
163 totalMassRatio = oneMass.totalMass / maxTotalMass;
164 maxTotalMass = oneMass.totalMass;
171 if (maxID > 10 && totalMassRatio < 1.5)
176 center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
177 center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
181 int maxR = qMin(subW - 1, subH - 1) / 2;
183 for (
int r = maxR; r > 1; r--)
187 for (
float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
189 int testX =
center->x + std::cos(theta) * r;
190 int testY =
center->y + std::sin(theta) * r;
193 if (testX < 0 || testX >= subW || testY < 0 || testY >= subH)
196 if (gradients[testX + testY * subW] > 0)
210 qCDebug(KSTARS_FITS) <<
"FITS: Weighted Center is X: " <<
center->x <<
" Y: " <<
center->y <<
" Width: " <<
center->width;
222 double FSum = 0, HF = 0, TF = 0;
223 const double resolution = 1.0 / 20.0;
225 int cen_y = qRound(
center->y);
233 const T * origBuffer =
reinterpret_cast<T
const *
>(m_ImageData->getImageBuffer()) + offset;
235 for (
double x = leftEdge; x <= rightEdge; x += resolution)
237 double slice = resolution * (origBuffer[
static_cast<int>(floor(x)) + cen_y * dataWidth]);
245 int subPixelCenter = (
center->width / resolution) / 2;
248 TF = subPixels[subPixelCenter];
252 for (
int k = 1; k < subPixelCenter; k++)
254 TF += subPixels[subPixelCenter + k];
255 TF += subPixels[subPixelCenter - k];
262 center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
277 starCenters.
append(center);
278 m_ImageData->setStarCenters(starCenters);
280 qCDebug(KSTARS_FITS) <<
"Flux: " << FSum <<
" Half-Flux: " << HF <<
" HFR: " <<
center->HFR;
295 FITSImage::Statistic
const &stats = data->getStatistics();
298 gradient.
resize(stats.samples_per_channel);
299 direction.
resize(stats.samples_per_channel);
301 for (
int y = 0; y < stats.height; y++)
303 size_t yOffset = y * stats.width;
304 const T * grayLine =
reinterpret_cast<T
const *
>(data->getImageBuffer()) + yOffset;
306 const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width;
307 const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width;
309 float * gradientLine = gradient.
data() + yOffset;
310 float * directionLine = direction.
data() + yOffset;
312 for (
int x = 0; x < stats.width; x++)
314 int x_m1 = x < 1 ? x : x - 1;
315 int x_p1 = x >= stats.width - 1 ? x : x + 1;
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];
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];
323 gradientLine[x] = qAbs(gradX) + qAbs(gradY);
351 if (gradX == 0 && gradY == 0)
352 directionLine[x] = 0;
354 directionLine[x] = 3;
357 qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
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;
366 directionLine[x] = 3;
376 for (
int y = 1; y < height - 1; y++)
378 for (
int x = 1; x < width - 1; x++)
380 int index = x + y * width;
381 float val = gradient[index];
382 if (val > 0 && ids[index] == 0)
384 trace(width, height, ++
id, gradient, ids, x, y);
396 int yOffset = y * width;
397 float * cannyLine = image.
data() + yOffset;
398 int * idLine = ids.
data() + yOffset;
405 for (
int j = -1; j < 2; j++)
409 if (nextY < 0 || nextY >= height)
412 float * cannyLineNext = cannyLine + j * width;
414 for (
int i = -1; i < 2; i++)
418 if (i == j || nextX < 0 || nextX >= width)
421 if (cannyLineNext[nextX] > 0)
424 trace(width, height,
id, image, ids, nextX, nextY);
void append(QList< T > &&value)
void reserve(qsizetype size)
void resize(qsizetype size)
qsizetype size() const const
T value(const Key &key, const T &defaultValue) const const
bool isNull() const const
QTextStream & center(QTextStream &stream)
QFuture< T > run(Function function,...)