circleOuterpix.hpp

Go to the documentation of this file.

/** \file circleOuterpix.hpp
 * \brief Declares and defines a class for finding the edge of a circle mask
 * \ingroup image_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2021 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 improc_circleOuterpix_hpp
#define improc_circleOuterpix_hpp
 
#include "eigenImage.hpp"
 
#include "imageUtils.hpp"
 
namespace mx
{
namespace improc
{
 
/// Find the center and the outermost pixels of a circular mask, giving an estimate of the radius
/** Takes a 1/0 mask of a filled circle or annulus, and calculates the center coordinate by
 * center-of-light, then finds the outermost 1 pixel in each row and column.  Calculates the average radius
 * of these pixels w.r.t. to the c.o.l. coords, as well as w.r.t. the average coords of the outer pixels
 * themselves.
 *
 * Also produces an image containing the outer pixels.
 *
 * \returns 0 on success
 *
 * \tparam realT the real floating-point type for calculations
 *
 * \ingroup image_utils
 *
 */
template <typename realT>
int circleOuterpix(
    realT &x0,               /// < [out] x coord of center of the mask by center-of-light
    realT &y0,               /// < [out] y coord of center of the mask by center-of-light
    realT &avgr0,            /// < [out] average radius of edge pixels w.r.t. center-of-light
    realT &avgx,             /// < [out] x coord of center of the mask by average edge pixel
    realT &avgy,             /// < [out] y coord of center of the mask by average edge pixel
    realT &avgr,             /// < [out] average radius of edge pixels w.r.t. to average edge pixel
    eigenImage<realT> &circ, /// < [out] image showing the edge pixels with value 1, all other pixels 0.  Resized.
    const eigenImage<realT> &masked /// < [in] an image with a (roughly) circular mask of 1s
)
{
 
    circ.resize( masked.rows(), masked.cols() );
    circ.setZero();
 
    mx::improc::imageCenterOfLight( x0, y0, masked );
 
    for( size_t i = 0; i < masked.rows(); ++i )
    {
        realT ndy = 0;
        realT pdy = 0;
 
        int nj = -1;
        int pj = -1;
 
        for( size_t j = 0; j < masked.cols(); ++j )
        {
            if( masked( i, j ) == 1 )
            {
                if( j - y0 < ndy )
                {
                    ndy = j - y0;
                    nj = j;
                }
                if( j - y0 > pdy )
                {
                    pdy = j - y0;
                    pj = j;
                }
            }
        }
 
        if( nj != -1 )
            circ( i, nj ) = 1;
        if( pj != -1 )
            circ( i, pj ) = 1;
    }
 
    for( size_t j = 0; j < masked.cols(); ++j )
    {
        realT ndx = 0;
        realT pdx = 0;
 
        int ni = -1;
        int pi = -1;
 
        for( size_t i = 0; i < masked.rows(); ++i )
        {
            if( masked( i, j ) == 1 )
            {
                if( i - x0 < ndx )
                {
                    ndx = i - x0;
                    ni = i;
                }
                if( i - x0 > pdx )
                {
                    pdx = i - x0;
                    pi = i;
                }
            }
        }
 
        if( ni != -1 )
            circ( ni, j ) = 1;
        if( pi != -1 )
            circ( pi, j ) = 1;
    }
 
    int npix = 0;
    avgx = 0;
    avgy = 0;
 
    for( size_t i = 0; i < masked.rows(); ++i )
    {
        for( size_t j = 0; j < (size_t)masked.cols(); ++j )
        {
            if( circ( i, j ) == 1 )
            {
                ++npix;
                avgx += i;
                avgy += j;
            }
        }
    }
 
    avgx /= npix;
    avgy /= npix;
 
    avgr = 0;
    avgr0 = 0;
    for( size_t i = 0; i < (size_t)circ.rows(); ++i )
    {
        for( size_t j = 0; j < (size_t)circ.cols(); ++j )
        {
            if( circ( i, j ) == 1 )
            {
                avgr += sqrt( pow( i - avgx, 2 ) + pow( j - avgy, 2 ) );
                avgr0 += sqrt( pow( i - x0, 2 ) + pow( j - y0, 2 ) );
            }
        }
    }
 
    avgr /= npix;
    avgr0 /= npix;
 
    return 0;
}
 
} // namespace improc
} // namespace mx
 
#endif // improc_circleOuterpix_hpp