imageFilters.hpp

Go to the documentation of this file.

/** \file imageFilters.hpp
 * \brief Image filters (smoothing, radial profiles, etc.)
 * \ingroup image_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
#ifndef __imageFilters_hpp__
#define __imageFilters_hpp__
 
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <format>
#include <limits>
 
#include "../math/floatUtils.hpp"
#include "../math/gslInterpolator.hpp"
#include "../math/vectorUtils.hpp"
#include "../math/geo.hpp"
 
#include "imageMasks.hpp"
 
namespace mx
{
namespace improc
{
 
/** \addtogroup image_filters_kernels
 * @{
 * The filter function use a kernel that specifies how to filter the image.  Filter kernels, usually
 * denoted as type kernelT below, must have the following interface:
 * \code
 * struct filterKernel
 * {
 *     typedef <eigen-like type> arrayT; // arrayT must have an eigen-like interface
 *
 *     typedef <type> arithT; // the type used for arithmetic, normally `typename arrayT::Scalar`
 *
 *     typedef <verbosity-type> verboseT; // the mxlib verbosity for error reports
 *
 *     //The maxWidth function returns the maximum possible full-width (in either direction) of the kernel
 *     // Called only once for each call to the filter function
 *     int maxWidth() const
 *     {
 *        //returns the maximum half-width given the configuration
 *     }
 *
 *     //The setKernel function is called for each pixel.
 *     void setKernel( arithT x,       // the pixel x-position relative to the image center (not the pixel index)
 *                     arithT y,       // the pixel y-position relative to the image center (not the pixel index)
 *                     arrayT & kernel // the array to be resized and populated
 *                   ) const
 *     {
 *        //This must resize and populate the passed in kernel array each time
 *        //so that it is re-entrant.
 *
 *        //Note: On output kernel array should be normalized so that sum() = 1.0
 *
 *        //Note: the width and height of the kernel array should always be odd
 *     }
 * };
 * \endcode
 * Additionally `kernelT` must be copyable.
 *
 */
/// Symetric Gaussian smoothing kernel
/** \ingroup image_filters_kernels
 *
 */
template <typename _arrayT, size_t _kernW = 4, class _verboseT = verbose::d>
struct gaussKernel
{
    typedef _arrayT arrayT;
    typedef typename _arrayT::Scalar arithT;
    static const int kernW = _kernW;
 
    typedef _verboseT verboseT;
 
    arrayT kernel;
 
    arithT _fwhm;
 
    explicit gaussKernel( arithT fwhm )
    {
        _fwhm = fwhm;
 
        int w = kernW * _fwhm;
 
        if( w % 2 == 0 )
            w++;
 
        kernel.resize( w, w );
 
        arithT kcen = 0.5 * ( w - 1.0 );
 
        arithT sig2 = _fwhm / 2.354820045030949327;
        sig2 *= sig2;
 
        arithT r2;
        for( int i = 0; i < w; ++i )
        {
            for( int j = 0; j < w; ++j )
            {
                r2 = pow( i - kcen, 2 ) + pow( j - kcen, 2 );
                kernel( i, j ) = exp( -r2 / ( 2.0 * sig2 ) );
            }
        }
 
        kernel /= kernel.sum();
    }
 
    int maxWidth() const
    {
        return _kernW * _fwhm;
    }
 
    error_t setKernel( arithT x,           /**< [in] x-coordinate relative to image center */
                       arithT y,           /**< [in] x-coordinate relative to image center */
                       arrayT &kernelArray /**< [in] the array to populate with the kernel, resized */
    ) const
    {
        // Unused parts of interface:
        static_cast<void>( x );
        static_cast<void>( y );
 
        kernelArray = kernel;
 
        return error_t::noerror;
    }
};
 
/// Azimuthally variable boxcar kernel.
/** Averages the image in a boxcar defined by a radial and azimuthal extent.
 *
 * \ingroup image_filters_kernels
 */
template <typename _arrayT, size_t _kernW = 2, class _verboseT = verbose::d>
struct azBoxKernel
{
    typedef _arrayT arrayT;
    typedef typename _arrayT::Scalar arithT;
 
    inline static constexpr int kernW = static_cast<int>( _kernW ); ///< kernel sampling factor.
 
    typedef _verboseT verboseT;
 
    arithT m_radWidth{ 0 }; ///< the half-width of the averaging box, in the radial direction, in pixels.
    arithT m_azWidth{ 0 };  ///< the half-width of the averaging box, in the azimuthal direction, in pixels.
    arithT m_maxAz{ 0 };    ///< maximum azimuthal half-width in radians; 0 means no angular limit.
 
    int m_maxWidth{ 0 };    ///< maximum kernel half-width needed to keep every generated kernel in bounds.
 
    /// Construct a kernel without an angular-position limit.
    azBoxKernel( arithT radWidth, ///< [in] the half-width of the averaging box, in the radial direction, in pixels.
                 arithT azWidth   ///< [in] the half-width of the averaging box, in the azimuthal direction, in pixels.
                 )
        : m_radWidth( fabs( radWidth ) ), m_azWidth( fabs( azWidth ) )
    {
        setMaxWidth();
    }
 
    /// Construct a kernel with an optional angular-position limit.
    azBoxKernel( arithT radWidth, ///< [in] the half-width of the averaging box, in the radial direction, in pixels.
                 arithT azWidth,  ///< [in] the half-width of the averaging box, in the azimuthal direction, in pixels.
                 arithT maxAz     ///< [in] the maximum half-width of the averaging box in the azimuthal direction, in
                                  ///< degrees. >= 0. If 0 or >= 180, then no maximum is enforced.
                 )
        : m_radWidth( fabs( radWidth ) ), m_azWidth( fabs( azWidth ) )
    {
        setMaxWidth();
 
        if( !math::isFinite( maxAz ) )
        {
            m_maxAz = maxAz;
            return;
        }
 
        maxAz = fabs( maxAz );
        if( maxAz >= 180 )
            maxAz = 0; // Larger than 180 means no limit.
 
        m_maxAz = math::dtor( maxAz );
    }
 
    /// Sets the max width based on the configured az and rad widths.
    void setMaxWidth()
    {
        m_maxWidth = 0;
 
        if( kernW <= 0 || !math::isFinite( m_radWidth ) || !math::isFinite( m_azWidth ) ||
            ( m_radWidth == 0 && m_azWidth == 0 ) )
        {
            return;
        }
 
        const arithT maximumDimension =
            kernW * ( std::floor( std::hypot( m_radWidth, m_azWidth ) ) + static_cast<arithT>( 1 ) );
 
        if( !math::isFinite( maximumDimension ) || maximumDimension > std::numeric_limits<int>::max() )
        {
            return;
        }
 
        m_maxWidth = static_cast<int>( maximumDimension / static_cast<arithT>( 2 ) );
    }
 
    /// Get the maximum kernel half-width in either image dimension.
    int maxWidth() const
    {
        return m_maxWidth;
    }
 
    /// Generate a normalized kernel at the requested image-relative coordinate.
    error_t setKernel( arithT x,      /**< [in] x-coordinate relative to image center */
                       arithT y,      /**< [in] x-coordinate relative to image center */
                       arrayT &kernel /**< [in] the array to populate with the kernel, resized */
    ) const
    {
        kernel.resize( 0, 0 );
 
        if( kernW <= 0 || !math::isFinite( m_radWidth ) || !math::isFinite( m_azWidth ) || !math::isFinite( m_maxAz ) ||
            m_maxWidth < 0 )
        {
            return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
                                                           "kernel widths and maxAz must be finite" );
        }
 
        if( !math::isFinite( x ) || !math::isFinite( y ) )
        {
            return internal::mxlib_error_report<verboseT>( error_t::invalidarg, "kernel coordinates must be finite" );
        }
 
        if( m_radWidth == 0 && m_azWidth == 0 )
        {
            kernel.resize( 1, 1 );
            kernel( 0, 0 ) = 1;
            return error_t::noerror;
        }
 
        const arithT rad0 = std::hypot( x, y );
 
        arithT sinq = 0;
        arithT cosq = 1;
        if( rad0 > 0 )
        {
            sinq = y / rad0;
            cosq = x / rad0;
        }
 
        // Only calc q if we're going to use it.
        arithT q = 0;
        if( m_maxAz > 0 && rad0 > 0 )
        {
            q = std::atan2( y, x );
        }
 
        const arithT width =
            kernW * ( std::floor( std::fabs( m_azWidth * sinq ) + std::fabs( m_radWidth * cosq ) ) + 1 );
        const arithT height =
            kernW * ( std::floor( std::fabs( m_azWidth * cosq ) + std::fabs( m_radWidth * sinq ) ) + 1 );
 
        if( !math::isFinite( width ) || !math::isFinite( height ) || width > std::numeric_limits<int>::max() ||
            height > std::numeric_limits<int>::max() )
        {
            return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
                                                           "kernel dimensions exceed the supported size" );
        }
 
        int w = static_cast<int>( width );
        int h = static_cast<int>( height );
 
        // Retain the output pixel itself so any valid angular limit has non-empty support. Increasing an even
        // dimension by one does not change its integer half-width.
        if( m_maxAz > 0 )
        {
            if( w % 2 == 0 )
            {
                ++w;
            }
            if( h % 2 == 0 )
            {
                ++h;
            }
        }
 
        if( w / 2 > m_maxWidth )
        {
            return internal::mxlib_error_report<verboseT>( error_t::sizeerr,
                                                           std::format( "Width half-width bigger than maxWidth. "
                                                                        "This is a bug. Details: "
                                                                        "|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|",
                                                                        kernW,
                                                                        m_radWidth,
                                                                        m_azWidth,
                                                                        m_maxWidth,
                                                                        x,
                                                                        y,
                                                                        rad0,
                                                                        sinq,
                                                                        cosq,
                                                                        w,
                                                                        h ) );
        }
 
        if( h / 2 > m_maxWidth )
        {
            return internal::mxlib_error_report<verboseT>( error_t::sizeerr,
                                                           std::format( "Height half-width bigger than maxWidth. "
                                                                        "This is a bug. Details: "
                                                                        "|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|",
                                                                        kernW,
                                                                        m_radWidth,
                                                                        m_azWidth,
                                                                        m_maxWidth,
                                                                        x,
                                                                        y,
                                                                        rad0,
                                                                        sinq,
                                                                        cosq,
                                                                        w,
                                                                        h ) );
        }
 
        kernel.resize( w, h );
 
        const arithT xcen = 0.5 * ( w - 1.0 );
        const arithT ycen = 0.5 * ( h - 1.0 );
 
        for( int j = 0; j < h; ++j )
        {
            for( int i = 0; i < w; ++i )
            {
                const arithT dx = i - xcen;
                const arithT dy = j - ycen;
                const arithT sampleX = x + dx;
                const arithT sampleY = y + dy;
                const arithT radP = std::hypot( sampleX, sampleY );
 
                if( std::fabs( radP - rad0 ) > m_radWidth )
                {
                    kernel( i, j ) = 0;
                    continue;
                }
 
                const arithT tangentialOffset = dy * cosq - dx * sinq;
                if( std::fabs( tangentialOffset ) <= m_azWidth )
                {
                    if( m_maxAz > 0 ) // Only check this if needed.
                    {
                        arithT q2 = 0;
                        if( radP > 0 )
                        {
                            q2 = std::atan2( sampleY, sampleX );
                        }
 
                        const arithT dq = math::angleDiff<math::radiansT<arithT>>( q, q2 );
                        if( std::fabs( dq ) > m_maxAz )
                        {
                            kernel( i, j ) = 0;
                            continue;
                        }
                    }
                    kernel( i, j ) = 1;
                }
                else
                    kernel( i, j ) = 0;
            }
        }
 
        const arithT ksum = kernel.sum();
        if( !math::isFinite( ksum ) || ksum <= 0 )
        {
            kernel.resize( 0, 0 );
            return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
                                                           std::format( "kernel sum 0 at {},{}", x, y ) );
        }
        kernel /= ksum;
 
        return error_t::noerror;
    }
};
 
/// A kernel that is pre-calculated for the entire image, useful for repeated applications
/** Use this for spatially variable kernels.  It is not needed for, e.g., the symmetric Gaussian.
 *
 * \tparam kernelT is the kernel type. See above for the requirements on kernelT.
 */
template <class kernelT>
struct precalcKernel
{
    typedef kernelT::arrayT arrayT;
    typedef kernelT::arrayT::Scalar arithT;
    typedef kernelT::verboseT verboseT;
 
    kernelT m_kernel;              ///< copied production kernel used to populate the cache.
 
    uint32_t m_rows;               ///< number of image rows represented by the cache.
    uint32_t m_cols;               ///< number of image columns represented by the cache.
    arithT m_xcen;                 ///< pixel x-coordinate of the image center.
    arithT m_ycen;                 ///< pixel y-coordinate of the image center.
 
    int m_maxWidth;                ///< maximum half-width reported by the copied production kernel.
 
    std::vector<arrayT> m_kernels; ///< generated kernels in column-major image-coordinate order.
 
    /// Disallow construction without a production kernel and image geometry.
    precalcKernel() = delete;
 
    /// Pre-calculate a production kernel for every coordinate in an image.
    /** \throws mx::exception when the production kernel cannot generate any cache entry. */
    precalcKernel( const kernelT &kernel, /**< [in] A fully initialized kernel.  Is copied.*/
                   uint32_t rows,         /**< [in] The rows in the images to be filtered*/
                   uint32_t cols,         /**< [in] The columns in the images to be filtered*/
                   arithT xcen,           /**< [in] The pixel x-coordinate of the center of the image*/
                   arithT ycen            /**< [in] The pixel y-coordinate of the center of the image*/
                   )
        : m_kernel( kernel ), m_rows( rows ), m_cols( cols ), m_xcen( xcen ), m_ycen( ycen )
    {
        m_maxWidth = kernel.maxWidth();
 
        m_kernels.resize( static_cast<size_t>( m_rows ) * static_cast<size_t>( m_cols ) );
        size_t n = 0;
        for( uint32_t cc = 0; cc < m_cols; ++cc )
        {
            for( uint32_t rr = 0; rr < m_rows; ++rr )
            {
                const error_t result = m_kernel.setKernel( rr - xcen, cc - ycen, m_kernels[n] );
                if( result != error_t::noerror )
                {
                    throw mx::exception<verboseT>(
                        result,
                        std::format( "failed to pre-calculate kernel at row {}, column {}", rr, cc ) );
                }
                ++n;
            }
        }
    }
 
    /// Get the maximum half-width reported by the cached production kernel.
    int maxWidth() const
    {
        return m_maxWidth;
    }
 
    /// Retrieve the cached kernel at an integral image-relative coordinate.
    error_t setKernel( arithT x,      /**< [in] x-coordinate relative to image center */
                       arithT y,      /**< [in] x-coordinate relative to image center */
                       arrayT &kernel /**< [in] the array to populate with the kernel, resized */
    ) const
    {
        if( m_kernels.size() == 0 )
        {
            return error_t::invalidconfig;
        }
 
        const arithT row = x + m_xcen;
        const arithT column = y + m_ycen;
 
        if( !math::isFinite( row ) || !math::isFinite( column ) || row < 0 || column < 0 || row >= m_rows ||
            column >= m_cols || std::floor( row ) != row || std::floor( column ) != column )
        {
            return error_t::invalidarg;
        }
 
        const size_t n = static_cast<size_t>( column ) * m_rows + static_cast<size_t>( row );
        kernel = m_kernels[n];
 
        return error_t::noerror;
    }
};
 
/// Filter an image with a mean kernel.
/** Applies the kernel to each pixel in the image and sums, storing the filtered result in the output image.
 *
 * \tparam imageOutT the type of the output image (must have an Eigen-like interface)
 * \tparam imageInT the type of the input image (must have an Eigen-like interface)
 * \tparam kernelT is the kernel type (see above for requirements)
 *
 * \ingroup image_filters_kernels
 *
 */
template <typename imageOutT, typename imageInT, typename kernelT>
error_t filterImage( imageOutT &fim,        /**< [out] Contains the filtered image, will be resized*/
                     imageInT im,           /**< [in] the image to be filtered*/
                     const kernelT &kernel, /**< [in] a fully configured kernel object*/
                     int maxr = 0           /**< [in] [opt] the maximum radius from the image center to
                                                            apply the kernel. pixels outside this radius are
                                                            set to 0.*/
)
{
    typedef typename kernelT::verboseT verboseT;
 
    fim.resize( im.rows(), im.cols() );
 
    float xcen = 0.5 * ( im.rows() - 1 );
    float ycen = 0.5 * ( im.cols() - 1 );
 
    if( maxr == 0 )
    {
        maxr = 0.5 * std::min( im.rows(), im.cols() ) - kernel.maxWidth();
    }
 
    int mini = 0.5 * im.rows() - maxr;
    int maxi = 0.5 * im.rows() + maxr;
    int minj = 0.5 * im.cols() - maxr;
    int maxj = 0.5 * im.cols() + maxr;
 
    typename kernelT::arrayT kernelArray;
 
    error_t top_errc = error_t::noerror;
 
    // clang-format off
    #pragma omp parallel private( kernelArray ) // clang-format on
    {
        int im_i, im_j, im_p, im_q;
        int kern_i, kern_j, kern_p, kern_q;
        typename imageOutT::Scalar norm;
 
        error_t errc = error_t::noerror;
 
        // clang-format off
        #pragma omp for // clang-format on
        for( int i = 0; i < im.rows(); ++i )
        {
            if( top_errc != error_t::noerror ) // can't return from omp loop
            {
                continue;
            }
 
            for( int j = 0; j < im.cols(); ++j )
            {
                if( errc != error_t::noerror ) // can't return from omp loop
                {
                    continue;
                }
 
                if( ( i >= mini && i < maxi ) && ( j >= minj && j < maxj ) )
                {
                    errc = kernel.setKernel( i - xcen, j - ycen, kernelArray );
 
                    fim( i, j ) = ( im.block( i - 0.5 * ( kernelArray.rows() - 1 ),
                                              j - 0.5 * ( kernelArray.cols() - 1 ),
                                              kernelArray.rows(),
                                              kernelArray.cols() ) *
                                    kernelArray )
                                      .sum();
                }
                else
                {
                    errc = kernel.setKernel( i - xcen, j - ycen, kernelArray );
 
                    im_i = i - 0.5 * ( kernelArray.rows() - 1 );
                    if( im_i < 0 )
                        im_i = 0;
 
                    im_j = j - 0.5 * ( kernelArray.cols() - 1 );
                    if( im_j < 0 )
                        im_j = 0;
 
                    im_p = im.rows() - im_i;
                    if( im_p > kernelArray.rows() )
                        im_p = kernelArray.rows();
 
                    im_q = im.cols() - im_j;
                    if( im_q > kernelArray.cols() )
                        im_q = kernelArray.cols();
 
                    kern_i = 0.5 * ( kernelArray.rows() - 1 ) - i;
                    if( kern_i < 0 )
                        kern_i = 0;
 
                    kern_j = 0.5 * ( kernelArray.cols() - 1 ) - j;
                    if( kern_j < 0 )
                        kern_j = 0;
 
                    kern_p = kernelArray.rows() - kern_i;
                    if( kern_p > kernelArray.rows() )
                        kern_p = kernelArray.rows();
 
                    kern_q = kernelArray.cols() - kern_j;
                    if( kern_q > kernelArray.cols() )
                        kern_q = kernelArray.cols();
 
                    // Pick only the smallest widths
                    if( im_p < kern_p )
                        kern_p = im_p;
                    if( im_q < kern_q )
                        kern_q = im_q;
 
                    norm = kernelArray.block( kern_i, kern_j, kern_p, kern_q ).sum();
 
                    fim( i, j ) =
                        ( im.block( im_i, im_j, kern_p, kern_q ) * kernelArray.block( kern_i, kern_j, kern_p, kern_q ) )
                            .sum() /
                        norm;
                    if( !std::isfinite( fim( i, j ) ) )
                        fim( i, j ) = 0.0;
                }
            } // for rows
 
            // clang-format off
            #pragma omp critical // clang-format on
            {
                if( errc != error_t::noerror && top_errc == error_t::noerror )
                {
                    top_errc = errc;
                }
            }
 
        } // for cols
    } // pragma omp parallel
 
    return top_errc;
}
 
/// Filter an image with a median kernel.
/** Calculates the median of all pixels corresponding to non-zero pixels in the kernel,
 *  storing the filtered result in the output image.
 *
 * \tparam imageOutT the type of the output image (must have an Eigen-like interface)
 * \tparam imageInT the type of the input image (must have an Eigen-like interface)
 * \tparam kernelT is the kernel type (see above for requirements)
 *
 * \ingroup image_filters_kernels
 *
 */
template <typename imageOutT, typename imageInT, typename kernelT>
void medianFilterImage( imageOutT &fim,        /**< [out] Contains the filtered image, will be resized*/
                        imageInT im,           /**< [in] the image to be filtered*/
                        const kernelT &kernel, /**< [in] a fully configured kernel object*/
                        int maxr = 0,          /**< [in] [opt] the maximum radius from the image center to apply
                                                               the kernel. Psixels outside this radius are set to 0.*/
                        int maxrproc = 1 )
{
    fim.resize( im.rows(), im.cols() );
 
    float xcen = 0.5 * ( im.rows() - 1 );
    float ycen = 0.5 * ( im.cols() - 1 );
 
    if( maxr == 0 )
        maxr = 0.5 * std::min( im.rows(), im.cols() ) - kernel.maxWidth();
 
    int mini = 0.5 * im.rows() - maxr;
    int maxi = 0.5 * im.rows() + maxr;
    int minj = 0.5 * im.cols() - maxr;
    int maxj = 0.5 * im.cols() + maxr;
 
    typename kernelT::arrayT kernelArray;
 
    // clang-format off
    #pragma omp parallel private( kernelArray ) // clang-format on
    {
        int im_i, im_j, im_p, im_q;
        int kern_i, kern_j, kern_p, kern_q;
        typename imageOutT::Scalar norm;
 
        std::vector<typename imageOutT::Scalar> pixels;
 
        // clang-format off
        #pragma omp for // clang-format on
        for( int i = 0; i < im.rows(); ++i )
        {
            for( int j = 0; j < im.cols(); ++j )
            {
                if( ( i >= mini && i < maxi ) && ( j >= minj && j < maxj ) )
                {
                    kernel.setKernel( i - xcen, j - ycen, kernelArray );
 
                    pixels.clear();
                    for( int cc = 0; cc < kernelArray.cols(); ++cc )
                    {
                        for( int rr = 0; rr < kernelArray.rows(); ++rr )
                        {
                            if( kernelArray( rr, cc ) != 0 )
                                pixels.push_back( im.block( i - 0.5 * ( kernelArray.rows() - 1 ),
                                                            j - 0.5 * ( kernelArray.cols() - 1 ),
                                                            kernelArray.rows(),
                                                            kernelArray.cols() )( rr, cc ) );
                        }
                    }
 
                    fim( i, j ) = math::vectorMedianInPlace( pixels );
                }
                else
                {
                    if( maxrproc == 2 )
                    {
                        fim( i, j ) = 0;
                        continue;
                    }
                    kernel.setKernel( i - xcen, j - ycen, kernelArray );
 
                    im_i = i - 0.5 * ( kernelArray.rows() - 1 );
                    if( im_i < 0 )
                        im_i = 0;
 
                    im_j = j - 0.5 * ( kernelArray.cols() - 1 );
                    if( im_j < 0 )
                        im_j = 0;
 
                    im_p = im.rows() - im_i;
                    if( im_p > kernelArray.rows() )
                        im_p = kernelArray.rows();
 
                    im_q = im.cols() - im_j;
                    if( im_q > kernelArray.cols() )
                        im_q = kernelArray.cols();
 
                    kern_i = 0.5 * ( kernelArray.rows() - 1 ) - i;
                    if( kern_i < 0 )
                        kern_i = 0;
 
                    kern_j = 0.5 * ( kernelArray.cols() - 1 ) - j;
                    if( kern_j < 0 )
                        kern_j = 0;
 
                    kern_p = kernelArray.rows() - kern_i;
                    if( kern_p > kernelArray.rows() )
                        kern_p = kernelArray.rows();
 
                    kern_q = kernelArray.cols() - kern_j;
                    if( kern_q > kernelArray.cols() )
                        kern_q = kernelArray.cols();
 
                    // Pick only the smallest widths
                    if( im_p < kern_p )
                        kern_p = im_p;
                    if( im_q < kern_q )
                        kern_q = im_q;
 
                    norm = kernelArray.block( kern_i, kern_j, kern_p, kern_q ).sum();
 
                    pixels.clear();
                    for( int cc = 0; cc < kern_q; ++cc )
                    {
                        for( int rr = 0; rr < kern_p; ++rr )
                        {
                            if( kernelArray.block( kern_i, kern_j, kern_p, kern_q )( rr, cc ) != 0 )
                            {
                                pixels.push_back( im.block( im_i, im_j, kern_p, kern_q )( rr, cc ) );
                            }
                        }
                    }
 
                    fim( i, j ) = math::vectorMedianInPlace( pixels );
 
                    // fim(i,j) = ( im.block(im_i, im_j, kern_p, kern_q) * kernelArray.block(kern_i, kern_j, kern_p,
                    // kern_q )).sum()/norm;
                }
            }
        }
    } // pragma omp parallel
}
 
///@}
 
/// Smooth an image using the mean in a rectangular box, optionally rejecting the highest and lowest values.
/** Calculates the mean value in a rectangular box of imIn, of size meanFullWidth X meanFullWidth and stores it in the
 * corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of the
 * central pair: it contains meanFullWidth/2 pixels before the output pixel and one fewer after it. Pixels without a
 * complete window are not changed in imOut (i.e. you should initialize them before the call).
 *
 * imOut is not re-allocated.
 *
 * If rejectMinMax is `true` then the minimum and maximum value in the box are not included in the mean.  Rejection
 * makes the algorithm somewhat slower, depending on box width.
 *
 * \tparam imageTout is an eigen-like image array
 * \tparam imageTin is an eigen-like image array
 *
 * \returns 0 on success
 * \returns -1 on error.
 *
 * \ingroup image_filters_average
 */
template <typename imageTout, typename imageTin>
int meanSmooth( imageTout &imOut,         /**< [out] the smoothed image. Not re-allocated, and the edge
                                                     pixels are not modified */
                const imageTin &imIn,     /**< [in] the image to smooth */
                int meanFullWidth,        /**< [in] the full width of the smoothing box */
                bool rejectMinMax = false /**< [in] whether to reject the minimum and maximum values */
)
{
    typedef typename imageTout::Scalar scalarT;
 
    if( meanFullWidth <= 0 || meanFullWidth > imIn.rows() || meanFullWidth > imIn.cols() ||
        imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
    {
        return -1;
    }
 
    const int before = meanFullWidth / 2;
    const int after = meanFullWidth - before - 1;
    int nPix = meanFullWidth * meanFullWidth;
 
    if( rejectMinMax && nPix <= 2 )
    {
        return -1;
    }
 
    if( rejectMinMax ) // avoid the branch on every pixel
    {
        nPix -= 2;
        for( int jj = before; jj < imIn.cols() - after; ++jj )
        {
            for( int ii = before; ii < imIn.rows() - after; ++ii )
            {
                scalarT sum = 0;
                scalarT max = std::numeric_limits<scalarT>::lowest();
                scalarT min = std::numeric_limits<scalarT>::max();
                for( int ll = 0; ll < meanFullWidth; ++ll )
                {
                    for( int kk = 0; kk < meanFullWidth; ++kk )
                    {
                        const scalarT value = imIn( ii - before + kk, jj - before + ll );
                        sum += value;
                        if( value > max )
                            max = value;
                        if( value < min )
                            min = value;
                    }
                }
                imOut( ii, jj ) = ( sum - max - min ) / nPix;
            }
        }
    }
    else
    {
        for( int jj = before; jj < imIn.cols() - after; ++jj )
        {
            for( int ii = before; ii < imIn.rows() - after; ++ii )
            {
                scalarT sum = 0;
                for( int ll = 0; ll < meanFullWidth; ++ll )
                {
                    for( int kk = 0; kk < meanFullWidth; ++kk )
                    {
                        sum += imIn( ii - before + kk, jj - before + ll );
                    }
                }
                imOut( ii, jj ) = sum / nPix;
            }
        }
    }
 
    return 0;
}
 
/** \brief Smooth an image using the mean in a rectangular box, optionally rejecting the highest and lowest values.
 * Determines  the location and value of the highest pixel.
 *
 * Calculates the mean value in a rectangular box of imIn, of size meanFullWidth X meanFullWidth and stores it in the
 * corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of the
 * central pair: it contains meanFullWidth/2 pixels before the output pixel and one fewer after it. Pixels without a
 * complete window are not changed in imOut (i.e. you should initialize them before the call).
 *
 * imOut is not re-allocated.
 *
 * If rejectMinMax is `true` then the minimum and maximum value in the box are not included in the mean.  Rejection
 * makes the somewhat slower, depending on box width.
 *
 * This version also determines the location and value of the maximum pixel.  This adds some overhead, maybe on the
 * order of 10% slower than without.
 *
 * \overload
 *
 * \tparam imageTout is an eigen-like image array
 * \tparam imageTin is an eigen-like image array
 *
 * \returns 0 on success
 * \returns -1 on error.
 *
 * \ingroup image_filters_average
 */
template <typename imageTout, typename imageTin>
int meanSmooth(
    imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
    int &xMax,        /**< [out] the x location of the maximum pixel */
    int &yMax,        /**< [out] the y location of the maximum pixel */
    typename imageTout::Scalar &pMax, /**< [out] the value of the maximum pixel */
    const imageTin &imIn,             /**< [in] the image to smooth */
    int meanFullWidth,                /**< [in] the full width of the smoothing box */
    bool rejectMinMax = false         /**< [in] whether to reject the minimum and maximum values */
)
{
    typedef typename imageTout::Scalar scalarT;
 
    xMax = -1;
    yMax = -1;
    pMax = std::numeric_limits<scalarT>::lowest();
 
    if( meanFullWidth <= 0 || meanFullWidth > imIn.rows() || meanFullWidth > imIn.cols() ||
        imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
    {
        return -1;
    }
 
    const int before = meanFullWidth / 2;
    const int after = meanFullWidth - before - 1;
    int nPix = meanFullWidth * meanFullWidth;
 
    if( rejectMinMax && nPix <= 2 )
    {
        return -1;
    }
 
    if( rejectMinMax ) // avoid the branch on every pixel.
    {
        nPix -= 2;
        for( int jj = before; jj < imIn.cols() - after; ++jj )
        {
            for( int ii = before; ii < imIn.rows() - after; ++ii )
            {
                scalarT sum = 0;
                scalarT max = std::numeric_limits<scalarT>::lowest();
                scalarT min = std::numeric_limits<scalarT>::max();
                for( int ll = 0; ll < meanFullWidth; ++ll )
                {
                    for( int kk = 0; kk < meanFullWidth; ++kk )
                    {
                        const scalarT value = imIn( ii - before + kk, jj - before + ll );
                        sum += value;
                        if( value > max )
                            max = value;
                        if( value < min )
                            min = value;
                    }
                }
                imOut( ii, jj ) = ( sum - max - min ) / nPix;
                if( imOut( ii, jj ) > pMax )
                {
                    pMax = imOut( ii, jj );
                    xMax = ii;
                    yMax = jj;
                }
            }
        }
    }
    else
    {
        for( int jj = before; jj < imIn.cols() - after; ++jj )
        {
            for( int ii = before; ii < imIn.rows() - after; ++ii )
            {
                scalarT sum = 0;
                for( int ll = 0; ll < meanFullWidth; ++ll )
                {
                    for( int kk = 0; kk < meanFullWidth; ++kk )
                    {
                        sum += imIn( ii - before + kk, jj - before + ll );
                    }
                }
                imOut( ii, jj ) = sum / nPix;
                if( imOut( ii, jj ) > pMax )
                {
                    pMax = imOut( ii, jj );
                    xMax = ii;
                    yMax = jj;
                }
            }
        }
    }
 
    return 0;
}
 
/// Smooth an image using the median in a rectangular box.  Also Determines the location and value of the highest pixel
/// in the smoothed image.
/** Calculates the median value in a rectangular box of imIn, of size medianFullWidth X medianFullWidth and stores it in
 * the corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of
 * the central pair: it contains medianFullWidth/2 pixels before the output pixel and one fewer after it. The median of
 * an even-sized window is the arithmetic mean of its two central samples. Pixels without a complete window are not
 * changed in imOut (i.e. you should initialize them before the call).
 *
 * imOut is not re-allocated.
 *
 * Also determines the location and value of the maximum pixel.  This is a negligble overhead compared to the median
 * operation.
 *
 *
 * \tparam imageTout is an eigen-like image array
 * \tparam imageTin is an eigen-like image array
 *
 * \returns 0 on success
 * \returns -1 on error.
 *
 * \ingroup  image_filters_average
 */
template <typename imageTout, typename imageTin>
int medianSmooth(
    imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
    int &xMax,        /**< [out] the x location of the maximum pixel */
    int &yMax,        /**< [out] the y location of the maximum pixel */
    typename imageTout::Scalar &pMax, /**< [out] the value of the maximum pixel */
    const imageTin &imIn,             /**< [in] the image to smooth */
    int medianFullWidth               /**< [in] the full width of the smoothing box */
)
{
    typedef typename imageTout::Scalar scalarT;
 
    xMax = -1;
    yMax = -1;
    pMax = std::numeric_limits<scalarT>::lowest();
 
    if( medianFullWidth <= 0 || medianFullWidth > imIn.rows() || medianFullWidth > imIn.cols() ||
        imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
    {
        return -1;
    }
 
    const int before = medianFullWidth / 2;
    const int after = medianFullWidth - before - 1;
    const size_t sampleCount = static_cast<size_t>( medianFullWidth ) * static_cast<size_t>( medianFullWidth );
    std::vector<scalarT> pixs( sampleCount );
 
    for( int jj = before; jj < imIn.cols() - after; ++jj )
    {
        for( int ii = before; ii < imIn.rows() - after; ++ii )
        {
            size_t n = 0;
            for( int ll = 0; ll < medianFullWidth; ++ll )
            {
                for( int kk = 0; kk < medianFullWidth; ++kk )
                {
                    pixs[n] = imIn( ii - before + kk, jj - before + ll );
                    ++n;
                }
            }
 
            imOut( ii, jj ) = math::vectorMedianInPlace( pixs );
            if( imOut( ii, jj ) > pMax )
            {
                pMax = imOut( ii, jj );
                xMax = ii;
                yMax = jj;
            }
        }
    }
 
    return 0;
}
 
/// Smooth an image using the median in a rectangular box.
/** Calculates the median value in a rectangular box of imIn, of size medianFullWidth X medianFullWidth and stores it in
 * the corresponding center pixel of imOut. Even widths use the higher-index member of the central pair as the output
 * pixel. Pixels without a complete window are not changed in imOut (i.e. you should initialize them before the call).
 *
 * imOut is not re-allocated.
 *
 * \overload
 *
 * \tparam imageTout is an eigen-like image array
 * \tparam imageTin is an eigen-like image array
 *
 * \returns 0 on success
 * \returns -1 on error.
 *
 * \ingroup image_filters_average
 */
template <typename imageTout, typename imageTin>
int medianSmooth(
    imageTout &imOut,     /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
    const imageTin &imIn, /**< [in] the image to smooth */
    int medianFullWidth   /**< [in] the full width of the smoothing box */
)
{
    int xMax, yMax;
    typename imageTout::Scalar pMax;
    return medianSmooth( imOut, xMax, yMax, pMax, imIn, medianFullWidth );
}
 
template <typename eigenImT>
void rowEdgeMedSubtract( eigenImT &im, ///< The image to filter
                         int ncols     ///< The number of columns on each side of the image to use as the reference
)
{
    typedef typename eigenImT::Scalar realT;
 
    std::vector<realT> edge( 2 * ncols );
    for( int rr = 0; rr < im.rows(); ++rr )
    {
        for( int cc = 0; cc < ncols; ++cc )
        {
            edge[cc] = im( rr, cc );
        }
 
        for( int cc = 0; cc < ncols; ++cc )
        {
            edge[ncols + cc] = im( rr, im.cols() - ncols + cc );
        }
 
        realT med = math::vectorMedian( edge );
        im.row( rr ) -= med;
    }
 
    return;
}
 
template <typename eigenImT>
void colEdgeMedSubtract( eigenImT &im, ///< The image to filter
                         int nrows     ///< The number of rows on each side of the image to use as the reference
)
{
    typedef typename eigenImT::Scalar realT;
 
    std::vector<realT> edge( 2 * nrows );
    for( int cc = 0; cc < im.cols(); ++cc )
    {
        for( int rr = 0; rr < nrows; ++rr )
        {
            edge[rr] = im( rr, cc );
        }
 
        for( int rr = 0; rr < nrows; ++rr )
        {
            edge[nrows + rr] = im( im.rows() - nrows + rr, cc );
        }
 
        realT med = math::vectorMedian( edge );
        im.col( cc ) -= med;
    }
 
    return;
}
 
//------------ Radial Profile --------------------//
 
template <typename floatT>
struct radval
{
    floatT r;
    floatT v;
};
 
template <typename floatT>
struct radvalRadComp
{
    bool operator()( radval<floatT> rv1, radval<floatT> rv2 )
    {
        return ( rv1.r < rv2.r );
    }
};
 
template <typename floatT>
struct radvalValComp
{
    bool operator()( radval<floatT> rv1, radval<floatT> rv2 )
    {
        return ( rv1.v < rv2.v );
    }
};
 
/// Calculate the the radial profile
/** The median radial profile is calculated by rebinning to a 1 pixel grid.
 *
 *
 * \tparam vecT the std::vector-like type to contain the profile
 * \tparam eigenImT1 the eigen-array-like type of the input image
 * \tparam eigenImT2 the eigen-array-like type of the radius and mask image
 * \tparam eigenImT3 the eigen-array-like type of the mask image
 *
 * \ingroup rad_prof
 */
template <typename vecT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
void radprof(
    vecT &rad,              ///< [out] the radius points for the profile. Should be empty.
    vecT &prof,             ///< [out] the median image value at the corresponding radius. Should be empty.
    const eigenImT1 &im,    ///< [in] the image of which to calculate the profile
    const eigenImT2 &radim, ///< [in] image of radius values per pixel
    const eigenImT3 *mask, ///< [in] [optional] 1/0 mask, only pixels with a value of 1 are included in the profile. Set
                           ///< to 0 to not use.
    bool mean = false,     ///< [in] [optional] set to true to use the mean.  If false (default) the median is used.
    typename eigenImT1::Scalar minr = 0 )
{
    typedef typename eigenImT1::Scalar floatT;
 
    int dim1 = im.rows();
    int dim2 = im.cols();
 
    floatT maxr;
    // floatT minr;
 
    size_t nPix;
    if( mask )
    {
        nPix = mask->sum();
    }
    else
    {
        nPix = dim1 * dim2;
    }
 
    /* A vector of radvals will be sorted, then binned*/
    std::vector<radval<floatT>> rv( nPix );
 
    size_t i = 0;
 
    for( int c = 0; c < im.cols(); ++c )
    {
        for( int r = 0; r < im.rows(); ++r )
        {
            if( mask )
            {
                if( ( *mask )( r, c ) == 0 )
                    continue;
            }
 
            rv[i].r = radim( r, c );
            rv[i].v = im( r, c );
            ++i;
        }
    }
 
    sort( rv.begin(), rv.end(), radvalRadComp<floatT>() );
 
    //    for(auto it=rv.begin(); it != rv.end(); ++it)
    //    {
    //       std::cout << it->r << " " << it->v << "\n";
    //    }
    //
    //    exit(0);
 
    if( minr == 0 )
        minr = rv[0].r;
    maxr = rv.back().r;
 
    /*Now bin*/
    floatT dr = 1;
    floatT r0 = minr;
    floatT r1 = minr + dr;
    int i1 = 0, i2, n;
 
    floatT med;
 
    while( r1 < maxr )
    {
        while( rv[i1].r < r0 )
            ++i1;
        i2 = i1;
        while( rv[i2].r <= r1 )
            ++i2;
 
        if( mean )
        {
            med = 0;
            for( int in = i1; in < i2; ++in )
                med += rv[in].v;
            med /= ( i2 - i1 );
        }
        else
        {
            n = 0.5 * ( i2 - i1 );
 
            std::nth_element( rv.begin() + i1, rv.begin() + i1 + n, rv.begin() + i2, radvalValComp<floatT>() );
 
            med = ( rv.begin() + i1 + n )->v;
 
            // Average two points if even number of points
            if( ( i2 - i1 ) % 2 == 0 )
            {
                med =
                    0.5 *
                    ( med + ( *std::max_element( rv.begin() + i1, rv.begin() + i1 + n, radvalValComp<floatT>() ) ).v );
            }
        }
 
        rad.push_back( .5 * ( r0 + r1 ) );
        prof.push_back( med );
        i1 = i2;
        r0 += dr;
        r1 += dr;
    }
}
 
/// Calculate the the radial profile
/** The median radial profile is calculated by rebinning to a 1 pixel grid.
 * This version calculates a centered radius image.
 *
 * \overload
 *
 * \tparam vecT the std::vector-like type to contain the profile
 * \tparam eigenImT1 the eigen-array-like type of the input image
 * \tparam eigenImT2 the eigen-array-like type of the radius and mask image
 * \tparam eigenImT3 the eigen-array-like type of the mask image
 *
 * \ingroup rad_prof
 */
template <typename vecT, typename eigenImT1, typename eigenImT2>
void radprof(
    vecT &rad,             ///< [out] the radius points for the profile. Should be empty.
    vecT &prof,            ///< [out] the median image value at the corresponding radius. Should be empty.
    const eigenImT1 &im,   ///< [in] the image of which to calculate the profile
    const eigenImT2 &mask, ///< [in] 1/0 mask, only pixels with a value of 1 are included in the profile
    bool mean = false      ///< [in] [optional] set to true to use the mean.  If false (default) the median is used.
)
{
    eigenImage<typename eigenImT1::Scalar> radim;
    radim.resize( im.cols(), im.rows() );
 
    radiusImage( radim );
 
    radprof( rad, prof, im, radim, &mask, mean );
}
 
/// Calculate the the radial profile
/** The median radial profile is calculated by rebinning to a 1 pixel grid.
 * This version calculates a centered radius image.
 *
 * \overload
 *
 * \tparam vecT the std::vector-like type to contain the profile
 * \tparam eigenImT1 the eigen-array-like type of the input image
 *
 * \ingroup rad_prof
 */
template <typename vecT, typename eigenImT1>
void radprof(
    vecT &rad,           ///< [out] the radius points for the profile. Should be empty.
    vecT &prof,          ///< [out] the median image value at the corresponding radius. Should be empty.
    const eigenImT1 &im, ///< [in] the image of which to calculate the profile
    bool mean = false,   ///< [in] [optional] set to true to use the mean.  If false (default) the median is used.
    double dr = 1 )
{
    eigenImage<typename eigenImT1::Scalar> radim;
    radim.resize( im.cols(), im.rows() );
 
    radiusImage( radim );
 
    radprof( rad, prof, im, radim, (eigenImage<typename eigenImT1::Scalar> *)nullptr, mean, dr );
}
 
/// Form a radial profile image, and optionally subtract it from the input
/** The radial profile is calculated using linear interpolation on a 1 pixel grid
 *
 *
 * \tparam radprofT the eigen array type of the output
 * \tparam eigenImT1 the eigen array type of the input image
 * \tparam eigenImT2 the eigen array type of the radius image
 * \tparam eigenImT3 the eigen array type of the mask image
 *
 * \ingroup rad_prof
 */
template <typename radprofT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
void radprofim( radprofT &radprofIm,   ///< [out] the radial profile image.  This will be resized.
                eigenImT1 &im,         ///< [in the image to form the profile of.
                const eigenImT2 &rad,  ///< [in] an array of radius values for each pixel
                const eigenImT3 *mask, /**< [in] [optional 1/0 mask, only pixels with a value of 1 are
                                                            included in the profile. Can be nullptr. */
                bool subtract,         /**< [in] if true, then on ouput im will have had its radial
                                                  profile subtracted. */
                bool mean = false      /**< [in] [optional] set to true to use the mean.
                                                            If false (default) the median is used. */)
{
 
    std::vector<double> med_r, med_v; // Must be double for GSL interpolator
 
    radprof( med_r, med_v, im, rad, mask );
 
    /* And finally, interpolate onto the radius image */
    radprofIm.resize( im.rows(), im.cols() );
 
    math::gslInterpolator<math::gsl_interp_linear<double>> interp( med_r, med_v );
 
    for( int c = 0; c < im.cols(); ++c )
    {
        for( int r = 0; r < im.rows(); ++r )
        {
            if( mask )
            {
                if( ( *mask )( r, c ) == 0 )
                {
                    radprofIm( r, c ) = 0;
                    continue;
                }
            }
 
            radprofIm( r, c ) = interp( ( (double)rad( r, c ) ) );
            if( subtract )
                im( r, c ) -= radprofIm( r, c );
        }
    }
}
 
/// Form a radial profile image, and optionally subtract it from the input
/** The radial profile is calculated using linear interpolation on a 1 pixel grid.
 * This version calculates a centered radius image.
 *
 * \tparam radprofT the eigen array type of the output
 * \tparam eigenImT the eigen array type of the input
 *
 * \ingroup rad_prof
 */
template <typename radprofT, typename eigenImT>
void radprofim(
    radprofT &radprof,     ///< [out] the radial profile image.  This will be resized.
    eigenImT &im,          ///< [in] the image to form the profile of.
    bool subtract = false, ///< [in] [optional] if true, then on ouput im will have had its radial profile subtracted.
    bool mean = false      ///< [in] [optional] set to true to use the mean.  If false (default) the median is used.
)
{
    eigenImage<typename eigenImT::Scalar> rad;
    rad.resize( im.rows(), im.cols() );
 
    radiusImage( rad );
 
    radprofim( radprof, im, rad, (eigenImage<typename eigenImT::Scalar> *)nullptr, subtract );
}
 
/** \ingroup std_prof
 * @{
 */
 
/// Form a standard deviation image, and optionally normalize the input relative to the local mean to form a S/N map.
/** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
 *
 * \tparam eigenImT the eigen array type of the output and non-reference images.  Each image input can be a different
 * type to allow references, etc.
 *
 */
template <typename eigenImT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
void stddevImage( eigenImT &stdIm,       ///< [out] the standard deviation image.  This will be resized.
                  const eigenImT1 &im,   ///< [in] the image to form the standard deviation profile of, never altered.
                  const eigenImT2 &rad,  ///< [in] array of radius values
                  const eigenImT3 &mask, ///< [in] a 1/0 mask.  0 pixels are excluded from the std-dev calculations.
                  typename eigenImT::Scalar minRad, ///< [in] the minimum radius to analyze
                  typename eigenImT::Scalar maxRad, ///< [in] the maximum radius to analyze
                  bool divide ///< [in] if true, the output is the input image minus the interpolated mean profile,
                              ///< divided by the std-dev profile, i.e. a S/N map. default is false.
)
{
    typedef typename eigenImT::Scalar floatT;
 
    int dim1 = im.cols();
    int dim2 = im.rows();
 
    floatT mr = rad.maxCoeff();
 
    /* A vector of radvals will be sorted, then binned*/
    std::vector<radval<floatT>> rv;
    rv.reserve( dim1 * dim2 );
 
    for( int i = 0; i < im.size(); ++i )
    {
        if( mask( i ) == 0 )
            continue;
 
        rv.push_back( { rad( i ), im( i ) } );
    }
 
    sort( rv.begin(), rv.end(), radvalRadComp<floatT>() );
 
    /*Now bin*/
    floatT dr = 1;
    floatT r0 = 0;
    floatT r1 = dr;
    int i1 = 0, i2;
 
    std::vector<double> std_r, std_v, mean_v;
 
    while( r1 < mr )
    {
        while( i1 < rv.size() && rv[i1].r < r0 )
        {
            ++i1;
        }
        if( i1 == rv.size() )
        {
            break;
        }
 
        i2 = i1;
        while( i2 < rv.size() && rv[i2].r <= r1 )
            ++i2;
 
        std::vector<double> vals;
 
        for( int i = i1; i < i2; ++i )
        {
            vals.push_back( rv[i].v );
        }
 
        const double mean = math::vectorMean( vals );
        std_r.push_back( .5 * ( r0 + r1 ) );
        mean_v.push_back( mean );
        std_v.push_back( std::sqrt( math::vectorVariance( vals, mean ) ) );
        i1 = i2;
        r0 += dr;
        r1 += dr;
    }
 
    /* And finally, interpolate onto the radius image */
    stdIm.resize( dim1, dim2 );
    math::gslInterpolator<math::gsl_interp_linear<double>> stdInterp( std_r, std_v );
    math::gslInterpolator<math::gsl_interp_linear<double>> meanInterp( std_r, mean_v );
 
    for( int i = 0; i < dim1; ++i )
    {
        for( int j = 0; j < dim2; ++j )
        {
            if( rad( i, j ) < minRad || rad( i, j ) > maxRad )
            {
                stdIm( i, j ) = 0;
            }
            else
            {
                stdIm( i, j ) = stdInterp( ( (double)rad( i, j ) ) );
                if( divide )
                    stdIm( i, j ) = ( im( i, j ) - meanInterp( ( (double)rad( i, j ) ) ) ) / stdIm( i, j );
            }
        }
    }
}
 
/// Form a standard deviation image, and optionally normalize the input relative to the local mean to form a S/N map.
/** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
 *
 * This version creates a radius map on each call, and calls the above version.  This should not
 * be used for repeated alls, rather create a radius map ahead of time.
 *
 * \overload
 *
 * \tparam eigenImT the eigen array type of the output and non-reference images
 *
 */
template <typename eigenImT, typename eigenImT1, typename eigenImT2>
void stddevImage(
    eigenImT &stdIm,                  ///< [out] the standard deviation image.  This will be resized.
    const eigenImT1 &im,              ///< [in] the image to form the standard deviation profile of, never altered.
    const eigenImT2 &mask,            ///< [in] a 1/0 mask.  0 pixels are excluded from the std-dev calculations.
    typename eigenImT::Scalar minRad, ///< [in] the minimum radius to analyze
    typename eigenImT::Scalar maxRad, ///< [in] the maximum radius to analyze
    bool divide = false               ///< [in] [optional] if true, the output is the input image minus the interpolated
                        ///< mean profile, divided by the std-dev profile, i.e. a S/N map. default is false.
)
{
    int dim1 = im.cols();
    int dim2 = im.rows();
 
    eigenImage<typename eigenImT::Scalar> rad;
    rad.resize( dim1, dim2 );
 
    radiusImage( rad );
    stddevImage( stdIm, im, rad, mask, minRad, maxRad, divide );
}
 
/// Form a standard deviation image for each image in a cube, and optionally normalize it relative to the local mean
/// to form a S/N map cube.
/** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
 *
 *
 * \tparam eigencubeT is the eigen cube type of the input and output cubes.
 * \tparam eigenImT the eigen array type of the output and non-reference images.
 *
 */
template <typename eigenCubeT, typename eigenCubeT1, typename eigenCubeT2, typename radT1, typename radT2>
void stddevImageCube(
    eigenCubeT &stdImc,          /**< [out] the standard deviation image cube.  This will be resized. */
    const eigenCubeT1 &imc,      /**< [in] the image cube to form the standard deviation profile of. */
    const eigenCubeT2 &maskCube, /**< [in] a 1/0 mask.  0 pixels are excluded from the std-dev calculations. */
    radT1 minRad,                /**< [in] the minimum radius to analyze */
    radT2 maxRad,                /**< [in] the maximum radius to analyze */
    bool divide = false          /**< [in] [optional] if true, the output is the input image minus the interpolated
                                           mean profile, divided by the std-dev profile, i.e. a S/N map. default is false.*/
)
{
    int dim1 = imc.cols();
    int dim2 = imc.rows();
 
    typename eigenCubeT::Scalar minRadF = minRad;
    typename eigenCubeT::Scalar maxRadF = maxRad;
 
    eigenImage<typename eigenCubeT::Scalar> rad;
    rad.resize( dim1, dim2 );
 
    radiusImage( rad );
 
    stdImc.resize( imc.rows(), imc.cols(), imc.planes() );
 
    // #pragma omp parallel for
    for( int i = 0; i < imc.planes(); ++i )
    {
        eigenImage<typename eigenCubeT::Scalar> im, stdIm, mask;
 
        im = imc.image( i );
        mask = maskCube.image( i );
 
        stddevImage( stdIm, im, rad, mask, minRadF, maxRadF, divide );
 
        stdImc.image( i ) = stdIm;
    }
}
 
///@}
 
} // namespace improc
} // namespace mx
 
#endif //__imageFilters_hpp__