File indexing completed on 2024-05-12 15:23:43
0001 /* 0002 SPDX-FileCopyrightText: 2020 Hy Murveit <hy@murveit.com> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "guidestars.h" 0008 0009 #include "ekos_guide_debug.h" 0010 #include "../guideview.h" 0011 #include "fitsviewer/fitsdata.h" 0012 #include "fitsviewer/fitssepdetector.h" 0013 #include "Options.h" 0014 0015 #include <math.h> 0016 #include <stellarsolver.h> 0017 #include "ekos/auxiliary/stellarsolverprofileeditor.h" 0018 #include <QTime> 0019 0020 #define DLOG if (false) qCDebug 0021 0022 // Then when looking for the guide star, gets this many candidates. 0023 #define STARS_TO_SEARCH 250 0024 0025 // Don't use stars with SNR lower than this when computing multi-star drift. 0026 #define MIN_DRIFT_SNR 8 0027 0028 // If there are fewer than this many stars, don't use this multi-star algorithm. 0029 // It will instead back-off to a reticle-based algorithm. 0030 #define MIN_STAR_CORRESPONDENCE_SIZE 5 0031 0032 // We limit the HFR for guide stars. When searching for the guide star, we relax this by the 0033 // margin below (e.g. if a guide star was selected that was near the max guide-star hfr, the later 0034 // the hfr increased a little, we still want to be able to find it. 0035 constexpr double HFR_MARGIN = 2.0; 0036 /* 0037 Start with a set of reference (x,y) positions from stars, where one is designated a guide star. 0038 Given these and a set of new input stars, determine a mapping of new stars to the references. 0039 As the star positions should mostly change via translation, the correspondences are determined by 0040 similar geometry. It is possible, though that there is noise in positions, and that 0041 some reference stars may not be in the new star group, or that the new stars may contain some 0042 stars not appearing in references. However, the guide star must appear in both for this to 0043 be successful. 0044 */ 0045 0046 namespace 0047 { 0048 0049 QString logHeader(const QString &label) 0050 { 0051 return QString("%1 %2 %3 %4 %5 %6 %7") 0052 .arg(label, -9).arg("#", 3).arg("x ", 6).arg("y ", 6).arg("flux", 6) 0053 .arg("HFR", 6).arg("SNR ", 5); 0054 } 0055 0056 QString logStar(const QString &label, int index, const SkyBackground &bg, const Edge &star) 0057 { 0058 double snr = bg.SNR(star.sum, star.numPixels); 0059 return QString("%1 %2 %3 %4 %5 %6 %7") 0060 .arg(label, -9) 0061 .arg(index, 3) 0062 .arg(star.x, 6, 'f', 1) 0063 .arg(star.y, 6, 'f', 1) 0064 .arg(star.sum, 6, 'f', 0) 0065 .arg(star.HFR, 6, 'f', 2) 0066 .arg(snr, 5, 'f', 1); 0067 } 0068 0069 void logStars(const QString &label, const QString &label2, const SkyBackground &bg, 0070 int size, 0071 std::function<const Edge & (int index)> stars, 0072 const QString &extraLabel, 0073 std::function<QString (int index)> extras) 0074 { 0075 QString header = logHeader(label); 0076 if (extraLabel.size() > 0) 0077 header.append(QString(" %1").arg(extraLabel)); 0078 qCDebug(KSTARS_EKOS_GUIDE) << header; 0079 for (int i = 0; i < size; ++i) 0080 { 0081 const auto &star = stars(i); 0082 QString line = logStar(label2, i, bg, star); 0083 if (extraLabel.size() > 0) 0084 line.append(QString(" %1").arg(extras(i))); 0085 qCDebug(KSTARS_EKOS_GUIDE) << line; 0086 } 0087 } 0088 } //namespace 0089 0090 GuideStars::GuideStars() 0091 { 0092 } 0093 0094 // It's possible that we don't map all the stars, if there are too many. 0095 int GuideStars::getStarMap(int index) 0096 { 0097 if (index >= starMap.size() || (index < 0)) 0098 return -1; 0099 return starMap[index]; 0100 } 0101 0102 void GuideStars::setupStarCorrespondence(const QList<Edge> &neighbors, int guideIndex) 0103 { 0104 qCDebug(KSTARS_EKOS_GUIDE) << "setupStarCorrespondence: neighbors " << neighbors.size() << "guide index" << guideIndex; 0105 if (neighbors.size() >= MIN_STAR_CORRESPONDENCE_SIZE) 0106 { 0107 starMap.clear(); 0108 for (int i = 0; i < neighbors.size(); ++i) 0109 { 0110 qCDebug(KSTARS_EKOS_GUIDE) << " adding ref: " << neighbors[i].x << neighbors[i].y; 0111 // We need to initialize starMap, because findGuideStar(), which normally 0112 // initializes starMap(), might call selectGuideStar() if it finds that 0113 // the starCorrespondence is empty. 0114 starMap.push_back(i); 0115 } 0116 starCorrespondence.initialize(neighbors, guideIndex); 0117 } 0118 else 0119 starCorrespondence.reset(); 0120 } 0121 0122 // Calls SEP to generate a set of star detections and score them, 0123 // then calls selectGuideStars(detections, scores) to select the guide star. 0124 QVector3D GuideStars::selectGuideStar(const QSharedPointer<FITSData> &imageData) 0125 { 0126 if (imageData == nullptr) 0127 return QVector3D(-1, -1, -1); 0128 0129 QList<double> sepScores; 0130 QList<double> minDistances; 0131 const double maxHFR = Options::guideMaxHFR(); 0132 findTopStars(imageData, Options::maxMultistarReferenceStars(), &detectedStars, maxHFR, 0133 nullptr, &sepScores, &minDistances); 0134 0135 int maxX = imageData->width(); 0136 int maxY = imageData->height(); 0137 return selectGuideStar(detectedStars, sepScores, maxX, maxY, minDistances); 0138 } 0139 0140 // This selects the guide star by using the input scores, and 0141 // downgrading stars too close to the edges of the image (so that calibration won't fail). 0142 // It also removes from consideration stars with huge HFR values (most likely not stars). 0143 QVector3D GuideStars::selectGuideStar(const QList<Edge> &stars, 0144 const QList<double> &sepScores, 0145 int maxX, int maxY, 0146 const QList<double> &minDistances) 0147 { 0148 constexpr int maxStarDiameter = 32; 0149 0150 if ((uint) stars.size() < Options::minDetectionsSEPMultistar()) 0151 { 0152 qCDebug(KSTARS_EKOS_GUIDE) << "Too few stars detected: " << stars.size(); 0153 return QVector3D(-1, -1, -1); 0154 } 0155 0156 int maxIndex = Options::maxMultistarReferenceStars() < stars.count() ? Options::maxMultistarReferenceStars() : 0157 stars.count(); 0158 QVector<int> scores(Options::maxMultistarReferenceStars()); 0159 0160 // Find the bottom 25% HFR value. 0161 QVector<float> hfrs; 0162 for (int i = 0; i < maxIndex; i++) 0163 hfrs.push_back(stars[i].HFR); 0164 std::sort(hfrs.begin(), hfrs.end()); 0165 float tooLowHFR = 0.0; 0166 if (maxIndex / 4 > 0) 0167 tooLowHFR = hfrs[maxIndex / 4]; 0168 0169 QList<Edge> guideStarNeighbors; 0170 for (int i = 0; i < maxIndex; i++) 0171 { 0172 int score = 100 + sepScores[i]; 0173 const Edge ¢er = stars.at(i); 0174 guideStarNeighbors.append(center); 0175 0176 // Severely reject stars close to edges 0177 // Worry about calibration? Make calibration distance parameter? 0178 constexpr int borderGuard = 35; 0179 if (center.x < borderGuard || center.y < borderGuard || 0180 center.x > (maxX - borderGuard) || center.y > (maxY - borderGuard)) 0181 score -= 1000; 0182 // Reject stars bigger than square 0183 if (center.HFR > float(maxStarDiameter) / 2) 0184 score -= 1000; 0185 0186 // Try not to use a star with an HFR in bottom 25%. 0187 if (center.HFR < tooLowHFR) 0188 score -= 1000; 0189 0190 // Add advantage to stars with SNRs between 40-100. 0191 auto bg = skybackground(); 0192 double snr = bg.SNR(center.sum, center.numPixels); 0193 if (snr >= 40 && snr <= 100) 0194 score += 75; 0195 if (snr >= 100) 0196 score -= 50; 0197 0198 // discourage stars that have close neighbors. 0199 // This isn't perfect, as we're not detecting all stars, just top 100 or so. 0200 const double neighborDistance = minDistances[i]; 0201 constexpr double closeNeighbor = 25; 0202 constexpr double veryCloseNeighbor = 15; 0203 if (neighborDistance < veryCloseNeighbor) 0204 score -= 100; 0205 else if (neighborDistance < closeNeighbor) 0206 score -= 50; 0207 0208 scores[i] = score; 0209 } 0210 0211 logStars("Select", "Star", skyBackground, 0212 maxIndex, 0213 [&](int i) -> const Edge& 0214 { 0215 return stars[i]; 0216 }, 0217 "score", [&](int i) -> QString 0218 { 0219 return QString("%1").arg(scores[i], 5); 0220 }); 0221 0222 int maxScore = -1; 0223 int maxScoreIndex = -1; 0224 for (int i = 0; i < maxIndex; i++) 0225 { 0226 if (scores[i] > maxScore) 0227 { 0228 maxScore = scores[i]; 0229 maxScoreIndex = i; 0230 } 0231 } 0232 0233 if (maxScoreIndex < 0) 0234 { 0235 qCDebug(KSTARS_EKOS_GUIDE) << "No suitable star detected."; 0236 return QVector3D(-1, -1, -1); 0237 } 0238 setupStarCorrespondence(guideStarNeighbors, maxScoreIndex); 0239 QVector3D newStarCenter(stars[maxScoreIndex].x, stars[maxScoreIndex].y, 0); 0240 qCDebug(KSTARS_EKOS_GUIDE) << "new star center: " << maxScoreIndex << " x: " 0241 << stars[maxScoreIndex].x << " y: " << stars[maxScoreIndex].y; 0242 return newStarCenter; 0243 } 0244 0245 // Find the current target positions for the guide-star neighbors, and add them 0246 // to the guideView. 0247 void GuideStars::plotStars(QSharedPointer<GuideView> &guideView, const QRect &trackingBox) 0248 { 0249 if (guideView == nullptr) return; 0250 guideView->clearNeighbors(); 0251 if (starCorrespondence.size() == 0) return; 0252 0253 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar()); 0254 0255 // Find the offset between the current reticle position and the original 0256 // guide star position (e.g. it might have moved due to dithering). 0257 double reticle_x = trackingBox.x() + trackingBox.width() / 2.0; 0258 double reticle_y = trackingBox.y() + trackingBox.height() / 2.0; 0259 0260 // See which neighbor stars were associated with detected stars. 0261 QVector<bool> found(starCorrespondence.size(), false); 0262 QVector<int> detectedIndex(starCorrespondence.size(), -1); 0263 0264 for (int i = 0; i < detectedStars.size(); ++i) 0265 { 0266 int refIndex = getStarMap(i); 0267 if (refIndex >= 0) 0268 { 0269 found[refIndex] = true; 0270 detectedIndex[refIndex] = i; 0271 } 0272 } 0273 // Add to the guideView the neighbor stars' target positions and whether they were found. 0274 for (int i = 0; i < starCorrespondence.size(); ++i) 0275 { 0276 bool isGuideStar = i == starCorrespondence.guideStar(); 0277 const QVector2D offset = starCorrespondence.offset(i); 0278 const double detected_x = found[i] ? detectedStars[detectedIndex[i]].x : 0; 0279 const double detected_y = found[i] ? detectedStars[detectedIndex[i]].y : 0; 0280 guideView->addGuideStarNeighbor(offset.x() + reticle_x, offset.y() + reticle_y, found[i], 0281 detected_x, detected_y, isGuideStar); 0282 } 0283 guideView->updateNeighbors(); 0284 } 0285 0286 // Find the guide star using the starCorrespondence algorithm (looking for 0287 // the other reference stars in the same relative position as when the guide star was selected). 0288 // If this method fails, it backs off to looking in the tracking box for the highest scoring star. 0289 GuiderUtils::Vector GuideStars::findGuideStar(const QSharedPointer<FITSData> &imageData, const QRect &trackingBox, 0290 QSharedPointer<GuideView> &guideView, bool firstFrame) 0291 { 0292 // We fall-back to single star guiding if multistar isn't working, 0293 // but limit this for a number of iterations. 0294 // If a user doesn't have multiple stars available, the user shouldn't be using multistar. 0295 constexpr int MAX_CONSECUTIVE_UNRELIABLE = 10; 0296 if (firstFrame) 0297 unreliableDectionCounter = 0; 0298 0299 // Don't accept reference stars whose position is more than this many pixels from expected. 0300 constexpr double maxStarAssociationDistance = 10; 0301 0302 if (imageData == nullptr) 0303 return GuiderUtils::Vector(-1, -1, -1); 0304 0305 // If the guide star has not yet been set up, then establish it here. 0306 // Not thrilled doing this, but this is the way the internal guider is setup 0307 // when guiding is restarted by the scheduler (the normal establish guide star 0308 // methods are not called). 0309 if (firstFrame && starCorrespondence.size() == 0) 0310 { 0311 QVector3D v = selectGuideStar(imageData); 0312 qCDebug(KSTARS_EKOS_GUIDE) << QString("findGuideStar: Called without starCorrespondence. Refound guide star at %1 %2") 0313 .arg(QString::number(v.x(), 'f', 1)).arg(QString::number(v.y(), 'f', 1)); 0314 return GuiderUtils::Vector(v.x(), v.y(), v.z()); 0315 } 0316 0317 // Allow a little margin above the max hfr for guide stars when searching for the guide star. 0318 const double maxHFR = Options::guideMaxHFR() + HFR_MARGIN; 0319 if (starCorrespondence.size() > 0) 0320 { 0321 0322 findTopStars(imageData, STARS_TO_SEARCH, &detectedStars, maxHFR); 0323 if (detectedStars.empty()) 0324 return GuiderUtils::Vector(-1, -1, -1); 0325 0326 // Allow it to guide even if the main guide star isn't detected (as long as enough reference stars are). 0327 starCorrespondence.setAllowMissingGuideStar(allowMissingGuideStar); 0328 0329 // Star correspondence can run quicker if it knows the image size. 0330 starCorrespondence.setImageSize(imageData->width(), imageData->height()); 0331 0332 // When using large star-correspondence sets and filtering with a StellarSolver profile, 0333 // the stars at the edge of detection can be lost. Best not to filter, but... 0334 double minFraction = 0.5; 0335 if (starCorrespondence.size() > 25) minFraction = 0.33; 0336 else if (starCorrespondence.size() > 15) minFraction = 0.4; 0337 0338 Edge foundStar = starCorrespondence.find(detectedStars, maxStarAssociationDistance, &starMap, true, minFraction); 0339 0340 // Is there a correspondence to the guide star 0341 // Should we also weight distance to the tracking box? 0342 for (int i = 0; i < detectedStars.size(); ++i) 0343 { 0344 if (getStarMap(i) == starCorrespondence.guideStar()) 0345 { 0346 auto &star = detectedStars[i]; 0347 double SNR = skyBackground.SNR(star.sum, star.numPixels); 0348 guideStarSNR = SNR; 0349 guideStarMass = star.sum; 0350 unreliableDectionCounter = 0; 0351 qCDebug(KSTARS_EKOS_GUIDE) << QString("StarCorrespondence found star %1 at %2 %3 SNR %4") 0352 .arg(i).arg(star.x, 0, 'f', 1).arg(star.y, 0, 'f', 1).arg(SNR, 0, 'f', 1); 0353 0354 if (guideView != nullptr) 0355 plotStars(guideView, trackingBox); 0356 return GuiderUtils::Vector(star.x, star.y, 0); 0357 } 0358 } 0359 // None of the stars matched the guide star, but it's possible star correspondence 0360 // invented a guide star position. 0361 if (foundStar.x >= 0 && foundStar.y >= 0) 0362 { 0363 guideStarSNR = skyBackground.SNR(foundStar.sum, foundStar.numPixels); 0364 guideStarMass = foundStar.sum; 0365 unreliableDectionCounter = 0; // debating this 0366 qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence invented at" << foundStar.x << foundStar.y << "SNR" << guideStarSNR; 0367 if (guideView != nullptr) 0368 plotStars(guideView, trackingBox); 0369 return GuiderUtils::Vector(foundStar.x, foundStar.y, 0); 0370 } 0371 } 0372 0373 qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence not used. It failed to find the guide star."; 0374 0375 if (++unreliableDectionCounter > MAX_CONSECUTIVE_UNRELIABLE) 0376 return GuiderUtils::Vector(-1, -1, -1); 0377 0378 logStars("findGuide", "Star", skyBackground, 0379 detectedStars.size(), 0380 [&](int i) -> const Edge& 0381 { 0382 return detectedStars[i]; 0383 }, 0384 "assoc", [&](int i) -> QString 0385 { 0386 return QString("%1").arg(getStarMap(i)); 0387 }); 0388 0389 if (trackingBox.isValid() == false) 0390 return GuiderUtils::Vector(-1, -1, -1); 0391 0392 // If we didn't find a star that way, then fall back 0393 findTopStars(imageData, 100, &detectedStars, maxHFR, &trackingBox); 0394 if (detectedStars.size() > 0) 0395 { 0396 // Find the star closest to the guide star position, if we have a position. 0397 // Otherwise use the center of the tracking box (which must be valid, see above if). 0398 // 1. Get the guide star position 0399 int best = 0; 0400 double refX = trackingBox.x() + trackingBox.width() / 2; 0401 double refY = trackingBox.y() + trackingBox.height() / 2; 0402 if (starCorrespondence.size() > 0 && starCorrespondence.guideStar() >= 0) 0403 { 0404 const auto gStar = starCorrespondence.reference(starCorrespondence.guideStar()); 0405 refX = gStar.x; 0406 refY = gStar.y; 0407 } 0408 // 2. Find the closest star to that position. 0409 double minDistSq = 1e8; 0410 for (int i = 0; i < detectedStars.size(); ++i) 0411 { 0412 const auto &dStar = detectedStars[i]; 0413 const double distSq = (dStar.x - refX) * (dStar.x - refX) + (dStar.y - refY) * (dStar.y - refY); 0414 if (distSq < minDistSq) 0415 { 0416 best = i; 0417 minDistSq = distSq; 0418 } 0419 } 0420 // 3. Return that star. 0421 auto &star = detectedStars[best]; 0422 double SNR = skyBackground.SNR(star.sum, star.numPixels); 0423 guideStarSNR = SNR; 0424 guideStarMass = star.sum; 0425 qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence. Standard method found at " << star.x << star.y << "SNR" << SNR; 0426 return GuiderUtils::Vector(star.x, star.y, 0); 0427 } 0428 return GuiderUtils::Vector(-1, -1, -1); 0429 } 0430 0431 SSolver::Parameters GuideStars::getStarExtractionParameters(int num) 0432 { 0433 SSolver::Parameters params; 0434 params.listName = "Guider"; 0435 params.apertureShape = SSolver::SHAPE_CIRCLE; 0436 params.minarea = 10; // changed from 5 --> 10 0437 params.deblend_thresh = 32; 0438 params.deblend_contrast = 0.005; 0439 params.initialKeep = num * 2; 0440 params.keepNum = num; 0441 params.removeBrightest = 0; 0442 params.removeDimmest = 0; 0443 params.saturationLimit = 0; 0444 return params; 0445 } 0446 0447 // This is the interface to star detection. 0448 int GuideStars::findAllSEPStars(const QSharedPointer<FITSData> &imageData, QList<Edge *> *sepStars, int num) 0449 { 0450 if (imageData == nullptr) 0451 return 0; 0452 0453 QVariantMap settings; 0454 settings["optionsProfileIndex"] = Options::guideOptionsProfile(); 0455 settings["optionsProfileGroup"] = static_cast<int>(Ekos::GuideProfiles); 0456 imageData->setSourceExtractorSettings(settings); 0457 imageData->findStars(ALGORITHM_SEP).waitForFinished(); 0458 skyBackground = imageData->getSkyBackground(); 0459 0460 QList<Edge *> edges = imageData->getStarCenters(); 0461 // Let's sort edges, starting with widest 0462 std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->HFR > edge2->HFR;}); 0463 0464 m_NumStarsDetected = edges.count(); 0465 // Take only the first num stars 0466 { 0467 int starCount = qMin(num, edges.count()); 0468 for (int i = 0; i < starCount; i++) 0469 sepStars->append(edges[i]); 0470 } 0471 qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: SEP detected" << edges.count() << "stars, keeping" << sepStars->size(); 0472 0473 edges.clear(); 0474 return sepStars->count(); 0475 } 0476 0477 double GuideStars::findMinDistance(int index, const QList<Edge*> &stars) 0478 { 0479 double closestDistanceSqr = 1e10; 0480 const Edge &star = *stars[index]; 0481 for (int i = 0; i < stars.size(); ++i) 0482 { 0483 if (i == index) continue; 0484 const Edge &thisStar = *stars[i]; 0485 const double xDist = star.x - thisStar.x; 0486 const double yDist = star.y - thisStar.y; 0487 const double distanceSqr = xDist * xDist + yDist * yDist; 0488 if (distanceSqr < closestDistanceSqr) 0489 closestDistanceSqr = distanceSqr; 0490 } 0491 return sqrt(closestDistanceSqr); 0492 } 0493 0494 // Returns a list of 'num' stars, sorted according to evaluateSEPStars(). 0495 // If the region-of-interest rectange is not null, it only returns scores in that area. 0496 void GuideStars::findTopStars(const QSharedPointer<FITSData> &imageData, int num, QList<Edge> *stars, 0497 const double maxHFR, const QRect *roi, 0498 QList<double> *outputScores, QList<double> *minDistances) 0499 { 0500 if (roi == nullptr) 0501 DLOG(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num; 0502 else 0503 DLOG(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num << 0504 QString("Roi(%1,%2 %3x%4)").arg(roi->x()).arg(roi->y()).arg(roi->width()).arg(roi->height()); 0505 0506 stars->clear(); 0507 if (imageData == nullptr) 0508 return; 0509 0510 QElapsedTimer timer; 0511 timer.restart(); 0512 QList<Edge*> sepStars; 0513 int count = findAllSEPStars(imageData, &sepStars, num * 2); 0514 if (count == 0) 0515 return; 0516 0517 if (sepStars.empty()) 0518 return; 0519 0520 QVector<double> scores; 0521 evaluateSEPStars(sepStars, &scores, roi, maxHFR); 0522 // Sort the sepStars by score, higher score to lower score. 0523 QVector<std::pair<int, double>> sc; 0524 for (int i = 0; i < scores.size(); ++i) 0525 sc.push_back(std::pair<int, double>(i, scores[i])); 0526 std::sort(sc.begin(), sc.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b) 0527 { 0528 return a.second > b.second; 0529 }); 0530 // Copy the top num results. 0531 for (int i = 0; i < std::min(num, scores.size()); ++i) 0532 { 0533 const int starIndex = sc[i].first; 0534 const double starScore = sc[i].second; 0535 if (starScore >= 0) 0536 { 0537 stars->append(*(sepStars[starIndex])); 0538 if (outputScores != nullptr) 0539 outputScores->append(starScore); 0540 if (minDistances != nullptr) 0541 minDistances->append(findMinDistance(starIndex, sepStars)); 0542 } 0543 } 0544 DLOG(KSTARS_EKOS_GUIDE) 0545 << QString("Multistar: findTopStars returning: %1 stars, %2s") 0546 .arg(stars->size()).arg(timer.elapsed() / 1000.0, 4, 'f', 2); 0547 } 0548 0549 // Scores star detection relative to each other. Uses the star's SNR as the main measure. 0550 void GuideStars::evaluateSEPStars(const QList<Edge *> &starCenters, QVector<double> *scores, 0551 const QRect *roi, const double maxHFR) const 0552 { 0553 auto centers = starCenters; 0554 scores->clear(); 0555 const int numDetections = centers.size(); 0556 for (int i = 0; i < numDetections; ++i) scores->push_back(0); 0557 if (numDetections == 0) return; 0558 0559 0560 // Rough constants used by this weighting. 0561 // If the center pixel is above this, assume it's clipped and don't emphasize. 0562 constexpr double snrWeight = 20; // Measure weight if/when multiple measures are used. 0563 auto bg = skybackground(); 0564 0565 // Sort by SNR in increasing order so the weighting goes up. 0566 // Assign score based on the sorted position. 0567 std::sort(centers.begin(), centers.end(), [&bg](const Edge * a, const Edge * b) 0568 { 0569 double snrA = bg.SNR(a->sum, a->numPixels); 0570 double snrB = bg.SNR(b->sum, b->numPixels); 0571 return snrA < snrB; 0572 }); 0573 0574 // See if the HFR restriction is too severe. 0575 int numRejectedByHFR = 0; 0576 for (int i = 0; i < numDetections; ++i) 0577 { 0578 if (centers.at(i)->HFR > maxHFR) 0579 numRejectedByHFR++; 0580 } 0581 const int starsRemaining = numDetections - numRejectedByHFR; 0582 const bool useHFRConstraint = 0583 starsRemaining > 5 || 0584 (starsRemaining >= 3 && numDetections <= 6) || 0585 (starsRemaining >= 2 && numDetections <= 4) || 0586 (starsRemaining >= 1 && numDetections <= 2); 0587 0588 for (int i = 0; i < numDetections; ++i) 0589 { 0590 for (int j = 0; j < numDetections; ++j) 0591 { 0592 if (starCenters.at(j) == centers.at(i)) 0593 { 0594 // Don't emphasize stars that are too wide. 0595 if (useHFRConstraint && centers.at(i)->HFR > maxHFR) 0596 (*scores)[j] = -1; 0597 else 0598 (*scores)[j] += snrWeight * i; 0599 break; 0600 } 0601 } 0602 } 0603 0604 // If we are insisting on a star in the tracking box. 0605 if (roi != nullptr) 0606 { 0607 // Reset the score of anything outside the tracking box. 0608 for (int j = 0; j < starCenters.size(); ++j) 0609 { 0610 const auto &star = starCenters.at(j); 0611 if (star->x < roi->x() || star->x >= roi->x() + roi->width() || 0612 star->y < roi->y() || star->y >= roi->y() + roi->height()) 0613 (*scores)[j] = -1; 0614 } 0615 } 0616 } 0617 0618 // Given a star detection and a reference star compute the RA and DEC drifts between the two. 0619 void GuideStars::computeStarDrift(const Edge &star, const Edge &reference, double *driftRA, double *driftDEC) const 0620 { 0621 if (!calibrationInitialized) return; 0622 0623 GuiderUtils::Vector position(star.x, star.y, 0); 0624 GuiderUtils::Vector reference_position(reference.x, reference.y, 0); 0625 GuiderUtils::Vector arc_position, arc_reference_position; 0626 arc_position = calibration.convertToArcseconds(position); 0627 arc_reference_position = calibration.convertToArcseconds(reference_position); 0628 0629 // translate into sky coords. 0630 GuiderUtils::Vector sky_coords = arc_position - arc_reference_position; 0631 sky_coords = calibration.rotateToRaDec(sky_coords); 0632 0633 // Save the drifts in RA and DEC. 0634 *driftRA = sky_coords.x; 0635 *driftDEC = sky_coords.y; 0636 } 0637 0638 // Computes the overall drift given the input image that was analyzed in findGuideStar(). 0639 // Returns true if this can be done--there are already guide stars and the current detections 0640 // can be compared to them. In that case, fills in RADRift and DECDrift with arc-second drifts. 0641 // Reticle_x,y is the target position. It may be different from when we 0642 // originally selected the guide star, e.g. due to dithering, etc. We compute an offset from the 0643 // original guide-star position and the reticle position and offset all the reference star 0644 // positions by that. 0645 bool GuideStars::getDrift(double oneStarDrift, double reticle_x, double reticle_y, 0646 double *RADrift, double *DECDrift) 0647 { 0648 if (starCorrespondence.size() == 0) 0649 return false; 0650 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar()); 0651 DLOG(KSTARS_EKOS_GUIDE) << "Multistar getDrift, reticle:" << reticle_x << reticle_y 0652 << "guidestar" << gStar.x << gStar.y 0653 << "so offsets:" << (reticle_x - gStar.x) << (reticle_y - gStar.y); 0654 // Revoke multistar if we're that far away. 0655 // Currently disabled by large constant. 0656 constexpr double maxDriftForMultistar = 4000000.0; // arc-seconds 0657 0658 if (!calibrationInitialized) 0659 return false; 0660 0661 // Don't run multistar if the 1-star approach is too far off. 0662 if (oneStarDrift > maxDriftForMultistar) 0663 { 0664 qCDebug(KSTARS_EKOS_GUIDE) << "MultiStar: 1-star drift too high: " << oneStarDrift; 0665 return false; 0666 } 0667 if (starCorrespondence.size() < 2) 0668 { 0669 qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: no guide stars"; 0670 return false; 0671 } 0672 0673 // The original guide star position has been offset. 0674 double offset_x = reticle_x - gStar.x; 0675 double offset_y = reticle_y - gStar.y; 0676 0677 int numStarsProcessed = 0; 0678 bool guideStarProcessed = false; 0679 double driftRASum = 0, driftDECSum = 0; 0680 double guideStarRADrift = 0, guideStarDECDrift = 0; 0681 QVector<double> raDrifts, decDrifts; 0682 0683 DLOG(KSTARS_EKOS_GUIDE) 0684 << QString("%1 %2 dRA dDEC").arg(logHeader("")).arg(logHeader(" Ref:")); 0685 for (int i = 0; i < detectedStars.size(); ++i) 0686 { 0687 const auto &star = detectedStars[i]; 0688 auto bg = skybackground(); 0689 double snr = bg.SNR(detectedStars[i].sum, detectedStars[i].numPixels); 0690 // Probably should test SNR on the reference as well. 0691 if (getStarMap(i) >= 0 && snr >= MIN_DRIFT_SNR) 0692 { 0693 auto ref = starCorrespondence.reference(getStarMap(i)); 0694 ref.x += offset_x; 0695 ref.y += offset_y; 0696 double driftRA, driftDEC; 0697 computeStarDrift(star, ref, &driftRA, &driftDEC); 0698 if (getStarMap(i) == starCorrespondence.guideStar()) 0699 { 0700 guideStarRADrift = driftRA; 0701 guideStarDECDrift = driftDEC; 0702 guideStarProcessed = true; 0703 } 0704 0705 raDrifts.push_back(driftRA); 0706 decDrifts.push_back(driftDEC); 0707 numStarsProcessed++; 0708 0709 DLOG(KSTARS_EKOS_GUIDE) 0710 << QString("%1 %2 %3 %4").arg(logStar("MultiStar", i, bg, star)) 0711 .arg(logStar(" Ref:", getStarMap(i), bg, ref)) 0712 .arg(driftRA, 5, 'f', 2).arg(driftDEC, 5, 'f', 2); 0713 } 0714 } 0715 0716 if (numStarsProcessed == 0 || (!allowMissingGuideStar && !guideStarProcessed)) 0717 return false; 0718 0719 // Filter out reference star drifts that are too different from the guide-star drift. 0720 QVector<double> raDriftsKeep, decDriftsKeep; 0721 for (int i = 0; i < raDrifts.size(); ++i) 0722 { 0723 if ((allowMissingGuideStar && !guideStarProcessed) || 0724 ((fabs(raDrifts[i] - guideStarRADrift) < 2.0) && 0725 (fabs(decDrifts[i] - guideStarDECDrift) < 2.0))) 0726 { 0727 driftRASum += raDrifts[i]; 0728 driftDECSum += decDrifts[i]; 0729 raDriftsKeep.push_back(raDrifts[i]); 0730 decDriftsKeep.push_back(decDrifts[i]); 0731 } 0732 } 0733 if (raDriftsKeep.empty() || raDriftsKeep.empty()) 0734 return false; // Shouldn't happen. 0735 0736 numStarsProcessed = raDriftsKeep.size(); 0737 0738 // Generate the drift either from the median or the average of the usable reference drifts. 0739 bool useMedian = true; 0740 if (useMedian) 0741 { 0742 std::sort(raDriftsKeep.begin(), raDriftsKeep.end(), [] (double d1, double d2) 0743 { 0744 return d1 < d2; 0745 }); 0746 std::sort(decDriftsKeep.begin(), decDriftsKeep.end(), [] (double d1, double d2) 0747 { 0748 return d1 < d2; 0749 }); 0750 if (numStarsProcessed % 2 == 0) 0751 { 0752 const int middle = numStarsProcessed / 2; 0753 *RADrift = ( raDriftsKeep[middle - 1] + raDriftsKeep[middle]) / 2.0; 0754 *DECDrift = (decDriftsKeep[middle - 1] + decDriftsKeep[middle]) / 2.0; 0755 } 0756 else 0757 { 0758 const int middle = numStarsProcessed / 2; // because it's 0-based 0759 *RADrift = raDriftsKeep[middle]; 0760 *DECDrift = decDriftsKeep[middle]; 0761 } 0762 DLOG(KSTARS_EKOS_GUIDE) << "MultiStar: Drift median " << *RADrift << *DECDrift << numStarsProcessed << " of " << 0763 detectedStars.size() << "#guide" << starCorrespondence.size(); 0764 } 0765 else 0766 { 0767 *RADrift = driftRASum / raDriftsKeep.size(); 0768 *DECDrift = driftDECSum / decDriftsKeep.size(); 0769 DLOG(KSTARS_EKOS_GUIDE) << "MultiStar: Drift " << *RADrift << *DECDrift << numStarsProcessed << " of " << 0770 detectedStars.size() << "#guide" << starCorrespondence.size(); 0771 } 0772 return true; 0773 } 0774 0775 void GuideStars::setCalibration(const Calibration &cal) 0776 { 0777 calibration = cal; 0778 calibrationInitialized = true; 0779 }