Functions
cpl_error_code	cpl_fit_image_gaussian (const cpl_image im, const cpl_image im_err, cpl_size xpos, cpl_size ypos, cpl_size xsize, cpl_size ysize, cpl_array parameters, cpl_array err_params, const cpl_array fit_params, double rms, double red_chisq, cpl_matrix covariance, double major, double minor, double angle, cpl_matrix **phys_cov)
	Fit a 2D gaussian to image values.

cpl_imagelist *	cpl_fit_imagelist_polynomial (const cpl_vector x_pos, const cpl_imagelist values, cpl_size mindeg, cpl_size maxdeg, cpl_boolean is_symsamp, cpl_type pixeltype, cpl_image *fiterror)
	Least-squares fit a polynomial to each pixel in a list of images.

cpl_imagelist *	cpl_fit_imagelist_polynomial_window (const cpl_vector x_pos, const cpl_imagelist values, cpl_size llx, cpl_size lly, cpl_size urx, cpl_size ury, cpl_size mindeg, cpl_size maxdeg, cpl_boolean is_symsamp, cpl_type pixeltype, cpl_image *fiterror)
	Least-squares fit a polynomial to each pixel in a list of images.

cpl_error_code	cpl_fit_lvmq (const cpl_matrix x, const cpl_matrix sigma_x, const cpl_vector y, const cpl_vector sigma_y, cpl_vector a, const int ia[], int(f)(const double x[], const double a[], double result), int(dfda)(const double x[], const double a[], double result[]), double relative_tolerance, int tolerance_count, int max_iterations, double mse, double red_chisq, cpl_matrix **covariance)
	Fit a function to a set of data.

double	cpl_gaussian_eval_2d (const cpl_array *self, double x, double y)
	Evaluate the Gaussian in a 2D-point.

Detailed Description

This module provides a routine for non-linear fitting.

Synopsis:: #include "cpl_fit.h"

Function Documentation

◆ cpl_fit_image_gaussian()

cpl_error_code cpl_fit_image_gaussian	(	const cpl_image *	im,
		const cpl_image *	im_err,
		cpl_size	xpos,
		cpl_size	ypos,
		cpl_size	xsize,
		cpl_size	ysize,
		cpl_array *	parameters,
		cpl_array *	err_params,
		const cpl_array *	fit_params,
		double *	rms,
		double *	red_chisq,
		cpl_matrix **	covariance,
		double *	major,
		double *	minor,
		double *	angle,
		cpl_matrix **	phys_cov
	)

Fit a 2D gaussian to image values.

Parameters

im	Input image with data values to fit.
im_err	Optional input image with statistical errors associated to data.
xpos	X position of center of fitting domain.
ypos	Y position of center of fitting domain.
xsize	X size of fitting domain. It must be at least 3 pixels.
ysize	Y size of fitting domain. It must be at least 3 pixels.
parameters	Preallocated array for returning the values of the best-fit gaussian parameters (the parametrisation of the fitted gaussian is described in the main documentation section, below). This array must be of type CPL_TYPE_DOUBLE, and it must have exactly 7 elements. Generally, when passed to this function, this array would not be initialised (all elements are "invalid"). A first-guess for the gaussian parameters is not mandatory: but it is possible to specify here a first-guess value for each parameter. First-guess values can also be specified just for a subset of parameters.
err_params	Optional preallocated array for returning the statistical error associated to each fitted parameter. This array must be of type CPL_TYPE_DOUBLE, and it must have exactly 7 elements. This makes mandatory to specify im_err. Note that the returned values are the square root of the diagonal elements (variances) of the covariance matrix (see ahead).
fit_params	Optional array, used for flagging the parameters to freeze. This array must be of type CPL_TYPE_INT, and it must have exactly 7 elements. If an array element is set to 0, the corresponding parameter will be frozen. Any other value (including an "invalid" array element) would indicate a free parameter. If a parameter is frozen, a first-guess value must be specified at the corresponding element of the parameters array. If no array is specified here (NULL pointer), all parameters are free.
rms	If not NULL, returned standard deviation of fit residuals.
red_chisq	If not NULL, returned reduced chi-squared of fit. This makes mandatory to specify im_err.
covariance	If not NULL, a newly allocated covariance matrix will be returned. This makes mandatory to specify im_err. On error it is not modified.
major	If not NULL, returned semi-major axis of ellipse at 1-sigma.
minor	If not NULL, returned semi-minor axis of ellipse at 1-sigma.
angle	If not NULL, returned angle between X axis and major axis of ellipse, counted counterclockwise (radians).
phys_cov	If not NULL, a newly allocated 3x3 covariance matrix for the derived physical parameters major, minor, and angle, will be returned. This makes mandatory to specify im_err. On error it is not modified.

Returns: CPL_ERROR_NONE on successful fit.

This function fits a 2d gaussian to pixel values within a specified region by minimizing $\chi^2$ using a Levenberg-Marquardt algorithm. The gaussian model adopted here is based on the well-known cartesian form

$z = B + \frac{A}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp\left({-\frac{1}{2\left(1-\rho^2\right)} \left(\left(\frac{x - \mu_x}{\sigma_x}\right)^2 -2\rho\left(\frac{x - \mu_x}{\sigma_x}\right) \left(\frac{y - \mu_y}{\sigma_y}\right) + \left(\frac{y - \mu_y}{\sigma_y}\right)^2\right)}\right)$

where is a background level and the volume of the gaussian (they both can be negative!), making 7 parameters altogether. Conventionally the parameters are indexed from 0 to 6 in the elements of the arrays parameters, err_params, fit_params, and of the 7x7 covariance matrix:

$\begin{eqnarray*} \mathrm{parameters[0]} &=& B \\ \mathrm{parameters[1]} &=& A \\ \mathrm{parameters[2]} &=& \rho \\ \mathrm{parameters[3]} &=& \mu_x \\ \mathrm{parameters[4]} &=& \mu_y \\ \mathrm{parameters[5]} &=& \sigma_x \\ \mathrm{parameters[6]} &=& \sigma_y \end{eqnarray*}$

The semi-axes and the orientation $\theta$ of the ellipse at 1-sigma level are finally derived from the fitting parameters as:

$\begin{eqnarray*} \theta &=& \frac{1}{2} \arctan \left(2 \rho \frac{\sigma_x \sigma_y} {\sigma_x^2 - \sigma_y^2}\right) \\ a &=& \sigma_x \sigma_y \sqrt{2(1-\rho^2) \frac{\cos 2\theta} {\left(\sigma_x^2 + \sigma_y^2\right) \cos 2\theta + \sigma_y^2 - \sigma_x^2}} \\ b &=& \sigma_x \sigma_y \sqrt{2(1-\rho^2) \frac{\cos 2\theta} {\left(\sigma_x^2 + \sigma_y^2\right) \cos 2\theta - \sigma_y^2 + \sigma_x^2}} \end{eqnarray*}$

Note that $\theta$ is counted counterclockwise starting from the positive direction of the axis, ranging bewteen $-\pi/2$ and $+\pi/2$ radians.

If the correlation $\rho = 0$ and $\sigma_x \geq \sigma_y$ (within uncertainties) the ellipse is either a circle or its major axis is aligned with the axis, so it is conventionally set

$\begin{eqnarray*} \theta &=& 0 \\ a &=& \sigma_x \\ b &=& \sigma_y \end{eqnarray*}$

If the correlation $\rho = 0$ and $\sigma_x < \sigma_y$ (within uncertainties) the major axis of the ellipse is aligned with the axis, so it is conventionally set

$\begin{eqnarray*} \theta &=& \frac{\pi}{2} \\ a &=& \sigma_y \\ b &=& \sigma_x \end{eqnarray*}$

If requested, the 3x3 covariance matrix G associated to the derived physical quantities is also computed, applying the usual

$\mathrm{G} = \mathrm{J} \mathrm{C} \mathrm{J}^\mathrm{T}$

where J is the Jacobian of the transformation $(B, A, \rho, \mu_x, \mu_y, \sigma_x, \sigma_y) \rightarrow (\theta, a, b)$ and C is the 7x7 matrix of the gaussian parameters.

References cpl_array_delete(), cpl_array_get_double(), cpl_array_get_int(), cpl_array_get_type(), cpl_array_new(), cpl_array_set_double(), cpl_array_set_invalid(), CPL_ERROR_ACCESS_OUT_OF_RANGE, CPL_ERROR_DATA_NOT_FOUND, CPL_ERROR_ILLEGAL_INPUT, CPL_ERROR_INCOMPATIBLE_INPUT, CPL_ERROR_INVALID_TYPE, CPL_ERROR_NONE, CPL_ERROR_NULL_INPUT, CPL_ERROR_SINGULAR_MATRIX, CPL_ERROR_TYPE_MISMATCH, cpl_errorstate_get(), cpl_errorstate_set(), cpl_free(), cpl_image_delete(), cpl_image_extract(), cpl_image_get(), cpl_image_get_bpm_const(), cpl_image_get_data_double_const(), cpl_image_get_maxpos(), cpl_image_get_mean(), cpl_image_get_median(), cpl_image_get_minpos(), cpl_image_get_size_x(), cpl_image_get_size_y(), cpl_image_get_type(), cpl_malloc(), cpl_mask_is_empty(), CPL_MATH_2PI, CPL_MATH_FWHM_SIG, CPL_MATH_PI_2, cpl_matrix_delete(), cpl_matrix_new(), cpl_matrix_set_(), cpl_matrix_unwrap(), cpl_matrix_wrap(), CPL_MAX, CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT, cpl_type_get_name(), CPL_TYPE_INT, cpl_vector_delete(), cpl_vector_unwrap(), and cpl_vector_wrap().

◆ cpl_fit_imagelist_polynomial()

cpl_imagelist * cpl_fit_imagelist_polynomial	(	const cpl_vector *	x_pos,
		const cpl_imagelist *	values,
		cpl_size	mindeg,
		cpl_size	maxdeg,
		cpl_boolean	is_symsamp,
		cpl_type	pixeltype,
		cpl_image *	fiterror
	)

Least-squares fit a polynomial to each pixel in a list of images.

Parameters

x_pos	The vector of positions to fit
values	The list of images with values to fit
mindeg	The smallest degree with a non-zero coefficient
maxdeg	The polynomial degree of the fit, at least mindeg
is_symsamp	True iff the x_pos values are symmetric around their mean
pixeltype	The pixel-type of the created image list
fiterror	When non-NULL, the error of the fit

Note: values and x_pos must have the same number of elements.; The created imagelist must be deallocated with cpl_imagelist_delete().; x_pos must have at least 1 + (maxdeg - mindeg) distinct values.

Returns: The image list of the fitted polynomial coefficients or NULL on error.

See also: cpl_fit_imagelist_polynomial_window()

Possible _cpl_error_code_ set in this function:

CPL_ERROR_NULL_INPUT if an input const pointer is NULL
CPL_ERROR_ILLEGAL_INPUT if mindeg is negative or maxdeg is less than mindeg.
CPL_ERROR_INCOMPATIBLE_INPUT if x_pos and values have different lengths, or if fiterror is non-NULL with a different size than that of values, or if the input images do not all have the same dimensions and pixel type.
CPL_ERROR_DATA_NOT_FOUND if x_pos contains less than nc values.
CPL_ERROR_SINGULAR_MATRIX if x_pos contains less than nc distinct values.
CPL_ERROR_UNSUPPORTED_MODE if the chosen pixel type is not one of CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT, CPL_TYPE_INT.

References cpl_fit_imagelist_polynomial_window(), cpl_image_get_size_x(), cpl_image_get_size_y(), and cpl_imagelist_get_const().

◆ cpl_fit_imagelist_polynomial_window()

cpl_imagelist * cpl_fit_imagelist_polynomial_window	(	const cpl_vector *	x_pos,
		const cpl_imagelist *	values,
		cpl_size	llx,
		cpl_size	lly,
		cpl_size	urx,
		cpl_size	ury,
		cpl_size	mindeg,
		cpl_size	maxdeg,
		cpl_boolean	is_symsamp,
		cpl_type	pixeltype,
		cpl_image *	fiterror
	)

Least-squares fit a polynomial to each pixel in a list of images.

Parameters

x_pos	The vector of positions to fit
values	The list of images with values to fit
llx	Lower left x position (FITS convention, 1 for leftmost)
lly	Lower left y position (FITS convention, 1 for lowest)
urx	Upper right x position (FITS convention)
ury	Upper right y position (FITS convention)
mindeg	The smallest degree with a non-zero coefficient
maxdeg	The polynomial degree of the fit, at least mindeg
is_symsamp	True iff the x_pos values are symmetric around their mean
pixeltype	The (non-complex) pixel-type of the created image list
fiterror	When non-NULL, the error of the fit. Must be non-complex

Note: values and x_pos must have the same number of elements.; The created imagelist must be deallocated with cpl_imagelist_delete().; x_pos must have at least 1 + (maxdeg - mindeg) distinct values.

Returns: The image list of the fitted polynomial coefficients or NULL on error.

See also: cpl_polynomial_fit()

For each pixel, a polynomial representing the relation value = P(x) is computed where: P(x) = x^{mindeg} * (a_0 + a_1 * x + ... + a_{nc-1} * x^{nc-1}), where mindeg >= 0 and maxdeg >= mindeg, and nc is the number of polynomial coefficients to determine, nc = 1 + (maxdeg - mindeg).

The returned image list thus contains nc coefficient images, a_0, a_1, ..., a_{nc-1}.

np is the number of sample points, i.e. the number of elements in x_pos and number of images in the input image list.

If mindeg is nonzero then is_symsamp is ignored, otherwise is_symsamp may to be set to CPL_TRUE if and only if the values in x_pos are known a-priori to be symmetric around their mean, e.g. (1, 2, 4, 6, 10, 14, 16, 18, 19), but not (1, 2, 4, 6, 10, 14, 16). Setting is_symsamp to CPL_TRUE while mindeg is zero eliminates certain round-off errors. For higher order fitting the fitting problem known as "Runge's phenomenon" is minimized using the socalled "Chebyshev nodes" as sampling points. For Chebyshev nodes is_symsamp can be set to CPL_TRUE.

Even though it is not an error, it is hardly useful to use an image of pixel type integer for the fitting error. An image of pixel type float should on the other hand be sufficient for most fitting errors.

The call requires the following number of FLOPs, where nz is the number of pixels in any one image in the imagelist:

2 * nz * nc * (nc + np) + np * nc^2 + nc^3/3 + O(nc * (nc + np)).

If mindeg is zero an additional nz * nc^2 FLOPs are required.

If fiterror is non-NULL an additional 2 * nz * nc * np FLOPs are required.

Bad pixels in the input is suported as follows: First all pixels are fitted ignoring any bad pixel maps in the input. If this succeeds then each fit, where bad pixel(s) are involved is redone. During this second pass all input pixels flagged as bad are ignored. For each pixel to be redone, the remaining good samples are passed to cpl_polynomial_fit(). The input is_symsamp is ignored in this second pass. The reduced number of samples may reduce the number of sampling points to equal the number of coefficients to fit. In this case the fit has another meaning (any non-zero residual is due to rounding errors, not a fitting error). If for a given fit bad pixels reduces the number of sampling points to less than the number of coefficients to fit, then as many coefficients are fit as there are sampling points. The higher order coefficients are set to zero and flagged as bad. If a given pixel has no good samples, then the resulting fit will consist of zeroes, all flagged as bad.

Possible _cpl_error_code_ set in this function:

CPL_ERROR_NULL_INPUT if an input const pointer is NULL
CPL_ERROR_ILLEGAL_INPUT if mindeg is negative or maxdeg is less than mindeg or if llx or lly are smaller than 1 or if urx or ury is smaller than llx and lly respectively.
CPL_ERROR_ACCESS_OUT_OF_RANGE if urx or ury exceed the size of values.
CPL_ERROR_INCOMPATIBLE_INPUT if x_pos and values have different lengths, or if fiterror is non-NULL with a different size than that of values, or if the input images do not all have the same dimensions and pixel type.
CPL_ERROR_DATA_NOT_FOUND if x_pos contains less than nc values.
CPL_ERROR_SINGULAR_MATRIX if x_pos contains less than nc distinct values.
CPL_ERROR_UNSUPPORTED_MODE if the chosen pixel type is not one of CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT, CPL_TYPE_INT.

References cpl_ensure, CPL_ERROR_ACCESS_OUT_OF_RANGE, CPL_ERROR_DATA_NOT_FOUND, CPL_ERROR_ILLEGAL_INPUT, CPL_ERROR_INCOMPATIBLE_INPUT, CPL_ERROR_NONE, CPL_ERROR_NULL_INPUT, CPL_ERROR_SINGULAR_MATRIX, CPL_ERROR_UNSUPPORTED_MODE, cpl_image_get_bpm_const(), cpl_image_get_size_x(), cpl_image_get_size_y(), cpl_imagelist_delete(), cpl_imagelist_get_const(), cpl_imagelist_get_size(), cpl_imagelist_is_uniform(), cpl_imagelist_new(), cpl_imagelist_set(), cpl_malloc(), cpl_mask_delete(), cpl_mask_or(), CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT, cpl_type_get_sizeof(), CPL_TYPE_INT, and cpl_vector_get_size().

Referenced by cpl_fit_imagelist_polynomial().

◆ cpl_fit_lvmq()

cpl_error_code cpl_fit_lvmq	(	const cpl_matrix *	x,
		const cpl_matrix *	sigma_x,
		const cpl_vector *	y,
		const cpl_vector *	sigma_y,
		cpl_vector *	a,
		const int	ia[],
		int()(const double x[], const double a[], double result)	f,
		int(*)(const double x[], const double a[], double result[])	dfda,
		double	relative_tolerance,
		int	tolerance_count,
		int	max_iterations,
		double *	mse,
		double *	red_chisq,
		cpl_matrix **	covariance
	)

Fit a function to a set of data.

Parameters

x	N x D matrix of the positions to fit. Each matrix row is a D-dimensional position.
sigma_x	Uncertainty (one sigma, gaussian errors assumed) assosiated with x. Taking into account the uncertainty of the independent variable is currently unsupported, and this parameter must therefore be set to NULL.
y	The N values to fit.
sigma_y	Vector of size N containing the uncertainties of the y-values. If this parameter is NULL, constant uncertainties are assumed.
a	Vector containing M fit parameters. Must contain a guess solution on input and contains the best fit parameters on output.
ia	Array of size M defining which fit parameters participate in the fit (non-zero) and which fit parameters are held constant (zero). At least one element must be non-zero. Alternatively, pass NULL to fit all parameters.
f	Function that evaluates the fit function at the position specified by the first argument (an array of size D) using the fit parameters specified by the second argument (an array of size M). The result must be output using the third parameter, and the function must return zero iff the evaluation succeded.
dfda	Function that evaluates the first order partial derivatives of the fit function with respect to the fit parameters at the position specified by the first argument (an array of size D) using the parameters specified by the second argument (an array of size M). The result must be output using the third parameter (array of size M), and the function must return zero iff the evaluation succeded.
relative_tolerance	The algorithm converges by definition if the relative decrease in chi squared is less than tolerance tolerance_count times in a row. Recommended default: CPL_FIT_LVMQ_TOLERANCE
tolerance_count	The algorithm converges by definition if the relative decrease in chi squared is less than tolerance tolerance_count times in a row. Recommended default: CPL_FIT_LVMQ_COUNT
max_iterations	If this number of iterations is reached without convergence, the algorithm diverges, by definition. Recommended default: CPL_FIT_LVMQ_MAXITER
mse	If non-NULL, the mean squared error of the best fit is computed.
red_chisq	If non-NULL, the reduced chi square of the best fit is computed. This requires sigma_y to be specified.
covariance	If non-NULL, the formal covariance matrix of the best fit parameters is computed (or NULL on error). On success the diagonal terms of the covariance matrix are guaranteed to be positive. However, terms that involve a constant parameter (as defined by the input array ia) are always set to zero. Computation of the covariacne matrix requires sigma_y to be specified.

Returns: CPL_ERROR_NONE iff OK.

This function makes a minimum chi squared fit of the specified function to the specified data set using a Levenberg-Marquardt algorithm.

Possible _cpl_error_code_ set in this function:

CPL_ERROR_NULL_INPUT if an input pointer other than sigma_x, sigma_y, mse, red_chisq or covariance is NULL.
CPL_ERROR_ILLEGAL_INPUT if an input matrix/vector is empty, if ia contains only zero values, if any of relative_tolerance, tolerance_count or max_iterations is non-positive, if N <= M and red_chisq is non-NULL, if any element of sigma_x or sigma_y is non-positive, or if evaluation of the fit function or its derivative failed.
CPL_ERROR_INCOMPATIBLE_INPUT if the dimensions of the input vectors/matrices do not match, or if chi square or covariance computation is requested and sigma_y is NULL.
CPL_ERROR_ILLEGAL_OUTPUT if memory allocation failed.
CPL_ERROR_CONTINUE if the Levenberg-Marquardt algorithm failed to converge.
CPL_ERROR_SINGULAR_MATRIX if the covariance matrix could not be computed.