Kstars

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

KDE's Doxygen guidelines are available online.