Fit spectral components to a casacore::Vector of data from an image. More...

#include <ImageFit1D.h>

Public Types
enum	AbcissaType { PIXEL, IM_NATIVE, VELOCITY, N_TYPES }

using	FitterType = typename casacore::NumericTraits< T >::PrecisionType

Public Member Functions
	ImageFit1D (std::shared_ptr< const casacore::ImageInterface< T > > image, casacore::uInt axis=0)
	Constructor. More...

	ImageFit1D (std::shared_ptr< const casacore::ImageInterface< T > > image, std::shared_ptr< const casacore::ImageInterface< T > > weights, casacore::uInt axis=0)
	Constructor with fitting weights image. More...

	~ImageFit1D ()
	Destructor. More...

	ImageFit1D (const ImageFit1D &other)
	Copy constructor. More...

ImageFit1D &	operator= (const ImageFit1D &other)
	Assignment operator. More...

void	setData (const casacore::IPosition &pos, casacore::Array< FitterType >(*yfunc)(const casacore::Array< FitterType > &)=0)
	Set the data to be fit. More...

void	setData (const casacore::IPosition &pos, const ImageFit1D< T >::AbcissaType type, const casacore::Bool doAbs=true, const casacore::Double const &abscissaDivisor=0, casacore::Array< casacore::Double >(xfunc)(const casacore::Array< casacore::Double > &)=0, casacore::Array< FitterType >(*yfunc)(const casacore::Array< FitterType > &)=0)

void	setElements (const SpectralList &list)
	Set a SpectralList of SpectralElements to fit for. More...

void	addElement (const SpectralElement &el)
	Add new SpectralElement(s) to the SpectralList (can be empty) of SpectralElements to be fit for. More...

void	addElements (const SpectralList &list)

void	setGaussianElements (casacore::uInt nGauss)
	Set a SpectralList of Gaussian SpectralElements to fit for. More...

void	clearList ()
	Clear the SpectralList of elements to be fit for. More...

casacore::Bool	fit ()
	Do the fit and return convergence status. More...

casacore::Double	getChiSquared () const
	Get Chi Squared of fit. More...

casacore::Double	getNumberIterations () const
	Get number of iterations for last fit. More...

const SpectralList &	getList (casacore::Bool fit=true) const
	Recover the list of elements. More...

casacore::Vector< T >	getEstimate (casacore::Int which=-1) const
	Recover vectors for the estimate, fit and residual. More...

casacore::Vector< T >	getFit (casacore::Int which=-1) const

casacore::Vector< T >	getResidual (casacore::Int which=-1, casacore::Bool fit=true) const

casacore::Bool	setXMask (const std::set< casacore::uInt > &indices, casacore::Bool specifiedPixelsAreGood)

casacore::Vector< casacore::Bool >	getDataMask () const
	get data mask More...

casacore::Vector< casacore::Bool >	getTotalMask () const
	Get Total Mask (data and range mask) More...

casacore::Bool	succeeded () const
	did the fit succeed? should only be called after fit(). More...

casacore::Bool	converged () const
	did the fit converge? should only be called after fit(). More...

void	invalidate ()
	flag the solution as invalid based on external criteria. More...

casacore::Bool	isValid () const
	is the solution valid? If false, some external logic has called invalidate() More...

void	setAbscissa (const casacore::Vector< casacore::Double > &x)
	Set the abscissa values prior to running setData. More...

casacore::Vector < casacore::Double >	makeAbscissa (ImageFit1D< T >::AbcissaType type, casacore::Bool doAbs, const casacore::Double *const &abscissaDivisor)
	make the abscissa values, `x`. More...

Static Public Member Functions
static casacore::Bool	setAbcissaState (casacore::String &errMsg, ImageFit1D< T >::AbcissaType &type, casacore::CoordinateSystem &cSys, const casacore::String &xUnit, const casacore::String &doppler, casacore::uInt pixelAxis)
	Helper function. More...

Private Member Functions
	ImageFit1D ()
	Disallow default constructor. More...

void	check () const

void	checkType () const

void	_construct ()

void	copy (const ImageFit1D< T > &other)

void	_resetFitter ()
	void setWeightsImage (const casacore::ImageInterface<T>& im); More...

Private Attributes
std::shared_ptr< const casacore::ImageInterface< T > >	_image

std::shared_ptr< const casacore::ImageInterface< T > >	_weights

casacore::uInt	_axis

ProfileFit1D< FitterType >	_fitter
	In the future I will be able to template the fitter on T. More...

casacore::Bool	_converged

casacore::Bool	_success

casacore::Bool	_isValid

casacore::Vector < casacore::Double >	_x

casacore::Vector < casacore::Double >	_unityWeights

casacore::Vector < casacore::Double >	_weightSlice

casacore::IPosition	_sliceShape

Detailed Description

template<class T>
class casa::ImageFit1D< T >

Fit spectral components to a casacore::Vector of data from an image.

Intended use:

Public interface

Review Status

Test programs:: tImageFit1D

Prerequisite

Synopsis

Fit lists (held in class SpectralList) of SpectralElements to a casacore::Vector of data from the image. Each SpectralElement can be one from a variety of types. The values of the parameters for each SpectralElement provide the initial starting guesses for the fitting process.

You specify the domain in which the fit is to be done via the enum AbcissaType. The casacore::CoordinateSystem in the image is used to convert the pixel coordinates to the desired abcissa. You can change the units of the casacore::CoordinateSystem if you want to fit in different units. If you set an estimate yourself (function setElements or addElement) it is the callers responsibility that the elements are in the correct abcissa domain. Function setGaussianElements will automatically make an estimate in the correct domain.

Also, a SpectralElement object holds a mask indicating whether a parameter should be held fixed or solved for. After the fitting is done, a new SpectralList holding SpectralElements with the fitted parameters is created.

For all the functions that return a status casacore::Bool, true is good. If false is returned, an error message can be recovered with function errorMessage, You should not proceed if false is returned.

Exceptions will be thrown if you do not set the Image and axis via the constructor or setImage function.

Example

casacore::PagedImage<casacore::Float> im("myimage");
casacore::Int axis = 2;
ImageFit1D<casacore::Float> fitter(image, axis);
casacore::IPosition pos(in.ndim(),0);
fitter.setData(pos, ImageFit1D<casacore::Float>::IM_NATIVE);     // Fit in native coordinate space
fitter.setGaussianElements(3);                      // FIt 3 Gaussians
if (fitter.fit()) {
cerr << fitter.getList() << endl;                // Print result
}

To Do

Add constraints

Definition at line 121 of file ImageFit1D.h.

Member Typedef Documentation

template<class T>

using casa::ImageFit1D< T >::FitterType = typename casacore::NumericTraits<T>::PrecisionType

Definition at line 124 of file ImageFit1D.h.

Member Enumeration Documentation

template<class T>

enum casa::ImageFit1D::AbcissaType

Enumerator
PIXEL
IM_NATIVE
VELOCITY
N_TYPES

Definition at line 126 of file ImageFit1D.h.

Constructor & Destructor Documentation

template<class T>

casa::ImageFit1D< T >::ImageFit1D	(	std::shared_ptr< const casacore::ImageInterface< T > >	image,
		casacore::uInt	axis = `0`
	)

Constructor.

Fitting weights are assumed all unity.

template<class T>

casa::ImageFit1D< T >::ImageFit1D	(	std::shared_ptr< const casacore::ImageInterface< T > >	image,
		std::shared_ptr< const casacore::ImageInterface< T > >	weights,
		casacore::uInt	axis = `0`
	)

Constructor with fitting weights image.

The data and weights images must be the same shape.

template<class T>

casa::ImageFit1D< T >::~ImageFit1D ( )

Destructor.

template<class T>

casa::ImageFit1D< T >::ImageFit1D ( const ImageFit1D< T > & other )

Copy constructor.

Uses reference semantics.

template<class T>

casa::ImageFit1D< T >::ImageFit1D ( )

inlineprivate

Disallow default constructor.

Definition at line 308 of file ImageFit1D.h.

Member Function Documentation

template<class T>

void casa::ImageFit1D< T >::_construct ( )

private

template<class T>

void casa::ImageFit1D< T >::_resetFitter ( )

private

void setWeightsImage (const casacore::ImageInterface<T>& im);

reset the fitter, for example if we've done a fit and want to move to the next position in the image

template<class T>

void casa::ImageFit1D< T >::addElement ( const SpectralElement & el )

inline

Add new SpectralElement(s) to the SpectralList (can be empty) of SpectralElements to be fit for.

You must have already called setData to call this function.

Definition at line 208 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::addElement().

template<class T>

void casa::ImageFit1D< T >::addElements ( const SpectralList & list )

inline

Definition at line 209 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::addElements().

template<class T>

void casa::ImageFit1D< T >::check ( ) const

private

template<class T>

void casa::ImageFit1D< T >::checkType ( ) const

private

template<class T>

void casa::ImageFit1D< T >::clearList ( )

inline

Clear the SpectralList of elements to be fit for.

Definition at line 221 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::clearList().

template<class T>

casacore::Bool casa::ImageFit1D< T >::converged ( ) const

did the fit converge? should only be called after fit().

template<class T>

void casa::ImageFit1D< T >::copy ( const ImageFit1D< T > & other )

private

template<class T>

casacore::Bool casa::ImageFit1D< T >::fit ( )

Do the fit and return convergence status.

Errors in the fitting process will generate an casacore::AipsError exception and you should catch these yourself.

Referenced by casa::ImageFit1D< T >::getList().

template<class T>

casacore::Double casa::ImageFit1D< T >::getChiSquared ( ) const

inline

Get Chi Squared of fit.

Definition at line 229 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::getChiSquared().

template<class T>

casacore::Vector<casacore::Bool> casa::ImageFit1D< T >::getDataMask ( ) const

inline

get data mask

Definition at line 256 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::getDataMask().

template<class T>

casacore::Vector<T> casa::ImageFit1D< T >::getEstimate ( casacore::Int which = -1 ) const

Recover vectors for the estimate, fit and residual.

If you don't specify which element, all elements are included If the Vectors are returned with zero length, it means an error condition exists (e.g. asking for fit before you do one). In this case an error message can be recovered with function errorMessage.

template<class T>

casacore::Vector<T> casa::ImageFit1D< T >::getFit ( casacore::Int which = -1 ) const

template<class T>

const SpectralList& casa::ImageFit1D< T >::getList ( casacore::Bool fit = true ) const

inline

Recover the list of elements.

You can get the elements as initially estimated (fit=false), or after fitting (fit=true). In the latter case, the SpectralElements hold the parameters and errors of the fit.

Definition at line 238 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, casa::ImageFit1D< T >::fit(), and casa::ProfileFit1D< T >::getList().

template<class T>

casacore::Double casa::ImageFit1D< T >::getNumberIterations ( ) const

inline

Get number of iterations for last fit.

Definition at line 232 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::getNumberIterations().

template<class T>

casacore::Vector<T> casa::ImageFit1D< T >::getResidual	(	casacore::Int	which = `-1`,
		casacore::Bool	fit = `true`
	)		const

template<class T>

casacore::Vector<casacore::Bool> casa::ImageFit1D< T >::getTotalMask ( ) const

inline

Get Total Mask (data and range mask)

Definition at line 259 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::getTotalMask().

template<class T>

void casa::ImageFit1D< T >::invalidate ( )

flag the solution as invalid based on external criteria.

template<class T>

casacore::Bool casa::ImageFit1D< T >::isValid ( ) const

is the solution valid? If false, some external logic has called invalidate()

template<class T>

casacore::Vector<casacore::Double> casa::ImageFit1D< T >::makeAbscissa	(	ImageFit1D< T >::AbcissaType	type,
		casacore::Bool	doAbs,
		const casacore::Double *const &	abscissaDivisor
	)

make the abscissa values, x.

If type=IN_NATIVE and abscissaDivisor is not null, then divide the native values by the value pointed to by <src>abscissaDivisor in making the abscissa values.

template<class T>

ImageFit1D& casa::ImageFit1D< T >::operator= ( const ImageFit1D< T > & other )

Assignment operator.

Uses reference semantics.

template<class T>

static casacore::Bool casa::ImageFit1D< T >::setAbcissaState	(	casacore::String &	errMsg,
		ImageFit1D< T >::AbcissaType &	type,
		casacore::CoordinateSystem &	cSys,
		const casacore::String &	xUnit,
		const casacore::String &	doppler,
		casacore::uInt	pixelAxis
	)

static

Helper function.

Sets up the casacore::CoordinateSystem to reflect the choice of abcissa unit and the doppler (if the axis is spectral).

template<class T>

void casa::ImageFit1D< T >::setAbscissa ( const casacore::Vector< casacore::Double > & x )

inline

Set the abscissa values prior to running setData.

If this is done, then the abscissa values will not be recomputed when setData is called. This can imporove performance if, for example, you are looping over several fitters for which you know the abscissa values do not change.

Definition at line 285 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_x, and casacore::Vector< T >::assign().

template<class T>

void casa::ImageFit1D< T >::setData	(	const casacore::IPosition &	pos,
		casacore::Array< FitterType >(*)(const casacore::Array< FitterType > &)	yfunc = `0`
	)

Set the data to be fit.

All non-profile axes data are averaged. For the profile axis, the full spectrum is taken. The abscissa world values are computed when you call these functions unless they have been set previously by a call to setAbscissa() in which case the values that were passed to that method are used. Use the first form of setData() in this case. The domain of the abscissa values is controlled by AbcissaType and doAbs (absolute coordinates). The casacore::CoordinateSystem in the image is used to convert from pixels to world values. If type=IN_NATIVE and abscissaDivisor is not null, the world abscissa values will be divided by the value pointed to by abscissaDivisor. This mitigates having very large or very small abscissa values when fitting. If xfunc and/or yfunc is not NULL, the x and/or y values are fed to the specified function and the resultant values are what are used for the x and/or y values in the fit. If xfunc is not NULL and setAbscissa values has been called prior, no abscissa value transformation occurs. Thus if you want to apply a function to the abscissa values, the caller should pass the result of that function into setAbscissaValues.

template<class T>

void casa::ImageFit1D< T >::setData	(	const casacore::IPosition &	pos,
		const ImageFit1D< T >::AbcissaType	type,
		const casacore::Bool	doAbs = `true`,
		const casacore::Double *const &	abscissaDivisor = `0`,
		casacore::Array< casacore::Double >(*)(const casacore::Array< casacore::Double > &)	xfunc = `0`,
		casacore::Array< FitterType >(*)(const casacore::Array< FitterType > &)	yfunc = `0`
	)

template<class T>

void casa::ImageFit1D< T >::setElements ( const SpectralList & list )

inline

Set a SpectralList of SpectralElements to fit for.

These elements

must be in the correct abcissa domain set in function setData. You must have already called setData to call this function. The SpectralElements in the list hold the initial estimates. They also contain the information about whether specific parameters are to be held fixed or allowed to vary in the fitting process. You can recover the list of elements with function getList.

Definition at line 202 of file ImageFit1D.h.

References casa::ImageFit1D< T >::_fitter, and casa::ProfileFit1D< T >::setElements().

template<class T>

void casa::ImageFit1D< T >::setGaussianElements ( casacore::uInt nGauss )

Set a SpectralList of Gaussian SpectralElements to fit for.

The initial estimates for the Gaussians will be automatically determined in the correct abcissa domain. All of the parameters created by this function will be solved for by default. You can recover the list of elements with function getList. Status is returned, if false, error message can be recovered with errorMessage

template<class T>