imageUtils.hpp

Go to the documentation of this file.

/** \file imageUtils.hpp
 * \author Jared R. Males
 * \brief  Header for the image processing utilities
 * \ingroup image_processing_files
 *
 */
 
//***********************************************************************//
// Copyright 2020 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_imageUtils_hpp
#define improc_imageUtils_hpp
 
#include <cmath>
#include <limits>
#include <type_traits>
 
#include "../math/floatUtils.hpp"
 
#include "imageTransforms.hpp"
 
namespace mx
{
namespace improc
{
 
/// Get the mxlib standard center coordinate of an image
/**
 * \returns 0.5*(rows_cols-1)
 */
template <typename realT = double>
realT imCen( int rows_cols /**< [in] the number of rows or columns in the image */ )
{
    return 0.5 * ( 1.0 * rows_cols - 1 );
}
 
/// Get the mxlib standard center x-coordinate of an image
/**
 * \returns 0.5*(rows-1)
 */
template <typename realT = double, typename imT>
realT imCenX( const imT &im /**< [in] the image to find the center of */ )
{
    return imCen<realT>( im.rows() );
}
 
/// Get the mxlib standard center y-coordinate of an image
/**
 * \returns 0.5*(cols-1)
 */
template <typename realT = double, typename imT>
realT imCenY( const imT &im /**< [in] the image to find the center of */ )
{
    return imCen<realT>( im.cols() );
}
 
/** \ingroup image_utils
 *@{
 */
 
/// Return the configured invalid image-pixel value.
/** The default is a quiet NaN so FITS and image-display tools recognize invalid pixels. Defining
 * `MXLIB_INVALID_NUMBER_VALUE` at compile time selects that finite value instead. CMake propagates its
 * `MXLIB_INVALID_NUMBER_VALUE` cache setting through the installed pkg-config metadata.
 *
 * \returns the invalid-pixel value for `realT`
 */
template <typename realT>
constexpr realT invalidNumber()
{
    static_assert( std::is_floating_point_v<realT>, "invalidNumber requires a floating-point type" );
 
#ifdef MXLIB_INVALID_NUMBER_VALUE
    return static_cast<realT>( MXLIB_INVALID_NUMBER_VALUE );
#else
    return std::numeric_limits<realT>::quiet_NaN();
#endif
}
 
/// Check whether a value represents an invalid image pixel.
/** Detects the configured mxlib invalid-number value as well as NaN and positive or negative infinity. The
 * non-finite classification remains reliable under finite-math-only optimization.
 *
 * \returns true if value is invalid, otherwise false.
 */
template <typename realT>
bool isInvalidPixel( realT value /**< [in] value to test */ )
{
    static_assert( std::is_floating_point_v<realT>, "isInvalidPixel requires a floating-point type" );
 
#ifdef MXLIB_INVALID_NUMBER_VALUE
    return value == invalidNumber<realT>() || !math::isFinite( value );
#else
    return !math::isFinite( value );
#endif
}
 
/// Reflect pixel coordinates across the given center pixel.
/**
 */
template <typename realT>
int reflectImageCoords( int &x1,  ///< [out] the reflected x coordinate
                        int &y1,  ///< [out] the reflected y coordinate
                        int x0,   ///< [in] the input x coordinate
                        int y0,   ///< [in] the input y coordinate
                        realT xc, ///< [in] the center pixel x coordinate
                        realT yc  ///< [in] the center pixel y coordinate
)
{
    x1 = xc - ( x0 - xc );
    y1 = yc - ( y0 - yc );
 
    return 0;
}
 
/// @}
 
/// Zero any NaNs in an image
/**
 * \overload
 *
 * \ingroup eigen_image_processing
 */
template <class imageT, typename valueT>
void zeroNaNs( imageT &im, ///< [in.out] image which will have any NaN pixels set to zero
               valueT val  ///< [in] [optional] The value to set NaN pixels too.  Default is 0.
)
{
    for( int c = 0; c < im.cols(); ++c )
    {
        for( int r = 0; r < im.rows(); ++r )
        {
            if( isInvalidPixel( im( r, c ) ) )
            {
                im( r, c ) = val;
            }
        }
    }
}
 
/// Zero any NaNs in an image
/**
 * \ingroup eigen_image_processing
 */
template <class imageT>
void zeroNaNs( imageT &im /**< [in.out] image which will have any NaN pixels set to zero */ )
{
    typename imageT::Scalar zero = 0;
    zeroNaNs<imageT, typename imageT::Scalar>( im, zero );
}
 
/// Zero any NaNs in an image cube
/** This version fills in a mask with 1s where there were nans, 0s elsewhere.
 *
 * \ingroup eigen_image_processing
 */
template <class cubeT, class maskCubeT>
void zeroNaNCube( cubeT &imc,     /**< [in.out] cube which will have any NaN pixels set to zero */
                  maskCubeT *mask /**< [out] a 1/0 mask with 1 indicating which pixels where nan */
)
{
    if( mask )
    {
        mask->resize( imc.rows(), imc.cols(), imc.planes() );
        mask->setZero();
    }
 
    for( int p = 0; p < imc.planes(); ++p )
    {
        for( int c = 0; c < imc.cols(); ++c )
        {
            for( int r = 0; r < imc.rows(); ++r )
            {
                if( isInvalidPixel( imc.image( p )( r, c ) ) )
                {
                    imc.image( p )( r, c ) = 0;
                    if( mask )
                    {
                        ( *mask ).image( p )( r, c ) = 1;
                    }
                }
            }
        }
    }
}
 
/// Zero any NaNs in an image cube
/**
 * \ingroup eigen_image_processing
 */
template <class cubeT>
void zeroNaNCube( cubeT &imc /**< [in.out] cube which will have any NaN pixels set to zero */ )
{
    return zeroNaNCube<cubeT, cubeT>( imc, nullptr );
}
 
/// Calculate the mean value of an image
/**
 *
 * \returns the mean of the input image
 *
 * \ingroup eigen_image_processing
 */
template <class calcT, class imageT>
calcT imageMean( imageT &im /**< [in] the image of which to calculate the mean*/ )
{
    return static_cast<calcT>( im.sum() ) / ( im.rows() * im.cols() );
}
 
/// Calculate the mean value of an image over a mask
/**
 * \returns the mean of the input image in the masked region
 *
 * \ingroup eigen_image_processing
 */
template <class calcT, class imageT, class maskT>
calcT imageMean( imageT &im, /**< [in] the image of which to calculate the mean*/
                 const maskT &mask /**< [in] a 1/0 mask where 1 defines the good pixels */ )
{
    return static_cast<calcT>( ( im * mask ).sum() ) / ( mask.sum() );
}
 
/// Calculate the variance of an image given its mean
/**
 *
 * \returns the variance of the input image
 *
 * \ingroup eigen_image_processing
 */
template <typename calcT, class imageT>
calcT imageVariance( imageT &im /**< [in] the image of which to calculate the variance*/,
                     calcT mn /**< [in] the mean value of the image w.r.t. which to calcualate the variance */
)
{
    return ( im.template cast<calcT>() - mn ).square().sum() / ( im.rows() * im.cols() );
}
 
/// Calculate the variance of an image given its mean
/**
 *
 * \returns the variance of the input image
 *
 * \ingroup eigen_image_processing
 */
template <typename calcT, class imageT, class maskT>
calcT imageVariance( imageT &im /**< [in] the image of which to calculate the variance*/,
                     calcT mn,         /**< [in] the mean value of the image w.r.t. which to calcualate the variance */
                     const maskT &mask /**< [in] a 1/0 mask where 1 defines the good pixels */
)
{
    return ( im.template cast<calcT>() * mask - mn ).square().sum() / ( mask.sum() );
}
 
/// Calculate the median of an Eigen-like array.
/** Calculates the median of the entire array, allowing for some pixels to be ignored using a mask.
 * Working memory can be retained between calls.
 *
 * \tparam imageT is an Eigen-like type
 * \tparam maskT is an Eigen-like type
 *
 * \returns the median of the unmasked pixels of mat, using \ref mx::math::vectorMedianInPlace().
 *
 * \ingroup eigen_image_processing
 *
 */
template <typename imageT, typename maskT = imageT>
typename imageT::Scalar
imageMedian( const imageT &mat,                             /**< [in] the image to take the median of*/
             const maskT *mask,                             /**< [in] if non-0, a 1/0 mask where 0 pixels are ignored.*/
             std::vector<typename imageT::Scalar> *work = 0 /**< [in] [optional] working memory
                                                                                 can be retained
                                                                       and re-passed. Is resized.*/
)
{
    typename imageT::Scalar med;
 
    bool localWork = false;
    if( work == 0 )
    {
        work = new std::vector<typename imageT::Scalar>;
        localWork = true;
    }
 
    int sz = mat.size();
 
    if( mask )
    {
        sz = mask->sum();
    }
 
    work->resize( sz );
 
    if( mask )
    {
        int ii = 0;
        for( int i = 0; i < mat.rows(); ++i )
        {
            for( int j = 0; j < mat.cols(); ++j )
            {
                if( ( *mask )( i, j ) == 0 )
                {
                    continue;
                }
                ( *work )[ii] = mat( i, j );
                ++ii;
            }
        }
    }
    else
    {
        int ii = 0;
        for( int i = 0; i < mat.rows(); ++i )
        {
            for( int j = 0; j < mat.cols(); ++j )
            {
                ( *work )[ii] = mat( i, j );
                ++ii;
            }
        }
    }
 
    med = math::vectorMedianInPlace( *work );
 
    if( localWork )
    {
        delete work;
    }
 
    return med;
}
 
/// Calculate the median of an Eigen-like array.
/** Calculates the median of the entire array.
 * Working memory can be retained between calls.
 *
 * \tparam imageT is an Eigen-like type
 *
 * \returns the median of the unmasked pixels of mat, using \ref mx::math::vectorMedianInPlace().
 *
 * \ingroup eigen_image_processing
 *
 */
template <typename imageT>
typename imageT::Scalar imageMedian( const imageT &mat, /**< [in] the image to take the median of*/
                                     std::vector<typename imageT::Scalar> *work = 0 /**< [in] [optional] working memory
                                                                                       can be retained and re-passed.*/
)
{
    return imageMedian( mat, static_cast<Eigen::Array<typename imageT::Scalar, -1, -1> *>( nullptr ), work );
}
 
/// Calculate the center of light of an image
/** Note that the sum of the image should be > 0.
 *
 * \ingroup eigen_image_processing
 */
template <typename imageT>
int imageCenterOfLight( typename imageT::Scalar &x, ///< [out] the x coordinate of the center of light [pixels]
                        typename imageT::Scalar &y, ///< [out] the y coordinate of hte center of light [pixels]
                        const imageT &im            ///< [in] the image to centroid
)
{
    x = 0;
    y = 0;
 
    typename imageT::Scalar sum = im.sum();
 
    if( sum == 0 )
    {
        x = 0.5 * ( im.rows() - 1.0 );
        y = 0.5 * ( im.cols() - 1.0 );
        return 0;
    }
 
    for( int j = 0; j < im.cols(); ++j )
    {
        for( int i = 0; i < im.rows(); ++i )
        {
            x += ( i + 1 ) * im( i, j );
            y += ( j + 1 ) * im( i, j );
        }
    }
 
    x = x / sum - 1;
    y = y / sum - 1;
 
    return 0;
}
 
/// Find the maximum in an image at sub-pixel resolution by interpolation
/** Uses imageMagnify() to expand the image to the desired scale.  Because of the
 * scaling used in imageMagnify, the desired scale may not be exact.  As a result
 * the actual scale is returned in scale_x and scale_y.
 *
 * \ingroup eigen_image_processing
 */
template <typename floatT, typename imageT, typename magImageT, typename transformT>
int imageMaxInterp( floatT &x,        ///< [out] the x-position of the maximum, in pixels of the input image
                    floatT &y,        ///< [out] the y-position of the maximum, in pixels of the input image
                    floatT &scale_x,  ///< [in.out] the desired scale or resolution, in pixels < 1, in the x direction.
                                      ///< On output contains the actual scale calculated.
                    floatT &scale_y,  ///< [in.out] the desired scale or resolution, in pixels < 1, in the y direction.
                                      ///< On output contains the actual scale calculated.
                    magImageT &magIm, ///< [in] the magnified image.  This is used as working memory, will be resized.
                    const imageT &im, ///< [in] the image to find the maximum of
                    transformT trans  ///< [in] the transform to use for interpolation
)
{
    floatT magSize_x = ceil( ( im.rows() - 1.0 ) / scale_x ) + 1;
    floatT magSize_y = ceil( ( im.cols() - 1.0 ) / scale_y ) + 1;
 
    floatT mag_x = ( (floatT)magSize_x - 1.0 ) / ( (floatT)im.rows() - 1.0 );
    floatT mag_y = ( (floatT)magSize_y - 1.0 ) / ( (floatT)im.cols() - 1.0 );
 
    scale_x = 1.0 / mag_x;
    scale_y = 1.0 / mag_y;
 
    magIm.resize( magSize_x, magSize_y );
 
    imageMagnify( magIm, im, trans );
 
    int ix, iy;
    magIm.maxCoeff( &ix, &iy );
    x = ix * scale_x;
    y = iy * scale_y;
 
    return 0;
}
 
/// Find the maximum in an image at sub-pixel resolution by cubic convolution interpolation
/** Uses imageMagnify() to expand the image to the desired scale.  Because of the
 * scaling used in imageMagnify, the desired scale may not be exact.  As a result
 * the actual scale is returned in scale_x and scale_y.
 *
 * \ingroup eigen_image_processing
 */
template <typename floatT, typename imageT, typename magImageT>
int imageMaxInterp( floatT &x,        ///< [out] the x-position of the maximum, in pixels of the input image
                    floatT &y,        ///< [out] the y-position of the maximum, in pixels of the input image
                    floatT &scale_x,  ///< [in.out] the desired scale or resolution, in pixels < 1, in the x direction.
                                      ///< On output contains the actual scale calculated.
                    floatT &scale_y,  ///< [in.out] the desired scale or resolution, in pixels < 1, in the y direction.
                                      ///< On output contains the actual scale calculated.
                    magImageT &magIm, ///< [in] the magnified image.  This is used as working memory, will be resized.
                    const imageT &im  ///< [in] the image to find the maximum of
)
{
    return imageMaxInterp( x, y, scale_x, scale_y, magIm, im, cubicConvolTransform<typename imageT::Scalar>() );
}
 
/// Combine two images, each with their own mask defining good pixels.
/** The combined image is made up of the pixels in im1 where mask1 is 1, and the pixels of im2 where mask2 is 1.
 * If a pixel in both mask1 and mask2 has a value of 1, that pixel in the combo is the average of im1 and im2.
 * All other pixels are set to 0 in the combined image.
 *
 * Separate template types are used for each argument to allow references, etc.
 *
 * \tparam imageT the eigen-like array type of the combined image
 * \tparam imageT1 the eigen-like array type of image 1
 * \tparam imageT2 the eigen-like array type of mask 1
 * \tparam imageT3 the eigen-like array type of image 2
 * \tparam imageT4 the eigen-like array type of mask 2
 *
 * \ingroup eigen_image_processing
 */
template <typename imageT, typename imageT1, typename imageT2, typename imageT3, typename imageT4>
void combine2ImagesMasked( imageT &combo,        ///< [out] the combined image.  will be resized.
                           const imageT1 &im1,   ///< [in] the first image
                           const imageT2 &mask1, ///< [in] the mask for the first image
                           const imageT3 &im2,   ///< [in] the second image
                           const imageT4 &mask2  ///< [in] the mask for the second image
)
{
    combo.resize( im1.rows(), im2.cols() );
 
    for( int c = 0; c < combo.cols(); ++c )
    {
        for( int r = 0; r < combo.rows(); ++r )
        {
            if( mask1( r, c ) == 1 && mask2( r, c ) == 0 )
                combo( r, c ) = im1( r, c );
            else if( mask2( r, c ) == 1 & mask1( r, c ) == 0 )
                combo( r, c ) = im2( r, c );
            else if( mask1( r, c ) == 1 && mask2( r, c ) == 1 )
                combo( r, c ) = 0.5 * ( im1( r, c ) + im2( r, c ) );
            else
                combo( r, c ) = 0;
        }
    }
}
 
/// Remove rows and columns
/**
 * \ingroup eigen_image_processing
 */
template <typename eigenT, typename eigenTin>
void removeRowsAndCols( eigenT &out, const eigenTin &in, int st, int w )
{
    out.resize( in.rows() - w, in.cols() - w );
 
    out.topLeftCorner( st, st ) = in.topLeftCorner( st, st );
 
    out.bottomLeftCorner( in.rows() - ( st + w ), st ) = in.bottomLeftCorner( in.rows() - ( st + w ), st );
 
    out.topRightCorner( st, in.cols() - ( st + w ) ) = in.topRightCorner( st, in.cols() - ( st + w ) );
 
    out.bottomRightCorner( in.rows() - ( st + w ), in.cols() - ( st + w ) ) =
        in.bottomRightCorner( in.rows() - ( st + w ), in.cols() - ( st + w ) );
}
 
/// Remove rows
/**
 * \ingroup eigen_image_processing
 */
template <typename eigenT, typename eigenTin>
void removeRows( eigenT &out, const eigenTin &in, int st, int w )
{
    out.resize( in.rows() - w, in.cols() );
 
    out.topLeftCorner( st, in.cols() ) = in.topLeftCorner( st, in.cols() );
 
    out.bottomLeftCorner( in.rows() - ( st + w ), in.cols() ) =
        in.bottomLeftCorner( in.rows() - ( st + w ), in.cols() );
}
 
/// Remove columns
/**
 * \ingroup eigen_image_processing
 */
template <typename eigenT, typename eigenTin>
void removeCols( eigenT &out, const eigenTin &in, int st, int w )
{
    out.resize( in.rows(), in.cols() - w );
 
    out.topLeftCorner( in.rows(), st ) = in.topLeftCorner( in.rows(), st );
 
    out.topRightCorner( in.rows(), in.cols() - ( st + w ) ) = in.topRightCorner( in.rows(), in.cols() - ( st + w ) );
}
 
/** \ingroup image_utils
 *@{
 */
 
/// Copy one image to another, with no transformation
/** This is merely memcpy
 *
 * \returns dest
 */
void *imcpy( void *dest,    ///< [out] the address of the first pixel in the destination image
             void *src,     ///< [in] the address of the first pixel in the source image
             size_t width,  ///< [in] the width in pixels of size szof
             size_t height, ///< [in] the height in pixels of size szof
             size_t szof    ///< [in] the size in bytes of a one pixel
);
 
/// Copy one image to another, flipping up-down
/** This is a reversed row-by-row memcpy
 *
 * \returns dest
 */
void *imcpy_flipUD( void *dest,    ///< [out] the address of the first pixel in the destination image
                    void *src,     ///< [in] the address of the first pixel in the source image
                    size_t width,  ///< [in] the width in pixels of size szof
                    size_t height, ///< [in] the height in pixels of size szof
                    size_t szof    ///< [in] the size in bytes of a one pixel
);
 
/// Copy one image to another, flipping left-right
/**
 *
 * \returns dest
 */
void *imcpy_flipLR( void *dest,    ///< [out] the address of the first pixel in the destination image
                    void *src,     ///< [in] the address of the first pixel in the source image
                    size_t width,  ///< [in] the width in pixels of size szof
                    size_t height, ///< [in] the height in pixels of size szof
                    size_t szof    ///< [in] the size in bytes of a one pixel
);
 
/// Copy one image to another, flipping up-down and left-right
/**
 *
 * \returns dest
 */
void *imcpy_flipUDLR( void *dest,    ///< [out] the address of the first pixel in the destination image
                      void *src,     ///< [in] the address of the first pixel in the source image
                      size_t width,  ///< [in] the width in pixels of size szof
                      size_t height, ///< [in] the height in pixels of size szof
                      size_t szof    ///< [in] the size in bytes of a one pixel
);
 
} // namespace improc
} // namespace mx
 
#endif // improc_imageUtils_hpp