# Kstars

curvefit.h
1 #pragma once
2
3 #include "../../auxiliary/robuststatistics.h"
4
5 #include <QVector>
6 #include <qcustomplot.h>
7 #include <gsl/gsl_vector.h>
8 #include <gsl/gsl_min.h>
9 #include <gsl/gsl_matrix.h>
10 #include <gsl/gsl_multifit.h>
11 #include <gsl/gsl_multifit_nlinear.h>
12 #include <gsl/gsl_blas.h>
13 #include <ekos_focus_debug.h>
14
15 namespace Ekos
16 {
17 // The curve fitting class provides curve fitting algorithms using the Lehvensberg-Marquart (LM)
18 // solver as provided the Gnu Science Library (GSL). See the comments at the start of curvefit.cpp
19 // for more details.
20 //
21 // 2 curves, hyperbola and parabola are provided.
22 // For compatibility with existing Ekos functionality a Quadratic option using the exising Ekos linear
23 // least squares solver (again provided by GSL) is supported. The Quadratic and Parabola curves are
24 // the same thing mathematically but Parabola uses the non-linear least squares LM solver whilst Quadratic
25 // uses the original Ekos linear least squares solver.
26 //
27 // Users of CurveFitting operate on focuser position and HFR. Within CurveFitting the curve uses the more
28 // usual mathematical notation of x, y
29 //
30 // Furture releases may migrate all curve fitting to the LM solver.
31 class CurveFitting
32 {
33  public:
34  typedef enum { FOCUS_QUADRATIC, FOCUS_HYPERBOLA, FOCUS_PARABOLA, FOCUS_GAUSSIAN } CurveFit;
35  typedef enum { OPTIMISATION_MINIMISE, OPTIMISATION_MAXIMISE } OptimisationDirection;
36  typedef enum { STANDARD, BEST, BEST_RETRY } FittingGoal;
37
38  // Data structures to hold datapoints for the class
39  struct DataPT
40  {
41  double x; // focuser position
42  double y; // focus measurement, e.g. HFR
43  double weight; // the measurement weight, e.g. inverse variance
44  };
45
46  struct DataPointT // This is the data strcture passed into GSL LM routines
47  {
48  bool useWeights; // Are we fitting the curve with or without weights
49  QVector<DataPT> dps; // Vector of DataPT
50  OptimisationDirection dir = OPTIMISATION_MINIMISE; // Are we minimising or maximising?
51  Mathematics::RobustStatistics::ScaleCalculation weightCalculation =
52  Mathematics::RobustStatistics::SCALE_VARIANCE; // How to compute weights
53
54  // Helper functions to operate on the data structure
55  void push_back(double x, double y, double weight = 1)
56  {
57  dps.push_back({x, y, weight});
58  }
59  };
60
61  struct DataPT3D
62  {
63  double x; // x position
64  double y; // y position
65  int z; // value
66  double weight; // the measurement weight, e.g. inverse variance
67  };
68
69  struct DataPoint3DT // This is the data strcture passed into GSL LM routines
70  {
71  bool useWeights; // Are we fitting the curve with or without weights
72  QVector<DataPT3D> dps; // Vector of DataPT3D
73  OptimisationDirection dir = OPTIMISATION_MAXIMISE; // Are we minimising or maximising?
74  Mathematics::RobustStatistics::ScaleCalculation weightCalculation =
75  Mathematics::RobustStatistics::SCALE_VARIANCE; // How to compute weights
76
77  // Helper function to operate on the data structure
78  void push_back(double x, double y, int z, double weight = 1)
79  {
80  dps.push_back({x, y, z, weight});
81  }
82  };
83
84  // Data structure to hold star parameters
85  struct StarParams
86  {
87  double background;
88  double peak;
89  double centroid_x;
90  double centroid_y;
91  double HFR;
92  double theta;
93  double FWHMx;
94  double FWHMy;
95  double FWHM;
96  };
97
98  // Constructor just initialises the object
99  CurveFitting();
100
101  // A constructor that reconstructs a previously built object.
102  // Does not implement getting the original data points.
103  CurveFitting(const QString &serialized);
104
105  // fitCurve takes in the vectors with the position, hfr and weight (e.g. variance in HFR) values
106  // along with the type of curve to use and whether or not to use weights in the calculation
107  // It fits the curve and solves for the coefficients.
108  void fitCurve(const FittingGoal goal, const QVector<int> &position, const QVector<double> &hfr,
109  const QVector<double> &weights, const QVector<bool> &outliers,
110  const CurveFit curveFit, const bool useWeights, const OptimisationDirection optDir);
111
112  // fitCurve3D 3-Dimensional version of fitCurve.
113  // Data is passed in in imageBuffer - a 2D array of width x height
114  // Approx star information is passed in to seed the LM solver initial parameters.
115  // Start and end define the x,y coordinates of a box around the star start is top left corner, end is bottom right
116  template <typename T>
117  void fitCurve3D(const T *imageBuffer, const int imageWidth, const QPair<int, int> start, const QPair<int, int> end,
118  const StarParams &starParams, const CurveFit curveFit, const bool useWeights)
119  {
120  if (imageBuffer == nullptr)
121  {
122  qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D null image ptr");
123  m_FirstSolverRun = true;
124  return;
125  }
126
127  if (imageWidth <= 0)
128  {
129  qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D imageWidth=%1").arg(imageWidth);
130  m_FirstSolverRun = true;
131  return;
132  }
133
134  if (useWeights)
135  {
136  qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with useWeights. Not yet implemented");
137  m_FirstSolverRun = true;
138  return;
139  }
140
141  m_dataPoints.dps.clear();
142  m_dataPoints.useWeights = useWeights;
143
144  // Load up the data structures for the solver.
145  // The pixel reference x, y refers to the top left corner of the pizel so add 0.5 to x and y to reference the
146  // centre of the pixel.
147  int width = end.first - start.first;
148  int height = end.second - start.second;
149
150  for (int j = 0; j < height; j++)
151  for (int i = 0; i < width; i++)
152  m_dataPoints.push_back(i + 0.5, j + 0.5, imageBuffer[start.first + i + ((start.second + j) * imageWidth)], 1.0);
153
154  m_CurveType = curveFit;
155  switch (m_CurveType)
156  {
157  case FOCUS_GAUSSIAN :
158  m_coefficients = gaussian_fit(m_dataPoints, starParams);
159  break;
160  default :
161  // Something went wrong, log an error and reset state so solver starts from scratch if called again
162  qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with curveFit=%1").arg(curveFit);
163  m_FirstSolverRun = true;
164  return;
165  }
166  m_LastCoefficients = m_coefficients;
167  m_LastCurveType = m_CurveType;
168  m_FirstSolverRun = false;
169  }
170
171  // Returns the minimum position and value in the pointers for the solved curve.
172  // Returns false if the curve couldn't be solved
173  bool findMinMax(double expected, double minPosition, double maxPosition, double *position, double *value, CurveFit curveFit,
174  const OptimisationDirection optDir);
175  // getStarParams returns the star parameters for the solved star
176  bool getStarParams(const CurveFit curveFit, StarParams *starParams);
177
178  // getCurveParams gets the coefficients of a curve solve
179  // setCurveParams sets the coefficients of a curve solve
180  // using get and set returns the solver to its state as it was when get was called
181  // This allows functions like "f" to be called to calculate curve values for a
182  // prior curve fit solution.
183  bool getCurveParams(const CurveFit curveType, QVector<double> &coefficients)
184  {
185  if (curveType != m_CurveType)
186  return false;
187  coefficients = m_coefficients;
188  return true;
189  }
190
191  bool setCurveParams(const CurveFit curveType, const QVector<double> coefficients)
192  {
193  if (curveType != m_CurveType)
194  return false;
195  m_coefficients = coefficients;
196  return true;
197  }
198
199  // f calculates the value of y for a given x using the appropriate curve algorithm
200  double f(double x);
201  // f calculates the value of z for a given x and y using the appropriate curve algorithm
202  double f3D(double x, double y);
203  // Calculates the R-squared which is a measure of how well the curve fits the datapoints
204  double calculateR2(CurveFit curveFit);
205  // Calculate the deltas of each datapoint from the curve
206  void calculateCurveDeltas(CurveFit curveFit, std::vector<std::pair<int, double>> &curveDeltas);
207
208  // Returns a QString which can be used to construct an identical object.
209  // Does not implement getting the original data points.
210  QString serialize() const;
211
212  private:
213  // TODO: This function will likely go when Linear and L1P merge to be closer.
214  // Calculates the value of the polynomial at x. Params will be cast to a CurveFit*.
215  static double curveFunction(double x, void *params);
216
217  // TODO: This function will likely go when Linear and L1P merge to be closer.
218  QVector<double> polynomial_fit(const double *const data_x, const double *const data_y, const int n, const int order);
219
220  QVector<double> hyperbola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
221  const QVector<double> weights, bool useWeights, const OptimisationDirection optDir);
222  QVector<double> parabola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
223  const QVector<double> data_weights,
224  bool useWeights, const OptimisationDirection optDir);
225  QVector<double> gaussian_fit(DataPoint3DT data, const StarParams &starParams);
226
227  bool minimumQuadratic(double expected, double minPosition, double maxPosition, double *position, double *value);
228  bool minMaxHyperbola(double expected, double minPosition, double maxPosition, double *position, double *value,
229  const OptimisationDirection optDir);
230  bool minMaxParabola(double expected, double minPosition, double maxPosition, double *position, double *value,
231  const OptimisationDirection optDir);
232  bool getGaussianParams(StarParams *starParams);
233
234  void hypMakeGuess(const int attempt, const QVector<double> inX, const QVector<double> inY,
235  const OptimisationDirection optDir, gsl_vector * guess);
236  void hypSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
237  double *ftol);
238
239  void parMakeGuess(const int attempt, const QVector<double> inX, const QVector<double> inY,
240  const OptimisationDirection optDir,
241  gsl_vector * guess);
242  void parSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
243  double *ftol);
244  void gauMakeGuess(const int attempt, const StarParams &starParams, gsl_vector * guess);
245  void gauSetupParams(gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol, double *ftol);
246
247  // Get the reason code from the passed in info
248  QString getLMReasonCode(int info);
249
250  // Calculation engine for the R-squared which is a measure of how well the curve fits the datapoints
251  double calcR2(const QVector<double> dataPoints, const QVector<double> curvePoints, const QVector<double> scale,
252  const bool useWeights);
253
254  // Used in the QString constructor.
255  bool recreateFromQString(const QString &serialized);
256
257  // Type of curve
258  CurveFit m_CurveType;
259  // The data values.
260  QVector<double> m_x, m_y, m_scale;
261  // Use weights or not
262  bool m_useWeights;
263  DataPoint3DT m_dataPoints;
264  // The solved parameters.
265  QVector<double> m_coefficients;
266  // State variables used by the LM solver. These variables provide a way of optimising the starting
267  // point for the solver by using the solution found by the previous run providing the relevant
268  // solver parameters are consistent between runs.
269  bool m_FirstSolverRun;
270  CurveFit m_LastCurveType;
271  QVector<double> m_LastCoefficients;
272 };
273
274 } //namespace
Ekos is an advanced Astrophotography tool for Linux. It is based on a modular extensible framework to...
Definition: align.cpp:69
void push_back(const T &value)
KIOFILEWIDGETS_EXPORT QString dir(const QString &fileClass)
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
const QList< QKeySequence > & end()
This file is part of the KDE documentation.