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
129 QList<double> sepScores;
130 QList<double> minDistances;
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,
144 const QList<double> &sepScores,
145 int maxX, int maxY,
146 const QList<double> &minDistances)
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() < static_cast<uint>(stars.count()) ?
157 Options::maxMultistarReferenceStars() :
158 stars.count();
159 QVector<int> scores(Options::maxMultistarReferenceStars());
160
161 // Find the bottom 25% HFR value.
162 QVector<float> hfrs;
163 for (int i = 0; i < maxIndex; i++)
164 hfrs.push_back(stars[i].HFR);
165 std::sort(hfrs.begin(), hfrs.end());
166 float tooLowHFR = 0.0;
167 if (maxIndex / 4 > 0)
168 tooLowHFR = hfrs[maxIndex / 4];
169
170 QList<Edge> guideStarNeighbors;
171 for (int i = 0; i < maxIndex; i++)
172 {
173 int score = 100 + sepScores[i];
174 const Edge &center = stars.at(i);
175 guideStarNeighbors.append(center);
176
177 // Severely reject stars close to edges
178 // Worry about calibration? Make calibration distance parameter?
179 constexpr int borderGuard = 35;
180 if (center.x < borderGuard || center.y < borderGuard ||
181 center.x > (maxX - borderGuard) || center.y > (maxY - borderGuard))
182 score -= 1000;
183 // Reject stars bigger than square
184 if (center.HFR > float(maxStarDiameter) / 2)
185 score -= 1000;
186
187 // Try not to use a star with an HFR in bottom 25%.
188 if (center.HFR < tooLowHFR)
189 score -= 1000;
190
191 // Add advantage to stars with SNRs between 40-100.
192 auto bg = skybackground();
193 double snr = bg.SNR(center.sum, center.numPixels);
194 if (snr >= 40 && snr <= 100)
195 score += 75;
196 if (snr >= 100)
197 score -= 50;
198
199 // discourage stars that have close neighbors.
200 // This isn't perfect, as we're not detecting all stars, just top 100 or so.
201 const double neighborDistance = minDistances[i];
202 constexpr double closeNeighbor = 25;
203 constexpr double veryCloseNeighbor = 15;
204 if (neighborDistance < veryCloseNeighbor)
205 score -= 100;
206 else if (neighborDistance < closeNeighbor)
207 score -= 50;
208
209 scores[i] = score;
210 }
211
212 logStars("Select", "Star", skyBackground,
213 maxIndex,
214 [&](int i) -> const Edge&
215 {
216 return stars[i];
217 },
218 "score", [&](int i) -> QString
219 {
220 return QString("%1").arg(scores[i], 5);
221 });
222
223 int maxScore = -1;
224 int maxScoreIndex = -1;
225 for (int i = 0; i < maxIndex; i++)
226 {
227 if (scores[i] > maxScore)
228 {
229 maxScore = scores[i];
230 maxScoreIndex = i;
231 }
232 }
233
234 if (maxScoreIndex < 0)
235 {
236 qCDebug(KSTARS_EKOS_GUIDE) << "No suitable star detected.";
237 return QVector3D(-1, -1, -1);
238 }
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;
244}
245
246// Find the current target positions for the guide-star neighbors, and add them
247// to the guideView.
248void GuideStars::plotStars(QSharedPointer<GuideView> &guideView, const QRect &trackingBox)
249{
250 if (guideView == nullptr) return;
251 guideView->clearNeighbors();
252 if (starCorrespondence.size() == 0) return;
253
254 const Edge gStar = starCorrespondence.reference(starCorrespondence.guideStar());
255
256 // Find the offset between the current reticle position and the original
257 // guide star position (e.g. it might have moved due to dithering).
258 double reticle_x = trackingBox.x() + trackingBox.width() / 2.0;
259 double reticle_y = trackingBox.y() + trackingBox.height() / 2.0;
260
261 // See which neighbor stars were associated with detected stars.
262 QVector<bool> found(starCorrespondence.size(), false);
263 QVector<int> detectedIndex(starCorrespondence.size(), -1);
264
265 for (int i = 0; i < detectedStars.size(); ++i)
266 {
267 int refIndex = getStarMap(i);
268 if (refIndex >= 0)
269 {
270 found[refIndex] = true;
271 detectedIndex[refIndex] = i;
272 }
273 }
274 // Add to the guideView the neighbor stars' target positions and whether they were found.
275 for (int i = 0; i < starCorrespondence.size(); ++i)
276 {
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);
283 }
284 guideView->updateNeighbors();
285}
286
287// Find the guide star using the starCorrespondence algorithm (looking for
288// the other reference stars in the same relative position as when the guide star was selected).
289// If this method fails, it backs off to looking in the tracking box for the highest scoring star.
290GuiderUtils::Vector GuideStars::findGuideStar(const QSharedPointer<FITSData> &imageData, const QRect &trackingBox,
291 QSharedPointer<GuideView> &guideView, bool firstFrame)
292{
293 QElapsedTimer timer;
294 timer.start();
295 // We fall-back to single star guiding if multistar isn't working,
296 // but limit this for a number of iterations.
297 // If a user doesn't have multiple stars available, the user shouldn't be using multistar.
298 constexpr int MAX_CONSECUTIVE_UNRELIABLE = 10;
299 if (firstFrame)
300 unreliableDectionCounter = 0;
301
302 // Don't accept reference stars whose position is more than this many pixels from expected.
303 constexpr double maxStarAssociationDistance = 10;
304
305 if (imageData == nullptr)
306 return GuiderUtils::Vector(-1, -1, -1);
307
308 // If the guide star has not yet been set up, then establish it here.
309 // Not thrilled doing this, but this is the way the internal guider is setup
310 // when guiding is restarted by the scheduler (the normal establish guide star
311 // methods are not called).
312 if (firstFrame && starCorrespondence.size() == 0)
313 {
314 QVector3D v = selectGuideStar(imageData);
315 qCDebug(KSTARS_EKOS_GUIDE) << QString("findGuideStar: Called without starCorrespondence. Refound guide star at %1 %2")
316 .arg(QString::number(v.x(), 'f', 1)).arg(QString::number(v.y(), 'f', 1));
317 return GuiderUtils::Vector(v.x(), v.y(), v.z());
318 }
319
320 // Allow a little margin above the max hfr for guide stars when searching for the guide star.
321 const double maxHFR = Options::guideMaxHFR() + HFR_MARGIN;
322 if (starCorrespondence.size() > 0)
323 {
324
325 findTopStars(imageData, STARS_TO_SEARCH, &detectedStars, maxHFR);
326 if (detectedStars.empty())
327 return GuiderUtils::Vector(-1, -1, -1);
328
329 // Allow it to guide even if the main guide star isn't detected (as long as enough reference stars are).
330 starCorrespondence.setAllowMissingGuideStar(allowMissingGuideStar);
331
332 // Star correspondence can run quicker if it knows the image size.
333 starCorrespondence.setImageSize(imageData->width(), imageData->height());
334
335 // When using large star-correspondence sets and filtering with a StellarSolver profile,
336 // the stars at the edge of detection can be lost. Best not to filter, but...
337 double minFraction = 0.5;
338 if (starCorrespondence.size() > 25) minFraction = 0.33;
339 else if (starCorrespondence.size() > 15) minFraction = 0.4;
340
341 Edge foundStar = starCorrespondence.find(detectedStars, maxStarAssociationDistance, &starMap, true, minFraction);
342
343 // Is there a correspondence to the guide star
344 // Should we also weight distance to the tracking box?
345 for (int i = 0; i < detectedStars.size(); ++i)
346 {
347 if (getStarMap(i) == starCorrespondence.guideStar())
348 {
349 auto &star = detectedStars[i];
350 double SNR = skyBackground.SNR(star.sum, star.numPixels);
351 guideStarSNR = SNR;
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);
356
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',
360 3);
361 return GuiderUtils::Vector(star.x, star.y, 0);
362 }
363 }
364 // None of the stars matched the guide star, but it's possible star correspondence
365 // invented a guide star position.
366 if (foundStar.x >= 0 && foundStar.y >= 0)
367 {
368 guideStarSNR = skyBackground.SNR(foundStar.sum, foundStar.numPixels);
369 guideStarMass = foundStar.sum;
370 unreliableDectionCounter = 0; // debating this
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,
375 'f', 3);
376 return GuiderUtils::Vector(foundStar.x, foundStar.y, 0);
377 }
378 }
379
380 qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence not used. It failed to find the guide star.";
381
382 if (++unreliableDectionCounter > MAX_CONSECUTIVE_UNRELIABLE)
383 return GuiderUtils::Vector(-1, -1, -1);
384
385 logStars("findGuide", "Star", skyBackground,
386 detectedStars.size(),
387 [&](int i) -> const Edge&
388 {
389 return detectedStars[i];
390 },
391 "assoc", [&](int i) -> QString
392 {
393 return QString("%1").arg(getStarMap(i));
394 });
395
396 if (trackingBox.isValid() == false)
397 return GuiderUtils::Vector(-1, -1, -1);
398
399 // If we didn't find a star that way, then fall back
400 findTopStars(imageData, 100, &detectedStars, maxHFR, &trackingBox);
401 if (detectedStars.size() > 0)
402 {
403 // Find the star closest to the guide star position, if we have a position.
404 // Otherwise use the center of the tracking box (which must be valid, see above if).
405 // 1. Get the guide star position
406 int best = 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)
410 {
411 const auto gStar = starCorrespondence.reference(starCorrespondence.guideStar());
412 refX = gStar.x;
413 refY = gStar.y;
414 }
415 // 2. Find the closest star to that position.
416 double minDistSq = 1e8;
417 for (int i = 0; i < detectedStars.size(); ++i)
418 {
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)
422 {
423 best = i;
424 minDistSq = distSq;
425 }
426 }
427 // 3. Return that star.
428 auto &star = detectedStars[best];
429 double SNR = skyBackground.SNR(star.sum, star.numPixels);
430 guideStarSNR = SNR;
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,
434 'f', 3);
435 return GuiderUtils::Vector(star.x, star.y, 0);
436 }
437 return GuiderUtils::Vector(-1, -1, -1);
438}
439
440SSolver::Parameters GuideStars::getStarExtractionParameters(int num)
441{
442 SSolver::Parameters params;
443 params.listName = "Guider";
444 params.apertureShape = SSolver::SHAPE_CIRCLE;
445 params.minarea = 10; // changed from 5 --> 10
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;
453 return params;
454}
455
456// This is the interface to star detection.
457int GuideStars::findAllSEPStars(const QSharedPointer<FITSData> &imageData, QList<Edge *> *sepStars, int num)
458{
459 if (imageData == nullptr)
460 return 0;
461
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();
468
469 QList<Edge *> edges = imageData->getStarCenters();
470 // Let's sort edges, starting with widest
471 std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->HFR > edge2->HFR;});
472
473 m_NumStarsDetected = edges.count();
474 // Take only the first num stars
475 {
476 int starCount = qMin(num, edges.count());
477 for (int i = 0; i < starCount; i++)
478 sepStars->append(edges[i]);
479 }
480 qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: SEP detected" << edges.count() << "stars, keeping" << sepStars->size();
481
482 edges.clear();
483 return sepStars->count();
484}
485
486double GuideStars::findMinDistance(int index, const QList<Edge*> &stars)
487{
488 double closestDistanceSqr = 1e10;
489 const Edge &star = *stars[index];
490 for (int i = 0; i < stars.size(); ++i)
491 {
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;
499 }
500 return sqrt(closestDistanceSqr);
501}
502
503// Returns a list of 'num' stars, sorted according to evaluateSEPStars().
504// If the region-of-interest rectange is not null, it only returns scores in that area.
505void GuideStars::findTopStars(const QSharedPointer<FITSData> &imageData, int num, QList<Edge> *stars,
506 const double maxHFR, const QRect *roi,
507 QList<double> *outputScores, QList<double> *minDistances)
508{
509 QElapsedTimer timer2;
510 timer2.start();
511 if (roi == nullptr)
512 {
513 DLOG(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num;
514 }
515 else
516 {
517 DLOG(KSTARS_EKOS_GUIDE) << "Multistar: findTopStars" << num <<
518 QString("Roi(%1,%2 %3x%4)").arg(roi->x()).arg(roi->y()).arg(roi->width()).arg(roi->height());
519 }
520
521 stars->clear();
522 if (imageData == nullptr)
523 return;
524
525 QElapsedTimer timer;
526 timer.restart();
527 QList<Edge*> sepStars;
528 int count = findAllSEPStars(imageData, &sepStars, num * 2);
529 if (count == 0)
530 return;
531
532 if (sepStars.empty())
533 return;
534
535 QVector<double> scores;
536 evaluateSEPStars(sepStars, &scores, roi, maxHFR);
537 // Sort the sepStars by score, higher score to lower score.
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)
542 {
543 return a.second > b.second;
544 });
545 // Copy the top num results.
546 for (int i = 0; i < std::min(num, (int)scores.size()); ++i)
547 {
548 const int starIndex = sc[i].first;
549 const double starScore = sc[i].second;
550 if (starScore >= 0)
551 {
552 stars->append(*(sepStars[starIndex]));
553 if (outputScores != nullptr)
554 outputScores->append(starScore);
555 if (minDistances != nullptr)
556 minDistances->append(findMinDistance(starIndex, sepStars));
557 }
558 }
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);
562
563 qCDebug(KSTARS_EKOS_GUIDE) << QString("FindTopStars. star detection took %1s").arg(timer2.elapsed() / 1000.0, 0, 'f', 3);
564}
565
566// Scores star detection relative to each other. Uses the star's SNR as the main measure.
567void GuideStars::evaluateSEPStars(const QList<Edge *> &starCenters, QVector<double> *scores,
568 const QRect *roi, const double maxHFR) const
569{
570 auto centers = starCenters;
571 scores->clear();
572 const int numDetections = centers.size();
573 for (int i = 0; i < numDetections; ++i) scores->push_back(0);
574 if (numDetections == 0) return;
575
576
577 // Rough constants used by this weighting.
578 // If the center pixel is above this, assume it's clipped and don't emphasize.
579 constexpr double snrWeight = 20; // Measure weight if/when multiple measures are used.
580 auto bg = skybackground();
581
582 // Sort by SNR in increasing order so the weighting goes up.
583 // Assign score based on the sorted position.
584 std::sort(centers.begin(), centers.end(), [&bg](const Edge * a, const Edge * b)
585 {
586 double snrA = bg.SNR(a->sum, a->numPixels);
587 double snrB = bg.SNR(b->sum, b->numPixels);
588 return snrA < snrB;
589 });
590
591 // See if the HFR restriction is too severe.
592 int numRejectedByHFR = 0;
593 for (int i = 0; i < numDetections; ++i)
594 {
595 if (centers.at(i)->HFR > maxHFR)
596 numRejectedByHFR++;
597 }
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);
604
605 for (int i = 0; i < numDetections; ++i)
606 {
607 for (int j = 0; j < numDetections; ++j)
608 {
609 if (starCenters.at(j) == centers.at(i))
610 {
611 // Don't emphasize stars that are too wide.
612 if (useHFRConstraint && centers.at(i)->HFR > maxHFR)
613 (*scores)[j] = -1;
614 else
615 (*scores)[j] += snrWeight * i;
616 break;
617 }
618 }
619 }
620
621 // If we are insisting on a star in the tracking box.
622 if (roi != nullptr)
623 {
624 // Reset the score of anything outside the tracking box.
625 for (int j = 0; j < starCenters.size(); ++j)
626 {
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())
630 (*scores)[j] = -1;
631 }
632 }
633}
634
635// Given a star detection and a reference star compute the RA and DEC drifts between the two.
636void GuideStars::computeStarDrift(const Edge &star, const Edge &reference, double *driftRA, double *driftDEC) const
637{
638 if (!calibrationInitialized) return;
639
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);
645
646 // translate into sky coords.
647 GuiderUtils::Vector sky_coords = arc_position - arc_reference_position;
648 sky_coords = calibration.rotateToRaDec(sky_coords);
649
650 // Save the drifts in RA and DEC.
651 *driftRA = sky_coords.x;
652 *driftDEC = sky_coords.y;
653}
654
655// Computes the overall drift given the input image that was analyzed in findGuideStar().
656// Returns true if this can be done--there are already guide stars and the current detections
657// can be compared to them. In that case, fills in RADRift and DECDrift with arc-second drifts.
658// Reticle_x,y is the target position. It may be different from when we
659// originally selected the guide star, e.g. due to dithering, etc. We compute an offset from the
660// original guide-star position and the reticle position and offset all the reference star
661// positions by that.
662bool GuideStars::getDrift(double oneStarDrift, double reticle_x, double reticle_y,
663 double *RADrift, double *DECDrift)
664{
665 if (starCorrespondence.size() == 0)
666 return false;
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);
671 // Revoke multistar if we're that far away.
672 // Currently disabled by large constant.
673 constexpr double maxDriftForMultistar = 4000000.0; // arc-seconds
674
675 if (!calibrationInitialized)
676 return false;
677
678 // Don't run multistar if the 1-star approach is too far off.
679 if (oneStarDrift > maxDriftForMultistar)
680 {
681 qCDebug(KSTARS_EKOS_GUIDE) << "MultiStar: 1-star drift too high: " << oneStarDrift;
682 return false;
683 }
684 if (starCorrespondence.size() < 2)
685 {
686 qCDebug(KSTARS_EKOS_GUIDE) << "Multistar: no guide stars";
687 return false;
688 }
689
690 // The original guide star position has been offset.
691 double offset_x = reticle_x - gStar.x;
692 double offset_y = reticle_y - gStar.y;
693
694 int numStarsProcessed = 0;
695 bool guideStarProcessed = false;
696 double driftRASum = 0, driftDECSum = 0;
697 double guideStarRADrift = 0, guideStarDECDrift = 0;
698 QVector<double> raDrifts, decDrifts;
699
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)
703 {
704 const auto &star = detectedStars[i];
705 auto bg = skybackground();
706 double snr = bg.SNR(detectedStars[i].sum, detectedStars[i].numPixels);
707 // Probably should test SNR on the reference as well.
708 if (getStarMap(i) >= 0 && snr >= MIN_DRIFT_SNR)
709 {
710 auto ref = starCorrespondence.reference(getStarMap(i));
711 ref.x += offset_x;
712 ref.y += offset_y;
713 double driftRA, driftDEC;
714 computeStarDrift(star, ref, &driftRA, &driftDEC);
715 if (getStarMap(i) == starCorrespondence.guideStar())
716 {
717 guideStarRADrift = driftRA;
718 guideStarDECDrift = driftDEC;
719 guideStarProcessed = true;
720 }
721
722 raDrifts.push_back(driftRA);
723 decDrifts.push_back(driftDEC);
724 numStarsProcessed++;
725
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);
730 }
731 }
732
733 if (numStarsProcessed == 0 || (!allowMissingGuideStar && !guideStarProcessed))
734 return false;
735
736 // Filter out reference star drifts that are too different from the guide-star drift.
737 QVector<double> raDriftsKeep, decDriftsKeep;
738 for (int i = 0; i < raDrifts.size(); ++i)
739 {
740 if ((allowMissingGuideStar && !guideStarProcessed) ||
741 ((fabs(raDrifts[i] - guideStarRADrift) < 2.0) &&
742 (fabs(decDrifts[i] - guideStarDECDrift) < 2.0)))
743 {
744 driftRASum += raDrifts[i];
745 driftDECSum += decDrifts[i];
746 raDriftsKeep.push_back(raDrifts[i]);
747 decDriftsKeep.push_back(decDrifts[i]);
748 }
749 }
750 if (raDriftsKeep.empty() || raDriftsKeep.empty())
751 return false; // Shouldn't happen.
752
753 numStarsProcessed = raDriftsKeep.size();
754
755 // Generate the drift either from the median or the average of the usable reference drifts.
756 bool useMedian = true;
757 if (useMedian)
758 {
759 std::sort(raDriftsKeep.begin(), raDriftsKeep.end(), [] (double d1, double d2)
760 {
761 return d1 < d2;
762 });
763 std::sort(decDriftsKeep.begin(), decDriftsKeep.end(), [] (double d1, double d2)
764 {
765 return d1 < d2;
766 });
767 if (numStarsProcessed % 2 == 0)
768 {
769 const int middle = numStarsProcessed / 2;
770 *RADrift = ( raDriftsKeep[middle - 1] + raDriftsKeep[middle]) / 2.0;
771 *DECDrift = (decDriftsKeep[middle - 1] + decDriftsKeep[middle]) / 2.0;
772 }
773 else
774 {
775 const int middle = numStarsProcessed / 2; // because it's 0-based
776 *RADrift = raDriftsKeep[middle];
777 *DECDrift = decDriftsKeep[middle];
778 }
779 DLOG(KSTARS_EKOS_GUIDE) << "MultiStar: Drift median " << *RADrift << *DECDrift << numStarsProcessed << " of " <<
780 detectedStars.size() << "#guide" << starCorrespondence.size();
781 }
782 else
783 {
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();
788 }
789 return true;
790}
791
792void GuideStars::setCalibration(const Calibration &cal)
793{
794 calibration = cal;
795 calibrationInitialized = true;
796}
qint64 elapsed() const const
qint64 restart()
void append(QList< T > &&value)
const_reference at(qsizetype i) const const
iterator begin()
void clear()
qsizetype count() const const
bool empty() const const
iterator end()
T & first()
void push_back(parameter_type value)
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)
qsizetype size() const const
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 Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.