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
19 #define SMART_THRESHOLD 0
20 #define SEP_THRESHOLD 1
21 #define CENTROID_THRESHOLD 2
22 #define AUTO_THRESHOLD 3
23 #define NO_THRESHOLD 4
24 #define SEP_MULTISTAR 5
25
26 // smart threshold algorithm param
27 // width of outer frame for background calculation
28 #define SMART_FRAME_WIDTH 4
29 // cut-factor above average threshold
30 #define SMART_CUT_FACTOR 0.1
31
32 struct Peak
33 {
34 int x;
35 int y;
36 float val;
37
38 Peak() = default;
PeakPeak39 Peak(int x_, int y_, float val_) : x(x_), y(y_), val(val_) { }
operator <Peak40 bool operator<(const Peak &rhs) const
41 {
42 return val < rhs.val;
43 }
44 };
45
46 // JM: Why not use QPoint?
47 typedef struct
48 {
49 int x, y;
50 } point_t;
51
52 namespace
53 {
54
psf_conv(float * dst,const float * src,int width,int height)55 void psf_conv(float *dst, const float *src, int width, int height)
56 {
57 //dst.Init(src.Size);
58
59 // A B1 B2 C1 C2 C3 D1 D2 D3
60 const double PSF[] = { 0.906, 0.584, 0.365, .117, .049, -0.05, -.064, -.074, -.094 };
61
62 //memset(dst.px, 0, src.NPixels * sizeof(float));
63
64 /* PSF Grid is:
65 D3 D3 D3 D3 D3 D3 D3 D3 D3
66 D3 D3 D3 D2 D1 D2 D3 D3 D3
67 D3 D3 C3 C2 C1 C2 C3 D3 D3
68 D3 D2 C2 B2 B1 B2 C2 D2 D3
69 D3 D1 C1 B1 A B1 C1 D1 D3
70 D3 D2 C2 B2 B1 B2 C2 D2 D3
71 D3 D3 C3 C2 C1 C2 C3 D3 D3
72 D3 D3 D3 D2 D1 D2 D3 D3 D3
73 D3 D3 D3 D3 D3 D3 D3 D3 D3
74
75 1@A
76 4@B1, B2, C1, C3, D1
77 8@C2, D2
78 44 * D3
79 */
80
81 int psf_size = 4;
82
83 for (int y = psf_size; y < height - psf_size; y++)
84 {
85 for (int x = psf_size; x < width - psf_size; x++)
86 {
87 float A, B1, B2, C1, C2, C3, D1, D2, D3;
88
89 #define PX(dx, dy) *(src + width * (y + (dy)) + x + (dx))
90 A = PX(+0, +0);
91 B1 = PX(+0, -1) + PX(+0, +1) + PX(+1, +0) + PX(-1, +0);
92 B2 = PX(-1, -1) + PX(+1, -1) + PX(-1, +1) + PX(+1, +1);
93 C1 = PX(+0, -2) + PX(-2, +0) + PX(+2, +0) + PX(+0, +2);
94 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);
95 C3 = PX(-2, -2) + PX(+2, -2) + PX(-2, +2) + PX(+2, +2);
96 D1 = PX(+0, -3) + PX(-3, +0) + PX(+3, +0) + PX(+0, +3);
97 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);
98 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,
99 +1) + PX(+4, +1) + PX(-4, +2) + PX(-3, +2) + PX(+3, +2) + PX(+4, +2);
100 #undef PX
101 int i;
102 const float *uptr;
103
104 uptr = src + width * (y - 4) + (x - 4);
105 for (i = 0; i < 9; i++)
106 D3 += *uptr++;
107
108 uptr = src + width * (y - 3) + (x - 4);
109 for (i = 0; i < 3; i++)
110 D3 += *uptr++;
111 uptr += 3;
112 for (i = 0; i < 3; i++)
113 D3 += *uptr++;
114
115 uptr = src + width * (y + 3) + (x - 4);
116 for (i = 0; i < 3; i++)
117 D3 += *uptr++;
118 uptr += 3;
119 for (i = 0; i < 3; i++)
120 D3 += *uptr++;
121
122 uptr = src + width * (y + 4) + (x - 4);
123 for (i = 0; i < 9; i++)
124 D3 += *uptr++;
125
126 double mean = (A + B1 + B2 + C1 + C2 + C3 + D1 + D2 + D3) / 81.0;
127 double PSF_fit = PSF[0] * (A - mean) + PSF[1] * (B1 - 4.0 * mean) + PSF[2] * (B2 - 4.0 * mean) +
128 PSF[3] * (C1 - 4.0 * mean) + PSF[4] * (C2 - 8.0 * mean) + PSF[5] * (C3 - 4.0 * mean) +
129 PSF[6] * (D1 - 4.0 * mean) + PSF[7] * (D2 - 8.0 * mean) + PSF[8] * (D3 - 44.0 * mean);
130
131 dst[width * y + x] = (float) PSF_fit;
132 }
133 }
134 }
135
GetStats(double * mean,double * stdev,int width,const float * img,const QRect & win)136 void GetStats(double *mean, double *stdev, int width, const float *img, const QRect &win)
137 {
138 // Determine the mean and standard deviation
139 double sum = 0.0;
140 double a = 0.0;
141 double q = 0.0;
142 double k = 1.0;
143 double km1 = 0.0;
144
145 const float *p0 = img + win.top() * width + win.left();
146 for (int y = 0; y < win.height(); y++)
147 {
148 const float *end = p0 + win.height();
149 for (const float *p = p0; p < end; p++)
150 {
151 double const x = (double) * p;
152 sum += x;
153 double const a0 = a;
154 a += (x - a) / k;
155 q += (x - a0) * (x - a);
156 km1 = k;
157 k += 1.0;
158 }
159 p0 += width;
160 }
161
162 *mean = sum / km1;
163 *stdev = sqrt(q / km1);
164 }
165
RemoveItems(std::set<Peak> & stars,const std::set<int> & to_erase)166 void RemoveItems(std::set<Peak> &stars, const std::set<int> &to_erase)
167 {
168 int n = 0;
169 for (std::set<Peak>::iterator it = stars.begin(); it != stars.end(); n++)
170 {
171 if (to_erase.find(n) != to_erase.end())
172 {
173 std::set<Peak>::iterator next = it;
174 ++next;
175 stars.erase(it);
176 it = next;
177 }
178 else
179 ++it;
180 }
181 }
182
createFloatImage(const QSharedPointer<FITSData> & target)183 float *createFloatImage(const QSharedPointer<FITSData> &target)
184 {
185 QSharedPointer<FITSData> imageData;
186
187 /*************
188 if (target.isNull())
189 imageData = m_ImageData;
190 else
191 *********/ // shouldn't be null
192 imageData = target;
193
194 // #1 Convert to float array
195 // We only process 1st plane if it is a color image
196 uint32_t imgSize = imageData->width() * imageData->height();
197 float *imgFloat = new float[imgSize];
198
199 if (imgFloat == nullptr)
200 {
201 qCritical() << "Not enough memory for float image array!";
202 return nullptr;
203 }
204
205 switch (imageData->getStatistics().dataType)
206 {
207 case TBYTE:
208 {
209 uint8_t const *buffer = imageData->getImageBuffer();
210 for (uint32_t i = 0; i < imgSize; i++)
211 imgFloat[i] = buffer[i];
212 }
213 break;
214
215 case TSHORT:
216 {
217 int16_t const *buffer = reinterpret_cast<int16_t const *>(imageData->getImageBuffer());
218 for (uint32_t i = 0; i < imgSize; i++)
219 imgFloat[i] = buffer[i];
220 }
221 break;
222
223 case TUSHORT:
224 {
225 uint16_t const *buffer = reinterpret_cast<uint16_t const*>(imageData->getImageBuffer());
226 for (uint32_t i = 0; i < imgSize; i++)
227 imgFloat[i] = buffer[i];
228 }
229 break;
230
231 case TLONG:
232 {
233 int32_t const *buffer = reinterpret_cast<int32_t const*>(imageData->getImageBuffer());
234 for (uint32_t i = 0; i < imgSize; i++)
235 imgFloat[i] = buffer[i];
236 }
237 break;
238
239 case TULONG:
240 {
241 uint32_t const *buffer = reinterpret_cast<uint32_t const*>(imageData->getImageBuffer());
242 for (uint32_t i = 0; i < imgSize; i++)
243 imgFloat[i] = buffer[i];
244 }
245 break;
246
247 case TFLOAT:
248 {
249 float const *buffer = reinterpret_cast<float const*>(imageData->getImageBuffer());
250 for (uint32_t i = 0; i < imgSize; i++)
251 imgFloat[i] = buffer[i];
252 }
253 break;
254
255 case TLONGLONG:
256 {
257 int64_t const *buffer = reinterpret_cast<int64_t const*>(imageData->getImageBuffer());
258 for (uint32_t i = 0; i < imgSize; i++)
259 imgFloat[i] = buffer[i];
260 }
261 break;
262
263 case TDOUBLE:
264 {
265 double const *buffer = reinterpret_cast<double const*>(imageData->getImageBuffer());
266 for (uint32_t i = 0; i < imgSize; i++)
267 imgFloat[i] = buffer[i];
268 }
269 break;
270
271 default:
272 delete[] imgFloat;
273 return nullptr;
274 }
275
276 return imgFloat;
277 }
278
279 } // namespace
280
281 // Based on PHD2 algorithm
detectStars(const QSharedPointer<FITSData> & imageData,const QRect & trackingBox)282 QList<Edge*> GuideAlgorithms::detectStars(const QSharedPointer<FITSData> &imageData, const QRect &trackingBox)
283 {
284 //Debug.Write(wxString::Format("Star::AutoFind called with edgeAllowance = %d searchRegion = %d\n", extraEdgeAllowance, searchRegion));
285
286 // run a 3x3 median first to eliminate hot pixels
287 //usImage smoothed;
288 //smoothed.CopyFrom(image);
289 //Median3(smoothed);
290 constexpr int extraEdgeAllowance = 0;
291 const QSharedPointer<FITSData> smoothed(new FITSData(imageData));
292 smoothed->applyFilter(FITS_MEDIAN);
293
294 int searchRegion = trackingBox.width();
295
296 int subW = smoothed->width();
297 int subH = smoothed->height();
298 int size = subW * subH;
299
300 // convert to floating point
301 float *conv = createFloatImage(smoothed);
302
303 // run the PSF convolution
304 {
305 float *tmp = new float[size];
306 memset(tmp, 0, size * sizeof(float));
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;
390 repeat:
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
485 template <typename T>
findLocalStarPosition(QSharedPointer<FITSData> & imageData,const int algorithmIndex,const int videoWidth,const int videoHeight,const QRect & trackingBox)486 GuiderUtils::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;
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 Edge *star = imageData->getSelectedHFRStar();
523 if (star)
524 ret = GuiderUtils::Vector(star->x, star->y, 0);
525 else
526 ret = GuiderUtils::Vector(-1, -1, -1);
527 }
528 else
529 ret = GuiderUtils::Vector(-1, -1, -1);
530
531 return ret;
532 }
533
534 T const *pdata = reinterpret_cast<T const*>(imageData->getImageBuffer());
535
536 qCDebug(KSTARS_EKOS_GUIDE) << "Tracking Square " << trackingBox;
537
538 double square_square = trackingBox.width() * trackingBox.width();
539
540 psrc = porigin = pdata + trackingBox.y() * videoWidth + trackingBox.x();
541
542 resx = resy = 0;
543 threshold = mass = 0;
544
545 // several threshold adaptive smart algorithms
546 switch (algorithmIndex)
547 {
548 case CENTROID_THRESHOLD:
549 {
550 int width = trackingBox.width();
551 int height = trackingBox.width();
552 float i0, i1, i2, i3, i4, i5, i6, i7, i8;
553 int ix = 0, iy = 0;
554 int xM4;
555 T const *p;
556 double average, fit, bestFit = 0;
557 int minx = 0;
558 int maxx = width;
559 int miny = 0;
560 int maxy = height;
561 for (int x = minx; x < maxx; x++)
562 for (int y = miny; y < maxy; y++)
563 {
564 i0 = i1 = i2 = i3 = i4 = i5 = i6 = i7 = i8 = 0;
565 xM4 = x - 4;
566 p = psrc + (y - 4) * videoWidth + xM4;
567 i8 += *p++;
568 i8 += *p++;
569 i8 += *p++;
570 i8 += *p++;
571 i8 += *p++;
572 i8 += *p++;
573 i8 += *p++;
574 i8 += *p++;
575 i8 += *p++;
576 p = psrc + (y - 3) * videoWidth + xM4;
577 i8 += *p++;
578 i8 += *p++;
579 i8 += *p++;
580 i7 += *p++;
581 i6 += *p++;
582 i7 += *p++;
583 i8 += *p++;
584 i8 += *p++;
585 i8 += *p++;
586 p = psrc + (y - 2) * videoWidth + xM4;
587 i8 += *p++;
588 i8 += *p++;
589 i5 += *p++;
590 i4 += *p++;
591 i3 += *p++;
592 i4 += *p++;
593 i5 += *p++;
594 i8 += *p++;
595 i8 += *p++;
596 p = psrc + (y - 1) * videoWidth + xM4;
597 i8 += *p++;
598 i7 += *p++;
599 i4 += *p++;
600 i2 += *p++;
601 i1 += *p++;
602 i2 += *p++;
603 i4 += *p++;
604 i8 += *p++;
605 i8 += *p++;
606 p = psrc + (y + 0) * videoWidth + xM4;
607 i8 += *p++;
608 i6 += *p++;
609 i3 += *p++;
610 i1 += *p++;
611 i0 += *p++;
612 i1 += *p++;
613 i3 += *p++;
614 i6 += *p++;
615 i8 += *p++;
616 p = psrc + (y + 1) * videoWidth + xM4;
617 i8 += *p++;
618 i7 += *p++;
619 i4 += *p++;
620 i2 += *p++;
621 i1 += *p++;
622 i2 += *p++;
623 i4 += *p++;
624 i8 += *p++;
625 i8 += *p++;
626 p = psrc + (y + 2) * videoWidth + xM4;
627 i8 += *p++;
628 i8 += *p++;
629 i5 += *p++;
630 i4 += *p++;
631 i3 += *p++;
632 i4 += *p++;
633 i5 += *p++;
634 i8 += *p++;
635 i8 += *p++;
636 p = psrc + (y + 3) * videoWidth + xM4;
637 i8 += *p++;
638 i8 += *p++;
639 i8 += *p++;
640 i7 += *p++;
641 i6 += *p++;
642 i7 += *p++;
643 i8 += *p++;
644 i8 += *p++;
645 i8 += *p++;
646 p = psrc + (y + 4) * videoWidth + xM4;
647 i8 += *p++;
648 i8 += *p++;
649 i8 += *p++;
650 i8 += *p++;
651 i8 += *p++;
652 i8 += *p++;
653 i8 += *p++;
654 i8 += *p++;
655 i8 += *p++;
656 average = (i0 + i1 + i2 + i3 + i4 + i5 + i6 + i7 + i8) / 85.0;
657 fit = P0 * (i0 - average) + P1 * (i1 - 4 * average) + P2 * (i2 - 4 * average) +
658 P3 * (i3 - 4 * average) + P4 * (i4 - 8 * average) + P5 * (i5 - 4 * average) +
659 P6 * (i6 - 4 * average) + P7 * (i7 - 8 * average) + P8 * (i8 - 48 * average);
660 if (bestFit < fit)
661 {
662 bestFit = fit;
663 ix = x;
664 iy = y;
665 }
666 }
667
668 if (bestFit > 50)
669 {
670 double sumX = 0;
671 double sumY = 0;
672 double total = 0;
673 for (int y = iy - 4; y <= iy + 4; y++)
674 {
675 p = psrc + y * width + ix - 4;
676 for (int x = ix - 4; x <= ix + 4; x++)
677 {
678 double w = *p++;
679 sumX += x * w;
680 sumY += y * w;
681 total += w;
682 }
683 }
684 if (total > 0)
685 {
686 ret = (GuiderUtils::Vector(trackingBox.x(), trackingBox.y(), 0) + GuiderUtils::Vector(sumX / total, sumY / total, 0));
687 return ret;
688 }
689 }
690
691 return GuiderUtils::Vector(-1, -1, -1);
692 }
693 break;
694 // Alexander's Stepanenko smart threshold algorithm
695 case SMART_THRESHOLD:
696 {
697 point_t bbox_lt = { trackingBox.x() - SMART_FRAME_WIDTH, trackingBox.y() - SMART_FRAME_WIDTH };
698 point_t bbox_rb = { trackingBox.x() + trackingBox.width() + SMART_FRAME_WIDTH,
699 trackingBox.y() + trackingBox.width() + SMART_FRAME_WIDTH
700 };
701 int offset = 0;
702
703 // clip frame
704 if (bbox_lt.x < 0)
705 bbox_lt.x = 0;
706 if (bbox_lt.y < 0)
707 bbox_lt.y = 0;
708 if (bbox_rb.x > videoWidth)
709 bbox_rb.x = videoWidth;
710 if (bbox_rb.y > videoHeight)
711 bbox_rb.y = videoHeight;
712
713 // calc top bar
714 int box_wd = bbox_rb.x - bbox_lt.x;
715 int box_ht = trackingBox.y() - bbox_lt.y;
716 int pix_cnt = 0;
717 if (box_wd > 0 && box_ht > 0)
718 {
719 pix_cnt += box_wd * box_ht;
720 for (j = bbox_lt.y; j < trackingBox.y(); ++j)
721 {
722 offset = j * videoWidth;
723 for (i = bbox_lt.x; i < bbox_rb.x; ++i)
724 {
725 pptr = pdata + offset + i;
726 threshold += *pptr;
727 }
728 }
729 }
730 // calc left bar
731 box_wd = trackingBox.x() - bbox_lt.x;
732 box_ht = trackingBox.width();
733 if (box_wd > 0 && box_ht > 0)
734 {
735 pix_cnt += box_wd * box_ht;
736 for (j = trackingBox.y(); j < trackingBox.y() + box_ht; ++j)
737 {
738 offset = j * videoWidth;
739 for (i = bbox_lt.x; i < trackingBox.x(); ++i)
740 {
741 pptr = pdata + offset + i;
742 threshold += *pptr;
743 }
744 }
745 }
746 // calc right bar
747 box_wd = bbox_rb.x - trackingBox.x() - trackingBox.width();
748 box_ht = trackingBox.width();
749 if (box_wd > 0 && box_ht > 0)
750 {
751 pix_cnt += box_wd * box_ht;
752 for (j = trackingBox.y(); j < trackingBox.y() + box_ht; ++j)
753 {
754 offset = j * videoWidth;
755 for (i = trackingBox.x() + trackingBox.width(); i < bbox_rb.x; ++i)
756 {
757 pptr = pdata + offset + i;
758 threshold += *pptr;
759 }
760 }
761 }
762 // calc bottom bar
763 box_wd = bbox_rb.x - bbox_lt.x;
764 box_ht = bbox_rb.y - trackingBox.y() - trackingBox.width();
765 if (box_wd > 0 && box_ht > 0)
766 {
767 pix_cnt += box_wd * box_ht;
768 for (j = trackingBox.y() + trackingBox.width(); j < bbox_rb.y; ++j)
769 {
770 offset = j * videoWidth;
771 for (i = bbox_lt.x; i < bbox_rb.x; ++i)
772 {
773 pptr = pdata + offset + i;
774 threshold += *pptr;
775 }
776 }
777 }
778 // find maximum
779 double max_val = 0;
780 for (j = 0; j < trackingBox.width(); ++j)
781 {
782 for (i = 0; i < trackingBox.width(); ++i)
783 {
784 pptr = psrc + i;
785 if (*pptr > max_val)
786 max_val = *pptr;
787 }
788 psrc += videoWidth;
789 }
790 if (pix_cnt != 0)
791 threshold /= (double)pix_cnt;
792
793 // cut by 10% higher then average threshold
794 if (max_val > threshold)
795 threshold += (max_val - threshold) * SMART_CUT_FACTOR;
796
797 //log_i("smart thr. = %f cnt = %d", threshold, pix_cnt);
798 break;
799 }
800 // simple adaptive threshold
801 case AUTO_THRESHOLD:
802 {
803 for (j = 0; j < trackingBox.width(); ++j)
804 {
805 for (i = 0; i < trackingBox.width(); ++i)
806 {
807 pptr = psrc + i;
808 threshold += *pptr;
809 }
810 psrc += videoWidth;
811 }
812 threshold /= square_square;
813 break;
814 }
815 // no threshold subtracion
816 default:
817 {
818 }
819 }
820
821 psrc = porigin;
822 for (j = 0; j < trackingBox.width(); ++j)
823 {
824 for (i = 0; i < trackingBox.width(); ++i)
825 {
826 pptr = psrc + i;
827 pval = *pptr - threshold;
828 pval = pval < 0 ? 0 : pval;
829
830 resx += (double)i * pval;
831 resy += (double)j * pval;
832
833 mass += pval;
834 }
835 psrc += videoWidth;
836 }
837
838 if (mass == 0)
839 mass = 1;
840
841 resx /= mass;
842 resy /= mass;
843
844 ret = GuiderUtils::Vector(trackingBox.x(), trackingBox.y(), 0) + GuiderUtils::Vector(resx, resy, 0);
845
846 return ret;
847 }
848
findLocalStarPosition(QSharedPointer<FITSData> & imageData,const int algorithmIndex,const int videoWidth,const int videoHeight,const QRect & trackingBox)849 GuiderUtils::Vector GuideAlgorithms::findLocalStarPosition(QSharedPointer<FITSData> &imageData,
850 const int algorithmIndex,
851 const int videoWidth,
852 const int videoHeight,
853 const QRect &trackingBox)
854 {
855 switch (imageData->dataType())
856 {
857 case TBYTE:
858 return findLocalStarPosition<uint8_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
859
860 case TSHORT:
861 return findLocalStarPosition<int16_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
862
863 case TUSHORT:
864 return findLocalStarPosition<uint16_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
865
866 case TLONG:
867 return findLocalStarPosition<int32_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
868
869 case TULONG:
870 return findLocalStarPosition<uint32_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
871
872 case TFLOAT:
873 return findLocalStarPosition<float>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
874
875 case TLONGLONG:
876 return findLocalStarPosition<int64_t>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
877
878 case TDOUBLE:
879 return findLocalStarPosition<double>(imageData, algorithmIndex, videoWidth, videoHeight, trackingBox);
880
881 default:
882 break;
883 }
884
885 return GuiderUtils::Vector(-1, -1, -1);
886 }
887