9#include "ekos_guide_debug.h"
10#include "../guideview.h"
11#include "fitsviewer/fitsdata.h"
12#include "fitsviewer/fitssepdetector.h"
16#include <stellarsolver.h>
17#include "ekos/auxiliary/stellarsolverprofileeditor.h"
20#define DLOG if (false) qCDebug
23#define STARS_TO_SEARCH 250
26#define MIN_DRIFT_SNR 8
30#define MIN_STAR_CORRESPONDENCE_SIZE 5
35constexpr double HFR_MARGIN = 2.0;
51 return QString(
"%1 %2 %3 %4 %5 %6 %7")
53 .
arg(
"HFR", 6).
arg(
"SNR ", 5);
56QString logStar(
const QString &label,
int index,
const SkyBackground &bg,
const Edge &star)
58 double snr = bg.SNR(star.sum, star.numPixels);
59 return QString(
"%1 %2 %3 %4 %5 %6 %7")
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)
69void logStars(
const QString &label,
const QString &label2,
const SkyBackground &bg,
71 std::function<
const Edge & (
int index)> stars,
73 std::function<
QString (
int index)> extras)
75 QString header = logHeader(label);
76 if (extraLabel.
size() > 0)
78 qCDebug(KSTARS_EKOS_GUIDE) << header;
79 for (
int i = 0; i < size; ++i)
81 const auto &star = stars(i);
82 QString line = logStar(label2, i, bg, star);
83 if (extraLabel.
size() > 0)
85 qCDebug(KSTARS_EKOS_GUIDE) << line;
90GuideStars::GuideStars()
95int GuideStars::getStarMap(
int index)
97 if (index >= starMap.size() || (index < 0))
99 return starMap[index];
102void GuideStars::setupStarCorrespondence(
const QList<Edge> &neighbors,
int guideIndex)
104 qCDebug(KSTARS_EKOS_GUIDE) <<
"setupStarCorrespondence: neighbors " << neighbors.
size() <<
"guide index" << guideIndex;
105 if (neighbors.
size() >= MIN_STAR_CORRESPONDENCE_SIZE)
108 for (
int i = 0; i < neighbors.
size(); ++i)
110 qCDebug(KSTARS_EKOS_GUIDE) <<
" adding ref: " << neighbors[i].x << neighbors[i].y;
114 starMap.push_back(i);
116 starCorrespondence.initialize(neighbors, guideIndex);
119 starCorrespondence.reset();
126 if (imageData ==
nullptr)
131 const double maxHFR = Options::guideMaxHFR();
132 findTopStars(imageData, Options::maxMultistarReferenceStars(), &detectedStars, maxHFR,
133 nullptr, &sepScores, &minDistances);
135 int maxX = imageData->width();
136 int maxY = imageData->height();
137 return selectGuideStar(detectedStars, sepScores, maxX, maxY, minDistances);
148 constexpr int maxStarDiameter = 32;
150 if ((uint) stars.
size() < Options::minDetectionsSEPMultistar())
152 qCDebug(KSTARS_EKOS_GUIDE) <<
"Too few stars detected: " << stars.
size();
156 int maxIndex = Options::maxMultistarReferenceStars() <
static_cast<uint
>(stars.
count()) ?
157 Options::maxMultistarReferenceStars() :
159 QVector<int> scores(Options::maxMultistarReferenceStars());
163 for (
int i = 0; i < maxIndex; i++)
165 std::sort(hfrs.
begin(), hfrs.
end());
166 float tooLowHFR = 0.0;
167 if (maxIndex / 4 > 0)
168 tooLowHFR = hfrs[maxIndex / 4];
171 for (
int i = 0; i < maxIndex; i++)
173 int score = 100 + sepScores[i];
175 guideStarNeighbors.
append(center);
179 constexpr int borderGuard = 35;
181 center.x > (maxX - borderGuard) ||
center.y > (maxY - borderGuard))
184 if (
center.HFR >
float(maxStarDiameter) / 2)
188 if (
center.HFR < tooLowHFR)
192 auto bg = skybackground();
194 if (snr >= 40 && snr <= 100)
201 const double neighborDistance = minDistances[i];
202 constexpr double closeNeighbor = 25;
203 constexpr double veryCloseNeighbor = 15;
204 if (neighborDistance < veryCloseNeighbor)
206 else if (neighborDistance < closeNeighbor)
212 logStars(
"Select",
"Star", skyBackground,
214 [&](
int i) ->
const Edge&
224 int maxScoreIndex = -1;
225 for (
int i = 0; i < maxIndex; i++)
227 if (scores[i] > maxScore)
229 maxScore = scores[i];
234 if (maxScoreIndex < 0)
236 qCDebug(KSTARS_EKOS_GUIDE) <<
"No suitable star detected.";
239 setupStarCorrespondence(guideStarNeighbors, maxScoreIndex);
240 QVector3D newStarCenter(stars[maxScoreIndex].x, stars[maxScoreIndex].y, 0);
241 qCDebug(KSTARS_EKOS_GUIDE) <<
"new star center: " << maxScoreIndex <<
" x: "
242 << stars[maxScoreIndex].x <<
" y: " << stars[maxScoreIndex].y;
243 return newStarCenter;
250 if (guideView ==
nullptr)
return;
251 guideView->clearNeighbors();
252 if (starCorrespondence.size() == 0)
return;
254 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
258 double reticle_x = trackingBox.
x() + trackingBox.
width() / 2.0;
259 double reticle_y = trackingBox.
y() + trackingBox.
height() / 2.0;
263 QVector<int> detectedIndex(starCorrespondence.size(), -1);
265 for (
int i = 0; i < detectedStars.
size(); ++i)
267 int refIndex = getStarMap(i);
270 found[refIndex] =
true;
271 detectedIndex[refIndex] = i;
275 for (
int i = 0; i < starCorrespondence.size(); ++i)
277 bool isGuideStar = i == starCorrespondence.guideStar();
278 const QVector2D offset = starCorrespondence.offset(i);
279 const double detected_x = found[i] ? detectedStars[detectedIndex[i]].x : 0;
280 const double detected_y = found[i] ? detectedStars[detectedIndex[i]].y : 0;
281 guideView->addGuideStarNeighbor(offset.
x() + reticle_x, offset.
y() + reticle_y, found[i],
282 detected_x, detected_y, isGuideStar);
284 guideView->updateNeighbors();
296 constexpr int MAX_CONSECUTIVE_UNRELIABLE = 10;
298 unreliableDectionCounter = 0;
301 constexpr double maxStarAssociationDistance = 10;
303 if (imageData ==
nullptr)
304 return GuiderUtils::Vector(-1, -1, -1);
310 if (firstFrame && starCorrespondence.size() == 0)
312 QVector3D v = selectGuideStar(imageData);
313 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"findGuideStar: Called without starCorrespondence. Refound guide star at %1 %2")
315 return GuiderUtils::Vector(v.
x(), v.
y(), v.
z());
319 const double maxHFR = Options::guideMaxHFR() + HFR_MARGIN;
320 if (starCorrespondence.size() > 0)
323 findTopStars(imageData, STARS_TO_SEARCH, &detectedStars, maxHFR);
324 if (detectedStars.
empty())
325 return GuiderUtils::Vector(-1, -1, -1);
328 starCorrespondence.setAllowMissingGuideStar(allowMissingGuideStar);
331 starCorrespondence.setImageSize(imageData->width(), imageData->height());
335 double minFraction = 0.5;
336 if (starCorrespondence.size() > 25) minFraction = 0.33;
337 else if (starCorrespondence.size() > 15) minFraction = 0.4;
339 Edge foundStar = starCorrespondence.find(detectedStars, maxStarAssociationDistance, &starMap,
true, minFraction);
343 for (
int i = 0; i < detectedStars.
size(); ++i)
345 if (getStarMap(i) == starCorrespondence.guideStar())
347 auto &star = detectedStars[i];
348 double SNR = skyBackground.SNR(star.sum, star.numPixels);
350 guideStarMass = star.sum;
351 unreliableDectionCounter = 0;
352 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"StarCorrespondence found star %1 at %2 %3 SNR %4")
353 .
arg(i).
arg(star.x, 0,
'f', 1).
arg(star.y, 0,
'f', 1).
arg(SNR, 0,
'f', 1);
355 if (guideView !=
nullptr)
356 plotStars(guideView, trackingBox);
357 return GuiderUtils::Vector(star.x, star.y, 0);
362 if (foundStar.x >= 0 && foundStar.y >= 0)
364 guideStarSNR = skyBackground.SNR(foundStar.sum, foundStar.numPixels);
365 guideStarMass = foundStar.sum;
366 unreliableDectionCounter = 0;
367 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence invented at" << foundStar.x << foundStar.y <<
"SNR" << guideStarSNR;
368 if (guideView !=
nullptr)
369 plotStars(guideView, trackingBox);
370 return GuiderUtils::Vector(foundStar.x, foundStar.y, 0);
374 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence not used. It failed to find the guide star.";
376 if (++unreliableDectionCounter > MAX_CONSECUTIVE_UNRELIABLE)
377 return GuiderUtils::Vector(-1, -1, -1);
379 logStars(
"findGuide",
"Star", skyBackground,
380 detectedStars.
size(),
381 [&](
int i) ->
const Edge&
383 return detectedStars[i];
387 return QString(
"%1").arg(getStarMap(i));
390 if (trackingBox.
isValid() ==
false)
391 return GuiderUtils::Vector(-1, -1, -1);
394 findTopStars(imageData, 100, &detectedStars, maxHFR, &trackingBox);
395 if (detectedStars.
size() > 0)
401 double refX = trackingBox.
x() + trackingBox.
width() / 2;
402 double refY = trackingBox.
y() + trackingBox.
height() / 2;
403 if (starCorrespondence.size() > 0 && starCorrespondence.guideStar() >= 0)
405 const auto gStar = starCorrespondence.reference(starCorrespondence.guideStar());
410 double minDistSq = 1e8;
411 for (
int i = 0; i < detectedStars.
size(); ++i)
413 const auto &dStar = detectedStars[i];
414 const double distSq = (dStar.x - refX) * (dStar.x - refX) + (dStar.y - refY) * (dStar.y - refY);
415 if (distSq < minDistSq)
422 auto &star = detectedStars[best];
423 double SNR = skyBackground.SNR(star.sum, star.numPixels);
425 guideStarMass = star.sum;
426 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence. Standard method found at " << star.x << star.y <<
"SNR" << SNR;
427 return GuiderUtils::Vector(star.x, star.y, 0);
429 return GuiderUtils::Vector(-1, -1, -1);
432SSolver::Parameters GuideStars::getStarExtractionParameters(
int num)
434 SSolver::Parameters params;
435 params.listName =
"Guider";
436 params.apertureShape = SSolver::SHAPE_CIRCLE;
438 params.deblend_thresh = 32;
439 params.deblend_contrast = 0.005;
440 params.initialKeep = num * 2;
441 params.keepNum = num;
442 params.removeBrightest = 0;
443 params.removeDimmest = 0;
444 params.saturationLimit = 0;
451 if (imageData ==
nullptr)
454 QVariantMap settings;
455 settings[
"optionsProfileIndex"] = Options::guideOptionsProfile();
456 settings[
"optionsProfileGroup"] =
static_cast<int>(Ekos::GuideProfiles);
457 imageData->setSourceExtractorSettings(settings);
458 imageData->findStars(ALGORITHM_SEP).waitForFinished();
459 skyBackground = imageData->getSkyBackground();
463 std::sort(edges.
begin(), edges.
end(), [](
const Edge * edge1,
const Edge * edge2) ->
bool { return edge1->HFR > edge2->HFR;});
465 m_NumStarsDetected = edges.
count();
468 int starCount = qMin(num, edges.
count());
469 for (
int i = 0; i < starCount; i++)
470 sepStars->
append(edges[i]);
472 qCDebug(KSTARS_EKOS_GUIDE) <<
"Multistar: SEP detected" << edges.
count() <<
"stars, keeping" << sepStars->
size();
475 return sepStars->
count();
478double GuideStars::findMinDistance(
int index,
const QList<Edge*> &stars)
480 double closestDistanceSqr = 1e10;
481 const Edge &star = *stars[index];
482 for (
int i = 0; i < stars.
size(); ++i)
484 if (i == index)
continue;
485 const Edge &thisStar = *stars[i];
486 const double xDist = star.x - thisStar.x;
487 const double yDist = star.y - thisStar.y;
488 const double distanceSqr = xDist * xDist + yDist * yDist;
489 if (distanceSqr < closestDistanceSqr)
490 closestDistanceSqr = distanceSqr;
492 return sqrt(closestDistanceSqr);
498 const double maxHFR,
const QRect *roi,
503 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar: findTopStars" << num;
507 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar: findTopStars" << num <<
512 if (imageData ==
nullptr)
518 int count = findAllSEPStars(imageData, &sepStars, num * 2);
522 if (sepStars.
empty())
526 evaluateSEPStars(sepStars, &scores, roi, maxHFR);
529 for (
int i = 0; i < scores.
size(); ++i)
530 sc.
push_back(std::pair<int, double>(i, scores[i]));
531 std::sort(sc.
begin(), sc.
end(), [](
const std::pair<int, double> &a,
const std::pair<int, double> &b)
533 return a.second > b.second;
536 for (
int i = 0; i < std::min(num, scores.
size()); ++i)
538 const int starIndex = sc[i].
first;
539 const double starScore = sc[i].second;
542 stars->
append(*(sepStars[starIndex]));
543 if (outputScores !=
nullptr)
544 outputScores->
append(starScore);
545 if (minDistances !=
nullptr)
546 minDistances->
append(findMinDistance(starIndex, sepStars));
549 DLOG(KSTARS_EKOS_GUIDE)
550 <<
QString(
"Multistar: findTopStars returning: %1 stars, %2s")
551 .
arg(stars->
size()).
arg(timer.elapsed() / 1000.0, 4,
'f', 2);
556 const QRect *roi,
const double maxHFR)
const
558 auto centers = starCenters;
560 const int numDetections = centers.size();
561 for (
int i = 0; i < numDetections; ++i) scores->
push_back(0);
562 if (numDetections == 0)
return;
567 constexpr double snrWeight = 20;
568 auto bg = skybackground();
572 std::sort(centers.begin(), centers.end(), [&bg](
const Edge * a,
const Edge * b)
574 double snrA = bg.SNR(a->sum, a->numPixels);
575 double snrB = bg.SNR(b->sum, b->numPixels);
580 int numRejectedByHFR = 0;
581 for (
int i = 0; i < numDetections; ++i)
583 if (centers.at(i)->HFR > maxHFR)
586 const int starsRemaining = numDetections - numRejectedByHFR;
587 const bool useHFRConstraint =
588 starsRemaining > 5 ||
589 (starsRemaining >= 3 && numDetections <= 6) ||
590 (starsRemaining >= 2 && numDetections <= 4) ||
591 (starsRemaining >= 1 && numDetections <= 2);
593 for (
int i = 0; i < numDetections; ++i)
595 for (
int j = 0; j < numDetections; ++j)
597 if (starCenters.
at(j) == centers.at(i))
600 if (useHFRConstraint && centers.at(i)->HFR > maxHFR)
603 (*scores)[j] += snrWeight * i;
613 for (
int j = 0; j < starCenters.
size(); ++j)
615 const auto &star = starCenters.
at(j);
616 if (star->x < roi->
x() || star->x >= roi->
x() + roi->
width() ||
617 star->y < roi->
y() || star->y >= roi->
y() + roi->
height())
624void GuideStars::computeStarDrift(
const Edge &star,
const Edge &reference,
double *driftRA,
double *driftDEC)
const
626 if (!calibrationInitialized)
return;
628 GuiderUtils::Vector position(star.x, star.y, 0);
629 GuiderUtils::Vector reference_position(reference.x, reference.y, 0);
630 GuiderUtils::Vector arc_position, arc_reference_position;
631 arc_position = calibration.convertToArcseconds(position);
632 arc_reference_position = calibration.convertToArcseconds(reference_position);
635 GuiderUtils::Vector sky_coords = arc_position - arc_reference_position;
636 sky_coords = calibration.rotateToRaDec(sky_coords);
639 *driftRA = sky_coords.x;
640 *driftDEC = sky_coords.y;
650bool GuideStars::getDrift(
double oneStarDrift,
double reticle_x,
double reticle_y,
651 double *RADrift,
double *DECDrift)
653 if (starCorrespondence.size() == 0)
655 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
656 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar getDrift, reticle:" << reticle_x << reticle_y
657 <<
"guidestar" << gStar.x << gStar.y
658 <<
"so offsets:" << (reticle_x - gStar.x) << (reticle_y - gStar.y);
661 constexpr double maxDriftForMultistar = 4000000.0;
663 if (!calibrationInitialized)
667 if (oneStarDrift > maxDriftForMultistar)
669 qCDebug(KSTARS_EKOS_GUIDE) <<
"MultiStar: 1-star drift too high: " << oneStarDrift;
672 if (starCorrespondence.size() < 2)
674 qCDebug(KSTARS_EKOS_GUIDE) <<
"Multistar: no guide stars";
679 double offset_x = reticle_x - gStar.x;
680 double offset_y = reticle_y - gStar.y;
682 int numStarsProcessed = 0;
683 bool guideStarProcessed =
false;
684 double driftRASum = 0, driftDECSum = 0;
685 double guideStarRADrift = 0, guideStarDECDrift = 0;
688 DLOG(KSTARS_EKOS_GUIDE)
689 <<
QString(
"%1 %2 dRA dDEC").
arg(logHeader(
"")).
arg(logHeader(
" Ref:"));
690 for (
int i = 0; i < detectedStars.
size(); ++i)
692 const auto &star = detectedStars[i];
693 auto bg = skybackground();
694 double snr = bg.SNR(detectedStars[i].sum, detectedStars[i].numPixels);
696 if (getStarMap(i) >= 0 && snr >= MIN_DRIFT_SNR)
698 auto ref = starCorrespondence.reference(getStarMap(i));
701 double driftRA, driftDEC;
702 computeStarDrift(star, ref, &driftRA, &driftDEC);
703 if (getStarMap(i) == starCorrespondence.guideStar())
705 guideStarRADrift = driftRA;
706 guideStarDECDrift = driftDEC;
707 guideStarProcessed =
true;
714 DLOG(KSTARS_EKOS_GUIDE)
715 <<
QString(
"%1 %2 %3 %4").
arg(logStar(
"MultiStar", i, bg, star))
716 .
arg(logStar(
" Ref:", getStarMap(i), bg, ref))
717 .
arg(driftRA, 5,
'f', 2).
arg(driftDEC, 5,
'f', 2);
721 if (numStarsProcessed == 0 || (!allowMissingGuideStar && !guideStarProcessed))
726 for (
int i = 0; i < raDrifts.
size(); ++i)
728 if ((allowMissingGuideStar && !guideStarProcessed) ||
729 ((fabs(raDrifts[i] - guideStarRADrift) < 2.0) &&
730 (fabs(decDrifts[i] - guideStarDECDrift) < 2.0)))
732 driftRASum += raDrifts[i];
733 driftDECSum += decDrifts[i];
738 if (raDriftsKeep.
empty() || raDriftsKeep.
empty())
741 numStarsProcessed = raDriftsKeep.
size();
744 bool useMedian =
true;
747 std::sort(raDriftsKeep.
begin(), raDriftsKeep.
end(), [] (
double d1,
double d2)
751 std::sort(decDriftsKeep.
begin(), decDriftsKeep.
end(), [] (
double d1,
double d2)
755 if (numStarsProcessed % 2 == 0)
757 const int middle = numStarsProcessed / 2;
758 *RADrift = ( raDriftsKeep[middle - 1] + raDriftsKeep[middle]) / 2.0;
759 *DECDrift = (decDriftsKeep[middle - 1] + decDriftsKeep[middle]) / 2.0;
763 const int middle = numStarsProcessed / 2;
764 *RADrift = raDriftsKeep[middle];
765 *DECDrift = decDriftsKeep[middle];
767 DLOG(KSTARS_EKOS_GUIDE) <<
"MultiStar: Drift median " << *RADrift << *DECDrift << numStarsProcessed <<
" of " <<
768 detectedStars.
size() <<
"#guide" << starCorrespondence.size();
772 *RADrift = driftRASum / raDriftsKeep.
size();
773 *DECDrift = driftDECSum / decDriftsKeep.
size();
774 DLOG(KSTARS_EKOS_GUIDE) <<
"MultiStar: Drift " << *RADrift << *DECDrift << numStarsProcessed <<
" of " <<
775 detectedStars.
size() <<
"#guide" << starCorrespondence.size();
780void GuideStars::setCalibration(
const Calibration &cal)
783 calibrationInitialized =
true;
void append(QList< T > &&value)
const_reference at(qsizetype i) const const
qsizetype count() const const
void push_back(parameter_type value)
qsizetype size() const const
bool isValid() const const
QString & append(QChar ch)
QString arg(Args &&... args) const const
QString number(double n, char format, int precision)
qsizetype size() const const
QTextStream & center(QTextStream &stream)