fitAiry.hpp

Go to the documentation of this file.

/** \file fitAiry.hpp
 * \author Jared R. Males
 * \brief Tools for fitting the Airy pattern to PSF images.
 * \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 fitAiry_hpp
#define fitAiry_hpp
 
#include "levmarInterface.hpp"
#include "../func/airyPattern.hpp"
 
namespace mx
{
namespace math
{
namespace fit
{
 
template <typename realT>
struct array2FitAiry;
 
template <typename _realT>
struct airy2D_obs_fitter;
 
template <typename _realT>
struct airy2D_obs_fitter_ps;
 
/** \defgroup airy_peak_fit Airy Patterns
 * \brief Fitting Airy patterns to data.
 *
 * The Airy pattern is fit to data, typically diffraction-limited images of point sources.
 *
 * \ingroup peak_fit
 */
 
/// Class to manage fitting a 2D Airy pattern to data via the \ref levmarInterface
/** In addition to the requirements on fitterT specified by \ref levmarInterface
 * this class also requires this definition in fitterT
 * \code
 * static const int nparams = 5;
 * \endcode
 * where the number 5 is replaced by the number of parameters that fitterT expects to fit.
 *
 * \tparam fitterT a type meeting the above requirements.
 *
 * \ingroup airy_peak_fit
 *
 */
template <typename fitterT>
class fitAiry2D : public levmarInterface<fitterT>
{
 
  public:
    typedef typename fitterT::realT realT;
 
    static const int nparams = fitterT::nparams;
 
    array2FitAiry<realT> arr;
 
    void initialize()
    {
        this->allocate_params( nparams );
        this->adata = &arr;
    }
 
    fitAiry2D()
    {
        initialize();
    }
 
    ~fitAiry2D()
    {
    }
 
    /// Set the initial guess when platescale and central obscuration are fixed.
    void setGuess( realT A0, ///< [in] the constant background
                   realT A,  ///< [in] the peak scaling
                   realT x0, ///< [in] the center x-coordinate [pixels]
                   realT y0  ///< [in] the center y-coordinate [pixels]
    )
    {
        static_assert( nparams == 4, "fitAiry2D: Wrong setGuess called for platescale and/or cen-obs variable." );
 
        this->p[0] = A0;
        this->p[1] = A;
        this->p[2] = x0;
        this->p[3] = y0;
    }
 
    /// Set the initial guess when platescale is variable, and central obscuration is fixed.
    void setGuess( realT A0, ///< [in] the constant background
                   realT A,  ///< [in] the peak scaling
                   realT x0, ///< [in] the center x-coordinate [pixels]
                   realT y0, ///< [in] the center y-coordinate [pixels]
                   realT ps  ///< [in] the platescale [ (lambda/D) / pixel]
    )
    {
        static_assert( nparams == 5, "fitAiry2D: Wrong setGuess called for only platescale variable." );
 
        this->p[0] = A0;
        this->p[1] = A;
        this->p[2] = x0;
        this->p[3] = y0;
        this->p[4] = ps;
    }
 
    /// Set the initial guess when central-obscuration is variable.
    void setGuess( realT A0, ///< [in] the constant background
                   realT A,  ///< [in] the peak scaling
                   realT x0, ///< [in] the center x-coordinate [pixels]
                   realT y0, ///< [in] the center y-coordinate [pixels]
                   realT ps, ///< [in] the platescale [ (lambda/D) / pixel]
                   realT co  ///< [in] the central obscuration [ratio]
    )
    {
        static_assert( nparams == 6, "fitAiry2D: Wrong setGuess called for cen-obs not variable." );
 
        this->p[0] = A0;
        this->p[1] = A;
        this->p[2] = x0;
        this->p[3] = y0;
        this->p[4] = ps;
        this->p[5] = co;
    }
 
    void setArray( realT *data, int nx, int ny )
    {
        arr.data = data;
        arr.nx = nx;
        arr.ny = ny;
 
        this->n = nx * ny;
    }
 
    void cenObs( realT co )
    {
        arr.cenObs = co;
    }
 
    void ps( realT ps )
    {
        arr.ps = ps;
    }
 
    int fit()
    {
        return levmarInterface<fitterT>::fit();
    }
 
    realT A0()
    {
        return this->p[0];
    }
 
    realT A()
    {
        return this->p[1];
    }
 
    realT x0()
    {
        return this->p[2];
    }
 
    realT y0()
    {
        return this->p[3];
    }
 
    realT ps()
    {
        if( nparams < 5 )
        {
            return arr.ps;
        }
        else
        {
            return this->p[4];
        }
    }
 
    realT cenObs()
    {
        if( nparams < 6 )
        {
            return arr.cenObs;
        }
        else
        {
            return this->p[5];
        }
    }
};
 
/// Wrapper for a native array to pass to \ref levmarInterface, with Airy details.
/** \ingroup airy_peak_fit
 */
template <typename realT>
struct array2FitAiry
{
    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 cenObs{ 0 }; ///< is the ratio of the circular central obscuration diameter to the diameter.
    realT ps{ 0 };     ///< the platescale in \f$ (\lambda/D)/pixel  \f$
};
 
///\ref levmarInterface fitter structure for the centrally obscured Airy pattern.
/**
 * Platescale and central obscuration are fixed.
 *
 * \ingroup airy_peak_fit
 *
 */
template <typename _realT>
struct airy2D_obs_fitter
{
    typedef _realT realT;
 
    static const int nparams = 4;
 
    static void func( realT *p, realT *hx, int m, int n, void *adata )
    {
        array2FitAiry<realT> *arr = (array2FitAiry<realT> *)adata;
 
        size_t idx_mat, idx_dat;
 
        idx_dat = 0;
 
        /*
          p[0] = floor
          p[1] = scale
          p[2] = x0
          p[3] = y0
        */
 
        realT r;
 
        for( int i = 0; i < arr->nx; i++ )
        {
            for( int j = 0; j < arr->ny; j++ )
            {
                idx_mat = i + j * arr->nx;
 
                r = sqrt( pow( i - p[2], 2 ) + pow( j - p[3], 2 ) );
 
                hx[idx_dat] = func::airyPattern( static_cast<realT>( i ),
                                                 static_cast<realT>( j ),
                                                 p[0],
                                                 p[1],
                                                 p[2],
                                                 p[3],
                                                 arr->ps,
                                                 arr->cenObs ) -
                              arr->data[idx_mat];
 
                // hx[idx_dat] *= fabs(arr->data[idx_mat]);
 
                idx_dat++;
            }
        }
    }
};
 
///\ref levmarInterface fitter structure for the obstructed Airy pattern, including platescale.
/** Central obscuration is fixed.
 *
 * \ingroup airy_peak_fit
 *
 */
template <typename _realT>
struct airy2D_obs_fitter_ps
{
    typedef _realT realT;
 
    static const int nparams = 5;
 
    static void func( realT *p, realT *hx, int m, int n, void *adata )
    {
        array2FitAiry<realT> *arr = (array2FitAiry<realT> *)adata;
 
        size_t idx_mat, idx_dat;
 
        idx_dat = 0;
 
        /*
          p[0] = floor
          p[1] = scale
          p[2] = x0
          p[3] = y0
          p[4] = ps [(lam/D)/pix]
        */
 
        realT r;
 
        for( int i = 0; i < arr->nx; i++ )
        {
            for( int j = 0; j < arr->ny; j++ )
            {
                idx_mat = i + j * arr->nx;
 
                r = sqrt( pow( i - p[2], 2 ) + pow( j - p[3], 2 ) );
 
                hx[idx_dat] =
                    func::airyPattern(
                        static_cast<realT>( i ), static_cast<realT>( j ), p[0], p[1], p[2], p[3], p[4], arr->cenObs ) -
                    arr->data[idx_mat];
 
                // hx[idx_dat] *= fabs( pow(arr->data[idx_mat],2));
 
                idx_dat++;
            }
        }
    }
};
 
///\ref levmarInterface fitter structure for the obstructed Airy pattern, including fitting platescale and central
///obscuration.
/** \ingroup airy_peak_fit
 *
 */
template <typename _realT>
struct airy2D_obs_fitter_ps_eps
{
    typedef _realT realT;
 
    static const int nparams = 6;
 
    static void func( realT *p, realT *hx, int m, int n, void *adata )
    {
        array2FitAiry<realT> *arr = (array2FitAiry<realT> *)adata;
 
        size_t idx_mat, idx_dat;
 
        idx_dat = 0;
 
        /*
          p[0] = floor
          p[1] = scale
          p[2] = x0
          p[3] = y0
          p[4] = ps [(lam/D)/pix]
          p[5] = cen-obs
        */
 
        // realT r;
 
        for( int i = 0; i < arr->nx; ++i )
        {
            for( int j = 0; j < arr->ny; ++j )
            {
                idx_mat = i + j * arr->nx;
 
                // r = sqrt( pow( i-p[2],2) + pow(j-p[3],2));
 
                hx[idx_dat] =
                    func::airyPattern(
                        static_cast<realT>( i ), static_cast<realT>( j ), p[0], p[1], p[2], p[3], p[4], p[5] ) -
                    arr->data[idx_mat];
 
                // hx[idx_dat] *= fabs( pow(arr->data[idx_mat],2));
 
                idx_dat++;
            }
        }
    }
};
 
/// Alias for the fitAiry2D type with both platescale and cen-obs fixed.
/** \ingroup airy_peak_fit
 */
template <typename realT>
using fitAiry2DbothFixed = mx::math::fit::fitAiry2D<mx::math::fit::airy2D_obs_fitter<realT>>;
 
/// Alias for the fitAiry2D type with cen-obs fixed.
/** \ingroup airy_peak_fit
 */
template <typename realT>
using fitAiry2DcenObsFixed = mx::math::fit::fitAiry2D<mx::math::fit::airy2D_obs_fitter_ps<realT>>;
 
/// Alias for the fitAiry2D type with none fixed.
/** \ingroup airy_peak_fit
 */
template <typename realT>
using fitAiry2DnoneFixed = mx::math::fit::fitAiry2D<mx::math::fit::airy2D_obs_fitter_ps_eps<realT>>;
 
} // namespace fit
} // namespace math
} // namespace mx
 
#endif // fitAiry_hpp