Kstars

guidealgorithms.cpp
1/* Algorithms used in the internal guider.
2 Copyright (C) 2012 Andrew Stepanenko
3
4 Modified by Jasem Mutlaq for KStars.
5 SPDX-FileCopyrightText: Jasem Mutlaq <mutlaqja@ikarustech.com>
6
7 SPDX-License-Identifier: GPL-2.0-or-later
8 */
9
10#include "guidealgorithms.h"
11
12#include <set>
13#include <QObject>
14
15#include "ekos_guide_debug.h"
16#include "ekos/auxiliary/stellarsolverprofileeditor.h"
17#include "Options.h"
18#include "fitsviewer/fitsdata.h"
19
20#define SMART_THRESHOLD 0
21#define SEP_THRESHOLD 1
22#define CENTROID_THRESHOLD 2
23#define AUTO_THRESHOLD 3
24#define NO_THRESHOLD 4
25#define SEP_MULTISTAR 5
26
27// smart threshold algorithm param
28// width of outer frame for background calculation
29#define SMART_FRAME_WIDTH 4
30// cut-factor above average threshold
31#define SMART_CUT_FACTOR 0.1
32
33struct Peak
34{
35 int x;
36 int y;
37 float val;
38
39 Peak() = default;
40 Peak(int x_, int y_, float val_) : x(x_), y(y_), val(val_) { }
41 bool operator<(const Peak &rhs) const
42 {
43 return val < rhs.val;
44 }
45};
46
47// JM: Why not use QPoint?
48typedef struct
49{
50 int x, y;
51} point_t;
52
53namespace
54{
55
56void psf_conv(float *dst, const float *src, int width, int height)
57{
58 //dst.Init(src.Size);
59
60 // A B1 B2 C1 C2 C3 D1 D2 D3
61 const double PSF[] = { 0.906, 0.584, 0.365, .117, .049, -0.05, -.064, -.074, -.094 };
62
63 //memset(dst.px, 0, src.NPixels * sizeof(float));
64
65 /* PSF Grid is:
66 D3 D3 D3 D3 D3 D3 D3 D3 D3
67 D3 D3 D3 D2 D1 D2 D3 D3 D3
68 D3 D3 C3 C2 C1 C2 C3 D3 D3
69 D3 D2 C2 B2 B1 B2 C2 D2 D3
70 D3 D1 C1 B1 A B1 C1 D1 D3
71 D3 D2 C2 B2 B1 B2 C2 D2 D3
72 D3 D3 C3 C2 C1 C2 C3 D3 D3
73 D3 D3 D3 D2 D1 D2 D3 D3 D3
74 D3 D3 D3 D3 D3 D3 D3 D3 D3
75
76 1@A
77 4@B1, B2, C1, C3, D1
78 8@C2, D2
79 44 * D3
80 */
81
82 int psf_size = 4;
83
84 for (int y = psf_size; y < height - psf_size; y++)
85 {
86 for (int x = psf_size; x < width - psf_size; x++)
87 {
88 float A, B1, B2, C1, C2, C3, D1, D2, D3;
89
90#define PX(dx, dy) *(src + width * (y + (dy)) + x + (dx))
91 A = PX(+0, +0);
92 B1 = PX(+0, -1) + PX(+0, +1) + PX(+1, +0) + PX(-1, +0);
93 B2 = PX(-1, -1) + PX(+1, -1) + PX(-1, +1) + PX(+1, +1);
94 C1 = PX(+0, -2) + PX(-2, +0) + PX(+2, +0) + PX(+0, +2);
95 C2 = PX(-1, -2) + PX(+1, -2) + PX(-2, -1) + PX(+2, -1) + PX(-2, +1) + PX(+2, +1) + PX(-1, +2) + PX(+1, +2);
96 C3 = PX(-2, -2) + PX(+2, -2) + PX(-2, +2) + PX(+2, +2);
97 D1 = PX(+0, -3) + PX(-3, +0) + PX(+3, +0) + PX(+0, +3);
98 D2 = PX(-1, -3) + PX(+1, -3) + PX(-3, -1) + PX(+3, -1) + PX(-3, +1) + PX(+3, +1) + PX(-1, +3) + PX(+1, +3);
99 D3 = PX(-4, -2) + PX(-3, -2) + PX(+3, -2) + PX(+4, -2) + PX(-4, -1) + PX(+4, -1) + PX(-4, +0) + PX(+4, +0) + PX(-4,
100 +1) + PX(+4, +1) + PX(-4, +2) + PX(-3, +2) + PX(+3, +2) + PX(+4, +2);
101#undef PX
102 int i;
103 const float *uptr;
104
105 uptr = src + width * (y - 4) + (x - 4);
106 for (i = 0; i < 9; i++)
107 D3 += *uptr++;
108
109 uptr = src + width * (y - 3) + (x - 4);
110 for (i = 0; i < 3; i++)
111 D3 += *uptr++;
112 uptr += 3;
113 for (i = 0; i < 3; i++)
114 D3 += *uptr++;
115
116 uptr = src + width * (y + 3) + (x - 4);
117 for (i = 0; i < 3; i++)
118 D3 += *uptr++;
119 uptr += 3;
120 for (i = 0; i < 3; i++)
121 D3 += *uptr++;
122
123 uptr = src + width * (y + 4) + (x - 4);
124 for (i = 0; i < 9; i++)
125 D3 += *uptr++;
126
127 double mean = (A + B1 + B2 + C1 + C2 + C3 + D1 + D2 + D3) / 81.0;
128 double PSF_fit = PSF[0] * (A - mean) + PSF[1] * (B1 - 4.0 * mean) + PSF[2] * (B2 - 4.0 * mean) +
129 PSF[3] * (C1 - 4.0 * mean) + PSF[4] * (C2 - 8.0 * mean) + PSF[5] * (C3 - 4.0 * mean) +
130 PSF[6] * (D1 - 4.0 * mean) + PSF[7] * (D2 - 8.0 * mean) + PSF[8] * (D3 - 44.0 * mean);
131
132 dst[width * y + x] = (float) PSF_fit;
133 }
134 }
135}
136
137void GetStats(double *mean, double *stdev, int width, const float *img, const QRect &win)
138{
139 // Determine the mean and standard deviation
140 double sum = 0.0;
141 double a = 0.0;
142 double q = 0.0;
143 double k = 1.0;
144 double km1 = 0.0;
145
146 const float *p0 = img + win.top() * width + win.left();
147 for (int y = 0; y < win.height(); y++)
148 {
149 const float *end = p0 + win.height();
150 for (const float *p = p0; p < end; p++)
151 {
152 double const x = (double) * p;
153 sum += x;
154 double const a0 = a;
155 a += (x - a) / k;
156 q += (x - a0) * (x - a);
157 km1 = k;
158 k += 1.0;
159 }
160 p0 += width;
161 }
162
163 *mean = sum / km1;
164 *stdev = sqrt(q / km1);
165}
166
167void RemoveItems(std::set<Peak> &stars, const std::set<int> &to_erase)
168{
169 int n = 0;
170 for (std::set<Peak>::iterator it = stars.begin(); it != stars.end(); n++)
171 {
172 if (to_erase.find(n) != to_erase.end())
173 {
174 std::set<Peak>::iterator next = it;
175 ++next;
176 stars.erase(it);
177 it = next;
178 }
179 else
180 ++it;
181 }
182}
183
184float *createFloatImage(const QSharedPointer<FITSData> &target)
185{
186 QSharedPointer<FITSData> imageData;
187
188 /*************
189 if (target.isNull())
190 imageData = m_ImageData;
191 else
192 *********/ // shouldn't be null
193 imageData = target;
194
195 // #1 Convert to float array
196 // We only process 1st plane if it is a color image
197 uint32_t imgSize = imageData->width() * imageData->height();
198 float *imgFloat = new float[imgSize];
199
200 if (imgFloat == nullptr)
201 {
202 qCritical() << "Not enough memory for float image array!";
203 return nullptr;
204 }
205
206 switch (imageData->getStatistics().dataType)
207 {
208 case TBYTE:
209 {
210 uint8_t const *buffer = imageData->getImageBuffer();
211 for (uint32_t i = 0; i < imgSize; i++)
212 imgFloat[i] = buffer[i];
213 }
214 break;
215
216 case TSHORT:
217 {
218 int16_t const *buffer = reinterpret_cast<int16_t const *>(imageData->getImageBuffer());
219 for (uint32_t i = 0; i < imgSize; i++)
220 imgFloat[i] = buffer[i];
221 }
222 break;
223
224 case TUSHORT:
225 {
226 uint16_t const *buffer = reinterpret_cast<uint16_t const*>(imageData->getImageBuffer());
227 for (uint32_t i = 0; i < imgSize; i++)
228 imgFloat[i] = buffer[i];
229 }
230 break;
231
232 case TLONG:
233 {
234 int32_t const *buffer = reinterpret_cast<int32_t const*>(imageData->getImageBuffer());
235 for (uint32_t i = 0; i < imgSize; i++)
236 imgFloat[i] = buffer[i];
237 }
238 break;
239
240 case TULONG:
241 {
242 uint32_t const *buffer = reinterpret_cast<uint32_t const*>(imageData->getImageBuffer());
243 for (uint32_t i = 0; i < imgSize; i++)
244 imgFloat[i] = buffer[i];
245 }
246 break;
247
248 case TFLOAT:
249 {
250 float const *buffer = reinterpret_cast<float const*>(imageData->getImageBuffer());
251 for (uint32_t i = 0; i < imgSize; i++)
252 imgFloat[i] = buffer[i];
253 }
254 break;
255
256 case TLONGLONG:
257 {
258 int64_t const *buffer = reinterpret_cast<int64_t const*>(imageData->getImageBuffer());
259 for (uint32_t i = 0; i < imgSize; i++)
260 imgFloat[i] = buffer[i];
261 }
262 break;
263
264 case TDOUBLE:
265 {
266 double const *buffer = reinterpret_cast<double const*>(imageData->getImageBuffer());
267 for (uint32_t i = 0; i < imgSize; i++)
268 imgFloat[i] = buffer[i];
269 }
270 break;
271
272 default:
273 delete[] imgFloat;
274 return nullptr;
275 }
276
277 return imgFloat;
278}
279
280} // namespace
281
282// Based on PHD2 algorithm
283QList<Edge*> GuideAlgorithms::detectStars(const QSharedPointer<FITSData> &imageData, const QRect &trackingBox)
284{
285 //Debug.Write(wxString::Format("Star::AutoFind called with edgeAllowance = %d searchRegion = %d\n", extraEdgeAllowance, searchRegion));
286
287 // run a 3x3 median first to eliminate hot pixels
288 //usImage smoothed;
289 //smoothed.CopyFrom(image);
290 //Median3(smoothed);
291 constexpr int extraEdgeAllowance = 0;
292 const QSharedPointer<FITSData> smoothed(new FITSData(imageData));
293 smoothed->applyFilter(FITS_MEDIAN);
294
295 int searchRegion = trackingBox.width();
296
297 int subW = smoothed->width();
298 int subH = smoothed->height();
299 int size = subW * subH;
300
301 // convert to floating point
302 float *conv = createFloatImage(smoothed);
303
304 // run the PSF convolution
305 {
306 float *tmp = new float[size] {};
307 psf_conv(tmp, conv, subW, subH);
308 delete [] conv;
309 // Swap
310 conv = tmp;
311 }
312
313 enum { CONV_RADIUS = 4 };
314 int dw = subW; // width of the downsampled image
315 int dh = subH; // height of the downsampled image
316 QRect convRect(CONV_RADIUS, CONV_RADIUS, dw - 2 * CONV_RADIUS, dh - 2 * CONV_RADIUS); // region containing valid data
317
318 enum { TOP_N = 100 }; // keep track of the brightest stars
319 std::set<Peak> stars; // sorted by ascending intensity
320
321 double global_mean, global_stdev;
322 GetStats(&global_mean, &global_stdev, subW, conv, convRect);
323
324 //Debug.Write(wxString::Format("AutoFind: global mean = %.1f, stdev %.1f\n", global_mean, global_stdev));
325
326 const double threshold = 0.1;
327 //Debug.Write(wxString::Format("AutoFind: using threshold = %.1f\n", threshold));
328
329 // find each local maximum
330 int srch = 4;
331 for (int y = convRect.top() + srch; y <= convRect.bottom() - srch; y++)
332 {
333 for (int x = convRect.left() + srch; x <= convRect.right() - srch; x++)
334 {
335 float val = conv[dw * y + x];
336 bool ismax = false;
337 if (val > 0.0)
338 {
339 ismax = true;
340 for (int j = -srch; j <= srch; j++)
341 {
342 for (int i = -srch; i <= srch; i++)
343 {
344 if (i == 0 && j == 0)
345 continue;
346 if (conv[dw * (y + j) + (x + i)] > val)
347 {
348 ismax = false;
349 break;
350 }
351 }
352 }
353 }
354 if (!ismax)
355 continue;
356
357 // compare local maximum to mean value of surrounding pixels
358 const int local = 7;
359 double local_mean, local_stdev;
360 QRect localRect(x - local, y - local, 2 * local + 1, 2 * local + 1);
361 localRect = localRect.intersected(convRect);
362 GetStats(&local_mean, &local_stdev, subW, conv, localRect);
363
364 // this is our measure of star intensity
365 double h = (val - local_mean) / global_stdev;
366
367 if (h < threshold)
368 {
369 // Debug.Write(wxString::Format("AG: local max REJECT [%d, %d] PSF %.1f SNR %.1f\n", imgx, imgy, val, SNR));
370 continue;
371 }
372
373 // coordinates on the original image
374 int downsample = 1;
375 int imgx = x * downsample + downsample / 2;
376 int imgy = y * downsample + downsample / 2;
377
378 stars.insert(Peak(imgx, imgy, h));
379 if (stars.size() > TOP_N)
380 stars.erase(stars.begin());
381 }
382 }
383
384 //for (std::set<Peak>::const_reverse_iterator it = stars.rbegin(); it != stars.rend(); ++it)
385 //qCDebug(KSTARS_EKOS_GUIDE) << "AutoFind: local max [" << it->x << "," << it->y << "]" << it->val;
386
387 // merge stars that are very close into a single star
388 {
389 const int minlimitsq = 5 * 5;
390repeat:
391 for (std::set<Peak>::const_iterator a = stars.begin(); a != stars.end(); ++a)
392 {
393 std::set<Peak>::const_iterator b = a;
394 ++b;
395 for (; b != stars.end(); ++b)
396 {
397 int dx = a->x - b->x;
398 int dy = a->y - b->y;
399 int d2 = dx * dx + dy * dy;
400 if (d2 < minlimitsq)
401 {
402 // very close, treat as single star
403 //Debug.Write(wxString::Format("AutoFind: merge [%d, %d] %.1f - [%d, %d] %.1f\n", a->x, a->y, a->val, b->x, b->y, b->val));
404 // erase the dimmer one
405 stars.erase(a);
406 goto repeat;
407 }
408 }
409 }
410 }
411
412 // exclude stars that would fit within a single searchRegion box
413 {
414 // build a list of stars to be excluded
415 std::set<int> to_erase;
416 const int extra = 5; // extra safety margin
417 const int fullw = searchRegion + extra;
418 for (std::set<Peak>::const_iterator a = stars.begin(); a != stars.end(); ++a)
419 {
420 std::set<Peak>::const_iterator b = a;
421 ++b;
422 for (; b != stars.end(); ++b)
423 {
424 int dx = abs(a->x - b->x);
425 int dy = abs(a->y - b->y);
426 if (dx <= fullw && dy <= fullw)
427 {
428 // stars closer than search region, exclude them both
429 // but do not let a very dim star eliminate a very bright star
430 if (b->val / a->val >= 5.0)
431 {
432 //Debug.Write(wxString::Format("AutoFind: close dim-bright [%d, %d] %.1f - [%d, %d] %.1f\n", a->x, a->y, a->val, b->x, b->y, b->val));
433 }
434 else
435 {
436 //Debug.Write(wxString::Format("AutoFind: too close [%d, %d] %.1f - [%d, %d] %.1f\n", a->x, a->y, a->val, b->x, b->y, b->val));
437 to_erase.insert(std::distance(stars.begin(), a));
438 to_erase.insert(std::distance(stars.begin(), b));
439 }
440 }
441 }
442 }
443 RemoveItems(stars, to_erase);
444 }
445
446 // exclude stars too close to the edge
447 {
448 enum { MIN_EDGE_DIST = 40 };
449 int edgeDist = MIN_EDGE_DIST;//pConfig->Profile.GetInt("/StarAutoFind/MinEdgeDist", MIN_EDGE_DIST);
450 if (edgeDist < searchRegion)
451 edgeDist = searchRegion;
452 edgeDist += extraEdgeAllowance;
453
454 std::set<Peak>::iterator it = stars.begin();
455 while (it != stars.end())
456 {
457 std::set<Peak>::iterator next = it;
458 ++next;
459 if (it->x <= edgeDist || it->x >= subW - edgeDist ||
460 it->y <= edgeDist || it->y >= subH - edgeDist)
461 {
462 //Debug.Write(wxString::Format("AutoFind: too close to edge [%d, %d] %.1f\n", it->x, it->y, it->val));
463 stars.erase(it);
464 }
465 it = next;
466 }
467 }
468
469 QList<Edge*> centers;
470 for (std::set<Peak>::reverse_iterator it = stars.rbegin(); it != stars.rend(); ++it)
471 {
472 Edge *center = new Edge;
473 center->x = it->x;
474 center->y = it->y;
475 center->val = it->val;
476 centers.append(center);
477 }
478
479 delete [] conv;
480
481 return centers;
482}
483
484
485template <typename T>
486GuiderUtils::Vector GuideAlgorithms::findLocalStarPosition(QSharedPointer<FITSData> &imageData,
487 const int algorithmIndex,
488 const int videoWidth,
489 const int videoHeight,
490 const QRect &trackingBox)
491{
492 static const double P0 = 0.906, P1 = 0.584, P2 = 0.365, P3 = 0.117, P4 = 0.049, P5 = -0.05, P6 = -0.064, P7 = -0.074,
493 P8 = -0.094;
494
495 GuiderUtils::Vector ret(-1, -1, -1);
496 int i, j;
497 double resx, resy, mass, threshold, pval;
498 T const *psrc = nullptr;
499 T const *porigin = nullptr;
500 T const *pptr;
501
502 if (trackingBox.isValid() == false)
503 return GuiderUtils::Vector(-1, -1, -1);
504
505 if (imageData.isNull())
506 {
507 qCWarning(KSTARS_EKOS_GUIDE) << "Cannot process a nullptr image.";
508 return GuiderUtils::Vector(-1, -1, -1);
509 }
510
511 if (algorithmIndex == SEP_THRESHOLD)
512 {
513 QVariantMap settings;
514 settings["optionsProfileIndex"] = Options::guideOptionsProfile();
515 settings["optionsProfileGroup"] = static_cast<int>(Ekos::GuideProfiles);
516 imageData->setSourceExtractorSettings(settings);
517 QFuture<bool> result = imageData->findStars(ALGORITHM_SEP, trackingBox);
518 result.waitForFinished();
519 if (result.result())
520 {
521 imageData->getHFR(HFR_MEDIAN);
522 const Edge star = imageData->getSelectedHFRStar();
523 ret = GuiderUtils::Vector(star.x, star.y, 0);
524 }
525 else
526 ret = GuiderUtils::Vector(-1, -1, -1);
527
528 return ret;
529 }
530
531 T const *pdata = reinterpret_cast<T const*>(imageData->getImageBuffer());
532
533 qCDebug(KSTARS_EKOS_GUIDE) << "Tracking Square " << trackingBox;
534
535 double square_square = trackingBox.width() * trackingBox.width();
536
537 psrc = porigin = pdata + trackingBox.y() * videoWidth + trackingBox.x();
538
539 resx = resy = 0;
540 threshold = mass = 0;
541
542 // several threshold adaptive smart algorithms
543 switch (algorithmIndex)
544 {
545 case CENTROID_THRESHOLD:
546 {
547 int width = trackingBox.width();
548 int height = trackingBox.width();
549 float i0, i1, i2, i3, i4, i5, i6, i7, i8;
550 int ix = 0, iy = 0;
551 int xM4;
552 T const *p;
553 double average, fit, bestFit = 0;
554 int minx = 0;
555 int maxx = width;
556 int miny = 0;
557 int maxy = height;
558 for (int x = minx; x < maxx; x++)
559 for (int y = miny; y < maxy; y++)
560 {
561 i0 = i1 = i2 = i3 = i4 = i5 = i6 = i7 = i8 = 0;
562 xM4 = x - 4;
563 p = psrc + (y - 4) * videoWidth + xM4;
564 i8 += *p++;
565 i8 += *p++;
566 i8 += *p++;
567 i8 += *p++;
568 i8 += *p++;
569 i8 += *p++;
570 i8 += *p++;
571 i8 += *p++;
572 i8 += *p++;
573 p = psrc + (y - 3) * videoWidth + xM4;
574 i8 += *p++;
575 i8 += *p++;
576 i8 += *p++;
577 i7 += *p++;
578 i6 += *p++;
579 i7 += *p++;
580 i8 += *p++;
581 i8 += *p++;
582 i8 += *p++;
583 p = psrc + (y - 2) * videoWidth + xM4;
584 i8 += *p++;
585 i8 += *p++;
586 i5 += *p++;
587 i4 += *p++;
588 i3 += *p++;
589 i4 += *p++;
590 i5 += *p++;
591 i8 += *p++;
592 i8 += *p++;
593 p = psrc + (y - 1) * videoWidth + xM4;
594 i8 += *p++;
595 i7 += *p++;
596 i4 += *p++;
597 i2 += *p++;
598 i1 += *p++;
599 i2 += *p++;
600 i4 += *p++;
601 i8 += *p++;
602 i8 += *p++;
603 p = psrc + (y + 0) * videoWidth + xM4;
604 i8 += *p++;
605 i6 += *p++;
606 i3 += *p++;
607 i1 += *p++;
608 i0 += *p++;
609 i1 += *p++;
610 i3 += *p++;
611 i6 += *p++;
612 i8 += *p++;
613 p = psrc + (y + 1) * videoWidth + xM4;
614 i8 += *p++;
615 i7 += *p++;
616 i4 += *p++;
617 i2 += *p++;
618 i1 += *p++;
619 i2 += *p++;
620 i4 += *p++;
621 i8 += *p++;
622 i8 += *p++;
623 p = psrc + (y + 2) * videoWidth + xM4;
624 i8 += *p++;
625 i8 += *p++;
626 i5 += *p++;
627 i4 += *p++;
628 i3 += *p++;
629 i4 += *p++;
630 i5 += *p++;
631 i8 += *p++;
632 i8 += *p++;
633 p = psrc + (y + 3) * videoWidth + xM4;
634 i8 += *p++;
635 i8 += *p++;
636 i8 += *p++;
637 i7 += *p++;
638 i6 += *p++;
639 i7 += *p++;
640 i8 += *p++;
641 i8 += *p++;
642 i8 += *p++;
643 p = psrc + (y + 4) * videoWidth + xM4;
644 i8 += *p++;
645 i8 += *p++;
646 i8 += *p++;
647 i8 += *p++;
648 i8 += *p++;
649 i8 += *p++;
650 i8 += *p++;
651 i8 += *p++;
652 i8 += *p++;
653 average = (i0 + i1 + i2 + i3 + i4 + i5 + i6 + i7 + i8) / 85.0;
654 fit = P0 * (i0 - average) + P1 * (i1 - 4 * average) + P2 * (i2 - 4 * average) +
655 P3 * (i3 - 4 * average) + P4 * (i4 - 8 * average) + P5 * (i5 - 4 * average) +
656 P6 * (i6 - 4 * average) + P7 * (i7 - 8 * average) + P8 * (i8 - 48 * average);
657 if (bestFit < fit)
658 {
659 bestFit = fit;
660 ix = x;
661 iy = y;
662 }
663 }
664
665 if (bestFit > 50)
666 {
667 double sumX = 0;
668 double sumY = 0;
669 double total = 0;
670 for (int y = iy - 4; y <= iy + 4; y++)
671 {
672 p = psrc + y * width + ix - 4;
673 for (int x = ix - 4; x <= ix + 4; x++)
674 {
675 double w = *p++;
676 sumX += x * w;
677 sumY += y * w;
678 total += w;
679 }
680 }
681 if (total > 0)
682 {
683 ret = (GuiderUtils::Vector(trackingBox.x(), trackingBox.y(), 0) + GuiderUtils::Vector(sumX / total, sumY / total, 0));
684 return ret;
685 }
686 }
687
688 return GuiderUtils::Vector(-1, -1, -1);
689 }
690 break;
691 // Alexander's Stepanenko smart threshold algorithm
692 case SMART_THRESHOLD:
693 {
694 point_t bbox_lt = { trackingBox.x() - SMART_FRAME_WIDTH, trackingBox.y() - SMART_FRAME_WIDTH };
695 point_t bbox_rb = { trackingBox.x() + trackingBox.width() + SMART_FRAME_WIDTH,
696 trackingBox.y() + trackingBox.width() + SMART_FRAME_WIDTH
697 };
698 int offset = 0;
699
700 // clip frame
701 if (bbox_lt.x < 0)
702 bbox_lt.x = 0;
703 if (bbox_lt.y < 0)
704 bbox_lt.y = 0;
705 if (bbox_rb.x > videoWidth)
706 bbox_rb.x = videoWidth;
707 if (bbox_rb.y > videoHeight)
708 bbox_rb.y = videoHeight;
709
710 // calc top bar
711 int box_wd = bbox_rb.x - bbox_lt.x;
712 int box_ht = trackingBox.y() - bbox_lt.y;
713 int pix_cnt = 0;
714 if (box_wd > 0 && box_ht > 0)
715 {
716 pix_cnt += box_wd * box_ht;
717 for (j = bbox_lt.y; j < trackingBox.y(); ++j)
718 {
719 offset = j * videoWidth;
720 for (i = bbox_lt.x; i < bbox_rb.x; ++i)
721 {
722 pptr = pdata + offset + i;
723 threshold += *pptr;
724 }
725 }
726 }
727 // calc left bar
728 box_wd = trackingBox.x() - bbox_lt.x;
729 box_ht = trackingBox.width();
730 if (box_wd > 0 && box_ht > 0)
731 {
732 pix_cnt += box_wd * box_ht;
733 for (j = trackingBox.y(); j < trackingBox.y() + box_ht; ++j)
734 {
735 offset = j * videoWidth;
736 for (i = bbox_lt.x; i < trackingBox.x(); ++i)
737 {
738 pptr = pdata + offset + i;
739 threshold += *pptr;
740 }
741 }
742 }
743 // calc right bar
744 box_wd = bbox_rb.x - trackingBox.x() - trackingBox.width();
745 box_ht = trackingBox.width();
746 if (box_wd > 0 && box_ht > 0)
747 {
748 pix_cnt += box_wd * box_ht;
749 for (j = trackingBox.y(); j < trackingBox.y() + box_ht; ++j)
750 {
751 offset = j * videoWidth;
752 for (i = trackingBox.x() + trackingBox.width(); i < bbox_rb.x; ++i)
753 {
754 pptr = pdata + offset + i;
755 threshold += *pptr;
756 }
757 }
758 }
759 // calc bottom bar
760 box_wd = bbox_rb.x - bbox_lt.x;
761 box_ht = bbox_rb.y - trackingBox.y() - trackingBox.width();
762 if (box_wd > 0 && box_ht > 0)
763 {
764 pix_cnt += box_wd * box_ht;
765 for (j = trackingBox.y() + trackingBox.width(); j < bbox_rb.y; ++j)
766 {
767 offset = j * videoWidth;
768 for (i = bbox_lt.x; i < bbox_rb.x; ++i)
769 {
770 pptr = pdata + offset + i;
771 threshold += *pptr;
772 }
773 }
774 }
775 // find maximum
776 double max_val = 0;
777 for (j = 0; j < trackingBox.width(); ++j)
778 {
779 for (i = 0; i < trackingBox.width(); ++i)
780 {
781 pptr = psrc + i;
782 if (*pptr > max_val)
783 max_val = *pptr;
784 }
785 psrc += videoWidth;
786 }
787 if (pix_cnt != 0)
788 threshold /= (double)pix_cnt;
789
790 // cut by 10% higher then average threshold
791 if (max_val > threshold)
792 threshold += (max_val - threshold) * SMART_CUT_FACTOR;
793
794 //log_i("smart thr. = %f cnt = %d", threshold, pix_cnt);
795 break;
796 }
797 // simple adaptive threshold
798 case AUTO_THRESHOLD:
799 {
800 for (j = 0; j < trackingBox.width(); ++j)
801 {
802 for (i = 0; i < trackingBox.width(); ++i)
803 {
804 pptr = psrc + i;
805 threshold += *pptr;
806 }
807 psrc += videoWidth;
808 }
809 threshold /= square_square;
810 break;
811 }
812 // no threshold subtracion
813 default:
814 {
815 }
816 }
817
818 psrc = porigin;
819 for (j = 0; j < trackingBox.width(); ++j)
820 {
821 for (i = 0; i < trackingBox.width(); ++i)
822 {
823 pptr = psrc + i;
824 pval = *pptr - threshold;
825 pval = pval < 0 ? 0 : pval;
826
827 resx += (double)i * pval;
828 resy += (double)j * pval;
829
830 mass += pval;
831 }
832 psrc += videoWidth;
833 }
834
835 if (mass == 0)
836 mass = 1;
837
838 resx /= mass;
839 resy /= mass;
840
841 ret = GuiderUtils::Vector(trackingBox.x(), trackingBox.y(), 0) + GuiderUtils::Vector(resx, resy, 0);
842
843 return ret;
844}
845
846GuiderUtils::Vector GuideAlgorithms::findLocalStarPosition(QSharedPointer<FITSData> &imageData,
847 const int algorithmIndex,
848 const int videoWidth,
849 const int videoHeight,
850 const QRect &trackingBox)
851{
852 switch (imageData->dataType())
853 {
854 case TBYTE:
855 return findLocalStarPosition<uint8_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
856
857 case TSHORT:
858 return findLocalStarPosition<int16_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
859
860 case TUSHORT:
861 return findLocalStarPosition<uint16_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
862
863 case TLONG:
864 return findLocalStarPosition<int32_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
865
866 case TULONG:
867 return findLocalStarPosition<uint32_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
868
869 case TFLOAT:
870 return findLocalStarPosition<float>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
871
872 case TLONGLONG:
873 return findLocalStarPosition<int64_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
874
875 case TDOUBLE:
876 return findLocalStarPosition<double>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
877
878 default:
879 break;
880 }
881
882 return GuiderUtils::Vector(-1, -1, -1);
883}
KGameStandardAction::repeat
QAction * repeat(const QObject *recvr, const char *slot, QObject *parent)
KGameStandardAction::end
QAction * end(const QObject *recvr, const char *slot, QObject *parent)
KStandardAction::next
QAction * next(const QObject *recvr, const char *slot, QObject *parent)
QFuture::result
T result() const const
QFuture::waitForFinished
void waitForFinished()
QList::append
void append(QList< T > &&value)
QRect::height
int height() const const
QRect::isValid
bool isValid() const const
QRect::left
int left() const const
QRect::top
int top() const const
QRect::width
int width() const const
QRect::x
int x() const const
QRect::y
int y() const const
QSharedPointer::isNull
bool isNull() const const
Qt::center
QTextStream & center(QTextStream &stream)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.