[Insight-developers] enhancement to itk::GDCMImageIO

Mathieu Malaterre mathieu.malaterre at kitware.com
Thu Feb 10 11:26:59 EST 2005


Kent,

	Good work, that looks great. But I was thinking shouldn't this be 
ultimately in GDCM itself. I am already concern about DICOM image which 
would be missing Image Orientation Patient, but would provide Patient 
Orientation (20,20).

	Or some images (old ACR/NEMA) also defined Image Orientation (RET) in 
(20,35) instead, and it has the same meaning as the one in DICOM v3.

	I can alreay patch gdcm 0.6 in ITK to do this, or GDCM CVS has a notion 
of Helper that provide /augmented/ methods like GetSpacing where in the 
case of the string is "0\0" (true!), it's able to default to (1,1).

	In both case I won't be commiting any changes until CVS is open again.

What do you think ?
Mathieu

Kent Williams wrote:
> 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 ;-)
> 
> 
> ------------------------------------------------------------------------
> 
> *** 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()
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers



More information about the Insight-developers mailing list