fitMoffat.hpp

Go to the documentation of this file.

/** \file fitMoffat.hpp
 * \author Jared R. Males
 * \brief Tools for fitting Moffat functions to data.
 * \ingroup fitting_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef fitMoffat_hpp
#define fitMoffat_hpp
 
#include "levmarInterface.hpp"
#include "../func/moffat.hpp"
 
#include "../../improc/eigenImage.hpp"
 
namespace mx
{
namespace math
{
namespace fit
{
 
/// Wrapper for a native array to pass to \ref levmarInterface, with Moffat details.
/** \ingroup moffat_peak_fit
 */
template <typename realT>
struct array2FitMoffat1D
{
    realT *data{ nullptr };   ///< Pointer to the array of y values
    realT *coords{ nullptr }; ///< Pointer to the array of x values (optional)
    size_t nx{ 0 };           ///< X dimension of the array
 
    realT m_I0{ 0 };
    realT m_I{ 0 };
    realT m_x0{ 0 };
    realT m_alpha{ 0 };
    realT m_beta{ 0 };
 
    int m_I0_idx{ 0 };
    int m_I_idx{ 1 };
    int m_x0_idx{ 2 };
    int m_alpha_idx{ 4 };
    int m_beta_idx{ 5 };
 
    int m_nparams = 5;
 
    /// Set whether each parameter is fixed.
    /** Sets the parameter indices appropriately.
     */
    void setFixed( bool I0, bool I, bool x0, bool alpha, bool beta )
    {
        int idx = 0;
 
        if( I0 )
            m_I0_idx = -1;
        else
            m_I0_idx = idx++;
 
        if( I )
            m_I_idx = -1;
        else
            m_I_idx = idx++;
 
        if( x0 )
            m_x0_idx = -1;
        else
            m_x0_idx = idx++;
 
        if( alpha )
            m_alpha_idx = -1;
        else
            m_alpha_idx = idx++;
 
        if( beta )
            m_beta_idx = -1;
        else
            m_beta_idx = idx++;
 
        m_nparams = idx;
    }
 
    realT I0( realT *p )
    {
        if( m_I0_idx < 0 )
        {
            return m_I0;
        }
        else
        {
            return p[m_I0_idx];
        }
    }
 
    void I0( realT *p, realT nI0 )
    {
        if( m_I0_idx < 0 )
        {
            m_I0 = nI0;
        }
        else
        {
            p[m_I0_idx] = nI0;
        }
    }
 
    realT I( realT *p )
    {
        if( m_I_idx < 0 )
        {
            return m_I;
        }
        else
        {
            return p[m_I_idx];
        }
    }
 
    void I( realT *p, realT nI )
    {
        if( m_I_idx < 0 )
        {
            m_I = nI;
        }
        else
        {
            p[m_I_idx] = nI;
        }
    }
 
    realT x0( realT *p )
    {
        if( m_x0_idx < 0 )
        {
            return m_x0;
        }
        else
        {
            return p[m_x0_idx];
        }
    }
 
    void x0( realT *p, realT nx0 )
    {
        if( m_x0_idx < 0 )
        {
            m_x0 = nx0;
        }
        else
        {
            p[m_x0_idx] = nx0;
        }
    }
 
    realT alpha( realT *p )
    {
        if( m_alpha_idx < 0 )
        {
            return m_alpha;
        }
        else
        {
            return p[m_alpha_idx];
        }
    }
 
    void alpha( realT *p, realT nalpha )
    {
        if( m_alpha_idx < 0 )
        {
            m_alpha = nalpha;
        }
        else
        {
            p[m_alpha_idx] = nalpha;
        }
    }
 
    realT beta( realT *p )
    {
        if( m_beta_idx < 0 )
        {
            return m_beta;
        }
        else
        {
            return p[m_beta_idx];
        }
    }
 
    void beta( realT *p, realT nbeta )
    {
        if( m_beta_idx < 0 )
        {
            m_beta = nbeta;
        }
        else
        {
            p[m_beta_idx] = nbeta;
        }
    }
 
    int nparams()
    {
        return m_nparams;
    }
};
 
/** \defgroup moffat_peak_fit Moffat Functions
  * \brief Fitting Moffat functions to data.
  *
  * The Moffat Function\cite moffat_1969, a.k.a. the Moffat Profile, a.k.a. the Moffat Distribution, has the form
  * \f[
    I(x) = I_{pk}\left[ 1 + \frac{x^2}{\alpha^2}\right]^{-\beta}
  * \f]
  * With \f$\beta=1\f$ it is the
  * Lorentzian or Cauchy distribution.  See also https://en.wikipedia.org/wiki/Moffat_distribution and
  * https://en.wikipedia.org/wiki/Cauchy_distribution.
  *
  * For utilities for working with Moffat functions see \ref gen_math_moffats.
  *
  * \ingroup peak_fit
  */
 
///\ref levmarInterface fitter structure for the 1D Moffat.
/** \ingroup moffat_peak_fit
 *
 */
template <typename _realT>
struct moffat1D_fitter
{
    typedef _realT realT;
 
    static const int nparams = 6;
 
    static void func( realT *p, realT *hx, int m, int n, void *adata )
    {
        array2FitMoffat1D<realT> *arr = (array2FitMoffat1D<realT> *)adata;
 
        size_t idx_mat, idx_dat;
 
        idx_dat = 0;
 
        realT I0 = arr->I0( p );
        realT I = arr->I( p );
        realT x0 = arr->x0( p );
        realT alpha = arr->alpha( p );
        realT beta = arr->beta( p );
 
        if( arr->coords )
        {
            for( int i = 0; i < arr->nx; ++i )
            {
                hx[i] = func::moffat<realT>( arr->coords[i], I0, I, x0, alpha, beta ) - arr->data[i];
            }
        }
        else
        {
            for( int i = 0; i < arr->nx; ++i )
            {
                hx[i] = func::moffat<realT>( i, I0, I, x0, alpha, beta ) - arr->data[i];
            }
        }
    }
};
 
/// Class to manage fitting a 1D Moffat to data via the \ref levmarInterface
/** Fits \ref gen_math_moffats to a 1-dimensional array of data.
 *
 * This class allows for treating any of the parameters as fixed.
 *
 * \tparam fitterT a type meeting the requirements specified in \ref levmarInterface.
 *
 * \ingroup moffat_peak_fit
 *
 */
template <typename _realT>
class fitMoffat1D : public levmarInterface<moffat1D_fitter<_realT>>
{
 
  public:
    typedef _realT realT;
    typedef moffat1D_fitter<realT> fitterT;
 
  protected:
    array2FitMoffat1D<realT> arr;
 
    void initialize()
    {
        this->allocate_params( arr.nparams() );
        this->adata = &arr;
    }
 
  public:
    fitMoffat1D()
    {
        initialize();
    }
 
    ~fitMoffat1D()
    {
    }
 
    /// Set whether each parameter is fixed.
    /** Sets the parameter indices appropriately.
     */
    void setFixed( bool I0,    ///< [in] if true, then I0 will be not be part of the fit
                   bool I,     ///< [in] if true, then I will be not be part of the fit
                   bool x0,    ///< [in] if true, then x0 will be not be part of the fit
                   bool alpha, ///< [in] if true, then alpha will be not be part of the fit
                   bool beta   ///< [in] if true, then beta will be not be part of the fit
    )
    {
        arr.setFixed( I0, I, x0, alpha, beta );
        this->allocate_params( arr.nparams() );
    }
 
    /// Set the initial guess for a symmetric Moffat.
    void setGuess( realT I0,    ///< [in] the constant background level
                   realT I,     ///< [in] the peak scaling
                   realT x0,    ///< [in] the center x-coordinate
                   realT alpha, ///< [in] the width parameter
                   realT beta   ///< [in] the shape parameter
    )
    {
        arr.I0( this->p, I0 );
        arr.I( this->p, I );
        arr.x0( this->p, x0 );
        arr.alpha( this->p, alpha );
        arr.beta( this->p, beta );
    }
 
    /// Set the data aray.
    void setArray( realT *data, ///< [in] pointer to an nx sized array of the y-values to be fit
                   int nx       ///< [in] the number of pixels in the x direction of the data array
    )
    {
        arr.data = data;
        arr.nx = nx;
 
        this->n = nx;
    }
 
    /// Set the data aray and the coordinates.
    void setArray( realT *data,   ///< [in] pointer to an nx sized array of the y-values to be fit
                   realT *coords, ///< [in] point to an nx sized array of the x-values of the data
                   int nx         ///< [in] the number of pixels in the x direction of the data array
    )
    {
        arr.data = data;
        arr.coords = coords;
        arr.nx = nx;
 
        this->n = nx;
    }
 
    /// Get the current value of I0, the constant.
    /**
     * \returns the current value of I0
     */
    realT I0()
    {
        return arr.I0( this->p );
    }
 
    /// Get the current value of I, the peak scaling.
    /**
     * \returns the current value of I
     */
    realT I()
    {
        return arr.I( this->p );
    }
 
    /// Get the center x-
    /**
     * \returns the current value of x0
     */
    realT x0()
    {
        return arr.x0( this->p );
    }
 
    /// Return the width parameter
    /**
     * \returns the current value of alpha
     */
    realT alpha()
    {
        return arr.alpha( this->p );
    }
 
    /// Return the shape parameter
    /**
     * \returns the current value of beta
     */
    realT beta()
    {
        return arr.beta( this->p );
    }
 
    /// Return the full-width at half maximum
    /**
     * \returns the FWHM calculated from alpha and beta
     */
    realT fwhm()
    {
        return func::moffatFWHM( alpha(), beta() );
    }
};
 
// forward
template <typename realT>
struct array2FitMoffat;
 
// forward
template <typename _realT>
struct moffat2D_sym_fitter;
 
/// Class to manage fitting a 2D Moffat to data via the \ref levmarInterface
/** Fits \ref gen_math_moffats to a 2-dimensional array of data.
 *
 * This class allows for treating any of the parameters as fixed.
 *
 * \tparam fitterT a type meeting the requirements specified in \ref levmarInterface.
 *
 * \ingroup moffat_peak_fit
 *
 */
template <typename fitterT>
class fitMoffat2D : public levmarInterface<fitterT>
{
 
  public:
    typedef typename fitterT::realT realT;
 
  protected:
    array2FitMoffat<realT> arr;
 
    void initialize()
    {
        this->allocate_params( arr.nparams() );
        this->adata = &arr;
    }
 
  public:
    fitMoffat2D()
    {
        initialize();
    }
 
    ~fitMoffat2D()
    {
    }
 
    /// Set whether each parameter is fixed.
    /** Sets the parameter indices appropriately.
     */
    void setFixed( bool I0,    ///< [in] if true, then I0 will be not be part of the fit
                   bool I,     ///< [in] if true, then I will be not be part of the fit
                   bool x0,    ///< [in] if true, then x0 will be not be part of the fit
                   bool y0,    ///< [in] if true, then y0 will be not be part of the fit
                   bool alpha, ///< [in] if true, then alpha will be not be part of the fit
                   bool beta   ///< [in] if true, then beta will be not be part of the fit
    )
    {
        arr.setFixed( I0, I, x0, y0, alpha, beta );
        this->allocate_params( arr.nparams() );
    }
 
    /// Set the initial guess for a symmetric Moffat.
    void setGuess( realT I0,    ///< [in] the constant background level
                   realT I,     ///< [in] the peak scaling
                   realT x0,    ///< [in] the center x-coordinate
                   realT y0,    ///< [in] the center y-coordinate
                   realT alpha, ///< [in] the width parameter
                   realT beta   ///< [in] the shape parameter
    )
    {
        arr.I0( this->p, I0 );
        arr.I( this->p, I );
        arr.x0( this->p, x0 );
        arr.y0( this->p, y0 );
        arr.alpha( this->p, alpha );
        arr.beta( this->p, beta );
    }
 
    /// Set the data aray.
    void setArray( realT *data, ///< [in] pointer to an nx X ny array of data to be fit
                   int nx,      ///< [in] the number of pixels in the x direction of the data array
                   int ny       ///< [in] the number of pixels in the y direction of the data array
    )
    {
        arr.data = data;
        arr.nx = nx;
        arr.ny = ny;
 
        this->n = nx * ny;
    }
 
    /// Get the current value of I0, the constant.
    /**
     * \returns the current value of I0
     */
    realT I0()
    {
        return arr.I0( this->p );
    }
 
    /// Get the current value of I, the peak scaling.
    /**
     * \returns the current value of I
     */
    realT I()
    {
        return arr.I( this->p );
    }
 
    /// Get the center x-
    /**
     * \returns the current value of x0
     */
    realT x0()
    {
        return arr.x0( this->p );
    }
 
    /// Get the center y-coordinate
    /**
     * \returns the current value of y0
     */
    realT y0()
    {
        return arr.y0( this->p );
    }
 
    /// Return the width parameter
    /**
     * \returns the current value of alpha
     */
    realT alpha()
    {
        return arr.alpha( this->p );
    }
 
    /// Return the shape parameter
    /**
     * \returns the current value of beta
     */
    realT beta()
    {
        return arr.beta( this->p );
    }
 
    /// Return the full-width at half maximum
    /**
     * \returns the FWHM calculated from alpha and beta
     */
    realT fwhm()
    {
        return func::moffatFWHM( alpha(), beta() );
    }
};
 
/// Wrapper for a native array to pass to \ref levmarInterface, with Moffat details.
/** \ingroup moffat_peak_fit
 */
template <typename realT>
struct array2FitMoffat
{
    realT *data{ nullptr }; ///< Pointer to the array
    size_t nx{ 0 };         ///< X dimension of the array
    size_t ny{ 0 };         ///< Y dimension of the array
 
    realT m_I0{ 0 };
    realT m_I{ 0 };
    realT m_x0{ 0 };
    realT m_y0{ 0 };
    realT m_alpha{ 0 };
    realT m_beta{ 0 };
 
    int m_I0_idx{ 0 };
    int m_I_idx{ 1 };
    int m_x0_idx{ 2 };
    int m_y0_idx{ 3 };
    int m_alpha_idx{ 4 };
    int m_beta_idx{ 5 };
 
    int m_nparams = 6;
 
    /// Set whether each parameter is fixed.
    /** Sets the parameter indices appropriately.
     */
    void setFixed( bool I0, bool I, bool x0, bool y0, bool alpha, bool beta )
    {
        int idx = 0;
 
        if( I0 )
            m_I0_idx = -1;
        else
            m_I0_idx = idx++;
 
        if( I )
            m_I_idx = -1;
        else
            m_I_idx = idx++;
 
        if( x0 )
            m_x0_idx = -1;
        else
            m_x0_idx = idx++;
 
        if( y0 )
            m_y0_idx = -1;
        else
            m_y0_idx = idx++;
 
        if( alpha )
            m_alpha_idx = -1;
        else
            m_alpha_idx = idx++;
 
        if( beta )
            m_beta_idx = -1;
        else
            m_beta_idx = idx++;
 
        m_nparams = idx;
    }
 
    realT I0( realT *p )
    {
        if( m_I0_idx < 0 )
        {
            return m_I0;
        }
        else
        {
            return p[m_I0_idx];
        }
    }
 
    void I0( realT *p, realT nI0 )
    {
        if( m_I0_idx < 0 )
        {
            m_I0 = nI0;
        }
        else
        {
            p[m_I0_idx] = nI0;
        }
    }
 
    realT I( realT *p )
    {
        if( m_I_idx < 0 )
        {
            return m_I;
        }
        else
        {
            return p[m_I_idx];
        }
    }
 
    void I( realT *p, realT nI )
    {
        if( m_I_idx < 0 )
        {
            m_I = nI;
        }
        else
        {
            p[m_I_idx] = nI;
        }
    }
 
    realT x0( realT *p )
    {
        if( m_x0_idx < 0 )
        {
            return m_x0;
        }
        else
        {
            return p[m_x0_idx];
        }
    }
 
    void x0( realT *p, realT nx0 )
    {
        if( m_x0_idx < 0 )
        {
            m_x0 = nx0;
        }
        else
        {
            p[m_x0_idx] = nx0;
        }
    }
 
    realT y0( realT *p )
    {
        if( m_y0_idx < 0 )
        {
            return m_y0;
        }
        else
        {
            return p[m_y0_idx];
        }
    }
 
    void y0( realT *p, realT ny0 )
    {
        if( m_y0_idx < 0 )
        {
            m_y0 = ny0;
        }
        else
        {
            p[m_y0_idx] = ny0;
        }
    }
 
    realT alpha( realT *p )
    {
        if( m_alpha_idx < 0 )
        {
            return m_alpha;
        }
        else
        {
            return p[m_alpha_idx];
        }
    }
 
    void alpha( realT *p, realT nalpha )
    {
        if( m_alpha_idx < 0 )
        {
            m_alpha = nalpha;
        }
        else
        {
            p[m_alpha_idx] = nalpha;
        }
    }
 
    realT beta( realT *p )
    {
        if( m_beta_idx < 0 )
        {
            return m_beta;
        }
        else
        {
            return p[m_beta_idx];
        }
    }
 
    void beta( realT *p, realT nbeta )
    {
        if( m_beta_idx < 0 )
        {
            m_beta = nbeta;
        }
        else
        {
            p[m_beta_idx] = nbeta;
        }
    }
 
    int nparams()
    {
        return m_nparams;
    }
};
 
///\ref levmarInterface fitter structure for the symmetric Moffat.
/** \ingroup moffat_peak_fit
 *
 */
template <typename _realT>
struct moffat2D_sym_fitter
{
    typedef _realT realT;
 
    static const int nparams = 6;
 
    static void func( realT *p, realT *hx, int m, int n, void *adata )
    {
        array2FitMoffat<realT> *arr = (array2FitMoffat<realT> *)adata;
 
        size_t idx_mat, idx_dat;
 
        idx_dat = 0;
 
        realT I0 = arr->I0( p );
        realT I = arr->I( p );
        realT x0 = arr->x0( p );
        realT y0 = arr->y0( p );
        realT alpha = arr->alpha( p );
        realT beta = arr->beta( p );
 
        for( int i = 0; i < arr->nx; ++i )
        {
            for( int j = 0; j < arr->ny; ++j )
            {
                idx_mat = i + j * arr->nx;
 
                hx[idx_dat] = func::moffat2D<realT>( i, j, I0, I, x0, y0, alpha, beta ) - arr->data[idx_mat];
 
                ++idx_dat;
            }
        }
    }
};
 
/// Alias for the fitMoffat2D type fitting the symmetric Moffat profile.
/** \ingroup moffat_peak_fit
 */
template <typename realT>
using fitMoffat2Dsym = mx::math::fit::fitMoffat2D<mx::math::fit::moffat2D_sym_fitter<realT>>;
 
} // namespace fit
} // namespace math
 
} // namespace mx
 
#endif // fitMoffat_hpp