Kstars

starcorrespondence.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 "starcorrespondence.h"
8
9#include <math.h>
10#include "ekos_guide_debug.h"
11#include "Options.h"
12
13// Finds the star in sortedStars that's closest to x,y and within maxDistance pixels.
14// Returns the index of the closest star in sortedStars, or -1 if none satisfies the criteria.
15// sortedStars must be sorted by x, from lower x to higher x.
16// Fills distance to the pixel distance to the closest star.
17int StarCorrespondence::findClosestStar(double x, double y, const QList<Edge> sortedStars,
18 double maxDistance, double *distance) const
19{
20 if (x < -maxDistance || y < -maxDistance ||
21 x > imageWidth + maxDistance || y > imageWidth + maxDistance)
22 return -1;
23
24 const int size = sortedStars.size();
25 int startPoint = 0;
26
27 // Efficiency improvement to quickly find the start point.
28 // Increments by size/10 to find a better start point.
29 const int coarseIncrement = size / 10;
30 if (coarseIncrement > 1)
31 {
32 for (int i = 0; i < size; i += coarseIncrement)
33 {
34 if (sortedStars[i].x < x - maxDistance)
35 startPoint = i;
36 else
37 break;
38 }
39 }
40
41 // Iterate through the stars starting with the star with lowest x where x >= x - maxDistance
42 // through the star with greatest x where x <= x + maxDistance. Find the closest star to x,y.
43 int bestIndex = -1;
44 double bestSquaredDistance = maxDistance * maxDistance;
45 for (int i = startPoint; i < size; ++i)
46 {
47 auto &star = sortedStars[i];
48 if (star.x < x - maxDistance) continue;
49 if (star.x > x + maxDistance) break;
50 const double xDiff = star.x - x;
51 const double yDiff = star.y - y;
52 const double squaredDistance = xDiff * xDiff + yDiff * yDiff;
53 if (squaredDistance <= bestSquaredDistance)
54 {
55 bestIndex = i;
56 bestSquaredDistance = squaredDistance;
57 }
58 }
59 if (distance != nullptr) *distance = sqrt(bestSquaredDistance);
60 return bestIndex;
61}
62
63namespace
64{
65// Sorts stars by their x values, places the sorted stars into sortedStars.
66// sortedToOriginal is a map whose index is the index of a star in sortedStars and whose
67// value is the index of a star in stars.
68// Therefore, sortedToOrigial[a] = b, means that sortedStars[a] corresponds to stars[b].
69void sortByX(const QList<Edge> &stars, QList<Edge> *sortedStars, QVector<int> *sortedToOriginal)
70{
71 sortedStars->clear();
72 sortedToOriginal->clear();
73
74 // Sort the stars by their x-positions, lower to higher.
76 for (int i = 0; i < stars.size(); ++i)
77 rx.push_back(std::pair<int, float>(i, stars[i].x));
78 std::sort(rx.begin(), rx.end(), [](const std::pair<int, float> &a, const std::pair<int, double> &b)
79 {
80 return a.second < b.second;
81 });
82 for (int i = 0; i < stars.size(); ++i)
83 {
84 sortedStars->append(stars[ rx[i].first ]);
85 sortedToOriginal->push_back(rx[i].first);
86 }
87}
88
89} // namespace
90
91StarCorrespondence::StarCorrespondence(const QList<Edge> &stars, int guideStar)
92{
93 initialize(stars, guideStar);
94}
95
96StarCorrespondence::StarCorrespondence()
97{
98}
99
100// Initialize with the reference stars and the guide-star index.
101void StarCorrespondence::initialize(const QList<Edge> &stars, int guideStar)
102{
103 qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence initialized with " << stars.size() << guideStar;
104 references = stars;
105 guideStarIndex = guideStar;
106 float guideX = references[guideStarIndex].x;
107 float guideY = references[guideStarIndex].y;
108 guideStarOffsets.clear();
109 referenceSums.clear();
110 referenceNumPixels.clear();
111
112 // Compute the x and y offsets from the guide star to all reference stars
113 // and store them in guideStarOffsets.
114 const int numRefs = references.size();
115 for (int i = 0; i < numRefs; ++i)
116 {
117 const auto &ref = references[i];
118 guideStarOffsets.push_back(Offsets(ref.x - guideX, ref.y - guideY));
119 referenceSums.push_back(ref.sum);
120 referenceNumPixels.push_back(ref.numPixels);
121 }
122
123 initializeAdaptation();
124
125 initialized = true;
126}
127
128void StarCorrespondence::reset()
129{
130 references.clear();
131 guideStarOffsets.clear();
132 referenceSums.clear();
133 referenceNumPixels.clear();
134 initialized = false;
135}
136
137int StarCorrespondence::findInternal(const QList<Edge> &stars, double maxDistance, QVector<int> *starMap,
138 int guideStarIndex, const QVector<Offsets> &offsets,
139 int *numFound, int *numNotFound, double minFraction) const
140{
141 // This is the cost of not finding one of the reference stars.
142 constexpr double missingRefStarCost = 100;
143 // This weight multiplies the distance between a reference star position and the closest star in stars.
144 constexpr double distanceWeight = 1.0;
145
146 // Initialize all stars to not-corresponding to any reference star.
147 *starMap = QVector<int>(stars.size(), -1);
148
149 // We won't accept a solution worse than bestCost.
150 // In the default case, we need to find about half the reference stars.
151 // Another possibility is a constraint on the input stars being mapped
152 // E.g. that the more likely solution is one where the stars are close to the references.
153 // This can be an issue if the number of input stars is way less than the number of reference stars
154 // but in that case we can fail and go to the default star-finding algorithm.
155 int bestCost = offsets.size() * missingRefStarCost * (1 - minFraction);
156 // Note that the above implies that if stars.size() < offsets.size() * minFraction
157 // then it is impossible to succeed.
158 if (stars.size() < minFraction * offsets.size())
159 return -1;
160
161 // Assume the guide star corresponds to each of the stars.
162 // Score the assignment, pick the best, and then assign the rest.
163 const int numStars = stars.size();
164 int bestStarIndex = -1, bestNumFound = 0, bestNumNotFound = 0;
165 for (int starIndex = 0; starIndex < numStars; ++starIndex)
166 {
167 const float starX = stars[starIndex].x;
168 const float starY = stars[starIndex].y;
169
170 double cost = 0.0;
171 QVector<int> mapping(stars.size(), -1);
172 int numFound = 0, numNotFound = 0;
173 for (int offsetIndex = 0; offsetIndex < offsets.size(); ++offsetIndex)
174 {
175 // Of course, the guide star has offsets 0 to itself.
176 if (offsetIndex == guideStarIndex) continue;
177
178 // We're already worse than the best cost. No need to search any more.
179 if (cost > bestCost) break;
180
181 // Look for an input star at the offset position.
182 // Note the index value returned is the sorted index, which is converted back later.
183 const auto &offset = offsets[offsetIndex];
184 double distance;
185 const int closestIndex = findClosestStar(starX + offset.x, starY + offset.y,
186 stars, maxDistance, &distance);
187 if (closestIndex < 0)
188 {
189 // This reference star position had no corresponding input star.
190 cost += missingRefStarCost;
191 numNotFound++;
192 continue;
193 }
194 numFound++;
195
196 // If starIndex is the star that corresponds to guideStarIndex, then
197 // stars[index] corresponds to references[offsetIndex]
198 mapping[closestIndex] = offsetIndex;
199 cost += distance * distanceWeight;
200 }
201 if (cost < bestCost)
202 {
203 bestCost = cost;
204 bestStarIndex = starIndex;
205 bestNumFound = numFound;
206 bestNumNotFound = numNotFound;
207
208 *starMap = QVector<int>(stars.size(), -1);
209 for (int i = 0; i < mapping.size(); ++i)
210 (*starMap)[i] = mapping[i];
211 (*starMap)[starIndex] = guideStarIndex;
212 }
213 }
214 *numFound = bestNumFound;
215 *numNotFound = bestNumNotFound;
216 return bestStarIndex;
217}
218
219// Converts offsets that are from the point-of-view of the guidestar, to offsets from targetStar.
220void StarCorrespondence::makeOffsets(const QVector<Offsets> &offsets, QVector<Offsets> *targetOffsets,
221 int targetStar) const
222{
223 targetOffsets->resize(offsets.size());
224 int xOffset = offsets[targetStar].x;
225 int yOffset = offsets[targetStar].y;
226 for (int i = 0; i < offsets.size(); ++i)
227 {
228 if (i == targetStar)
229 {
230 (*targetOffsets)[i] = Offsets(0, 0);
231 }
232 else
233 {
234 const double x = offsets[i].x - xOffset;
235 const double y = offsets[i].y - yOffset;
236 (*targetOffsets)[i] = Offsets(x, y);
237 }
238 }
239
240}
241
242// We create an imaginary star from the ones we did find.
243Edge StarCorrespondence::inventStarPosition(const QList<Edge> &stars, const QVector<int> &starMap,
244 const QVector<Offsets> &offsets, const Offsets &offset) const
245{
246 Edge inventedStar;
247 inventedStar.invalidate();
248
249 QVector<double> xPositions, yPositions;
250 float refSum = 0, origSum = 0, refNumPixels = 0, origNumPixels = 0;
251 for (int i = 0; i < starMap.size(); ++i)
252 {
253 int reference = starMap[i];
254 if (reference >= 0)
255 {
256 xPositions.push_back(stars[i].x - offsets[reference].x);
257 yPositions.push_back(stars[i].y - offsets[reference].y);
258
259 // These are used to help invent the SNR.
260 refSum += stars[i].sum;
261 refNumPixels += stars[i].numPixels;
262 origSum += referenceSums[reference];
263 origNumPixels += referenceNumPixels[reference];
264 }
265 }
266 float inventedSum = 0;
267 float inventedNumPixels = 0;
268 if (origSum > 0 && origNumPixels > 0)
269 {
270 // If possible, interpolate to estimate the invented star's flux and number of pixels.
271 inventedSum = (refSum / origSum) * referenceSums[guideStarIndex];
272 inventedNumPixels = (refNumPixels / origNumPixels) * referenceNumPixels[guideStarIndex];
273 }
274
275 if (xPositions.size() == 0)
276 return inventedStar;
277
278 // Compute the median x and y values. After gathering the values above,
279 // we sort them and use the middle positions.
280 std::sort(xPositions.begin(), xPositions.end(), [] (double d1, double d2)
281 {
282 return d1 < d2;
283 });
284 std::sort(yPositions.begin(), yPositions.end(), [] (double d1, double d2)
285 {
286 return d1 < d2;
287 });
288 const int num = xPositions.size();
289 double xVal = 0, yVal = 0;
290 if (num % 2 == 0)
291 {
292 const int middle = num / 2;
293 xVal = (xPositions[middle - 1] + xPositions[middle]) / 2.0;
294 yVal = (yPositions[middle - 1] + yPositions[middle]) / 2.0;
295 }
296 else
297 {
298 const int middle = num / 2; // because it's 0-based
299 xVal = xPositions[middle];
300 yVal = yPositions[middle];
301 }
302 inventedStar.x = xVal - offset.x;
303 inventedStar.y = yVal - offset.y;
304 inventedStar.sum = inventedSum;
305 inventedStar.numPixels = inventedNumPixels;
306 return inventedStar;
307}
308
309namespace
310{
311// This utility converts the starMap to refer to indexes in the original star QVector
312// instead of the sorted-by-x version.
313void unmapStarMap(const QVector<int> &sortedStarMap, const QVector<int> &sortedToOriginal, QVector<int> *starMap)
314{
315 for (int i = 0; i < sortedStarMap.size(); ++i)
316 (*starMap)[sortedToOriginal[i]] = sortedStarMap[i];
317}
318}
319
320Edge StarCorrespondence::find(const QList<Edge> &stars, double maxDistance,
321 QVector<int> *starMap, bool adapt, double minFraction)
322{
323 m_NumReferencesFound = 0;
324 *starMap = QVector<int>(stars.size(), -1);
325 Edge foundStar;
326 foundStar.invalidate();
327 if (!initialized) return foundStar;
328 int numFound = 0, numNotFound = 0;
329
330 // findClosestStar needs an input with stars sorted by their x.
331 // Do this outside of the loops.
332 QList<Edge> sortedStars;
333 QVector<int> sortedToOriginal;
334 sortByX(stars, &sortedStars, &sortedToOriginal);
335
336 const bool alwaysInvent = Options::alwaysInventGuideStar() &&
337 stars.size() >= minFraction * guideStarOffsets.size();
338 if (!alwaysInvent && Options::alwaysInventGuideStar())
339 qCDebug(KSTARS_EKOS_GUIDE)
340 << QString("Could not use AlwaysInvent--too few stars: %1 < %2 * %3 (%4)").arg(stars.size())
341 .arg(minFraction, 0, 'f', 2).arg(guideStarOffsets.size()).arg(minFraction * guideStarOffsets.size(), 0, 'f', 1);
342 const bool canInvent = allowMissingGuideStar && stars.size() >= minFraction * guideStarOffsets.size();
343
344 QVector<int> sortedStarMap;
345 int bestStarIndex =
346 alwaysInvent ? -1 : findInternal(sortedStars, maxDistance, &sortedStarMap, guideStarIndex,
347 guideStarOffsets, &numFound, &numNotFound, minFraction);
348
349 if (!alwaysInvent && bestStarIndex > -1)
350 {
351 // Convert back to the unsorted index value.
352 bestStarIndex = sortedToOriginal[bestStarIndex];
353 unmapStarMap(sortedStarMap, sortedToOriginal, starMap);
354
355 foundStar = stars[bestStarIndex];
356 qCDebug(KSTARS_EKOS_GUIDE)
357 << "StarCorrespondence found guideStar at " << bestStarIndex << "found/not"
358 << numFound << numNotFound;
359 m_NumReferencesFound = numFound;
360 }
361 else if (alwaysInvent || canInvent)
362 {
363 // If we didn't find a good solution in findInternal() above (with alwaysInvent false), perhaps
364 // the guide star was not detected. We try to get around that by getting an estimate of its position
365 // from many reference stars. If alwaysInvent is true, we always wind up here, as that estimate is
366 // likely more robust anyway than consistently finding the one guide star.
367 if (!alwaysInvent)
368 qCDebug(KSTARS_EKOS_GUIDE)
369 << "Trying to invent. StarCorrespondence did NOT find guideStar. Found/not"
370 << numFound << numNotFound << "minFraction" << minFraction << "maxDistance" << maxDistance;
371 else
372 qCDebug(KSTARS_EKOS_GUIDE)
373 << "Using AlwaysInvent. #stars" << stars.size() << "minFraction" << minFraction << "maxDistance" << maxDistance;
374
375 int bestNumFound = 0;
376 int bestNumNotFound = 0;
377 Edge bestInvented;
378 bestInvented.invalidate();
379 // For each reference star, pretend it was the guide star by using makeOffsets() to convert the
380 // original guide-star's offsets to new offsets relative to that reference star.
381 // Then, using findInternal(), map the detected stars to the reference stars.
382 // If we were successful in doing that (i.e. the reference star was actually detected)
383 // then we should know the positions of all (detected) reference stars.
384 // We then use inventStarPosition() to make an "invented guide star position" by reversing the
385 // offsets and translating each found ref star position back to a guide-star position and finding
386 // the median position of those. We only need to do this for one (actually detected) reference star.
387 for (int gStarIndex = 0; gStarIndex < guideStarOffsets.size(); gStarIndex++)
388 {
389 if (!alwaysInvent && gStarIndex == guideStarIndex)
390 continue;
391 QVector<Offsets> gStarOffsets;
392 makeOffsets(guideStarOffsets, &gStarOffsets, gStarIndex);
393 QVector<int> newStarMap;
394 int detectedStarIndex = findInternal(sortedStars, maxDistance, &newStarMap,
395 gStarIndex, gStarOffsets,
396 &numFound, &numNotFound, minFraction);
397 if (detectedStarIndex >= 0 && numFound > bestNumFound)
398 {
399 Edge invented = inventStarPosition(sortedStars, newStarMap, gStarOffsets,
400 guideStarOffsets[gStarIndex]);
401 if (invented.x < 0 || invented.y < 0)
402 continue;
403
404 bestInvented = invented;
405 bestNumFound = numFound;
406 bestNumNotFound = numNotFound;
407 sortedStarMap = newStarMap;
408
409 // If enough of the references were found to reliably estimate the guide-star position
410 // then we can break out of the loop. Nothing special about 7, just a guess.
411 if (numNotFound <= 1 || bestNumFound >= 7)
412 break;
413 }
414 }
415 if (bestNumFound > 0)
416 {
417 // Convert back to the unsorted index value.
418 unmapStarMap(sortedStarMap, sortedToOriginal, starMap);
419 qCDebug(KSTARS_EKOS_GUIDE)
420 << "StarCorrespondence found guideStar (invented) at "
421 << bestInvented.x << bestInvented.y << "found/not" << bestNumFound << bestNumNotFound;
422 m_NumReferencesFound = bestNumFound;
423
424 if (alwaysInvent && adapt)
425 adaptOffsets(stars, *starMap, bestInvented.x, bestInvented.y);
426
427 return bestInvented;
428 }
429 else qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence could not invent guideStar.";
430 }
431 else qCDebug(KSTARS_EKOS_GUIDE) << "StarCorrespondence could not find guideStar.";
432
433 if (adapt && (bestStarIndex != -1))
434 adaptOffsets(stars, *starMap, foundStar.x, foundStar.y);
435 return foundStar;
436}
437
438// Compute the alpha cooeficient for a simple IIR filter used in adaptOffsets()
439// that causes the offsets to be, roughtly, the moving average of the last timeConstant samples.
440// See the discussion here:
441// https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order-iir-ar-filter-approximation-to-a-moving-average-f
442void StarCorrespondence::initializeAdaptation()
443{
444 // Approximately average the last 25 samples.
445 const double timeConstant = 25.0;
446 alpha = 1.0 / pow(timeConstant, 0.865);
447
448 // Don't let the adapted offsets move far from the original ones.
449 originalGuideStarOffsets = guideStarOffsets;
450}
451
452void StarCorrespondence::adaptOffsets(const QList<Edge> &stars, const QVector<int> &starMap, double x, double y)
453{
454 const int numStars = stars.size();
455 if (starMap.size() != numStars)
456 {
457 qCDebug(KSTARS_EKOS_GUIDE) << "Adapt error: StarMap size != stars.size()" << starMap.size() << numStars;
458 return;
459 }
460
461 for (int i = 0; i < numStars; ++i)
462 {
463 // We don't adapt if the ith star doesn't correspond to any of the references,
464 // or if it corresponds to the guide star (whose offsets are, of course, always 0).
465 const int refIndex = starMap[i];
466 if (refIndex == -1 || refIndex == guideStarIndex)
467 continue;
468
469 if ((refIndex >= references.size()) || (refIndex < 0))
470 {
471 qCDebug(KSTARS_EKOS_GUIDE) << "Adapt error: bad refIndex[" << i << "] =" << refIndex;
472 return;
473 }
474 // Adapt the x and y offsets using the following IIR filter:
475 // output[t] = alpha * offset[t] + (1-alpha) * output[t-1]
476 const double xOffset = stars[i].x - x;
477 const double yOffset = stars[i].y - y;
478 const double currentXOffset = guideStarOffsets[refIndex].x;
479 const double currentYOffset = guideStarOffsets[refIndex].y;
480 const double newXOffset = alpha * xOffset + (1 - alpha) * currentXOffset;
481 const double newYOffset = alpha * yOffset + (1 - alpha) * currentYOffset;
482
483 // The adaptation is meant for small movements of at most a few pixels.
484 // We don't want it to move enough to find a new star.
485 constexpr double maxMovement = 2.5; // pixels
486 if ((fabs(newXOffset - originalGuideStarOffsets[refIndex].x) < maxMovement) &&
487 (fabs(newYOffset - originalGuideStarOffsets[refIndex].y) < maxMovement))
488 {
489 guideStarOffsets[refIndex].x = newXOffset;
490 guideStarOffsets[refIndex].y = newYOffset;
491 }
492 }
493}
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
void append(QList< T > &&value)
iterator begin()
void clear()
iterator end()
void push_back(parameter_type value)
void resize(qsizetype size)
qsizetype size() const const
QString arg(Args &&... args) const const
float x() const const
float y() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Dec 6 2024 12:09:52 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.