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