Kstars

fitsgradientdetector.cpp
1/*
2 SPDX-FileCopyrightText: 2004 Jasem Mutlaq
3 SPDX-FileCopyrightText: 2020 Eric Dejouhanet <eric.dejouhanet@gmail.com>
4 SPDX-FileCopyrightText: 2015 Gonzalo Exequiel Pedone <hipersayan.x@gmail.com>
5
6 SPDX-License-Identifier: GPL-2.0-or-later AND GPL-3.0-or-later
7
8 Some code fragments were adapted from Peter Kirchgessner's FITS plugin
9 SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
10*/
11
12#include <math.h>
13#include <cmath>
14#include <QtConcurrent>
15
16#include "fits_debug.h"
17#include "fitsgradientdetector.h"
18#include "fitsdata.h"
19
20QFuture<bool> FITSGradientDetector::findSources(const QRect &boundary)
21{
22 FITSImage::Statistic const &stats = m_ImageData->getStatistics();
23 switch (stats.dataType)
24 {
25
26 case TBYTE:
27 default:
28 return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint8_t>, boundary);
29
30 case TSHORT:
31 return QtConcurrent::run(this, &FITSGradientDetector::findSources<int16_t>, boundary);
32
33 case TUSHORT:
34 return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint16_t>, boundary);
35
36 case TLONG:
37 return QtConcurrent::run(this, &FITSGradientDetector::findSources<int32_t>, boundary);
38
39 case TULONG:
40 return QtConcurrent::run(this, &FITSGradientDetector::findSources<uint16_t>, boundary);
41
42 case TFLOAT:
43 return QtConcurrent::run(this, &FITSGradientDetector::findSources<float>, boundary);
44
45 case TLONGLONG:
46 return QtConcurrent::run(this, &FITSGradientDetector::findSources<int64_t>, boundary);
47
48 case TDOUBLE:
49 return QtConcurrent::run(this, &FITSGradientDetector::findSources<double>, boundary);
50 }
51}
52
53template <typename T>
54bool FITSGradientDetector::findSources(const QRect &boundary)
55{
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());
60
61 int BBP = m_ImageData->getBytesPerPixel();
62
63 uint16_t dataWidth = m_ImageData->width();
64
65 // #1 Find offsets
66 uint32_t size = subW * subH;
67 uint32_t offset = subX + subY * dataWidth;
68
69 // #2 Create new buffer
70 auto * buffer = new uint8_t[size * BBP];
71 // If there is no offset, copy whole buffer in one go
72 if (offset == 0)
73 memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP);
74 else
75 {
76 uint8_t * dataPtr = buffer;
77 uint8_t const * origDataPtr = m_ImageData->getImageBuffer();
78 uint32_t lineOffset = 0;
79 // Copy data line by line
80 for (int height = subY; height < (subY + subH); height++)
81 {
82 lineOffset = (subX + height * dataWidth) * BBP;
84 dataPtr += (subW * BBP);
85 }
86 }
87
88 // #3 Create new FITSData to hold it
89 auto * boundedImage = new FITSData();
90 FITSImage::Statistic stats;
91 stats.width = subW;
92 stats.height = subH;
93 stats.dataType = m_ImageData->getStatistics().dataType;
94 stats.bytesPerPixel = m_ImageData->getStatistics().bytesPerPixel;
95 stats.samples_per_channel = size;
96 stats.ndim = 2;
97 boundedImage->restoreStatistics(stats);
98
99 // #4 Set image buffer and calculate stats.
100 boundedImage->setImageBuffer(buffer);
101
102 boundedImage->calculateStats(true);
103
104 // #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain
105 boundedImage->applyFilter(FITS_MEDIAN);
106 boundedImage->applyFilter(FITS_HIGH_CONTRAST);
107
108 // #6 Perform Sobel to find gradients and their directions
111
112 // TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed
113 // and discarded whenever necessary. It won't work on noisy images unless this is done.
115
116 QVector<int> ids(gradients.size());
117
118 int maxID = partition(subW, subH, gradients, ids);
119
120 // Not needed anymore
121 delete boundedImage;
122
123 if (maxID == 0)
124 return 0;
125
126 typedef struct massInfo
127 {
128 float massX = 0;
129 float massY = 0;
130 float totalMass = 0;
131 } massInfo;
132
134
135 // #7 Calculate center of mass for all detected regions
136 for (int y = 0; y < subH; y++)
137 {
138 for (int x = 0; x < subW; x++)
139 {
140 int index = x + y * subW;
141
142 int regionID = ids[index];
143 if (regionID > 0)
144 {
145 float pixel = gradients[index];
146
147 masses[regionID].totalMass += pixel;
148 masses[regionID].massX += x * pixel;
149 masses[regionID].massY += y * pixel;
150 }
151 }
152 }
153
154 // Compare multiple masses, and only select the highest total mass one as the desired star
155 int maxRegionID = 1;
156 int maxTotalMass = masses[1].totalMass;
157 double totalMassRatio = 1e6;
158 for (auto key : masses.keys())
159 {
160 massInfo oneMass = masses.value(key);
161 if (oneMass.totalMass > maxTotalMass)
162 {
163 totalMassRatio = oneMass.totalMass / maxTotalMass;
164 maxTotalMass = oneMass.totalMass;
165 maxRegionID = key;
166 }
167 }
168
169 // If image has many regions and there is no significant relative center of mass then it's just noise and no stars
170 // are probably there above a useful threshold.
171 if (maxID > 10 && totalMassRatio < 1.5)
172 return false;
173
174 auto * center = new Edge;
175 center->width = -1;
176 center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
177 center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
178 center->HFR = 1;
179
180 // Maximum Radius
181 int maxR = qMin(subW - 1, subH - 1) / 2;
182
183 for (int r = maxR; r > 1; r--)
184 {
185 int pass = 0;
186
187 for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
188 {
189 int testX = center->x + std::cos(theta) * r;
190 int testY = center->y + std::sin(theta) * r;
191
192 // if out of bound, break;
194 break;
195
196 if (gradients[testX + testY * subW] > 0)
197 //if (thresholded[testX + testY * subW] > 0)
198 {
199 if (++pass >= 24)
200 {
201 center->width = r * 2;
202 // Break of outer loop
203 r = 0;
204 break;
205 }
206 }
207 }
208 }
209
210 qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width;
211
212 // If no stars were detected
213 if (center->width == -1)
214 {
215 delete center;
216 return false;
217 }
218
219 // 30% fuzzy
220 //center->width += center->width*0.3 * (running_threshold / threshold);
221
222 double FSum = 0, HF = 0, TF = 0;
223 const double resolution = 1.0 / 20.0;
224
225 int cen_y = qRound(center->y);
226
227 double rightEdge = center->x + center->width / 2.0;
228 double leftEdge = center->x - center->width / 2.0;
229
231 subPixels.reserve(center->width / resolution);
232
233 const T * origBuffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer()) + offset;
234
235 for (double x = leftEdge; x <= rightEdge; x += resolution)
236 {
237 double slice = resolution * (origBuffer[static_cast<int>(floor(x)) + cen_y * dataWidth]);
238 FSum += slice;
239 subPixels.append(slice);
240 }
241
242 // Half flux
243 HF = FSum / 2.0;
244
245 int subPixelCenter = (center->width / resolution) / 2;
246
247 // Start from center
249 double lastTF = TF;
250 // Integrate flux along radius axis until we reach half flux
251 //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
252 for (int k = 1; k < subPixelCenter; k++)
253 {
256
257 if (TF >= HF)
258 {
259 // We overpassed HF, let's calculate from last TF how much until we reach HF
260
261 // #1 Accurate calculation, but very sensitive to small variations of flux
262 center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
263
264 // #2 Less accurate calculation, but stable against small variations of flux
265 //center->HFR = (k - 1) * resolution;
266 break;
267 }
268
269 lastTF = TF;
270 }
271
272 // Correct center for subX and subY
273 center->x += subX;
274 center->y += subY;
275
276 QList<Edge*> starCenters;
277 starCenters.append(center);
278 m_ImageData->setStarCenters(starCenters);
279
280 qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR;
281
282 return true;
283}
284
285/* CannyDetector, Implementation of Canny edge detector in Qt/C++.
286 * Web-Site: https://github.com/hipersayanX/CannyDetector
287 */
288
289template <typename T>
290void FITSGradientDetector::sobel(FITSData const *data, QVector<float> &gradient, QVector<float> &direction) const
291{
292 if (data == nullptr)
293 return;
294
295 FITSImage::Statistic const &stats = data->getStatistics();
296
297 //int size = image.width() * image.height();
298 gradient.resize(stats.samples_per_channel);
299 direction.resize(stats.samples_per_channel);
300
301 for (int y = 0; y < stats.height; y++)
302 {
303 size_t yOffset = y * stats.width;
304 const T * grayLine = reinterpret_cast<T const *>(data->getImageBuffer()) + yOffset;
305
306 const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width;
307 const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width;
308
309 float * gradientLine = gradient.data() + yOffset;
310 float * directionLine = direction.data() + yOffset;
311
312 for (int x = 0; x < stats.width; x++)
313 {
314 int x_m1 = x < 1 ? x : x - 1;
315 int x_p1 = x >= stats.width - 1 ? x : x + 1;
316
319
321 2 * grayLine_p1[x] - grayLine_p1[x_p1];
322
324
325 /* Gradient directions are classified in 4 possible cases
326 *
327 * dir 0
328 *
329 * x x x
330 * - - -
331 * x x x
332 *
333 * dir 1
334 *
335 * x x /
336 * x / x
337 * / x x
338 *
339 * dir 2
340 *
341 * \ x x
342 * x \ x
343 * x x \
344 *
345 * dir 3
346 *
347 * x | x
348 * x | x
349 * x | x
350 */
351 if (gradX == 0 && gradY == 0)
352 directionLine[x] = 0;
353 else if (gradX == 0)
354 directionLine[x] = 3;
355 else
356 {
357 qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
358
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;
365 else
366 directionLine[x] = 3;
367 }
368 }
369 }
370}
371
372int FITSGradientDetector::partition(int width, int height, QVector<float> &gradient, QVector<int> &ids) const
373{
374 int id = 0;
375
376 for (int y = 1; y < height - 1; y++)
377 {
378 for (int x = 1; x < width - 1; x++)
379 {
380 int index = x + y * width;
381 float val = gradient[index];
382 if (val > 0 && ids[index] == 0)
383 {
384 trace(width, height, ++id, gradient, ids, x, y);
385 }
386 }
387 }
388
389 // Return max id
390 return id;
391}
392
393void FITSGradientDetector::trace(int width, int height, int id, QVector<float> &image, QVector<int> &ids, int x,
394 int y) const
395{
396 int yOffset = y * width;
397 float * cannyLine = image.data() + yOffset;
398 int * idLine = ids.data() + yOffset;
399
400 if (idLine[x] != 0)
401 return;
402
403 idLine[x] = id;
404
405 for (int j = -1; j < 2; j++)
406 {
407 int nextY = y + j;
408
409 if (nextY < 0 || nextY >= height)
410 continue;
411
412 float * cannyLineNext = cannyLine + j * width;
413
414 for (int i = -1; i < 2; i++)
415 {
416 int nextX = x + i;
417
418 if (i == j || nextX < 0 || nextX >= width)
419 continue;
420
421 if (cannyLineNext[nextX] > 0)
422 {
423 // Trace neighbors.
424 trace(width, height, id, image, ids, nextX, nextY);
425 }
426 }
427 }
428}
int height() const const
bool isNull() const const
int width() const const
QTextStream & center(QTextStream &stream)
QFuture< T > run(Function function,...)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Tue Mar 26 2024 11:19:03 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.