# Kstars

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