src/gdcmOrientation.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmOrientation.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/09/20 14:13:41 $
   7   Version:   $Revision: 1.9 $
   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 } // end namespace gdcm
 207
 208
 209
 210
 211 // ---------------------------------------------------------------------------
 212 // Here is the original Python code, kindly supplied by THERALYS
 213 //
 214 // C++ code doesn't give good results
 215 // --> FIXME
 216
 217 /*
 218
 219 def TypeOrientation(self,file0):
 220 """
 221 # ------------------------- Purpose : -----------------------------------
 222 # - This function compare the orientation of the given image and the
 223 #   basics orientations (Axial, Cornal, Sagital)
 224 # ------------------------- Parameters : --------------------------------
 225 # - <file0> : - type : string
 226 #             - The name of the first image file of the serie
 227 # ------------------------- Return : ------------------------------------
 228 #   1 :   Axial
 229 #  -1 :   Axial invert
 230 #   2 :   Coronal
 231 #  -2 :   Coronal invert
 232 #   3 :   Sagital
 233 #  -3 :   Sagital invert
 234 #   4 :   Heart Axial
 235 #  -4 :   Heart Axial invert
 236 #   5 :   Heart Coronal
 237 #  -5 :   Heart Coronal invert
 238 #   6 :   Heart Sagital
 239 #  -6 :   Heart Sagital invert
 240 #
 241    # ------------------------- Other : -------------------------------------
 242 # This method finds the most similar basic orientation.
 243 """
 244 try:
 245    toRead = gdcm.File(file0)
 246    ValDict = GetValuesDict(toRead)
 247    try:
 248       imageOrientation=ValDict["Image Orientation (Patient)"]
 249    except KeyError:
 250       imageOrientation=ValDict["Image Orientation"]
 251
 252    ori1=[float(split(imageOrientation,"\\")[0]),\
 253       float(split(imageOrientation,"\\")[1]),\
 254       float(split(imageOrientation,"\\")[2])]
 255    ori2=[float(split(imageOrientation,"\\")[3]),\
 256       float(split(imageOrientation,"\\")[4]),\
 257       float(split(imageOrientation,"\\")[5])]
 258
 259 ## two vectors perpendicular describe one plane
 260    dicPlane=[ [  [1,0,0],[0,1,0]   ],  ## Axial
 261             [  [1,0,0],[0,0,-1]  ],  ## Coronal
 262             [  [0,1,0],[0,0,-1]  ],  ## Sagittal
 263             [  [ 0.8 , 0.5 ,  0.0 ],[-0.1 , 0.1 , -0.95]        ],## Axial - HEART
 264             [  [ 0.8 , 0.5 ,  0.0 ],[-0.6674 , 0.687 , 0.1794]  ],## Coronal - HEART
 265             [  [-0.1 , 0.1 , -0.95],[-0.6674 , 0.687 , 0.1794]  ] ] ## Sagittal - HEART
 266
 267    i=0
 268    res=[0,99999]  ## [ <result> , <memory of the last succes calcule> ]
 269    for plane in dicPlane:
 270       i=i+1
 271       refA=plane[0]
 272       refB=plane[1]
 273       res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
 274       res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
 275    return res[0]
 276
 277    except KeyError:
 278    return 0
 279
 280
 281    def VerfCriterion(self,typeCriterion,criterionNew,res):
 282       type = res[0]
 283       criterion = res[1]
 284 #     if criterionNew<0.1 and criterionNew<criterion:
 285       if criterionNew<criterion:
 286          criterion=criterionNew
 287          type=typeCriterion
 288       return [ type , criterion ]
 289
 290
 291    def CalculLikelyhood2Vec(self,refA,refB,ori1,ori2):
 292 """
 293    # ------------------------- Purpose : -----------------------------------
 294    # - This function determine the orientation similarity of two planes.
 295    #   Each plane is described by two vector.
 296    # ------------------------- Parameters : --------------------------------
 297    # - <refA>  : - type : vector 3D (float)
 298    # - <refB>  : - type : vector 3D (float)
 299    #             - Description of the first plane
 300    # - <ori1>  : - type : vector 3D (float)
 301    # - <ori2>  : - type : vector 3D (float)
 302    #             - Description of the second plane
 303    # ------------------------- Return : ------------------------------------
 304    #  float :   0 if the planes are perpendicular.
 305    # While the difference of the orientation between the planes
 306    # are big more enlarge is
 307    # the criterion.
 308    # ------------------------- Other : -------------------------------------
 309    #  The calculus is based with vectors normalice
 310    """
 311
 312       ori3=self.ProductVectorial(ori1,ori2)
 313       refC=self.ProductVectorial(refA,refB)
 314       res=math.pow(refC[0]-ori3[0],2) + math.pow(refC[1]-ori3[1],2) + math.pow(refC[2]-ori3[2],2)
 315       return math.sqrt(res)
 316
 317    def ProductVectorial(self,vec1,vec2):
 318       """
 319       # ------------------------- Purpose : -----------------------------------
 320       # - Calculus of the poduct vectorial between two vectors 3D
 321       # ------------------------- Parameters : --------------------------------
 322       # - <vec1>  : - type : vector 3D (float)
 323       # - <vec2>  : - type : vector 3D (float)
 324       # ------------------------- Return : ------------------------------------
 325       #  (vec) :    - Vector 3D
 326       # ------------------------- Other : -------------------------------------
 327       """
 328       vec3=[0,0,0]
 329       vec3[0]=vec1[1]*vec2[2] - vec1[2]*vec2[1]
 330       vec3[1]=-( vec1[0]*vec2[2] - vec1[2]*vec2[0])
 331       vec3[2]=vec1[0]*vec2[1] - vec1[1]*vec2[0]
 332       return vec3
 333
 334    def GetValuesDict(image):
 335       """
 336       Returns a dictionnary containing values associated with Field Names
 337       dict["Dicom Field Name"]="Dicom field value"
 338       """
 339       val=image.GetFirstEntry()
 340       dic={}
 341       while(val):
 342          if isinstance(val,gdcm.ValEntryPtr):
 343             dic[val.GetName()]=val.GetValue()
 344          val=image.GetNextEntry()
 345       return dic
 346
 347 */