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

KDE's Doxygen guidelines are available online.