NiftiImage.h

#ifndef _NIFTI_IMAGE_H_
#define _NIFTI_IMAGE_H_
 
 
#ifdef USING_R
 
#include <Rcpp.h>
 
#else
 
#define R_NegInf -INFINITY
 
#include <stdint.h>
#include <cstddef>
#include <cmath>
#include <string>
#include <sstream>
#include <vector>
#include <list>
#include <complex>
#include <stdexcept>
#include <algorithm>
#include <map>
#include <locale>
#include <limits>
 
#endif

 
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;
            rValue.r = value.real();
            rValue.i = 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
    {
    private:
        // NB: "parent" cannot be a reference because reference members are immutable. That renders
        // the class non-copy-assignable, which is a requirement for iterators (issue #31)
        const NiftiImageData *parent;
        void *ptr;
        size_t step;
        
    public:
        // Standard iterator typedefs
        typedef std::random_access_iterator_tag iterator_category;
        typedef Element value_type;
        typedef std::ptrdiff_t difference_type;
        typedef Element* pointer;
        typedef Element& reference;
        
        Iterator (const NiftiImageData *parent = NULL, 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
 

template <typename ElementType, int Length>
class Vector
{
protected:
    ElementType elements[Length];
    
public:
    Vector (const ElementType value = 0.0)
    {
        std::fill(elements, elements + Length, value);
    }
    
    Vector (const ElementType * source)
    {
        std::copy(source, source + Length, this->elements);
    }
    
    Vector<ElementType,Length> operator- () const
    {
        Vector<ElementType,Length> result;
        for (int i=0; i<Length; i++)
            result.elements[i] = -elements[i];
        return result;
    }
    
    const ElementType & operator[] (const size_t i) const { return elements[i]; }
    
    ElementType & operator[] (const size_t i) { return elements[i]; }
};
 

template <class NiftiType, typename ElementType, int Order>
class SquareMatrix
{
protected:
    ElementType elements[Order*Order];          
    
    NiftiType * niftiPointer () const { return (NiftiType *) elements; }
    
    NiftiType niftiCopy () const
    {
        NiftiType value;
        std::copy(elements, elements + Order*Order, *value.m);
        return value;
    }
    
public:
    typedef NiftiType NativeType;                                       
    typedef SquareMatrix<NiftiType,ElementType,Order> MatrixType;       
    typedef Vector<ElementType,Order> VectorType;                       
    
    SquareMatrix (const ElementType value = 0.0)
    {
        std::fill(elements, elements + Order*Order, value);
    }
    
    SquareMatrix (const ElementType * source)
    {
        std::copy(source, source + Order*Order, this->elements);
    }
    
    SquareMatrix (const NiftiType &source)
    {
        const ElementType *castSource = (const ElementType *) *source.m;
        std::copy(castSource, castSource + Order*Order, this->elements);
    }
    
#ifdef USING_R
    SquareMatrix (SEXP source)
    {
        Rcpp::NumericMatrix matrix(source);
        if (matrix.cols() != Order && matrix.rows() != Order)
            throw std::runtime_error("Matrix does not have the expected dimensions");
        for (int i=0; i<Order; i++)
        {
            for (int j=0; j<Order; j++)
                elements[j + i*Order] = matrix(i,j);
        }
    }
#endif
    
    operator const NiftiType () const { return niftiCopy(); }
    
    operator NiftiType () { return niftiCopy(); }
    
    const ElementType * begin () const { return elements; }
    
    ElementType * begin () { return elements; }
    
    const ElementType * end () const { return elements + Order*Order; }
    
    ElementType * end () { return elements + Order*Order; }
    
    static MatrixType eye ()
    {
        MatrixType matrix;
        for (int i=0; i<Order; i++)
            matrix.elements[i + i*Order] = 1.0;
        return matrix;
    }
    
    MatrixType inverse () const;                                
    MatrixType polar () const;                                  
    ElementType colnorm () const;                               
    ElementType rownorm () const;                               
    ElementType determ () const;                                
    MatrixType multiply (const MatrixType &other) const;        
    VectorType multiply (const VectorType &vec) const;          
    
    MatrixType operator* (const MatrixType &other) const { return multiply(other); }    
    VectorType operator* (const VectorType &vec) const { return multiply(vec); }        
    
    const ElementType & operator() (const int i, const int j) const { return elements[j + i*Order]; }
    
    ElementType & operator() (const int i, const int j) { return elements[j + i*Order]; }
    
#ifdef USING_R
    operator SEXP () const
    {
        Rcpp::NumericMatrix result(Order, Order);
        for (int i=0; i<Order; i++)
        {
            for (int j=0; j<Order; j++)
                result(i,j) = elements[j + i*Order];
        }
        return result;
    }
#endif
};
 
 
// Include matrix implementations
#include "RNifti/NiftiImage_matrix.h"
 

class NiftiImage
{
public:
#if RNIFTI_NIFTILIB_VERSION == 1
    typedef int dim_t;                  
    typedef float pixdim_t;             
    typedef float scale_t;              
#elif RNIFTI_NIFTILIB_VERSION == 2
    typedef int64_t dim_t;              
    typedef double pixdim_t;            
    typedef double scale_t;             
#endif
    
    struct Block
    {
        const NiftiImage &image;        
        const int dimension;            
        const dim_t index;              
        
        Block (const NiftiImage &image, const int dimension, const dim_t 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 != size_t(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;
    };
    
    class Extension
    {
    protected:
        nifti1_extension *ext;          
        
        void copy (const nifti1_extension *source);
        
        template <typename SourceType>
        void copy (const SourceType *data, const size_t length, const int code);
        
    public:
        Extension ()
            : ext(NULL) {}
        
        Extension (nifti1_extension * const extension, const bool copy = false)
        {
            if (!copy || extension == NULL)
                this->ext = extension;
            else
                this->copy(extension);
        }
        
        Extension (const Extension &source)
        {
            copy(source.ext);
        }
        
        template <typename SourceType>
        Extension (const SourceType *data, const size_t length, const int code)
        {
            copy(data, length, code);
        }
        
#ifdef USING_R
        Extension (SEXP source, int code = -1)
        {
            const Rcpp::RObject object(source);
            if (code == -1 && object.hasAttribute("code"))
                code = Rcpp::as<int>(object.attr("code"));
            
            switch (object.sexp_type())
            {
                case RAWSXP:  copy(RAW(source), Rf_length(source), code);       break;
                case REALSXP: copy(REAL(source), Rf_length(source), code);      break;
                case CPLXSXP: copy(COMPLEX(source), Rf_length(source), code);   break;
                case INTSXP:  copy(INTEGER(source), Rf_length(source), code);   break;
                case LGLSXP:  copy(LOGICAL(source), Rf_length(source), code);   break;
                case STRSXP:
                {
                    if (Rf_length(source) > 1)
                        Rf_warning("Character vector elements after the first will not be stored in a NIfTI extension");
                    const char *string = CHAR(STRING_ELT(source, 0));
                    copy(string, strlen(string), code);
                    break;
                }
                default: Rf_error("Unable to convert SEXP type %d to NIfTI extension", object.sexp_type());
            }
        }
#endif
        
        int code () const { return (ext == NULL ? -1 : ext->ecode); }
        
        const char * data () const { return (ext == NULL ? NULL : ext->edata); }
        
        size_t length () const { return (ext == NULL || ext->esize < 8 ? 0 : size_t(ext->esize - 8)); }
        
        size_t size () const { return (ext == NULL || ext->esize < 8 ? 0 : size_t(ext->esize - 8)); }
        
#ifdef USING_R
        operator SEXP () const
        {
            if (ext == NULL || ext->esize < 8)
                return R_NilValue;
            
            const int length = ext->esize - 8;
            Rcpp::RawVector result(length);
            const Rbyte *source = (const Rbyte *) ext->edata;
            std::copy(source, source+length, result.begin());
            result.attr("code") = ext->ecode;
            return result;
        }
#endif
    };
    
    class Xform
    {
    public:
#if RNIFTI_NIFTILIB_VERSION == 1
        typedef float Element;                                      
        typedef Vector<float,4> Vector4;                            
        typedef Vector<float,3> Vector3;                            
        typedef SquareMatrix<mat44,float,4> Matrix;                 
        typedef SquareMatrix<mat33,float,3> Submatrix;              
#elif RNIFTI_NIFTILIB_VERSION == 2
        typedef double Element;                                     
        typedef Vector<double,4> Vector4;                           
        typedef Vector<double,3> Vector3;                           
        typedef SquareMatrix<nifti_dmat44,double,4> Matrix;         
        typedef SquareMatrix<nifti_dmat33,double,3> Submatrix;      
#endif
        
    protected:
        Element *forward, *inverse, *qparams;                       
        Matrix mat;                                                 
        
        void replace (const Matrix &source);
        
    public:
        Xform ()
            : forward(NULL), inverse(NULL), qparams(NULL), mat() {}
        
        Xform (const Matrix &source)
            : forward(NULL), inverse(NULL), qparams(NULL), mat(source) {}
        
        Xform (const Matrix::NativeType &source)
            : forward(NULL), inverse(NULL), qparams(NULL), mat(source) {}
        
        Xform (Matrix::NativeType &source)
            : forward(*source.m), inverse(NULL), qparams(NULL), mat(source) {}
        
        Xform (Matrix::NativeType &source, Matrix::NativeType &inverse, Element *qparams = NULL)
            : forward(*source.m), inverse(*inverse.m), qparams(qparams), mat(source) {}
        
#ifdef USING_R
        Xform (SEXP source)
            : forward(NULL), inverse(NULL), qparams(NULL), mat(Matrix(source)) {}
#endif
        
        operator const Matrix::NativeType () const { return mat; }
        
        operator Matrix::NativeType () { return mat; }
        
        Xform & operator= (const Xform &source)
        {
            replace(source.mat);
            return *this;
        }
        
        Xform & operator= (const Matrix &source)
        {
            replace(source);
            return *this;
        }
        
#ifdef USING_R
        Xform & operator= (SEXP source)
        {
            replace(Matrix(source));
            return *this;
        }
#endif
        
        const Matrix & matrix () const { return mat; }
        
        Submatrix submatrix () const;
        
        Submatrix rotation () const;
        
        Element handedness () const;
        
        Vector4 quaternion () const;
        
        Vector3 offset () const;
        
        Vector3 spacing () const;
        
        std::string orientation () 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 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<dim_t> &dim, const int datatype);

    void updatePixdim (const std::vector<pixdim_t> &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 (v%d) with pointer %p (from NiftiImage)\n", RNIFTI_NIFTILIB_VERSION, (void *) this->image);
#endif
    }
    
    NiftiImage (const Block &source)
        : image(NULL), refCount(NULL)
    {
        this->copy(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage (v%d) with pointer %p (from Block)\n", RNIFTI_NIFTILIB_VERSION, (void *) 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 (v%d) with pointer %p (from pointer)\n", RNIFTI_NIFTILIB_VERSION, (void *) this->image);
#endif
    }
    
    NiftiImage (const nifti_1_header &header)
        : image(NULL), refCount(NULL)
    {
#if RNIFTI_NIFTILIB_VERSION == 1
        acquire(nifti_convert_nhdr2nim(header, NULL));
#elif RNIFTI_NIFTILIB_VERSION == 2
        acquire(nifti_convert_n1hdr2nim(header, NULL));
#endif
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage (v%d) with pointer %p (from header)\n", RNIFTI_NIFTILIB_VERSION, (void *) this->image);
#endif
    }
    
#if RNIFTI_NIFTILIB_VERSION == 2
    NiftiImage (const nifti_2_header &header)
        : image(NULL), refCount(NULL)
    {
        acquire(nifti_convert_n2hdr2nim(header, NULL));
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage (v%d) with pointer %p (from header)\n", RNIFTI_NIFTILIB_VERSION, (void *) this->image);
#endif
    }
#endif
    
    NiftiImage (const std::vector<dim_t> &dim, const int datatype);
    
    NiftiImage (const std::vector<dim_t> &dim, const std::string &datatype);
    
    NiftiImage (const std::string &path, const bool readData = true);
    
    NiftiImage (const std::string &path, const std::vector<dim_t> &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 (v%d), with pointer %p (from NiftiImage)\n", RNIFTI_NIFTILIB_VERSION, (void *) this->image);
#endif
        return *this;
    }
    
    NiftiImage & operator= (const Block &source)
    {
        copy(source);
#ifndef NDEBUG
        Rc_printf("Creating NiftiImage (v%d) with pointer %p (from Block)\n", RNIFTI_NIFTILIB_VERSION, (void *) 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<dim_t> dim () const
    {
        if (image == NULL)
            return std::vector<dim_t>();
        else
            return std::vector<dim_t>(image->dim+1, image->dim+image->ndim+1);
    }
    
    std::vector<pixdim_t> pixdim () const
    {
        if (image == NULL)
            return std::vector<pixdim_t>();
        else
            return std::vector<pixdim_t>(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 ()
    {
#if RNIFTI_NIFTILIB_VERSION == 1
        nifti_image_unload(image);
#elif RNIFTI_NIFTILIB_VERSION == 2
        nifti2_image_unload(image);
#endif
        return *this;
    }
    
    NiftiImage & rescale (const std::vector<pixdim_t> &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
    
    const Xform xform (const bool preferQuaternion = true) const;
    
    const Xform qform () const { return (image == NULL ? Xform() : Xform(image->qto_xyz)); }
    
    Xform qform () { return (image == NULL ? Xform() : Xform(image->qto_xyz, image->qto_ijk, &image->quatern_b)); }
    
    const Xform sform () const { return (image == NULL ? Xform() : Xform(image->sto_xyz)); }
    
    Xform sform () { return (image == NULL ? Xform() : Xform(image->sto_xyz, image->sto_ijk)); }
    
    dim_t 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); }
    
    int nExtensions () const { return (image == NULL ? 0 : image->num_ext); }
    
    std::list<Extension> extensions (const int code = -1) const
    {
        if (image == NULL)
            return std::list<Extension>();
        else
        {
            std::list<Extension> result;
            for (int i=0; i<image->num_ext; i++)
            {
                const Extension extension(image->ext_list + i);
                if (code < 0 || code == extension.code())
                    result.push_back(extension);
            }
            return result;
        }
    }
    
    NiftiImage & addExtension (const Extension &extension)
    {
        if (image != NULL)
#if RNIFTI_NIFTILIB_VERSION == 1
            nifti_add_extension(image, extension.data(), int(extension.length()), extension.code());
#elif RNIFTI_NIFTILIB_VERSION == 2
            nifti2_add_extension(image, extension.data(), int(extension.length()), extension.code());
#endif
        return *this;
    }
    
    NiftiImage & replaceExtensions (const std::list<Extension> extensions)
    {
        dropExtensions();
        for (std::list<Extension>::const_iterator it=extensions.begin(); it!=extensions.end(); ++it)
            addExtension(*it);
        return *this;
    }
    
    NiftiImage & dropExtensions ()
    {
        if (image != NULL)
#if RNIFTI_NIFTILIB_VERSION == 1
            nifti_free_extensions(image);
#elif RNIFTI_NIFTILIB_VERSION == 2
            nifti2_free_extensions(image);
#endif
        return *this;
    }
    
    std::pair<std::string,std::string> toFile (const std::string fileName, const int datatype = DT_NONE, const int filetype = -1, const int compression = 6) const;
    
    std::pair<std::string,std::string> toFile (const std::string fileName, const std::string &datatype, const int filetype = -1, const int compression = 6) 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;
    
#endif
 
};
 
// Include image implementations
#include "RNifti/NiftiImage_impl.h"
 
} // main namespace
 
#endif