1 /*=========================================================================
4 Module: $RCSfile: gdcmOrientation.cxx,v $
6 Date: $Date: 2005/10/03 14:54:16 $
7 Version: $Revision: 1.16 $
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.
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.
17 =========================================================================*/
19 #include "gdcmOrientation.h"
21 #include "gdcmDebug.h"
22 #include <math.h> // for sqrt
26 //--------------------------------------------------------------------
27 // THERALYS Algorithm to determine the most similar basic orientation
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 //-----------------------------------------------------------------------
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)
44 * # -2 : Coronal invert
46 * # -3 : Sagital invert
48 * # -4 : Heart Axial invert
50 * # -5 : Heart Coronal invert
52 * # -6 : Heart Sagital invert
55 static const char *OrientationTypeStrings[] = {
67 "Heart Coronal invert",
68 "Heart Sagital invert",
72 const char* Orientation::GetOrientationTypeString(OrientationType const o)
78 return OrientationTypeStrings[k];
81 OrientationType Orientation::GetOrientationType( File *f )
84 bool succ = f->GetImageOrientationPatient( iop );
87 gdcmErrorMacro( "No Image Orientation (0020,0037)/(0020,0032) found in the file, cannot proceed." )
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];
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
109 Res res; // [ <result> , <memory of the last succes calcule> ]
113 for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
117 refA.x = dicPlane[numDicPlane][0][0];
118 refA.y = dicPlane[numDicPlane][0][1];
119 refA.z = dicPlane[numDicPlane][0][2];
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 );
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;
135 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
138 double type = in.first;
139 double criterion = in.second;
140 if (/*criterionNew < 0.1 && */criterionNew < criterion)
142 type = typeCriterion;
143 criterion = criterionNew;
146 res.second = criterion;
150 inline double square_dist(vector3D const &v1, vector3D const &v2)
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);
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
173 //------------------------- Other : -------------------------------------
174 // The calculus is based with vectors normalice
176 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB,
177 vector3D const &ori1, vector3D const &ori2 )
180 vector3D ori3 = ProductVectorial(ori1,ori2);
181 vector3D refC = ProductVectorial(refA,refB);
182 double res = square_dist(refC, ori3);
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 : -------------------------------------
196 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2)
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;
208 // ---------------------------------------------------------------------------
209 // Here is the original Python code, kindly supplied by THERALYS
211 // C++ code doesn't give good results
216 def TypeOrientation(self,file0):
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 : ------------------------------------
228 # -2 : Coronal invert
230 # -3 : Sagital invert
232 # -4 : Heart Axial invert
234 # -5 : Heart Coronal invert
236 # -6 : Heart Sagital invert
238 # ------------------------- Other : -------------------------------------
239 # This method finds the most similar basic orientation.
242 toRead = gdcm.File(file0)
243 ValDict = GetValuesDict(toRead)
245 imageOrientation=ValDict["Image Orientation (Patient)"]
247 imageOrientation=ValDict["Image Orientation"]
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])]
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
265 res=[0,99999] ## [ <result> , <memory of the last succes calcule> ]
266 for plane in dicPlane:
270 res=self.VerfCriterion( i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
271 res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
278 def VerfCriterion(self,typeCriterion,criterionNew,res):
281 # if criterionNew<0.1 and criterionNew<criterion:
282 if criterionNew<criterion:
283 criterion=criterionNew
285 return [ type , criterion ]
288 def CalculLikelyhood2Vec(self,refA,refB,ori1,ori2):
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
305 # ------------------------- Other : -------------------------------------
306 # The calculus is based with vectors normalice
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)
314 def ProductVectorial(self,vec1,vec2):
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 : -------------------------------------
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]
331 def GetValuesDict(image):
333 Returns a dictionnary containing values associated with Field Names
334 dict["Dicom Field Name"]="Dicom field value"
336 val=image.GetFirstEntry()
339 if isinstance(val,gdcm.ValEntryPtr):
340 dic[val.GetName()]=val.GetValue()
341 val=image.GetNextEntry()
347 // ------------------------------------------------------------------------
349 2.2.2 Orientation of DICOM images
352 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation
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:
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
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),
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."
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
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."
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):
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
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"
429 std::string Orientation::GetOrientation ( File *f )
432 if ( !f->GetImageOrientationPatient( iop ) )
435 std::string orientation;
436 orientation = GetSingleOrientation ( iop )
438 + GetSingleOrientation ( iop + 3 );
443 std::string Orientation::GetSingleOrientation ( float *iop)
445 std::string orientation;
447 char orientationX = iop[0] < 0 ? 'R' : 'L';
448 char orientationY = iop[1] < 0 ? 'A' : 'P';
449 char orientationZ = iop[2] < 0 ? 'F' : 'H';
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;
458 for (int i=0; i<3; ++i)
460 if (absX>.0001 && absX>absY && absX>absZ)
462 orientation = orientation + orientationX;
465 else if (absY>.0001 && absY>absX && absY>absZ)
467 orientation = orientation + orientationY;
470 else if (absZ>.0001 && absZ>absX && absZ>absY)
472 orientation = orientation + orientationZ;
482 /*-------------------------------------------------------------------
484 Some more stuff, from XMedcon
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 !
490 enum patient_slice_orientation_type
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
507 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
509 h.GetImageOrientationPatient(image_orientation_patient);
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)
519 F32 image_orientation_patient[6];
521 // protected, do it the hard way
522 // h.GetImageOrientationPatient(image_orientation_patient);
523 GetImageOrientationPatient(h,image_orientation_patient);
525 enum { headfirst, feetfirst } patient_orientation;
526 enum { supine, prone } patient_rotation;
527 enum { transaxial, sagittal, coronal } slice_orientation;
529 std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
530 if(patient_position == "gdcm::Unfound")
532 patient_position = "HF";
534 if(patient_position.find("HF") != std::string::npos)
536 patient_orientation = headfirst;
538 else if(patient_position.find("FF") != std::string::npos)
540 patient_orientation = feetfirst;
543 if(patient_position.find("S") != std::string::npos)
545 patient_rotation = supine;
547 else if(patient_position.find("P") != std::string::npos)
549 patient_rotation = prone;
552 if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
553 (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
555 slice_orientation = transaxial;
557 else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
558 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
560 slice_orientation = sagittal;
562 else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
563 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
565 slice_orientation = coronal;
569 patient_slice_orientation_type patient_slice_orientation =
570 patient_slice_orientation_unknown;
571 switch (patient_rotation)
574 switch (patient_orientation)
577 switch (slice_orientation)
580 patient_slice_orientation = supine_headfirst_transaxial;
583 patient_slice_orientation = supine_headfirst_sagittal;
586 patient_slice_orientation = supine_headfirst_coronal;
591 switch (slice_orientation)
594 patient_slice_orientation = supine_feetfirst_transaxial;
597 patient_slice_orientation = supine_feetfirst_sagittal;
600 patient_slice_orientation = supine_feetfirst_coronal;
607 switch (patient_orientation)
610 switch (slice_orientation)
613 patient_slice_orientation = prone_headfirst_transaxial;
616 patient_slice_orientation = prone_headfirst_sagittal;
619 patient_slice_orientation = prone_headfirst_coronal;
624 switch (slice_orientation)
627 patient_slice_orientation = prone_feetfirst_transaxial;
630 patient_slice_orientation = prone_feetfirst_sagittal;
633 patient_slice_orientation = prone_feetfirst_coronal;
640 if(patient_slice_orientation != patient_slice_orientation_unknown)
641 return patient_slice_orientation;
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;
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;
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;
686 -------------------------------------------------------------------------*/
688 } // end namespace gdcm