src/gdcmOrientation.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmOrientation.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/09/19 09:48:27 $
   7   Version:   $Revision: 1.8 $
   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    for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
  92    {
  93        ++i;
  94        // refA=plane[0]
  95        refA.x = dicPlane[numDicPlane][0][0];
  96        refA.y = dicPlane[numDicPlane][0][1];
  97        refA.z = dicPlane[numDicPlane][0][2];
  98        // refB=plane[1]
  99        refB.x = dicPlane[numDicPlane][1][0];
 100        refB.y = dicPlane[numDicPlane][1][1];
 101        refB.z = dicPlane[numDicPlane][1][2];
 102        res=VerfCriterion(  i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
 103        res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
 104    }
 105    return res.first;
 106 /*
 107 //  i=0
 108 //  res=[0,99999]  ## [ <result> , <memory of the last succes calculus> ]
 109 //  for plane in dicPlane:
 110 //      i=i+1
 111 //      refA=plane[0]
 112 //      refB=plane[1]
 113 //      res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
 114 //      res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
 115 //  return res[0]
 116 */
 117 }
 118
 119 Res
 120 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
 121 {
 122    Res res;
 123    double criterion = in.second;
 124    if (criterionNew < criterion)
 125    {
 126       res.first  = typeCriterion;;
 127       res.second = criterionNew;
 128    }
 129 /*
 130 //   type = res[0]
 131 //   criterion = res[1]
 132 // #     if criterionNew<0.1 and criterionNew<criterion:
 133 //   if criterionNew<criterion:
 134 //      criterion=criterionNew
 135 //      type=typeCriterion
 136 //   return [ type , criterion ]
 137 */
 138    return res;
 139 }
 140
 141 inline double square_dist(vector3D const &v1, vector3D const &v2)
 142 {
 143   double res;
 144   res = (v1.x - v2.x)*(v1.x - v2.x) +
 145         (v1.y - v2.y)*(v1.y - v2.y) +
 146         (v1.z - v2.z)*(v1.z - v2.z);
 147   return res;
 148 }
 149
 150 //------------------------- Purpose : -----------------------------------
 151 //- This function determines the orientation similarity of two planes.
 152 //  Each plane is described by two vectors.
 153 //------------------------- Parameters : --------------------------------
 154 //- <refA>  : - type : vector 3D (double)
 155 //- <refB>  : - type : vector 3D (double)
 156 //            - Description of the first plane
 157 //- <ori1>  : - type : vector 3D (double)
 158 //- <ori2>  : - type : vector 3D (double)
 159 //            - Description of the second plane
 160 //------------------------- Return : ------------------------------------
 161 // double :   0 if the planes are perpendicular. While the difference of
 162 //            the orientation between the planes are big more enlarge is
 163 //            the criterion.
 164 //------------------------- Other : -------------------------------------
 165 // The calculus is based with vectors normalice
 166 double
 167 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB,
 168                                   vector3D const &ori1, vector3D const &ori2 )
 169 {
 170
 171    vector3D ori3 = ProductVectorial(ori1,ori2);
 172    vector3D refC = ProductVectorial(refA,refB);
 173    double res = square_dist(refC, ori3);
 174
 175    return sqrt(res);
 176 }
 177
 178 //------------------------- Purpose : -----------------------------------
 179 //- Calculus of the poduct vectorial between two vectors 3D
 180 //------------------------- Parameters : --------------------------------
 181 //- <vec1>  : - type : vector 3D (double)
 182 //- <vec2>  : - type : vector 3D (double)
 183 //------------------------- Return : ------------------------------------
 184 // (vec) :    - Vector 3D
 185 //------------------------- Other : -------------------------------------
 186 vector3D
 187 Orientation::ProductVectorial(vector3D const & vec1, vector3D const & vec2)
 188 {
 189    vector3D vec3;
 190    vec3.x =    vec1.y*vec2.z - vec1.z*vec2.y;
 191    vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x);
 192    vec3.z =    vec1.x*vec2.y - vec1.y*vec2.x;
 193
 194    return vec3;
 195 }
 196
 197 } // end namespace gdcm
 198
 199
 200
 201
 202 // ---------------------------------------------------------------------------
 203 // Here is the original Python code, kindly supplied by THERALYS
 204 //
 205 // C++ code doesn't give good results
 206 // --> FIXME
 207
 208 /*
 209
 210 def TypeOrientation(self,file0):
 211 """
 212 # ------------------------- Purpose : -----------------------------------
 213 # - This function compare the orientation of the given image and the
 214 #   basics orientations (Axial, Cornal, Sagital)
 215 # ------------------------- Parameters : --------------------------------
 216 # - <file0> : - type : string
 217 #             - The name of the first image file of the serie
 218 # ------------------------- Return : ------------------------------------
 219 #   1 :   Axial
 220 #  -1 :   Axial invert
 221 #   2 :   Coronal
 222 #  -2 :   Coronal invert
 223 #   3 :   Sagital
 224 #  -3 :   Sagital invert
 225 #   4 :   Heart Axial
 226 #  -4 :   Heart Axial invert
 227 #   5 :   Heart Coronal
 228 #  -5 :   Heart Coronal invert
 229 #   6 :   Heart Sagital
 230 #  -6 :   Heart Sagital invert
 231 #
 232    # ------------------------- Other : -------------------------------------
 233 # This method finds the most similar basic orientation.
 234 """
 235 try:
 236    toRead = gdcm.File(file0)
 237    ValDict = GetValuesDict(toRead)
 238    try:
 239       imageOrientation=ValDict["Image Orientation (Patient)"]
 240    except KeyError:
 241       imageOrientation=ValDict["Image Orientation"]
 242
 243    ori1=[float(split(imageOrientation,"\\")[0]),\
 244       float(split(imageOrientation,"\\")[1]),\
 245       float(split(imageOrientation,"\\")[2])]
 246    ori2=[float(split(imageOrientation,"\\")[3]),\
 247       float(split(imageOrientation,"\\")[4]),\
 248       float(split(imageOrientation,"\\")[5])]
 249
 250 ## two vectors perpendicular describe one plane
 251    dicPlane=[ [  [1,0,0],[0,1,0]   ],  ## Axial
 252             [  [1,0,0],[0,0,-1]  ],  ## Coronal
 253             [  [0,1,0],[0,0,-1]  ],  ## Sagittal
 254             [  [ 0.8 , 0.5 ,  0.0 ],[-0.1 , 0.1 , -0.95]        ],## Axial - HEART
 255             [  [ 0.8 , 0.5 ,  0.0 ],[-0.6674 , 0.687 , 0.1794]  ],## Coronal - HEART
 256             [  [-0.1 , 0.1 , -0.95],[-0.6674 , 0.687 , 0.1794]  ] ] ## Sagittal - HEART
 257
 258    i=0
 259    res=[0,99999]  ## [ <result> , <memory of the last succes calcule> ]
 260    for plane in dicPlane:
 261       i=i+1
 262       refA=plane[0]
 263       refB=plane[1]
 264       res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
 265       res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
 266    return res[0]
 267
 268    except KeyError:
 269    return 0
 270
 271
 272    def VerfCriterion(self,typeCriterion,criterionNew,res):
 273       type = res[0]
 274       criterion = res[1]
 275 #     if criterionNew<0.1 and criterionNew<criterion:
 276       if criterionNew<criterion:
 277          criterion=criterionNew
 278          type=typeCriterion
 279       return [ type , criterion ]
 280
 281
 282    def CalculLikelyhood2Vec(self,refA,refB,ori1,ori2):
 283 """
 284    # ------------------------- Purpose : -----------------------------------
 285    # - This function determine the orientation similarity of two planes.
 286    #   Each plane is described by two vector.
 287    # ------------------------- Parameters : --------------------------------
 288    # - <refA>  : - type : vector 3D (float)
 289    # - <refB>  : - type : vector 3D (float)
 290    #             - Description of the first plane
 291    # - <ori1>  : - type : vector 3D (float)
 292    # - <ori2>  : - type : vector 3D (float)
 293    #             - Description of the second plane
 294    # ------------------------- Return : ------------------------------------
 295    #  float :   0 if the planes are perpendicular.
 296    # While the difference of the orientation between the planes
 297    # are big more enlarge is
 298    # the criterion.
 299    # ------------------------- Other : -------------------------------------
 300    #  The calculus is based with vectors normalice
 301    """
 302
 303       ori3=self.ProductVectorial(ori1,ori2)
 304       refC=self.ProductVectorial(refA,refB)
 305       res=math.pow(refC[0]-ori3[0],2) + math.pow(refC[1]-ori3[1],2) + math.pow(refC[2]-ori3[2],2)
 306       return math.sqrt(res)
 307
 308    def ProductVectorial(self,vec1,vec2):
 309       """
 310       # ------------------------- Purpose : -----------------------------------
 311       # - Calculus of the poduct vectorial between two vectors 3D
 312       # ------------------------- Parameters : --------------------------------
 313       # - <vec1>  : - type : vector 3D (float)
 314       # - <vec2>  : - type : vector 3D (float)
 315       # ------------------------- Return : ------------------------------------
 316       #  (vec) :    - Vector 3D
 317       # ------------------------- Other : -------------------------------------
 318       """
 319       vec3=[0,0,0]
 320       vec3[0]=vec1[1]*vec2[2] - vec1[2]*vec2[1]
 321       vec3[1]=-( vec1[0]*vec2[2] - vec1[2]*vec2[0])
 322       vec3[2]=vec1[0]*vec2[1] - vec1[1]*vec2[0]
 323       return vec3
 324
 325    def GetValuesDict(image):
 326       """
 327       Returns a dictionnary containing values associated with Field Names
 328       dict["Dicom Field Name"]="Dicom field value"
 329       """
 330       val=image.GetFirstEntry()
 331       dic={}
 332       while(val):
 333          if isinstance(val,gdcm.ValEntryPtr):
 334             dic[val.GetName()]=val.GetValue()
 335          val=image.GetNextEntry()
 336       return dic
 337
 338 */