NiftiImage.h

#ifndef _NIFTI_IMAGE_H_
#define _NIFTI_IMAGE_H_


#ifdef USING_R

#include <Rcpp.h>

// Defined since R 3.1.0, according to Tomas Kalibera, but there's no reason to break compatibility with 3.0.x
#ifndef MAYBE_SHARED
#define MAYBE_SHARED(x) (NAMED(x) > 1)
#endif

#else

#define R_NegInf -INFINITY

#include <stdint.h>
#include <cstddef>
#include <cmath>
#include <string>
#include <sstream>
#include <vector>
#include <complex>
#include <stdexcept>
#include <algorithm>
#include <map>
#include <locale>
#include <limits>

#endif


#include "niftilib/nifti1_io.h"

namespace RNifti {

typedef std::complex<float> complex64_t;
typedef std::complex<double> complex128_t;

struct rgba32_t
{
    union ValueType {
        int packed;
        unsigned char bytes[4];
    };
    ValueType value;
    rgba32_t () { value.packed = 0; }
};

class NiftiImageData
{
public:
    double slope;                       
    double intercept;                   
protected:
    struct TypeHandler
    {
        virtual ~TypeHandler() {}
        virtual size_t size () const { return 0; }
        virtual bool hasNaN () const { return false; }
        virtual complex128_t getComplex (void *ptr) const { return complex128_t(0.0, 0.0); }
        virtual double getDouble (void *ptr) const { return 0.0; }
        virtual int getInt (void *ptr) const { return 0; }
        virtual rgba32_t getRgb (void *ptr) const { return rgba32_t(); }
        virtual void setComplex (void *ptr, const complex128_t value) const {}
        virtual void setDouble (void *ptr, const double value) const {}
        virtual void setInt (void *ptr, const int value) const {}
        virtual void setRgb (void *ptr, const rgba32_t value) const {}
        virtual void minmax (void *ptr, const size_t length, double *min, double *max) const { *min = 0.0; *max = 0.0; }
    };
    
    template <typename Type, bool alpha = false>
    struct ConcreteTypeHandler : public TypeHandler
    {
        size_t size () const { return (sizeof(Type)); }
        bool hasNaN () const { return std::numeric_limits<Type>::has_quiet_NaN; }
        complex128_t getComplex (void *ptr) const { return complex128_t(static_cast<double>(*static_cast<Type*>(ptr)), 0.0); }
        double getDouble (void *ptr) const { return static_cast<double>(*static_cast<Type*>(ptr)); }
        int getInt (void *ptr) const { return static_cast<int>(*static_cast<Type*>(ptr)); }
        void setComplex (void *ptr, const complex128_t value) const
        {
            *(static_cast<Type*>(ptr)) = Type(value.real());
            *(static_cast<Type*>(ptr) + 1) = Type(0);
        }
        void setDouble (void *ptr, const double value) const { *(static_cast<Type*>(ptr)) = Type(value); }
        void setInt (void *ptr, const int value) const { *(static_cast<Type*>(ptr)) = Type(value); }
        void minmax (void *ptr, const size_t length, double *min, double *max) const;
    };
    
    template <typename ElementType>
    struct ConcreteTypeHandler<std::complex<ElementType>,false> : public TypeHandler
    {
        size_t size () const { return (sizeof(ElementType) * 2); }
        bool hasNaN () const { return std::numeric_limits<ElementType>::has_quiet_NaN; }
        std::complex<ElementType> getNative (void *ptr) const
        {
            const ElementType real = *static_cast<ElementType*>(ptr);
            const ElementType imag = *(static_cast<ElementType*>(ptr) + 1);
            return std::complex<ElementType>(real, imag);
        }
        void setNative (void *ptr, const std::complex<ElementType> native) const
        {
            *(static_cast<ElementType*>(ptr)) = native.real();
            *(static_cast<ElementType*>(ptr) + 1) = native.imag();
        }
        complex128_t getComplex (void *ptr) const { return complex128_t(getNative(ptr)); }
        double getDouble (void *ptr) const { return static_cast<double>(getNative(ptr).real()); }
        int getInt (void *ptr) const { return static_cast<int>(getNative(ptr).real()); }
        void setComplex (void *ptr, const complex128_t value) const { setNative(ptr, std::complex<ElementType>(value)); }
        void setDouble (void *ptr, const double value) const { setNative(ptr, std::complex<ElementType>(value, 0.0)); }
        void setInt (void *ptr, const int value) const { setNative(ptr, std::complex<ElementType>(static_cast<ElementType>(value), 0.0)); }
        void minmax (void *ptr, const size_t length, double *min, double *max) const;
    };
    
    template <bool alpha>
    struct ConcreteTypeHandler<rgba32_t,alpha> : public TypeHandler
    {
        size_t size () const { return alpha ? 4 : 3; }
        int getInt (void *ptr) const { return getRgb(ptr).value.packed; }
        rgba32_t getRgb (void *ptr) const
        {
            rgba32_t value;
            unsigned char *source = static_cast<unsigned char *>(ptr);
            std::copy(source, source + (alpha ? 4 : 3), value.value.bytes);
            return value;
        }
        void setInt (void *ptr, const int value) const
        {
            rgba32_t native;
            native.value.packed = value;
            setRgb(ptr, native);
        }
        void setRgb (void *ptr, const rgba32_t value) const
        {
            unsigned char *target = static_cast<unsigned char *>(ptr);
            std::copy(value.value.bytes, value.value.bytes + (alpha ? 4 : 3), target);
        }
        void minmax (void *ptr, const size_t length, double *min, double *max) const { *min = 0.0; *max = 255.0; }
    };
    
    TypeHandler * createHandler ()
    {
        if (_datatype == DT_NONE)
            return NULL;
        
        switch (_datatype)
        {
            case DT_UINT8:      return new ConcreteTypeHandler<uint8_t>();          break;
            case DT_INT16:      return new ConcreteTypeHandler<int16_t>();          break;
            case DT_INT32:      return new ConcreteTypeHandler<int32_t>();          break;
            case DT_FLOAT32:    return new ConcreteTypeHandler<float>();            break;
            case DT_FLOAT64:    return new ConcreteTypeHandler<double>();           break;
            case DT_INT8:       return new ConcreteTypeHandler<int8_t>();           break;
            case DT_UINT16:     return new ConcreteTypeHandler<uint16_t>();         break;
            case DT_UINT32:     return new ConcreteTypeHandler<uint32_t>();         break;
            case DT_INT64:      return new ConcreteTypeHandler<int64_t>();          break;
            case DT_UINT64:     return new ConcreteTypeHandler<uint64_t>();         break;
            case DT_COMPLEX64:  return new ConcreteTypeHandler<complex64_t>();      break;
            case DT_COMPLEX128: return new ConcreteTypeHandler<complex128_t>();     break;
            case DT_RGB24:      return new ConcreteTypeHandler<rgba32_t,false>();   break;
            case DT_RGBA32:     return new ConcreteTypeHandler<rgba32_t,true>();    break;
            
            default:
            throw std::runtime_error("Unsupported data type (" + std::string(nifti_datatype_string(_datatype)) + ")");
        }
    }
    
    void *dataPtr;                      
    int _datatype;                      
    TypeHandler *handler;               
    size_t _length;                     
    bool owner;                         
    void init (void *data, const size_t length, const int datatype, const double slope, const double intercept, const bool alloc = true)
    {
        this->_length = length;
        this->_datatype = datatype;
        this->slope = slope;
        this->intercept = intercept;
        
        owner = false;
        handler = createHandler();
        if (handler == NULL)
            dataPtr = NULL;
        else if (alloc && data == NULL)
        {
            dataPtr = calloc(length, handler->size());
            owner = true;
        }
        else
            dataPtr = data;
    }
    
    void calibrateFrom (const NiftiImageData &data)
    {
        slope = 1.0;
        intercept = 0.0;
        
        if (this->isInteger())
        {
            double dataMin, dataMax, typeMin, typeMax;
            data.minmax(&dataMin, &dataMax);
            handler->minmax(NULL, 0, &typeMin, &typeMax);
            
            // If the source type is floating-point but values are in range, we will just round them
            if (dataMin < typeMin || dataMax > typeMax)
            {
                slope = (dataMax - dataMin) / (typeMax - typeMin);
                intercept = dataMin - (slope) * typeMin;
            }
        }
    }
    
public:
    struct Element
    {
    private:
        const NiftiImageData &parent;
        void *ptr;
        
    public:
        Element (const NiftiImageData &parent, void *ptr = NULL)
            : parent(parent)
        {
            this->ptr = (ptr == NULL ? parent.dataPtr : ptr);
        }
        
        template <typename SourceType>
        Element & operator= (const SourceType &value);
        
        Element & operator= (const Element &other);
        
        template <typename TargetType>
        operator TargetType() const
        {
            if (parent.isScaled())
                return TargetType(parent.handler->getDouble(ptr) * parent.slope + parent.intercept);
            else if (std::numeric_limits<TargetType>::is_integer)
                return TargetType(parent.handler->getInt(ptr));
            else
                return TargetType(parent.handler->getDouble(ptr));
        }
        
        template <typename ElementType>
        operator std::complex<ElementType>() const
        {
            if (parent.isScaled())
                return std::complex<ElementType>(parent.handler->getComplex(ptr) * parent.slope + complex128_t(parent.intercept, parent.intercept));
            else
                return std::complex<ElementType>(parent.handler->getComplex(ptr));
        }
        
#ifdef USING_R

        operator Rcomplex() const
        {
            const complex128_t value = parent.handler->getComplex(ptr);
            Rcomplex rValue = { value.real(), value.imag() };
            if (parent.isScaled())
            {
                rValue.r = rValue.r * parent.slope + parent.intercept;
                rValue.i = rValue.i * parent.slope + parent.intercept;
            }
            return rValue;
        }
#endif
        
        operator rgba32_t() const
        {
            return parent.handler->getRgb(ptr);
        }
    };
    
    class Iterator : public std::iterator<std::random_access_iterator_tag, Element>
    {
    private:
        const NiftiImageData &parent;
        void *ptr;
        size_t step;
        
    public:
        Iterator (const NiftiImageData &parent, void *ptr = NULL, const size_t step = 0)
            : parent(parent)
        {
            this->ptr = (ptr == NULL ? parent.dataPtr : ptr);
            this->step = (step == 0 ? parent.handler->size() : step);
        }
        
        Iterator (const Iterator &other)
            : parent(other.parent), ptr(other.ptr), step(other.step) {}
        
        Iterator & operator++ () { ptr = static_cast<char*>(ptr) + step; return *this; }
        Iterator operator++ (int) { Iterator copy(*this); ptr = static_cast<char*>(ptr) + step; return copy; }
        Iterator operator+ (ptrdiff_t n) const
        {
            void *newptr = static_cast<char*>(ptr) + (n * step);
            return Iterator(parent, newptr, step);
        }
        Iterator & operator-- () { ptr = static_cast<char*>(ptr) - step; return *this; }
        Iterator operator-- (int) { Iterator copy(*this); ptr = static_cast<char*>(ptr) - step; return copy; }
        Iterator operator- (ptrdiff_t n) const
        {
            void *newptr = static_cast<char*>(ptr) - (n * step);
            return Iterator(parent, newptr, step);
        }
        
        ptrdiff_t operator- (const Iterator &other) const
        {
            const ptrdiff_t difference = static_cast<char*>(ptr) - static_cast<char*>(other.ptr);
            return difference / step;
        }
        
        bool operator== (const Iterator &other) const { return (ptr==other.ptr && step==other.step); }
        bool operator!= (const Iterator &other) const { return (ptr!=other.ptr || step!=other.step); }
        bool operator> (const Iterator &other) const { return (ptr > other.ptr); }
        bool operator< (const Iterator &other) const { return (ptr < other.ptr); }
        
        const Element operator* () const { return Element(parent, ptr); }
        Element operator* () { return Element(parent, ptr); }
        const Element operator[] (const size_t i) const { return Element(parent, static_cast<char*>(ptr) + (i * step)); }
        Element operator[] (const size_t i) { return Element(parent, static_cast<char*>(ptr) + (i * step)); }
    };
    
    NiftiImageData ()
        : slope(1.0), intercept(0.0), dataPtr(NULL), _datatype(DT_NONE), handler(NULL), _length(0), owner(false) {}
    
    NiftiImageData (void *data, const size_t length, const int datatype, const double slope = 1.0, const double intercept = 0.0)
    {
        init(data, length, datatype, slope, intercept);
    }
    
    NiftiImageData (nifti_image *image)
    {
        if (image == NULL)
            init(NULL, 0, DT_NONE, 0.0, 0.0, false);
        else
            init(image->data, image->nvox, image->datatype, static_cast<double>(image->scl_slope), static_cast<double>(image->scl_inter), false);
    }
    
    NiftiImageData (const NiftiImageData &source, const int datatype = DT_NONE)
    {
        init(NULL, source.length(), datatype == DT_NONE ? source.datatype() : datatype, source.slope, source.intercept);
        
        if (datatype == DT_NONE || datatype == source.datatype())
            memcpy(dataPtr, source.dataPtr, source.totalBytes());
        else
        {
            calibrateFrom(source);
            std::copy(source.begin(), source.end(), this->begin());
        }
    }
    
    template <class InputIterator>
    NiftiImageData (InputIterator from, InputIterator to, const int datatype)
    {
        const size_t length = static_cast<size_t>(std::distance(from, to));
        init(NULL, length, datatype, 1.0, 0.0);
        std::copy(from, to, this->begin());
    }
    
    virtual ~NiftiImageData ()
    {
        delete handler;
        if (owner)
            free(dataPtr);
    }
    
    NiftiImageData & operator= (const NiftiImageData &source)
    {
        if (source.dataPtr != NULL)
        {
            // Free the old data, if we allocated it
            if (owner)
                free(dataPtr);
            init(NULL, source.length(), source.datatype(), source.slope, source.intercept);
            memcpy(dataPtr, source.dataPtr, source.totalBytes());
        }
        return *this;
    }
    
    void * blob () const             { return dataPtr; }                
    int datatype () const            { return _datatype; }              
    size_t length () const           { return _length; }                
    size_t size () const             { return _length; }                
    size_t bytesPerPixel () const    { return (handler == NULL ? 0 : handler->size()); }
    
    size_t totalBytes () const       { return _length * bytesPerPixel(); }
    
    bool isEmpty () const            { return (dataPtr == NULL); }
    
    bool isScaled () const           { return (slope != 0.0 && (slope != 1.0 || intercept != 0.0)); }
    
    bool isComplex () const          { return (_datatype == DT_COMPLEX64 || _datatype == DT_COMPLEX128); }
    
    bool isFloatingPoint () const    { return (_datatype == DT_FLOAT32 || _datatype == DT_FLOAT64); }
    
    bool isInteger () const          { return nifti_is_inttype(_datatype); }
    
    bool isRgb () const              { return (_datatype == DT_RGB24 || _datatype == DT_RGBA32); }
    
    NiftiImageData unscaled () const { return NiftiImageData(dataPtr, _length, _datatype); }
    
    NiftiImageData & disown ()       { this->owner = false; return *this; }
    
    const Iterator begin () const { return Iterator(*this); }
    
    const Iterator end () const { return Iterator(*this, static_cast<char*>(dataPtr) + totalBytes()); }
    
    Iterator begin () { return Iterator(*this); }
    
    Iterator end () { return Iterator(*this, static_cast<char*>(dataPtr) + totalBytes()); }
    
    const Element operator[] (const size_t i) const { return Element(*this, static_cast<char*>(dataPtr) + (i * bytesPerPixel())); }
    
    Element operator[] (const size_t i) { return Element(*this, static_cast<char*>(dataPtr) + (i * bytesPerPixel())); }
    
    void minmax (double *min, double *max) const
    {
        if (handler == NULL)
        {
            *min = 0.0;
            *max = 0.0;
        }
        else
            handler->minmax(dataPtr, _length, min, max);
    }
};


// R provides an NaN (NA) value for integers
#ifdef USING_R
template <>
inline bool NiftiImageData::ConcreteTypeHandler<int>::hasNaN () const { return true; }
#endif


class NiftiImage
{
public:
    struct Block
    {
        const NiftiImage &image;        
        const int dimension;            
        const int index;                
        Block (const NiftiImage &image, const int dimension, const int index)
            : image(image), dimension(dimension), index(index)
        {
            if (dimension != image->ndim)
                throw std::runtime_error("Blocks must be along the last dimension in the image");
        }
        
        Block & operator= (const NiftiImage &source)
        {
            if (source->datatype != image->datatype)
                throw std::runtime_error("New data does not have the same datatype as the target block");
            if (source->scl_slope != image->scl_slope || source->scl_inter != image->scl_inter)
                throw std::runtime_error("New data does not have the same scale parameters as the target block");
            
            size_t blockSize = 1;
            for (int i=1; i<dimension; i++)
                blockSize *= image->dim[i];
            
            if (blockSize != source->nvox)
                throw std::runtime_error("New data does not have the same size as the target block");
            
            blockSize *= image->nbyper;
            memcpy(static_cast<char*>(image->data) + blockSize*index, source->data, blockSize);
            return *this;
        }
        
        NiftiImageData data () const
        {
            if (image.isNull())
                return NiftiImageData();
            else
            {
                size_t blockSize = 1;
                for (int i=1; i<dimension; i++)
                    blockSize *= image->dim[i];
                return NiftiImageData(static_cast<char*>(image->data) + blockSize * index * image->nbyper, blockSize, image->datatype, static_cast<double>(image->scl_slope), static_cast<double>(image->scl_inter));
            }
        }
        
        template <typename TargetType>
        std::vector<TargetType> getData (const bool useSlope = true) const;
    };
    
#ifdef USING_R

    static int sexpTypeToNiftiType (const int sexpType)
    {
        if (sexpType == INTSXP || sexpType == LGLSXP)
            return DT_INT32;
        else if (sexpType == REALSXP)
            return DT_FLOAT64;
        else if (sexpType == CPLXSXP)
            return DT_COMPLEX128;
        else
            throw std::runtime_error("Array elements must be numeric");
    }
#endif
    
    static mat33 xformToRotation (const mat44 matrix);
    
    static std::string xformToString (const mat44 matrix);
    
    static int fileVersion (const std::string &path);
    

protected:
    nifti_image *image;         
    int *refCount;              
    void acquire (nifti_image * const image);
    
    void acquire (const NiftiImage &source)
    {
        refCount = source.refCount;
        acquire(source.image);
    }
    
    void release ();
    
    void copy (const nifti_image *source);
    
    void copy (const NiftiImage &source);
    
    void copy (const Block &source);


#ifdef USING_R

    void initFromNiftiS4 (const Rcpp::RObject &object, const bool copyData = true);
    
    void initFromMriImage (const Rcpp::RObject &object, const bool copyData = true);
    
    void initFromList (const Rcpp::RObject &object);
    
    void initFromArray (const Rcpp::RObject &object, const bool copyData = true);
   
#endif
    
    void initFromDims (const std::vector<int> &dim, const int datatype);

    void updatePixdim (const std::vector<float> &pixdim);
    
    void setPixunits (const std::vector<std::string> &pixunits);
    
public:
    NiftiImage ()
        : image(NULL), refCount(NULL) {}
    
    NiftiImage (const NiftiImage &source, const bool copy = true)
        : image(NULL), refCount(NULL)
    {
        if (copy)
            this->copy(source);
        else
            acquire(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
#endif
    }
    
    NiftiImage (const Block &source)
        : image(NULL), refCount(NULL)
    {
        this->copy(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
#endif
    }
    
    NiftiImage (nifti_image * const image, const bool copy = false)
        : image(NULL), refCount(NULL)
    {
        if (copy)
            this->copy(image);
        else
            acquire(image);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage with pointer %p (from pointer)\n", this->image);
#endif
    }
    
    NiftiImage (const std::vector<int> &dim, const int datatype);
    
    NiftiImage (const std::vector<int> &dim, const std::string &datatype);
    
    NiftiImage (const std::string &path, const bool readData = true);
    
    NiftiImage (const std::string &path, const std::vector<int> &volumes);
    
#ifdef USING_R

    NiftiImage (const SEXP object, const bool readData = true, const bool readOnly = false);
#endif
    
    virtual ~NiftiImage () { release(); }
    
    operator const nifti_image* () const { return image; }
    
    operator nifti_image* () { return image; }
   
    const nifti_image * operator-> () const { return image; }
    
    nifti_image * operator-> () { return image; }
    
    NiftiImage & operator= (const NiftiImage &source)
    {
        copy(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
#endif
        return *this;
    }
    
    NiftiImage & operator= (const Block &source)
    {
        copy(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
#endif
        return *this;
    }
    
    NiftiImage & setPersistence (const bool persistent) { return *this; }
    
    bool isNull () const { return (image == NULL); }
    
    bool isShared () const { return (refCount != NULL && *refCount > 1); }
    
    bool isPersistent () const { return false; }
    
    bool isDataScaled () const { return (image != NULL && image->scl_slope != 0.0 && (image->scl_slope != 1.0 || image->scl_inter != 0.0)); }
    
    int nDims () const
    {
        if (image == NULL)
            return 0;
        else
            return image->ndim;
    }
    
    std::vector<int> dim () const
    {
        if (image == NULL)
            return std::vector<int>();
        else
            return std::vector<int>(image->dim+1, image->dim+image->ndim+1);
    }
    
    std::vector<float> pixdim () const
    {
        if (image == NULL)
            return std::vector<float>();
        else
            return std::vector<float>(image->pixdim+1, image->pixdim+image->ndim+1);
    }
    
    NiftiImage & drop ()
    {
        int ndim = image->ndim;
        while (image->dim[ndim] < 2)
            ndim--;
        image->dim[0] = image->ndim = ndim;
        
        return *this;
    }
    
    const NiftiImageData data () const { return NiftiImageData(image); }
    
    NiftiImageData data () { return NiftiImageData(image); }
    
    template <typename TargetType>
    std::vector<TargetType> getData (const bool useSlope = true) const;
    
    NiftiImage & changeDatatype (const int datatype, const bool useSlope = false);
    
    NiftiImage & changeDatatype (const std::string &datatype, const bool useSlope = false);
    
    template <typename SourceType>
    NiftiImage & replaceData (const std::vector<SourceType> &data, const int datatype = DT_NONE);
    
    NiftiImage & replaceData (const NiftiImageData &data);
    
    NiftiImage & dropData ()
    {
        nifti_image_unload(image);
        return *this;
    }
    
    NiftiImage & rescale (const std::vector<float> &scales);
    
    NiftiImage & reorient (const int i, const int j, const int k);
    
    NiftiImage & reorient (const std::string &orientation);
    
#ifdef USING_R

    NiftiImage & update (const Rcpp::RObject &object);
#endif
    
    mat44 xform (const bool preferQuaternion = true) const;
    
    int nBlocks () const { return (image == NULL ? 0 : image->dim[image->ndim]); }
    
    const Block block (const int i) const { return Block(*this, nDims(), i); }
    
    Block block (const int i) { return Block(*this, nDims(), i); }
    
    const Block slice (const int i) const { return Block(*this, 3, i); }
    
    Block slice (const int i) { return Block(*this, 3, i); }
    
    const Block volume (const int i) const { return Block(*this, 4, i); }
    
    Block volume (const int i) { return Block(*this, 4, i); }
    
    int nChannels () const
    {
        if (image == NULL)
            return 0;
        else
        {
            switch (image->datatype)
            {
                case DT_NONE:   return 0;
                case DT_RGB24:  return 3;
                case DT_RGBA32: return 4;
                default:        return 1;
            }
        }
    }
    
    size_t nVoxels () const { return (image == NULL ? 0 : image->nvox); }
    
    std::pair<std::string,std::string> toFile (const std::string fileName, const int datatype = DT_NONE) const;
    
    std::pair<std::string,std::string> toFile (const std::string fileName, const std::string &datatype) const;
    
#ifdef USING_R
    
    Rcpp::RObject toArray () const;
    
    Rcpp::RObject toPointer (const std::string label) const;
    
    Rcpp::RObject toArrayOrPointer (const bool internal, const std::string label) const;
    
    Rcpp::RObject headerToList () const;
    
#endif

};

// Include implementations
#include "RNifti/NiftiImage_impl.h"

} // main namespace

#endif