src/gdcmOrientation.cxx

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