Grid3_i.h

#ifndef GRID3_I_H
#define GRID3_I_H

#include "Endian.h"
#include "FileName.h"
#include <typeinfo>
#include <sstream>
#include <fstream>
using namespace std;

// Configure gdcm support
// Define 'grid3_use_gdcm' to enable gdcm support.

#ifdef grid3_use_gdcm
#include "gdcmImageReader.h"
#endif

inline string ToLower(string s){
    string result=s;
    for(int i=0; i<s.size(); i++){
        result[i]=tolower(s[i]);
    }
    return result;

}

template <class T>
inline Grid3<T>::Grid3(int width, int height, int depth, double x_spacing, double y_spacing, double z_spacing){
    _width=width;
    _height=height;
    _depth=depth;
    _spacingX=x_spacing;
    _spacingY=y_spacing;
    _spacingZ=z_spacing;
    data=std::vector<T>(NumberOfElements());        
}

template <class T>
inline Grid3<T>::Grid3(Vec3i dimensions, Vec3d spacing){
    _width=dimensions.x;
    _height=dimensions.y;
    _depth=dimensions.z;
    _spacingX=spacing.x;
    _spacingY=spacing.y;
    _spacingZ=spacing.z;
    data=std::vector<T>(NumberOfElements());   
}

template <class T>
inline void Grid3<T>::Init(int width, int height, int depth, double x_spacing, double y_spacing, double z_spacing) {
    _width=width;
    _height=height;
    _depth=depth;
    _spacingX=x_spacing;
    _spacingY=y_spacing;
    _spacingZ=z_spacing;
    data=std::vector<T>(NumberOfElements());        
}


template <class T>
inline void Grid3<T>::SetSpacing(double x_spacing, double y_spacing, double z_spacing){
    _spacingX=x_spacing;
    _spacingY=y_spacing;    
    _spacingZ=z_spacing;
}
    
template <class T>
inline void Grid3<T>::SetSpacing(Vec3d spacing){
    _spacingX=spacing.x;
    _spacingY=spacing.y;
    _spacingZ=spacing.z;
}

template <class T>
inline Vec3d Grid3<T>::WorldToGrid(Vec3d p){
    Vec3d result;
    result.x = p.x / GetSpacingX();
    result.y = p.y / GetSpacingY();
    result.z = p.z / GetSpacingZ();
    return result;
}

template <class T>
inline Vec3d Grid3<T>::GridToWorld(Vec3d p){
    Vec3d result;
    result.x = p.x * GetSpacingX();
    result.y = p.y * GetSpacingY();
    result.z = p.z * GetSpacingZ();
    return result;
}

template <class T>
inline int Grid3<T>::NumberOfElements(){ 
    return _width*_height*_depth;
}

template <class T>
inline T* Grid3<T>::Data(){ 
    return &data[0]; 
}

template <>
inline std::string Grid3<uint8_t>::GetDataType(){
    return "UINT8";
}

template <>
inline std::string Grid3<int8_t>::GetDataType(){
    return "INT8";
}

template <>
inline std::string Grid3<uint16_t>::GetDataType(){
    return "UINT16";
}

template <>
inline std::string Grid3<int16_t>::GetDataType(){
    return "INT16";
}

template <>
inline std::string Grid3<uint32_t>::GetDataType(){
    return "UINT32";
}

template <>
inline std::string Grid3<int32_t>::GetDataType(){
    return "INT32";
}

template <class T>
inline std::string Grid3<T>::GetDataType(){
    return std::string( typeid(T).name() );
}


template <class T>
inline int Grid3<T>::Offset(int x, int y, int z){ 
    return x+ y*_width + z*_width*_height; 
}

template <class T>
inline Vec3i Grid3<T>::Coordinate(long index){
    Vec3i c;        
    c.z = index/(_width*_height);   
    index = index%(_width*_height);
    c.y = index/_width;
    c.x = index%_width;
    return c;
}

template <class T>
inline bool Grid3<T>::IsInside(int x, int y, int z){
    return (x>=0 && x<_width && y>=0 && y<_height && z>=0 && z<_depth);
}

template <class T>
inline bool Grid3<T>::IsInside(Vec3i p){ 
    return IsInside(p.x,p.y,p.z); 
}

template <class T>
inline T& Grid3<T>::operator[](long index){
    return data[index];
}

template <class T>
inline T& Grid3<T>::operator()(int x,int y, int z, bool Safe){
    if(Safe){
        x=max(x,0);
        x=min(x,_width-1);
        y=max(y,0);
        y=min(y,_height-1);
        z=max(z,0);
        z=min(z,_depth-1);
    }
    return data[Offset(x,y,z)];
}

template <class T>
inline T& Grid3<T>::operator()(Vec3i p, bool Safe){
    if(Safe){
        p.x=max(p.x,0);
        p.x=min(p.x,_width-1);
        p.y=max(p.y,0);
        p.y=min(p.y,_height-1);
        p.z=max(p.z,0);
        p.z=min(p.z,_depth-1);
    }
    return data[Offset(p.x,p.y,p.z)];
}

template <class T>
inline T Grid3<T>::LinearAt(double x, double y, double z){
    double xt=x-floor(x);
    double yt=y-floor(y);
    double zt=z-floor(z);

    int x1=int(floor(x));
    int x2=int(ceil(x));
    int y1=int(floor(y));
    int y2=int(ceil(y));
    int z1=int(floor(z));
    int z2=int(ceil(z));

    return (1-zt)*((1-yt)*((1-xt)*operator()(x1,y1,z1,true)+
        (xt)*operator()(x2,y1,z1,true))+ 
        (yt)*((1-xt)*operator()(x1,y2,z1,true)+
        (xt)*operator()(x2,y2,z1,true)))+
        (zt)*((1-yt)*((1-xt)*operator()(x1,y1,z2,true)+
        (xt)*operator()(x2,y1,z2,true))+ 
        (yt)*((1-xt)*operator()(x1,y2,z2,true)+
        (xt)*operator()(x2,y2,z2,true)));
}

template <class T>
inline T Grid3<T>::LinearAt(Vec3d p){
    return LinearAt(p.x,p.y,p.z);
}

template <class T>
inline void Grid3<T>::Get6Neighbors(long index, vector<long>& n){
    n.clear();

    Vec3i coord=Coordinate(index);
    if(coord.x+1<_width){
        n.push_back(Offset(coord.x+1, coord.y, coord.z));
    }

    if(coord.x-1>=0){
        n.push_back(Offset(coord.x-1, coord.y, coord.z));
    }

    if(coord.y+1<_height){
        n.push_back(Offset(coord.x, coord.y+1, coord.z));
    }

    if(coord.y- 1>=0){
        n.push_back(Offset(coord.x, coord.y-1, coord.z));
    }

    if(coord.z+1<_depth){
        n.push_back(Offset(coord.x, coord.y, coord.z+1));
    }

    if(coord.z- 1>=0){
        n.push_back(Offset(coord.x, coord.y, coord.z-1));
    }

}

template <class T>
inline void Grid3<T>::Get6Neighbors(Vec3i p, vector<Vec3i>& n){
    n.clear();

    if(p.x+1<_width){
        n.push_back(Vec3i(p.x+1, p.y, p.z));
    }

    if(p.x-1>=0){
        n.push_back(Vec3i(p.x-1, p.y, p.z));
    }

    if(p.y+1<_height){
        n.push_back(Vec3i(p.x, p.y+1, p.z));
    }

    if(p.y-1>=0){
        n.push_back(Vec3i(p.x, p.y-1, p.z));
    }

    if(p.z+1<_depth){
        n.push_back(Vec3i(p.x, p.y, p.z+1));
    }

    if(p.z-1>=0){
        n.push_back(Vec3i(p.x, p.y, p.z-1));
    }
}

template <class T>
inline void Grid3<T>::Fill(T value){
    std::fill(data.begin(),data.end(),value);
}

template <class T>
inline T Grid3<T>::GetMaxValue(){
    return *std::max_element(data.begin(), data.end());
}

template <class T>
inline T Grid3<T>::GetMinValue(){
    return *std::min_element(data.begin(), data.end());
}


/*template <class T>
inline Grid3<T> Grid3<T>::operator*(double x){
    Grid3<T> result(this->GetDimensions(), this->GetSpacing());

    #pragma omp parallel for
    for (long i = 0; i < data.size(); i++){
        result.data[i] = x*data[i];
    }
    return result;
}

template <class T>
inline Grid3<T> Grid3<T>::operator/ (double x){
    Grid3<T> result(this->GetDimensions(), this->GetSpacing());
    double x_inv = 1 / x;
    #pragma omp parallel for
    for (long i = 0; i < data.size(); i++){
        result.data[i] = x_inv*data[i];
    }
    return result;
}

template <class T>
inline Grid3<T> Grid3<T>::operator+(T x){
    Grid3<T> result(this->GetDimensions(), this->GetSpacing());
#pragma omp parallel for
    for (long i = 0; i < data.size(); i++){
        result.data[i] += x;
    }
    return result;
}

template <class T>
inline Grid3<T> Grid3<T>::operator-(T x){
    Grid3<T> result(this->GetDimensions(), this->GetSpacing());
#pragma omp parallel for
    for (long i = 0; i < data.size(); i++){
        result.data[i] -= x;
    }
    return result;
}*/

template <class T>
inline void Grid3<T>::CopyDataFrom(Grid3<T>& src){
    if(src.NumberOfElements() == this->NumberOfElements()){
        std::copy(src.data.begin(), src.data.end(), data.begin());
    }
}

template <class T>
inline void Grid3<T>::SwapContentsWith(Grid3<T>& v){
    if(this->NumberOfElements() == v.NumberOfElements()){
        std::swap_ranges(v.data.begin(), v.data.end(), data.begin());
    }
}

template <class T>
inline bool Grid3<T>::CreateFromFile(std::string filename){
    if(filename.find(".vtk")!=std::string::npos){
        return ReadVTK(filename);
    }
    else if(filename.find(".mhd")!=string::npos){
        return ReadMetaImage(filename);
    }
    else if(filename.find(".mha")!=string::npos){
        return ReadMetaImage(filename);
    }
    else if(filename.find(".nii")!=string::npos){
        return ReadNIFTI(filename);
    }
    else if(filename.find(".dcm")!=string::npos){
        return ReadDICOM(filename);
    }
    return false;
}

template <class T>
inline bool Grid3<T>::ReadVTK(std::string filename){

    //   # vtk DataFile Version x.x
    //   Some information about the file
    //   BINARY
    //   DATASET STRUCTURED_POINTS
    //   DIMENSIONS 128 128 128
    //   ORIGIN 0.0 0.0 0.0
    //   SPACING 1.0 1.0 1.0
    //   POINT_DATA 2097152
    //   SCALARS image_data unsigned_char
    //   LOOKUP_TABLE default
    //   raw data........

    //Note: Binary VTK legacy files are assumed to ALWAYS be encoded in big endian format.

    // open file, use binary mode (necessary on windows)
    ifstream is;
    is.open( filename.c_str(), ios::in | ios::binary );
    
    if( !is ) { 
        //cout<<"Stream invalid"<<endl;
        return false; 
    }

    std::stringstream ss;
    std::string line;

    // check if it is a VTK file
    getline(is,line);
    ss.str(line);

    std::string checkvtk;
    ss >> checkvtk; // '#'
    ss >> checkvtk; // vtk
    if( !(ToLower(checkvtk) == "vtk" )){ return false; }

    bool done=false;
    string keyword;

    bool read_dimensions=false;
    bool read_lut=false;
    bool read_type=false;
    bool read_asciibinary=false;

    bool binary=true;
    string dataset_type;
    string typestring;
    long W,H,D;
    double spacing_x=1,spacing_y=1,spacing_z=1;
    double origin_x=0,origin_y=0,origin_z=0;

    while(!done && getline(is,line)){
        //ss=std::stringstream(line);
        std::stringstream ss2(line);

        ss2>>keyword;
        if(ToLower(keyword)=="dimensions"){

            ss2>>W;
            ss2>>H;
            ss2>>D;
            read_dimensions=true;
        }
        else if(ToLower(keyword)=="binary"){
            read_asciibinary=true;
        }
        else if(ToLower(keyword)=="ascii"){
            binary=false;
            read_asciibinary=true;
        }
        else if(ToLower(keyword)=="dataset"){
            ss2>>dataset_type;
            if(!(ToLower(dataset_type)=="structured_points")){ 
                return false;
            }
        }
        else if(ToLower(keyword)=="spacing" ||ToLower(keyword)=="aspect_ratio" ){
            ss2>>spacing_x;
            ss2>>spacing_y;
            ss2>>spacing_z;
        }
        else if(ToLower(keyword)=="origin"){
            ss2>>origin_x;
            ss2>>origin_y;
            ss2>>origin_z;
        }
        else if(ToLower(keyword)=="scalars"){
            ss2>>typestring;
            ss2>>typestring;
            read_type=true;
        }
        else if(ToLower(keyword)=="lookup_table"){
            read_lut=true;
        }

        done=read_dimensions && read_lut && read_type && read_asciibinary;
    }

    if(!done){
        if(read_dimensions){
            cout<<"dim"<<endl;
        }
        cout<<read_lut<<endl;
        cout<<read_type<<endl;
        cout<<read_asciibinary<<endl;
        return false;
    }

    long nelem=W*H*D;

    Init(W,H,D);
    T* data=Data();

    if(typestring=="unsigned_char"){ 
        uint8_t* tmpdata=new uint8_t[nelem];
        is.read((char*)tmpdata, nelem*sizeof(uint8_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT8");
    }

    else if(typestring=="unsigned_short"){ 
        uint16_t* tmpdata=new uint16_t[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint16_t), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT16");
    }

    else if(typestring=="unsigned_int"){ 
        uint32_t* tmpdata=new uint32_t[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint32_t), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT32");
    }

    else if(typestring=="char"){ 
        int8_t* tmpdata=new int8_t[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int8_t), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT8");
    }

    else if(typestring=="short"){ 
        int16_t* tmpdata=new int16_t[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int16_t), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT16");
    }

    else if(typestring=="int"){ 
        int32_t* tmpdata=new int32_t[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int32_t), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT32");
    }

    else if(typestring=="float"){ 
        float* tmpdata=new float[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(float), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "FLOAT");
    }

    else if(typestring=="double"){ 
        double* tmpdata=new double[nelem];
        ReadFromBigEndian((char*)tmpdata, nelem, sizeof(double), is);
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "DOUBLE");
    }

    _spacingX=spacing_x;
    _spacingY=spacing_y;
    _spacingZ=spacing_z;

    return true;
}

template <class T>
inline bool Grid3<T>::WriteVTK(std::string filename){
    //   # vtk DataFile Version x.x
    //   Some information about the file
    //   BINARY
    //   DATASET STRUCTURED_POINTS
    //   DIMENSIONS 128 128 128
    //   ORIGIN 0.0 0.0 0.0
    //   SPACING 1.0 1.0 1.0
    //   POINT_DATA 2097152
    //   SCALARS image_data unsigned_char
    //   LOOKUP_TABLE default
    //   raw data........

    ofstream outfile;
    outfile.open(filename,ofstream::out|ofstream::binary);
    if(!outfile) { return false; }

    long nelem= _width*_height*_depth;

    outfile<<"# vtk DataFile Version 3.0"<<endl;
    outfile<<"Created using the Grid3 library"<<endl;

    string Datamode="BINARY";
    //if(mode=="ASCII"){Datamode="ASCII";}
    outfile<<Datamode<<endl;

    outfile<<"DATASET STRUCTURED_POINTS"<<endl;

    outfile<<"DIMENSIONS "<<_width<<" "<<_height<<" "<<_depth<<endl;
    outfile<<"ORIGIN 0.0 0.0 0.0"<<endl;
    outfile<<"SPACING "<<_spacingX<<" "<<_spacingY<<" "<<_spacingZ<<endl;
    outfile<<"POINT_DATA "<<nelem<<endl;

    string datatype=GetDataType();
    string vtk_datatype="";

    if(datatype=="UINT8"){vtk_datatype="unsigned_char";}
    else if(datatype=="INT8"){vtk_datatype="char";}
    else if(datatype=="UINT16"){vtk_datatype="unsigned_short";}
    else if(datatype=="INT16"){vtk_datatype="short";}
    else if(datatype=="UINT32"){vtk_datatype="unsigned_int";}
    else if(datatype=="INT32"){vtk_datatype="int";}
    else if(datatype=="float"){vtk_datatype="float";}
    else if(datatype=="double"){vtk_datatype="double";}

    outfile<<"SCALARS image_data "<<vtk_datatype<<endl;
    outfile<<"LOOKUP_TABLE default"<<endl;

    //if(Datamode=="BINARY"){
    //ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint16_t), is);

    T* data=Data();
    WriteToBigEndian(reinterpret_cast<char*>(data), nelem, sizeof(T), outfile);

    outfile.close();

    return true;
}

/*template <class T>
inline bool Grid3<T>::ReadAmira(std::string filename){
    return false;
}*/

/*template <class T>
inline bool Grid3<T>::WriteAmira(std::string filename){
    return false;
}*/

template <class T>
inline bool Grid3<T>::ReadMetaImage(std::string filename){

    //See http://www.itk.org/Wiki/ITK/MetaIO/Documentation for fileformat specification

    map<string,string> header_info;

    ifstream f;
    f.open(filename);
    if(!f){return false;}

    string line,key, value;

    //Read key/value pairs from header file and store in header_info
    while(!f.eof()){
        getline(f,line);
        size_t pos=line.find("=");
        if(pos!=string::npos){
            line.erase(pos,1);
            key=line.substr(0,pos);
            value=line.substr(pos);

            //remove trailing and leading whitespaces
            size_t first_non_space=key.find_first_not_of(" ");
            size_t last_non_space=key.find_last_not_of(" ");
            key=key.substr(first_non_space,last_non_space+1);

            //remove trailing and leading whitespaces
            first_non_space=value.find_first_not_of(" ");
            last_non_space=value.find_last_not_of(" ");
            value=value.substr(first_non_space,last_non_space+1);

            header_info[key]=value;
        }
    }

    f.close();

    map<string,string>::iterator it;
    stringstream ss;

    //Verify that the object is an image
    it=header_info.find("ObjectType");
    if(it==header_info.end() || it->second!="Image"){ return false;}

    //Read the NDims field
    it=header_info.find("NDims");
    if(it==header_info.end()){return false;}
    int ndims;
    ss=stringstream(it->second);
    ss>>ndims;
    if(!(ndims==2 || ndims==3)){return false;}


    //Read the dimsize field
    it=header_info.find("DimSize");
    if(it==header_info.end()){return false;}
    ss=stringstream(it->second);
    int dimx, dimy, dimz;
    dimx=1;
    dimy=1;
    dimz=1;
    ss>>dimx>>dimy;;
    if(ndims==3){ss>>dimz;}


    //Read the ElementDataFile field
    it=header_info.find("ElementDataFile");
    if(it==header_info.end()){return false;}
    string datafile=it->second;

    //Read the ElementType field
    it=header_info.find("ElementType");
    if(it==header_info.end()){return false;}
    string typestring=it->second;

    //Read the binarydata field
    bool BinaryData=true;
    it=header_info.find("BinaryData");
    if(it!=header_info.end()){
        BinaryData=(it->second=="True");
    }

    //Read the BinaryDataByteOrderMSB field
    bool BinaryDataByteOrderMSB=IsBigEndian(); //default byte ordering is the same as the native ordering
    it=header_info.find("BinaryDataByteOrderMSB");
    if(it!=header_info.end()){
        BinaryDataByteOrderMSB=(it->second=="True");
    }

    //Read the HeaderSize field
    int HeaderSize=-1; 
    it=header_info.find("HeaderSize");
    if(it!=header_info.end()){
        ss=stringstream(it->second);
        ss>>HeaderSize;
    }

    //Read the ElementSíze field
    it=header_info.find("ElementSize");
    double voxelsize_x=1, voxelsize_y=1, voxelsize_z=1;
    if(it!=header_info.end()){
        ss=stringstream(it->second);
        ss>>voxelsize_x>>voxelsize_y;
        if(ndims==3){ ss>>voxelsize_z;}
    }   

    //Read the ElementSpacing field
    it=header_info.find("ElementSpacing");
    _spacingX=voxelsize_x;
    _spacingY=voxelsize_y;
    _spacingZ=voxelsize_z;
    if(it!=header_info.end()){
        ss=stringstream(it->second);
        ss>>_spacingX>>_spacingY;
        if(ndims==3){ ss>>_spacingZ;}
        
    }   

    long nelem=dimx*dimy*dimz;
    Init(dimx,dimy,dimz,_spacingX,_spacingY,_spacingZ);

    T* data=Data();


    /*typedef enum
    {
    MET_NONE,
    MET_ASCII_CHAR, 
    MET_CHAR, 
    MET_UCHAR, 
    MET_SHORT,
    MET_USHORT,
    MET_INT, 
    MET_UINT,
    MET_LONG,
    MET_ULONG,
    MET_LONG_LONG,
    MET_ULONG_LONG,
    MET_FLOAT, 
    MET_DOUBLE,
    MET_STRING, 
    MET_CHAR_ARRAY, 
    MET_UCHAR_ARRAY, 
    MET_SHORT_ARRAY,
    MET_USHORT_ARRAY,
    MET_INT_ARRAY, 
    MET_UINT_ARRAY,
    MET_LONG_ARRAY,
    MET_ULONG_ARRAY,
    MET_LONG_LONG_ARRAY,
    MET_ULONG_LONG_ARRAY,
    MET_FLOAT_ARRAY, 
    MET_DOUBLE_ARRAY,
    MET_FLOAT_MATRIX,
    MET_OTHER
    } MET_ValueEnumType;*/


    string Datafile_fullpath="";

    size_t found=filename.find_last_of("/\\");
    if(found!=string::npos){
        Datafile_fullpath.append(filename.substr(0,found+1));
    }

    Datafile_fullpath.append(datafile);


    //Determine file size automatically
    f.open(Datafile_fullpath,ios::binary | ios::ate);
    if(!f.good()){ 

        return false;
    }
    int fileSize=f.tellg();
    f.close();

    f.open(Datafile_fullpath,ios::in|ios::binary);

    if(!f.good()){ 
        cout << Datafile_fullpath << endl;
        return false;
    }

    if(typestring=="MET_UCHAR"){ 
        uint8_t* tmpdata=new uint8_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint8_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(uint8_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT8");


    }

    else if(typestring=="MET_CHAR"){ 
        int8_t* tmpdata=new int8_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int8_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(int8_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT8");

    }

    else if(typestring=="MET_USHORT"){ 
        uint16_t* tmpdata=new uint16_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint16_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(uint16_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT16");

    }

    else if(typestring=="MET_SHORT"){ 
        int16_t* tmpdata=new int16_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int16_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(int16_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT16");

    }
    else if(typestring=="MET_UINT"){ 
        uint32_t* tmpdata=new uint32_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(uint32_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(uint32_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "UINT32");

    }

    else if(typestring=="MET_INT"){ 
        int32_t* tmpdata=new int32_t[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(int32_t), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(int32_t), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "INT32");

    }

    else if(typestring=="MET_FLOAT"){ 
        float* tmpdata=new float[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(float), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(float), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "FLOAT");

    }
    else if(typestring=="MET_DOUBLE"){ 
        double* tmpdata=new double[nelem];
        if(BinaryDataByteOrderMSB){
            ReadFromBigEndian((char*)tmpdata, nelem, sizeof(double), f);
        }
        else{
            ReadFromLittleEndian((char*)tmpdata, nelem, sizeof(double), f);
        }

        #pragma omp parallel for
        for(long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
        delete[] tmpdata;
        metadata.SetString("ORIGINAL_TYPE", "DOUBLE");

    }

    f.close();


    //Set metadata

    //Read the Comment field
    it=header_info.find("Comment");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the Name field
    it=header_info.find("Name");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the ID field
    it=header_info.find("ID");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the Position field
    it=header_info.find("Position");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the Orientation field
    it=header_info.find("Orientation");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the AnatomicalOrientation field
    it=header_info.find("AnatomicalOrientation");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }

    //Read the Modality field
    it=header_info.find("Modality");
    if(it!=header_info.end()){
        metadata.SetString(it->first,it->second);
    }


    return true;
}

template <class T>
inline bool Grid3<T>::WriteMetaImage(std::string filename){

    //string datatype=GetDataType();


    FileName f(filename);
    f.SetExtension("mhd");
    string headerfile_name=f.GetAsString();

    f.SetExtension("raw");
    string datafile=f.GetAsString();

    f.ClearPath();
    string relative_datafile=f.GetAsString();
    
    ofstream outfile;
    string datatype=GetDataType();
    string metaimage_datatype="";
    if(datatype=="UINT8"){metaimage_datatype="MET_UCHAR";}
    else if(datatype=="UINT16"){metaimage_datatype="MET_USHORT";}
    else if(datatype=="UINT32"){metaimage_datatype="MET_UINT";}
    else if(datatype=="INT8"){metaimage_datatype="MET_CHAR";}
    else if(datatype=="INT16"){metaimage_datatype="MET_SHORT";}
    else if(datatype=="INT32"){metaimage_datatype="MET_INT";}
    else{return false;} //Unsupported datatype for this format

    outfile.open(headerfile_name,ofstream::out|ofstream::binary);
    if(!outfile) { return false; }
    
    outfile<<"ObjectType = Image"<<endl;
    
    int ndims=3;
    if( GetDepth()==1){ ndims=2; }
    outfile<<"NDims = "<<ndims<<endl;

    outfile<<"DimSize = "<<GetWidth()<<" "<<GetHeight(); 
    if(ndims==3){outfile<<" "<<GetDepth();} 
    outfile<<endl;

    outfile<<"HeaderSize = 0"<<endl;
    outfile<<"ElementSize = 1 1 1"<<endl;
    outfile<<"ElementSpacing = "<<GetSpacingX()<<" "<<GetSpacingY()<<" "<<GetSpacingZ()<<" "<<endl;

    outfile<<"ElementType = "<<metaimage_datatype<<endl;

    if(IsBigEndian()){
        outfile<<"ElementByteOrderMSB = True"<<endl;
    }
    else{
        outfile<<"ElementByteOrderMSB = False"<<endl;
    }

    outfile<<"ElementDataFile = "<<relative_datafile<<endl;

    outfile.close();
    
    long nelem= _width*_height*_depth;

    T* data=Data();
    outfile.open(datafile,ofstream::out|ofstream::binary);
    outfile.write((char*)data, nelem*sizeof(T));
    //WriteToBigEndian((char*)&data[0], nelem, sizeof(T), outfile);
    outfile.close();

    return true;
}

template <class T>
inline bool Grid3<T>::ReadNIFTI(std::string filename){
    //TODO: Handle Endianness
    //It is assumed that sizeof(float)==4

    //Volume* result = NULL;

    ifstream f;

    f.open( filename.c_str(), std::ios::in | std::ios::binary );
    if( !f ) { return false; }

    //Read size of the header. Must be 348 (bytes).
    int32_t sizeof_hdr=0;
    f.read((char*)&sizeof_hdr,4);
    if(sizeof_hdr!=348){ 
        
        f.close(); 
        return false; 
    }


    //Read data array dimensions. Grid3 only supports 3 dimensions.
    f.seekg(40);
    int16_t dim[8];
    f.read((char*)&dim[0],16);

    if(!(dim[0]>0 && dim[0]<5)){
        f.close(); 
        return false;
    }

    int dimx=dim[1];
    int dimy=1;
    if(dim[0]>=2){ dimy=dim[2]; } 
    int dimz=1;
    if(dim[0]>=3){ dimz=dim[3]; } 
    int nchannels=1;
    if(dim[0]>=4){ nchannels=dim[4]; } 

    //Read data type
    f.seekg(70);
    int16_t datatype=0;
    f.read((char*)&datatype,2);
    if(datatype==0){
        f.close();
        return false;
    }

    //Read grid spacings (unit per dimension)
    f.seekg(76);
    float pixdim[8]; 
    f.read((char*)&pixdim[0],32);
    _spacingX=pixdim[1];
    _spacingY=pixdim[2];
    _spacingZ=pixdim[3];

    //Read offset into a .nii file
    float vox_offset=0;
    f.read((char*)&vox_offset,4);

    //Read data
    f.seekg(std::streamoff(vox_offset));

    long nelem=dimx*dimy*dimz;
    Init(dimx,dimy,dimz);

    T* data=Data();

    switch (datatype){
    case 1: //bool, unsupported
        f.close(); 
        return false;
        break;

    case 2: //unsigned char
        {uint8_t* tmpdata=new uint8_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT8");
        f.read((char*)tmpdata, nelem*sizeof(uint8_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;}
        break;

    case 4: // signed short
        {int16_t* tmpdata=new int16_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "INT16");
        f.read((char*)tmpdata, nelem*sizeof(int16_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 8: // signed int
        {int32_t* tmpdata=new int32_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "INT32");
        f.read((char*)tmpdata, nelem*sizeof(int32_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 16: // float
        {float* tmpdata=new float[nelem];
        metadata.SetString("ORIGINAL_TYPE", "FLOAT");
        f.read((char*)tmpdata, nelem*sizeof(float));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 32: // complex, unsupported
        f.close();
        return false;
        break;
    case 64: // double
        {double* tmpdata=new double[nelem];
        metadata.SetString("ORIGINAL_TYPE", "DOUBLE");
        f.read((char*)tmpdata, nelem*sizeof(double));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 256: // signed char
        {int8_t* tmpdata=new int8_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "INT8");
        f.read((char*)tmpdata, nelem*sizeof(int8_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 512: // unsigned short
        {uint16_t* tmpdata=new uint16_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT16");
        f.read((char*)tmpdata, nelem*sizeof(uint16_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    case 768: // unsigned int
        {uint32_t* tmpdata=new uint32_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT32");
        f.read((char*)tmpdata, nelem*sizeof(uint32_t));
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
        }
        break;
    default: //unsupported types
        f.close();
        return false;
        break;
    }
    f.close();
    return true;
}

template <class T>
inline bool Grid3<T>::WriteNIFTI(std::string filename){
    //Not done yet
    ofstream outfile;
    outfile.open(filename);
    if( !outfile ) { return false; }

    //Write header size, always 348
    int32_t sizeof_hdr=348;
    outfile.write((char*)&sizeof_hdr,4);

    //Write image dimensions
    int16_t dim[8];
    dim[0]=3;
    dim[1]= GetWidth();
    dim[2]=GetHeight();
    dim[3]=GetDepth();
    dim[4]=0;
    dim[5]=0;
    dim[6]=0;
    dim[7]=0;

    outfile.seekp(40);
    outfile.write((char*)&dim[0],16);

    //Read data type
    int16_t datatype=0;

    string datatype_name=GetDataType();

    if(datatype_name=="UINT8"){datatype=2;}
    else if(datatype_name=="INT8"){datatype=256;}
    else if(datatype_name=="UINT16"){datatype=512;}
    else if(datatype_name=="INT16"){datatype=4;}
    else if(datatype_name=="UINT32"){datatype=768;}
    else if(datatype_name=="INT32"){datatype=8;}
    else if(datatype_name=="float"){datatype=16;}
    else if(datatype_name=="double"){datatype=64;}

    outfile.seekp(70);
    outfile.write((char*)&datatype,2);

    //Write grid spacing
    outfile.seekp(76);
    float pixdim[8]; 
    pixdim[0]=1; //This entry has a special meaning, unused by Grid3
    pixdim[1]=_spacingX;
    pixdim[2]=_spacingY;
    pixdim[3]=_spacingZ;
    outfile.write((char*)&pixdim[0],32);

    //Write offset to data
    float vox_offset=352;
    outfile.seekp(108);
    outfile.write((char*)&vox_offset,4);

    //Write "magic" string
    char magic[4]={'n','+','1','\0'}; //Single .nii file
    outfile.seekp(344);
    outfile.write((char*)&magic[0],4);

    //Write data
    outfile.seekp(352);
    T* data=Data();
    outfile.write((char*)data,sizeof(T)*NumberOfElements());

    outfile.close();
    return true;
}


template <class T>
inline bool Grid3<T>::ReadDICOM(std::string filename){
    #ifdef grid3_use_gdcm
    
    gdcm::ImageReader ir;
    ir.SetFileName(filename.c_str());
    if(!ir.Read()){ return false; }

    const gdcm::Image &image = ir.GetImage();

    const unsigned int* dim = image.GetDimensions();
    Init(dim[0], dim[1], dim[2], image.GetSpacing(0), image.GetSpacing(1), image.GetSpacing(2));

    gdcm::PixelFormat pf= image.GetPixelFormat();
    if(pf.GetSamplesPerPixel()!=1){return false;} //Only scalar images supported

    long nelem=NumberOfElements();
    string type= string(pf.GetScalarTypeAsString());
    
    //Defined data types for DICOM:
    //UINT8, 
    //INT8, 
    //UINT12, 
    //INT12, 
    //UINT16, 
    //INT16, 
    //UINT32, 
    //INT32, 
    //FLOAT16, 
    //FLOAT32, 
    //FLOAT64, 
    //SINGLEBIT, 
    //UNKNOWN 


    if(type=="UINT8"){
        uint8_t* tmpdata=new uint8_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT8");
        if( !image.GetBuffer( (char*)tmpdata) ){
            delete[] tmpdata; 
            return false; 
        }
        #pragma omp parallel for
        for(long i=0; i<nelem; i++){
            data[i]=T(tmpdata[i]); 
        }
        delete[] tmpdata;
    }
    else if(type=="UINT16"){
        uint16_t* tmpdata = new uint16_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT16");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
        
    }
    else if (type == "UINT32"){
        uint32_t* tmpdata = new uint32_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "UINT32");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
    }
    else if (type == "INT8"){
        int8_t* tmpdata = new int8_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "INT8");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;

    }
    else if (type == "INT16"){
        int16_t* tmpdata = new int16_t[nelem];
        metadata.SetString("ORIGINAL_TYPE", "INT16");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
    }
    else if (type == "INT32"){
        metadata.SetString("ORIGINAL_TYPE", "INT32");
        int32_t* tmpdata = new int32_t[nelem];
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
    }
    else if (type == "FLOAT16"){
        return false;   //unsupported
    }
    else if (type == "FLOAT32"){
        //We require the float type to have 32 bits
        if (CHAR_BIT*sizeof(float) != 32){  
            return false;
        }
        float* tmpdata = new float[nelem];
        metadata.SetString("ORIGINAL_TYPE", "FLOAT");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
    }
    else if (type == "FLOAT64"){
        //We require the double type to have 64 bits
        if (CHAR_BIT*sizeof(double) != 64){
            return false;
        }
        double* tmpdata = new double[nelem];
        metadata.SetString("ORIGINAL_TYPE", "DOUBLE");
        if (!image.GetBuffer((char*)tmpdata)){
            delete[] tmpdata;
            return false;
        }
        #pragma omp parallel for
        for (long i = 0; i<nelem; i++){
            data[i] = T(tmpdata[i]);
        }
        delete[] tmpdata;
    }

    return true;
    #endif
    
    return false;
}


#endif