[Insight-developers] enhancement to itk::GDCMImageIO

[Kent Williams]

Enclosed is a context diff for a proposed change to itkGDCMImageIO.cxx, 
which sets the ITK_CoordinateOrientation attribute in the 
MetadataDictionary for images read in via GDCM.

I didn't just check it in, partly because it's relatively late in the 
day, and partly because I would like feedback about image orientation in 
ImageIO in general, and the way I've added it to the GDCMImageIO class 
in particular.

I don't think what I've done would break anything but I've been wrong 
before ;-)
-------------- next part --------------
*** itkGDCMImageIO.cxx.~1.54.~	2005-01-25 09:45:48.000000000 -0600
--- itkGDCMImageIO.cxx	2005-02-09 13:39:37.000000000 -0600
***************
*** 18,24 ****
  
  =========================================================================*/
  #include "itkGDCMImageIO.h"
! 
  #include "itkMetaDataObject.h"
  #include <itksys/Base64.h>
  #include "gdcm/src/gdcmValEntry.h" //internal of gdcm
--- 18,24 ----
  
  =========================================================================*/
  #include "itkGDCMImageIO.h"
! #include "itkIOCommon.h"
  #include "itkMetaDataObject.h"
  #include <itksys/Base64.h>
  #include "gdcm/src/gdcmValEntry.h" //internal of gdcm
***************
*** 258,263 ****
--- 258,413 ----
    file.close();
  }
  
+ namespace gdcmOrient
+ {
+ inline float Abs(float x)
+ {
+   if(x < 0)
+     return -x;
+   return x;
+ }
+ 
+ inline unsigned Max3(float x, float y, float z)
+ {
+   if(Abs(x) > Abs(y) && Abs(x) > Abs(z))
+     {
+     return 0;
+     }
+   if(Abs(y) > Abs(x) && Abs(y) > Abs(z))
+     {
+     return 1;
+     }
+   if(Abs(z) > Abs(x) && Abs(z) > Abs(y))
+     {
+     return 2;
+     }
+   // they must all be equal, so just say x
+   return 0;
+ }
+ 
+ inline int Sign(float x)
+ {
+   if(x < 0)
+     return -1;
+   return 1;
+ }
+ 
+ static void getcosines(const std::string &orientstring,float *cosines)
+ {
+   int i = 0;
+   std::string::size_type numstart, numend;
+   numstart = 0;
+   while(numstart != std::string::npos)
+     {
+     numend = orientstring.find("\\",numstart);
+     while(numstart != numend && orientstring[numstart] == ' ')
+       {
+       numstart++;
+       }
+     cosines[i++] = atof(orientstring.substr(numstart,numend-numstart).c_str());
+     numstart = (numend == std::string::npos ? numend : numend + 1);
+     }
+ }
+ //
+ // this is all based on this from the DICOM standard
+ //
+ // C.7.6.2.1.1
+ // Image Orientation (0020,0037) specifies the direction 
+ // cosines of the first row and the first column with 
+ // respect to the patient. These Attributes shall be 
+ // provide as a pair. Row value for the x, y, and z 
+ // axes respectively followed by the Column value for 
+ // the x, y, and z axes respectively. The direction 
+ // of the axes is defined fully by the patient's 
+ // orientation. The x-axis is increasing to the 
+ // left hand side of the patient. The y-axis is increasing 
+ // to the posterior side of the patient. The z-axis is 
+ // increasing toward the head of the patient. The patient 
+ // based coordinate system is a right handed system, 
+ // i.e. the vector cross product of a unit vector 
+ // along the positive x-axis and a unit vector 
+ // along the positive y-axis is equal to a unit 
+ // vector along the positive z-axis.
+ 
+ static SpatialOrientation::ValidCoordinateOrientationFlags 
+ dicomOrient(const std::string &orientstring)
+ {
+   float cosines[6];
+   int axes[9] = {0,0,0,0,0,0,0,0,0};
+   int dominant_axis;
+ 
+   // convert orientation string to floats
+   getcosines(orientstring,cosines);
+ 
+   // figure out the dominant axis of the row
+   dominant_axis = Max3(cosines[0],
+                        cosines[1],
+                        cosines[2]);
+   axes[dominant_axis] = 
+     Sign(cosines[dominant_axis]);
+ 
+   // figure the dominant axis of the column dimension
+   dominant_axis = Max3(cosines[3],
+                        cosines[4],
+                        cosines[5]);
+   axes[dominant_axis + 3] =
+     Sign(cosines[dominant_axis + 3]);
+ 
+   // cross product to get slice direction
+   axes[6] = axes[1] * axes[5] - axes[2] * axes[4];
+   axes[7] = axes[2] * axes[3] - axes[0] * axes[5];
+   axes[8] = axes[0] * axes[4] - axes[1] * axes[3];
+ 
+   //
+   // swizzle up the spatial orientation.
+   // this is based on the idea that SpatialOrientation in the patient's
+   // frame of reference.
+ #define SICTZero static_cast<SpatialOrientation::CoordinateTerms>(0)
+   SpatialOrientation::CoordinateTerms terms[3] = {SICTZero,SICTZero,SICTZero};
+   for(unsigned i = 0; i < 3; i++)
+     {
+     if(axes[(i*3)] == 1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Right;
+       }
+     else if(axes[(i*3)] == -1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Left;
+       }
+     else if(axes[(i*3)+1] == 1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Anterior;
+       }
+     else if(axes[(i*3)+1] == -1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Posterior;
+       }
+     else if(axes[(i*3)+2] == 1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Inferior;
+       }
+     else if(axes[(i*3)+2] == -1)
+       {
+       terms[i] = SpatialOrientation::ITK_COORDINATE_Superior;
+       }
+     }
+   //
+   // all terms must be defined, otherwise just punt
+   if(terms[0] == SICTZero || terms[1] == SICTZero || terms[2] == SICTZero)
+     {
+     return SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP;
+     }
+   return static_cast<SpatialOrientation::ValidCoordinateOrientationFlags>
+     ((terms[0] << 
+      SpatialOrientation::ITK_COORDINATE_PrimaryMinor) +
+     (terms[1] << 
+      SpatialOrientation::ITK_COORDINATE_SecondaryMinor) +
+     (terms[2] << 
+      SpatialOrientation::ITK_COORDINATE_TertiaryMinor));
+     
+ }
+ }
+ 
  void GDCMImageIO::InternalReadImageInformation(std::ifstream& file)
  {
    //read header
***************
*** 395,400 ****
--- 545,564 ----
  
      d = header.GetNextEntry();
    }
+   //
+   // encode the image orientation, based on Image Orientation (Patient)
+   SpatialOrientation::ValidCoordinateOrientationFlags coord_orient;
+   std::string iop = header.GetEntryByNumber(0x0020,0x0037);
+   if(iop == "gdcm::Unfound")
+     {
+     // total punt, but any value has a 1 in 48 chance of being right.
+     coord_orient = SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP;
+     }
+   else
+     {
+     coord_orient = gdcmOrient::dicomOrient(iop);
+     }
+     itk::EncapsulateMetaData<SpatialOrientation::ValidCoordinateOrientationFlags>(dico,ITK_CoordinateOrientation, coord_orient);
  }
  
  void GDCMImageIO::ReadImageInformation()