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