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 <[email protected]>
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;
39  Peak(int x_, int y_, float val_) : x(x_), y(y_), val(val_) { }
40  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 
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 
76  [email protected], B2, C1, C3, D1
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 
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 
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 
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
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>
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(-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 
846 GuiderUtils::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 }
void append(const T &value)
bool isNull() const const
int width() const const
int x() const const
int y() const const
int left() const const
QAction * repeat(const QObject *recvr, const char *slot, QObject *parent)
int top() const const
bool isValid() const const
void waitForFinished()
int height() const const
QTextStream & center(QTextStream &stream)
QAction * next(const QObject *recvr, const char *slot, QObject *parent)
const QList< QKeySequence > & end()
T result() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Tue Aug 9 2022 04:06:03 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.