Kstars

focusalgorithms.cpp
1 /*
2  SPDX-FileCopyrightText: 2019 Hy Murveit <[email protected]>
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 
17 namespace 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  */
26 class 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,
49  const QList<Edge*> *stars) override;
50 
51  FocusAlgorithmInterface *Copy() override;
52 
53  void getMeasurements(QVector<int> *pos, QVector<double> *hfrs, QVector<double> *sds) const override
54  {
55  *pos = positions;
56  *hfrs = values;
57  *sds = sigmas;
58  }
59 
60  void getPass1Measurements(QVector<int> *pos, QVector<double> *hfrs, QVector<double> *sds) const override
61  {
62  *pos = pass1Positions;
63  *hfrs = pass1Values;
64  *sds = pass1Sigmas;
65  }
66 
67  double latestHFR() const override
68  {
69  if (values.size() > 0)
70  return values.last();
71  return -1;
72  }
73 
74  QString getTextStatus(double R2 = 0) const override;
75 
76  private:
77 
78  // Called in newMeasurement. Sets up the next iteration.
79  int completeIteration(int step, bool foundFit, double minPos, double minVal);
80 
81  // Does the bookkeeping for the final focus solution.
82  int setupSolution(int position, double value, double sigma);
83 
84  // Called when we've found a solution, e.g. the HFR value is within tolerance of the desired value.
85  // It it returns true, then it's decided tht we should try one more sample for a possible improvement.
86  // If it returns false, then "play it safe" and use the current position as the solution.
87  bool setupPendingSolution(int position);
88 
89  // Determines the desired focus position for the first sample.
90  void computeInitialPosition();
91 
92  // Sets the internal state for re-finding the minimum, and returns the requested next
93  // position to sample.
94  // Called by Linear. L1P calls finishFirstPass instead
95  int setupSecondPass(int position, double value, double margin = 2.0);
96 
97  // Called by L1P to finish off the first pass, setting state variables.
98  int finishFirstPass(int position, double value);
99 
100  // Used in the 2nd pass. Focus is getting worse. Requires several consecutive samples getting worse.
101  bool gettingWorse();
102 
103  // If one of the last 2 samples are as good or better than the 2nd best so far, return true.
104  bool bestSamplesHeuristic();
105 
106  // Adds to the debug log a line summarizing the result of running this algorithm.
107  void debugLog();
108 
109  // Use the detailed star measurements to adjust the HFR value.
110  // See comment above method definition.
111  double relativeHFR(double origHFR, const QList<Edge*> *stars);
112 
113  // Calculate the standard deviation (=sigma) of the star HFRs
114  double calculateStarSigma(const bool useWeights, const QList<Edge*> *stars);
115 
116  // Used to time the focus algorithm.
117  QTime stopWatch;
118 
119  // A vector containing the HFR values sampled by this algorithm so far.
121  QVector<double> pass1Values;
122  // A vector containing star measurement standard deviation = sigma
123  QVector<double> sigmas;
124  QVector<double> pass1Sigmas;
125  // A vector containing the focus positions corresponding to the HFR values stored above.
126  QVector<int> positions;
127  QVector<int> pass1Positions;
128 
129  // Focus position requested by this algorithm the previous step.
130  int requestedPosition;
131  // The position where the first pass is begun. Usually requestedPosition unless there's a restart.
132  int passStartPosition;
133  // Number of iterations processed so far. Used to see it doesn't exceed params.maxIterations.
134  int numSteps;
135  // The best value in the first pass. The 2nd pass attempts to get within
136  // tolerance of this value.
137  double firstPassBestValue;
138  // The position of the minimum found in the first pass.
139  int firstPassBestPosition;
140  // The sampling interval--the recommended number of focuser steps moved inward each iteration
141  // of the first pass.
142  int stepSize;
143  // The minimum focus position to use. Computed from the focuser limits and maxTravel.
144  int minPositionLimit;
145  // The maximum focus position to use. Computed from the focuser limits and maxTravel.
146  int maxPositionLimit;
147  // Counter for the number of times a v-curve minimum above the current position was found,
148  // which implies the initial focus sweep has passed the minimum and should be terminated.
149  int numPolySolutionsFound;
150  // Counter for the number of times a v-curve minimum above the passStartPosition was found,
151  // which implies the sweep should restart at a higher position.
152  int numRestartSolutionsFound;
153  // The index (into values and positions) when the most recent 2nd pass started.
154  int secondPassStartIndex;
155  // True if performing the first focus sweep.
156  bool inFirstPass;
157  // True if the 2nd pass has found a solution, and it's now optimizing the solution.
158  bool solutionPending;
159  // When we're near done, but the HFR just got worse, we may retry the current position
160  // in case the HFR value was noisy.
161  int retryNumber = 0;
162 
163  // Keep the solution-pending position/value for the status messages.
164  double solutionPendingValue = 0;
165  int solutionPendingPosition = 0;
166 
167  // RelativeHFR variables
168  bool relativeHFREnabled = false;
169  QVector<FocusStars> starLists;
170  double bestRelativeHFR = 0;
171  int bestRelativeHFRIndex = 0;
172 };
173 
174 // Copies the object. Used in testing to examine alternate possible inputs given
175 // the current state.
176 FocusAlgorithmInterface *LinearFocusAlgorithm::Copy()
177 {
178  LinearFocusAlgorithm *alg = new LinearFocusAlgorithm(params);
179  *alg = *this;
180  return dynamic_cast<FocusAlgorithmInterface*>(alg);
181 }
182 
183 FocusAlgorithmInterface *MakeLinearFocuser(const FocusAlgorithmInterface::FocusParams &params)
184 {
185  return new LinearFocusAlgorithm(params);
186 }
187 
188 LinearFocusAlgorithm::LinearFocusAlgorithm(const FocusParams &focusParams)
189  : FocusAlgorithmInterface(focusParams)
190 {
191  stopWatch.start();
192  // These variables don't get reset if we restart the algorithm.
193  numSteps = 0;
194  maxPositionLimit = std::min(params.maxPositionAllowed, params.startPosition + params.maxTravel);
195  minPositionLimit = std::max(params.minPositionAllowed, params.startPosition - params.maxTravel);
196  computeInitialPosition();
197 }
198 
199 QString LinearFocusAlgorithm::getTextStatus(double R2) const
200 {
201  if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
202  {
203  if (done)
204  {
205  if (focusSolution > 0)
206  return QString("Solution: %1").arg(focusSolution);
207  else
208  return "Failed";
209  }
210  if (solutionPending)
211  return QString("Pending: %1, %2").arg(solutionPendingPosition).arg(solutionPendingValue, 0, 'f', 2);
212  if (inFirstPass)
213  return "1st Pass";
214  else
215  return QString("2nd Pass. 1st: %1, %2")
216  .arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 2);
217  }
218  else // Linear 1 Pass
219  {
220  QString text = "L1P: ";
221  if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
222  text.append("Quadratic");
223  else if (params.curveFit == CurveFitting::FOCUS_HYPERBOLA)
224  {
225  text.append("Hyperbola");
226  text.append(params.useWeights ? " (W)" : " (U)");
227  }
228  else if (params.curveFit == CurveFitting::FOCUS_PARABOLA)
229  {
230  text.append("Parabola");
231  text.append(params.useWeights ? " (W)" : " (U)");
232  }
233 
234  if (inFirstPass)
235  return text;
236  else if (!done)
237  return text.append(". Moving to Solution");
238  else
239  {
240  if (focusSolution > 0)
241  {
242  if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
243  return text.append(" Solution: %1").arg(focusSolution);
244  else
245  // Add R2 to 2 decimal places. Round down to be conservative
246  return text.append(" Solution: %1, R²=%2").arg(focusSolution).arg(trunc(R2 * 100.0) / 100.0, 0, 'f', 2);
247  }
248  else
249  return text.append(" Failed");
250  }
251  }
252 }
253 
254 void LinearFocusAlgorithm::computeInitialPosition()
255 {
256  // These variables get reset if the algorithm is restarted.
257  stepSize = params.initialStepSize;
258  inFirstPass = true;
259  solutionPending = false;
260  retryNumber = 0;
261  firstPassBestValue = -1;
262  firstPassBestPosition = 0;
263  numPolySolutionsFound = 0;
264  numRestartSolutionsFound = 0;
265  secondPassStartIndex = -1;
266 
267  qCDebug(KSTARS_EKOS_FOCUS)
268  << 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")
269  .arg(params.maxTravel).arg(params.initialStepSize).arg(params.startPosition).arg(params.minPositionAllowed)
270  .arg(params.maxPositionAllowed).arg(params.maxIterations).arg(params.focusTolerance).arg(minPositionLimit)
271  .arg(maxPositionLimit).arg(params.initialOutwardSteps).arg(params.focusAlgorithm).arg(params.backlash)
272  .arg(params.curveFit).arg(params.useWeights);
273 
274  requestedPosition = std::min(maxPositionLimit,
275  static_cast<int>(params.startPosition + params.initialOutwardSteps * params.initialStepSize));
276  passStartPosition = requestedPosition;
277  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: initialPosition %1 sized %2")
278  .arg(requestedPosition).arg(params.initialStepSize);
279 }
280 
281 // The Problem:
282 // The HFR values passed in to these focus algorithms are simply the mean or median HFR values
283 // for all stars detected (with some other constraints). However, due to randomness in the
284 // star detection scheme, two sets of star measurements may use different sets of actual stars to
285 // compute their mean HFRs. This adds noise to determining best focus (e.g. including a
286 // wider star in one set but not the other will tend to increase the HFR for that star set).
287 // relativeHFR tries to correct for this by comparing two sets of stars only using their common stars.
288 //
289 // The Solution:
290 // We maintain a reference set of stars, along with the HFR computed for those reference stars (bestRelativeHFR).
291 // We find the common set of stars between an input star set and those reference stars. We compute
292 // HFRs for the input stars in that common set, and for the reference common set as well.
293 // The ratio of those common-set HFRs multiply bestRelativeHFR to generate the relative HFR for this
294 // input star set--that is, it's the HFR "relative to the reference set". E.g. if the common-set HFR
295 // for the input stars is twice worse than that from the reference common stars, then the relative HFR
296 // would be 2 * bestRelativeHFR.
297 // The reference set is the best HFR set of stars seen so far, but only from the 1st pass.
298 // Thus, in the 2nd pass, the reference stars would be the best HFR set from the 1st pass.
299 //
300 // In unusual cases, e.g. when the common set can't be computed, we just return the input origHFR.
301 double LinearFocusAlgorithm::relativeHFR(double origHFR, const QList<Edge*> *stars)
302 {
303  constexpr double minHFR = 0.3;
304  if (origHFR < minHFR) return origHFR;
305 
306  const int currentIndex = values.size();
307  const bool isFirstSample = (currentIndex == 0);
308  double relativeHFR = origHFR;
309 
310  if (isFirstSample && stars != nullptr)
311  relativeHFREnabled = true;
312 
313  if (relativeHFREnabled && stars == nullptr)
314  {
315  // Error--we have been getting relativeHFR stars, and now we're not.
316  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Error: Inconsistent relativeHFR, disabling.");
317  relativeHFREnabled = false;
318  }
319 
320  if (!relativeHFREnabled)
321  return origHFR;
322 
323  if (starLists.size() != positions.size())
324  {
325  // Error--inconsistent data structures!
326  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Error: Inconsistent relativeHFR data structures, disabling.");
327  relativeHFREnabled = false;
328  return origHFR;
329  }
330 
331  // Copy the input stars.
332  // FocusStars enables us to measure HFR relatively consistently,
333  // by using the same stars when comparing two measurements.
334  constexpr int starDistanceThreshold = 5; // Max distance (pixels) for two positions to be considered the same star.
335  FocusStars fs(*stars, starDistanceThreshold);
336  starLists.push_back(fs);
337  auto &currentStars = starLists.back();
338 
339  if (isFirstSample)
340  {
341  relativeHFR = currentStars.getHFR();
342  if (relativeHFR <= 0)
343  {
344  // Fall back to the simpler scheme.
345  relativeHFREnabled = false;
346  return origHFR;
347  }
348  // 1st measurement. By definition this is the best HFR.
349  bestRelativeHFR = relativeHFR;
350  bestRelativeHFRIndex = currentIndex;
351  }
352  else
353  {
354  // HFR computed relative to the best measured so far.
355  auto &bestStars = starLists[bestRelativeHFRIndex];
356  double hfr = currentStars.relativeHFR(bestStars, bestRelativeHFR);
357  if (hfr > 0)
358  {
359  relativeHFR = hfr;
360  }
361  else
362  {
363  relativeHFR = currentStars.getHFR();
364  if (relativeHFR <= 0)
365  return origHFR;
366  }
367 
368  // In the 1st pass we compute the current HFR relative to the best HFR measured yet.
369  // In the 2nd pass we compute the current HFR relative to the best HFR in the 1st pass.
370  if (inFirstPass && relativeHFR <= bestRelativeHFR)
371  {
372  bestRelativeHFR = relativeHFR;
373  bestRelativeHFRIndex = currentIndex;
374  }
375  }
376 
377  qCDebug(KSTARS_EKOS_FOCUS) << QString("RelativeHFR: orig %1 computed %2 relative %3")
378  .arg(origHFR).arg(currentStars.getHFR()).arg(relativeHFR);
379 
380  return relativeHFR;
381 }
382 
383 // Calculate the standard deviation (=sigma) of the star HFR measurements
384 // on the stars used in getHFR. Sigma is then used in the curve fitting process.
385 double LinearFocusAlgorithm::calculateStarSigma(const bool useWeights, const QList<Edge*> *stars)
386 {
387  if (stars == nullptr || stars->size() <= 0 || !useWeights)
388  // If we can't calculate sigma set to 1 - which will stop curve fit weightings being used
389  return 1.0f;
390 
391  const int n = stars->size();
392  double sum = 0;
393  for (int i = 0; i < n; ++i)
394  sum += stars->value(i)->HFR;
395  const double mean = sum / n;
396  double variance = 0;
397  for(int i = 0; i < n; ++i)
398  variance += pow(stars->value(i)->HFR - mean, 2.0);
399  variance = variance / n;
400  return sqrt(variance);
401 }
402 
403 int LinearFocusAlgorithm::newMeasurement(int position, double value, const QList<Edge*> *stars)
404 {
405  double minPos = -1.0, minVal = 0;
406  bool foundFit = false;
407 
408  double starSigma = 1.0, origValue = value;
409  // For QUADRATIC continue to use the relativeHFR functionality
410  // For HYPERBOLA and PARABOLA the stars used for the HDR calculation and the sigma calculation
411  // should be the same. For now, we will use the full set of stars and therefore not adjust the HFR
412  // JEE TODO: revisit adding relativeHFR functionality for FOCUS_PARABOLA and FOCUS_HYPERBOLA
413  if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
414  value = relativeHFR(value, stars);
415  else
416  // Calculate the standard deviation (=sigma) of the star measurements to be used by the curve fitting process.
417  // Currently only available for Hyperbola and Parabola
418  starSigma = calculateStarSigma(params.useWeights, stars);
419 
420  // For LINEAR 1 PASS don't bother with a full 2nd pass just jump to the solution
421  // Do the step out and back in to deal with backlash
422  if (params.focusAlgorithm == Focus::FOCUS_LINEAR1PASS && !inFirstPass)
423  {
424  return setupSolution(position, value, starSigma);
425  }
426 
427  int thisStepSize = stepSize;
428  ++numSteps;
429  if (inFirstPass)
430  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: step %1, newMeasurement(%2, %3 -> %4, %5)")
431  .arg(numSteps).arg(position).arg(origValue).arg(value)
432  .arg(stars == nullptr ? 0 : stars->size());
433  else
434  qCDebug(KSTARS_EKOS_FOCUS)
435  << QString("Linear: step %1, newMeasurement(%2, %3 -> %4, %5) 1stBest %6 %7").arg(numSteps)
436  .arg(position).arg(origValue).arg(value).arg(stars == nullptr ? 0 : stars->size())
437  .arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
438 
439  const int LINEAR_POSITION_TOLERANCE = params.initialStepSize;
440  if (abs(position - requestedPosition) > LINEAR_POSITION_TOLERANCE)
441  {
442  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error didn't get the requested position");
443  return requestedPosition;
444  }
445  // Have we already found a solution?
446  if (focusSolution != -1)
447  {
448  doneString = i18n("Called newMeasurement after a solution was found.");
449  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
450  debugLog();
451  return -1;
452  }
453 
454  // Store the sample values.
455  values.push_back(value);
456  positions.push_back(position);
457  sigmas.push_back(starSigma);
458  if (inFirstPass)
459  {
460  pass1Values.push_back(value);
461  pass1Positions.push_back(position);
462  pass1Sigmas.push_back(starSigma);
463  }
464 
465  // If we're in the 2nd pass and either
466  // - the current value is within tolerance, or
467  // - we're optimizing because we've previously found a within-tolerance solution,
468  // then setupPendingSolution decides whether to continue optimizing or to complete.
469  if (solutionPending ||
470  (!inFirstPass && (value < firstPassBestValue * (1.0 + params.focusTolerance))))
471  {
472  if (setupPendingSolution(position))
473  // We can continue to look a little further.
474  return completeIteration(retryNumber > 0 ? 0 : stepSize, false, -1.0f, -1.0f);
475  else
476  // Finish now
477  return setupSolution(position, value, starSigma);
478  }
479 
480  if (inFirstPass)
481  {
482  constexpr int kMinPolynomialPoints = 5;
483  constexpr int kNumPolySolutionsRequired = 2;
484  constexpr int kNumRestartSolutionsRequired = 3;
485  constexpr double kDecentValue = 2.5;
486 
487  if (values.size() >= kMinPolynomialPoints)
488  {
489  if (params.curveFit == CurveFitting::FOCUS_QUADRATIC)
490  {
491  PolynomialFit fit(2, positions, values);
492  foundFit = fit.findMinimum(position, 0, params.maxPositionAllowed, &minPos, &minVal);
493  if (!foundFit)
494  {
495  // I've found that the first sample can be odd--perhaps due to backlash.
496  // Try again skipping the first sample, if we have sufficient points.
497  if (values.size() > kMinPolynomialPoints)
498  {
499  PolynomialFit fit2(2, positions.mid(1), values.mid(1));
500  foundFit = fit2.findMinimum(position, 0, params.maxPositionAllowed, &minPos, &minVal);
501  minPos = minPos + 1;
502  }
503  }
504  }
505  else // Hyperbola or Parabola so use the LM solver
506  {
507  curveFit.fitCurve(positions, values, sigmas, params.curveFit, params.useWeights);
508  foundFit = curveFit.findMin(position, 0, params.maxPositionAllowed, &minPos, &minVal,
509  static_cast<CurveFitting::CurveFit>(params.curveFit));
510  }
511 
512  if (foundFit)
513  {
514  const int distanceToMin = static_cast<int>(position - minPos);
515  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: fit(%1): %2 = %3 @ %4 distToMin %5")
516  .arg(positions.size()).arg(minPos).arg(minVal).arg(position).arg(distanceToMin);
517  if (distanceToMin >= 0)
518  {
519  // The minimum is further inward.
520  numPolySolutionsFound = 0;
521  numRestartSolutionsFound = 0;
522  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Solutions reset %1 = %2").arg(minPos).arg(minVal);
523  if (value > kDecentValue)
524  {
525  // Only skip samples if the HFV values aren't very good.
526  const int stepsToMin = distanceToMin / stepSize;
527  // Temporarily increase the step size if the minimum is very far inward.
528  if (stepsToMin >= 8)
529  thisStepSize = stepSize * 4;
530  else if (stepsToMin >= 4)
531  thisStepSize = stepSize * 2;
532  }
533  }
534  else if (!bestSamplesHeuristic())
535  {
536  // We have potentially passed the bottom of the curve,
537  // but it's possible it is further back than the start of our sweep.
538  if (minPos > passStartPosition)
539  {
540  numRestartSolutionsFound++;
541  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: RESTART Solution #%1 %2 = %3 @ %4")
542  .arg(numRestartSolutionsFound).arg(minPos).arg(minVal).arg(position);
543  }
544  else
545  {
546  numPolySolutionsFound++;
547  numRestartSolutionsFound = 0;
548  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Solution #%1: %2 = %3 @ %4")
549  .arg(numPolySolutionsFound).arg(minPos).arg(minVal).arg(position);
550  }
551  }
552 
553  // The LINEAR algo goes >2 steps past the minimum. The LINEAR 1 PASS algo uses the InitialOutwardSteps parameter
554  // in order to have some user control over the number of steps past the minimum when fitting the curve.
555  // JEE TODO: Are we stepping out too far and wasting time? Needs more analysis.
556  // With LINEAR 1 PASS setupSecondPass with a margin of 0.0 as no need to move the focuser further
557  // This will result in a step out, past the solution of either "backlash" id set, or 5 * step size, followed
558  // by a step in to the solution. By stepping out and in, backlash will be neutralised.
559  int howFarToGo;
560  if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
561  howFarToGo = kNumPolySolutionsRequired;
562  else
563  howFarToGo = std::max(1, static_cast<int> (params.initialOutwardSteps) - 1);
564 
565  if (numPolySolutionsFound >= howFarToGo)
566  {
567  // We found a minimum. Setup the 2nd pass. We could use either the polynomial min or the
568  // min measured star as the target HFR. I've seen using the polynomial minimum to be
569  // sometimes too conservative, sometimes too low. For now using the min sample.
570  double minMeasurement = *std::min_element(values.begin(), values.end());
571  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: 1stPass solution @ %1: pos %2 val %3, min measurement %4")
572  .arg(position).arg(minPos).arg(minVal).arg(minMeasurement);
573 
574  // For Linear setup the second pass with a margin of 2
575  if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
576  return setupSecondPass(static_cast<int>(minPos), minMeasurement, 2.0);
577  else
578  // For Linear 1 Pass do a round on the double minPos parameter before casting to an int
579  // as the cast will truncate and potentially change the position down by 1
580  // The code to display v-curve rounds so this keeps things consistent.
581  // Could make this change for Linear as well but this will then need testfocus to change
582  return finishFirstPass(static_cast<int>(round(minPos)), minMeasurement);
583  }
584  else if (numRestartSolutionsFound >= kNumRestartSolutionsRequired)
585  {
586  // We need to restart--that is the error is rising and thus our initial position
587  // was already past the minimum. Restart the algorithm with a greater initial position.
588  // If the min position from the polynomial solution is not too far from the current start position
589  // use that, but don't let it go too far away.
590  const double highestStartPosition = params.startPosition + params.initialOutwardSteps * params.initialStepSize;
591  params.startPosition = std::min(minPos, highestStartPosition);
592  computeInitialPosition();
593  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: Restart. Current pos %1, min pos %2, min val %3, re-starting at %4")
594  .arg(position).arg(minPos).arg(minVal).arg(requestedPosition);
595  return requestedPosition;
596  }
597 
598  }
599  else
600  {
601  // Minimum failed indicating the 2nd-order polynomial is an inverted U--it has a maximum,
602  // but no minimum. This is, of course, not a sensible solution for the focuser, but can
603  // happen with noisy data and perhaps small step sizes. We still might be able to take advantage,
604  // and notice whether the polynomial is increasing or decreasing locally. For now, do nothing.
605  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: ******** No poly min: Poly must be inverted");
606  }
607  }
608  else
609  {
610  // Don't have enough samples to reliably fit a polynomial.
611  // Simply step the focus in one more time and iterate.
612  }
613  }
614  else if (gettingWorse())
615  {
616  // Doesn't look like we'll find something close to the min. Retry the 2nd pass.
617  // NOTE: this branch will only be executed for the 2 pass version of Linear
618  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: getting worse, re-running 2nd pass");
619  return setupSecondPass(firstPassBestPosition, firstPassBestValue);
620  }
621 
622  return completeIteration(thisStepSize, foundFit, minPos, minVal);
623 }
624 
625 int LinearFocusAlgorithm::setupSolution(int position, double value, double sigma)
626 {
627  focusSolution = position;
628  focusHFR = value;
629  done = true;
630  doneString = i18n("Solution found.");
631  if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
632  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: solution @ %1 = %2 (best %3)")
633  .arg(position).arg(value).arg(firstPassBestValue);
634  else // Linear 1 Pass
635  {
636  double delta = fabs(value - firstPassBestValue);
637  QString str("Linear: solution @ ");
638  str.append(QString("%1 = %2 (expected %3) delta=%4").arg(position).arg(value).arg(firstPassBestValue).arg(delta));
639 
640  if (params.useWeights && sigma > 0.0)
641  {
642  double numSigmas = delta / sigma;
643  str.append(QString(" or %1 sigma %2 than expected")
644  .arg(numSigmas).arg(value <= firstPassBestValue ? "better" : "worse"));
645  if (value <= firstPassBestValue || numSigmas < 1)
646  str.append(QString(". GREAT result"));
647  else if (numSigmas < 3)
648  str.append(QString(". OK result"));
649  else
650  str.append(QString(". POOR result. If this happens repeatedly it may be a sign of poor backlash compensation."));
651  }
652 
653  qCInfo(KSTARS_EKOS_FOCUS) << str;
654  }
655  debugLog();
656  return -1;
657 }
658 
659 int LinearFocusAlgorithm::completeIteration(int step, bool foundFit, double minPos, double minVal)
660 {
661  if (params.focusAlgorithm == Focus::FOCUS_LINEAR && numSteps == params.maxIterations - 2)
662  {
663  // If we're close to exceeding the iteration limit, retry this pass near the old minimum position.
664  // For Linear 1 Pass we are still in the first pass so no point retrying the 2nd pass. Just carry on and
665  // fail in 2 more datapoints if necessary.
666  const int minIndex = static_cast<int>(std::min_element(values.begin(), values.end()) - values.begin());
667  return setupSecondPass(positions[minIndex], values[minIndex], 0.5);
668  }
669  else if (numSteps > params.maxIterations)
670  {
671  // Fail. Exceeded our alloted number of iterations.
672  done = true;
673  doneString = i18n("Too many steps.");
674  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
675  debugLog();
676  return -1;
677  }
678 
679  // Setup the next sample.
680  requestedPosition = requestedPosition - step;
681 
682  // Make sure the next sample is within bounds.
683  if (requestedPosition < minPositionLimit)
684  {
685  // The position is too low. Pick the min value and go to (or retry) a 2nd iteration.
686  if (params.focusAlgorithm == Focus::FOCUS_LINEAR)
687  {
688  const int minIndex = static_cast<int>(std::min_element(values.begin(), values.end()) - values.begin());
689  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: reached end without Vmin. Restarting %1 pos %2 value %3")
690  .arg(minIndex).arg(positions[minIndex]).arg(values[minIndex]);
691  return setupSecondPass(positions[minIndex], values[minIndex]);
692  }
693  else // Linear 1 Pass
694  {
695  if (foundFit && minPos >= minPositionLimit)
696  // Although we ran out of room, we have got a viable, although not perfect, solution
697  return finishFirstPass(static_cast<int>(round(minPos)), minVal);
698  else
699  {
700  // We can't travel far enough to find a solution so fail as focuser parameters require user attention
701  done = true;
702  doneString = i18n("Solution lies outside max travel.");
703  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: error %1").arg(doneString);
704  debugLog();
705  return -1;
706  }
707  }
708  }
709  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: requesting position %1").arg(requestedPosition);
710  return requestedPosition;
711 }
712 
713 // This is called in the 2nd pass, when we've found a sample point that's within tolerance
714 // of the 1st pass' best value. Since it's within tolerance, the 2nd pass might finish, using
715 // the current value as the final solution. However, this method checks to see if we might go
716 // a bit further. Once solutionPending is true, it's committed to finishing soon.
717 // It goes further if:
718 // - the current HFR value is an improvement on the previous 2nd-pass value,
719 // - the number of steps taken so far is not close the max number of allowable steps.
720 // If instead the HFR got worse, we retry the current position a few times to see if it might
721 // get better before giving up.
722 bool LinearFocusAlgorithm::setupPendingSolution(int position)
723 {
724  const int length = values.size();
725  const int secondPassIndex = length - secondPassStartIndex;
726  const double thisValue = values[length - 1];
727 
728  // ReferenceValue is the value used in the notGettingWorse computation.
729  // Basically, it's the previous value, or the value before retries.
730  double referenceValue = (secondPassIndex <= 1) ? 1e6 : values[length - 2];
731  if (retryNumber > 0 && length - (2 + retryNumber) >= 0)
732  referenceValue = values[length - (2 + retryNumber)];
733 
734  // NotGettingWorse: not worse than the previous (non-repeat) 2nd pass sample, or the 1st 2nd pass sample.
735  const bool notGettingWorse = (secondPassIndex <= 1) || (thisValue <= referenceValue);
736 
737  // CouldGoFurther: Not near a boundary in position or number of steps.
738  const bool couldGoFather = (numSteps < params.maxIterations - 2) && (position - stepSize > minPositionLimit);
739 
740  // Allow passing the 1st pass' minimum HFR position, but take smaller steps and don't retry as much.
741  int maxNumRetries = 3;
742  if (position - stepSize / 2 < firstPassBestPosition)
743  {
744  stepSize = std::max(2, params.initialStepSize / 4);
745  maxNumRetries = 1;
746  }
747 
748  if (notGettingWorse && couldGoFather)
749  {
750  qCDebug(KSTARS_EKOS_FOCUS)
751  << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, solution pending, searching further")
752  .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
753  solutionPending = true;
754  solutionPendingPosition = position;
755  solutionPendingValue = thisValue;
756  retryNumber = 0;
757  return true;
758  }
759  else if (solutionPending && couldGoFather && retryNumber < maxNumRetries &&
760  (secondPassIndex > 1) && (thisValue >= referenceValue))
761  {
762  qCDebug(KSTARS_EKOS_FOCUS)
763  << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, solution pending, got worse, retrying")
764  .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
765  // Try this poisition again.
766  retryNumber++;
767  return true;
768  }
769  qCDebug(KSTARS_EKOS_FOCUS)
770  << QString("Linear: %1: Position(%2) & HFR(%3) -- Pass1: %4 %5, finishing, can't go further")
771  .arg(length).arg(position).arg(thisValue, 0, 'f', 3).arg(firstPassBestPosition).arg(firstPassBestValue, 0, 'f', 3);
772  retryNumber = 0;
773  return false;
774 }
775 
776 void LinearFocusAlgorithm::debugLog()
777 {
778  QString str("Linear: points=[");
779  for (int i = 0; i < positions.size(); ++i)
780  {
781  str.append(QString("(%1, %2, %3)").arg(positions[i]).arg(values[i]).arg(sigmas[i]));
782  if (i < positions.size() - 1)
783  str.append(", ");
784  }
785  str.append(QString("];iterations=%1").arg(numSteps));
786  str.append(QString(";duration=%1").arg(stopWatch.elapsed() / 1000));
787  str.append(QString(";solution=%1").arg(focusSolution));
788  str.append(QString(";HFR=%1").arg(focusHFR));
789  str.append(QString(";filter='%1'").arg(params.filterName));
790  str.append(QString(";temperature=%1").arg(params.temperature));
791  str.append(QString(";focusalgorithm=%1").arg(params.focusAlgorithm));
792  str.append(QString(";backlash=%1").arg(params.backlash));
793  str.append(QString(";curvefit=%1").arg(params.curveFit));
794  str.append(QString(";useweights=%1").arg(params.useWeights));
795 
796  qCDebug(KSTARS_EKOS_FOCUS) << str;
797 }
798 
799 // Called by Linear to set state for the second pass
800 int LinearFocusAlgorithm::setupSecondPass(int position, double value, double margin)
801 {
802  firstPassBestPosition = position;
803  firstPassBestValue = value;
804  inFirstPass = false;
805  solutionPending = false;
806  secondPassStartIndex = values.size();
807 
808  // Arbitrarily go back "margin" steps above the best position.
809  // Could be a problem if backlash were worse than that many steps.
810  requestedPosition = std::min(static_cast<int>(firstPassBestPosition + stepSize * margin), maxPositionLimit);
811  stepSize = params.initialStepSize / 2;
812  qCDebug(KSTARS_EKOS_FOCUS) << QString("Linear: 2ndPass starting at %1 step %2").arg(requestedPosition).arg(stepSize);
813  return requestedPosition;
814 }
815 
816 // Called by L1P to finish off the first pass by setting state variables
817 int LinearFocusAlgorithm::finishFirstPass(int position, double value)
818 {
819  firstPassBestPosition = position;
820  firstPassBestValue = value;
821  inFirstPass = false;
822  requestedPosition = position;
823  return requestedPosition;
824 }
825 
826 // Return true if one of the 2 recent samples is among the best 2 samples so far.
827 bool LinearFocusAlgorithm::bestSamplesHeuristic()
828 {
829  const int length = values.size();
830  if (length < 5) return true;
831  QVector<double> tempValues = values;
832  std::nth_element(tempValues.begin(), tempValues.begin() + 2, tempValues.end());
833  double secondBest = tempValues[1];
834  if ((values[length - 1] <= secondBest) || (values[length - 2] <= secondBest))
835  return true;
836  return false;
837 }
838 
839 // Return true if there are "streak" consecutive values which are successively worse.
840 bool LinearFocusAlgorithm::gettingWorse()
841 {
842  // Must have this many consecutive values getting worse.
843  constexpr int streak = 3;
844  const int length = values.size();
845  if (secondPassStartIndex < 0)
846  return false;
847  if (length < streak + 1)
848  return false;
849  // This insures that all the values we're checking are in the latest 2nd pass.
850  if (length - secondPassStartIndex < streak + 1)
851  return false;
852  for (int i = length - 1; i >= length - streak; --i)
853  if (values[i] <= values[i - 1])
854  return false;
855  return true;
856 }
857 
858 }
859 
Ekos is an advanced Astrophotography tool for Linux. It is based on a modular extensible framework to...
Definition: align.cpp:70
QVector::iterator begin()
int size() const const
QString i18n(const char *text, const TYPE &arg...)
QVector::iterator end()
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
QVector< V > values(const QMultiHash< K, V > &c)
T value(int i) const const
QString & append(QChar ch)
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:54 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.