/scr/tee1/isis/lib/Core/DataStorage/image.hpp

Go to the documentation of this file.

//
// C++ Interface: image
//
// Description:
//
//
// Author: Enrico Reimer<reimer@cbs.mpg.de>, (C) 2009
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef IMAGE_H
#define IMAGE_H

#include "chunk.hpp"

#include <set>
#include <boost/shared_ptr.hpp>
#include <vector>
#include <boost/foreach.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <stack>
#include "sortedchunklist.hpp"
#include "common.hpp"

namespace isis
{
namespace data
{
namespace _internal
{
template<typename CHUNK_TYPE> class ImageIteratorTemplate: public std::iterator <
    std::random_access_iterator_tag,
    typename boost::mpl::if_<boost::is_const<CHUNK_TYPE>, typename CHUNK_TYPE::const_iterator, typename CHUNK_TYPE::iterator>::type::value_type,
    typename boost::mpl::if_<boost::is_const<CHUNK_TYPE>, typename CHUNK_TYPE::const_iterator, typename CHUNK_TYPE::iterator>::type::difference_type,
    typename boost::mpl::if_<boost::is_const<CHUNK_TYPE>, typename CHUNK_TYPE::const_iterator, typename CHUNK_TYPE::iterator>::type::pointer,
    typename boost::mpl::if_<boost::is_const<CHUNK_TYPE>, typename CHUNK_TYPE::const_iterator, typename CHUNK_TYPE::iterator>::type::reference
    >
{
protected:
    typedef typename boost::mpl::if_<boost::is_const<CHUNK_TYPE>, typename CHUNK_TYPE::const_iterator, typename CHUNK_TYPE::iterator>::type inner_iterator;
    typedef CHUNK_TYPE chunk_type;
    typedef ImageIteratorTemplate<CHUNK_TYPE> ThisType;

    std::vector<chunk_type *> chunks;
    size_t ch_idx;
    inner_iterator current_it;
    typename inner_iterator::difference_type ch_len;

    typename inner_iterator::difference_type currentDist() const {
        if ( ch_idx >= chunks.size() )
            return 0; // if we're behind the last chunk assume we are at the "start" of the "end"-chunk
        else {
            const inner_iterator chit_begin = chunks[ch_idx]->begin(); // cast in a const or cast out a non existing one
            return std::distance ( chit_begin, current_it ); // so we use same iterators here
        }
    }
    friend class ImageIteratorTemplate<const CHUNK_TYPE>; //yes, I'm my own friend, sometimes :-) (enables the constructor below)
public:

    //will become additional constructor from non const if this is const, otherwise overrride the default copy contructor
    ImageIteratorTemplate ( const ImageIteratorTemplate<typename boost::remove_const<CHUNK_TYPE>::type > &src ) :
        chunks ( src.chunks.begin(), src.chunks.end() ), ch_idx ( src.ch_idx ),
        current_it ( src.current_it ),
        ch_len ( src.ch_len )
    {}

    // empty constructor
    ImageIteratorTemplate() : ch_idx ( 0 ), ch_len ( 0 ) {}


    // normal conytructor
    explicit ImageIteratorTemplate ( const std::vector<chunk_type *>& _chunks ) :
        chunks ( _chunks ), ch_idx ( 0 ),
        current_it ( chunks[0]->begin() ),
        ch_len ( std::distance ( current_it, chunks[0]->end() ) )
    {}

    ThisType &operator++() {
        return operator+= ( 1 );
    }
    ThisType &operator--() {
        return operator-= ( 1 );
    }

    ThisType operator++ ( int ) {
        ThisType tmp = *this;
        operator++();
        return tmp;
    }
    ThisType operator-- ( int ) {
        ThisType tmp = *this;
        operator--();
        return tmp;
    }

    typename inner_iterator::reference operator*() const {
        return current_it.operator * ();
    }
    typename inner_iterator::pointer  operator->() const {
        return current_it.operator->();
    }

    bool operator== ( const ThisType &cmp ) const {
        return ch_idx == cmp.ch_idx && current_it == cmp.current_it;
    }
    bool operator!= ( const ThisType &cmp ) const {
        return !operator== ( cmp );
    }

    bool operator> ( const ThisType &cmp ) const {
        return ch_idx > cmp.ch_idx || ( ch_idx == cmp.ch_idx && current_it > cmp.current_it );
    }
    bool operator< ( const ThisType &cmp ) const {
        return ch_idx < cmp.ch_idx || ( ch_idx == cmp.ch_idx && current_it < cmp.current_it );
    }


    bool operator>= ( const ThisType &cmp ) const {
        return operator> ( cmp ) || operator== ( cmp );
    }
    bool operator<= ( const ThisType &cmp ) const {
        return operator< ( cmp ) || operator== ( cmp );
    }

    typename inner_iterator::difference_type operator- ( const ThisType &cmp ) const {
        typename inner_iterator::difference_type dist = ( ch_idx - cmp.ch_idx ) * ch_len; // get the (virtual) distance from my current block to cmp's current block

        if ( ch_idx >= cmp.ch_idx ) { //if I'm beyond cmp add my current pos to the distance, and substract his
            dist += currentDist() - cmp.currentDist();
        } else {
            dist += cmp.currentDist() - currentDist();
        }

        return dist;
    }

    ThisType operator+ ( typename ThisType::difference_type n ) const {
        return ThisType ( *this ) += n;
    }
    ThisType operator- ( typename ThisType::difference_type n ) const {
        return ThisType ( *this ) -= n;
    }

    ThisType &operator+= ( typename inner_iterator::difference_type n ) {
        n += currentDist(); //start from current begin (add current_it-(begin of the current chunk) to n)
        assert ( ( n / ch_len + static_cast<typename ThisType::difference_type> ( ch_idx ) ) >= 0 );
        ch_idx += n / ch_len; //if neccesary jump to next chunk

        if ( ch_idx < chunks.size() )
            current_it = chunks[ch_idx]->begin() + n % ch_len; //set new current iterator in new chunk plus the "rest"
        else
            current_it = ( * ( chunks.end() - 1 ) )->end() ; //set current_it to the last chunks end iterator if we are behind it

        return *this;
    }
    ThisType &operator-= ( typename inner_iterator::difference_type n ) {
        return operator+= ( -n );
    }

    typename ThisType::reference operator[] ( typename inner_iterator::difference_type n ) const {
        return * ( ThisType ( chunks ) += n );
    }

};
}

class ChunkOp : std::unary_function<Chunk &, bool>
{
public:
    virtual bool operator() ( Chunk &, util::vector4<size_t> posInImage ) = 0;
    virtual ~ChunkOp();
};

class Image:
    public _internal::NDimensional<4>,
    public util::PropertyMap
{
    dimensions minIndexingDim;
public:
    void setIndexingDim ( dimensions d = rowDim );
    enum orientation {axial, reversed_axial, sagittal, reversed_sagittal, coronal, reversed_coronal};

    typedef _internal::ImageIteratorTemplate<Chunk> iterator;
    typedef _internal::ImageIteratorTemplate<const Chunk> const_iterator;
    typedef iterator::reference reference;
    typedef const_iterator::reference const_reference;
    static const char *neededProperties;
protected:
    _internal::SortedChunkList set;
    std::vector<boost::shared_ptr<Chunk> > lookup;
private:
    size_t chunkVolume;

    void deduplicateProperties();

    const boost::shared_ptr<Chunk> &chunkPtrAt ( size_t at ) const;

    inline std::pair<size_t, size_t> commonGet ( size_t first, size_t second, size_t third, size_t fourth ) const {
        const size_t idx[] = {first, second, third, fourth};
        LOG_IF ( ! clean, Debug, error )
                << "Getting data from a non indexed image will result in undefined behavior. Run reIndex first.";
        LOG_IF ( set.isEmpty(), Debug, error )
                << "Getting data from a empty image will result in undefined behavior.";
        LOG_IF ( !isInRange ( idx ), Debug, isis::error )
                << "Index " << util::vector4<size_t> ( idx ) << " is out of range (" << getSizeAsString() << ")";
        const size_t index = getLinearIndex ( idx );
        return std::make_pair ( index / chunkVolume, index % chunkVolume );
    }


protected:
    bool clean;
    static const char *defaultChunkEqualitySet;

    size_t getChunkStride ( size_t base_stride = 1 );
    Chunk &chunkAt ( size_t at );
    Image();

    util::fvector3 m_RowVec;
    util::fvector3 m_RowVecInv;
    util::fvector3 m_ColumnVec;
    util::fvector3 m_ColumnVecInv;
    util::fvector3 m_SliceVec;
    util::fvector3 m_SliceVecInv;
    util::fvector3 m_Offset;

public:
    Image ( const Image &ref );

    template<typename T> Image ( std::list<T> &chunks, dimensions min_dim = rowDim ) :
        _internal::NDimensional<4>(), util::PropertyMap(), minIndexingDim ( min_dim ),
        set ( defaultChunkEqualitySet ),
        clean ( false ) {
        util::Singletons::get<NeededsList<Image>, 0>().applyTo( *this );
        set.addSecondarySort ( "acquisitionNumber" );
        insertChunksFromList ( chunks );
    }
    template<typename T> Image ( std::vector<T> &chunks, dimensions min_dim = rowDim ) :
        _internal::NDimensional<4>(), util::PropertyMap(), minIndexingDim ( min_dim ),
        set ( defaultChunkEqualitySet ),
        clean ( false ) {
        util::Singletons::get<NeededsList<Image>, 0>().applyTo( *this );
        set.addSecondarySort ( "acquisitionNumber" );
        std::list<T> tmp( chunks.begin(), chunks.end() );
        insertChunksFromList ( tmp );
        chunks.assign( tmp.begin(), tmp.end() );
    }

    template<typename T> size_t insertChunksFromList ( std::list<T> &chunks ) {
        BOOST_MPL_ASSERT ( ( boost::is_base_of<Chunk, T> ) );
        size_t cnt = 0;

        for ( typename std::list<T>::iterator i = chunks.begin(); i != chunks.end(); ) { // for all remaining chunks
            if ( insertChunk ( *i ) ) {
                chunks.erase ( i++ );
                cnt++;
            } else {
                i++;
            }
        }

        if ( ! isEmpty() ) {
            LOG ( Debug, info ) << "Reindexing image with " << cnt << " chunks.";

            if ( !reIndex() ) {
                LOG ( Runtime, error ) << "Failed to create image from " << cnt << " chunks.";
            } else {
                LOG_IF ( !getMissing().empty(), Debug, warning )
                        << "The created image is missing some properties: " << getMissing() << ". It will be invalid.";
            }
        } else {
            LOG ( Debug, warning ) << "Image is empty after inserting chunks.";
        }

        return cnt;
    }


    Image ( const Chunk &chunk, dimensions min_dim = rowDim );

    Image &operator= ( const Image &ref );

    bool checkMakeClean();
    bool isClean() const;
    template <typename T> T &voxel ( size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0 ) {
        checkMakeClean();
        const std::pair<size_t, size_t> index = commonGet ( first, second, third, fourth );
        ValueArray<T> &data = chunkAt ( index.first ).asValueArray<T>();
        return data[index.second];
    }

    template <typename T> const T &voxel ( size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0 ) const {
        const std::pair<size_t, size_t> index = commonGet ( first, second, third, fourth );
        const ValueArray<T> &data = chunkPtrAt ( index.first )->getValueArray<T>();
        return data[index.second];
    }

    const util::ValueReference getVoxelValue ( size_t nrOfColumns, size_t nrOfRows = 0, size_t nrOfSlices = 0, size_t nrOfTimesteps = 0 ) const;
    void setVoxelValue ( const util::ValueReference &val, size_t nrOfColumns, size_t nrOfRows = 0, size_t nrOfSlices = 0, size_t nrOfTimesteps = 0 );

    unsigned short getMajorTypeID() const;
    std::string getMajorTypeName() const;

    iterator begin();
    iterator end();
    const_iterator begin() const;
    const_iterator end() const;

    Chunk getChunkAt ( size_t at, bool copy_metadata = true ) const;

    const Chunk getChunk ( size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0, bool copy_metadata = true ) const;

    Chunk getChunk ( size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0, bool copy_metadata = true );

    template<typename TYPE> Chunk getChunkAs ( size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0, bool copy_metadata = true ) const {
        return getChunkAs<TYPE> ( getScalingTo ( ValueArray<TYPE>::staticID ), first, second, third, fourth, copy_metadata );
    }
    template<typename TYPE> Chunk getChunkAs ( const scaling_pair &scaling, size_t first, size_t second = 0, size_t third = 0, size_t fourth = 0, bool copy_metadata = true ) const {
        Chunk ret = getChunk ( first, second, third, fourth, copy_metadata ); // get a cheap copy
        ret.convertToType ( ValueArray<TYPE>::staticID, scaling ); // make it of type T
        return ret; //return that
    }

    scaling_pair getScalingTo ( unsigned short typeID, autoscaleOption scaleopt = autoscale ) const;


    bool insertChunk ( const Chunk &chunk );
    bool reIndex();

    bool isEmpty() const;

    std::list<util::PropertyValue> getChunksProperties ( const util::PropertyMap::KeyType &key, bool unique = false ) const;

    size_t getMaxBytesPerVoxel() const;

    template<typename T> std::pair<T, T> getMinMaxAs() const {
        util::checkType<T>();// works only for T from _internal::types
        std::pair<util::ValueReference, util::ValueReference> minmax = getMinMax();
        return std::make_pair ( minmax.first->as<T>(), minmax.second->as<T>() );
    }

    std::pair<util::ValueReference, util::ValueReference> getMinMax() const;

    size_t compare ( const Image &comp ) const;

    orientation getMainOrientation() const;

    bool transformCoords ( boost::numeric::ublas::matrix<float> transform_matrix, bool transformCenterIsImageCenter = false );

    dimensions mapScannerAxisToImageDimension ( scannerAxis scannerAxes );

    util::fvector3 getPhysicalCoordsFromIndex ( const util::ivector4 &index ) const;


    util::ivector4 getIndexFromPhysicalCoords ( const util::fvector3 &physicalCoords ) const;

    template<typename T> void copyToMem ( T *dst, size_t len,  scaling_pair scaling = scaling_pair() ) const {
        if ( clean ) {
            if ( scaling.first.isEmpty() || scaling.second.isEmpty() ) {
                scaling = getScalingTo ( ValueArray<T>::staticID );
            }

            // we could do this using convertToType - but this solution does not need any additional temporary memory
            BOOST_FOREACH ( const boost::shared_ptr<Chunk> &ref, lookup ) {
                const size_t cSize = ref->getSizeAsVector().product();

                if ( !ref->copyToMem<T> ( dst, len, scaling ) ) {
                    LOG ( Runtime, error ) << "Failed to copy raw data of type " << ref->getTypeName() << " from image into memory of type " << ValueArray<T>::staticName();
                } else {
                    if ( len < cSize ) {
                        LOG ( Runtime, error ) << "Abborting copy, because there is no space left in the target";
                        break;
                    }

                    len -= cSize;
                }

                dst += ref->getVolume(); // increment the cursor
            }
        } else {
            LOG ( Runtime, error ) << "Cannot copy from non clean images. Run reIndex first";
        }
    }

    template<typename T> MemChunk<T> copyAsMemChunk() const {
        const util::vector4<size_t> size = getSizeAsVector();
        data::MemChunk<T> ret ( size[0], size[1], size[2], size[3] );
        copyToMem<T> ( &ret.voxel<T>(0,0), ret.getVolume() );
        static_cast<util::PropertyMap &>( ret ) = static_cast<const util::PropertyMap &>( getChunkAt( 0 ) );
        return ret;
    }

    template<typename T> ValueArray<T> copyAsValueArray() const {
        data::ValueArray<T> ret ( getVolume() );
        copyToMem<T> ( ret.begin().operator->(), ret.getLength() );
        return ret;
    }
    void copyToValueArray ( data::ValueArrayBase &dst,  scaling_pair scaling = scaling_pair() ) const;

    Image copyByID( unsigned short ID = 0, scaling_pair scaling = scaling_pair() )const;


    std::vector<isis::data::Chunk> copyChunksToVector ( bool copy_metadata = true ) const;

    bool convertToType ( short unsigned int ID, isis::data::autoscaleOption scaleopt = autoscale );

    size_t spliceDownTo ( dimensions dim );

    size_t foreachChunk ( ChunkOp &op, bool copyMetaData = false );


    template <typename TYPE> size_t foreachVoxel ( VoxelOp<TYPE> &op ) {
        class _proxy: public ChunkOp
        {
            VoxelOp<TYPE> &op;
        public:
            _proxy ( VoxelOp<TYPE> &_op ) : op ( _op ) {}
            bool operator() ( Chunk &ch, util::vector4<size_t> posInImage ) {
                return ch.foreachVoxel<TYPE> ( op, posInImage ) == 0;
            }
        };
        _proxy prx ( op );
        return convertToType ( data::ValueArray<TYPE>::staticID ) && foreachChunk ( prx, false );
    }

    size_t getNrOfRows() const;
    size_t getNrOfColumns() const;
    size_t getNrOfSlices() const;
    size_t getNrOfTimesteps() const;

    util::fvector3 getFoV() const;
    bool updateOrientationMatrices();

    std::string identify( bool withpath = true )const;
};

template<typename T> class TypedImage: public Image
{
protected:
    TypedImage() {} // to be used only by inheriting classes
public:
    typedef _internal::ImageIteratorTemplate<data::ValueArray<T> > iterator;
    typedef _internal::ImageIteratorTemplate<const data::ValueArray<T> > const_iterator;
    typedef typename iterator::reference reference;
    typedef typename const_iterator::reference const_reference;
    TypedImage ( const Image &src ) : Image ( src ) { // ok we just copied the whole image
        //but we want it to be of type T
        convertToType ( ValueArray<T>::staticID );
    }
    TypedImage &operator= ( const TypedImage &ref ) { //its already of the given type - so just copy it
        Image::operator= ( ref );
        return *this;
    }
    TypedImage &operator= ( const Image &ref ) { // copy the image, and make sure its of the given type
        Image::operator= ( ref );
        convertToType ( ValueArray<T>::staticID );
        return *this;
    }
    void copyToMem ( void *dst ) {
        Image::copyToMem<T> ( ( T * ) dst );
    }
    void copyToMem ( void *dst ) const {
        Image::copyToMem<T> ( ( T * ) dst );
    }
    iterator begin() {
        if ( checkMakeClean() ) {
            std::vector<data::ValueArray<T>*> vec ( lookup.size() );

            for ( size_t i = 0; i < lookup.size(); i++ )
                vec[i] = &lookup[i]->template asValueArray<T>();

            return iterator ( vec );
        } else {
            LOG ( Debug, error )  << "Image is not clean. Returning empty iterator ...";
            return iterator();
        }
    }
    iterator end() {
        return begin() + getVolume();
    };
    const_iterator begin() const {
        if ( isClean() ) {
            std::vector<const data::ValueArray<T>*> vec ( lookup.size() );

            for ( size_t i = 0; i < lookup.size(); i++ )
                vec[i] = &lookup[i]->template asValueArray<T>();

            return const_iterator ( vec );
        } else {
            LOG ( Debug, error )  << "Image is not clean. Returning empty iterator ...";
            return const_iterator();
        }
    }
    const_iterator end() const {
        return begin() + getVolume();
    };
};

template<typename T> class MemImage: public TypedImage<T>
{
public:
    MemImage ( const Image &src ) {
        operator= ( src );
    }

    MemImage &operator= ( const Image &ref ) { // copy the image, and make sure its of the given type

        Image::operator= ( ref ); // ok we just copied the whole image

        //we want deep copies of the chunks, and we want them to be of type T
        struct : _internal::SortedChunkList::chunkPtrOperator {
            std::pair<util::ValueReference, util::ValueReference> scale;
            boost::shared_ptr<Chunk> operator() ( const boost::shared_ptr< Chunk >& ptr ) {
                return boost::shared_ptr<Chunk> ( new MemChunk<T> ( *ptr, scale ) );
            }
        } conv_op;
        conv_op.scale = ref.getScalingTo ( ValueArray<T>::staticID );
        LOG ( Debug, info ) << "Computed scaling for conversion from source image: [" << conv_op.scale << "]";

        this->set.transform ( conv_op );

        if ( ref.isClean() ) {
            this->lookup = this->set.getLookup(); // the lookup table still points to the old chunks
        } else {
            LOG ( Debug, info ) << "Copied unclean image. Running reIndex on the copy.";
            this->reIndex();
        }

        return *this;
    }
};

}
}

#endif // IMAGE_H