imCenterCircleSym.hpp

Go to the documentation of this file.

/** \file imCenterCircleSym.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 __imCenterCircleSym_hpp__
#define __imCenterCircleSym_hpp__
 
#include "../math/vectorUtils.hpp"
#include "../math/fit/fitGaussian.hpp"
#include "imageFilters.hpp"
#include "eigenImage.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 realT>
struct imCenterCircleSym
{
 
    eigenImage<realT> *mask; ///< An optional pointer to a mask.  If not NULL and the same size as the input image, any
                             ///< 0 pixels in this image are ignored.
 
    eigenImage<realT> grid;   ///< The cost grid.
    eigenImage<realT> smGrid; ///< The smoothed cost grid.
 
    realT maxPix; ///< The maximum range of the grid search, the grid will be +/- maxPix.  Default: 5.0.
    realT dPix;   ///< The spacing of the points in the grid search.  Default: 0.1.
 
    realT minRad; ///< The minimum radius to include in the cost function.   Default: 1.0.
    realT maxRad; ///< The maximum radius to include in the cost function  Default: 20.0.
    realT dRad;   ///< The spacing of the radii.   Default: 1.0.
 
    realT smWidth; ///< FWHM of the Guassian smoothing kernel.  Default is 10 (seems to work well regardless of dPix).
                   ///< Set to 0 for no smoothing.
    realT guessWidth; ///< 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> fit;
 
    realT _x0; ///< Internal storage of the guess center position
    realT _y0; ///< Internal storage of the guess center position.
 
    imCenterCircleSym()
    {
 
        mask = 0;
 
        init( 5, 0.1, 1.0, 20.0, 1.0 );
 
        smWidth = 10;
        guessWidth = 0.5;
    }
 
    /// Initialize the grid.
    int init( realT maxP, ///< [in] the new value of maxPix
              realT dP,   ///< [in] the new value of dPix
              realT minR, ///< [in] the new value of minRad
              realT maxR, ///< [in] the new value of maxRad
              realT dR    ///< [in] the new value of dRad
    )
    {
        maxPix = maxP;
        dPix = dP;
        minRad = minR;
        maxRad = maxR;
        dRad = dR;
        return 0;
    }
 
    /// 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.
    )
    {
        _x0 = x0;
        _y0 = y0;
 
        bool doMask = false;
        if( mask )
        {
            if( mask->rows() == im.rows() && mask->cols() == im.cols() )
                doMask = true;
        }
 
        grid.resize( 2 * ( maxPix / dPix ) + 1, 2 * ( maxPix / dPix ) + 1 );
 
#ifdef ICCS_OMP
#pragma omp parallel
        {
#endif
            int nRad = ( maxRad - minRad ) / dRad + 1;
            std::vector<std::vector<realT>> rads;
            rads.resize( nRad );
 
            int N = 2 * ( maxPix / dPix ) + 1;
 
            realT startX, startY;
#ifdef ICCS_OMP
#pragma omp for
#endif
            for( int i = 0; i < N; ++i )
            {
                startX = x0 - maxPix + i * dPix;
                for( int j = 0; j < N; ++j )
                {
                    startY = y0 - maxPix + j * dPix;
 
                    for( int k = 0; k < rads.size(); ++k )
                    {
                        rads[k].clear();
                    }
 
                    for( int ii = startX - maxRad; ii <= startX + maxRad; ++ii )
                    {
                        for( int jj = startY - maxRad; jj <= startY + maxRad; ++jj )
                        {
 
                            if( doMask )
                            {
                                if( ( *mask )( ii, jj ) == 0 )
                                    continue;
                            }
 
                            realT rad;
                            rad = sqrt( pow( ii - startX, 2 ) + pow( jj - startY, 2 ) );
 
                            if( rad < minRad || rad >= maxRad )
                                continue;
 
                            rads[( rad - minRad ) / dRad].push_back( im( ii, jj ) );
                        } // jj
                    } // ii
 
                    grid( i, j ) = 0;
 
                    realT mean, var;
                    for( int k = 0; k < rads.size(); ++k )
                    {
                        if( rads[k].size() <= 1 )
                            continue;
                        mean = math::vectorMean( rads[k] );
                        var = math::vectorVariance( rads[k], mean );
 
                        grid( i, j ) += var / pow( mean, 2 );
                    }
                } // j
            } // i
#ifdef ICCS_OMP
        } // omp parallel
#endif
        //-------- Now Smooth the grid, and multiply by -1 for fitting -------------
 
        if( smWidth > 0 )
        {
            filterImage( smGrid, grid, gaussKernel<eigenImage<realT>, 2>( 10 ), 0 );
        }
        else
        {
            smGrid = grid;
        }
 
        smGrid *= -1;
 
        // Get initial guess for the fit.
        realT mx, mn, gx, gy;
 
        mx = smGrid.maxCoeff( &gx, &gy );
        mn = smGrid.minCoeff();
 
        //------------- Now fit the smoothed cost grid to a 2D elliptical Gaussian ------------
        /** \todo This needs to be able to handle highly non-symmetric peaks much better.
         */
        fit.setArray( smGrid.data(), smGrid.rows(), smGrid.cols() );
        fit.setGuess( mn, mx - mn, gx, gy, guessWidth / dPix, guessWidth / dPix, 0 );
 
        fit.fit();
 
        // The estimated enter of circular symmetry from the fit.
        x = x0 - maxPix + fit.x0() * dPix;
        y = y0 - maxPix + fit.y0() * 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::imCenterCircleSym Results \n";
        ios << c << "--------------------------------------\n";
        ios << c << "Estimated x-center: " << _x0 - maxPix + fit.x0() * dPix << "\n";
        ios << c << "Estimated y-center: " << _y0 - maxPix + fit.y0() * dPix << "\n";
        ios << c << "Cost half-width: "
            << 0.5 * sigma2fwhm( sqrt( pow( fit.sigma_x(), 2 ) + pow( fit.sigma_y(), 2 ) ) * dPix ) << " pixels\n";
        ios << c << "--------------------------------------\n";
        ios << c << "Setup:\n";
        ios << c << "--------------------------------------\n";
        ios << c << "maxPix: " << maxPix << "\n";
        ios << c << "dPix: " << dPix << "\n";
        ios << c << "minRad: " << minRad << "\n";
        ios << c << "maxRad: " << maxRad << "\n";
        ios << c << "dRad: " << dRad << "\n";
        ios << c << "smWidth: " << smWidth << "\n";
        ios << c << "guessWidth: " << guessWidth << "\n";
        ios << c << "--------------------------------------\n";
        ios << c << "Fit results:\n";
        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 //__imCenterCircleSym_hpp__