File: fileio_mhd.cpp

package info (click to toggle)
odin 2.0.5-9
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 9,224 kB
  • sloc: cpp: 62,639; sh: 4,541; makefile: 779
file content (131 lines) | stat: -rw-r--r-- 4,645 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include "fileio.h"

/////////////////////////////////////////////////////////////

struct MhdFormat : public FileFormat {
  STD_string description() const {return "MetaImage";}
  svector suffix() const  {
    svector result; result.resize(1);
    result[0]="mhd";
    return result;
  }
  svector dialects() const {return svector();}

  int read(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol& prot) {
    Log<FileIO> odinlog("MhdFormat","read");
    STD_string mhdfile;
    if(::load(mhdfile,filename)<0) return -1;
    mhdfile=replaceStr(mhdfile,"=", " = "); // make sure key, value and '=' are separated by whitespaces
    svector toks=tokens(mhdfile);
    int ndims=-1;
    int i;
    int ntoks=toks.size();

    for(i=0; i<ntoks; i++) {
      if(toks[i]=="NDims" && i<(ntoks-2) && toks[i+1]=="=") {
        ndims=atoi(toks[i+2].c_str());
        break;
      }
    }
    ODINLOG(odinlog,normalDebug) << "ndims=" << ndims << STD_endl;
    if(ndims<0 || ndims>4) {
      ODINLOG(odinlog,errorLog) << "Invalid NDims=" << ndims << STD_endl;
      return -1;
    }

    fvector spacing(3);
    TinyVector<int,4> shape; shape=1;
    STD_string format;
    STD_string rawfname;
    for(i=0; i<ntoks; i++) {
      if(toks[i]=="DimSize" && i<(ntoks-1-ndims) && toks[i+1]=="=") {
        for(int idim=0; idim<ndims; idim++) {
          shape(3-idim)=atoi(toks[i+2+idim].c_str());
        }
      }
      if(toks[i]=="ElementSpacing" && i<(ntoks-1-ndims) && toks[i+1]=="=") {
        for(int idim=0; idim<ndims; idim++) {
          spacing[2-idim]=atof(toks[i+2+idim].c_str());
        }
      }
      if(toks[i]=="ElementType" && i<(ntoks-2) && toks[i+1]=="=") {
        STD_string eltype=toks[i+2];
        if(eltype=="MET_FLOAT")  format="float";
        if(eltype=="MET_DOUBLE") format="double";
        if(eltype=="MET_SHORT")  format="short";
        if(eltype=="MET_LONG")   format="long";
        if(format=="") {
          ODINLOG(odinlog,errorLog) << "Unrecognized ElementType=" << eltype << STD_endl;
          return -1;
        }
      }
      if(toks[i]=="ElementDataFile" && i<(ntoks-2) && toks[i+1]=="=") rawfname=toks[i+2];
    }
    ODINLOG(odinlog,normalDebug) << "spacing/format/rawfname=" << spacing.printbody() << "/" << format << "/" << rawfname << STD_endl;


    LDRfileName fname(filename);
    data.resize(shape);
    if(data.read(format, fname.get_dirname()+rawfname)<0) {
      ODINLOG(odinlog,errorLog) << "Unable to ElementDataFile =" << fname.get_dirname()+rawfname << STD_endl;
      return -1;
    }

    prot.geometry.set_sliceThickness(spacing[0]).set_sliceDistance(spacing[0]);
    prot.geometry.set_FOV(phaseDirection,float(shape(2))*spacing[1]);
    prot.geometry.set_FOV(readDirection, float(shape(3))*spacing[2]);
    return shape(0)*shape(1);
  }


  int write(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol& prot) {
    Log<FileIO> odinlog("MhdFormat","write");
    TinyVector<int,4> shape=data.shape();

    farray fa(data);
    ndim nn=fa.get_extent();
    nn.autosize(); // get rid of 1-size dims

    fvector res(nn.dim()); res=1.0;
    ODINLOG(odinlog,normalDebug) << "fa.dim()=" << fa.dim() << STD_endl;
    if(fa.dim()>2) {
      dvector soffset=prot.geometry.get_sliceOffsetVector();
      ODINLOG(odinlog,normalDebug) << "soffset=" << soffset.printbody() << STD_endl;
      if(soffset.size()>1) res[nn.dim()-3]=fabs(soffset[1]-soffset[0]);
      else res[nn.dim()-3]=prot.geometry.get_sliceThickness();
    }
    if(nn.dim()>1) res[nn.dim()-2]=secureDivision(prot.geometry.get_FOV(phaseDirection),prot.seqpars.get_MatrixSize(phaseDirection));
    if(nn.dim()>0) res[nn.dim()-1]=secureDivision(prot.geometry.get_FOV(readDirection), prot.seqpars.get_MatrixSize(readDirection));

    LDRfileName fname(filename);

    STD_string rawfname=fname.get_basename_nosuffix()+".raw";

    STD_string header;
    header+="NDims = "+itos(nn.dim())+"\n";
    header+="DimSize =";
    for(unsigned int idim=0; idim<nn.dim(); idim++) header+=" "+itos(nn[nn.dim()-idim-1]);
    header+="\n";
    header+="ElementType = MET_FLOAT\n";
    header+="ElementSpacing =";
    for(unsigned int idim=0; idim<nn.dim(); idim++) header+=" "+ftos(res[nn.dim()-idim-1]);
    header+="\n";
    header+="ElementByteOrderMSB = False\n";
    header+="ElementDataFile = "+rawfname+"\n";

    if(::write(header,filename)<0) return -1;
    return data.write<float>(fname.get_dirname()+rawfname);
  }


};



//////////////////////////////////////////////////////////////

void register_mhd_format() {
 static MhdFormat mf;
 mf.register_format();
}