Kstars

guidestars.cpp
1 /*
2  SPDX-FileCopyrightText: 2020 Hy Murveit <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "guidestars.h"
8 
9 #include "ekos_guide_debug.h"
10 #include "../guideview.h"
11 #include "Options.h"
12 
13 #include <math.h>
14 #include <stellarsolver.h>
15 #include "ekos/auxiliary/stellarsolverprofileeditor.h"
16 #include <QTime>
17 
18 // Then when looking for the guide star, gets this many candidates.
19 #define STARS_TO_SEARCH 250
20 
21 // Don't use stars with SNR lower than this when computing multi-star drift.
22 #define MIN_DRIFT_SNR 8
23 
24 // If there are fewer than this many stars, don't use this multi-star algorithm.
25 // It will instead back-off to a reticle-based algorithm.
26 #define MIN_STAR_CORRESPONDENCE_SIZE 5
27 
28 // We limit the HFR for guide stars. When searching for the guide star, we relax this by the
29 // margin below (e.g. if a guide star was selected that was near the max guide-star hfr, the later
30 // the hfr increased a little, we still want to be able to find it.
31 constexpr double HFR_MARGIN = 2.0;
32 /*
33  Start with a set of reference (x,y) positions from stars, where one is designated a guide star.
34  Given these and a set of new input stars, determine a mapping of new stars to the references.
35  As the star positions should mostly change via translation, the correspondences are determined by
36  similar geometry. It is possible, though that there is noise in positions, and that
37  some reference stars may not be in the new star group, or that the new stars may contain some
38  stars not appearing in references. However, the guide star must appear in both for this to
39  be successful.
40  */
41 
42 namespace
43 {
44 
45 QString logHeader(const QString &label)
46 {
47  return QString("%1 %2 %3 %4 %5 %6 %7")
48  .arg(label, -9).arg("#", 3).arg("x ", 6).arg("y ", 6).arg("flux", 6)
49  .arg("HFR", 6).arg("SNR ", 5);
50 }
51 
52 QString logStar(const QString &label, int index, const SkyBackground &bg, const Edge &star)
53 {
54  double snr = bg.SNR(star.sum, star.numPixels);
55  return QString("%1 %2 %3 %4 %5 %6 %7")
56  .arg(label, -9)
57  .arg(index, 3)
58  .arg(star.x, 6, 'f', 1)
59  .arg(star.y, 6, 'f', 1)
60  .arg(star.sum, 6, 'f', 0)
61  .arg(star.HFR, 6, 'f', 2)
62  .arg(snr, 5, 'f', 1);
63 }
64 
65 void logStars(const QString &label, const QString &label2, const SkyBackground &bg,
66  int size,
67  std::function<const Edge & (int index)> stars,
68  const QString &extraLabel,
69  std::function<QString (int index)> extras)
70 {
71  QString header = logHeader(label);
72  if (extraLabel.size() > 0)
73  header.append(QString(" %1").arg(extraLabel));
74  qCDebug(KSTARS_EKOS_GUIDE) << header;
75  for (int i = 0; i < size; ++i)
76  {
77  const auto &star = stars(i);
78  QString line = logStar(label2, i, bg, star);
79  if (extraLabel.size() > 0)
80  line.append(QString(" %1").arg(extras(i)));
81  qCDebug(KSTARS_EKOS_GUIDE) << line;
82  }
83 }
84 } //namespace
85 
86 GuideStars::GuideStars()
87 {
88 }
89 
90 // It's possible that we don't map all the stars, if there are too many.
91 int GuideStars::getStarMap(int index)
92 {
93  if (index >= starMap.size() || (index < 0))
94  return -1;
95  return starMap[index];
96 }
97 
98 void GuideStars::setupStarCorrespondence(const QList<Edge> &neighbors, int guideIndex)
99 {
100  qCDebug(KSTARS_EKOS_GUIDE) << "setupStarCorrespondence: neighbors " << neighbors.size() << "guide index" << guideIndex;
101  if (neighbors.size() >= MIN_STAR_CORRESPONDENCE_SIZE)
102  {
103  starMap.clear();
104  for (int i = 0; i < neighbors.size(); ++i)
105  {
106  qCDebug(KSTARS_EKOS_GUIDE) << " adding ref: " << neighbors[i].x << neighbors[i].y;
107  // We need to initialize starMap, because findGuideStar(), which normally
108  // initializes starMap(), might call selectGuideStar() if it finds that
109  // the starCorrespondence is empty.
110  starMap.push_back(i);
111  }
112  starCorrespondence.initialize(neighbors, guideIndex);
113  }
114  else
115  starCorrespondence.reset();
116 }
117 
118 // Calls SEP to generate a set of star detections and score them,
119 // then calls selectGuideStars(detections, scores) to select the guide star.
120 QVector3D GuideStars::selectGuideStar(const QSharedPointer<FITSData> &imageData)
121 {
122  if (imageData == nullptr)
123  return QVector3D(-1, -1, -1);
124 
125  QList<double> sepScores;
126  QList<double> minDistances;
127  const double maxHFR = Options::guideMaxHFR();
128  findTopStars(imageData, Options::maxMultistarReferenceStars(), &detectedStars, maxHFR,
129  nullptr, &sepScores, &minDistances);
130 
131  int maxX = imageData->width();
132  int maxY = imageData->height();
133  return selectGuideStar(detectedStars, sepScores, maxX, maxY, minDistances);
134 }
135 
136 // This selects the guide star by using the input scores, and
137 // downgrading stars too close to the edges of the image (so that calibration won't fail).
138 // It also removes from consideration stars with huge HFR values (most likely not stars).
139 QVector3D GuideStars::selectGuideStar(const QList<Edge> &stars,
140  const QList<double> &sepScores,
141  int maxX, int maxY,
142  const QList<double> &minDistances)
143 {
144  constexpr int maxStarDiameter = 32;
145 
146  if ((uint) stars.size() < Options::minDetectionsSEPMultistar())
147  {
148  qCDebug(KSTARS_EKOS_GUIDE) << "Too few stars detected: " << stars.size();
149  return QVector3D(-1, -1, -1);
150  }
151 
152  int maxIndex = Options::maxMultistarReferenceStars() < stars.count() ? Options::maxMultistarReferenceStars() :
153  stars.count();
154  QVector<int> scores(Options::maxMultistarReferenceStars());
155  QList<Edge> guideStarNeighbors;
156  for (int i = 0; i < maxIndex; i++)
157  {
158  int score = 100 + sepScores[i];
159  const Edge &center = stars.at(i);
160  guideStarNeighbors.append(center);
161 
162  // Severely reject stars close to edges
163  // Worry about calibration? Make calibration distance parameter?
164  constexpr int borderGuard = 35;
165  if (center.x < borderGuard || center.y < borderGuard ||
166  center.x > (maxX - borderGuard) || center.y > (maxY - borderGuard))
167  score -= 1000;
168  // Reject stars bigger than square
169  if (center.HFR > float(maxStarDiameter) / 2)
170  score -= 1000;
171 
172  // Add advantage to stars with SNRs between 40-100.
173  auto bg = skybackground();
174  double snr = bg.SNR(center.sum, center.numPixels);
175  if (snr >= 40 && snr <= 100)
176  score += 75;
177  if (snr >= 100)
178  score -= 50;
179 
180  // discourage stars that have close neighbors.
181  // This isn't perfect, as we're not detecting all stars, just top 100 or so.
182  const double neighborDistance = minDistances[i];
183  constexpr double closeNeighbor = 25;
184  constexpr double veryCloseNeighbor = 15;
185  if (neighborDistance < veryCloseNeighbor)
186  score -= 100;
187  else if (neighborDistance < closeNeighbor)
188  score -= 50;
189 
190  scores[i] = score;
191  }
192 
193  logStars("Select", "Star", skyBackground,
194  maxIndex,
195  [&](int i) -> const Edge&
196  {
197  return stars[i];
198  },
199  "score", [&](int i) -> QString
200  {
201  return QString("%1").arg(scores[i], 5);
202  });
203 
204  int maxScore = -1;
205  int maxScoreIndex = -1;
206  for (int i = 0; i < maxIndex; i++)
207  {
208  if (scores[i] > maxScore)
209  {
210  maxScore = scores[i];
211  maxScoreIndex = i;
212  }
213  }
214 
215  if (maxScoreIndex < 0)
216  {
217  qCDebug(KSTARS_EKOS_GUIDE) << "No suitable star detected.";
218  return QVector3D(-1, -1, -1);
219  }
220  setupStarCorrespondence(guideStarNeighbors, maxScoreIndex);
221  QVector3D newStarCenter(stars[maxScoreIndex].x, stars[maxScoreIndex].y, 0);
222  qCDebug(KSTARS_EKOS_GUIDE) << "new star center: " << maxScoreIndex << " x: "
223  << stars[maxScoreIndex].x << " y: " << stars[maxScoreIndex].y;
224  return newStarCenter;
225 }
226 
227 // Find the current target positions for the guide-star neighbors, and add them
228 // to the guideView.
229 void GuideStars::plotStars(QSharedPointer<GuideView> &guideView, const QRect &trackingBox)
230 {
231  if (guideView == nullptr) return;
232  guideView->clearNeighbors();
233  if (starCorrespondence.size() == 0) return;
234 
235  const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
236 
237  // Find the offset between the current reticle position and the original
238  // guide star position (e.g. it might have moved due to dithering).
239  double reticle_x = trackingBox.x() + trackingBox.width() / 2.0;
240  double reticle_y = trackingBox.y() + trackingBox.height() / 2.0;
241 
242  // See which neighbor stars were associated with detected stars.
243  QVector<bool> found(starCorrespondence.size(), false);
244  QVector<int> detectedIndex(starCorrespondence.size(), -1);
245 
246  for (int i = 0; i < detectedStars.size(); ++i)
247  {
248  int refIndex = getStarMap(i);
249  if (refIndex >= 0)
250  {
251  found[refIndex] = true;
252  detectedIndex[refIndex] = i;
253  }
254  }
255  // Add to the guideView the neighbor stars' target positions and whether they were found.
256  for (int i = 0; i < starCorrespondence.size(); ++i)
257  {
258  bool isGuideStar = i == starCorrespondence.guideStar();
259  const QVector2D offset = starCorrespondence.offset(i);
260  const double detected_x = found[i] ? detectedStars[detectedIndex[i]].x : 0;
261  const double detected_y = found[i] ? detectedStars[detectedIndex[i]].y : 0;
262  guideView->addGuideStarNeighbor(offset.x() + reticle_x, offset.y() + reticle_y, found[i],
263  detected_x, detected_y, isGuideStar);
264  }
265  guideView->updateNeighbors();
266 }
267 
268 // Find the guide star using the starCorrespondence algorithm (looking for
269 // the other reference stars in the same relative position as when the guide star was selected).
270 // If this method fails, it backs off to looking in the tracking box for the highest scoring star.
271 GuiderUtils::Vector GuideStars::findGuideStar(const QSharedPointer<FITSData> &imageData, const QRect &trackingBox,
272  QSharedPointer<GuideView> &guideView, bool firstFrame)
273 {
274  // We fall-back to single star guiding if multistar isn't working,
275  // but limit this for a number of iterations.
276  // If a user doesn't have multiple stars available, the user shouldn't be using multistar.
277  constexpr int MAX_CONSECUTIVE_UNRELIABLE = 10;
278  if (firstFrame)
279  unreliableDectionCounter = 0;
280 
281  // Don't accept reference stars whose position is more than this many pixels from expected.
282  constexpr double maxStarAssociationDistance = 10;
283 
284  if (imageData == nullptr)
285  return GuiderUtils::Vector(-1, -1, -1);
286 
287  // If the guide star has not yet been set up, then establish it here.
288  // Not thrilled doing this, but this is the way the internal guider is setup
289  // when guiding is restarted by the scheduler (the normal establish guide star
290  // methods are not called).
291  if (firstFrame && starCorrespondence.size() == 0)
292  {
293  QVector3D v = selectGuideStar(imageData);
294  qCDebug(KSTARS_EKOS_GUIDE) << QString("findGuideStar: Called without starCorrespondence. Refound guide star at %1 %2")
295  .arg(QString::number(v.x(), 'f', 1)).arg(QString::number(v.y(), 'f', 1));
296  return GuiderUtils::Vector(v.x(), v.y(), v.z());
297  }
298 
299  // Allow a little margin above the max hfr for guide stars when searching for the guide star.
300  const double maxHFR = Options::guideMaxHFR() + HFR_MARGIN;
301  if (starCorrespondence.size() > 0)
302  {
303 
304  findTopStars(imageData, STARS_TO_SEARCH, &detectedStars, maxHFR);
305  if (detectedStars.empty())
306  return GuiderUtils::Vector(-1, -1, -1);
307 
308  // Allow it to guide even if the main guide star isn't detected (as long as enough reference stars are).
309  starCorrespondence.setAllowMissingGuideStar(allowMissingGuideStar);
310 
311  // Star correspondence can run quicker if it knows the image size.
312  starCorrespondence.setImageSize(imageData->width(), imageData->height());
313 
314  // When using large star-correspondence sets and filtering with a StellarSolver profile,
315  // the stars at the edge of detection can be lost. Best not to filter, but...
316  double minFraction = 0.5;
317  if (starCorrespondence.size() > 25) minFraction = 0.33;
318  else if (starCorrespondence.size() > 15) minFraction = 0.4;
319 
320  Edge foundStar = starCorrespondence.find(detectedStars, maxStarAssociationDistance, &starMap, true, minFraction);
321 
322  // Is there a correspondence to the guide star
323  // Should we also weight distance to the tracking box?
324  for (int i = 0; i < detectedStars.size(); ++i)
325  {
326  if (getStarMap(i) == starCorrespondence.guideStar())
327  {
328  auto &star = detectedStars[i];
329  double SNR = skyBackground.SNR(star.sum, star.numPixels);
330  guideStarSNR = SNR;
331  guideStarMass = star.sum;
332  unreliableDectionCounter = 0;
333  qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence found " << i << "at" << star.x << star.y << "SNR" << SNR;
334  if (guideView != nullptr)
335  plotStars(guideView, trackingBox);
336  return GuiderUtils::Vector(star.x, star.y, 0);
337  }
338  }
339  // None of the stars matched the guide star, but it's possible star correspondence
340  // invented a guide star position.
341  if (foundStar.x >= 0 && foundStar.y >= 0)
342  {
343  guideStarSNR = skyBackground.SNR(foundStar.sum, foundStar.numPixels);
344  guideStarMass = foundStar.sum;
345  unreliableDectionCounter = 0; // debating this
346  qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence invented at" << foundStar.x << foundStar.y << "SNR" << guideStarSNR;
347  if (guideView != nullptr)
348  plotStars(guideView, trackingBox);
349  return GuiderUtils::Vector(foundStar.x, foundStar.y, 0);
350  }
351  }
352 
353  qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence not used. It failed to find the guide star.";
354 
355  if (++unreliableDectionCounter > MAX_CONSECUTIVE_UNRELIABLE)
356  return GuiderUtils::Vector(-1, -1, -1);
357 
358  logStars("findGuide", "Star", skyBackground,
359  detectedStars.size(),
360  [&](int i) -> const Edge&
361  {
362  return detectedStars[i];
363  },
364  "assoc", [&](int i) -> QString
365  {
366  return QString("%1").arg(getStarMap(i));
367  });
368 
369  if (trackingBox.isValid() == false)
370  return GuiderUtils::Vector(-1, -1, -1);
371 
372  // If we didn't find a star that way, then fall back
373  findTopStars(imageData, 100, &detectedStars, maxHFR, &trackingBox);
374  if (detectedStars.size() > 0)
375  {
376  // Find the star closest to the guide star position, if we have a position.
377  // Otherwise use the center of the tracking box (which must be valid, see above if).
378  // 1. Get the guide star position
379  int best = 0;
380  double refX = trackingBox.x() + trackingBox.width() / 2;
381  double refY = trackingBox.y() + trackingBox.height() / 2;
382  if (starCorrespondence.size() > 0 && starCorrespondence.guideStar() >= 0)
383  {
384  const auto gStar = starCorrespondence.reference(starCorrespondence.guideStar());
385  refX = gStar.x;
386  refY = gStar.y;
387  }
388  // 2. Find the closest star to that position.
389  double minDistSq = 1e8;
390  for (int i = 0; i < detectedStars.size(); ++i)
391  {
392  const auto &dStar = detectedStars[i];
393  const double distSq = (dStar.x - refX) * (dStar.x - refX) + (dStar.y - refY) * (dStar.y - refY);
394  if (distSq < minDistSq)
395  {
396  best = i;
397  minDistSq = distSq;
398  }
399  }
400  // 3. Return that star.
401  auto &star = detectedStars[best];
402  double SNR = skyBackground.SNR(star.sum, star.numPixels);
403  guideStarSNR = SNR;
404  guideStarMass = star.sum;
405  qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence. Standard method found at " << star.x << star.y << "SNR" << SNR;
406  return GuiderUtils::Vector(star.x, star.y, 0);
407  }
408  return GuiderUtils::Vector(-1, -1, -1);
409 }
410 
411 SSolver::Parameters GuideStars::getStarExtractionParameters(int num)
412 {
413  SSolver::Parameters params;
414  params.listName = "Guider";
415  params.apertureShape = SSolver::SHAPE_CIRCLE;
416  params.minarea = 10; // changed from 5 --> 10
417  params.deblend_thresh = 32;
418  params.deblend_contrast = 0.005;
419  params.initialKeep = num * 2;
420  params.keepNum = num;
421  params.removeBrightest = 0;
422  params.removeDimmest = 0;
423  params.saturationLimit = 0;
424  return params;
425 }
426 
427 // This is the interface to star detection.
428 int GuideStars::findAllSEPStars(const QSharedPointer<FITSData> &imageData, QList<Edge *> *sepStars, int num)
429 {
430  if (imageData == nullptr)
431  return 0;
432 
433  QVariantMap settings;
434  settings["optionsProfileIndex"] = Options::guideOptionsProfile();
435  settings["optionsProfileGroup"] = static_cast<int>(Ekos::GuideProfiles);
436  imageData->setSourceExtractorSettings(settings);
437  imageData->findStars(ALGORITHM_SEP).waitForFinished();
438  skyBackground = imageData->getSkyBackground();
439 
440  QList<Edge *> edges = imageData->getStarCenters();
441  // Let's sort edges, starting with widest
442  std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->HFR > edge2->HFR;});
443 
444  // Take only the first num stars
445  {
446  int starCount = qMin(num, edges.count());
447  for (int i = 0; i < starCount; i++)
448  sepStars->append(edges[i]);
449  }
450  qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: SEP detected" << edges.count() << "stars, keeping" << sepStars->size();
451 
452  edges.clear();
453  return sepStars->count();
454 }
455 
456 double GuideStars::findMinDistance(int index, const QList<Edge*> &stars)
457 {
458  double closestDistanceSqr = 1e10;
459  const Edge &star = *stars[index];
460  for (int i = 0; i < stars.size(); ++i)
461  {
462  if (i == index) continue;
463  const Edge &thisStar = *stars[i];
464  const double xDist = star.x - thisStar.x;
465  const double yDist = star.y - thisStar.y;
466  const double distanceSqr = xDist * xDist + yDist * yDist;
467  if (distanceSqr < closestDistanceSqr)
468  closestDistanceSqr = distanceSqr;
469  }
470  return sqrt(closestDistanceSqr);
471 }
472 
473 // Returns a list of 'num' stars, sorted according to evaluateSEPStars().
474 // If the region-of-interest rectange is not null, it only returns scores in that area.
475 void GuideStars::findTopStars(const QSharedPointer<FITSData> &imageData, int num, QList<Edge> *stars,
476  const double maxHFR, const QRect *roi,
477  QList<double> *outputScores, QList<double> *minDistances)
478 {
479  if (roi == nullptr)
480  qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num;
481  else
482  qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num <<
483  QString("Roi(%1,%2 %3x%4)").arg(roi->x()).arg(roi->y()).arg(roi->width()).arg(roi->height());
484 
485  stars->clear();
486  if (imageData == nullptr)
487  return;
488 
489  QElapsedTimer timer;
490  timer.restart();
491  QList<Edge*> sepStars;
492  int count = findAllSEPStars(imageData, &sepStars, num * 2);
493  if (count == 0)
494  return;
495 
496  if (sepStars.empty())
497  return;
498 
499  QVector<double> scores;
500  evaluateSEPStars(sepStars, &scores, roi, maxHFR);
501  // Sort the sepStars by score, higher score to lower score.
503  for (int i = 0; i < scores.size(); ++i)
504  sc.push_back(std::pair<int, double>(i, scores[i]));
505  std::sort(sc.begin(), sc.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
506  {
507  return a.second > b.second;
508  });
509  // Copy the top num results.
510  for (int i = 0; i < std::min(num, scores.size()); ++i)
511  {
512  const int starIndex = sc[i].first;
513  const double starScore = sc[i].second;
514  if (starScore >= 0)
515  {
516  stars->append(*(sepStars[starIndex]));
517  if (outputScores != nullptr)
518  outputScores->append(starScore);
519  if (minDistances != nullptr)
520  minDistances->append(findMinDistance(starIndex, sepStars));
521  }
522  }
523  qCDebug(KSTARS_EKOS_GUIDE)
524  << QString("Multistar: findTopStars returning: %1 stars, %2s")
525  .arg(stars->size()).arg(timer.elapsed() / 1000.0, 4, 'f', 2);
526 }
527 
528 // Scores star detection relative to each other. Uses the star's SNR as the main measure.
529 void GuideStars::evaluateSEPStars(const QList<Edge *> &starCenters, QVector<double> *scores,
530  const QRect *roi, const double maxHFR) const
531 {
532  auto centers = starCenters;
533  scores->clear();
534  const int numDetections = centers.size();
535  for (int i = 0; i < numDetections; ++i) scores->push_back(0);
536  if (numDetections == 0) return;
537 
538 
539  // Rough constants used by this weighting.
540  // If the center pixel is above this, assume it's clipped and don't emphasize.
541  constexpr double snrWeight = 20; // Measure weight if/when multiple measures are used.
542  auto bg = skybackground();
543 
544  // Sort by SNR in increasing order so the weighting goes up.
545  // Assign score based on the sorted position.
546  std::sort(centers.begin(), centers.end(), [&bg](const Edge * a, const Edge * b)
547  {
548  double snrA = bg.SNR(a->sum, a->numPixels);
549  double snrB = bg.SNR(b->sum, b->numPixels);
550  return snrA < snrB;
551  });
552 
553  // See if the HFR restriction is too severe.
554  int numRejectedByHFR = 0;
555  for (int i = 0; i < numDetections; ++i)
556  {
557  if (centers.at(i)->HFR > maxHFR)
558  numRejectedByHFR++;
559  }
560  const int starsRemaining = numDetections - numRejectedByHFR;
561  const bool useHFRConstraint =
562  starsRemaining > 5 ||
563  (starsRemaining >= 3 && numDetections <= 6) ||
564  (starsRemaining >= 2 && numDetections <= 4) ||
565  (starsRemaining >= 1 && numDetections <= 2);
566 
567  for (int i = 0; i < numDetections; ++i)
568  {
569  for (int j = 0; j < numDetections; ++j)
570  {
571  if (starCenters.at(j) == centers.at(i))
572  {
573  // Don't emphasize stars that are too wide.
574  if (useHFRConstraint && centers.at(i)->HFR > maxHFR)
575  (*scores)[j] = -1;
576  else
577  (*scores)[j] += snrWeight * i;
578  break;
579  }
580  }
581  }
582 
583  // If we are insisting on a star in the tracking box.
584  if (roi != nullptr)
585  {
586  // Reset the score of anything outside the tracking box.
587  for (int j = 0; j < starCenters.size(); ++j)
588  {
589  const auto &star = starCenters.at(j);
590  if (star->x < roi->x() || star->x >= roi->x() + roi->width() ||
591  star->y < roi->y() || star->y >= roi->y() + roi->height())
592  (*scores)[j] = -1;
593  }
594  }
595 }
596 
597 // Given a star detection and a reference star compute the RA and DEC drifts between the two.
598 void GuideStars::computeStarDrift(const Edge &star, const Edge &reference, double *driftRA, double *driftDEC) const
599 {
600  if (!calibrationInitialized) return;
601 
602  GuiderUtils::Vector position(star.x, star.y, 0);
603  GuiderUtils::Vector reference_position(reference.x, reference.y, 0);
604  GuiderUtils::Vector arc_position, arc_reference_position;
605  arc_position = calibration.convertToArcseconds(position);
606  arc_reference_position = calibration.convertToArcseconds(reference_position);
607 
608  // translate into sky coords.
609  GuiderUtils::Vector sky_coords = arc_position - arc_reference_position;
610  sky_coords = calibration.rotateToRaDec(sky_coords);
611 
612  // Save the drifts in RA and DEC.
613  *driftRA = sky_coords.x;
614  *driftDEC = sky_coords.y;
615 }
616 
617 // Computes the overall drift given the input image that was analyzed in findGuideStar().
618 // Returns true if this can be done--there are already guide stars and the current detections
619 // can be compared to them. In that case, fills in RADRift and DECDrift with arc-second drifts.
620 // Reticle_x,y is the target position. It may be different from when we
621 // originally selected the guide star, e.g. due to dithering, etc. We compute an offset from the
622 // original guide-star position and the reticle position and offset all the reference star
623 // positions by that.
624 bool GuideStars::getDrift(double oneStarDrift, double reticle_x, double reticle_y,
625  double *RADrift, double *DECDrift)
626 {
627  if (starCorrespondence.size() == 0)
628  return false;
629  const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
630  qCDebug(KSTARS_EKOS_GUIDE) << "Multistar getDrift, reticle:" << reticle_x << reticle_y
631  << "guidestar" << gStar.x << gStar.y
632  << "so offsets:" << (reticle_x - gStar.x) << (reticle_y - gStar.y);
633  // Revoke multistar if we're that far away.
634  constexpr double maxDriftForMultistar = 4.0; // arc-seconds
635 
636  if (!calibrationInitialized)
637  return false;
638 
639  // Don't run multistar if the 1-star approach is too far off.
640  if (oneStarDrift > maxDriftForMultistar)
641  {
642  qCDebug(KSTARS_EKOS_GUIDE) << "MultiStar: 1-star drift too high: " << oneStarDrift;
643  return false;
644  }
645  if (starCorrespondence.size() < 2)
646  {
647  qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: no guide stars";
648  return false;
649  }
650 
651  // The original guide star position has been offset.
652  double offset_x = reticle_x - gStar.x;
653  double offset_y = reticle_y - gStar.y;
654 
655  int numStarsProcessed = 0;
656  bool guideStarProcessed = false;
657  double driftRASum = 0, driftDECSum = 0;
658  double guideStarRADrift = 0, guideStarDECDrift = 0;
659  QVector<double> raDrifts, decDrifts;
660  qCDebug(KSTARS_EKOS_GUIDE)
661  << QString("%1 %2 dRA dDEC").arg(logHeader("")).arg(logHeader(" Ref:"));
662  for (int i = 0; i < detectedStars.size(); ++i)
663  {
664  const auto &star = detectedStars[i];
665  auto bg = skybackground();
666  double snr = bg.SNR(detectedStars[i].sum, detectedStars[i].numPixels);
667  // Probably should test SNR on the reference as well.
668  if (getStarMap(i) >= 0 && snr >= MIN_DRIFT_SNR)
669  {
670  auto ref = starCorrespondence.reference(getStarMap(i));
671  ref.x += offset_x;
672  ref.y += offset_y;
673  double driftRA, driftDEC;
674  computeStarDrift(star, ref, &driftRA, &driftDEC);
675  if (getStarMap(i) == starCorrespondence.guideStar())
676  {
677  guideStarRADrift = driftRA;
678  guideStarDECDrift = driftDEC;
679  guideStarProcessed = true;
680  }
681 
682  raDrifts.push_back(driftRA);
683  decDrifts.push_back(driftDEC);
684  numStarsProcessed++;
685 
686  qCDebug(KSTARS_EKOS_GUIDE)
687  << QString("%1 %2 %3 %4").arg(logStar("MultiStar", i, bg, star))
688  .arg(logStar(" Ref:", getStarMap(i), bg, ref))
689  .arg(driftRA, 5, 'f', 2).arg(driftDEC, 5, 'f', 2);
690  }
691  }
692 
693  if (numStarsProcessed == 0 || (!allowMissingGuideStar && !guideStarProcessed))
694  return false;
695 
696  // Filter out reference star drifts that are too different from the guide-star drift.
697  QVector<double> raDriftsKeep, decDriftsKeep;
698  for (int i = 0; i < raDrifts.size(); ++i)
699  {
700  if ((allowMissingGuideStar && !guideStarProcessed) ||
701  ((fabs(raDrifts[i] - guideStarRADrift) < 2.0) &&
702  (fabs(decDrifts[i] - guideStarDECDrift) < 2.0)))
703  {
704  driftRASum += raDrifts[i];
705  driftDECSum += decDrifts[i];
706  raDriftsKeep.push_back(raDrifts[i]);
707  decDriftsKeep.push_back(decDrifts[i]);
708  }
709  }
710  if (raDriftsKeep.empty() || raDriftsKeep.empty())
711  return false; // Shouldn't happen.
712 
713  numStarsProcessed = raDriftsKeep.size();
714 
715  // Generate the drift either from the median or the average of the usable reference drifts.
716  bool useMedian = true;
717  if (useMedian)
718  {
719  std::sort(raDriftsKeep.begin(), raDriftsKeep.end(), [] (double d1, double d2)
720  {
721  return d1 < d2;
722  });
723  std::sort(decDriftsKeep.begin(), decDriftsKeep.end(), [] (double d1, double d2)
724  {
725  return d1 < d2;
726  });
727  if (numStarsProcessed % 2 == 0)
728  {
729  const int middle = numStarsProcessed / 2;
730  *RADrift = ( raDriftsKeep[middle - 1] + raDriftsKeep[middle]) / 2.0;
731  *DECDrift = (decDriftsKeep[middle - 1] + decDriftsKeep[middle]) / 2.0;
732  }
733  else
734  {
735  const int middle = numStarsProcessed / 2; // because it's 0-based
736  *RADrift = raDriftsKeep[middle];
737  *DECDrift = decDriftsKeep[middle];
738  }
739  qCDebug(KSTARS_EKOS_GUIDE) << "MultiStar: Drift median " << *RADrift << *DECDrift << numStarsProcessed << " of " <<
740  detectedStars.size() << "#guide" << starCorrespondence.size();
741  }
742  else
743  {
744  *RADrift = driftRASum / raDriftsKeep.size();
745  *DECDrift = driftDECSum / decDriftsKeep.size();
746  qCDebug(KSTARS_EKOS_GUIDE) << "MultiStar: Drift " << *RADrift << *DECDrift << numStarsProcessed << " of " <<
747  detectedStars.size() << "#guide" << starCorrespondence.size();
748  }
749  return true;
750 }
751 
752 void GuideStars::setCalibration(const Calibration &cal)
753 {
754  calibration = cal;
755  calibrationInitialized = true;
756 }
void append(const T &value)
QString number(int n, int base)
int size() const const
QVector::iterator begin()
int count(const T &value) const const
void ref()
void push_back(const T &value)
int width() const const
int x() const const
int y() const const
void push_back(const T &value)
T & first()
bool empty() const const
int size() const const
void clear()
qint64 restart()
bool isValid() const const
const T & at(int i) const const
qint64 elapsed() const const
QVector::iterator end()
int height() const const
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
void clear()
QList::iterator begin()
float x() const const
float y() const const
float z() const const
QTextStream & center(QTextStream &stream)
int size() const const
float x() const const
float y() const const
QList::iterator end()
QString & append(QChar ch)
bool empty() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Mon Aug 15 2022 04:04:01 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.