eigenCube.hpp

Go to the documentation of this file.

/** \file eigenCube.hpp
 * \brief An image cube with an Eigen API
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup image_processing_files
 *
 */
 
#ifndef eigenCube_hpp
#define eigenCube_hpp
 
#include <vector>
 
#pragma GCC system_header
#include <Eigen/Dense>
 
#include "../error/exception.hpp"
#include "../math/vectorUtils.hpp"
#include "eigenImage.hpp"
#include "imageUtils.hpp"
 
namespace mx
{
namespace improc
{
 
/// An image cube with an Eigen-like API
/** \ingroup eigen_image_processing
 */
template <typename dataT>
class eigenCube
{
  public:
    typedef bool is_eigenCube;
    typedef dataT Scalar;
 
    typedef typename Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>::Index Index;
 
    typedef Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> imageRef;
 
    typedef Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic> imageT;
 
  protected:
    Index _rows;
    Index _cols;
    Index _planes;
 
    dataT *m_data{ nullptr };
 
    bool _owner;
 
    /// Validate a mask and good-pixel threshold for a masked combination.
    template <typename eigenCubeT>
    void
    validateMaskedCombination( const eigenCubeT &mask, /**< [in] mask cube whose dimensions must match this cube. */
                               double minGoodFract     /**< [in] inclusive minimum fraction of good pixels in [0, 1]. */
    ) const;
 
    /// Validate the number of weights for a weighted combination.
    void
    validateWeights( const std::vector<dataT> &weights /**< [in] weights whose count must match the number of planes. */
    ) const;
 
    /// Test whether a nonempty set of good pixels meets the requested fraction.
    bool hasEnoughGoodPixels( size_t goodPixels,  /**< [in] number of unmasked pixels. */
                              double minGoodFract /**< [in] inclusive minimum fraction of good pixels. */
    ) const;
 
  public:
    eigenCube();
 
    /// C'tor which will allocate space..
    /**
     */
    eigenCube( Index nrows,  ///< [in] Number of rows in the cube
               Index ncols,  ///< [in] Number of columns the cube
               Index nplanes ///< [in] Number of planes in the cube
    );
 
    /// C'tor taking an existing array as an argument.
    /** The existing array is used as is, as in a map, and ownership is not taken.  You are
     * responsible for memory management, e.g. free-ing this array.
     */
    eigenCube( dataT *ndata,  ///< [in] Allocated array with a cube of nrows x ncols x nplanes
               size_t nrows,  ///< [in] Number of rows in the cube
               size_t ncols,  ///< [in] Number of columns the cube
               size_t nplanes ///< [in] Number of planes in the cube
    );
 
    /// Copy a cube into independently owned storage.
    eigenCube( const eigenCube<dataT> &ec /**< [in] cube to copy. */ );
 
    /// Move a cube, transferring any owned storage.
    eigenCube( eigenCube<dataT> &&ec /**< [in,out] cube to move from. */ ) noexcept;
 
    ~eigenCube();
 
    void setZero();
 
    /// Copy a cube into independently owned storage.
    eigenCube<dataT> &operator=( const eigenCube<dataT> &ec /**< [in] cube to copy. */ );
 
    /// Move a cube, transferring any owned storage.
    eigenCube<dataT> &operator=( eigenCube<dataT> &&ec /**< [in,out] cube to move from. */ ) noexcept;
 
    void shallowCopy( eigenCube<dataT> &src, bool takeOwner = false );
 
    /// De-allocate and set all sizes to 0
    void clear();
 
    void resize( int r, int c, int p );
 
    void resize( int r, int c );
 
    dataT *data();
 
    const dataT *data() const;
 
    const Index rows() const;
 
    const Index cols() const;
 
    const Index planes() const;
 
    /// Returns a 2D Eigen::Eigen::Map pointed at the entire cube.
    Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> cube();
 
    /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
    Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> image( Index n /**< [in] the image number */ );
 
    /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
    const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>
    image( Index n /**< [in] the image number */ ) const;
 
    /// Returns an Eigen::Eigen::Map-ed vector of the pixels at the given coordinate
    Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
               Eigen::Unaligned,
               Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
    pixel( Index i, Index j );
 
    /// Return an Eigen::Eigen::Map of the cube where each image is a vector
    Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> asVectors();
 
    /// Calculate the covariance matrix of the images in the cube
    void Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv );
 
    /// Calculate the sum image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     */
    template <typename eigenT>
    void sum( eigenT &mim /**< [out] the resultant sum image.  Is resized. */ );
 
    /// Calculate the mean image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     */
    template <typename eigenT>
    void mean( eigenT &mim /**< [out] the resultant mean image.  Is resized. */ );
 
    /// Calculate the mean image of the cube with a mask.
    /**
     * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
     * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
     *
     * \tparam eigenT an Eigen-like type.
     * \tparam eigenCubeT an eigenCube type.
     *
     * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
     * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
     */
    template <typename eigenT, typename eigenCubeT>
    void mean( eigenT &mim,              ///< [out] the resultant mean image.  Is resized.
               eigenCubeT &mask,         /**< [in] a mask cube.  Only pixels with value 1 are included in the
                                                   mean calculation. */
               double minGoodFract = 0.0 /**< [in] inclusive minimum fraction of good pixels in [0, 1]. */
    );
 
    /// Calculate the weighted mean image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     *
     * \throws mx::exception with error_t::sizeerr if the weight count does not match the cube planes.
     */
    template <typename eigenT>
    void mean( eigenT &mim,                ///< [out] the resultant mean image. Is resized.
               std::vector<dataT> &weights ///< [in] a vector of weights to use for calculating the mean
    );
 
    /// Calculate the weighted mean image of the cube, with a mask cube.
    /**
     * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
     * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
     *
     * \tparam eigenT an Eigen-like type.
     * \tparam eigenCubeT an eigenCube type.
     *
     * \throws mx::exception with error_t::sizeerr if the mask dimensions or weight count do not match the cube.
     * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
     */
    template <typename eigenT, typename eigenCubeT>
    void mean( eigenT &mim,                 ///< [out] the resultant mean image. Is resized.
               std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
               eigenCubeT &mask, ///< [in] a mask cube.  Only pixels with value 1 are included in the mean calculation.
               double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
    );
 
    /// Calculate the median image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     */
    template <typename eigenT>
    void median( eigenT &mim /**< [out] the resultant median image. Is resized. */ );
 
    /// Calculate the median image of the cube with a mask.
    /**
     * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
     * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
     *
     * \tparam eigenT an Eigen-like type.
     * \tparam eigenCubeT an eigenCube type.
     *
     * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
     * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
     */
    template <typename eigenT, typename eigenCubeT>
    void median( eigenT &mim,              ///< [out] the resultant median image. Is resized.
                 eigenCubeT &mask,         ///< [in] a mask cube. Only pixels with value 1 are included.
                 double minGoodFract = 0.0 /**< [in] the inclusive minimum fraction of good pixels. */
    );
 
    /// Calculate the sigma clipped mean image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     */
    template <typename eigenT>
    void sigmaMean( eigenT &mim, ///< [out] the resultant mean image
                    Scalar sigma ///< [in] the sigma value at which to clip.
    );
 
    /// Calculate the sigma clipped mean image of the cube, with a mask cube
    /**
     * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
     * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
     *
     * \tparam eigenT an Eigen-like type.
     * \tparam eigenCubeT an eigenCube type.
     *
     * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
     * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
     */
    template <typename eigenT, typename eigenCubeT>
    void
    sigmaMean( eigenT &mim,      ///< [out] the resultant mean image.  Is resized.
               eigenCubeT &mask, ///< [in] a mask cube.  Only pixels with value 1 are included in the mean calculation.
               Scalar sigma,     ///< [in] the sigma value at which to clip.
               double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
    );
 
    /// Calculate the sigma clipped weighted mean image of the cube
    /**
     * \tparam eigenT an Eigen-like type.
     *
     * \throws mx::exception with error_t::sizeerr if the weight count does not match the cube planes.
     */
    template <typename eigenT>
    void sigmaMean( eigenT &mim,                 ///< [out] the resultant mean image. Is resized.
                    std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
                    Scalar sigma                 ///< [in] the sigma value at which to clip.
    );
 
    /// Calculate the sigma clipped weighted mean image of the cube, with a mask cube.
    /**
     * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
     * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
     *
     * \tparam eigenT an Eigen-like type.
     * \tparam eigenCubeT an eigenCube type.
     *
     * \throws mx::exception with error_t::sizeerr if the mask dimensions or weight count do not match the cube.
     * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
     */
    template <typename eigenT, typename eigenCubeT>
    void
    sigmaMean( eigenT &mim,                 ///< [out] the resultant mean image. Is resized.
               std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
               eigenCubeT &mask, ///< [in] a mask cube.  Only pixels with value 1 are included in the mean calculation.
               Scalar sigma,     ///< [in] the sigma value at which to clip.
               double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
    );
};
 
template <typename dataT>
template <typename eigenCubeT>
void eigenCube<dataT>::validateMaskedCombination( const eigenCubeT &mask, double minGoodFract ) const
{
    if( mask.rows() != _rows || mask.cols() != _cols || mask.planes() != _planes )
    {
        throw mx::exception( error_t::sizeerr, "mask dimensions must match cube dimensions" );
    }
 
    if( !math::isFinite( minGoodFract ) || minGoodFract < 0.0 || minGoodFract > 1.0 )
    {
        throw mx::exception( error_t::invalidarg, "minGoodFract must be finite and in [0, 1]" );
    }
}
 
template <typename dataT>
void eigenCube<dataT>::validateWeights( const std::vector<dataT> &weights ) const
{
    if( weights.size() != static_cast<size_t>( _planes ) )
    {
        throw mx::exception( error_t::sizeerr, "weights size must match cube planes" );
    }
}
 
template <typename dataT>
bool eigenCube<dataT>::hasEnoughGoodPixels( size_t goodPixels, double minGoodFract ) const
{
    return goodPixels > 0 && static_cast<double>( goodPixels ) >= minGoodFract * static_cast<double>( _planes );
}
 
template <typename dataT>
eigenCube<dataT>::eigenCube()
{
    _rows = 0;
    _cols = 0;
    _planes = 0;
    _owner = false;
}
 
template <typename dataT>
eigenCube<dataT>::eigenCube( Index nrows, Index ncols, Index nplanes )
{
    _rows = nrows;
    _cols = ncols;
    _planes = nplanes;
 
    m_data = new dataT[_rows * _cols * _planes];
    _owner = true;
}
 
template <typename dataT>
eigenCube<dataT>::eigenCube( dataT *ndata, size_t nrows, size_t ncols, size_t nplanes )
{
    _rows = nrows;
    _cols = ncols;
    _planes = nplanes;
 
    m_data = ndata;
    _owner = false;
}
 
template <typename dataT>
eigenCube<dataT>::eigenCube( const eigenCube<dataT> &ec ) : eigenCube()
{
    *this = ec;
}
 
template <typename dataT>
eigenCube<dataT>::eigenCube( eigenCube<dataT> &&ec ) noexcept
    : _rows( ec._rows ), _cols( ec._cols ), _planes( ec._planes ), m_data( ec.m_data ), _owner( ec._owner )
{
    ec._rows = 0;
    ec._cols = 0;
    ec._planes = 0;
    ec.m_data = nullptr;
    ec._owner = false;
}
 
template <typename dataT>
eigenCube<dataT>::~eigenCube()
{
    if( _owner && m_data )
    {
        delete[] m_data;
    }
}
 
template <typename dataT>
void eigenCube<dataT>::setZero()
{
    int N = _rows * _cols * _planes;
 
    for( int i = 0; i < N; ++i )
        m_data[i] = ( (dataT)0 );
}
 
template <typename dataT>
eigenCube<dataT> &eigenCube<dataT>::operator=( const eigenCube<dataT> &ec )
{
    if( this == &ec )
    {
        return *this;
    }
 
    resize( ec.rows(), ec.cols(), ec.planes() );
 
    int N = _rows * _cols * _planes;
 
    for( int i = 0; i < N; ++i )
        m_data[i] = ec.m_data[i];
 
    return *this;
}
 
template <typename dataT>
eigenCube<dataT> &eigenCube<dataT>::operator=( eigenCube<dataT> &&ec ) noexcept
{
    if( this == &ec )
    {
        return *this;
    }
 
    clear();
 
    _rows = ec._rows;
    _cols = ec._cols;
    _planes = ec._planes;
    m_data = ec.m_data;
    _owner = ec._owner;
 
    ec._rows = 0;
    ec._cols = 0;
    ec._planes = 0;
    ec.m_data = nullptr;
    ec._owner = false;
 
    return *this;
}
 
template <typename dataT>
void eigenCube<dataT>::shallowCopy( eigenCube<dataT> &src, bool takeOwner )
{
    if( this == &src )
    {
        return;
    }
 
    clear();
 
    _rows = src._rows;
    _cols = src._cols;
    _planes = src._planes;
    m_data = src.m_data;
 
    if( takeOwner == true )
    {
        _owner = src._owner;
        src._rows = 0;
        src._cols = 0;
        src._planes = 0;
        src.m_data = nullptr;
        src._owner = false;
    }
    else
    {
        _owner = false;
    }
}
 
template <typename dataT>
void eigenCube<dataT>::clear()
{
    if( _owner && m_data )
    {
        delete[] m_data;
    }
 
    _rows = 0;
    _cols = 0;
    _planes = 0;
    m_data = nullptr;
    _owner = false;
}
 
template <typename dataT>
void eigenCube<dataT>::resize( int r, int c, int p )
{
    clear();
 
    _rows = r;
    _cols = c;
    _planes = p;
 
    m_data = new dataT[_rows * _cols * _planes];
    _owner = true;
}
 
template <typename dataT>
void eigenCube<dataT>::resize( int r, int c )
{
    resize( r, c, 1 );
}
 
template <typename dataT>
dataT *eigenCube<dataT>::data()
{
    return m_data;
}
 
template <typename dataT>
const dataT *eigenCube<dataT>::data() const
{
    return m_data;
}
 
template <typename dataT>
const typename eigenCube<dataT>::Index eigenCube<dataT>::rows() const
{
    return _rows;
}
 
template <typename dataT>
const typename eigenCube<dataT>::Index eigenCube<dataT>::cols() const
{
    return _cols;
}
 
template <typename dataT>
const typename eigenCube<dataT>::Index eigenCube<dataT>::planes() const
{
    return _planes;
}
 
template <typename dataT>
Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::cube()
{
    return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
}
 
template <typename dataT>
Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n )
{
    return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
}
 
template <typename dataT>
const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n ) const
{
    return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
}
 
template <typename dataT>
Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
           Eigen::Unaligned,
           Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
eigenCube<dataT>::pixel( Index i, Index j )
{
    return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
                      Eigen::Unaligned,
                      Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(
        m_data + j * _rows + i,
        _planes,
        1,
        Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>( 0, _rows * _cols ) );
}
 
template <typename dataT>
Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::asVectors()
{
    return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
}
 
template <typename dataT>
void eigenCube<dataT>::Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv )
{
    cv = asVectors().matrix().transpose() * asVectors().matrix();
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::sum( eigenT &mim )
{
    mim.resize( _rows, _cols );
 
#pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
    for( Index i = 0; i < _rows; ++i )
    {
        for( Index j = 0; j < _cols; ++j )
        {
            mim( i, j ) = pixel( i, j ).sum();
        }
    }
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::mean( eigenT &mim )
{
    mim.resize( _rows, _cols );
 
#pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
    for( Index i = 0; i < _rows; ++i )
    {
        for( Index j = 0; j < _cols; ++j )
        {
            mim( i, j ) = pixel( i, j ).mean();
        }
    }
}
 
template <typename dataT>
template <typename eigenT, typename eigenCubeT>
void eigenCube<dataT>::mean( eigenT &mim, eigenCubeT &mask, double minGoodFract )
{
    validateMaskedCombination( mask, minGoodFract );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel
    {
        std::vector<dataT> work;
 
#pragma omp for
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
 
                work.clear();
 
                for( Index k = 0; k < _planes; ++k )
                {
                    if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
                    {
                        work.push_back( ( pixel( i, j ) )( k, 0 ) );
                    }
                }
 
                if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
                {
                    mim( i, j ) = math::vectorMean( work );
                }
                else
                {
                    mim( i, j ) = invalidNumber<Scalar>();
                }
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights )
{
    validateWeights( weights );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel num_threads( Eigen::nbThreads() )
    {
        std::vector<Scalar> work;
 
#pragma omp for schedule( static, 10 )
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.resize( _planes ); // work could be smaller after sigmaMean
                // int ii=0;
                for( int k = 0; k < _planes; ++k )
                    work[k] = ( pixel( i, j ) )( k, 0 );
 
                mim( i, j ) = math::vectorMean( work, weights );
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT, typename eigenCubeT>
void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, double minGoodFract )
{
    validateMaskedCombination( mask, minGoodFract );
    validateWeights( weights );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel num_threads( Eigen::nbThreads() )
    {
        std::vector<Scalar> work, wwork;
 
#pragma omp for schedule( static, 10 )
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.clear();
                wwork.clear();
 
                // int ii=0;
 
                for( Index k = 0; k < _planes; ++k )
                {
                    if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
                    {
                        work.push_back( ( pixel( i, j ) )( k, 0 ) );
                        wwork.push_back( weights[k] );
                    }
                }
                if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
                {
                    mim( i, j ) = math::vectorMean( work, wwork );
                }
                else
                    mim( i, j ) = invalidNumber<Scalar>();
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::median( eigenT &mim )
{
    mim.resize( _rows, _cols );
 
#pragma omp parallel for schedule( static, 10 ) num_threads( Eigen::nbThreads() )
    for( Index i = 0; i < _rows; ++i )
    {
        // Eigen::Eigen::Array<Scalar, Eigen::Eigen::Dynamic, Eigen::Eigen::Dynamic> work;
        std::vector<Scalar> work;
        for( Index j = 0; j < _cols; ++j )
        {
            mim( i, j ) = imageMedian( pixel( i, j ), &work );
        }
    }
}
 
template <typename dataT>
template <typename eigenT, typename eigenCubeT>
void eigenCube<dataT>::median( eigenT &mim, eigenCubeT &mask, double minGoodFract )
{
    validateMaskedCombination( mask, minGoodFract );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel
    {
        std::vector<Scalar> work;
 
#pragma omp for
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.clear();
 
                for( Index k = 0; k < _planes; ++k )
                {
                    if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
                    {
                        work.push_back( ( pixel( i, j ) )( k, 0 ) );
                    }
                }
 
                if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
                {
                    mim( i, j ) = math::vectorMedianInPlace( work );
                }
                else
                {
                    mim( i, j ) = invalidNumber<Scalar>();
                }
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::sigmaMean( eigenT &mim, dataT sigma )
{
    mim.resize( _rows, _cols );
 
    // clang-format off
    #pragma omp parallel //clang-format on
    {
        std::vector<Scalar> work;
 
        // clang-format off
        #pragma omp for // clang-format off
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.resize( _planes );
                for( int k = 0; k < _planes; ++k )
                {
                    work[k] = ( pixel( i, j ) )( k, 0 );
                }
 
                mim( i, j ) = math::vectorSigmaMean( work, sigma );
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT, typename eigenCubeT>
void eigenCube<dataT>::sigmaMean( eigenT &mim, eigenCubeT &mask, dataT sigma, double minGoodFract )
{
    validateMaskedCombination( mask, minGoodFract );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel
    {
        std::vector<Scalar> work;
 
#pragma omp for
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.clear();
                for( Index k = 0; k < _planes; ++k )
                {
                    if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
                    {
                        work.push_back( ( pixel( i, j ) )( k, 0 ) );
                    }
                }
 
                if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
                {
                    mim( i, j ) = math::vectorSigmaMean( work, sigma );
                }
                else
                {
                    mim( i, j ) = invalidNumber<Scalar>();
                }
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT>
void eigenCube<dataT>::sigmaMean( eigenT &mim, std::vector<dataT> &weights, dataT sigma )
{
    validateWeights( weights );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel num_threads( Eigen::nbThreads() )
    {
        std::vector<Scalar> work;
 
#pragma omp for schedule( static, 10 )
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.resize( _planes ); // work could be smaller after sigmaMean
                for( int k = 0; k < _planes; ++k )
                    work[k] = ( pixel( i, j ) )( k, 0 );
 
                mim( i, j ) = math::vectorSigmaMean( work, weights, sigma );
            }
        }
    }
}
 
template <typename dataT>
template <typename eigenT, typename eigenCubeT>
void eigenCube<dataT>::sigmaMean(
    eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, dataT sigma, double minGoodFract )
{
    validateMaskedCombination( mask, minGoodFract );
    validateWeights( weights );
 
    mim.resize( _rows, _cols );
 
#pragma omp parallel num_threads( Eigen::nbThreads() )
    {
        std::vector<Scalar> work, wwork;
 
#pragma omp for schedule( static, 10 )
        for( Index i = 0; i < _rows; ++i )
        {
            for( Index j = 0; j < _cols; ++j )
            {
                work.clear();
                wwork.clear();
 
                for( Index k = 0; k < _planes; ++k )
                {
                    if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
                    {
                        work.push_back( ( pixel( i, j ) )( k, 0 ) );
                        wwork.push_back( weights[k] );
                    }
                }
                if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
                {
                    mim( i, j ) = math::vectorSigmaMean( work, wwork, sigma );
                }
                else
                    mim( i, j ) = invalidNumber<Scalar>();
            }
        }
    }
}
 
} // namespace improc
 
} // namespace mx
 
#endif // eigenCube_hpp