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();
298 constexpr int MAX_CONSECUTIVE_UNRELIABLE = 10;
300 unreliableDectionCounter = 0;
303 constexpr double maxStarAssociationDistance = 10;
305 if (imageData ==
nullptr)
306 return GuiderUtils::Vector(-1, -1, -1);
312 if (firstFrame && starCorrespondence.size() == 0)
314 QVector3D v = selectGuideStar(imageData);
315 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"findGuideStar: Called without starCorrespondence. Refound guide star at %1 %2")
317 return GuiderUtils::Vector(v.
x(), v.
y(), v.
z());
321 const double maxHFR = Options::guideMaxHFR() + HFR_MARGIN;
322 if (starCorrespondence.size() > 0)
325 findTopStars(imageData, STARS_TO_SEARCH, &detectedStars, maxHFR);
326 if (detectedStars.
empty())
327 return GuiderUtils::Vector(-1, -1, -1);
330 starCorrespondence.setAllowMissingGuideStar(allowMissingGuideStar);
333 starCorrespondence.setImageSize(imageData->width(), imageData->height());
337 double minFraction = 0.5;
338 if (starCorrespondence.size() > 25) minFraction = 0.33;
339 else if (starCorrespondence.size() > 15) minFraction = 0.4;
341 Edge foundStar = starCorrespondence.find(detectedStars, maxStarAssociationDistance, &starMap,
true, minFraction);
345 for (
int i = 0; i < detectedStars.
size(); ++i)
347 if (getStarMap(i) == starCorrespondence.guideStar())
349 auto &star = detectedStars[i];
350 double SNR = skyBackground.SNR(star.sum, star.numPixels);
352 guideStarMass = star.sum;
353 unreliableDectionCounter = 0;
354 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"StarCorrespondence found star %1 at %2 %3 SNR %4")
355 .
arg(i).
arg(star.x, 0,
'f', 1).
arg(star.y, 0,
'f', 1).
arg(SNR, 0,
'f', 1);
357 if (guideView !=
nullptr)
358 plotStars(guideView, trackingBox);
359 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"StarCorrespondence. findGuideStar took %1s").
arg(timer.elapsed() / 1000.0, 0,
'f',
361 return GuiderUtils::Vector(star.x, star.y, 0);
366 if (foundStar.x >= 0 && foundStar.y >= 0)
368 guideStarSNR = skyBackground.SNR(foundStar.sum, foundStar.numPixels);
369 guideStarMass = foundStar.sum;
370 unreliableDectionCounter = 0;
371 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence invented at" << foundStar.x << foundStar.y <<
"SNR" << guideStarSNR;
372 if (guideView !=
nullptr)
373 plotStars(guideView, trackingBox);
374 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"StarCorrespondence. findGuideStar/invent took %1s").
arg(timer.elapsed() / 1000.0, 0,
376 return GuiderUtils::Vector(foundStar.x, foundStar.y, 0);
380 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence not used. It failed to find the guide star.";
382 if (++unreliableDectionCounter > MAX_CONSECUTIVE_UNRELIABLE)
383 return GuiderUtils::Vector(-1, -1, -1);
385 logStars(
"findGuide",
"Star", skyBackground,
386 detectedStars.
size(),
387 [&](
int i) ->
const Edge&
389 return detectedStars[i];
393 return QString(
"%1").arg(getStarMap(i));
396 if (trackingBox.
isValid() ==
false)
397 return GuiderUtils::Vector(-1, -1, -1);
400 findTopStars(imageData, 100, &detectedStars, maxHFR, &trackingBox);
401 if (detectedStars.
size() > 0)
407 double refX = trackingBox.
x() + trackingBox.
width() / 2;
408 double refY = trackingBox.
y() + trackingBox.
height() / 2;
409 if (starCorrespondence.size() > 0 && starCorrespondence.guideStar() >= 0)
411 const auto gStar = starCorrespondence.reference(starCorrespondence.guideStar());
416 double minDistSq = 1e8;
417 for (
int i = 0; i < detectedStars.
size(); ++i)
419 const auto &dStar = detectedStars[i];
420 const double distSq = (dStar.x - refX) * (dStar.x - refX) + (dStar.y - refY) * (dStar.y - refY);
421 if (distSq < minDistSq)
428 auto &star = detectedStars[best];
429 double SNR = skyBackground.SNR(star.sum, star.numPixels);
431 guideStarMass = star.sum;
432 qCDebug(KSTARS_EKOS_GUIDE) <<
"StarCorrespondence. Standard method found at " << star.x << star.y <<
"SNR" << SNR;
433 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"StarCorrespondence. findGuideStar backoff took %1s").
arg(timer.elapsed() / 1000.0, 0,
435 return GuiderUtils::Vector(star.x, star.y, 0);
437 return GuiderUtils::Vector(-1, -1, -1);
440SSolver::Parameters GuideStars::getStarExtractionParameters(
int num)
442 SSolver::Parameters params;
443 params.listName =
"Guider";
444 params.apertureShape = SSolver::SHAPE_CIRCLE;
446 params.deblend_thresh = 32;
447 params.deblend_contrast = 0.005;
448 params.initialKeep = num * 2;
449 params.keepNum = num;
450 params.removeBrightest = 0;
451 params.removeDimmest = 0;
452 params.saturationLimit = 0;
459 if (imageData ==
nullptr)
462 QVariantMap settings;
463 settings[
"optionsProfileIndex"] = Options::guideOptionsProfile();
464 settings[
"optionsProfileGroup"] =
static_cast<int>(Ekos::GuideProfiles);
465 imageData->setSourceExtractorSettings(settings);
466 imageData->findStars(ALGORITHM_SEP).waitForFinished();
467 skyBackground = imageData->getSkyBackground();
471 std::sort(edges.
begin(), edges.
end(), [](
const Edge * edge1,
const Edge * edge2) ->
bool { return edge1->HFR > edge2->HFR;});
473 m_NumStarsDetected = edges.
count();
476 int starCount = qMin(num, edges.
count());
477 for (
int i = 0; i < starCount; i++)
478 sepStars->
append(edges[i]);
480 qCDebug(KSTARS_EKOS_GUIDE) <<
"Multistar: SEP detected" << edges.
count() <<
"stars, keeping" << sepStars->
size();
483 return sepStars->
count();
486double GuideStars::findMinDistance(
int index,
const QList<Edge*> &stars)
488 double closestDistanceSqr = 1e10;
489 const Edge &star = *stars[index];
490 for (
int i = 0; i < stars.
size(); ++i)
492 if (i == index)
continue;
493 const Edge &thisStar = *stars[i];
494 const double xDist = star.x - thisStar.x;
495 const double yDist = star.y - thisStar.y;
496 const double distanceSqr = xDist * xDist + yDist * yDist;
497 if (distanceSqr < closestDistanceSqr)
498 closestDistanceSqr = distanceSqr;
500 return sqrt(closestDistanceSqr);
506 const double maxHFR,
const QRect *roi,
513 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar: findTopStars" << num;
517 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar: findTopStars" << num <<
522 if (imageData ==
nullptr)
528 int count = findAllSEPStars(imageData, &sepStars, num * 2);
532 if (sepStars.
empty())
536 evaluateSEPStars(sepStars, &scores, roi, maxHFR);
539 for (
int i = 0; i < scores.
size(); ++i)
540 sc.
push_back(std::pair<int, double>(i, scores[i]));
541 std::sort(sc.
begin(), sc.
end(), [](
const std::pair<int, double> &a,
const std::pair<int, double> &b)
543 return a.second > b.second;
546 for (
int i = 0; i < std::min(num, (
int)scores.
size()); ++i)
548 const int starIndex = sc[i].
first;
549 const double starScore = sc[i].second;
552 stars->
append(*(sepStars[starIndex]));
553 if (outputScores !=
nullptr)
554 outputScores->
append(starScore);
555 if (minDistances !=
nullptr)
556 minDistances->
append(findMinDistance(starIndex, sepStars));
559 DLOG(KSTARS_EKOS_GUIDE)
560 <<
QString(
"Multistar: findTopStars returning: %1 stars, %2s")
561 .
arg(stars->
size()).
arg(timer.elapsed() / 1000.0, 4,
'f', 2);
563 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"FindTopStars. star detection took %1s").
arg(timer2.
elapsed() / 1000.0, 0,
'f', 3);
568 const QRect *roi,
const double maxHFR)
const
570 auto centers = starCenters;
572 const int numDetections = centers.size();
573 for (
int i = 0; i < numDetections; ++i) scores->
push_back(0);
574 if (numDetections == 0)
return;
579 constexpr double snrWeight = 20;
580 auto bg = skybackground();
584 std::sort(centers.begin(), centers.end(), [&bg](
const Edge * a,
const Edge * b)
586 double snrA = bg.SNR(a->sum, a->numPixels);
587 double snrB = bg.SNR(b->sum, b->numPixels);
592 int numRejectedByHFR = 0;
593 for (
int i = 0; i < numDetections; ++i)
595 if (centers.at(i)->HFR > maxHFR)
598 const int starsRemaining = numDetections - numRejectedByHFR;
599 const bool useHFRConstraint =
600 starsRemaining > 5 ||
601 (starsRemaining >= 3 && numDetections <= 6) ||
602 (starsRemaining >= 2 && numDetections <= 4) ||
603 (starsRemaining >= 1 && numDetections <= 2);
605 for (
int i = 0; i < numDetections; ++i)
607 for (
int j = 0; j < numDetections; ++j)
609 if (starCenters.
at(j) == centers.at(i))
612 if (useHFRConstraint && centers.at(i)->HFR > maxHFR)
615 (*scores)[j] += snrWeight * i;
625 for (
int j = 0; j < starCenters.
size(); ++j)
627 const auto &star = starCenters.
at(j);
628 if (star->x < roi->
x() || star->x >= roi->
x() + roi->
width() ||
629 star->y < roi->
y() || star->y >= roi->
y() + roi->
height())
636void GuideStars::computeStarDrift(
const Edge &star,
const Edge &reference,
double *driftRA,
double *driftDEC)
const
638 if (!calibrationInitialized)
return;
640 GuiderUtils::Vector position(star.x, star.y, 0);
641 GuiderUtils::Vector reference_position(reference.x, reference.y, 0);
642 GuiderUtils::Vector arc_position, arc_reference_position;
643 arc_position = calibration.convertToArcseconds(position);
644 arc_reference_position = calibration.convertToArcseconds(reference_position);
647 GuiderUtils::Vector sky_coords = arc_position - arc_reference_position;
648 sky_coords = calibration.rotateToRaDec(sky_coords);
651 *driftRA = sky_coords.x;
652 *driftDEC = sky_coords.y;
662bool GuideStars::getDrift(
double oneStarDrift,
double reticle_x,
double reticle_y,
663 double *RADrift,
double *DECDrift)
665 if (starCorrespondence.size() == 0)
667 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
668 DLOG(KSTARS_EKOS_GUIDE) <<
"Multistar getDrift, reticle:" << reticle_x << reticle_y
669 <<
"guidestar" << gStar.x << gStar.y
670 <<
"so offsets:" << (reticle_x - gStar.x) << (reticle_y - gStar.y);
673 constexpr double maxDriftForMultistar = 4000000.0;
675 if (!calibrationInitialized)
679 if (oneStarDrift > maxDriftForMultistar)
681 qCDebug(KSTARS_EKOS_GUIDE) <<
"MultiStar: 1-star drift too high: " << oneStarDrift;
684 if (starCorrespondence.size() < 2)
686 qCDebug(KSTARS_EKOS_GUIDE) <<
"Multistar: no guide stars";
691 double offset_x = reticle_x - gStar.x;
692 double offset_y = reticle_y - gStar.y;
694 int numStarsProcessed = 0;
695 bool guideStarProcessed =
false;
696 double driftRASum = 0, driftDECSum = 0;
697 double guideStarRADrift = 0, guideStarDECDrift = 0;
700 DLOG(KSTARS_EKOS_GUIDE)
701 <<
QString(
"%1 %2 dRA dDEC").
arg(logHeader(
"")).
arg(logHeader(
" Ref:"));
702 for (
int i = 0; i < detectedStars.
size(); ++i)
704 const auto &star = detectedStars[i];
705 auto bg = skybackground();
706 double snr = bg.SNR(detectedStars[i].sum, detectedStars[i].numPixels);
708 if (getStarMap(i) >= 0 && snr >= MIN_DRIFT_SNR)
710 auto ref = starCorrespondence.reference(getStarMap(i));
713 double driftRA, driftDEC;
714 computeStarDrift(star, ref, &driftRA, &driftDEC);
715 if (getStarMap(i) == starCorrespondence.guideStar())
717 guideStarRADrift = driftRA;
718 guideStarDECDrift = driftDEC;
719 guideStarProcessed =
true;
726 DLOG(KSTARS_EKOS_GUIDE)
727 <<
QString(
"%1 %2 %3 %4").
arg(logStar(
"MultiStar", i, bg, star))
728 .
arg(logStar(
" Ref:", getStarMap(i), bg, ref))
729 .
arg(driftRA, 5,
'f', 2).
arg(driftDEC, 5,
'f', 2);
733 if (numStarsProcessed == 0 || (!allowMissingGuideStar && !guideStarProcessed))
738 for (
int i = 0; i < raDrifts.
size(); ++i)
740 if ((allowMissingGuideStar && !guideStarProcessed) ||
741 ((fabs(raDrifts[i] - guideStarRADrift) < 2.0) &&
742 (fabs(decDrifts[i] - guideStarDECDrift) < 2.0)))
744 driftRASum += raDrifts[i];
745 driftDECSum += decDrifts[i];
750 if (raDriftsKeep.
empty() || raDriftsKeep.
empty())
753 numStarsProcessed = raDriftsKeep.
size();
756 bool useMedian =
true;
759 std::sort(raDriftsKeep.
begin(), raDriftsKeep.
end(), [] (
double d1,
double d2)
763 std::sort(decDriftsKeep.
begin(), decDriftsKeep.
end(), [] (
double d1,
double d2)
767 if (numStarsProcessed % 2 == 0)
769 const int middle = numStarsProcessed / 2;
770 *RADrift = ( raDriftsKeep[middle - 1] + raDriftsKeep[middle]) / 2.0;
771 *DECDrift = (decDriftsKeep[middle - 1] + decDriftsKeep[middle]) / 2.0;
775 const int middle = numStarsProcessed / 2;
776 *RADrift = raDriftsKeep[middle];
777 *DECDrift = decDriftsKeep[middle];
779 DLOG(KSTARS_EKOS_GUIDE) <<
"MultiStar: Drift median " << *RADrift << *DECDrift << numStarsProcessed <<
" of " <<
780 detectedStars.
size() <<
"#guide" << starCorrespondence.size();
784 *RADrift = driftRASum / raDriftsKeep.
size();
785 *DECDrift = driftDECSum / decDriftsKeep.
size();
786 DLOG(KSTARS_EKOS_GUIDE) <<
"MultiStar: Drift " << *RADrift << *DECDrift << numStarsProcessed <<
" of " <<
787 detectedStars.
size() <<
"#guide" << starCorrespondence.size();
792void GuideStars::setCalibration(
const Calibration &cal)
795 calibrationInitialized =
true;
qint64 elapsed() const const
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)