File: fileio_ismrmrd.cpp

package info (click to toggle)
odin 2.0.5-8
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 9,196 kB
  • sloc: cpp: 62,638; sh: 4,541; makefile: 779
file content (160 lines) | stat: -rw-r--r-- 4,781 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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#include "fileio.h"

#ifdef ISMRMRDSUPPORT

#include <ismrmrd/dataset.h>
#include <ismrmrd/ismrmrd.h>

#define HDF_DATASET_STR "dataset"
#define HDF_IMAGE_STR "image"

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

struct IsmrmrdFormat : public FileFormat{

  STD_string description() const {return "ISMRMRD Image";}

  svector suffix() const {
    svector result; result.resize(1);
    result[0]="h5";
    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("IsmrmrdFormat","read");

    ISMRMRD::Dataset dataset(filename.c_str(), HDF_DATASET_STR, false);

    int nimages=dataset.getNumberOfImages(HDF_IMAGE_STR);
    ODINLOG(odinlog,normalDebug) << "nimages=" << nimages << STD_endl;

    int nread=-1;
    int nphase=-1;
    int nslice=-1;

    Geometry& geometry=prot.geometry;

    ISMRMRD::Image<float> image;
    for(int i=0; i<nimages; i++) {
      dataset.readImage(HDF_IMAGE_STR, i, image);

      if(i==0) {
        nread =image.getMatrixSizeX();
        nphase=image.getMatrixSizeY();
        nslice=image.getMatrixSizeZ();
        ODINLOG(odinlog,normalDebug) << "nslice/nphase/nread=" << nslice << "/" << nphase << "/" << nread << STD_endl;

        data.resize(nimages, nslice, nphase, nread);

        geometry.set_FOV(readDirection,image.getFieldOfViewX());
        geometry.set_FOV(phaseDirection,image.getFieldOfViewY());
        geometry.set_FOV(sliceDirection,image.getFieldOfViewZ());

        float slicethick=secureDivision(image.getFieldOfViewZ(), nslice);
        geometry.set_sliceThickness(slicethick);
        geometry.set_sliceDistance(slicethick);

        dvector readvec(3);
        readvec[0]=image.getReadDirectionX();
        readvec[1]=image.getReadDirectionY();
        readvec[2]=image.getReadDirectionZ();

        dvector phasevec(3);
        phasevec[0]=image.getPhaseDirectionX();
        phasevec[1]=image.getPhaseDirectionY();
        phasevec[2]=image.getPhaseDirectionZ();

        dvector slicevec(3);
        slicevec[0]=image.getSliceDirectionX();
        slicevec[1]=image.getSliceDirectionY();
        slicevec[2]=image.getSliceDirectionZ();

        dvector centervec(3);
        centervec[0]=image.getPositionX();
        centervec[1]=image.getPositionY();
        centervec[2]=image.getPositionZ();

        geometry.set_orientation_and_offset(readvec, phasevec, slicevec, centervec);

      } else {
        if(nread!=image.getMatrixSizeX() || nphase!=image.getMatrixSizeY() || nslice!=image.getMatrixSizeZ()) {
          ODINLOG(odinlog,errorLog) << "size mismatch nslice/nphase/nread=" << nslice << "/" << nphase << "/" << nread << STD_endl;
          return -1;
        }
      }

      for(int iz=0; iz<nslice; iz++) {
        for(int iy=0; iy<nphase; iy++) {
          for(int ix=0; ix<nread; ix++) {
            data(i,iz,iy,ix)=image(ix,iy,iz);
          }
        }
      }

    }

    return nimages*nslice;
  }



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

    rmfile(filename.c_str()); // remove old file

    ISMRMRD::Dataset dataset(filename.c_str(), HDF_DATASET_STR, true);

    const TinyVector<int, 4> shape(data.shape());

    ISMRMRD::Image<float> image(shape(readDim), shape(phaseDim), shape(sliceDim), 1); // only one channel

    image.setImageType(ISMRMRD::ISMRMRD_IMTYPE_MAGNITUDE);

    const Geometry& geometry=prot.geometry;

    image.setFieldOfView(geometry.get_FOV(readDirection), geometry.get_FOV(phaseDirection), geometry.get_FOV(sliceDirection));

    dvector center(geometry.get_center());
    image.setPosition(center[0], center[1], center[2]);

    dvector readvec(geometry.get_readVector());
    image.setReadDirection(readvec[0], readvec[1], readvec[2]);

    dvector phasevec(geometry.get_phaseVector());
    image.setPhaseDirection(phasevec[0], phasevec[1], phasevec[2]);

    dvector slicevec(geometry.get_sliceVector());
    image.setSliceDirection(slicevec[0], slicevec[1], slicevec[2]);

    for(int itime=0; itime<shape(timeDim); itime++) {

      for(int iz=0; iz<shape(sliceDim); iz++) {
        for(int iy=0; iy<shape(phaseDim); iy++) {
          for(int ix=0; ix<shape(readDim); ix++) {
            image(ix,iy,iz) = data(itime,iz,iy,ix);
          }
        }
      }
      image.setRepetition(itime);

      dataset.appendImage(HDF_IMAGE_STR, image);

    }

    return 1;
  }

};
#endif // ISMRMRDSUPPORT

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

void register_ismrmrd_format() {
#ifdef ISMRMRDSUPPORT
  static IsmrmrdFormat irf;
  irf.register_format();
#endif
}