gdcmOrientation.cxx

Go to the documentation of this file.
00001 /*=========================================================================
00002                                                                                 
00003   Program:   gdcm
00004   Module:    $RCSfile: gdcmOrientation.cxx,v $
00005   Language:  C++
00006   Date:      $Date: 2007/05/23 14:18:10 $
00007   Version:   $Revision: 1.25 $
00008                                                                                 
00009   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
00010   l'Image). All rights reserved. See Doc/License.txt or
00011   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
00012                                                                                 
00013      This software is distributed WITHOUT ANY WARRANTY; without even
00014      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00015      PURPOSE.  See the above copyright notices for more information.
00016                                                                                 
00017 =========================================================================*/
00018 
00019 #include "gdcmOrientation.h"
00020 #include "gdcmFile.h"
00021 #include "gdcmDebug.h"
00022 #include <math.h> // for sqrt
00023 
00024 namespace GDCM_NAME_SPACE 
00025 {
00026 //--------------------------------------------------------------------
00027 //  THERALYS Algorithm to determine the most similar basic orientation
00028 //
00029 //  Transliterated from original Python code.
00030 //  Kept as close as possible to the original code
00031 //  in order to speed up any further modif of Python code :-(
00032 //-----------------------------------------------------------------------
00033 
00055 static const char  *OrientationTypeStrings[] = { 
00056   "Not Applicable",
00057   "Axial",
00058   "Coronal",
00059   "Sagital",
00060   "Heart Axial",
00061   "Heart Coronal",
00062   "Heart Sagital",
00063   "Axial invert",
00064   "Coronal invert",
00065   "Sagital invert",
00066   "Heart Axial invert",
00067   "Heart Coronal invert",
00068   "Heart Sagital invert",
00069   NULL
00070 };
00071 
00074 const char* Orientation::GetOrientationTypeString(OrientationType const o)
00075 {
00076   int k = (int)o;
00077   if (k < 0) 
00078        k = -k + 6;
00079 
00080   return OrientationTypeStrings[k];
00081 }
00082 
00085 OrientationType Orientation::GetOrientationType( File *f )
00086 {
00087    float iop[6];
00088    bool succ = f->GetImageOrientationPatient( iop );
00089    if ( !succ )
00090    {
00091       gdcmWarningMacro( "No Image Orientation (0020,0037)/(0020,0032) found in the file, cannot proceed." )
00092       return NotApplicable;
00093    }
00094    vector3D ori1;
00095    vector3D ori2;
00096 
00097    ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2]; 
00098    ori2.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5];
00099 
00100    // two perpendicular vectors describe one plane
00101    double dicPlane[6][2][3] =
00102    { {  { 1,   0,    0   },{ 0,      1,     0     }  }, // Axial
00103      {  { 1,   0,    0   },{ 0,      0,    -1     }  }, // Coronal
00104      {  { 0,   1,    0   },{ 0,      0,    -1     }  }, // Sagittal
00105      {  { 0.8, 0.5,  0.0 },{-0.1,    0.1 , -0.95  }  }, // Axial - HEART
00106      {  { 0.8, 0.5,  0.0 },{-0.6674, 0.687, 0.1794}  }, // Coronal - HEART
00107      {  {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794}  }  // Sagittal - HEART
00108    };
00109 
00110    vector3D refA;
00111    vector3D refB;
00112    int i = 0;
00113    Res res;   // [ <result> , <memory of the last succes calcule> ]
00114    res.first = 0;
00115    res.second = 99999;
00116 
00117    for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
00118    {
00119        ++i;
00120        // refA=plane[0]
00121        refA.x = dicPlane[numDicPlane][0][0]; 
00122        refA.y = dicPlane[numDicPlane][0][1]; 
00123        refA.z = dicPlane[numDicPlane][0][2];
00124        // refB=plane[1]
00125        refB.x = dicPlane[numDicPlane][1][0]; 
00126        refB.y = dicPlane[numDicPlane][1][1]; 
00127        refB.z = dicPlane[numDicPlane][1][2];
00128        res=VerfCriterion(  i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
00129        res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
00130    }
00131    // res thought looks like is a float value, but is indeed an int
00132    // try casting it to int first then enum value to please VS7:
00133    int int_res = (int)res.first;
00134    gdcmAssertMacro( int_res <= 6 && int_res >= -6);
00135    return (OrientationType)int_res;
00136 }
00137 
00138 Res 
00139 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
00140 {
00141    Res res;
00142    double type = in.first;
00143    double criterion = in.second;
00144    if (/*criterionNew < 0.1 && */criterionNew < criterion)
00145    {
00146       type      = typeCriterion;
00147       criterion = criterionNew;
00148    }
00149    res.first  = type;
00150    res.second = criterion;
00151    return res;
00152 } 
00153 
00154 inline double square_dist(vector3D const &v1, vector3D const &v2)
00155 {
00156   double res;
00157   res = (v1.x - v2.x)*(v1.x - v2.x) +
00158         (v1.y - v2.y)*(v1.y - v2.y) +
00159         (v1.z - v2.z)*(v1.z - v2.z);
00160   return res;
00161 }
00162 
00163 //------------------------- Purpose : -----------------------------------
00164 //- This function determines the orientation similarity of two planes.
00165 //  Each plane is described by two vectors.
00166 //------------------------- Parameters : --------------------------------
00167 //- <refA>  : - type : vector 3D (double)
00168 //- <refB>  : - type : vector 3D (double)
00169 //            - Description of the first plane
00170 //- <ori1>  : - type : vector 3D (double)
00171 //- <ori2>  : - type : vector 3D (double)
00172 //            - Description of the second plane
00173 //------------------------- Return : ------------------------------------
00174 // double :   0 if the planes are perpendicular. While the difference of
00175 //            the orientation between the planes are big more enlarge is
00176 //            the criterion.
00177 //------------------------- Other : -------------------------------------
00178 // The calculus is based with vectors normalice
00179 double
00180 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB, 
00181                                   vector3D const &ori1, vector3D const &ori2 )
00182 {
00183 
00184    vector3D ori3 = ProductVectorial(ori1,ori2);
00185    vector3D refC = ProductVectorial(refA,refB);
00186    double res = square_dist(refC, ori3);
00187 
00188    return sqrt(res);
00189 }
00190 
00191 //------------------------- Purpose : -----------------------------------
00192 //- Calculus of the poduct vectorial between two vectors 3D
00193 //------------------------- Parameters : --------------------------------
00194 //- <vec1>  : - type : vector 3D (double)
00195 //- <vec2>  : - type : vector 3D (double)
00196 //------------------------- Return : ------------------------------------
00197 // (vec) :    - Vector 3D
00198 //------------------------- Other : -------------------------------------
00199 vector3D
00200 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2)
00201 {
00202    vector3D vec3;
00203    vec3.x =    vec1.y*vec2.z - vec1.z*vec2.y;
00204    vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x);
00205    vec3.z =    vec1.x*vec2.y - vec1.y*vec2.x;
00206 
00207    return vec3;
00208 }
00209 
00210 
00211 // ------------------------------------------------------------------------
00212 /*
00213 2.2.2 Orientation of DICOM images
00214 
00215 
00216 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation
00217 says :
00218 
00219 A question that is frequently asked in comp.protocols.dicom is how to determine
00220  which side of an image is which (e.g. left, right) and so on. 
00221  The short answer is that for projection radiographs this is specified
00222  explicitly using the "Patient Orientation" attribute, and for cross-sectional
00223  images it needs to be derived from the "Image Orientation (Patient)" direction
00224  cosines. In the standard these are explained as follows:
00225 
00226     * "C.7.6.1.1.1 Patient Orientation. 
00227                 The Patient Orientation (0020,0020) relative to the image
00228                 plane shall be specified by two values that designate the 
00229                 anatomical direction of the positive row axis (left to right)
00230                 and the positive column axis (top to bottom). 
00231                 The first entry is the direction of the rows, given by the 
00232                 direction of the last pixel in the first row from the first 
00233                 pixel in that row. 
00234                 The second entry is the direction of the columns, given by 
00235                 the direction of the last pixel in the first column from the
00236                 first pixel in that column. 
00237                 Anatomical direction shall be designated by the capital 
00238                 letters: A (anterior), P (posterior), R (right),L (left), 
00239                 H (head), F (foot). 
00240                 Each value of the orientation attribute shall contain at 
00241                 least one of these characters. 
00242                 If refinements in the orientation descriptions are to be 
00243                 specified, then they shall be designated by one or two 
00244                 additional letters in each value. 
00245                 Within each value, the letters shall be ordered with the 
00246                 principal orientation designated in the first character."
00247  
00248     * "C.7.6.2.1.1 Image Position And Image Orientation. 
00249                 The "Image Position (Patient)" (0020,0032) specifies the x, y, 
00250                 and z coordinates of the upper left hand corner of the image; 
00251                 it is the center of the first voxel transmitted. 
00252                 The "Image Orientation (Patient)" (0020,0037) specifies the 
00253                 direction cosines of the first row and the first column 
00254                 with respect to the patient. 
00255                 These Attributes shall be provided as a pair. 
00256                 Row value for the x, y, and z axes respectively followed by 
00257                 the Column value for the x, y, and z axes respectively. 
00258                 The direction of the axes is defined fully by the patient's 
00259                 orientation. 
00260                 The x-axis is increasing to the left hand side of the patient.
00261                 The y-axis is increasing to the posterior side of the patient.
00262                 The z-axis is increasing toward the head of the patient. 
00263                 The patient based coordinate system is a right handed system,
00264                 i.e. the vector cross product of a unit vector along the 
00265                 positive x-axis and a unit vector along the positive y-axis
00266                 is equal to a unit vector along the positive z-axis." 
00267 
00268 Some simple code to take one of the direction cosines (vectors) from the 
00269 Image Orientation (Patient) attribute and generate strings equivalent to one 
00270 of the values of Patient Orientation looks like this (noting that if the vector
00271 is not aligned exactly with one of the major axes, the resulting string will 
00272 have multiple letters in as described under "refinements" in C.7.6.1.1.1): 
00273 
00274 */
00275 
00294 std::string Orientation::GetOrientation ( File *f )
00295 {
00296    float iop[6];
00297    if ( !f->GetImageOrientationPatient( iop ) )
00298    return GDCM_UNFOUND;
00299 
00300    std::string orientation;
00301    orientation = GetSingleOrientation ( iop ) 
00302                + "\\" 
00303                + GetSingleOrientation ( iop + 3 );
00304    return orientation;
00305 }
00306 
00307 
00308 std::string Orientation::GetSingleOrientation ( float *iop)
00309 {
00310    std::string orientation;
00311 
00312    char orientationX = iop[0] < 0 ? 'R' : 'L';
00313    char orientationY = iop[1] < 0 ? 'A' : 'P';
00314    char orientationZ = iop[2] < 0 ? 'F' : 'H';
00315 
00316    double absX = iop[0];
00317    if (absX < 0) absX = -absX;
00318       double absY = iop[1];
00319    if (absY < 0) absY = -absY;
00320       double absZ = iop[2];
00321    if (absZ < 0) absZ = -absZ;
00322 
00323    for (int i=0; i<3; ++i) 
00324    {
00325       if (absX>.0001 && absX>absY && absX>absZ) 
00326       {
00327          orientation = orientation + orientationX;
00328          absX=0;
00329        }
00330        else if (absY>.0001 && absY>absX && absY>absZ) 
00331        {
00332           orientation = orientation + orientationY;
00333           absY=0;
00334        }
00335        else if (absZ>.0001 && absZ>absX && absZ>absY) 
00336        {
00337            orientation = orientation + orientationZ;
00338            absZ=0;
00339        }
00340        else 
00341           break;
00342      }
00343    return orientation;
00344 } 
00345 
00346 
00347 /*-------------------------------------------------------------------
00348 
00349 Some more stuff, from XMedcon
00350 
00351 ---> Personal remark from JPRx :
00352 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP " 
00353 --> or, not so common, 
00354 // HFDR = Head First-Decubitus Right
00355 // HFDL = Head First-Decubitus Left
00356 // FFDR = Feet First-Decubitus Right
00357 // FFDL = Feet First-Decubitus Left
00358 --> the cosines may have any value -1.< <+1., for MR images !
00359 
00360 enum patient_slice_orientation_type
00361   {
00362     patient_slice_orientation_unknown = 0,
00363     supine_headfirst_transaxial,
00364     supine_headfirst_sagittal,
00365     supine_headfirst_coronal,
00366     supine_feetfirst_transaxial,
00367     supine_feetfirst_sagittal,
00368     supine_feetfirst_coronal,
00369     prone_headfirst_transaxial,
00370     prone_headfirst_sagittal,
00371     prone_headfirst_coronal,
00372     prone_feetfirst_transaxial,
00373     prone_feetfirst_sagittal,
00374     prone_feetfirst_coronal
00375   };
00376 
00377 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
00378 {
00379   h.GetImageOrientationPatient(image_orientation_patient);
00380 }
00381 
00382 #if 0
00383 //
00384 // this is all completely cribbed from the xmedcon library, since
00385 // we're trying to do what it does, mostly.
00386 patient_slice_orientation_type
00387 GetPatSliceOrient(gdcm::File &h)
00388 {
00389   F32 image_orientation_patient[6];
00390 
00391   // protected, do it the hard way
00392   //  h.GetImageOrientationPatient(image_orientation_patient);
00393   GetImageOrientationPatient(h,image_orientation_patient);
00394 
00395   enum { headfirst, feetfirst } patient_orientation;
00396   enum { supine, prone } patient_rotation;
00397   enum { transaxial, sagittal, coronal } slice_orientation;
00398 
00399   std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
00400   if(patient_position == "gdcm::Unfound")
00401     {
00402     patient_position = "HF";
00403     }
00404   if(patient_position.find("HF") != std::string::npos)
00405     {
00406     patient_orientation = headfirst;
00407     }
00408   else if(patient_position.find("FF") != std::string::npos)
00409     {
00410     patient_orientation = feetfirst;
00411     }
00412 
00413   if(patient_position.find("S") != std::string::npos)
00414     {
00415     patient_rotation = supine;
00416     }
00417   else if(patient_position.find("P") != std::string::npos)
00418     {
00419     patient_rotation = prone;
00420     }
00421 
00422   if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
00423      (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
00424     {
00425     slice_orientation = transaxial;
00426     }
00427   else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
00428           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
00429     {
00430     slice_orientation = sagittal;
00431     }
00432   else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
00433           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
00434     {
00435     slice_orientation = coronal;
00436     }
00437   //
00438   // combine results
00439   patient_slice_orientation_type patient_slice_orientation = 
00440     patient_slice_orientation_unknown;
00441   switch (patient_rotation) 
00442     {
00443     case supine:
00444       switch (patient_orientation) 
00445         {
00446         case headfirst:
00447           switch (slice_orientation) 
00448             {
00449             case transaxial:
00450               patient_slice_orientation = supine_headfirst_transaxial;
00451               break;
00452             case sagittal:
00453               patient_slice_orientation = supine_headfirst_sagittal;
00454               break;
00455             case coronal:
00456               patient_slice_orientation = supine_headfirst_coronal;
00457               break;
00458             }
00459           break;
00460         case feetfirst:
00461           switch (slice_orientation) 
00462             {
00463             case transaxial:
00464               patient_slice_orientation = supine_feetfirst_transaxial;
00465               break;
00466             case sagittal:
00467               patient_slice_orientation = supine_feetfirst_sagittal;
00468               break;
00469             case coronal:
00470               patient_slice_orientation = supine_feetfirst_coronal;
00471               break;
00472             }
00473           break;
00474         }
00475       break;
00476     case prone:
00477       switch (patient_orientation) 
00478         {
00479         case headfirst:
00480           switch (slice_orientation) 
00481             {
00482             case transaxial:
00483               patient_slice_orientation = prone_headfirst_transaxial;
00484               break;
00485             case sagittal:
00486               patient_slice_orientation = prone_headfirst_sagittal;
00487               break;
00488             case coronal:
00489               patient_slice_orientation = prone_headfirst_coronal;
00490               break;
00491             }
00492           break;
00493         case feetfirst:
00494           switch (slice_orientation) 
00495             {
00496             case transaxial:
00497               patient_slice_orientation = prone_feetfirst_transaxial;
00498               break;
00499             case sagittal:
00500               patient_slice_orientation = prone_feetfirst_sagittal;
00501               break;
00502             case coronal:
00503               patient_slice_orientation = prone_feetfirst_coronal;
00504               break;
00505             }
00506           break;
00507         }
00508       break;
00509     }
00510   if(patient_slice_orientation != patient_slice_orientation_unknown)
00511     return patient_slice_orientation;
00512   //
00513   // this is what xmedcon does
00514   if ((image_orientation_patient[0] == +1)   &&
00515       (image_orientation_patient[4] == +1))   
00516     patient_slice_orientation = supine_headfirst_transaxial;
00517   else if ((image_orientation_patient[0] == -1)   &&
00518            (image_orientation_patient[4] == +1))   
00519     patient_slice_orientation = supine_feetfirst_transaxial;
00520   else if ((image_orientation_patient[0] == -1)   &&
00521            (image_orientation_patient[4] == -1))   
00522     patient_slice_orientation = prone_headfirst_transaxial;
00523   else if ((image_orientation_patient[0] == +1)   &&
00524            (image_orientation_patient[4] == -1))   
00525     patient_slice_orientation = prone_feetfirst_transaxial;
00526 
00527   else if ((image_orientation_patient[1] == +1)   &&
00528            (image_orientation_patient[5] == -1))   
00529     patient_slice_orientation = supine_headfirst_sagittal;
00530   else if ((image_orientation_patient[1] == +1)   &&
00531            (image_orientation_patient[5] == +1))   
00532     patient_slice_orientation = supine_feetfirst_sagittal;
00533   else if ((image_orientation_patient[1] == -1)   &&
00534            (image_orientation_patient[5] == -1))   
00535     patient_slice_orientation = prone_headfirst_sagittal;
00536   else if ((image_orientation_patient[1] == -1)   &&
00537            (image_orientation_patient[5] == +1))   
00538     patient_slice_orientation = prone_feetfirst_sagittal;
00539 
00540   else if ((image_orientation_patient[0] == +1)   &&
00541            (image_orientation_patient[5] == -1))   
00542     patient_slice_orientation = supine_headfirst_coronal;
00543   else if ((image_orientation_patient[0] == -1)   &&
00544            (image_orientation_patient[5] == +1))   
00545     patient_slice_orientation = supine_feetfirst_coronal;
00546   else if ((image_orientation_patient[0] == -1)   &&
00547            (image_orientation_patient[5] == -1))   
00548     patient_slice_orientation = prone_headfirst_coronal;
00549   else if ((image_orientation_patient[0] == +1)   &&
00550            (image_orientation_patient[5] == +1))   
00551     patient_slice_orientation = prone_feetfirst_coronal;
00552   return patient_slice_orientation;
00553 }
00554 #else
00555 
00556 -------------------------------------------------------------------------*/
00557 
00558 } // end namespace gdcm

Generated on Fri Aug 24 12:53:17 2007 for gdcm by  doxygen 1.4.6