Kstars

focusalgorithms.cpp
1/*
2 SPDX-FileCopyrightText: 2019 Hy Murveit <hy-1@murveit.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5*/
6
7#include "focusalgorithms.h"
8#include "curvefit.h"
9
10#include "polynomialfit.h"
11#include <QVector>
12#include "kstars.h"
13#include "fitsviewer/fitsstardetector.h"
14#include "focusstars.h"
15#include <ekos_focus_debug.h>
16
17namespace Ekos
18{
19/**
20 * @class LinearFocusAlgorithm
21 * @short Autofocus algorithm that linearly samples HFR values.
22 *
23 * @author Hy Murveit
24 * @version 1.0
25 */
26class LinearFocusAlgorithm : public FocusAlgorithmInterface
27{
28 public:
29
30 // Constructor initializes a linear autofocus algorithm, starting at some initial position,
31 // sampling HFR values, decreasing the position by some step size, until the algorithm believe's
32 // it's seen a minimum. It then either:
33 // a) tries to find that minimum in a 2nd pass. This is the LINEAR algorithm, or
34 // b) moves to the minimum. This is the LINEAR 1 PASS algorithm.
35 LinearFocusAlgorithm(const FocusParams &params);
36
37 // After construction, this should be called to get the initial position desired by the
38 // focus algorithm. It returns the initial position passed to the constructor if
39 // it has no movement request.
40 int initialPosition() override
41 {
42 return requestedPosition;
43 }
44
45 // Pass in the measurement for the last requested position. Returns the position for the next
46 // requested measurement, or -1 if the algorithm's done or if there's an error.
47 // If stars is not nullptr, then the relativeHFR scheme is used to modify the HFR value.
48 int newMeasurement(int position, double value, const double starWeight, const QList<Edge*> *stars) override;
49
50 FocusAlgorithmInterface *Copy() override;
51
52 void getMeasurements(QVector<int> *pos, QVector<double> *val, QVector<double> *sds) const override
53 {
54 *pos = positions;
55 *val = values;
56 *sds = weights;
57 }
58
59 void getPass1Measurements(QVector<int> *pos, QVector<double> *val, QVector<double> *sds, QVector<bool> *out) const override
60 {
61 *pos = pass1Positions;
62 *val = pass1Values;
63 *sds = pass1Weights;
64 *out = pass1Outliers;
65 }
66
67 double latestValue() const override
68 {
69 if (values.size() > 0)
70 return values.last();
71 return -1;
72 }
73
74 bool isInFirstPass() const override
75 {
76 if (params.focusAlgorithm == Focus::FOCUS_LINEAR || params.focusAlgorithm == Focus::FOCUS_LINEAR1PASS)
77 return inFirstPass;
78 else
79 return true;
80 }
81
82 QString getTextStatus(double R2 = 0) const override;
83
84 int currentStep() const override
85 {
86 return numSteps;
87 }
88
89 private:
90
91 // Called in newMeasurement. Sets up the next iteration.
92 int completeIteration(int step, bool foundFit, double minPos, double minVal);
93
94 // Does the bookkeeping for the final focus solution.
95 int setupSolution(int position, double value, double sigma);
96
97 // Perform the linear walk for FOCUS_WALK_FIXED_STEPS and FOCUS_WALK_CFZ_SHUFFLE
98 int linearWalk(int position, double value, const double starWeight);
99
100 // Calculate the curve fitting goal based on the algorithm parameters and progress
101 CurveFitting::FittingGoal getGoal(int numSteps);
102
103 // Get the minimum number of datapoints before attempting the first curve fit
104 int getCurveMinPoints();
105
106 // Analyze data for donuts are remove
107 void removeDonuts();
108
109 // Calc the next step size for Linear1Pass for FOCUS_WALK_FIXED_STEPS and FOCUS_WALK_CFZ_SHUFFLE
110 int getNextStepSize();
111
112 // Called when we've found a solution, e.g. the HFR value is within tolerance of the desired value.
113 // It it returns true, then it's decided tht we should try one more sample for a possible improvement.
114 // If it returns false, then "play it safe" and use the current position as the solution.
115 bool setupPendingSolution(int position);
116
117 // Determines the desired focus position for the first sample.
118 void computeInitialPosition();
119
120 // Get the first step focus position
121 int getFirstPosition();
122
123 // Sets the internal state for re-finding the minimum, and returns the requested next
124 // position to sample.
125 // Called by Linear. L1P calls finishFirstPass instead
126 int setupSecondPass(int position, double value, double margin = 2.0);
127
128 // Called by L1P to finish off the first pass, setting state variables.
129 int finishFirstPass(int position, double value);
130
131 // Peirce's criterion
132 // Returns the squared threshold error deviation for outlier identification
133 // using Peirce's criterion based on Gould's methodology.
134 // Arguments:
135 // N = total number of observations
136 // n = max number of outliers to be removed
137 // m = number of model unknowns
138 // Returns: squared error threshold (x2)
139 double peirce_criterion(double N, double n, double m);
140
141 // Used in the 2nd pass. Focus is getting worse. Requires several consecutive samples getting worse.
142 bool gettingWorse();
143
144 // If one of the last 2 samples are as good or better than the 2nd best so far, return true.
145 bool bestSamplesHeuristic();
146
147 // Adds to the debug log a line summarizing the result of running this algorithm.
148 void debugLog();
149
150 // Use the detailed star measurements to adjust the HFR value.
151 // See comment above method definition.
152 double relativeHFR(double origHFR, const QList<Edge*> *stars);
153
154 // Calculate the standard deviation (=sigma) of the star HFRs
155 double calculateStarSigma(const bool useWeights, const QList<Edge*> *stars);
156
157 // Used to time the focus algorithm.
158 QElapsedTimer stopWatch;
159
160 // A vector containing the HFR values sampled by this algorithm so far.
161 QVector<double> values;
162 QVector<double> pass1Values;
163 // A vector containing star weights
164 QVector<double> weights;
165 QVector<double> pass1Weights;
166 // A vector containing the focus positions corresponding to the HFR values stored above.
167 QVector<int> positions;
168 QVector<int> pass1Positions;
169 // A vector containing whether or not the datapoint has been identified as an outlier.
170 QVector<bool> pass1Outliers;
171
172 // Focus position requested by this algorithm the previous step.
173 int requestedPosition;
174 // The position where the first pass is begun. Usually requestedPosition unless there's a restart.
175 int passStartPosition;
176 // Number of iterations processed so far. Used to see it doesn't exceed params.maxIterations.
177 int numSteps;
178 // The best value in the first pass. The 2nd pass attempts to get within
179 // tolerance of this value.
180 double firstPassBestValue;
181 // The position of the minimum found in the first pass.
182 int firstPassBestPosition;
183 // The sampling interval--the recommended number of focuser steps moved inward each iteration
184 // of the first pass.
185 int stepSize;
186 // The minimum focus position to use. Computed from the focuser limits and maxTravel.
187 int minPositionLimit;
188 // The maximum focus position to use. Computed from the focuser limits and maxTravel.
189 int maxPositionLimit;
190 // Counter for the number of times a v-curve minimum above the current position was found,
191 // which implies the initial focus sweep has passed the minimum and should be terminated.
192 int numPolySolutionsFound;
193 // Counter for the number of times a v-curve minimum above the passStartPosition was found,
194 // which implies the sweep should restart at a higher position.
195 int numRestartSolutionsFound;
196 // The index (into values and positions) when the most recent 2nd pass started.
197 int secondPassStartIndex;
198 // True if performing the first focus sweep.
199 bool inFirstPass;
200 // True if the 2nd pass has found a solution, and it's now optimizing the solution.
201 bool solutionPending;
202 // When we're near done, but the HFR just got worse, we may retry the current position
203 // in case the HFR value was noisy.
204 int retryNumber = 0;
205
206 // Keep the solution-pending position/value for the status messages.
207 double solutionPendingValue = 0;
208 int solutionPendingPosition = 0;
209
210 // RelativeHFR variables
211 bool relativeHFREnabled = false;
212 QVector<FocusStars> starLists;
213 double bestRelativeHFR = 0;
214 int bestRelativeHFRIndex = 0;
215};
216
217// Copies the object. Used in testing to examine alternate possible inputs given
218// the current state.
219FocusAlgorithmInterface *LinearFocusAlgorithm::Copy()
220{
221 LinearFocusAlgorithm *alg = new LinearFocusAlgorithm(params);
222 *alg = *this;
223 return dynamic_cast<FocusAlgorithmInterface*>(alg);
224}
225
226FocusAlgorithmInterface *MakeLinearFocuser(const FocusAlgorithmInterface::FocusParams &params)
227{
228 return new LinearFocusAlgorithm(params);
229}
230
231LinearFocusAlgorithm::LinearFocusAlgorithm(const FocusParams &focusParams)
232 : FocusAlgorithmInterface(focusParams)
233{
234 stopWatch.start();
235 // These variables don't get reset if we restart the algorithm.
236 numSteps = 0;
237 maxPositionLimit = std::min(params.maxPositionAllowed, params.startPosition + params.maxTravel);
238 minPositionLimit = std::max(params.minPositionAllowed, params.startPosition - params.maxTravel);
239 computeInitialPosition();
240}
241
242QString LinearFocusAlgorithm::getTextStatus(double R2) const
243{
244 if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
245 {
246 if (done)
247 {
248 if (focusSolution > 0)
249 return QString("Solution: %1").arg(focusSolution);
250 else
251 return "Failed";
252 }
253 if (solutionPending)
254 return QString("Pending: %1, %2").arg(solutionPendingPosition).arg(solutionPendingValue, 0, 'f', 2);
255 if (inFirstPass)
256 return "1st Pass";
257 else
258 return QString("2nd Pass. 1st: %1, %2")
259 .arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 2);
260 }
261 else // Linear 1 Pass
262 {
263 QString text = "L1P";
264 if (params.filterName != "")
265 text.append(QString(" [%1]").arg(params.filterName));
266 if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
267 text.append(": Quadratic");
268 else if (params.curveFit == CurveFitting::FOCUS_HYPERBOLA)
269 {
270 text.append(": Hyperbola");
271 text.append(params.useWeights ? " (W)" : " (U)");
272 }
273 else if (params.curveFit == CurveFitting::FOCUS_PARABOLA)
274 {
275 text.append(": Parabola");
276 text.append(params.useWeights ? " (W)" : " (U)");
277 }
278 else if (params.curveFit == CurveFitting::FOCUS_2DGAUSSIAN)
279 {
280 text.append(": Gaussian");
281 text.append(params.useWeights ? " (W)" : " (U)");
282 }
283
284 if (inFirstPass)
285 return text;
286 else if (!done)
287 return text.append(". Moving to Solution");
288 else
289 {
290 if (focusSolution > 0)
291 {
292 if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
293 return text.append(" Solution: %1").arg(focusSolution);
294 else
295 // Add R2 to 2 decimal places. Round down to be conservative
296 return text.append(" Solution: %1, R²=%2").arg(focusSolution).arg(trunc(R2 * 100.0) / 100.0, 0, 'f', 2);
297 }
298 else
299 return text.append(" Failed");
300 }
301 }
302}
303
304void LinearFocusAlgorithm::computeInitialPosition()
305{
306 // These variables get reset if the algorithm is restarted.
307 stepSize = params.initialStepSize;
308 inFirstPass = true;
309 solutionPending = false;
310 retryNumber = 0;
311 firstPassBestValue = -1;
312 firstPassBestPosition = 0;
313 numPolySolutionsFound = 0;
314 numRestartSolutionsFound = 0;
315 secondPassStartIndex = -1;
316
317 qCDebug(KSTARS_EKOS_FOCUS)
318 << QString("Linear: Travel %1 initStep %2 pos %3 min %4 max %5 maxIters %6 tolerance %7 minlimit %8 maxlimit %9 init#steps %10 algo %11 backlash %12 curvefit %13 useweights %14")
319 .arg(params.maxTravel).arg(params.initialStepSize).arg(params.startPosition).arg(params.minPositionAllowed)
320 .arg(params.maxPositionAllowed).arg(params.maxIterations).arg(params.focusTolerance).arg(minPositionLimit)
321 .arg(maxPositionLimit).arg(params.initialOutwardSteps).arg(params.focusAlgorithm).arg(params.backlash)
322 .arg(params.curveFit).arg(params.useWeights);
323
324 requestedPosition = std::min(maxPositionLimit, getFirstPosition());
325 passStartPosition = requestedPosition;
326 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: initialPosition %1 sized %2")
327 .arg(requestedPosition).arg(params.initialStepSize);
328}
329
330// Calculate the first position relative to the start position
331// Generally this is just the number of outward points * step size.
332// For LINEAR1PASS most walks have a constant step size so again it is the number of outward points * step size.
333// For LINEAR1PASS and FOCUS_WALK_CFZ_SHUFFLE the middle 3 or 4 points are separated by half step size
334int LinearFocusAlgorithm::getFirstPosition()
335{
336 if (params.focusAlgorithm != Focus::FOCUS_LINEAR1PASS)
337 return static_cast<int>(params.startPosition + params.initialOutwardSteps * params.initialStepSize);
338
339 int firstPosition;
340 double outSteps, numFullStepsOut, numHalfStepsOut;
341
342 switch (params.focusWalk)
343 {
344 case Focus::FOCUS_WALK_CLASSIC:
345 firstPosition = static_cast<int>(params.startPosition + params.initialOutwardSteps * params.initialStepSize);
346 break;
347
348 case Focus::FOCUS_WALK_FIXED_STEPS:
349 outSteps = (params.numSteps - 1) / 2.0f;
350 firstPosition = params.startPosition + (outSteps * params.initialStepSize);
351 break;
352
353 case Focus::FOCUS_WALK_CFZ_SHUFFLE:
354 if (params.numSteps % 2)
355 {
356 // Odd number of steps
357 numHalfStepsOut = 1.0f;
358 numFullStepsOut = ((params.numSteps - 1) / 2.0f) - numHalfStepsOut;
359 }
360 else
361 {
362 // Even number of steps
363 numHalfStepsOut = 1.5f;
364 numFullStepsOut = (params.numSteps / 2.0f) - 2.0f;
365 }
366 firstPosition = params.startPosition + (numFullStepsOut * params.initialStepSize) + (numHalfStepsOut *
367 (params.initialStepSize / 2.0f));
368
369 break;
370
371 default:
372 firstPosition = static_cast<int>(params.startPosition + params.initialOutwardSteps * params.initialStepSize);
373 break;
374 }
375
376 return firstPosition;
377}
378
379// The Problem:
380// The HFR values passed in to these focus algorithms are simply the mean or median HFR values
381// for all stars detected (with some other constraints). However, due to randomness in the
382// star detection scheme, two sets of star measurements may use different sets of actual stars to
383// compute their mean HFRs. This adds noise to determining best focus (e.g. including a
384// wider star in one set but not the other will tend to increase the HFR for that star set).
385// relativeHFR tries to correct for this by comparing two sets of stars only using their common stars.
386//
387// The Solution:
388// We maintain a reference set of stars, along with the HFR computed for those reference stars (bestRelativeHFR).
389// We find the common set of stars between an input star set and those reference stars. We compute
390// HFRs for the input stars in that common set, and for the reference common set as well.
391// The ratio of those common-set HFRs multiply bestRelativeHFR to generate the relative HFR for this
392// input star set--that is, it's the HFR "relative to the reference set". E.g. if the common-set HFR
393// for the input stars is twice worse than that from the reference common stars, then the relative HFR
394// would be 2 * bestRelativeHFR.
395// The reference set is the best HFR set of stars seen so far, but only from the 1st pass.
396// Thus, in the 2nd pass, the reference stars would be the best HFR set from the 1st pass.
397//
398// In unusual cases, e.g. when the common set can't be computed, we just return the input origHFR.
399double LinearFocusAlgorithm::relativeHFR(double origHFR, const QList<Edge*> *stars)
400{
401 constexpr double minHFR = 0.3;
402 if (origHFR < minHFR) return origHFR;
403
404 const int currentIndex = values.size();
405 const bool isFirstSample = (currentIndex == 0);
406 double relativeHFR = origHFR;
407
408 if (isFirstSample && stars != nullptr)
409 relativeHFREnabled = true;
410
411 if (relativeHFREnabled && stars == nullptr)
412 {
413 // Error--we have been getting relativeHFR stars, and now we're not.
414 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Error: Inconsistent relativeHFR, disabling.");
415 relativeHFREnabled = false;
416 }
417
418 if (!relativeHFREnabled)
419 return origHFR;
420
421 if (starLists.size() != positions.size())
422 {
423 // Error--inconsistent data structures!
424 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Error: Inconsistent relativeHFR data structures, disabling.");
425 relativeHFREnabled = false;
426 return origHFR;
427 }
428
429 // Copy the input stars.
430 // FocusStars enables us to measure HFR relatively consistently,
431 // by using the same stars when comparing two measurements.
432 constexpr int starDistanceThreshold = 5; // Max distance (pixels) for two positions to be considered the same star.
433 FocusStars fs(*stars, starDistanceThreshold);
434 starLists.push_back(fs);
435 auto &currentStars = starLists.back();
436
437 if (isFirstSample)
438 {
439 relativeHFR = currentStars.getHFR();
440 if (relativeHFR <= 0)
441 {
442 // Fall back to the simpler scheme.
443 relativeHFREnabled = false;
444 return origHFR;
445 }
446 // 1st measurement. By definition this is the best HFR.
447 bestRelativeHFR = relativeHFR;
448 bestRelativeHFRIndex = currentIndex;
449 }
450 else
451 {
452 // HFR computed relative to the best measured so far.
453 auto &bestStars = starLists[bestRelativeHFRIndex];
454 double hfr = currentStars.relativeHFR(bestStars, bestRelativeHFR);
455 if (hfr > 0)
456 {
457 relativeHFR = hfr;
458 }
459 else
460 {
461 relativeHFR = currentStars.getHFR();
462 if (relativeHFR <= 0)
463 return origHFR;
464 }
465
466 // In the 1st pass we compute the current HFR relative to the best HFR measured yet.
467 // In the 2nd pass we compute the current HFR relative to the best HFR in the 1st pass.
468 if (inFirstPass && relativeHFR <= bestRelativeHFR)
469 {
470 bestRelativeHFR = relativeHFR;
471 bestRelativeHFRIndex = currentIndex;
472 }
473 }
474
475 qCDebug(KSTARS_EKOS_FOCUS) << QString("RelativeHFR: orig %1 computed %2 relative %3")
476 .arg(origHFR).arg(currentStars.getHFR()).arg(relativeHFR);
477
478 return relativeHFR;
479}
480
481int LinearFocusAlgorithm::newMeasurement(int position, double value, const double starWeight, const QList<Edge*> *stars)
482{
483 if (params.focusAlgorithm == Focus::FOCUS_LINEAR1PASS && (params.focusWalk == Focus::FOCUS_WALK_FIXED_STEPS
484 || params.focusWalk == Focus::FOCUS_WALK_CFZ_SHUFFLE))
485 return linearWalk(position, value, starWeight);
486
487 double minPos = -1.0, minVal = 0;
488 bool foundFit = false;
489
490 double origValue = value;
491 // For QUADRATIC continue to use the relativeHFR functionality
492 // For HYPERBOLA and PARABOLA the stars used for the HDR calculation and the sigma calculation
493 // should be the same. For now, we will use the full set of stars and therefore not adjust the HFR
494 if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
495 value = relativeHFR(value, stars);
496
497 // For LINEAR 1 PASS don't bother with a full 2nd pass just jump to the solution
498 // Do the step out and back in to deal with backlash
499 if (params.focusAlgorithm == Focus::FOCUS_LINEAR1PASS && !inFirstPass)
500 {
501 return setupSolution(position, value, starWeight);
502 }
503
504 int thisStepSize = stepSize;
505 ++numSteps;
506 if (inFirstPass)
507 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: step %1, newMeasurement(%2, %3 -> %4, %5)")
508 .arg(numSteps).arg(position).arg(origValue).arg(value)
509 .arg(stars == nullptr ? 0 : stars->size());
510 else
511 qCDebug(KSTARS_EKOS_FOCUS)
512 << QString("Linear: step %1, newMeasurement(%2, %3 -> %4, %5) 1stBest %6 %7").arg(numSteps)
513 .arg(position).arg(origValue).arg(value).arg(stars == nullptr ? 0 : stars->size())
514 .arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
515
516 const int LINEAR_POSITION_TOLERANCE = params.initialStepSize;
517 if (abs(position - requestedPosition) > LINEAR_POSITION_TOLERANCE)
518 {
519 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error didn't get the requested position");
520 return requestedPosition;
521 }
522 // Have we already found a solution?
523 if (focusSolution != -1)
524 {
525 doneString = i18n("Called newMeasurement after a solution was found.");
526 failCode = Ekos::FOCUS_FAIL_INTERNAL;
527 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
528 debugLog();
529 return -1;
530 }
531
532 // Store the sample values.
533 positions.push_back(position);
534 values.push_back(value);
535 weights.push_back(starWeight);
536 if (inFirstPass)
537 {
538 pass1Positions.push_back(position);
539 pass1Values.push_back(value);
540 pass1Weights.push_back(starWeight);
541 pass1Outliers.push_back(false);
542 }
543
544 // If we're in the 2nd pass and either
545 // - the current value is within tolerance, or
546 // - we're optimizing because we've previously found a within-tolerance solution,
547 // then setupPendingSolution decides whether to continue optimizing or to complete.
548 if (solutionPending ||
549 (!inFirstPass && (value < firstPassBestValue * (1.0 + params.focusTolerance))))
550 {
551 if (setupPendingSolution(position))
552 // We can continue to look a little further.
553 return completeIteration(retryNumber > 0 ? 0 : stepSize, false, -1.0f, -1.0f);
554 else
555 // Finish now
556 return setupSolution(position, value, starWeight);
557 }
558
559 if (inFirstPass)
560 {
561 constexpr int kNumPolySolutionsRequired = 2;
562 constexpr int kNumRestartSolutionsRequired = 3;
563 //constexpr double kDecentValue = 3.0;
564
565 if (values.size() >= getCurveMinPoints())
566 {
567 if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
568 {
569 PolynomialFit fit(2, positions, values);
570 foundFit = fit.findMinimum(position, 0, params.maxPositionAllowed, &minPos, &minVal);
571 if (!foundFit)
572 {
573 // I've found that the first sample can be odd--perhaps due to backlash.
574 // Try again skipping the first sample, if we have sufficient points.
575 if (values.size() > getCurveMinPoints())
576 {
577 PolynomialFit fit2(2, positions.mid(1), values.mid(1));
578 foundFit = fit2.findMinimum(position, 0, params.maxPositionAllowed, &minPos, &minVal);
579 minPos = minPos + 1;
580 }
581 }
582 }
583 else // Hyperbola or Parabola so use the LM solver
584 {
585 auto goal = getGoal(numSteps);
586 params.curveFitting->fitCurve(goal, positions, values, weights, pass1Outliers, params.curveFit, params.useWeights,
587 params.optimisationDirection);
588
589 foundFit = params.curveFitting->findMinMax(position, 0, params.maxPositionAllowed, &minPos, &minVal,
590 static_cast<CurveFitting::CurveFit>(params.curveFit),
591 params.optimisationDirection);
592 }
593
594 if (foundFit)
595 {
596 const int distanceToMin = static_cast<int>(position - minPos);
597 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: fit(%1): %2 = %3 @ %4 distToMin %5")
598 .arg(positions.size()).arg(minPos).arg(minVal).arg(position).arg(distanceToMin);
599 if (distanceToMin >= 0)
600 {
601 // The minimum is further inward.
602 numPolySolutionsFound = 0;
603 numRestartSolutionsFound = 0;
604 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Solutions reset %1 = %2").arg(minPos).arg(minVal);
605
606 // The next code block is an optimisation where the user starts focus off a long way from prime focus.
607 // The idea is that by detecting this situation we can warp quickly to the solution.
608 // However, when the solver is solving on a small number of points (e.g. 6 or 7) its possible that the
609 // best fit curve produces a solution much further in than what turns out to be the actual focus point.
610 // This code therefore results in inappropriate warping. So this code is now disabled.
611 /*if (value / minVal > kDecentValue)
612 {
613 // Only skip samples if the value is a long way from the minimum.
614 const int stepsToMin = distanceToMin / stepSize;
615 // Temporarily increase the step size if the minimum is very far inward.
616 if (stepsToMin >= 8)
617 thisStepSize = stepSize * 4;
618 else if (stepsToMin >= 4)
619 thisStepSize = stepSize * 2;
620 }*/
621 }
622 else if (!bestSamplesHeuristic())
623 {
624 // We have potentially passed the bottom of the curve,
625 // but it's possible it is further back than the start of our sweep.
626 if (minPos > passStartPosition)
627 {
628 numRestartSolutionsFound++;
629 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: RESTART Solution #%1 %2 = %3 @ %4")
630 .arg(numRestartSolutionsFound).arg(minPos).arg(minVal).arg(position);
631 }
632 else
633 {
634 numPolySolutionsFound++;
635 numRestartSolutionsFound = 0;
636 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Solution #%1: %2 = %3 @ %4")
637 .arg(numPolySolutionsFound).arg(minPos).arg(minVal).arg(position);
638 }
639 }
640
641 // The LINEAR algo goes >2 steps past the minimum. The LINEAR 1 PASS algo uses the InitialOutwardSteps parameter
642 // in order to have some user control over the number of steps past the minimum when fitting the curve.
643 // With LINEAR 1 PASS setupSecondPass with a margin of 0.0 as no need to move the focuser further
644 // This will result in a step out, past the solution of either "backlash" id set, or 5 * step size, followed
645 // by a step in to the solution. By stepping out and in, backlash will be neutralised.
646 int howFarToGo;
647 if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
648 howFarToGo = kNumPolySolutionsRequired;
649 else
650 howFarToGo = std::max(1, static_cast<int> (params.initialOutwardSteps) - 1);
651
652 if (numPolySolutionsFound >= howFarToGo)
653 {
654 // We found a minimum. Setup the 2nd pass. We could use either the polynomial min or the
655 // min measured star as the target HFR. I've seen using the polynomial minimum to be
656 // sometimes too conservative, sometimes too low. For now using the min sample.
657 double minMeasurement = *std::min_element(values.begin(), values.end());
658 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: 1stPass solution @ %1: pos %2 val %3, min measurement %4")
659 .arg(position).arg(minPos).arg(minVal).arg(minMeasurement);
660
661 // For Linear setup the second pass with a margin of 2
662 if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
663 return setupSecondPass(static_cast<int>(minPos), minMeasurement, 2.0);
664 else
665 // For Linear 1 Pass do a round on the double minPos parameter before casting to an int
666 // as the cast will truncate and potentially change the position down by 1
667 // The code to display v-curve rounds so this keeps things consistent.
668 // Could make this change for Linear as well but this will then need testfocus to change
669 return finishFirstPass(static_cast<int>(round(minPos)), minMeasurement);
670 }
671 else if (numRestartSolutionsFound >= kNumRestartSolutionsRequired)
672 {
673 // We need to restart--that is the error is rising and thus our initial position
674 // was already past the minimum. Restart the algorithm with a greater initial position.
675 // If the min position from the polynomial solution is not too far from the current start position
676 // use that, but don't let it go too far away.
677 const double highestStartPosition = params.startPosition + params.initialOutwardSteps * params.initialStepSize;
678 params.startPosition = std::min(minPos, highestStartPosition);
679 computeInitialPosition();
680 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Restart. Current pos %1, min pos %2, min val %3, re-starting at %4")
681 .arg(position).arg(minPos).arg(minVal).arg(requestedPosition);
682 return requestedPosition;
683 }
684
685 }
686 else
687 {
688 // Minimum failed indicating the 2nd-order polynomial is an inverted U--it has a maximum,
689 // but no minimum. This is, of course, not a sensible solution for the focuser, but can
690 // happen with noisy data and perhaps small step sizes. We still might be able to take advantage,
691 // and notice whether the polynomial is increasing or decreasing locally. For now, do nothing.
692 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: ******** No poly min: Poly must be inverted");
693 }
694 }
695 else
696 {
697 // Don't have enough samples to reliably fit a polynomial.
698 // Simply step the focus in one more time and iterate.
699 }
700 }
701 else if (gettingWorse())
702 {
703 // Doesn't look like we'll find something close to the min. Retry the 2nd pass.
704 // NOTE: this branch will only be executed for the 2 pass version of Linear
705 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: getting worse, re-running 2nd pass");
706 return setupSecondPass(firstPassBestPosition, firstPassBestValue);
707 }
708
709 return completeIteration(thisStepSize, foundFit, minPos, minVal);
710}
711
712// Get the curve fitting goal, based on the "walk".
713// Start with STANDARD but at the end of the pass use BEST
714CurveFitting::FittingGoal LinearFocusAlgorithm::getGoal(int numSteps)
715{
716 if (params.focusWalk == Focus::FOCUS_WALK_CLASSIC)
717 {
718 if (numSteps > 2.0 * params.initialOutwardSteps)
719 return CurveFitting::FittingGoal::BEST;
720 else
721 return CurveFitting::FittingGoal::STANDARD;
722 }
723 else // FOCUS_WALK_FIXED_STEPS or FOCUS_WALK_CFZ_SHUFFLE
724 {
725 if (numSteps >= params.numSteps)
726 return CurveFitting::FittingGoal::BEST;
727 else
728 return CurveFitting::FittingGoal::STANDARD;
729 }
730}
731
732// Get the minimum number of points to start fitting a curve
733// The default is 5. Try to get past the minimum so the curve shape is hyperbolic/parabolic
734int LinearFocusAlgorithm::getCurveMinPoints()
735{
736 if (params.focusAlgorithm != Focus::FOCUS_LINEAR1PASS)
737 return 5;
738
739 if (params.focusWalk == Focus::FOCUS_WALK_CLASSIC)
740 return params.initialOutwardSteps + 2;
741
742 // For donut busting don't try drawing an intermediate curve... as we need to scrub data first
743 if (params.donutBuster)
744 return params.numSteps;
745
746 return params.numSteps / 2 + 2;
747}
748
749// Process next step for LINEAR1PASS for walks: FOCUS_WALK_FIXED_STEPS and FOCUS_WALK_CFZ_SHUFFLE
750int LinearFocusAlgorithm::linearWalk(int position, double value, const double starWeight)
751{
752 if (!inFirstPass)
753 {
754 return setupSolution(position, value, starWeight);
755 }
756
757 bool foundFit = false;
758 double minPos = -1.0, minVal = 0;
759 ++numSteps;
760 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear1Pass: step %1, linearWalk %2, %3")
761 .arg(numSteps).arg(position).arg(value);
762
763 // If we are within stepsize ticks of where we should be then fine, otherwise try again
764 if (abs(position - requestedPosition) > params.initialStepSize)
765 {
766 qCDebug(KSTARS_EKOS_FOCUS) << QString("linearWalk error: position %1, requested position %2").arg(position)
767 .arg(requestedPosition);
768 return requestedPosition;
769 }
770
771 // Store the sample values.
772 positions.push_back(position);
773 values.push_back(value);
774 weights.push_back(starWeight);
775 pass1Positions.push_back(position);
776 pass1Values.push_back(value);
777 pass1Weights.push_back(starWeight);
778 pass1Outliers.push_back(false);
779
780 if (values.size() < getCurveMinPoints())
781 {
782 // Don't have enough samples to reliably fit a curve.
783 // Simply step the focus in one more time and iterate.
784 }
785 else
786 {
787 if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
788 qCDebug(KSTARS_EKOS_FOCUS) << QString("linearWalk called incorrectly for FOCUS_QUADRATIC");
789 else
790 {
791 // Hyperbola or Parabola so use the LM solver
792 auto goal = getGoal(numSteps);
793 // Check for Donuts
794 if (params.donutBuster)
795 removeDonuts();
796 params.curveFitting->fitCurve(goal, positions, values, weights, pass1Outliers, params.curveFit, params.useWeights,
797 params.optimisationDirection);
798
799 foundFit = params.curveFitting->findMinMax(position, 0, params.maxPositionAllowed, &minPos, &minVal,
800 static_cast<CurveFitting::CurveFit>(params.curveFit), params.optimisationDirection);
801
802 if (numSteps >= params.numSteps)
803 {
804 // linearWalk has now completed either successfully (foundFit=true) or not
805 if (foundFit)
806 return finishFirstPass(static_cast<int>(round(minPos)), minVal);
807 else
808 {
809 done = true;
810 doneString = i18n("Failed to fit curve to data.");
811 failCode = Ekos::FOCUS_FAIL_CURVEFIT;
812 qCDebug(KSTARS_EKOS_FOCUS) << QString("LinearWalk: %1").arg(doneString);
813 debugLog();
814 return -1;
815 }
816 }
817 }
818 }
819
820 int nextStepSize = getNextStepSize();
821 return completeIteration(nextStepSize, foundFit, minPos, minVal);
822}
823
824// Check for donuts
825// Assumption is that we are starting near to focus so that we should have a reasonable v-curve
826// near the central datapoints but the wings may contain donuts=poor quality data.
827// 1. Make sure we have enough data to work with
828// 2. Find the central minimum datapoint
829// 3. Move outwards on each side of the minimum checking the next datapoint. If value drops this could be bad data
830// Only remove the wings - bad data near the centre won't be filtered out.
831void LinearFocusAlgorithm::removeDonuts()
832{
833 // 1. We need at least 7 datapoints to do any analysis
834 if (values.size() <= 7)
835 return;
836
837 QVector<bool> queryOutliers = pass1Outliers;
838
839 // 2. Firstly find the central minimum datapoint
840 int centre = values.size() / 2;
841 int minElement = centre;
842 double minValue = values[centre];
843 for (int i = 1; i < 4; i++)
844 {
845 if (values[centre + i] < minValue)
846 {
847 minValue = values[centre + i];
848 minElement = centre + i;
849 }
850 if (values[centre - i] < minValue)
851 {
852 minValue = values[centre - i];
853 minElement = centre - i;
854 }
855 }
856
857 // 3a. Move outward from the central minimum checking for increasing values
858 double threshold = values[minElement];
859 double nextValue, deltaValue;
860 double tolerance = 0.5; //
861 for (int i = minElement; i > 0; i--)
862 {
863 nextValue = values[i - 1];
864 if (nextValue < threshold)
865 // Smaller value indicates bad data
866 queryOutliers[i - 1] = true;
867 else
868 {
869 deltaValue = (nextValue - values[i]) * tolerance;
870 threshold = nextValue + deltaValue;
871 }
872 }
873
874 // Starting at the furthest outward point, move inward to minimum point, marking bad data (donuts) as outliers
875 for (int i = 0; i < minElement - 3; i++)
876 {
877 if (queryOutliers[i])
878 pass1Outliers[i] = true;
879 else
880 // Stop the process at the first good datapoint - don't discard all bad data
881 break;
882 }
883
884 // 3b. Move inward from the central minimum checking for increasing values
885 threshold = values[minElement];
886 for (int i = minElement; i < values.size() - 1; i++)
887 {
888 nextValue = values[i + 1];
889 if (nextValue < threshold)
890 // Smaller value indicates bad data
891 queryOutliers[i + 1] = true;
892 else
893 {
894 deltaValue = (nextValue - values[i]) * tolerance;
895 threshold = nextValue + deltaValue;
896 }
897 }
898
899 // Starting at the furthest outward point, move inward to minimum point, marking bad data (donuts) as outliers
900 for (int i = values.size() - 1; i > minElement + 3; i--)
901 {
902 if (queryOutliers[i])
903 pass1Outliers[i] = true;
904 else
905 // Stop the process at the first good datapoint - don't discard all bad data
906 break;
907 }
908
909 QString str;
910 for (int i = pass1Outliers.size() - 1; i >= 0; i--)
911 str.append((pass1Outliers[i]) ? 'Y' : 'N');
912 qCDebug(KSTARS_EKOS_FOCUS) << QString("Donut analysis: %1").arg(str);
913
914}
915
916// Function to calculate the next step size for LINEAR1PASS for walks: FOCUS_WALK_FIXED_STEPS and FOCUS_WALK_CFZ_SHUFFLE
917int LinearFocusAlgorithm::getNextStepSize()
918{
919 int nextStepSize, lower, upper;
920
921 switch (params.focusWalk)
922 {
923 case Focus::FOCUS_WALK_CLASSIC:
924 nextStepSize = stepSize;
925 break;
926 case Focus::FOCUS_WALK_FIXED_STEPS:
927 nextStepSize = stepSize;
928 break;
929 case Focus::FOCUS_WALK_CFZ_SHUFFLE:
930 if (params.numSteps % 2)
931 {
932 // Odd number of steps
933 lower = (params.numSteps - 3) / 2;
934 upper = params.numSteps - lower;
935 }
936 else
937 {
938 // Evem number of steps
939 lower = (params.numSteps - 4) / 2;
940 upper = (params.numSteps - lower);
941 }
942
943 if (numSteps <= lower)
944 nextStepSize = stepSize;
945 else if (numSteps >= upper)
946 nextStepSize = stepSize;
947 else
948 nextStepSize = stepSize / 2;
949 break;
950 default:
951 nextStepSize = stepSize;
952 break;
953 }
954
955 return nextStepSize;
956}
957
958int LinearFocusAlgorithm::setupSolution(int position, double value, double weight)
959{
960 focusSolution = position;
961 focusValue = value;
962 focusWeight = weight;
963 done = true;
964 doneString = i18n("Solution found.");
965 failCode = Ekos::FOCUS_FAIL_NONE;
966 if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
967 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: solution @ %1 = %2 (best %3)")
968 .arg(position).arg(value).arg(firstPassBestValue);
969 else // Linear 1 Pass
970 {
971 double delta = fabs(value - firstPassBestValue);
972 QString str("Linear: solution @ ");
973 str.append(QString("%1 = %2 (expected %3) delta=%4").arg(position).arg(value).arg(firstPassBestValue).arg(delta));
974
975 if (params.useWeights && weight > 0.0)
976 {
977 // TODO fix this for weights
978 double numSigmas = delta * weight;
979 str.append(QString(" or %1 sigma %2 than expected")
980 .arg(numSigmas).arg(value <= firstPassBestValue ? "better" : "worse"));
981 if (value <= firstPassBestValue || numSigmas < 1)
982 str.append(QString(". GREAT result"));
983 else if (numSigmas < 3)
984 str.append(QString(". OK result"));
985 else
986 str.append(QString(". POOR result. If this happens repeatedly it may be a sign of poor backlash compensation."));
987 }
988
989 qCInfo(KSTARS_EKOS_FOCUS) << str;
990 }
991 debugLog();
992 return -1;
993}
994
995int LinearFocusAlgorithm::completeIteration(int step, bool foundFit, double minPos, double minVal)
996{
997 if (params.focusAlgorithm == Focus::FOCUS_LINEAR && numSteps == params.maxIterations - 2)
998 {
999 // If we're close to exceeding the iteration limit, retry this pass near the old minimum position.
1000 // For Linear 1 Pass we are still in the first pass so no point retrying the 2nd pass. Just carry on and
1001 // fail in 2 more datapoints if necessary.
1002 const int minIndex = static_cast<int>(std::min_element(values.begin(), values.end()) - values.begin());
1003 return setupSecondPass(positions[minIndex], values[minIndex], 0.5);
1004 }
1005 else if (numSteps > params.maxIterations)
1006 {
1007 // Fail. Exceeded our alloted number of iterations.
1008 done = true;
1009 doneString = i18n("Too many steps.");
1010 failCode = Ekos::FOCUS_FAIL_MAX_ITERS;
1011 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
1012 debugLog();
1013 return -1;
1014 }
1015
1016 // Setup the next sample.
1017 requestedPosition = requestedPosition - step;
1018
1019 // Make sure the next sample is within bounds.
1020 if (requestedPosition < minPositionLimit)
1021 {
1022 // The position is too low. Pick the min value and go to (or retry) a 2nd iteration.
1023 if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
1024 {
1025 const int minIndex = static_cast<int>(std::min_element(values.begin(), values.end()) - values.begin());
1026 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: reached end without Vmin. Restarting %1 pos %2 value %3")
1027 .arg(minIndex).arg(positions[minIndex]).arg(values[minIndex]);
1028 return setupSecondPass(positions[minIndex], values[minIndex]);
1029 }
1030 else // Linear 1 Pass
1031 {
1032 if (foundFit && minPos >= minPositionLimit)
1033 // Although we ran out of room, we have got a viable, although not perfect, solution
1034 return finishFirstPass(static_cast<int>(round(minPos)), minVal);
1035 else
1036 {
1037 // We can't travel far enough to find a solution so fail as focuser parameters require user attention
1038 done = true;
1039 doneString = i18n("Solution lies outside max travel.");
1040 failCode = Ekos::FOCUS_FAIL_FOCUSER_OOB;
1041 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
1042 debugLog();
1043 return -1;
1044 }
1045 }
1046 }
1047 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: requesting position %1").arg(requestedPosition);
1048 return requestedPosition;
1049}
1050
1051// This is called in the 2nd pass, when we've found a sample point that's within tolerance
1052// of the 1st pass' best value. Since it's within tolerance, the 2nd pass might finish, using
1053// the current value as the final solution. However, this method checks to see if we might go
1054// a bit further. Once solutionPending is true, it's committed to finishing soon.
1055// It goes further if:
1056// - the current HFR value is an improvement on the previous 2nd-pass value,
1057// - the number of steps taken so far is not close the max number of allowable steps.
1058// If instead the HFR got worse, we retry the current position a few times to see if it might
1059// get better before giving up.
1060bool LinearFocusAlgorithm::setupPendingSolution(int position)
1061{
1062 const int length = values.size();
1063 const int secondPassIndex = length - secondPassStartIndex;
1064 const double thisValue = values[length - 1];
1065
1066 // ReferenceValue is the value used in the notGettingWorse computation.
1067 // Basically, it's the previous value, or the value before retries.
1068 double referenceValue = (secondPassIndex <= 1) ? 1e6 : values[length - 2];
1069 if (retryNumber > 0 && length - (2 + retryNumber) >= 0)
1070 referenceValue = values[length - (2 + retryNumber)];
1071
1072 // NotGettingWorse: not worse than the previous (non-repeat) 2nd pass sample, or the 1st 2nd pass sample.
1073 const bool notGettingWorse = (secondPassIndex <= 1) || (thisValue <= referenceValue);
1074
1075 // CouldGoFurther: Not near a boundary in position or number of steps.
1076 const bool couldGoFather = (numSteps < params.maxIterations - 2) && (position - stepSize > minPositionLimit);
1077
1078 // Allow passing the 1st pass' minimum HFR position, but take smaller steps and don't retry as much.
1079 int maxNumRetries = 3;
1080 if (position - stepSize / 2 < firstPassBestPosition)
1081 {
1082 stepSize = std::max(2, params.initialStepSize / 4);
1083 maxNumRetries = 1;
1084 }
1085
1086 if (notGettingWorse && couldGoFather)
1087 {
1088 qCDebug(KSTARS_EKOS_FOCUS)
1089 << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, solution pending, searching further")
1090 .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
1091 solutionPending = true;
1092 solutionPendingPosition = position;
1093 solutionPendingValue = thisValue;
1094 retryNumber = 0;
1095 return true;
1096 }
1097 else if (solutionPending && couldGoFather && retryNumber < maxNumRetries &&
1098 (secondPassIndex > 1) && (thisValue >= referenceValue))
1099 {
1100 qCDebug(KSTARS_EKOS_FOCUS)
1101 << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, solution pending, got worse, retrying")
1102 .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
1103 // Try this poisition again.
1104 retryNumber++;
1105 return true;
1106 }
1107 qCDebug(KSTARS_EKOS_FOCUS)
1108 << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, finishing, can't go further")
1109 .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
1110 retryNumber = 0;
1111 return false;
1112}
1113
1114void LinearFocusAlgorithm::debugLog()
1115{
1116 QString str("Linear: points=[");
1117 for (int i = 0; i < positions.size(); ++i)
1118 {
1119 str.append(QString("(%1, %2, %3)").arg(positions[i]).arg(values[i]).arg(weights[i]));
1120 if (i < positions.size() - 1)
1121 str.append(", ");
1122 }
1123 str.append(QString("];iterations=%1").arg(numSteps));
1124 str.append(QString(";duration=%1").arg(stopWatch.elapsed() / 1000));
1125 str.append(QString(";solution=%1").arg(focusSolution));
1126 str.append(QString(";value=%1").arg(focusValue));
1127 str.append(QString(";filter='%1'").arg(params.filterName));
1128 str.append(QString(";temperature=%1").arg(params.temperature));
1129 str.append(QString(";focusalgorithm=%1").arg(params.focusAlgorithm));
1130 str.append(QString(";backlash=%1").arg(params.backlash));
1131 str.append(QString(";curvefit=%1").arg(params.curveFit));
1132 str.append(QString(";useweights=%1").arg(params.useWeights));
1133
1134 qCDebug(KSTARS_EKOS_FOCUS) << str;
1135}
1136
1137// Called by Linear to set state for the second pass
1138int LinearFocusAlgorithm::setupSecondPass(int position, double value, double margin)
1139{
1140 firstPassBestPosition = position;
1141 firstPassBestValue = value;
1142 inFirstPass = false;
1143 solutionPending = false;
1144 secondPassStartIndex = values.size();
1145
1146 // Arbitrarily go back "margin" steps above the best position.
1147 // Could be a problem if backlash were worse than that many steps.
1148 requestedPosition = std::min(static_cast<int>(firstPassBestPosition + stepSize * margin), maxPositionLimit);
1149 stepSize = params.initialStepSize / 2;
1150 qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: 2ndPass starting at %1 step %2").arg(requestedPosition).arg(stepSize);
1151 return requestedPosition;
1152}
1153
1154// Called by L1P to finish off the first pass by setting state variables
1155int LinearFocusAlgorithm::finishFirstPass(int position, double value)
1156{
1157 firstPassBestPosition = position;
1158 firstPassBestValue = value;
1159 inFirstPass = false;
1160 requestedPosition = position;
1161
1162 if (params.refineCurveFit)
1163 {
1164 // We've completed pass1 so, if required, reprocess all the datapoints excluding outliers
1165 // First, we need the distance of each datapoint from the curve
1166 std::vector<std::pair<int, double>> curveDeltas;
1167
1168 params.curveFitting->calculateCurveDeltas(params.curveFit, curveDeltas);
1169
1170 // Check for outliers if there are enough datapoints
1171 if (curveDeltas.size() > 6)
1172 {
1173 // Build a vector of the square of the curve deltas
1174 std::vector<double> curveDeltasX2Vec;
1175 for (ulong i = 0; i < curveDeltas.size(); i++)
1176 curveDeltasX2Vec.push_back(pow(curveDeltas[i].second, 2.0));
1177
1178 auto curveDeltasX2Mean = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_MEAN,
1179 curveDeltasX2Vec);
1180
1181 // Order the curveDeltas, highest first, then check against the limit
1182 // Remove points over the limit, but don't remove too many points to compromise the curve.
1183 double outlierRejection = params.donutBuster ? params.outlierRejection : 0.2;
1184 int maxOutliers = curveDeltas.size() * outlierRejection;
1185
1186 // Parabola has 3 degrees of freedom; Hyperbola and Gaussian 4
1187 double modelUnknowns = params.curveFit == CurveFitting::FOCUS_PARABOLA ? 3.0 : 4.0;
1188
1189 // Use Peirce's Criterion to get the outlier threshold
1190 // Note this operates on the square of the curve deltas
1191 double pc = peirce_criterion(static_cast<double>(curveDeltas.size()), static_cast<double>(maxOutliers), modelUnknowns);
1192 double pc_threshold = sqrt(pc * curveDeltasX2Mean);
1193
1194 // Sort the curve deltas, largest first
1195 std::sort(curveDeltas.begin(), curveDeltas.end(),
1196 [](const std::pair<int, double> &a, const std::pair<int, double> &b)
1197 {
1198 return ( a.second > b.second );
1199 });
1200
1201 // Save off pass1 variables in case they need to be reset later
1202 auto bestPass1Outliers = pass1Outliers;
1203
1204 // Work out how many outliers to exclude
1205 int outliers = 0;
1206 for (ulong i = 0; i < curveDeltas.size(); i++)
1207 {
1208 if(curveDeltas[i].second <= pc_threshold)
1209 break;
1210 else
1211 {
1212 pass1Outliers[curveDeltas[i].first] = true;
1213 outliers++;
1214 }
1215 if (outliers >= maxOutliers)
1216 break;
1217 }
1218
1219 // If there are any outliers then try to improve the curve fit by removing these
1220 // datapoints and reruning the curve fitting process. Note that this may not improve
1221 // the fit so we need to be prepared to reinstate the original solution.
1222 QVector<double> currentCoefficients;
1223 if (outliers > 0 && params.curveFitting->getCurveParams(params.curveFit, currentCoefficients))
1224 {
1225 double currentR2 = params.curveFitting->calculateR2(static_cast<CurveFitting::CurveFit>(params.curveFit));
1226
1227 // Try another curve fit on the data without outliers
1228 params.curveFitting->fitCurve(CurveFitting::FittingGoal::BEST, pass1Positions, pass1Values, pass1Weights,
1229 pass1Outliers, params.curveFit, params.useWeights, params.optimisationDirection);
1230 double minPos, minVal;
1231 bool foundFit = params.curveFitting->findMinMax(position, 0, params.maxPositionAllowed, &minPos, &minVal,
1232 static_cast<CurveFitting::CurveFit>(params.curveFit),
1233 params.optimisationDirection);
1234 if (foundFit)
1235 {
1236 double newR2 = params.curveFitting->calculateR2(static_cast<CurveFitting::CurveFit>(params.curveFit));
1237 if (newR2 > currentR2)
1238 {
1239 // We have a better solution so we'll use it
1240 requestedPosition = minPos;
1241 qCDebug(KSTARS_EKOS_FOCUS) << QString("Refine Curve Fit improved focus from %1 to %2. R2 improved from %3 to %4")
1242 .arg(position).arg(minPos).arg(currentR2).arg(newR2);
1243 }
1244 else
1245 {
1246 // New solution is not better than previous one so go with the previous one
1247 qCDebug(KSTARS_EKOS_FOCUS) <<
1248 QString("Refine Curve Fit could not improve on the original solution");
1249 pass1Outliers = bestPass1Outliers;
1250 params.curveFitting->setCurveParams(params.curveFit, currentCoefficients);
1251 }
1252 }
1253 else
1254 {
1255 qCDebug(KSTARS_EKOS_FOCUS) <<
1256 QString("Refine Curve Fit failed to fit curve whilst refining curve fit. Running with original solution");
1257 pass1Outliers = bestPass1Outliers;
1258 params.curveFitting->setCurveParams(params.curveFit, currentCoefficients);
1259 }
1260 }
1261 }
1262 }
1263 return requestedPosition;
1264}
1265
1266// Peirce's criterion based on Gould's methodology for outlier identification
1267// See https://en.wikipedia.org/wiki/Peirce%27s_criterion for details incl Peirce's original paper
1268// and Gould's paper both referenced in the notes. The wikipedia entry also includes some code fragments
1269// that form the basis of this function.
1270//
1271// Arguments:
1272// N = total number of observations
1273// n = number of outliers to be removed
1274// m = number of model unknowns
1275// Returns: squared error threshold (x2)
1276double LinearFocusAlgorithm::peirce_criterion(double N, double n, double m)
1277{
1278 double x2 = 0.0;
1279 if (N > 1)
1280 {
1281 // Calculate Q (Nth root of Gould's equation B):
1282 double Q = (pow(n, (n / N)) * pow((N - n), ((N - n) / N))) / N;
1283
1284 // Initialize R values
1285 double r_new = 1.0;
1286 double r_old = 0.0;
1287
1288 // Start iteration to converge on R:
1289 while (abs(r_new - r_old) > (N * 2.0e-16))
1290 {
1291 // Calculate Lamda
1292 // (1 / (N - n) th root of Gould's equation B
1293 double ldiv = pow(r_new, n);
1294 if (ldiv == 0)
1295 ldiv = 1.0e-6;
1296
1297 double Lamda = pow((pow(Q, N) / (ldiv)), (1.0 / (N - n)));
1298
1299 // Calculate x -squared(Gould's equation C)
1300 x2 = 1.0 + (N - m - n) / n * (1.0 - pow(Lamda, 2.0));
1301
1302 if (x2 < 0.0)
1303 {
1304 // If x2 goes negative, return 0
1305 x2 = 0.0;
1306 r_old = r_new;
1307 }
1308 else
1309 {
1310 // Use x-squared to update R(Gould 's equation D)
1311 r_old = r_new;
1312 r_new = exp((x2 - 1) / 2.0) * gsl_sf_erfc(sqrt(x2) / sqrt(2.0));
1313 }
1314 }
1315 }
1316
1317 return x2;
1318}
1319
1320// Return true if one of the 2 recent samples is among the best 2 samples so far.
1321bool LinearFocusAlgorithm::bestSamplesHeuristic()
1322{
1323 const int length = values.size();
1324 if (length < 5) return true;
1325 QVector<double> tempValues = values;
1326 std::nth_element(tempValues.begin(), tempValues.begin() + 2, tempValues.end());
1327 double secondBest = tempValues[1];
1328 if ((values[length - 1] <= secondBest) || (values[length - 2] <= secondBest))
1329 return true;
1330 return false;
1331}
1332
1333// Return true if there are "streak" consecutive values which are successively worse.
1334bool LinearFocusAlgorithm::gettingWorse()
1335{
1336 // Must have this many consecutive values getting worse.
1337 constexpr int streak = 3;
1338 const int length = values.size();
1339 if (secondPassStartIndex < 0)
1340 return false;
1341 if (length < streak + 1)
1342 return false;
1343 // This insures that all the values we're checking are in the latest 2nd pass.
1344 if (length - secondPassStartIndex < streak + 1)
1345 return false;
1346 for (int i = length - 1; i >= length - streak; --i)
1347 if (values[i] <= values[i - 1])
1348 return false;
1349 return true;
1350}
1351
1352}
1353
QString i18n(const char *text, const TYPE &arg...)
Ekos is an advanced Astrophotography tool for Linux.
Definition align.cpp:83
iterator begin()
iterator end()
qsizetype size() const const
QString & append(QChar ch)
QString arg(Args &&... args) const const
void push_back(QChar ch)
qsizetype size() 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.