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 "itkGE5ImageIO.h"
19 #include "itkByteSwapper.h"
20 #include "itksys/SystemTools.hxx"
21 #include "Ge5xHdr.h"
22 #include <algorithm>
23 #include <iostream>
24 #include <fstream>
25
26 #include "vnl/vnl_cross.h"
27
28
29 //From uiig library "The University of Iowa Imaging Group-UIIG"
30
31 namespace itk
32 {
33 static constexpr char GE_PROD_STR[] = "SIGNA";
34 // Default constructor
35 GE5ImageIO::GE5ImageIO() = default;
36
~GE5ImageIO()37 GE5ImageIO::~GE5ImageIO()
38 {
39 //Purposefully left blank
40 }
41
42 int GE5ImageIO
CheckGE5xImages(char const * const imageFileTemplate,std::string & reason)43 ::CheckGE5xImages(char const *const imageFileTemplate, std::string & reason)
44 {
45 //
46 // Does it exist?
47 if ( !itksys::SystemTools::FileExists(imageFileTemplate) )
48 {
49 reason = "File does not exist";
50 return -1;
51 }
52 //
53 // is it at least 5000 bytes?
54 if ( itksys::SystemTools::FileLength(imageFileTemplate) < 5000 )
55 {
56 reason = "File size is less than 5000 bytes";
57 return -1;
58 }
59
60 std::ifstream f;
61 try
62 {
63 this->OpenFileForReading( f, imageFileTemplate );
64 }
65 catch( ExceptionObject & )
66 {
67 reason = "File could not be opened for read";
68 return -1;
69 }
70
71 Ge5xPixelHeader imageHdr; /* Header Structure for GE 5x images
72 */
73 char hdr[GENESIS_SU_HDR_LEN]; /* Header to hold GE Suite header */
74 char prod[16]; /* Product name from Suite Header */
75
76 // First pass see if image is a raw MR extracted via ximg
77 if ( !this->ReadBufferAsBinary( f, (void *)&imageHdr, sizeof( imageHdr ) ) )
78 {
79 f.close();
80 return -1;
81 }
82 ByteSwapper< int >::SwapFromSystemToBigEndian(&imageHdr.GENESIS_IH_img_magic);
83 if ( imageHdr.GENESIS_IH_img_magic == GE_5X_MAGIC_NUMBER )
84 {
85 f.close();
86 return 0;
87 }
88 f.seekg(0, std::ios::beg);
89
90 //
91 // Second pass see if image was extracted via tape by Gene's tape
92 // reading software.
93 //
94 if ( !this->ReadBufferAsBinary(f, (void *)hdr, GENESIS_SU_HDR_LEN) )
95 {
96 reason = "Failed to read study header";
97 f.close();
98 return -1;
99 }
100 strncpy (prod, hdr + GENESIS_SU_PRODID, 13);
101 prod[13] = '\0';
102 if ( strcmp (prod, GE_PROD_STR) == 0 )
103 {
104 f.close();
105 return 0;
106 }
107
108 reason = "Failed to find string SIGNA";
109 f.close();
110 return -1;
111 }
112
CanReadFile(const char * FileNameToRead)113 bool GE5ImageIO::CanReadFile(const char *FileNameToRead)
114 {
115 std::string reason;
116
117 return this->CheckGE5xImages(FileNameToRead, reason) == 0 ? true : false;
118 }
119
120 namespace
121 {
122 void
SwapPixHdr(Ge5xPixelHeader * hdr)123 SwapPixHdr(Ge5xPixelHeader *hdr)
124 {
125 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_magic ) );
126 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_hdr_length ) );
127 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_width ) );
128 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_height ) );
129 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_depth ) );
130 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_compress ) );
131 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_dwindow ) );
132 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_dlevel ) );
133 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_bgshade ) );
134 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_ovrflow ) );
135 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_undflow ) );
136 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_top_offset ) );
137 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_bot_offset ) );
138 ByteSwapper< short >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_version ) );
139 ByteSwapper< unsigned short >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_checksum ) );
140 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_id ) );
141 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_id ) );
142 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_unpack ) );
143 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_unpack ) );
144 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_compress ) );
145 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_compress ) );
146 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_histo ) );
147 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_histo ) );
148 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_text ) );
149 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_text ) );
150 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_graphics ) );
151 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_graphics ) );
152 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_dbHdr ) );
153 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_dbHdr ) );
154 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_levelOffset ) );
155 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_user ) );
156 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_user ) );
157 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_suite ) );
158 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_suite ) );
159 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_exam ) );
160 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_exam ) );
161 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_series ) );
162 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_series ) );
163 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_p_image ) );
164 ByteSwapper< int >::SwapFromSystemToBigEndian ( &( hdr->GENESIS_IH_img_l_image ) );
165 }
166 }
167
168 GEImageHeader *
ReadHeader(const char * FileNameToRead)169 GE5ImageIO::ReadHeader(const char *FileNameToRead)
170 {
171 //#define VERBOSE_DEBUGGING
172 #if defined( VERBOSE_DEBUGGING )
173 #define RGEDEBUG(x) x
174 #else
175 #define RGEDEBUG(x)
176 #endif
177
178 Ge5xPixelHeader imageHdr; // GE 5x Header
179 GEImageHeader * curImage;
180 bool pixelHdrFlag;
181 std::string reason;
182 if ( this->CheckGE5xImages(FileNameToRead, reason) != 0 )
183 {
184 itkExceptionMacro(
185 "GE5ImageIO could not open file "
186 << FileNameToRead << " for reading."
187 << std::endl
188 << "Reason: "
189 << reason
190 );
191 }
192
193 curImage = new GEImageHeader;
194 if ( curImage == nullptr )
195 {
196 itkExceptionMacro(
197 "GE5ImageIO failed to create a GEImageHeader while reading "
198 << FileNameToRead << " ." << std::endl << "Reason: "
199 << "new GEImageHeader failed."
200 );
201 }
202 memset( curImage, 0, sizeof( GEImageHeader ) );
203
204 std::ifstream f;
205 this->OpenFileForReading( f, FileNameToRead );
206
207 f.read( (char *)&imageHdr, sizeof( imageHdr ) );
208 if ( f.fail() )
209 {
210 itkExceptionMacro(
211 "GE5ImageIO IO error while reading "
212 << FileNameToRead << " ." << std::endl
213 << "Reason: " << itksys::SystemTools::GetLastSystemError()
214 );
215 }
216 SwapPixHdr(&imageHdr);
217
218 // NOTE: The handling of version 2 vs Version3 is modelled after
219 // the sivic GE5Signa5x reader -- found here
220 // http://sivic.svn.sourceforge.net
221 //
222 // if i didn't see it with my own eyes I wouldn't believe it.
223 // Signa 5x either have a proper header or they don't!
224 // Appparently Version 2 files do always have a header or
225 // we'd be totally lost!
226 // if they don't have a header, we have to make assumptions
227 // about where they start and hope we're right; below, the offset
228 // is computed once the X & Y dims are known
229 if ( imageHdr.GENESIS_IH_img_magic == GE_5X_MAGIC_NUMBER )
230 {
231 pixelHdrFlag = true;
232 curImage->offset = imageHdr.GENESIS_IH_img_hdr_length;
233 }
234 else
235 {
236 pixelHdrFlag = false;
237 }
238 strncpy (curImage->filename, FileNameToRead, IOCommon::ITK_MAXPATHLEN);
239
240 //
241 // if there's no GE5 header on the file we have to assume
242 // it's Version 3. If there is a header, it could be version2
243 // in which case we need to fill out the fields in the header
244 // that are defined in version 3 and not version 2.
245 if(pixelHdrFlag && imageHdr.GENESIS_IH_img_version == 2)
246 {
247 imageHdr.GENESIS_IH_img_p_suite = 124; // Version 3 is 2304
248 imageHdr.GENESIS_IH_img_l_suite = 116; // Version 3 is 114
249 imageHdr.GENESIS_IH_img_p_exam = 240; // Version 3 is 2418
250 imageHdr.GENESIS_IH_img_l_exam = 1040; // Version 3 is 1024
251 imageHdr.GENESIS_IH_img_p_series = 1280; // Version 3 is 3442
252 imageHdr.GENESIS_IH_img_l_series = 1028; // Version 3 is 1020
253 imageHdr.GENESIS_IH_img_p_image = 2308; // Version 3 is 4462
254 imageHdr.GENESIS_IH_img_l_image = 1044; // Don't know for sure?
255 }
256
257 // if we have a version2 file, most of the fields are offset from
258 // their version3 positions.
259 #define VOff(a,b) (imageHdr.GENESIS_IH_img_version != 2 ? a : b)
260 // Create a buffer to read the exam header.
261 // Now seek to the exam header and read the data into the buffer.
262 char *buffer = nullptr;
263 if(pixelHdrFlag)
264 {
265 buffer = new char[imageHdr.GENESIS_IH_img_l_exam];
266 if ( buffer == nullptr )
267 {
268 f.close();
269 itkExceptionMacro("GE5ImageIO:Unable to allocate memory for exam header!");
270 }
271 f.seekg(imageHdr.GENESIS_IH_img_p_exam,std::ios::beg);
272 f.read(buffer,imageHdr.GENESIS_IH_img_l_exam);
273 }
274 else
275 {
276 buffer = new char[GENESIS_EX_HDR_LEN];
277 if ( buffer == nullptr )
278 {
279 f.close();
280 itkExceptionMacro("GE5ImageIO:Unable to allocate memory for exam header!");
281 }
282 f.seekg(GENESIS_EX_HDR_START,std::ios::beg);
283 f.read(buffer,GENESIS_EX_HDR_LEN);
284 }
285 if(f.fail())
286 {
287 f.close();
288 delete[] buffer;
289 itkExceptionMacro("GE5ImageIO:Could not read exam header!");
290 }
291
292 // Now extract the exam information from the buffer.
293 curImage->examNumber = hdr2Short(buffer+8);
294
295 strncpy(curImage->hospital,buffer+10,34);
296 curImage->hospital[34] = '\0';
297
298 // patient id
299 std::string tmpId(buffer+VOff(84,88), 13);
300 std::remove(tmpId.begin(), tmpId.end(), '-');
301 strncpy(curImage->patientId, tmpId.c_str(), sizeof(curImage->patientId)-1);
302 curImage->patientId[sizeof(curImage->patientId)-1] = '\0';
303
304 strncpy(curImage->name,buffer+VOff(97,101),25);
305 curImage->name[24] = '\0';
306
307 // Need to know modality as well.
308 strncpy(curImage->modality,buffer+VOff(305,309),3);
309 curImage->modality[3] = '\0';
310 bool isCT = false;
311 if( strncmp(curImage->modality, "CT", 2) == 0)
312 {
313 isCT = true;
314 }
315
316 // Done with exam, delete buffer.
317 delete[] buffer;
318 buffer = nullptr;
319
320 // Allocate buffer for series header.
321 // Now seek to the series header and read the data into the buffer.
322 if(pixelHdrFlag)
323 {
324 buffer = new char[imageHdr.GENESIS_IH_img_l_series];
325 if ( buffer == nullptr )
326 {
327 f.close();
328 itkExceptionMacro
329 ("GE5ImageIO:Unable to allocate memory for series header!");
330 }
331 f.seekg(imageHdr.GENESIS_IH_img_p_series, std::ios::beg);
332 f.read(buffer, imageHdr.GENESIS_IH_img_l_series);
333 }
334 else
335 {
336 buffer = new char[GENESIS_SE_HDR_LEN];
337 if ( buffer == nullptr )
338 {
339 f.close();
340 itkExceptionMacro
341 ("GE5ImageIO:Unable to allocate memory for series header!");
342 }
343 f.seekg(GENESIS_SE_HDR_START);
344 f.read(buffer,GENESIS_SE_HDR_LEN);
345 }
346 if(f.fail())
347 {
348 f.close();
349 itkExceptionMacro("GE5ImageIO:Could not read exam header!");
350 }
351
352 // Now extract the series information from the buffer.
353 curImage->seriesNumber = hdr2Short(buffer+10);
354
355 int timeStamp = hdr2Int(buffer+12);
356 statTimeToAscii(&timeStamp,curImage->date,sizeof(curImage->date));
357
358 // Done with series, delete buffer and allocate for MR header.
359 delete[] buffer;
360 buffer = nullptr;
361 if(pixelHdrFlag)
362 {
363 buffer = new char[imageHdr.GENESIS_IH_img_l_image];
364 if ( buffer == nullptr )
365 {
366 f.close();
367 itkExceptionMacro("GE5ImageIO:Unable to allocate memory for MR header!");
368 }
369 // Now seek to the MR header and read the data into the buffer.
370 f.seekg(imageHdr.GENESIS_IH_img_p_image, std::ios::beg);
371 f.read(buffer,imageHdr.GENESIS_IH_img_l_image);
372 }
373 else
374 {
375 buffer = new char[GENESIS_MR_HDR_LEN];
376 if ( buffer == nullptr )
377 {
378 f.close();
379 itkExceptionMacro("GE5ImageIO:Unable to allocate memory for MR header!");
380 }
381 f.seekg(GENESIS_IM_HDR_START, std::ios::beg);
382 f.read(buffer,GENESIS_MR_HDR_LEN);
383 }
384 if(f.fail())
385 {
386 itkExceptionMacro("GE5ImageIOCould not read exam header!");
387 }
388 // Won't need anymore info from the file after this, so close file.
389 f.close();
390
391 // Now extract the MR information from the buffer.
392 // This is the largest header!
393 curImage->imageNumber = hdr2Short(buffer+12);
394
395 curImage->sliceThickness = hdr2Float(buffer+VOff(26,28));
396
397 curImage->imageXsize = hdr2Short(buffer+VOff(30,32));
398 curImage->imageYsize = hdr2Short(buffer+VOff(32,34));
399 //
400 // if this a headerless flag, we don't know until now
401 // where to begin reading image data
402 if(!pixelHdrFlag)
403 {
404 curImage->offset =
405 itksys::SystemTools::FileLength(FileNameToRead)
406 - ( curImage->imageXsize * curImage->imageYsize * 2 );
407 }
408
409 curImage->xFOV = hdr2Float(buffer+VOff(34,36));
410 curImage->yFOV = hdr2Float(buffer+VOff(38,40));
411
412 curImage->acqXsize = hdr2Short(buffer+VOff(42,44));
413 curImage->acqYsize = hdr2Short(buffer+VOff(46,48));
414
415 curImage->imageXres = hdr2Float(buffer+VOff(50,52));
416 curImage->imageYres = hdr2Float(buffer+VOff(54,56));
417
418 short int GE_Plane(hdr2Short(buffer+VOff(114,116)));
419 switch ( GE_Plane )
420 {
421 case GE_CORONAL:
422 curImage->coordinateOrientation =
423 itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
424 break;
425 case GE_SAGITTAL:
426 curImage->coordinateOrientation =
427 itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR;
428 break;
429 case GE_AXIAL:
430 curImage->coordinateOrientation =
431 itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;
432 break;
433 default:
434 curImage->coordinateOrientation =
435 itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP;
436 break;
437 }
438
439
440 curImage->sliceLocation = hdr2Float(buffer+VOff(126,132));
441
442 curImage->centerR = hdr2Float(buffer+VOff(130,136));
443 curImage->centerA = hdr2Float(buffer+VOff(134,140));
444 curImage->centerS = hdr2Float(buffer+VOff(138,144));
445 curImage->normR = hdr2Float(buffer+VOff(142,146));
446 curImage->normA = hdr2Float(buffer+VOff(146,152));
447 curImage->normS = hdr2Float(buffer+VOff(150,156));
448 curImage->tlhcR = hdr2Float(buffer+VOff(154,160));
449 curImage->tlhcA = hdr2Float(buffer+VOff(158,164));
450 curImage->tlhcS = hdr2Float(buffer+VOff(162,168));
451 curImage->trhcR = hdr2Float(buffer+VOff(166,172));
452 curImage->trhcA = hdr2Float(buffer+VOff(170,176));
453 curImage->trhcS = hdr2Float(buffer+VOff(174,180));
454 curImage->brhcR = hdr2Float(buffer+VOff(178,184));
455 curImage->brhcA = hdr2Float(buffer+VOff(182,188));
456 curImage->brhcS = hdr2Float(buffer+VOff(186,192));
457
458 // These values are all MR specific
459 if(!isCT)
460 {
461 curImage->TR = hdr2Int(buffer+VOff(194,200));
462 curImage->TI = hdr2Int(buffer+VOff(198,204));
463 curImage->TE = hdr2Int(buffer+VOff(202,208));
464 curImage->TE2 = hdr2Int(buffer+VOff(206,212));
465
466 if((curImage->numberOfEchoes = hdr2Short(buffer+VOff(210,216))) == 0)
467 {
468 curImage->numberOfEchoes = 1;
469 }
470
471 curImage->echoNumber = hdr2Short(buffer+VOff(212,218));
472
473 curImage->NEX = hdr2Int(buffer+VOff(218,224));
474
475 curImage->flipAngle = hdr2Short(buffer+VOff(254,260));
476
477 strncpy(curImage->pulseSequence,
478 buffer+VOff(308,320),
479 34);
480 curImage->pulseSequence[33] = '\0';
481
482 curImage->numberOfSlices = hdr2Short(buffer+VOff(398,416));
483 }
484 else
485 {
486 curImage->TR = 0;
487 curImage->TI = 0;
488 curImage->TE = 0;
489 curImage->TE2 = 0;
490 curImage->numberOfEchoes = 1;
491 curImage->echoNumber = 1;
492 curImage->NEX = 1;
493 curImage->flipAngle = 0;
494 curImage->pulseSequence[0] = '\0';
495 curImage->numberOfSlices = 1;
496 }
497
498 // Delete the buffer and return the pointer to the header.
499 // The function that receives the pointer must do memory
500 // cleanup or a memory leak will occur.
501 delete[] buffer;
502
503 return ( curImage );
504 }
505
506 void
ModifyImageInformation()507 GE5ImageIO::ModifyImageInformation()
508 {
509 vnl_vector< double > dirx(3), diry(3), dirz(3);
510
511 // NOTE: itk use LPS coordinates while the GE system uses RAS
512 // coordinates. Consequently, the R and A coordinates must be negated
513 // to convert them to L and P.
514
515 dirx[0] = -( m_ImageHeader->trhcR - m_ImageHeader->tlhcR );
516 dirx[1] = -( m_ImageHeader->trhcA - m_ImageHeader->tlhcA );
517 dirx[2] = ( m_ImageHeader->trhcS - m_ImageHeader->tlhcS );
518 dirx.normalize();
519
520 diry[0] = -( m_ImageHeader->brhcR - m_ImageHeader->trhcR );
521 diry[1] = -( m_ImageHeader->brhcA - m_ImageHeader->trhcA );
522 diry[2] = ( m_ImageHeader->brhcS - m_ImageHeader->trhcS );
523 diry.normalize();
524
525 dirz[0] = -m_ImageHeader->normR;
526 dirz[1] = -m_ImageHeader->normA;
527 dirz[2] = m_ImageHeader->normS;
528 dirz.normalize();
529
530 // Set the directions
531 this->SetDirection(0, dirx);
532 this->SetDirection(1, diry);
533 this->SetDirection(2, dirz);
534
535 // See if slices need to be reversed. itk uses a right hand
536 // coordinate system. If the computed slice direction is opposite
537 // the direction in the header, the files have to be read in reverse
538 // order.
539 vnl_vector< double > sliceDirection = vnl_cross_3d(dirx, diry);
540 if ( dot_product(sliceDirection, dirz) < 0 )
541 {
542 // Use the computed direction
543 this->SetDirection(2, sliceDirection);
544
545 // Sort image list in reverse order
546 m_FilenameList->SetSortOrder(IPLFileNameList::SortGlobalDescend);
547 m_FilenameList->sortImageList();
548 }
549
550 // Compute the spacing between two slices from the origins of the
551 // first two files in the study
552 if ( m_FilenameList->NumFiles() > 1 )
553 {
554 auto it = m_FilenameList->begin();
555
556 // The first file
557 std::string file1 = ( *it )->GetImageFileName();
558
559 // The second file
560 it++;
561 std::string file2 = ( *it )->GetImageFileName();
562
563 GEImageHeader *hdr1 = this->ReadHeader( file1.c_str() );
564 GEImageHeader *hdr2 = this->ReadHeader( file2.c_str() );
565
566 float origin1[3], origin2[3];
567 origin1[0] = hdr1->tlhcR;
568 origin1[1] = hdr1->tlhcA;
569 origin1[2] = hdr1->tlhcS;
570
571 // Origin shopuld always come from the first slice
572 this->SetOrigin(0, -hdr1->tlhcR);
573 this->SetOrigin(1, -hdr1->tlhcA);
574 this->SetOrigin(2, hdr1->tlhcS);
575
576 origin2[0] = hdr2->tlhcR;
577 origin2[1] = hdr2->tlhcA;
578 origin2[2] = hdr2->tlhcS;
579
580 float distanceBetweenTwoSlices = std::sqrt(
581 ( origin1[0] - origin2[0] ) * ( origin1[0] - origin2[0] )
582 + ( origin1[1] - origin2[1] ) * ( origin1[1] - origin2[1] )
583 + ( origin1[2] - origin2[2] ) * ( origin1[2] - origin2[2] ) );
584
585 this->SetSpacing(2, distanceBetweenTwoSlices);
586
587 // Cleanup
588 delete hdr1;
589 delete hdr2;
590 }
591 else
592 // If there is only one slice, the use it's origin
593 {
594 this->SetOrigin(0, -m_ImageHeader->tlhcR);
595 this->SetOrigin(1, -m_ImageHeader->tlhcA);
596 this->SetOrigin(2, m_ImageHeader->tlhcS);
597 }
598 }
599 } // end namespace itk
600