imagingArray.hpp

Go to the documentation of this file.

/** \file imagingArray.hpp
 * \brief Declares and defines a class for managing images
 * \ingroup imaging
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
#ifndef imagingArray_hpp
#define imagingArray_hpp
 
#pragma GCC system_header
#include <Eigen/Dense>
 
#include "../math/templateBLAS.hpp"
#include "../math/fft/fft.hpp"
 
namespace mx
{
namespace wfp
{
 
template <typename Scalar>
struct fftwAllocator
{
    Scalar *alloc( size_t sz )
    {
        return (Scalar *)::fftw_malloc( sz * sizeof( Scalar ) );
    }
 
    void free( Scalar *ptr )
    {
        fftw_free( ptr );
    }
};
 
/// Test whether a type is an imagingArray by testing whether it has a typedef of "is_imagingArray"
/** Used for compile-time determination of type
 * Example usage:
 * \code
 * bool is_eC = is_imagingArray<eigenCube<float> >; //Evaluates to true
 * bool is_not_eC = is_imagingArray<float>; //Evaluates to false
 * \endcode
 */
// This was taken directly from the example at http://en.wikipedia.org/wiki/Substitution_failure_is_not_an_error
template <typename T>
struct is_imagingArray
{
    // Types "yes" and "no" are guaranteed to have different sizes,
    // specifically sizeof(yes) == 1 and sizeof(no) == 2.
    typedef char yes[1];
    typedef char no[2];
 
    template <typename imageT>
    static yes &test( typename imageT::imagingArrayT * );
 
    template <typename>
    static no &test( ... );
 
    // If the "sizeof" of the result of calling test<T>(0) would be equal to sizeof(yes),
    // the first overload worked and T has a nested type named "is_mmatrix".
    static const bool value = sizeof( test<T>( 0 ) ) == sizeof( yes );
};
 
template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
struct imagingArrayInplaceProduct
{
    imagingArrayInplaceProduct()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceProduct template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT &im, typeT &alpha )
    {
        math::scal<typename imageT::Scalar>( im.cols() * im.rows(), (typeT)alpha, im.data(), 1 );
    }
};
 
template <typename imageT, typename typeT>
struct imagingArrayInplaceProduct<imageT, typeT, true>
{
    imagingArrayInplaceProduct()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceProduct template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT &im, typeT &alpha )
    {
        typedef typename imageT::Scalar scalarT;
        // hadp<scalarT>(im.cols()*im.rows(), alpha.data(), im.data());
 
        Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigY( alpha.data(), alpha.cols(), alpha.rows() );
        Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
 
        eigX *= eigY;
    }
};
 
template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
struct imagingArrayInplaceDivision
{
    imagingArrayInplaceDivision()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceDivision template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT im, typeT alpha )
    {
        math::scal<typename imageT::Scalar>( im.cols() * im.rows(), ( (typeT)( 1 ) ) / alpha, im.data(), 1 );
    }
};
 
template <typename imageT, typename typeT>
struct imagingArrayInplaceDivision<imageT, typeT, true>
{
    imagingArrayInplaceDivision()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceDivision template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT im, typeT alpha )
    {
        typedef typename imageT::Scalar scalarT;
        math::hadd<scalarT>( im.cols() * im.rows(), (scalarT)1.0, alpha.data(), 1, im.data(), 1 );
    }
};
 
template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
struct imagingArrayInplaceAdd
{
    imagingArrayInplaceAdd()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceAdd template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT &im, typeT &alpha )
    {
        typedef typename imageT::Scalar scalarT;
        Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
        eigX += alpha;
    }
};
 
template <typename imageT, typename typeT>
struct imagingArrayInplaceAdd<imageT, typeT, true>
{
    imagingArrayInplaceAdd()
    {
        static_assert( is_imagingArray<imageT>::value,
                       "imagingArrayInplaceAdd template parameter 1 must be an imagingArray type" );
    }
 
    void operator()( imageT &im, typeT &alpha )
    {
        typedef typename imageT::Scalar scalarT;
        // hadp<scalarT>(im.cols()*im.rows(), alpha.data(), im.data());
 
        Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigY( alpha.data(), alpha.cols(), alpha.rows() );
        Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
 
        eigX += eigY;
    }
};
 
template <typename _Scalar, class _allocatorT, int cudaGPU>
class imagingArray;
 
template <typename _Scalar, class _allocatorT>
class imagingArray<_Scalar, _allocatorT, 0>
{
  public:
    typedef bool imagingArrayT;
    typedef _Scalar Scalar;
    typedef _allocatorT allocatorT;
 
  protected:
    allocatorT allocator;
 
    Scalar *_data;
 
    bool _owner;
 
    int _cols;
    int _rows;
 
  public:
    imagingArray()
    {
        _data = 0;
        _owner = false;
        _cols = 0;
        _rows = 0;
    }
 
    imagingArray( int x, int y )
    {
        _data = 0;
        _owner = false;
        _cols = 0;
        _rows = 0;
 
        resize( x, y );
    }
 
    template <typename eigenT>
    explicit imagingArray( eigenT &im )
    {
        _data = im.data();
        _owner = false;
        _cols = im.cols();
        _rows = im.rows();
    }
 
    ~imagingArray()
    {
        free();
    }
 
    void resize( int szx, int szy )
    {
        if( szx == _cols && szy == _rows && _owner )
            return;
 
        free();
 
        _cols = szx;
        _rows = szy;
 
        _data = allocator.alloc( _cols * _rows );
        _owner = true;
    }
 
    void free()
    {
        if( _data && _owner )
        {
            allocator.free( _data );
            _data = 0;
            _cols = 0;
            _rows = 0;
        }
    }
 
    int cols() const
    {
        return _cols;
    }
 
    int rows() const
    {
        return _rows;
    }
 
    Scalar *data()
    {
        return _data;
    }
 
    Scalar *data() const
    {
        return _data;
    }
 
    Scalar operator()( int x, int y ) const
    {
        return _data[x + y * _cols];
    }
 
    Scalar &operator()( int x, int y )
    {
        return _data[x + y * _cols];
    }
 
    template <typename argT>
    void set( argT val )
    {
        int N = _cols * _rows;
 
        Scalar dval = (Scalar)val;
        for( int i = 0; i < N; ++i )
            _data[i] = dval;
    }
 
    void setZero()
    {
        set( (Scalar)0 );
    }
 
    imagingArray &operator=( const imagingArray &arr )
    {
        resize( arr.rows(), arr.cols() );
 
        for( int i = 0; i < _rows; ++i )
        {
            for( int j = 0; j < _cols; ++j )
            {
                operator()( i, j ) = arr( i, j );
            }
        }
 
        return ( *this );
    }
 
    /// In-place product.  Works for both scalars and imagingArray types.
    template <typename typeT>
    imagingArray &operator*=( typeT &alpha )
    {
        imagingArrayInplaceProduct<imagingArray, typeT> prod;
        prod( *this, alpha );
        return *this;
    }
 
    /// In-place division.  Works for both scalars and imagingArray types.
    template <typename typeT>
    imagingArray &operator/=( typeT alpha )
    {
        imagingArrayInplaceDivision<imagingArray, typeT> div;
        div( *this, alpha );
        return *this;
    }
 
    /// In-place add.  Works for both scalars and imagingArray types.
    template <typename typeT>
    imagingArray &operator+=( typeT &alpha )
    {
        imagingArrayInplaceAdd<imagingArray, typeT> prod;
        prod( *this, alpha );
        return *this;
    }
 
    Scalar sum()
    {
        int N = _cols * _rows;
        Scalar _sum = 0;
 
        for( int i = 0; i < N; ++i )
            _sum += _data[i];
 
        return _sum;
    }
 
    Eigen::Map<Eigen::Array<Scalar, -1, -1>> eigenMap()
    {
        return Eigen::Map<Eigen::Array<Scalar, -1, -1>>( _data, _rows, _cols );
    }
};
 
template <typename scalarT>
using imagingArrayT = imagingArray<scalarT, fftwAllocator<scalarT>, 0>;
 
// typedef imagingArray<float,mx::fftwAllocator<float>,0> imagingArrayf;
// typedef imagingArray<double,mx::fftwAllocator<double>,0> imagingArrayd;
// typedef imagingArray<std::complex<float>,mx::fftwAllocator<std::complex<float> >,0> imagingArraycf;
 
} // namespace wfp
} // namespace mx
 
#endif // imagingArray_hpp