src/gdcmOrientation.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmOrientation.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2007/05/23 14:18:10 $
   7   Version:   $Revision: 1.25 $
   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_NAME_SPACE
  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
  55 static const char  *OrientationTypeStrings[] = {
  56   "Not Applicable",
  57   "Axial",
  58   "Coronal",
  59   "Sagital",
  60   "Heart Axial",
  61   "Heart Coronal",
  62   "Heart Sagital",
  63   "Axial invert",
  64   "Coronal invert",
  65   "Sagital invert",
  66   "Heart Axial invert",
  67   "Heart Coronal invert",
  68   "Heart Sagital invert",
  69   NULL
  70 };
  71
  72 /// \brief returns human readable interpretation of the most
  73 ///        similar basic orientation (Axial, Coronal, Sagital, ...) of the image
  74 const char* Orientation::GetOrientationTypeString(OrientationType const o)
  75 {
  76   int k = (int)o;
  77   if (k < 0)
  78        k = -k + 6;
  79
  80   return OrientationTypeStrings[k];
  81 }
  82
  83 /// \brief returns of the most similar basic orientation
  84 ///        (Axial, Coronal, Sagital, ...) of the image
  85 OrientationType Orientation::GetOrientationType( File *f )
  86 {
  87    float iop[6];
  88    bool succ = f->GetImageOrientationPatient( iop );
  89    if ( !succ )
  90    {
  91       gdcmWarningMacro( "No Image Orientation (0020,0037)/(0020,0032) found in the file, cannot proceed." )
  92       return NotApplicable;
  93    }
  94    vector3D ori1;
  95    vector3D ori2;
  96
  97    ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2];
  98    ori2.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5];
  99
 100    // two perpendicular vectors describe one plane
 101    double dicPlane[6][2][3] =
 102    { {  { 1,   0,    0   },{ 0,      1,     0     }  }, // Axial
 103      {  { 1,   0,    0   },{ 0,      0,    -1     }  }, // Coronal
 104      {  { 0,   1,    0   },{ 0,      0,    -1     }  }, // Sagittal
 105      {  { 0.8, 0.5,  0.0 },{-0.1,    0.1 , -0.95  }  }, // Axial - HEART
 106      {  { 0.8, 0.5,  0.0 },{-0.6674, 0.687, 0.1794}  }, // Coronal - HEART
 107      {  {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794}  }  // Sagittal - HEART
 108    };
 109
 110    vector3D refA;
 111    vector3D refB;
 112    int i = 0;
 113    Res res;   // [ <result> , <memory of the last succes calcule> ]
 114    res.first = 0;
 115    res.second = 99999;
 116
 117    for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
 118    {
 119        ++i;
 120        // refA=plane[0]
 121        refA.x = dicPlane[numDicPlane][0][0];
 122        refA.y = dicPlane[numDicPlane][0][1];
 123        refA.z = dicPlane[numDicPlane][0][2];
 124        // refB=plane[1]
 125        refB.x = dicPlane[numDicPlane][1][0];
 126        refB.y = dicPlane[numDicPlane][1][1];
 127        refB.z = dicPlane[numDicPlane][1][2];
 128        res=VerfCriterion(  i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
 129        res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
 130    }
 131    // res thought looks like is a float value, but is indeed an int
 132    // try casting it to int first then enum value to please VS7:
 133    int int_res = (int)res.first;
 134    gdcmAssertMacro( int_res <= 6 && int_res >= -6);
 135    return (OrientationType)int_res;
 136 }
 137
 138 Res
 139 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
 140 {
 141    Res res;
 142    double type = in.first;
 143    double criterion = in.second;
 144    if (/*criterionNew < 0.1 && */criterionNew < criterion)
 145    {
 146       type      = typeCriterion;
 147       criterion = criterionNew;
 148    }
 149    res.first  = type;
 150    res.second = criterion;
 151    return res;
 152 }
 153
 154 inline double square_dist(vector3D const &v1, vector3D const &v2)
 155 {
 156   double res;
 157   res = (v1.x - v2.x)*(v1.x - v2.x) +
 158         (v1.y - v2.y)*(v1.y - v2.y) +
 159         (v1.z - v2.z)*(v1.z - v2.z);
 160   return res;
 161 }
 162
 163 //------------------------- Purpose : -----------------------------------
 164 //- This function determines the orientation similarity of two planes.
 165 //  Each plane is described by two vectors.
 166 //------------------------- Parameters : --------------------------------
 167 //- <refA>  : - type : vector 3D (double)
 168 //- <refB>  : - type : vector 3D (double)
 169 //            - Description of the first plane
 170 //- <ori1>  : - type : vector 3D (double)
 171 //- <ori2>  : - type : vector 3D (double)
 172 //            - Description of the second plane
 173 //------------------------- Return : ------------------------------------
 174 // double :   0 if the planes are perpendicular. While the difference of
 175 //            the orientation between the planes are big more enlarge is
 176 //            the criterion.
 177 //------------------------- Other : -------------------------------------
 178 // The calculus is based with vectors normalice
 179 double
 180 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB,
 181                                   vector3D const &ori1, vector3D const &ori2 )
 182 {
 183
 184    vector3D ori3 = ProductVectorial(ori1,ori2);
 185    vector3D refC = ProductVectorial(refA,refB);
 186    double res = square_dist(refC, ori3);
 187
 188    return sqrt(res);
 189 }
 190
 191 //------------------------- Purpose : -----------------------------------
 192 //- Calculus of the poduct vectorial between two vectors 3D
 193 //------------------------- Parameters : --------------------------------
 194 //- <vec1>  : - type : vector 3D (double)
 195 //- <vec2>  : - type : vector 3D (double)
 196 //------------------------- Return : ------------------------------------
 197 // (vec) :    - Vector 3D
 198 //------------------------- Other : -------------------------------------
 199 vector3D
 200 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2)
 201 {
 202    vector3D vec3;
 203    vec3.x =    vec1.y*vec2.z - vec1.z*vec2.y;
 204    vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x);
 205    vec3.z =    vec1.x*vec2.y - vec1.y*vec2.x;
 206
 207    return vec3;
 208 }
 209
 210
 211 // ------------------------------------------------------------------------
 212 /*
 213 2.2.2 Orientation of DICOM images
 214
 215
 216 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation
 217 says :
 218
 219 A question that is frequently asked in comp.protocols.dicom is how to determine
 220  which side of an image is which (e.g. left, right) and so on.
 221  The short answer is that for projection radiographs this is specified
 222  explicitly using the "Patient Orientation" attribute, and for cross-sectional
 223  images it needs to be derived from the "Image Orientation (Patient)" direction
 224  cosines. In the standard these are explained as follows:
 225
 226     * "C.7.6.1.1.1 Patient Orientation.
 227                 The Patient Orientation (0020,0020) relative to the image
 228                 plane shall be specified by two values that designate the
 229                 anatomical direction of the positive row axis (left to right)
 230                 and the positive column axis (top to bottom).
 231                 The first entry is the direction of the rows, given by the
 232                 direction of the last pixel in the first row from the first
 233                 pixel in that row.
 234                 The second entry is the direction of the columns, given by
 235                 the direction of the last pixel in the first column from the
 236                 first pixel in that column.
 237                 Anatomical direction shall be designated by the capital
 238                 letters: A (anterior), P (posterior), R (right),L (left),
 239                 H (head), F (foot).
 240                 Each value of the orientation attribute shall contain at
 241                 least one of these characters.
 242                 If refinements in the orientation descriptions are to be
 243                 specified, then they shall be designated by one or two
 244                 additional letters in each value.
 245                 Within each value, the letters shall be ordered with the
 246                 principal orientation designated in the first character."
 247
 248     * "C.7.6.2.1.1 Image Position And Image Orientation.
 249                 The "Image Position (Patient)" (0020,0032) specifies the x, y,
 250                 and z coordinates of the upper left hand corner of the image;
 251                 it is the center of the first voxel transmitted.
 252                 The "Image Orientation (Patient)" (0020,0037) specifies the
 253                 direction cosines of the first row and the first column
 254                 with respect to the patient.
 255                 These Attributes shall be provided as a pair.
 256                 Row value for the x, y, and z axes respectively followed by
 257                 the Column value for the x, y, and z axes respectively.
 258                 The direction of the axes is defined fully by the patient's
 259                 orientation.
 260                 The x-axis is increasing to the left hand side of the patient.
 261                 The y-axis is increasing to the posterior side of the patient.
 262                 The z-axis is increasing toward the head of the patient.
 263                 The patient based coordinate system is a right handed system,
 264                 i.e. the vector cross product of a unit vector along the
 265                 positive x-axis and a unit vector along the positive y-axis
 266                 is equal to a unit vector along the positive z-axis."
 267
 268 Some simple code to take one of the direction cosines (vectors) from the
 269 Image Orientation (Patient) attribute and generate strings equivalent to one
 270 of the values of Patient Orientation looks like this (noting that if the vector
 271 is not aligned exactly with one of the major axes, the resulting string will
 272 have multiple letters in as described under "refinements" in C.7.6.1.1.1):
 273
 274 */
 275
 276 /**
 277  * \brief Computes the Patient Orientation relative to the image plane
 278  *          from the 'Image Orientation (Patient)'
 279  *          - or from 0020 0035Image Orientation (RET) -
 280  *          - The first entry is the direction of the rows, given by the
 281  *          direction of the last pixel in the first row from the first
 282  *          pixel in that row.
 283  *          - The second entry is the direction of the columns, given by
 284  *          the direction of the last pixel in the first column from the
 285  *          first pixel in that column.
 286  *          Anatomical direction is designated by the capital
 287  *          letters: A (anterior), P (posterior), R (right),L (left),
 288  *          H (head), F (foot).
 289  *          - Refinements in the orientation descriptions are designated
 290  *          by one or two additional letters in each value.
 291  *          Use it when "Patient Orientation" (0020,0020) is not found
 292  * @return orientation string as "rowsOrientation\columnsOrientation"
 293  */
 294 std::string Orientation::GetOrientation ( File *f )
 295 {
 296    float iop[6];
 297    if ( !f->GetImageOrientationPatient( iop ) )
 298    return GDCM_UNFOUND;
 299
 300    std::string orientation;
 301    orientation = GetSingleOrientation ( iop )
 302                + "\\"
 303                + GetSingleOrientation ( iop + 3 );
 304    return orientation;
 305 }
 306
 307
 308 std::string Orientation::GetSingleOrientation ( float *iop)
 309 {
 310    std::string orientation;
 311
 312    char orientationX = iop[0] < 0 ? 'R' : 'L';
 313    char orientationY = iop[1] < 0 ? 'A' : 'P';
 314    char orientationZ = iop[2] < 0 ? 'F' : 'H';
 315
 316    double absX = iop[0];
 317    if (absX < 0) absX = -absX;
 318       double absY = iop[1];
 319    if (absY < 0) absY = -absY;
 320       double absZ = iop[2];
 321    if (absZ < 0) absZ = -absZ;
 322
 323    for (int i=0; i<3; ++i)
 324    {
 325       if (absX>.0001 && absX>absY && absX>absZ)
 326       {
 327          orientation = orientation + orientationX;
 328          absX=0;
 329        }
 330        else if (absY>.0001 && absY>absX && absY>absZ)
 331        {
 332           orientation = orientation + orientationY;
 333           absY=0;
 334        }
 335        else if (absZ>.0001 && absZ>absX && absZ>absY)
 336        {
 337            orientation = orientation + orientationZ;
 338            absZ=0;
 339        }
 340        else
 341           break;
 342      }
 343    return orientation;
 344 }
 345
 346
 347 /*-------------------------------------------------------------------
 348
 349 Some more stuff, from XMedcon
 350
 351 ---> Personal remark from JPRx :
 352 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP "
 353 --> or, not so common,
 354 // HFDR = Head First-Decubitus Right
 355 // HFDL = Head First-Decubitus Left
 356 // FFDR = Feet First-Decubitus Right
 357 // FFDL = Feet First-Decubitus Left
 358 --> the cosines may have any value -1.< <+1., for MR images !
 359
 360 enum patient_slice_orientation_type
 361   {
 362     patient_slice_orientation_unknown = 0,
 363     supine_headfirst_transaxial,
 364     supine_headfirst_sagittal,
 365     supine_headfirst_coronal,
 366     supine_feetfirst_transaxial,
 367     supine_feetfirst_sagittal,
 368     supine_feetfirst_coronal,
 369     prone_headfirst_transaxial,
 370     prone_headfirst_sagittal,
 371     prone_headfirst_coronal,
 372     prone_feetfirst_transaxial,
 373     prone_feetfirst_sagittal,
 374     prone_feetfirst_coronal
 375   };
 376
 377 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
 378 {
 379   h.GetImageOrientationPatient(image_orientation_patient);
 380 }
 381
 382 #if 0
 383 //
 384 // this is all completely cribbed from the xmedcon library, since
 385 // we're trying to do what it does, mostly.
 386 patient_slice_orientation_type
 387 GetPatSliceOrient(gdcm::File &h)
 388 {
 389   F32 image_orientation_patient[6];
 390
 391   // protected, do it the hard way
 392   //  h.GetImageOrientationPatient(image_orientation_patient);
 393   GetImageOrientationPatient(h,image_orientation_patient);
 394
 395   enum { headfirst, feetfirst } patient_orientation;
 396   enum { supine, prone } patient_rotation;
 397   enum { transaxial, sagittal, coronal } slice_orientation;
 398
 399   std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
 400   if(patient_position == "gdcm::Unfound")
 401     {
 402     patient_position = "HF";
 403     }
 404   if(patient_position.find("HF") != std::string::npos)
 405     {
 406     patient_orientation = headfirst;
 407     }
 408   else if(patient_position.find("FF") != std::string::npos)
 409     {
 410     patient_orientation = feetfirst;
 411     }
 412
 413   if(patient_position.find("S") != std::string::npos)
 414     {
 415     patient_rotation = supine;
 416     }
 417   else if(patient_position.find("P") != std::string::npos)
 418     {
 419     patient_rotation = prone;
 420     }
 421
 422   if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
 423      (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
 424     {
 425     slice_orientation = transaxial;
 426     }
 427   else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
 428           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
 429     {
 430     slice_orientation = sagittal;
 431     }
 432   else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
 433           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
 434     {
 435     slice_orientation = coronal;
 436     }
 437   //
 438   // combine results
 439   patient_slice_orientation_type patient_slice_orientation =
 440     patient_slice_orientation_unknown;
 441   switch (patient_rotation)
 442     {
 443     case supine:
 444       switch (patient_orientation)
 445         {
 446         case headfirst:
 447           switch (slice_orientation)
 448             {
 449             case transaxial:
 450               patient_slice_orientation = supine_headfirst_transaxial;
 451               break;
 452             case sagittal:
 453               patient_slice_orientation = supine_headfirst_sagittal;
 454               break;
 455             case coronal:
 456               patient_slice_orientation = supine_headfirst_coronal;
 457               break;
 458             }
 459           break;
 460         case feetfirst:
 461           switch (slice_orientation)
 462             {
 463             case transaxial:
 464               patient_slice_orientation = supine_feetfirst_transaxial;
 465               break;
 466             case sagittal:
 467               patient_slice_orientation = supine_feetfirst_sagittal;
 468               break;
 469             case coronal:
 470               patient_slice_orientation = supine_feetfirst_coronal;
 471               break;
 472             }
 473           break;
 474         }
 475       break;
 476     case prone:
 477       switch (patient_orientation)
 478         {
 479         case headfirst:
 480           switch (slice_orientation)
 481             {
 482             case transaxial:
 483               patient_slice_orientation = prone_headfirst_transaxial;
 484               break;
 485             case sagittal:
 486               patient_slice_orientation = prone_headfirst_sagittal;
 487               break;
 488             case coronal:
 489               patient_slice_orientation = prone_headfirst_coronal;
 490               break;
 491             }
 492           break;
 493         case feetfirst:
 494           switch (slice_orientation)
 495             {
 496             case transaxial:
 497               patient_slice_orientation = prone_feetfirst_transaxial;
 498               break;
 499             case sagittal:
 500               patient_slice_orientation = prone_feetfirst_sagittal;
 501               break;
 502             case coronal:
 503               patient_slice_orientation = prone_feetfirst_coronal;
 504               break;
 505             }
 506           break;
 507         }
 508       break;
 509     }
 510   if(patient_slice_orientation != patient_slice_orientation_unknown)
 511     return patient_slice_orientation;
 512   //
 513   // this is what xmedcon does
 514   if ((image_orientation_patient[0] == +1)   &&
 515       (image_orientation_patient[4] == +1))
 516     patient_slice_orientation = supine_headfirst_transaxial;
 517   else if ((image_orientation_patient[0] == -1)   &&
 518            (image_orientation_patient[4] == +1))
 519     patient_slice_orientation = supine_feetfirst_transaxial;
 520   else if ((image_orientation_patient[0] == -1)   &&
 521            (image_orientation_patient[4] == -1))
 522     patient_slice_orientation = prone_headfirst_transaxial;
 523   else if ((image_orientation_patient[0] == +1)   &&
 524            (image_orientation_patient[4] == -1))
 525     patient_slice_orientation = prone_feetfirst_transaxial;
 526
 527   else if ((image_orientation_patient[1] == +1)   &&
 528            (image_orientation_patient[5] == -1))
 529     patient_slice_orientation = supine_headfirst_sagittal;
 530   else if ((image_orientation_patient[1] == +1)   &&
 531            (image_orientation_patient[5] == +1))
 532     patient_slice_orientation = supine_feetfirst_sagittal;
 533   else if ((image_orientation_patient[1] == -1)   &&
 534            (image_orientation_patient[5] == -1))
 535     patient_slice_orientation = prone_headfirst_sagittal;
 536   else if ((image_orientation_patient[1] == -1)   &&
 537            (image_orientation_patient[5] == +1))
 538     patient_slice_orientation = prone_feetfirst_sagittal;
 539
 540   else if ((image_orientation_patient[0] == +1)   &&
 541            (image_orientation_patient[5] == -1))
 542     patient_slice_orientation = supine_headfirst_coronal;
 543   else if ((image_orientation_patient[0] == -1)   &&
 544            (image_orientation_patient[5] == +1))
 545     patient_slice_orientation = supine_feetfirst_coronal;
 546   else if ((image_orientation_patient[0] == -1)   &&
 547            (image_orientation_patient[5] == -1))
 548     patient_slice_orientation = prone_headfirst_coronal;
 549   else if ((image_orientation_patient[0] == +1)   &&
 550            (image_orientation_patient[5] == +1))
 551     patient_slice_orientation = prone_feetfirst_coronal;
 552   return patient_slice_orientation;
 553 }
 554 #else
 555
 556 -------------------------------------------------------------------------*/
 557
 558 } // end namespace gdcm