src/gdcmOrientation.cxx

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