mrimage.h

Go to the documentation of this file.

/******************************************************************************
 * Copyright 2014 Micah C Chambers (micahc.vt@gmail.com)
 *
 * NPL is free software: you can redistribute it and/or modify it under the
 * terms of the BSD 2-Clause License available in LICENSE or at
 * http://opensource.org/licenses/BSD-2-Clause
 *
 * @file mrimage.h
 * @brief This file contains the definition for MRImage and its derived types.
 * The derived types are templated over dimensionality and pixel type.
 ******************************************************************************/

#ifndef NDIMAGE_H
#define NDIMAGE_H

#include "npltypes.h"
#include "ndarray.h"
#include "nifti.h"

#include "zlib.h"

#include <string>
#include <iomanip>
#include <cassert>
#include <memory>
#include <map>

namespace npl {

using std::vector;

enum SliceOrderT {UNKNOWN_SLICE=0, SEQ=1, RSEQ=2, ALT=3, RALT=4, ALT_SHFT=5,
    RALT_SHFT=6};

enum CoordinateT {NOFORM=0, QFORM=1, SFORM=2};

enum BoundaryConditionT {ZEROFLUX=0, CONSTZERO=1, WRAP=2};

class MRImage;

/*****************************************************************************
 * Helper Function to Create MRImage's
 ****************************************************************************/
ptr<MRImage> createMRImage(size_t ndim, const size_t* size, PixelT type);

ptr<MRImage> createMRImage(const std::vector<size_t>& dim, PixelT type);

ptr<MRImage> createMRImage(size_t ndim, const size_t* size, PixelT type,
                void* ptr, std::function<void(void*)> deleter);

ptr<MRImage> createMRImage(const std::vector<size_t>& dim, PixelT type,
                void* ptr, std::function<void(void*)> deleter);

std::ostream& operator<<(std::ostream &out, const MRImage& img);

/******************************************************************************
 * Classes.
 ******************************************************************************/

class MRImage : public virtual NDArray
{
public:
    MRImage() : m_freqdim(-1), m_phasedim(-1), m_slicedim(-1),
    m_coordinate(NOFORM), m_slice_duration(0), m_slice_start(-1),
    m_slice_end(-1), m_slice_order(UNKNOWN_SLICE) {} ;

    /********************************************
     * Orientation Functions/Variables
     *******************************************/

    void orientDefault();

    bool isOriented() { return m_coordinate != NOFORM; };

    void setOrient(const VectorXd& neworig, const VectorXd& newspace,
                    const MatrixXd& newdir, bool reinit = true,
                    CoordinateT coord = QFORM);

    const double& direction(int64_t row, int64_t col) const;

    const double& invdirection(int64_t row, int64_t col) const;

    const MatrixXd& getDirection() const;

    void setDirection(const MatrixXd& newdir, bool reinit, CoordinateT coord = QFORM);

    double& origin(int64_t row);

    const double& origin(int64_t row) const;

    const VectorXd& getOrigin() const;

    void setOrigin(const VectorXd& neworigin, bool reinit);

    double& spacing(int64_t row);

    const double& spacing(int64_t row) const;

    const VectorXd& getSpacing() const;

    void setSpacing(const VectorXd& newspacing, bool reinit);

    /**********************************
     * Orientation Transform Functions
     *********************************/
    /*
     * @brief Converts an index in pixel space to RAS, physical/time coordinates.
     * If len < dimensions, additional dimensions are assumed to be 0. If len >
     * dimensions then additional values are ignored, and only the first DIM
     * values will be transformed and written to ras.
     *
     * @param len Length of xyz/ras arrays.
     * @param xyz Array in xyz... coordinates (maybe as long as you want).
     * @param ras Corresponding coordinate
     *
     * @return 0
     */
    virtual int indexToPoint(size_t len, const int64_t* xyz, double* ras) const=0;

    virtual int indexToPoint(size_t len, const double* xyz, double* ras) const=0;

    virtual int pointToIndex(size_t len, const double* ras, double* xyz) const=0;

    virtual int pointToIndex(size_t len, const double* ras, int64_t* index) const=0;

    virtual int orientVector(size_t len, const double* xyz, double* ras) const=0;

    virtual int disOrientVector(size_t len, const double* ras, double* xyz) const=0;

    virtual bool pointInsideFOV(size_t len, const double* ras) const=0;

    virtual bool indexInsideFOV(size_t len, const double* xyz) const=0;

    virtual bool indexInsideFOV(size_t len, const int64_t* xyz) const=0;

    virtual bool matchingOrient(ptr<const MRImage> other, bool checkdim,
            bool checksize, double tol = 0.01) const;

    virtual bool isIsotropic(bool only3d = true, double tol = 0.01) const;

    /********************************************
     * Output Functions
     *******************************************/

    virtual int write(std::string filename, double version = 1) const = 0;

    /********************************************
     * Copying/Pointer Functions
     *******************************************/

    ptr<MRImage> getPtr()  {
            return dPtrCast<MRImage>(shared_from_this());
    };

    ptr<const MRImage> getConstPtr() const {
            return dPtrCast<const MRImage>(shared_from_this());
    };


    virtual ptr<MRImage> cloneImage() const = 0;


    virtual ptr<NDArray> copy() const = 0;

    virtual ptr<NDArray> createAnother() const = 0;

    virtual ptr<NDArray> createAnother(size_t newdims, const size_t* newsize,
                    PixelT newtype) const = 0;

    virtual ptr<NDArray> createAnother(PixelT newtype) const = 0;

    virtual ptr<NDArray> createAnother(size_t newdims,
                    const size_t* newsize) const = 0;


    virtual ptr<NDArray> copyCast(size_t newdims, const size_t* newsize,
                    PixelT newtype) const = 0;

    virtual ptr<NDArray> copyCast(PixelT newtype) const = 0;

    virtual ptr<NDArray> copyCast(size_t newdims,
                    const size_t* newsize) const = 0;

    virtual ptr<NDArray> extractCast(size_t len, const int64_t* index,
                    const size_t* size) const = 0;

    virtual ptr<NDArray> extractCast(size_t len, const size_t* size) const = 0;

    virtual ptr<NDArray> extractCast(size_t len,
                    const int64_t* index, const size_t* size, PixelT newtype) const = 0;

    virtual ptr<NDArray> extractCast(size_t len, const size_t* size,
                    PixelT newtype) const = 0;


    virtual void copyMetadata(ptr<const MRImage> src);


    //  virtual int unary(double(*func)(double,double)) const = 0;
    //  virtual int binOp(const MRImage* right, double(*func)(double,double), bool elevR) const = 0;

    /**********************************************************************
     * Medical Image Specific
     *********************************************************************/

    // < 0 indicate unset variables
    int m_freqdim;
    int m_phasedim;
    int m_slicedim;
    CoordinateT m_coordinate;

    /*
     * nifti specific stuff, eventually these should be moved to a nifti
     * image subclass
     */

    void updateSliceTiming(double duration, int start, int end, SliceOrderT order);

    // raw values for slice data, < 0 indicate unset
    double m_slice_duration;
    int m_slice_start;
    int m_slice_end;

    // SEQ, RSEQ, ALT, RALT, ALT_SHFT, RALT_SHFT
    // SEQ (sequential):    slice_start .. slice_end
    // RSEQ (reverse seq):  slice_end .. slice_start
    // ALT (alternated):    slice_start, slice_start+2, .. slice_end|slice_end-1,
    //                      slice_start+1 .. slice_end|slice_end-1
    // RALT (reverse alt):  slice_end, slice_end-2, .. slice_start|slice_start+1,
    //                      slice_end-1 .. slice_start|slice_start+1
    // ALT_SHFT (siemens alt):slice_start+1, slice_start+3, .. slice_end|slice_end-1,
    //                      slice_start .. slice_end|slice_end-1
    // RALT (reverse alt):  slice_end-1, slice_end-3, .. slice_start|slice_start+1,
    //                      slice_end-2 .. slice_start|slice_start+1
    SliceOrderT m_slice_order;

    // each slice is given its relative time, with 0 as the first
    std::map<int64_t,double> m_slice_timing;

protected:
    virtual int writeNifti1Image(gzFile file) const = 0;
    virtual int writeNifti2Image(gzFile file) const = 0;

    MatrixXd m_direction;

    MatrixXd m_inv_direction;

    VectorXd m_spacing;

    VectorXd m_origin;

    friend ptr<MRImage> readNiftiImage(gzFile file, bool verbose);
};

template <size_t D, typename T>
class MRImageStore :  public virtual NDArrayStore<D,T>, public virtual MRImage
{

public:

    /*****************************************
     * Constructors
     ****************************************/
    MRImageStore(std::initializer_list<size_t> a_args);

    MRImageStore(const std::vector<size_t>& a_args);

    MRImageStore(size_t len, const size_t* size);

    MRImageStore(size_t len, const size_t* size, T* ptr,
                    const std::function<void(void*)>& deleter);

    MRImageStore();

    /*************************************************************************
     * Coordinate Transform Functions
     ************************************************************************/

    virtual int indexToPoint(size_t len, const int64_t* xyz, double* ras) const;

    virtual int indexToPoint(size_t len, const double* xyz, double* ras) const;

    virtual int pointToIndex(size_t len, const double* ras, double* xyz) const;

    virtual int pointToIndex(size_t len, const double* ras, int64_t* index) const;

    virtual int orientVector(size_t len, const double* xyz, double* ras) const;

    virtual int disOrientVector(size_t len, const double* ras, double* xyz) const;

    virtual bool pointInsideFOV(size_t len, const double* ras) const;

    virtual bool indexInsideFOV(size_t len, const double* xyz) const;

    virtual bool indexInsideFOV(size_t len, const int64_t* xyz) const;

    std::string getUnits(size_t d) { return m_units[d]; };

    /********************************************
     * Output Functions
     *******************************************/
    int write(std::string filename, double version) const;

    void printSelf();

    /********************************************
     * Copying/Pointer Functions
     *******************************************/

    virtual ptr<NDArray> copy() const;

    virtual ptr<NDArray> createAnother() const;

    virtual ptr<NDArray> createAnother(size_t newdims, const size_t* newsize,
                    PixelT newtype) const;

    virtual ptr<NDArray> createAnother(PixelT newtype) const;

    virtual ptr<NDArray> createAnother(size_t newdims,
                    const size_t* newsize) const;

    virtual ptr<NDArray> copyCast(size_t newdims, const size_t* newsize,
                    PixelT newtype) const;

    virtual ptr<NDArray> copyCast(PixelT newtype) const;

    virtual ptr<NDArray> copyCast(size_t newdims, const size_t* newsize) const;

    virtual ptr<NDArray> extractCast(size_t len, const int64_t* index,
                    const size_t* size) const;

    virtual ptr<NDArray> extractCast(size_t len, const size_t* size) const;

    virtual ptr<NDArray> extractCast(size_t len,
                    const int64_t* index, const size_t* size, PixelT newtype) const;

    virtual ptr<NDArray> extractCast(size_t len, const size_t* size,
                    PixelT newtype) const;

    ptr<MRImage> cloneImage() const;

protected:

    std::string m_units[D];

    int writeNifti1Image(gzFile file) const;
    int writeNifti2Image(gzFile file) const;
    int writeNifti1Header(gzFile file) const;
    int writeNifti2Header(gzFile file) const;
//  int writePixels(gzFile file) const;
    int writeJSON(gzFile file) const;
};
} // npl
#endif