imCenterRadon.hpp

Go to the documentation of this file.

/** \file imCenterRadon.hpp
 * \brief A class to find PSF centers using circular symmetry.
 * \ingroup image_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// 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 __imCenterRadon_hpp__
#define __imCenterRadon_hpp__
 
#include "../math/vectorUtils.hpp"
#include "../math/fit/fitGaussian.hpp"
#include "imageFilters.hpp"
#include "eigenImage.hpp"
#include "imagePeakInterp.hpp"
 
// #define ICCS_OMP
namespace mx
{
namespace improc
{
 
/// Find the center of the image of a point source using circular symmetry of the PSF.
/** Conducts a grid search through possible center positions.  Finds the point with the minimum
 * total variance around concentric circles about that point.  The cost grid is fit with a Gaussian.
 *
 * \tparam realT is a real floating point type.
 *
 * \ingroup image_reg
 */
template <typename transformT>
struct imCenterRadon
{
 
    typedef typename transformT::arithT realT;
 
    eigenImage<realT> m_grid; ///< The cost grid.
 
    realT m_maxPix{ 2 };    ///< The maximum range of the grid search, the grid will be +/- maxPix.  Default: 5.0.
    realT m_dPix{ 0.1 };    ///< The spacing of the points in the grid search.  Default: 0.1.
    realT m_peakTol{ 0.1 }; ///< The tolerance for the peak finding.  This is multiplicative with m_dPix.
 
    realT m_minRad{ 1.0 }; ///< The minimum radius to include in the transform.   Default: 1.0.
    realT m_maxRad{ 10 };  ///< The maximum radius to include in the transform.  Default: 10.0.
    realT m_dRad{ 0.1 };   ///< The spacing of the radii.   Default: 0.1.
 
    realT m_angTolPix{
        1.0 }; ///< The angle tolerance in pixels at the maximum radius, sets the angular step size. Default is 1.0.
 
    realT m_guessWidth{ 0.5 }; ///< Guess in original pixels for the FWHM of the cost grid.  Default is 0.5, which seems
                               ///< to work well in most cases.
 
    /// The fitter for finding the centroid of the grid.  You can use this to inspect the details of the fit if desired.
    // mx::math::fit::fitGaussian2D<mx::gaussian2D_gen_fitter<realT>> fit;
    mx::math::fit::fitGaussian2Dgen<realT> m_fit;
 
    realT m_x0; ///< Internal storage of the guess center position
    realT m_y0; ///< Internal storage of the guess center position.
 
    imCenterRadon()
    {
    }
 
    /// Peform the grid search and fit.
    int center( realT &x,              ///< [out] the x-coordinate of the estimated center of rotational symmetry
                realT &y,              ///< [out] the y-coordinate estimated center of rotational symmetry
                eigenImage<realT> &im, ///< [in] the image to centroid
                realT x0, ///< [in] an initial guess at the x-coordinate, the grid is centered about this point.
                realT y0  ///< [in] an initial guess at the y-coordinate, the grid is centered about this point.
    )
    {
        m_x0 = x0;
        m_y0 = y0;
 
        const int lbuff = transformT::lbuff;
        const int width = transformT::width;
 
        transformT trans;
 
        m_grid.resize( 2 * ( m_maxPix / m_dPix ) + 1, 2 * ( m_maxPix / m_dPix ) + 1 );
 
        int N = 2 * ( m_maxPix / m_dPix ) + 1;
 
        realT startX, startY;
 
        cubicConvolTransform<realT> cct;
 
        for( int i = 0; i < N; ++i )
        {
            startX = x0 - m_maxPix + i * m_dPix;
            for( int j = 0; j < N; ++j )
            {
                startY = y0 - m_maxPix + j * m_dPix;
 
                m_grid( i, j ) = 0;
 
                realT dang = math::rtod( atan( 1.0 / m_maxRad ) );
                int Nang = 360. / dang + 1;
                dang = 360. / Nang;
 
#ifdef MX_IMPROC_IMCENTERRADON_OMP
#pragma omp parallel for
#endif
                for( int iang = 0; iang < Nang; ++iang )
                {
                    realT ang = iang * dang;
 
                    eigenImage<realT> kern;
                    kern.resize( width, width );
 
                    realT cosang = cos( ang * 3.14159 / 180. );
                    realT sinang = sin( ang * 3.14159 / 180. );
                    realT lineint = 0;
                    for( realT rad = m_minRad; rad < m_maxRad; rad += m_dRad )
                    {
                        realT xx = startX + rad * cosang;
                        realT yy = startY + rad * sinang;
 
                        int i0 = xx;
                        int j0 = yy;
                        realT dx = xx - i0;
                        realT dy = yy - j0;
 
                        trans( kern, dx, dy );
                        lineint += ( im.block( i0 - lbuff, j0 - lbuff, width, width ) * kern ).sum();
                    }
#pragma omp critical
                    m_grid( i, j ) += lineint * lineint * pow( m_dRad, 2 );
                }
                m_grid( i, j ) /= Nang;
            } // j
        } // i
 
        return fitGrid( x, y, m_grid );
    }
 
    /// Fit the transform grid
    /** A separate function for testing, does not normally need to be called independently
     */
    int fitGrid( realT &x, ///< [out] the x-coordinate of the estimated center of rotational symmetry
                 realT &y, ///< [out] the y-coordinate estimated center of rotational symmetry
                 eigenImage<realT> grid )
    {
        imagePeakInterp<cubicConvolTransform<realT>> ipi( m_peakTol );
 
        ipi( x, y, grid );
 
        x = m_x0 - m_maxPix + x * m_dPix;
        y = m_y0 - m_maxPix + y * m_dPix;
 
        return 0;
 
        // Get initial guess for the fit.
        realT mx, mn, gx, gy;
 
        mx = grid.maxCoeff( &gx, &gy );
        mn = grid.minCoeff();
 
        //------------- Now fit the grid to a 2D elliptical Gaussian ------------
        /** \todo This needs to be able to handle highly non-symmetric peaks much better.
         */
        m_fit.setArray( grid.data(), grid.rows(), grid.cols() );
        m_fit.setGuess( mn, mx - mn, gx, gy, m_guessWidth / m_dPix, m_guessWidth / m_dPix, 0 );
 
        m_fit.fit();
 
        // The estimated enter of circular symmetry from the fit.
        x = m_x0 - m_maxPix + m_fit.x0() * m_dPix;
        y = m_y0 - m_maxPix + m_fit.y0() * m_dPix;
 
        return 0;
    }
 
    /// Output the results to a stream
    /** Prints a result summary to the input stream.
     *
     * \tparam iosT is a std::ostream-like type.
     * \tparam comment is a comment character to start each line.  Can be '\0'.
     */
    template <typename iosT, char comment = '#'>
    iosT &dumpResults( iosT &ios )
    {
        char c[] = { comment, '\0' };
 
        ios << c << "--------------------------------------\n";
        ios << c << "mx::improc::imCenterRadon Results \n";
        ios << c << "--------------------------------------\n";
        ios << c << "Estimated x-center: " << m_x0 - m_maxPix + m_fit.x0() * m_dPix << "\n";
        ios << c << "Estimated y-center: " << m_y0 - m_maxPix + m_fit.y0() * m_dPix << "\n";
        ios << c << "Cost half-width: "
            << 0.5 * math::func::sigma2fwhm<realT>( sqrt( pow( m_fit.sigma_x(), 2 ) + pow( m_fit.sigma_y(), 2 ) ) *
                                                    m_dPix )
            << " pixels\n";
        ios << c << "--------------------------------------\n";
        ios << c << "Setup:\n";
        ios << c << "--------------------------------------\n";
        ios << c << "maxPix: " << m_maxPix << "\n";
        ios << c << "dPix: " << m_dPix << "\n";
        ios << c << "minRad: " << m_minRad << "\n";
        ios << c << "maxRad: " << m_maxRad << "\n";
        ios << c << "dRad: " << m_dRad << "\n";
        ios << c << "angTolPix: " << m_angTolPix << "\n";
        ios << c << "guessWidth: " << m_guessWidth << "\n";
        ios << c << "--------------------------------------\n";
        ios << c << "Fit results:\n";
        m_fit.dumpReport( ios );
        ios << c << "--------------------------------------\n";
 
        return ios;
    }
 
    /// Output the results to std::cout
    /** Prints a result summary to std::cout
     *
     * \tparam comment is a comment character to start each line.  Can be '\0'.
     */
    template <char comment = '#'>
    std::ostream &dumpResults()
    {
        return dumpResults<std::ostream, comment>( std::cout );
    }
};
 
} // namespace improc
} // namespace mx
 
#endif //__imCenterRadon_hpp__