src/gdcmOrientation.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmOrientation.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/10/03 14:54:16 $
   7   Version:   $Revision: 1.16 $
   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 // 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 (Patient)" (0020,0032) specifies the x, y,
 386                 and z coordinates of the upper left hand corner of the image;
 387                 it is the center of the first voxel transmitted.
 388                 The "Image Orientation (Patient)" (0020,0037) specifies the
 389                 direction cosines of the first row and the first column
 390                 with respect to the patient.
 391                 These Attributes shall be provided as a pair.
 392                 Row value for the x, y, and z axes respectively followed by
 393                 the Column value for the x, y, and z axes respectively.
 394                 The direction of the axes is defined fully by the patient's
 395                 orientation.
 396                 The x-axis is increasing to the left hand side of the patient.
 397                 The y-axis is increasing to the posterior side of the patient.
 398                 The z-axis is increasing toward the head of the patient.
 399                 The patient based coordinate system is a right handed system,
 400                 i.e. the vector cross product of a unit vector along the
 401                 positive x-axis and a unit vector along the positive y-axis
 402                 is equal to a unit vector along the positive z-axis."
 403
 404 Some simple code to take one of the direction cosines (vectors) from the
 405 Image Orientation (Patient) attribute and generate strings equivalent to one
 406 of the values of Patient Orientation looks like this (noting that if the vector
 407 is not aligned exactly with one of the major axes, the resulting string will
 408 have multiple letters in as described under "refinements" in C.7.6.1.1.1):
 409
 410 */
 411
 412 /**
 413  * \brief computes the Patient Orientation relative to the image plane
 414  *          from the 'Image Orientation (Patient)'
 415  *          The first entry is the direction of the rows, given by the
 416  *          direction of the last pixel in the first row from the first
 417  *          pixel in that row.
 418  *          The second entry is the direction of the columns, given by
 419  *          the direction of the last pixel in the first column from the
 420  *          first pixel in that column.
 421  *          Anatomical direction is designated by the capital
 422  *          letters: A (anterior), P (posterior), R (right),L (left),
 423  *          H (head), F (foot).
 424  *          Refinements in the orientation descriptions are designated
 425  *          by one or two additional letters in each value.
 426  *          Use it when "Patient Orientation" (0020,0020) is not found
 427  * @return orientation string as "rawOrientation\columnsOrientation"
 428  */
 429 std::string Orientation::GetOrientation ( File *f )
 430 {
 431    float iop[6];
 432    if ( !f->GetImageOrientationPatient( iop ) )
 433    return GDCM_UNFOUND;
 434
 435    std::string orientation;
 436    orientation = GetSingleOrientation ( iop )
 437                + "\\"
 438                + GetSingleOrientation ( iop + 3 );
 439    return orientation;
 440 }
 441
 442
 443 std::string Orientation::GetSingleOrientation ( float *iop)
 444 {
 445    std::string orientation;
 446
 447    char orientationX = iop[0] < 0 ? 'R' : 'L';
 448    char orientationY = iop[1] < 0 ? 'A' : 'P';
 449    char orientationZ = iop[2] < 0 ? 'F' : 'H';
 450
 451    double absX = iop[0];
 452    if (absX < 0) absX = -absX;
 453       double absY = iop[1];
 454    if (absY < 0) absY = -absY;
 455       double absZ = iop[2];
 456    if (absZ < 0) absZ = -absZ;
 457
 458    for (int i=0; i<3; ++i)
 459    {
 460       if (absX>.0001 && absX>absY && absX>absZ)
 461       {
 462          orientation = orientation + orientationX;
 463          absX=0;
 464        }
 465        else if (absY>.0001 && absY>absX && absY>absZ)
 466        {
 467           orientation = orientation + orientationY;
 468           absY=0;
 469        }
 470        else if (absZ>.0001 && absZ>absX && absZ>absY)
 471        {
 472            orientation = orientation + orientationZ;
 473            absZ=0;
 474        }
 475        else
 476           break;
 477      }
 478    return orientation;
 479 }
 480
 481
 482 /*-------------------------------------------------------------------
 483
 484 Some more stuff, from XMedcon
 485
 486 ---> Personal remark from JPRx :
 487 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP "
 488 --> the cosines may ahave any value -1< <+1, for MR images !
 489
 490 enum patient_slice_orientation_type
 491   {
 492     patient_slice_orientation_unknown = 0,
 493     supine_headfirst_transaxial,
 494     supine_headfirst_sagittal,
 495     supine_headfirst_coronal,
 496     supine_feetfirst_transaxial,
 497     supine_feetfirst_sagittal,
 498     supine_feetfirst_coronal,
 499     prone_headfirst_transaxial,
 500     prone_headfirst_sagittal,
 501     prone_headfirst_coronal,
 502     prone_feetfirst_transaxial,
 503     prone_feetfirst_sagittal,
 504     prone_feetfirst_coronal
 505   };
 506
 507 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
 508 {
 509   h.GetImageOrientationPatient(image_orientation_patient);
 510 }
 511
 512 #if 0
 513 //
 514 // this is all completely cribbed from the xmedcon library, since
 515 // we're trying to do what it does, mostly.
 516 patient_slice_orientation_type
 517 GetPatSliceOrient(gdcm::File &h)
 518 {
 519   F32 image_orientation_patient[6];
 520
 521   // protected, do it the hard way
 522   //  h.GetImageOrientationPatient(image_orientation_patient);
 523   GetImageOrientationPatient(h,image_orientation_patient);
 524
 525   enum { headfirst, feetfirst } patient_orientation;
 526   enum { supine, prone } patient_rotation;
 527   enum { transaxial, sagittal, coronal } slice_orientation;
 528
 529   std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
 530   if(patient_position == "gdcm::Unfound")
 531     {
 532     patient_position = "HF";
 533     }
 534   if(patient_position.find("HF") != std::string::npos)
 535     {
 536     patient_orientation = headfirst;
 537     }
 538   else if(patient_position.find("FF") != std::string::npos)
 539     {
 540     patient_orientation = feetfirst;
 541     }
 542
 543   if(patient_position.find("S") != std::string::npos)
 544     {
 545     patient_rotation = supine;
 546     }
 547   else if(patient_position.find("P") != std::string::npos)
 548     {
 549     patient_rotation = prone;
 550     }
 551
 552   if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
 553      (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
 554     {
 555     slice_orientation = transaxial;
 556     }
 557   else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
 558           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
 559     {
 560     slice_orientation = sagittal;
 561     }
 562   else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
 563           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
 564     {
 565     slice_orientation = coronal;
 566     }
 567   //
 568   // combine results
 569   patient_slice_orientation_type patient_slice_orientation =
 570     patient_slice_orientation_unknown;
 571   switch (patient_rotation)
 572     {
 573     case supine:
 574       switch (patient_orientation)
 575         {
 576         case headfirst:
 577           switch (slice_orientation)
 578             {
 579             case transaxial:
 580               patient_slice_orientation = supine_headfirst_transaxial;
 581               break;
 582             case sagittal:
 583               patient_slice_orientation = supine_headfirst_sagittal;
 584               break;
 585             case coronal:
 586               patient_slice_orientation = supine_headfirst_coronal;
 587               break;
 588             }
 589           break;
 590         case feetfirst:
 591           switch (slice_orientation)
 592             {
 593             case transaxial:
 594               patient_slice_orientation = supine_feetfirst_transaxial;
 595               break;
 596             case sagittal:
 597               patient_slice_orientation = supine_feetfirst_sagittal;
 598               break;
 599             case coronal:
 600               patient_slice_orientation = supine_feetfirst_coronal;
 601               break;
 602             }
 603           break;
 604         }
 605       break;
 606     case prone:
 607       switch (patient_orientation)
 608         {
 609         case headfirst:
 610           switch (slice_orientation)
 611             {
 612             case transaxial:
 613               patient_slice_orientation = prone_headfirst_transaxial;
 614               break;
 615             case sagittal:
 616               patient_slice_orientation = prone_headfirst_sagittal;
 617               break;
 618             case coronal:
 619               patient_slice_orientation = prone_headfirst_coronal;
 620               break;
 621             }
 622           break;
 623         case feetfirst:
 624           switch (slice_orientation)
 625             {
 626             case transaxial:
 627               patient_slice_orientation = prone_feetfirst_transaxial;
 628               break;
 629             case sagittal:
 630               patient_slice_orientation = prone_feetfirst_sagittal;
 631               break;
 632             case coronal:
 633               patient_slice_orientation = prone_feetfirst_coronal;
 634               break;
 635             }
 636           break;
 637         }
 638       break;
 639     }
 640   if(patient_slice_orientation != patient_slice_orientation_unknown)
 641     return patient_slice_orientation;
 642   //
 643   // this is what xmedcon does
 644   if ((image_orientation_patient[0] == +1)   &&
 645       (image_orientation_patient[4] == +1))
 646     patient_slice_orientation = supine_headfirst_transaxial;
 647   else if ((image_orientation_patient[0] == -1)   &&
 648            (image_orientation_patient[4] == +1))
 649     patient_slice_orientation = supine_feetfirst_transaxial;
 650   else if ((image_orientation_patient[0] == -1)   &&
 651            (image_orientation_patient[4] == -1))
 652     patient_slice_orientation = prone_headfirst_transaxial;
 653   else if ((image_orientation_patient[0] == +1)   &&
 654            (image_orientation_patient[4] == -1))
 655     patient_slice_orientation = prone_feetfirst_transaxial;
 656
 657   else if ((image_orientation_patient[1] == +1)   &&
 658            (image_orientation_patient[5] == -1))
 659     patient_slice_orientation = supine_headfirst_sagittal;
 660   else if ((image_orientation_patient[1] == +1)   &&
 661            (image_orientation_patient[5] == +1))
 662     patient_slice_orientation = supine_feetfirst_sagittal;
 663   else if ((image_orientation_patient[1] == -1)   &&
 664            (image_orientation_patient[5] == -1))
 665     patient_slice_orientation = prone_headfirst_sagittal;
 666   else if ((image_orientation_patient[1] == -1)   &&
 667            (image_orientation_patient[5] == +1))
 668     patient_slice_orientation = prone_feetfirst_sagittal;
 669
 670   else if ((image_orientation_patient[0] == +1)   &&
 671            (image_orientation_patient[5] == -1))
 672     patient_slice_orientation = supine_headfirst_coronal;
 673   else if ((image_orientation_patient[0] == -1)   &&
 674            (image_orientation_patient[5] == +1))
 675     patient_slice_orientation = supine_feetfirst_coronal;
 676   else if ((image_orientation_patient[0] == -1)   &&
 677            (image_orientation_patient[5] == -1))
 678     patient_slice_orientation = prone_headfirst_coronal;
 679   else if ((image_orientation_patient[0] == +1)   &&
 680            (image_orientation_patient[5] == +1))
 681     patient_slice_orientation = prone_feetfirst_coronal;
 682   return patient_slice_orientation;
 683 }
 684 #else
 685
 686 -------------------------------------------------------------------------*/
 687
 688 } // end namespace gdcm