1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #include "itkSiemensVisionImageIO.h"
19 #include "itksys/SystemTools.hxx"
20 #include <iostream>
21 #include <fstream>
22 #include <algorithm>
23
24 //From uiig library "The University of Iowa Imaging Group-UIIG"
25
26 namespace itk
27 {
28 // Default constructor
SiemensVisionImageIO()29 SiemensVisionImageIO::SiemensVisionImageIO()
30 {
31 //Purposefully left blank
32 }
33
~SiemensVisionImageIO()34 SiemensVisionImageIO::~SiemensVisionImageIO()
35 {
36 //Purposefully left blank
37 }
38
CanReadFile(const char * FileNameToRead)39 bool SiemensVisionImageIO::CanReadFile(const char *FileNameToRead)
40 {
41 this->SetFileName(FileNameToRead);
42 //
43 // Can you open it?
44 std::ifstream f;
45 try
46 {
47 this->OpenFileForReading( f, FileNameToRead );
48 }
49 catch( ExceptionObject & )
50 {
51 return false;
52 }
53 int matrixX;
54 //
55 // another lame heuristic, check the actual file size against
56 // the image size suggested in header + the header size.
57 if ( this->GetIntAt(f, HDR_DISPLAY_SIZE, &matrixX, false) != 0 )
58 {
59 return false;
60 }
61
62 if ( ( HDR_TOTAL_LENGTH + ( matrixX * matrixX * 2 ) ) !=
63 (int)itksys::SystemTools::FileLength(FileNameToRead) )
64 {
65 return false;
66 }
67
68 return true;
69 }
70
ReadHeader(const char * FileNameToRead)71 GEImageHeader * SiemensVisionImageIO::ReadHeader(const char *FileNameToRead)
72 {
73 if ( !this->CanReadFile(FileNameToRead) )
74 {
75 RAISE_EXCEPTION();
76 }
77 int tmpInt;
78 double tmpDble;
79
80 // #define DEBUGHEADER
81 #if defined( DEBUGHEADER )
82 #define DB(x) std::cerr << #x << " " << x << std::endl
83 #else
84 #define DB(x)
85 #endif
86
87 #define GE_PROD_STR "SIEMENS"
88 #define TEMPLEN 2048
89 auto * hdr = new GEImageHeader;
90 if ( hdr == nullptr )
91 {
92 RAISE_EXCEPTION();
93 }
94 #if defined( DEBUGHEADER )
95 std::cerr << "----------------------" << FileNameToRead << "----------------------" << std::endl;
96 #endif
97
98 std::ifstream f;
99 this->OpenFileForReading( f, FileNameToRead );
100
101 sprintf (hdr->scanner, "GE-ADW");
102
103 // Set modality to UNKNOWN
104 strcpy(hdr->modality, "UNK");
105
106 strncpy(hdr->filename, FileNameToRead, itk::IOCommon::ITK_MAXPATHLEN);
107
108 // Get VITAL Information from the header
109 this->GetStringAt(f, HDR_PAT_ID, hdr->patientId, HDR_PAT_ID_LEN);
110 hdr->patientId[HDR_PAT_ID_LEN] = '\0';
111 DB(hdr->patientId);
112 // fprintf(stderr, "Patient %s\n", hdr->patientId);a
113
114 this->GetStringAt(f, HDR_PAT_NAME, hdr->name, HDR_PAT_NAME_LEN);
115 hdr->name[HDR_PAT_NAME_LEN] = '\0';
116 DB(hdr->name);
117
118 int year, month, day, hour, minute, second;
119
120 this->GetIntAt(f, HDR_REG_YEAR, &year);
121
122 this->GetIntAt(f, HDR_REG_MONTH, &month);
123
124 this->GetIntAt(f, HDR_REG_DAY, &day);
125
126 this->GetIntAt(f, HDR_REG_HOUR, &hour);
127
128 this->GetIntAt(f, HDR_REG_MIN, &minute);
129
130 this->GetIntAt(f, HDR_REG_SEC, &second);
131
132 sprintf (hdr->date, "%d/%d/%d %d:%d:%d", year, month, day, hour, minute, second);
133 DB(hdr->date);
134
135 this->GetStringAt(f, HDR_INSTUTE_NAME, hdr->hospital, HDR_INSTUTE_NAME_LEN);
136 hdr->hospital[HDR_INSTUTE_NAME_LEN] = '\0';
137 DB(hdr->hospital);
138
139 this->GetStringAt(f, HDR_MODEL_NAME, hdr->scanner, HDR_MODEL_NAME_LEN);
140 hdr->scanner[HDR_MODEL_NAME_LEN] = '\0';
141 DB(hdr->scanner);
142 for ( unsigned int i = 0; i < strlen(hdr->scanner); i++ )
143 {
144 if ( hdr->scanner[i] == ' ' ) { hdr->scanner[i] = '-'; }
145 }
146
147 char tmpStr[TEMPLEN];
148 char tmpStr2[TEMPLEN];
149
150 this->GetStringAt(f, TEXT_STUDY_NUM2, tmpStr, TEXT_STUDY_NUM2_LEN);
151 tmpStr[TEXT_STUDY_NUM2_LEN] = '\0';
152 hdr->seriesNumber = std::stoi(tmpStr);
153 DB(hdr->seriesNumber);
154
155 this->GetStringAt(f, TEXT_IMG_NUMBER, tmpStr, TEXT_IMG_NUMBER_LEN);
156 tmpStr[TEXT_IMG_NUMBER_LEN] = '\0';
157 hdr->imageNumber = std::stoi(tmpStr);
158 DB(hdr->imageNumber);
159
160 this->GetStringAt(f, TEXT_SLICE_THCK, tmpStr, TEXT_SLICE_THCK_LEN);
161 tmpStr[TEXT_SLICE_THCK_LEN] = '\0';
162 hdr->sliceThickness = std::stoi(tmpStr);
163 hdr->sliceGap = 0.0f;
164
165 DB(hdr->sliceThickness);
166
167 this->GetIntAt( f, HDR_DISPLAY_SIZE, &tmpInt, sizeof( int ) );
168 hdr->imageXsize = (int)tmpInt;
169 DB(hdr->imageXsize);
170 hdr->imageYsize = (int)tmpInt;
171 DB(hdr->imageYsize);
172
173 this->GetStringAt(f, TEXT_ACQ_MTRX_FREQ, tmpStr, TEXT_ACQ_MTRX_FREQ_LEN);
174 tmpStr[TEXT_ACQ_MTRX_FREQ_LEN] = '\0';
175 hdr->acqXsize = std::stoi(tmpStr);
176 DB(hdr->acqXsize);
177
178 this->GetStringAt(f, TEXT_ACQ_MTRX_PHASE, tmpStr, TEXT_ACQ_MTRX_PHASE_LEN);
179 tmpStr[TEXT_ACQ_MTRX_PHASE_LEN] = '\0';
180 hdr->acqYsize = std::stoi(tmpStr);
181 DB(hdr->acqYsize);
182
183 this->GetStringAt(f, TEXT_FOVH, tmpStr, TEXT_FOVH_LEN);
184 tmpStr[TEXT_FOVH_LEN] = '\0';
185 hdr->xFOV = static_cast< float >( std::stod(tmpStr) );
186 DB(hdr->xFOV);
187
188 this->GetStringAt(f, TEXT_FOVV, tmpStr, TEXT_FOVV_LEN);
189 tmpStr[TEXT_FOVV_LEN] = '\0';
190 hdr->yFOV = static_cast< float >( std::stod(tmpStr) );
191 DB(hdr->yFOV);
192
193 this->GetDoubleAt( f, HDR_PIXELSIZE_ROW, &tmpDble, sizeof( double ) );
194 hdr->imageXres = (float)tmpDble;
195 DB(hdr->imageXres);
196
197 this->GetDoubleAt( f, HDR_PIXELSIZE_CLMN, &tmpDble, sizeof( double ) );
198 hdr->imageYres = (float)tmpDble;
199 DB(hdr->imageYres);
200
201 this->GetStringAt(f, TEXT_ANGLE_FLAG1, tmpStr, TEXT_ANGLE_FLAG1_LEN);
202 tmpStr[TEXT_ANGLE_FLAG1_LEN] = '\0';
203
204 this->GetStringAt(f, TEXT_ANGLE_FLAG3, tmpStr2, TEXT_ANGLE_FLAG3_LEN);
205 tmpStr2[TEXT_ANGLE_FLAG3_LEN] = '\0';
206
207 std::string text_angle_len;
208 {
209 char tmpStr3[TEMPLEN];
210 this->GetStringAt(f, TEXT_ANGLE, tmpStr3, TEXT_ANGLE_LEN);
211 tmpStr3[TEXT_ANGLE_LEN] = '\0';
212 text_angle_len = tmpStr3;
213 }
214 // An empty string implies an angle less than 45 degrees for backwards compatibility
215 text_angle_len.erase(
216 std::remove_if(text_angle_len.begin(), text_angle_len.end(),isspace),
217 text_angle_len.end() ); //Remove all whitespace
218 if ( strcmp(tmpStr, "Cor") == 0 )
219 {
220 if ( text_angle_len.empty() || std::fabs( std::stod(text_angle_len) ) <= 45.0 )
221 {
222 //hdr->imagePlane = itk::IOCommon::ITK_ANALYZE_ORIENTATION_IRP_CORONAL;
223 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
224 }
225 else
226 {
227 if ( strcmp(tmpStr2, "Sag") == 0 )
228 {
229 //hdr->imagePlane =
230 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_SAGITTAL;
231 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
232 }
233 else
234 {
235 //hdr->imagePlane =
236 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_TRANSVERSE;
237 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
238 }
239 }
240 }
241 else if ( strcmp(tmpStr, "Sag") == 0 )
242 {
243 if ( text_angle_len.empty() || std::fabs( std::stod(text_angle_len) ) <= 45.0 )
244 {
245 //hdr->imagePlane =
246 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_SAGITTAL;
247 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
248 }
249 else
250 {
251 if ( strcmp(tmpStr2, "Cor") == 0 )
252 {
253 //hdr->imagePlane =
254 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_CORONAL;
255 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
256 }
257 else
258 {
259 //hdr->imagePlane =
260 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_TRANSVERSE;
261 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
262 }
263 }
264 }
265 else
266 {
267 if ( text_angle_len.empty() || std::fabs( std::stod(text_angle_len) ) <= 45.0 )
268 {
269 //hdr->imagePlane =
270 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_TRANSVERSE;
271 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
272 }
273 else
274 {
275 if ( strcmp(tmpStr2, "Cor") == 0 )
276 {
277 //hdr->imagePlane =
278 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_CORONAL;
279 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
280 }
281 else
282 {
283 //hdr->imagePlane =
284 // itk::SpatialOrientation::ITK_ANALYZE_ORIENTATION_IRP_SAGITTAL;
285 hdr->coordinateOrientation = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
286 }
287 }
288 }
289
290 /* fprintf(stderr, "Plane %d\n", hdr->imagePlane); */
291 this->GetStringAt(f, TEXT_SLICE_POS, tmpStr, TEXT_SLICE_POS_LEN);
292 tmpStr[TEXT_SLICE_POS_LEN] = '\0';
293 hdr->sliceLocation = static_cast< float >( std::stod(tmpStr) );
294 DB(hdr->sliceLocation);
295
296 /* fprintf(stderr, "Slice Location %f\n", hdr->sliceLocation); */
297 this->GetDoubleAt( f, HDR_TR, &tmpDble, sizeof( double ) );
298 hdr->TR = (float)tmpDble / 1000.0f;
299 DB(hdr->TR);
300
301 this->GetDoubleAt( f, HDR_TE + 8, &tmpDble, sizeof( double ) );
302 hdr->TI = (float)tmpDble / 1000.0f;
303 DB(hdr->TI);
304
305 this->GetDoubleAt( f, HDR_TE, &tmpDble, sizeof( double ) );
306 hdr->TE = (float)tmpDble / 1000.0f;
307 DB(hdr->TE);
308
309 this->GetStringAt(f, TEXT_ECHO_NUM, tmpStr, TEXT_ECHO_NUM_LEN);
310 tmpStr[TEXT_ECHO_NUM_LEN] = '\0';
311 hdr->echoNumber = (int)std::stoi(tmpStr);
312 DB(hdr->echoNumber);
313
314 this->GetDoubleAt( f, HDR_FLIP_ANGLE, &tmpDble, sizeof( double ) );
315 hdr->flipAngle = (int)tmpDble;
316 DB(hdr->flipAngle);
317
318 this->GetStringAt(f, HDR_SEQPROG_NAME, hdr->pulseSequence, HDR_SEQPROG_NAME_LEN);
319 hdr->pulseSequence[HDR_SEQPROG_NAME_LEN] = '\0';
320
321 hdr->offset = HDR_TOTAL_LENGTH;
322 return hdr;
323 }
324 } // end namespace itk
325