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