1 /*=========================================================================
4 Module: $RCSfile: gdcmOrientation.cxx,v $
6 Date: $Date: 2005/09/22 14:41:25 $
7 Version: $Revision: 1.12 $
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
54 double Orientation::TypeOrientation( File *f )
57 bool succ = f->GetImageOrientationPatient( iop );
60 gdcmErrorMacro( "No Image Orientation (0020,0037) found in the file, cannot proceed." )
64 std::cout << " iop : ";
66 std::cout << iop[i] << " ";
67 std::cout << std::endl;
72 ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2];
73 ori2.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5];
75 // two perpendicular vectors describe one plane
76 double dicPlane[6][2][3] =
77 { { { 1, 0, 0 },{ 0, 1, 0 } }, // Axial
78 { { 1, 0, 0 },{ 0, 0, -1 } }, // Coronal
79 { { 0, 1, 0 },{ 0, 0, -1 } }, // Sagittal
80 { { 0.8, 0.5, 0.0 },{-0.1, 0.1 , -0.95 } }, // Axial - HEART
81 { { 0.8, 0.5, 0.0 },{-0.6674, 0.687, 0.1794} }, // Coronal - HEART
82 { {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794} } // Sagittal - HEART
88 Res res; // [ <result> , <memory of the last succes calcule> ]
92 std::cout << "-------------- res : " << res.first << "|" << res.second
95 for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
99 refA.x = dicPlane[numDicPlane][0][0];
100 refA.y = dicPlane[numDicPlane][0][1];
101 refA.z = dicPlane[numDicPlane][0][2];
103 refB.x = dicPlane[numDicPlane][1][0];
104 refB.y = dicPlane[numDicPlane][1][1];
105 refB.z = dicPlane[numDicPlane][1][2];
106 res=VerfCriterion( i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
107 std::cout << "-------------- res : " << res.first << "|" << res.second
109 res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
110 std::cout << "-------------- res : " << res.first << "|" << res.second
116 // res=[0,99999] ## [ <result> , <memory of the last succes calculus> ]
117 // for plane in dicPlane:
121 // res=self.VerfCriterion( i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
122 // res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
128 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
131 // double type = in.first;
132 double criterion = in.second;
133 if (/*criterionNew < 0.1 && */criterionNew < criterion)
135 res.first = typeCriterion;
136 res.second = criterionNew;
140 // criterion = res[1]
141 // # if criterionNew<0.1 and criterionNew<criterion:
142 // if criterionNew<criterion:
143 // criterion=criterionNew
144 // type=typeCriterion
145 // return [ type , 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;
483 } // end namespace gdcm