Peak Centroiding

Macros
#define	_POSIX_C_SOURCE   200112L

Functions
cpl_error_code	clipm_centroiding_gauss (const cpl_image *image, const cpl_size window_xxyy[4], double *xy_centre, double *xy_centre_err, double *xy_sigma, double *xy_sigma_err, double *xy_fwhm, double *xy_fwhm_err, double *centre_intensity, int robustness)
	Determine the position of an object in an image window, for x and y separately, by collapsing the image and fitting a gaussian. More...
cpl_error_code	clipm_centroiding_moment (const cpl_image *image, const cpl_size window_xxyy[4], int allow_wdw_enlarge, double gain, double *out_xy_centre, double *out_xy_centre_err, double *out_xy_sigma, double *out_xy_fwhm, double *centre_intensity)
	Determine the barycenter of an object in an image window. More...
cpl_error_code	clipm_centroiding_multi_gauss (const cpl_image *image, const cpl_matrix *locations, cpl_size areasize, cpl_matrix **xy_centre, cpl_matrix **xy_centre_err, cpl_matrix **xy_sigma, cpl_matrix **xy_sigma_err, cpl_matrix **xy_fwhm, cpl_matrix **xy_fwhm_err, cpl_matrix **centre_intensities, cpl_array **all_error_codes, int robustness)
	Invoke clipm_centroiding_gauss() at several locations in an image. More...

Detailed Description

This module provides functions for centroiding inside image windows.

Synopsis:

#include "clipm_centroiding.h"

Macro Definition Documentation

_POSIX_C_SOURCE

#define _POSIX_C_SOURCE   200112L

Function Documentation

clipm_centroiding_gauss()

cpl_error_code clipm_centroiding_gauss	(	const cpl_image *	image,
		const cpl_size	window_xxyy[4],
		double *	xy_centre,
		double *	xy_centre_err,
		double *	xy_sigma,
		double *	xy_sigma_err,
		double *	xy_fwhm,
		double *	xy_fwhm_err,
		double *	centre_intensity,
		int	robustness
	)

Determine the position of an object in an image window, for x and y separately, by collapsing the image and fitting a gaussian.

Parameters

image	Input image (FITS convention)
window_xxyy	Coordinate buffer of the form {xa, xb, ya, yb}, can be NULL, minimum/maximum order is irrelevant
xy_centre	Output centre coordinate buffer of size 2 ([x, y]), can be NULL
xy_centre_err	Output centre error buffer (1 sigma) of size 2 ([x_err, y_err]), can be NULL
xy_sigma	Output sigma buffer of size 2 ([sx, sy]), can be NULL
xy_sigma_err	Output sigma error buffer (1 sigma) of size 2 ([sx_err, sy_err]), can be NULL
xy_fwhm	Output full-width-half-maximum buffer of size 2 ([fwhm_x, fwhm_y]), can be NULL
xy_fwhm_err	Output FWHM error buffer (1 sigma) of size 2 ([fwhm_x_err, fwhm_y_err]), can be NULL
centre_intensity	Output intensity of the (closest) centre pixel, can be NULL
robustness	Maximum number of retries (>= 0)

Returns

CPL error code

Overview:

The central position of an enclosed object is computed by fitting a gaussian to the marginal distributions respectively in both dimensions.

The centre (in x and y) is returned in xy_centre, the sigma of the gaussian in xy_sigma, and the full-width-half-maximum in xy_fwhm.

The measurement errors of each of these values are returned as 1 sigma in xy_centre_err, xy_sigma_err, and xy_fwhm_err respectively. To get the error variance, these values must be squared by the user.

If the fitting was successful, then the intensity of the pixel which is closest to the center is returned in centre_intensity.

window contains at the end the last used window coordinates, if robustness > 0 (see below).

Robustness:

If robustness > 0 and the fitting fails in the first try, then the window is resized automatically (around a centre guess, see below) and the fitting is retried. This is repeated a maximum of robustness + 1 times.
In this case, window contains at the end the last used window coordinates.

The fitting is generally first repeated maximal robustness times until it is successful. Then (if robustness > 0):

• if it was a failure, then just one other try is added,
• if it was successful, one final fitting process is again performed in the +- 3 sigma region in x and y around the found centroid for improving the result.

Reasonable values for robustness are:

• 0: for one pure gaussian fit (with exactly the given window, if provided), and no retries
• 1..2: for good signal/noise ratio
• < ca. 10: for weak signals, eventually side peaks
• >= ca. 10: other strange cases? None known so far.

If there are other peaks than the main object to be centered, it might be possible that the fitting locks onto them. With enough difference in weight, this should seldom happen. For the case that there are two peaks with diagonal offset, the x and y fit could lock onto the different peaks. This would be detected by a not successful final fitting process or strange results.

Fitting Constraints:

• It is possible that a successful fit is returned. This does not necessarily mean that the correct centroid has been found. Therefore, always the errors (uncertainties) should be evaluated and checked, whether these are in a tolerable range(!!!).
• If the sigma-error (sigma-uncertainty) is greater than half the sigma, then the fitting is interpreted as a failure (and evt. retried depending on robustness).

What It Does Not:

This function was designed to fit a gaussian object with a flat background and limited noise. Searching gaussian peaks in difficult patterns is not what it does and can.

Note

As a rule of thumb, to get low measurement errors, the window or image should enclose at least the 3 sigma region of the fitted gaussian, i.e. for example the "whole" airy disk of a star.

Bad Pixel Handling:

Bad pixel maps are supported. This means that bad pixels are omitted during computation of the marginal distributions, this means during averaging. If this is not desired, bad pixels should be interpolated before calling this function.

Error Handling:

• The following error codes can be set and returned:
  • CPL_ERROR_NULL_INPUT if image is NULL
  • CPL_ERROR_ILLEGAL_INPUT if:
    • the size of window or image in any dimension is less then 5
    • robustness < 0
  • CPL_ERROR_ACCESS_OUT_OF_RANGE if window contains outlying coordinates
  • CPL_ERROR_INVALID_TYPE if image is not of type int, float, double
  • If the fitting fails finally, one of the following codes is returned:
    • CPL_ERROR_ILLEGAL_OUTPUT if the centre is outside window or image (this function is intended for centroiding on existing objects, otherwise have a look at pure fitting functions like e.g. cpl_vector_fit_gaussian())
    • CPL_ERROR_CONTINUE:
      • if obviously a local but not global minimum was found, or
      • if the fitting algorithm could not converge
    • CPL_ERROR_SINGULAR_MATRIX if the covariance matrix could not be computed
• If a fitting of at least one dimension fails to converge (i.e. the returned error is CPL_ERROR_CONTINUE), then the following values are returned for the respective dimension (x or y):
  • centre: the median position after estimating the signal background
  • sigma: the median of the absolute residuals multiplied by 1.4828
  • fwhm: sigma * 2.35482 ( )
  • centre_err: -1
  • sigma_err: -1
  • fwhm_err: -1
  • centre_intensity: not computed
• In all other error cases (wrong input parameters), the output values are undefined.

Todo:

• flag to allow enlargement of window

Example:

#include <cpl.h>
#include <clipm.h>
 
cpl_error_code  myfitting(  cpl_image *image,
                            int x1,
                            int y1,
                            int x2,
                            int y2)
{ 
    double  icent;
    int     window[4] = {x1, x2, y1, y2};
    double  center[2],
            sigma[2],
            center_errorsigma[2],
            sigma_errorsigma[2];
    cpl_error_code
            error;
    const double
            max_center_error = 0.5,
            max_sigma_error = 3;
    
    error = clipm_centroiding_gauss(
                            image,
                            window,
                            center,
                            center_errorsigma,
                            sigma,
                            sigma_errorsigma,
                            NULL,
                            NULL,
                            &icent,
                            7);
    if (error != CPL_ERROR_NONE ||
        center_errorsigma[0] > max_center_error ||
        center_errorsigma[1] > max_center_error ||
        sigma_errorsigma[0] > max_sigma_error ||
        sigma_errorsigma[1] > max_sigma_error)
        printf("Failure\n");
    else {
        printf("Success\n");
        printf("Pixel %f at [%.2f, %.2f]\n",
                            icent,
                            center[0],
                            center[1]);
    }
    
    return error;
}

clipm_centroiding_moment()

cpl_error_code clipm_centroiding_moment	(	const cpl_image *	image,
		const cpl_size	window_xxyy[4],
		int	allow_wdw_enlarge,
		double	gain,
		double *	out_xy_centre,
		double *	out_xy_centre_err,
		double *	out_xy_sigma,
		double *	out_xy_fwhm,
		double *	centre_intensity
	)

Determine the barycenter of an object in an image window.

Parameters

image	Input image (FITS convention)
window_xxyy	Coordinate buffer of the form {xa, xb, ya, yb}, can be NULL, minimum/maximum order is irrelevant
allow_wdw_enlarge	Flag to allow the function to increase the window size when it detects that the window cuts the wings of the centered object, default should be 0
gain	(Optional) detector gain input in unit [1/electrons], required for computation of out_xy_centre_err, ignored if out_xy_centre_err is NULL
out_xy_centre	(Output) centre coordinate buffer of size 2 ([x, y])
out_xy_centre_err	(Optional output) centre uncertainty buffer of size 2 ([x, y]), can be NULL
out_xy_sigma	(Optional output) sigma buffer of size 2 ([x, y]), can be NULL
out_xy_fwhm	(Optional Output) FWHM (buffer of size 2 [x, y]), measured (not derived from sigma)
centre_intensity	(Optional output) intensity of the (closest) centre pixel, can be NULL

Returns

CPL error code

Principle:

• The background is subtracted, and the barycentre position of an enclosed object is computed, and returned in out_xy_centre.
• If out_xy_centre_err is not NULL, the uncertainty of the centre position is computed. This is done by statistical analysis, and so the gain is required therefore.
• If out_xy_sigma is not NULL, the standard deviation is computed in a region of +/- 1.1*FWHM.
• If out_xy_fwhm is not NULL, the measured full-width-half-maximum is returned.

Bad Pixel Handling:

Bad pixel maps are supported. This means that bad pixels are ignored during the computation. But be aware that this may lead to a shift of the barycenter, so a proper way of interpolation should better be applied before centroiding.

Constraints:

• The peak must be positive, i.e. brighter than the background.
• image must be of size 3x3 or greater.

Algorithm

• The following process is iterated, until the window is stable:
  • The background and its noise are computed using iterative kappa-sigma clipping, using values inside +/- 3 sigma,
  • the barycentre is computed for values above background + 3*sigma.
  • in the case of failure, the window is enlarged,
  • the FWHM is measured,
  • and the window is refined to cover the +/- 2*FWHM region.
• If required, the PSF sigma out_xy_sigma and the centre uncertainty out_xy_centre_err are computed in the +/- 1.1*FWHM region.
• If required, the intensity of the pixel closest to the centre (and not flagged as bad) is returned.

Error Handling:

In the case of error, all output values are set to -1.0.
The following errors can be set and returned:

• CPL_ERROR_NULL_INPUT: if image or xy_centre is NULL
• CPL_ERROR_INVALID_TYPE: image is not of type CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE
• CPL_ERROR_ACCESS_OUT_OF_RANGE: window specifies coordinates outside the image plane
• CPL_ERROR_ILLEGAL_INPUT:
  • the image (window) size is not >= 3x3
• CPL_ERROR_DATA_NOT_FOUND:
  • less than 2 pixels are good, or
  • around the centre of the estimated peak are (too many) bad pixels
• CPL_ERROR_CONTINUE: the barycentre computation could not succeed.

Todo:

• CLIPM_CENTROIDING_MIN_WINDOWSIZE
• replace +/- 2 fwhm by _clipm_priv_image_get_optimal_box_width() (or better, use the output from clipm_priv_image_estimate_fwhm_xy())

clipm_centroiding_multi_gauss()

cpl_error_code clipm_centroiding_multi_gauss	(	const cpl_image *	image,
		const cpl_matrix *	locations,
		cpl_size	areasize,
		cpl_matrix **	xy_centre,
		cpl_matrix **	xy_centre_err,
		cpl_matrix **	xy_sigma,
		cpl_matrix **	xy_sigma_err,
		cpl_matrix **	xy_fwhm,
		cpl_matrix **	xy_fwhm_err,
		cpl_matrix **	centre_intensities,
		cpl_array **	all_error_codes,
		int	robustness
	)

Invoke clipm_centroiding_gauss() at several locations in an image.

Parameters

image	Input image (FITS convention)
locations	Center coordinates of sub-regions (windows) as 2 x N matrix
areasize	Edge-length of the square sub-regions (windows), at least 5
xy_centre	Output centre coordinate matrix of size 2 x N ([x, y]), can be NULL
xy_centre_err	Output centre error matrix (1 sigma) of size 2 x N, can be NULL
xy_sigma	Output sigma matrix of size 2 x N, can be NULL
xy_sigma_err	Output sigma error matrix (1 sigma) of size 2 x N, can be NULL
xy_fwhm	Output full-width-half-maximum matrix of size 2 x N, can be NULL
xy_fwhm_err	Output FWHM error matrix (1 sigma) of size 2 x N, can be NULL
centre_intensities	Output intensity matrix of the (closest) centre pixels of size 1 x N, can be NULL
all_error_codes	Output array of type CPL_TYPE_INT and size N, can be NULL
robustness	Maximum number of respective retries (>= 0)

Returns

CPL error code

Description:

• The function clipm_centroiding_gauss() is invoked N times, for N different square windows of image. The windows are defined by their centers, given by locations, which must be of the form:

  

• The edge-length of the square windows is given by areasize.
• all_error_codes returns all N returned error codes of clipm_centroiding_gauss().
• centre_intensities is returned as an 1 x N matrix containing the N intensities of the pixels closest to the N detected centres.
• All other output parameters are returned as 2 x N matrices, containing in each column the corresponding values, for x in row 0, and for y in row 1; For further details, please refer to clipm_centroiding_gauss().
• Any output parameter can be NULL.
• The window corner indices are rounded.

Error Handling:

The error handling is divided into two sections, the return value of this function, and the returned centroiding error states in all_error_codes.

1. The following error codes are set and returned by this function:
   • CPL_ERROR_NULL_INPUT: image or locations are NULL
   • CPL_ERROR_ILLEGAL_INPUT:
     • the number of rows in ref_locations is not 2
     • area_size < 5
     • robustness < 0
   • CPL_ERROR_INVALID_TYPE: image is not of type CPL_TYPE_INT, CPL_TYPE_FLOAT, or CPL_TYPE_DOUBLE
     
     
2. These error codes are stored in all_error_codes if either the horizontal or vertical centroiding failed:
   • 0: success
   • CPL_ERROR_ACCESS_OUT_OF_RANGE: the specified area exceeds the boundaries of the respective image
   • CPL_ERROR_ILLEGAL_OUTPUT: the determined centroid seems to be outside the window
   • CPL_ERROR_CONTINUE: the centroiding failed
   • CPL_ERROR_SINGULAR_MATRIX: the uncertainties could not be computed
   • -1: this value just catches potential library version conflicts