src/gdcmOrientation.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmOrientation.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/09/29 14:23:58 $
   7   Version:   $Revision: 1.13 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 #include "gdcmOrientation.h"
  20 #include "gdcmFile.h"
  21 #include "gdcmDebug.h"
  22 #include <math.h> // for sqrt
  23
  24 namespace gdcm
  25 {
  26 //--------------------------------------------------------------------
  27 //  THERALYS Algorithm to determine the most similar basic orientation
  28 //
  29 //  Transliterated from original Python code.
  30 //  Kept as close as possible to the original code
  31 //  in order to speed up any further modif of Python code :-(
  32 //-----------------------------------------------------------------------
  33
  34 /**
  35  * \brief  THERALYS' Algorithm to determine the most similar basic orientation
  36  *           (Axial, Coronal, Sagital) of the image
  37  * \note Should be run on the first gdcm::File of a 'coherent' Serie
  38  * @return orientation code
  39  *   #   0 : Not Applicable (neither 0020,0037 Image Orientation Patient
  40  *   #                       nor     0020,0032 Image Position           found)
  41  *   #   1 : Axial
  42  *   #  -1 : Axial invert
  43  *   #   2 : Coronal
  44  *   #  -2 : Coronal invert
  45  *   #   3 : Sagital
  46  *   #  -3 : Sagital invert
  47  *   #   4 : Heart Axial
  48  *   #  -4 : Heart Axial invert
  49  *   #   5 : Heart Coronal
  50  *   #  -5 : Heart Coronal invert
  51  *   #   6 : Heart Sagital
  52  *   #  -6 : Heart Sagital invert
  53  */
  54 double Orientation::TypeOrientation( File *f )
  55 {
  56    float iop[6];
  57    bool succ = f->GetImageOrientationPatient( iop );
  58    if ( !succ )
  59    {
  60       gdcmErrorMacro( "No Image Orientation (0020,0037) found in the file, cannot proceed." )
  61       return 0;
  62    }
  63 /*
  64 std::cout << " iop : ";
  65 for(int i=0;i<6;i++)
  66    std::cout << iop[i] << "  ";
  67 std::cout << std::endl;
  68 */
  69    vector3D ori1;
  70    vector3D ori2;
  71
  72    ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2];
  73    ori2.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5];
  74
  75    // two perpendicular vectors describe one plane
  76    double dicPlane[6][2][3] =
  77    { {  { 1,   0,    0   },{ 0,      1,     0     }  }, // Axial
  78      {  { 1,   0,    0   },{ 0,      0,    -1     }  }, // Coronal
  79      {  { 0,   1,    0   },{ 0,      0,    -1     }  }, // Sagittal
  80      {  { 0.8, 0.5,  0.0 },{-0.1,    0.1 , -0.95  }  }, // Axial - HEART
  81      {  { 0.8, 0.5,  0.0 },{-0.6674, 0.687, 0.1794}  }, // Coronal - HEART
  82      {  {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794}  }  // Sagittal - HEART
  83    };
  84
  85    vector3D refA;
  86    vector3D refB;
  87    int i = 0;
  88    Res res;   // [ <result> , <memory of the last succes calcule> ]
  89    res.first = 0;
  90    res.second = 99999;
  91
  92    for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
  93    {
  94        ++i;
  95        // refA=plane[0]
  96        refA.x = dicPlane[numDicPlane][0][0];
  97        refA.y = dicPlane[numDicPlane][0][1];
  98        refA.z = dicPlane[numDicPlane][0][2];
  99        // refB=plane[1]
 100        refB.x = dicPlane[numDicPlane][1][0];
 101        refB.y = dicPlane[numDicPlane][1][1];
 102        refB.z = dicPlane[numDicPlane][1][2];
 103        res=VerfCriterion(  i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
 104        res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
 105    }
 106    return res.first;
 107 /*
 108 //  i=0
 109 //  res=[0,99999]  ## [ <result> , <memory of the last succes calculus> ]
 110 //  for plane in dicPlane:
 111 //      i=i+1
 112 //      refA=plane[0]
 113 //      refB=plane[1]
 114 //      res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
 115 //      res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
 116 //  return res[0]
 117 */
 118 }
 119
 120 Res
 121 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
 122 {
 123    Res res;
 124    double type = in.first;
 125    double criterion = in.second;
 126    if (/*criterionNew < 0.1 && */criterionNew < criterion)
 127    {
 128       type      = typeCriterion;
 129       criterion = criterionNew;
 130    }
 131    res.first  = type;
 132    res.second = criterion;
 133    return res;
 134 /*
 135 //   type = res[0]
 136 //   criterion = res[1]
 137 // #     if criterionNew<0.1 and criterionNew<criterion:
 138 //   if criterionNew<criterion:
 139 //      criterion=criterionNew
 140 //      type=typeCriterion
 141 //   return [ type , criterion ]
 142 */
 143
 144 }
 145
 146 inline double square_dist(vector3D const &v1, vector3D const &v2)
 147 {
 148   double res;
 149   res = (v1.x - v2.x)*(v1.x - v2.x) +
 150         (v1.y - v2.y)*(v1.y - v2.y) +
 151         (v1.z - v2.z)*(v1.z - v2.z);
 152   return res;
 153 }
 154
 155 //------------------------- Purpose : -----------------------------------
 156 //- This function determines the orientation similarity of two planes.
 157 //  Each plane is described by two vectors.
 158 //------------------------- Parameters : --------------------------------
 159 //- <refA>  : - type : vector 3D (double)
 160 //- <refB>  : - type : vector 3D (double)
 161 //            - Description of the first plane
 162 //- <ori1>  : - type : vector 3D (double)
 163 //- <ori2>  : - type : vector 3D (double)
 164 //            - Description of the second plane
 165 //------------------------- Return : ------------------------------------
 166 // double :   0 if the planes are perpendicular. While the difference of
 167 //            the orientation between the planes are big more enlarge is
 168 //            the criterion.
 169 //------------------------- Other : -------------------------------------
 170 // The calculus is based with vectors normalice
 171 double
 172 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB,
 173                                   vector3D const &ori1, vector3D const &ori2 )
 174 {
 175
 176    vector3D ori3 = ProductVectorial(ori1,ori2);
 177    vector3D refC = ProductVectorial(refA,refB);
 178    double res = square_dist(refC, ori3);
 179
 180    return sqrt(res);
 181 }
 182
 183 //------------------------- Purpose : -----------------------------------
 184 //- Calculus of the poduct vectorial between two vectors 3D
 185 //------------------------- Parameters : --------------------------------
 186 //- <vec1>  : - type : vector 3D (double)
 187 //- <vec2>  : - type : vector 3D (double)
 188 //------------------------- Return : ------------------------------------
 189 // (vec) :    - Vector 3D
 190 //------------------------- Other : -------------------------------------
 191 vector3D
 192 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2)
 193 {
 194    vector3D vec3;
 195    vec3.x =    vec1.y*vec2.z - vec1.z*vec2.y;
 196    vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x);
 197    vec3.z =    vec1.x*vec2.y - vec1.y*vec2.x;
 198
 199    return vec3;
 200 }
 201
 202
 203
 204 // ---------------------------------------------------------------------------
 205 // Here is the original Python code, kindly supplied by THERALYS
 206 //
 207 // C++ code doesn't give good results
 208 // --> FIXME
 209
 210 /*
 211
 212 def TypeOrientation(self,file0):
 213 """
 214 # ------------------------- Purpose : -----------------------------------
 215 # - This function compare the orientation of the given image and the
 216 #   basics orientations (Axial, Cornal, Sagital)
 217 # ------------------------- Parameters : --------------------------------
 218 # - <file0> : - type : string
 219 #             - The name of the first image file of the serie
 220 # ------------------------- Return : ------------------------------------
 221 #   1 :   Axial
 222 #  -1 :   Axial invert
 223 #   2 :   Coronal
 224 #  -2 :   Coronal invert
 225 #   3 :   Sagital
 226 #  -3 :   Sagital invert
 227 #   4 :   Heart Axial
 228 #  -4 :   Heart Axial invert
 229 #   5 :   Heart Coronal
 230 #  -5 :   Heart Coronal invert
 231 #   6 :   Heart Sagital
 232 #  -6 :   Heart Sagital invert
 233 #
 234    # ------------------------- Other : -------------------------------------
 235 # This method finds the most similar basic orientation.
 236 """
 237 try:
 238    toRead = gdcm.File(file0)
 239    ValDict = GetValuesDict(toRead)
 240    try:
 241       imageOrientation=ValDict["Image Orientation (Patient)"]
 242    except KeyError:
 243       imageOrientation=ValDict["Image Orientation"]
 244
 245    ori1=[float(split(imageOrientation,"\\")[0]),\
 246       float(split(imageOrientation,"\\")[1]),\
 247       float(split(imageOrientation,"\\")[2])]
 248    ori2=[float(split(imageOrientation,"\\")[3]),\
 249       float(split(imageOrientation,"\\")[4]),\
 250       float(split(imageOrientation,"\\")[5])]
 251
 252 ## two vectors perpendicular describe one plane
 253    dicPlane=[ [  [1,0,0],[0,1,0]   ],  ## Axial
 254             [  [1,0,0],[0,0,-1]  ],  ## Coronal
 255             [  [0,1,0],[0,0,-1]  ],  ## Sagittal
 256             [  [ 0.8 , 0.5 ,  0.0 ],[-0.1 , 0.1 , -0.95]        ],## Axial - HEART
 257             [  [ 0.8 , 0.5 ,  0.0 ],[-0.6674 , 0.687 , 0.1794]  ],## Coronal - HEART
 258             [  [-0.1 , 0.1 , -0.95],[-0.6674 , 0.687 , 0.1794]  ] ] ## Sagittal - HEART
 259
 260    i=0
 261    res=[0,99999]  ## [ <result> , <memory of the last succes calcule> ]
 262    for plane in dicPlane:
 263       i=i+1
 264       refA=plane[0]
 265       refB=plane[1]
 266       res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
 267       res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
 268    return res[0]
 269
 270    except KeyError:
 271    return 0
 272
 273
 274    def VerfCriterion(self,typeCriterion,criterionNew,res):
 275       type = res[0]
 276       criterion = res[1]
 277 #     if criterionNew<0.1 and criterionNew<criterion:
 278       if criterionNew<criterion:
 279          criterion=criterionNew
 280          type=typeCriterion
 281       return [ type , criterion ]
 282
 283
 284    def CalculLikelyhood2Vec(self,refA,refB,ori1,ori2):
 285 """
 286    # ------------------------- Purpose : -----------------------------------
 287    # - This function determine the orientation similarity of two planes.
 288    #   Each plane is described by two vector.
 289    # ------------------------- Parameters : --------------------------------
 290    # - <refA>  : - type : vector 3D (float)
 291    # - <refB>  : - type : vector 3D (float)
 292    #             - Description of the first plane
 293    # - <ori1>  : - type : vector 3D (float)
 294    # - <ori2>  : - type : vector 3D (float)
 295    #             - Description of the second plane
 296    # ------------------------- Return : ------------------------------------
 297    #  float :   0 if the planes are perpendicular.
 298    # While the difference of the orientation between the planes
 299    # are big more enlarge is
 300    # the criterion.
 301    # ------------------------- Other : -------------------------------------
 302    #  The calculus is based with vectors normalice
 303    """
 304
 305       ori3=self.ProductVectorial(ori1,ori2)
 306       refC=self.ProductVectorial(refA,refB)
 307       res=math.pow(refC[0]-ori3[0],2) + math.pow(refC[1]-ori3[1],2) + math.pow(refC[2]-ori3[2],2)
 308       return math.sqrt(res)
 309
 310    def ProductVectorial(self,vec1,vec2):
 311       """
 312       # ------------------------- Purpose : -----------------------------------
 313       # - Calculus of the poduct vectorial between two vectors 3D
 314       # ------------------------- Parameters : --------------------------------
 315       # - <vec1>  : - type : vector 3D (float)
 316       # - <vec2>  : - type : vector 3D (float)
 317       # ------------------------- Return : ------------------------------------
 318       #  (vec) :    - Vector 3D
 319       # ------------------------- Other : -------------------------------------
 320       """
 321       vec3=[0,0,0]
 322       vec3[0]=vec1[1]*vec2[2] - vec1[2]*vec2[1]
 323       vec3[1]=-( vec1[0]*vec2[2] - vec1[2]*vec2[0])
 324       vec3[2]=vec1[0]*vec2[1] - vec1[1]*vec2[0]
 325       return vec3
 326
 327    def GetValuesDict(image):
 328       """
 329       Returns a dictionnary containing values associated with Field Names
 330       dict["Dicom Field Name"]="Dicom field value"
 331       """
 332       val=image.GetFirstEntry()
 333       dic={}
 334       while(val):
 335          if isinstance(val,gdcm.ValEntryPtr):
 336             dic[val.GetName()]=val.GetValue()
 337          val=image.GetNextEntry()
 338       return dic
 339
 340 */
 341
 342
 343 // ------------------------------------------------------------------------
 344 /*
 345 2.2.2 Orientation of DICOM images
 346
 347
 348 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation
 349 says :
 350
 351 A question that is frequently asked in comp.protocols.dicom is how to determine
 352  which side of an image is which (e.g. left, right) and so on.
 353  The short answer is that for projection radiographs this is specified
 354  explicitly using the "Patient Orientation" attribute, and for cross-sectional
 355  images it needs to be derived from the "Image Orientation (Patient)" direction
 356  cosines. In the standard these are explained as follows:
 357
 358     * "C.7.6.1.1.1 Patient Orientation.
 359                 The Patient Orientation (0020,0020) relative to the image
 360                 plane shall be specified by two values that designate the
 361                 anatomical direction of the positive row axis (left to right)
 362                 and the positive column axis (top to bottom).
 363                 The first entry is the direction of the rows, given by the
 364                 direction of the last pixel in the first row from the first
 365                 pixel in that row.
 366                 The second entry is the direction of the columns, given by
 367                 the direction of the last pixel in the first column from the
 368                 first pixel in that column.
 369                 Anatomical direction shall be designated by the capital
 370                 letters: A (anterior), P (posterior), R (right),L (left),
 371                 H (head), F (foot).
 372                 Each value of the orientation attribute shall contain at
 373                 least one of these characters.
 374                 If refinements in the orientation descriptions are to be
 375                 specified, then they shall be designated by one or two
 376                 additional letters in each value.
 377                 Within each value, the letters shall be ordered with the
 378                 principal orientation designated in the first character."
 379
 380     * "C.7.6.2.1.1 Image Position And Image Orientation.
 381                 The "Image Position (Patient)" (0020,0032) specifies the x, y,
 382                 and z coordinates of the upper left hand corner of the image;
 383                 it is the center of the first voxel transmitted.
 384                 The "Image Orientation (Patient)" (0020,0037) specifies the
 385                 direction cosines of the first row and the first column
 386                 with respect to the patient.
 387                 These Attributes shall be provided as a pair.
 388                 Row value for the x, y, and z axes respectively followed by
 389                 the Column value for the x, y, and z axes respectively.
 390                 The direction of the axes is defined fully by the patient's
 391                 orientation.
 392                 The x-axis is increasing to the left hand side of the patient.
 393                 The y-axis is increasing to the posterior side of the patient.
 394                 The z-axis is increasing toward the head of the patient.
 395                 The patient based coordinate system is a right handed system,
 396                 i.e. the vector cross product of a unit vector along the
 397                 positive x-axis and a unit vector along the positive y-axis
 398                 is equal to a unit vector along the positive z-axis."
 399
 400 Some simple code to take one of the direction cosines (vectors) from the
 401 Image Orientation (Patient) attribute and generate strings equivalent to one
 402 of the values of Patient Orientation looks like this (noting that if the vector
 403 is not aligned exactly with one of the major axes, the resulting string will
 404 have multiple letters in as described under "refinements" in C.7.6.1.1.1):
 405
 406 */
 407
 408 /**
 409  * \brief computes the Patient Orientation relative to the image plane
 410  *          from the 'Image Orientation (Patient)'
 411  *          The first entry is the direction of the rows, given by the
 412  *          direction of the last pixel in the first row from the first
 413  *          pixel in that row.
 414  *          The second entry is the direction of the columns, given by
 415  *          the direction of the last pixel in the first column from the
 416  *          first pixel in that column.
 417  *          Anatomical direction is designated by the capital
 418  *          letters: A (anterior), P (posterior), R (right),L (left),
 419  *          H (head), F (foot).
 420  *          Refinements in the orientation descriptions are designated
 421  *          by one or two additional letters in each value.
 422  *          Use it when "Patient Orientation" (0020,0020) is not found
 423  * @return orientation string as "rawOrientation\columnsOrientation"
 424  */
 425 std::string Orientation::GetOrientation ( File *f )
 426 {
 427    float iop[6];
 428    if ( !f->GetImageOrientationPatient( iop ) )
 429    return GDCM_UNFOUND;
 430
 431    std::string orientation;
 432    orientation = GetSingleOrientation ( iop )
 433                + "\\"
 434                + GetSingleOrientation ( iop + 3 );
 435    return orientation;
 436 }
 437
 438
 439 std::string Orientation::GetSingleOrientation ( float *iop)
 440 {
 441    std::string orientation;
 442
 443    char orientationX = iop[0] < 0 ? 'R' : 'L';
 444    char orientationY = iop[1] < 0 ? 'A' : 'P';
 445    char orientationZ = iop[2] < 0 ? 'F' : 'H';
 446
 447    double absX = iop[0];
 448    if (absX < 0) absX = -absX;
 449       double absY = iop[1];
 450    if (absY < 0) absY = -absY;
 451       double absZ = iop[2];
 452    if (absZ < 0) absZ = -absZ;
 453
 454    for (int i=0; i<3; ++i)
 455    {
 456       if (absX>.0001 && absX>absY && absX>absZ)
 457       {
 458          orientation = orientation + orientationX;
 459          absX=0;
 460        }
 461        else if (absY>.0001 && absY>absX && absY>absZ)
 462        {
 463           orientation = orientation + orientationY;
 464           absY=0;
 465        }
 466        else if (absZ>.0001 && absZ>absX && absZ>absY)
 467        {
 468            orientation = orientation + orientationZ;
 469            absZ=0;
 470        }
 471        else
 472           break;
 473      }
 474    return orientation;
 475 }
 476
 477
 478
 479 } // end namespace gdcm