]> Creatis software - gdcm.git/blob - gdcmOrientation.cxx
7e54d554e5518e01359c97dd5333b6183d12659a
[gdcm.git] / gdcmOrientation.cxx
1 /*=========================================================================
2                                                                                 
3   Program:   gdcm
4   Module:    $RCSfile: gdcmOrientation.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/11/28 11:54:51 $
7   Version:   $Revision: 1.21 $
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 /// \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       gdcmErrorMacro( "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  *          - The first entry is the direction of the rows, given by the 
280  *          direction of the last pixel in the first row from the first 
281  *          pixel in that row. 
282  *          - The second entry is the direction of the columns, given by 
283  *          the direction of the last pixel in the first column from the
284  *          first pixel in that column. 
285  *          Anatomical direction is designated by the capital 
286  *          letters: A (anterior), P (posterior), R (right),L (left), 
287  *          H (head), F (foot).
288  *          - Refinements in the orientation descriptions are designated 
289  *          by one or two additional letters in each value.
290  *          Use it when "Patient Orientation" (0020,0020) is not found 
291  * @return orientation string as "rawOrientation\columnsOrientation"
292  */
293 std::string Orientation::GetOrientation ( File *f )
294 {
295    float iop[6];
296    if ( !f->GetImageOrientationPatient( iop ) )
297    return GDCM_UNFOUND;
298
299    std::string orientation;
300    orientation = GetSingleOrientation ( iop ) 
301                + "\\" 
302                + GetSingleOrientation ( iop + 3 );
303    return orientation;
304 }
305
306
307 std::string Orientation::GetSingleOrientation ( float *iop)
308 {
309    std::string orientation;
310
311    char orientationX = iop[0] < 0 ? 'R' : 'L';
312    char orientationY = iop[1] < 0 ? 'A' : 'P';
313    char orientationZ = iop[2] < 0 ? 'F' : 'H';
314
315    double absX = iop[0];
316    if (absX < 0) absX = -absX;
317       double absY = iop[1];
318    if (absY < 0) absY = -absY;
319       double absZ = iop[2];
320    if (absZ < 0) absZ = -absZ;
321
322    for (int i=0; i<3; ++i) 
323    {
324       if (absX>.0001 && absX>absY && absX>absZ) 
325       {
326          orientation = orientation + orientationX;
327          absX=0;
328        }
329        else if (absY>.0001 && absY>absX && absY>absZ) 
330        {
331           orientation = orientation + orientationY;
332           absY=0;
333        }
334        else if (absZ>.0001 && absZ>absX && absZ>absY) 
335        {
336            orientation = orientation + orientationZ;
337            absZ=0;
338        }
339        else 
340           break;
341      }
342    return orientation;
343
344
345
346 /*-------------------------------------------------------------------
347
348 Some more stuff, from XMedcon
349
350 ---> Personal remark from JPRx :
351 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP " 
352 --> or, not so common, 
353 // HFDR = Head First-Decubitus Right
354 // HFDL = Head First-Decubitus Left
355 // FFDR = Feet First-Decubitus Right
356 // FFDL = Feet First-Decubitus Left
357 --> the cosines may have any value -1.< <+1., for MR images !
358
359 enum patient_slice_orientation_type
360   {
361     patient_slice_orientation_unknown = 0,
362     supine_headfirst_transaxial,
363     supine_headfirst_sagittal,
364     supine_headfirst_coronal,
365     supine_feetfirst_transaxial,
366     supine_feetfirst_sagittal,
367     supine_feetfirst_coronal,
368     prone_headfirst_transaxial,
369     prone_headfirst_sagittal,
370     prone_headfirst_coronal,
371     prone_feetfirst_transaxial,
372     prone_feetfirst_sagittal,
373     prone_feetfirst_coronal
374   };
375
376 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
377 {
378   h.GetImageOrientationPatient(image_orientation_patient);
379 }
380
381 #if 0
382 //
383 // this is all completely cribbed from the xmedcon library, since
384 // we're trying to do what it does, mostly.
385 patient_slice_orientation_type
386 GetPatSliceOrient(gdcm::File &h)
387 {
388   F32 image_orientation_patient[6];
389
390   // protected, do it the hard way
391   //  h.GetImageOrientationPatient(image_orientation_patient);
392   GetImageOrientationPatient(h,image_orientation_patient);
393
394   enum { headfirst, feetfirst } patient_orientation;
395   enum { supine, prone } patient_rotation;
396   enum { transaxial, sagittal, coronal } slice_orientation;
397
398   std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
399   if(patient_position == "gdcm::Unfound")
400     {
401     patient_position = "HF";
402     }
403   if(patient_position.find("HF") != std::string::npos)
404     {
405     patient_orientation = headfirst;
406     }
407   else if(patient_position.find("FF") != std::string::npos)
408     {
409     patient_orientation = feetfirst;
410     }
411
412   if(patient_position.find("S") != std::string::npos)
413     {
414     patient_rotation = supine;
415     }
416   else if(patient_position.find("P") != std::string::npos)
417     {
418     patient_rotation = prone;
419     }
420
421   if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
422      (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
423     {
424     slice_orientation = transaxial;
425     }
426   else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
427           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
428     {
429     slice_orientation = sagittal;
430     }
431   else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
432           (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
433     {
434     slice_orientation = coronal;
435     }
436   //
437   // combine results
438   patient_slice_orientation_type patient_slice_orientation = 
439     patient_slice_orientation_unknown;
440   switch (patient_rotation) 
441     {
442     case supine:
443       switch (patient_orientation) 
444         {
445         case headfirst:
446           switch (slice_orientation) 
447             {
448             case transaxial:
449               patient_slice_orientation = supine_headfirst_transaxial;
450               break;
451             case sagittal:
452               patient_slice_orientation = supine_headfirst_sagittal;
453               break;
454             case coronal:
455               patient_slice_orientation = supine_headfirst_coronal;
456               break;
457             }
458           break;
459         case feetfirst:
460           switch (slice_orientation) 
461             {
462             case transaxial:
463               patient_slice_orientation = supine_feetfirst_transaxial;
464               break;
465             case sagittal:
466               patient_slice_orientation = supine_feetfirst_sagittal;
467               break;
468             case coronal:
469               patient_slice_orientation = supine_feetfirst_coronal;
470               break;
471             }
472           break;
473         }
474       break;
475     case prone:
476       switch (patient_orientation) 
477         {
478         case headfirst:
479           switch (slice_orientation) 
480             {
481             case transaxial:
482               patient_slice_orientation = prone_headfirst_transaxial;
483               break;
484             case sagittal:
485               patient_slice_orientation = prone_headfirst_sagittal;
486               break;
487             case coronal:
488               patient_slice_orientation = prone_headfirst_coronal;
489               break;
490             }
491           break;
492         case feetfirst:
493           switch (slice_orientation) 
494             {
495             case transaxial:
496               patient_slice_orientation = prone_feetfirst_transaxial;
497               break;
498             case sagittal:
499               patient_slice_orientation = prone_feetfirst_sagittal;
500               break;
501             case coronal:
502               patient_slice_orientation = prone_feetfirst_coronal;
503               break;
504             }
505           break;
506         }
507       break;
508     }
509   if(patient_slice_orientation != patient_slice_orientation_unknown)
510     return patient_slice_orientation;
511   //
512   // this is what xmedcon does
513   if ((image_orientation_patient[0] == +1)   &&
514       (image_orientation_patient[4] == +1))   
515     patient_slice_orientation = supine_headfirst_transaxial;
516   else if ((image_orientation_patient[0] == -1)   &&
517            (image_orientation_patient[4] == +1))   
518     patient_slice_orientation = supine_feetfirst_transaxial;
519   else if ((image_orientation_patient[0] == -1)   &&
520            (image_orientation_patient[4] == -1))   
521     patient_slice_orientation = prone_headfirst_transaxial;
522   else if ((image_orientation_patient[0] == +1)   &&
523            (image_orientation_patient[4] == -1))   
524     patient_slice_orientation = prone_feetfirst_transaxial;
525
526   else if ((image_orientation_patient[1] == +1)   &&
527            (image_orientation_patient[5] == -1))   
528     patient_slice_orientation = supine_headfirst_sagittal;
529   else if ((image_orientation_patient[1] == +1)   &&
530            (image_orientation_patient[5] == +1))   
531     patient_slice_orientation = supine_feetfirst_sagittal;
532   else if ((image_orientation_patient[1] == -1)   &&
533            (image_orientation_patient[5] == -1))   
534     patient_slice_orientation = prone_headfirst_sagittal;
535   else if ((image_orientation_patient[1] == -1)   &&
536            (image_orientation_patient[5] == +1))   
537     patient_slice_orientation = prone_feetfirst_sagittal;
538
539   else if ((image_orientation_patient[0] == +1)   &&
540            (image_orientation_patient[5] == -1))   
541     patient_slice_orientation = supine_headfirst_coronal;
542   else if ((image_orientation_patient[0] == -1)   &&
543            (image_orientation_patient[5] == +1))   
544     patient_slice_orientation = supine_feetfirst_coronal;
545   else if ((image_orientation_patient[0] == -1)   &&
546            (image_orientation_patient[5] == -1))   
547     patient_slice_orientation = prone_headfirst_coronal;
548   else if ((image_orientation_patient[0] == +1)   &&
549            (image_orientation_patient[5] == +1))   
550     patient_slice_orientation = prone_feetfirst_coronal;
551   return patient_slice_orientation;
552 }
553 #else
554
555 -------------------------------------------------------------------------*/
556
557 } // end namespace gdcm