1 /*=========================================================================
4 Module: $RCSfile: gdcmOrientation.cxx,v $
6 Date: $Date: 2008/08/18 12:27:10 $
7 Version: $Revision: 1.27 $
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
24 namespace GDCM_NAME_SPACE
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, Sagittal) 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 : Sagittal invert
48 * # -4 : Heart Axial invert
50 * # -5 : Heart Coronal invert
51 * # 6 : Heart Sagittal
52 * # -6 : Heart Sagittal invert
55 static const char *OrientationTypeStrings[] = {
67 "Heart Coronal invert",
68 "Heart Sagittal invert",
72 /// \brief returns human readable interpretation of the most
73 /// similar basic orientation (Axial, Coronal, Sagittal, ...) of the image
74 const char* Orientation::GetOrientationTypeString(OrientationType const o)
80 return OrientationTypeStrings[k];
83 /// \brief returns of the most similar basic orientation
84 /// (Axial, Coronal, Sagittal, ...) of the image
85 OrientationType Orientation::GetOrientationType( File *f )
88 bool succ = f->GetImageOrientationPatient( iop );
91 gdcmWarningMacro( "No Image Orientation (0020,0037)/(0020,0032) found in the file, cannot proceed." )
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];
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
113 Res res; // [ <result> , <memory of the last succes calcule> ]
117 for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
121 refA.x = dicPlane[numDicPlane][0][0];
122 refA.y = dicPlane[numDicPlane][0][1];
123 refA.z = dicPlane[numDicPlane][0][2];
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 );
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;
139 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in)
142 double type = in.first;
143 double criterion = in.second;
144 if (/*criterionNew < 0.1 && */criterionNew < criterion)
146 type = typeCriterion;
147 criterion = criterionNew;
150 res.second = criterion;
154 inline double square_dist(vector3D const &v1, vector3D const &v2)
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);
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
177 //------------------------- Other : -------------------------------------
178 // The calculus is based with vectors normalice
180 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB,
181 vector3D const &ori1, vector3D const &ori2 )
184 vector3D ori3 = ProductVectorial(ori1,ori2);
185 vector3D refC = ProductVectorial(refA,refB);
186 double res = square_dist(refC, ori3);
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 : -------------------------------------
200 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2)
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;
211 // ------------------------------------------------------------------------
213 2.2.2 Orientation of DICOM images
216 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation
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:
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
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),
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."
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
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."
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):
277 * \brief Computes the Patient Orientation relative to the image plane
278 * from the 'Image Orientation (Patient)'
279 * - or from 0020 0035Image Orientation (RET) -
280 * - The first entry is the direction of the rows, given by the
281 * direction of the last pixel in the first row from the first
283 * - The second entry is the direction of the columns, given by
284 * the direction of the last pixel in the first column from the
285 * first pixel in that column.
286 * Anatomical direction is designated by the capital
287 * letters: A (anterior), P (posterior), R (right),L (left),
288 * H (head), F (foot).
289 * - Refinements in the orientation descriptions are designated
290 * by one or two additional letters in each value.
291 * Use it when "Patient Orientation" (0020,0020) is not found
292 * @return orientation string as "rowsOrientation\columnsOrientation"
294 std::string Orientation::GetOrientation ( File *f )
297 if ( !f->GetImageOrientationPatient( iop ) )
300 std::string orientation;
301 orientation = GetSingleOrientation ( iop )
303 + GetSingleOrientation ( iop + 3 );
308 std::string Orientation::GetSingleOrientation ( float *iop)
310 std::string orientation;
312 char orientationX = iop[0] < 0 ? 'R' : 'L';
313 char orientationY = iop[1] < 0 ? 'A' : 'P';
314 char orientationZ = iop[2] < 0 ? 'F' : 'H';
316 double absX = iop[0];
317 if (absX < 0) absX = -absX;
318 double absY = iop[1];
319 if (absY < 0) absY = -absY;
320 double absZ = iop[2];
321 if (absZ < 0) absZ = -absZ;
323 for (int i=0; i<3; ++i)
325 if (absX>.0001 && absX>absY && absX>absZ)
327 orientation = orientation + orientationX;
330 else if (absY>.0001 && absY>absX && absY>absZ)
332 orientation = orientation + orientationY;
335 else if (absZ>.0001 && absZ>absX && absZ>absY)
337 orientation = orientation + orientationZ;
347 /*-------------------------------------------------------------------
349 Some more stuff, from XMedcon
351 ---> Personal remark from JPRx :
352 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP "
353 --> or, not so common,
354 // HFDR = Head First-Decubitus Right
355 // HFDL = Head First-Decubitus Left
356 // FFDR = Feet First-Decubitus Right
357 // FFDL = Feet First-Decubitus Left
358 --> the cosines may have any value -1.< <+1., for MR images !
360 enum patient_slice_orientation_type
362 patient_slice_orientation_unknown = 0,
363 supine_headfirst_transaxial,
364 supine_headfirst_sagittal,
365 supine_headfirst_coronal,
366 supine_feetfirst_transaxial,
367 supine_feetfirst_sagittal,
368 supine_feetfirst_coronal,
369 prone_headfirst_transaxial,
370 prone_headfirst_sagittal,
371 prone_headfirst_coronal,
372 prone_feetfirst_transaxial,
373 prone_feetfirst_sagittal,
374 prone_feetfirst_coronal
377 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6])
379 h.GetImageOrientationPatient(image_orientation_patient);
384 // this is all completely cribbed from the xmedcon library, since
385 // we're trying to do what it does, mostly.
386 patient_slice_orientation_type
387 GetPatSliceOrient(gdcm::File &h)
389 F32 image_orientation_patient[6];
391 // protected, do it the hard way
392 // h.GetImageOrientationPatient(image_orientation_patient);
393 GetImageOrientationPatient(h,image_orientation_patient);
395 enum { headfirst, feetfirst } patient_orientation;
396 enum { supine, prone } patient_rotation;
397 enum { transaxial, sagittal, coronal } slice_orientation;
399 std::string patient_position = h.GetEntryByNumber(0x0018,0x5100);
400 if(patient_position == GDCM_UNFOUND)
402 patient_position = "HF";
404 if(patient_position.find("HF") != std::string::npos)
406 patient_orientation = headfirst;
408 else if(patient_position.find("FF") != std::string::npos)
410 patient_orientation = feetfirst;
413 if(patient_position.find("S") != std::string::npos)
415 patient_rotation = supine;
417 else if(patient_position.find("P") != std::string::npos)
419 patient_rotation = prone;
422 if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
423 (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1))
425 slice_orientation = transaxial;
427 else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) &&
428 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
430 slice_orientation = sagittal;
432 else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) &&
433 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1))
435 slice_orientation = coronal;
439 patient_slice_orientation_type patient_slice_orientation =
440 patient_slice_orientation_unknown;
441 switch (patient_rotation)
444 switch (patient_orientation)
447 switch (slice_orientation)
450 patient_slice_orientation = supine_headfirst_transaxial;
453 patient_slice_orientation = supine_headfirst_sagittal;
456 patient_slice_orientation = supine_headfirst_coronal;
461 switch (slice_orientation)
464 patient_slice_orientation = supine_feetfirst_transaxial;
467 patient_slice_orientation = supine_feetfirst_sagittal;
470 patient_slice_orientation = supine_feetfirst_coronal;
477 switch (patient_orientation)
480 switch (slice_orientation)
483 patient_slice_orientation = prone_headfirst_transaxial;
486 patient_slice_orientation = prone_headfirst_sagittal;
489 patient_slice_orientation = prone_headfirst_coronal;
494 switch (slice_orientation)
497 patient_slice_orientation = prone_feetfirst_transaxial;
500 patient_slice_orientation = prone_feetfirst_sagittal;
503 patient_slice_orientation = prone_feetfirst_coronal;
510 if(patient_slice_orientation != patient_slice_orientation_unknown)
511 return patient_slice_orientation;
513 // this is what xmedcon does
514 if ((image_orientation_patient[0] == +1) &&
515 (image_orientation_patient[4] == +1))
516 patient_slice_orientation = supine_headfirst_transaxial;
517 else if ((image_orientation_patient[0] == -1) &&
518 (image_orientation_patient[4] == +1))
519 patient_slice_orientation = supine_feetfirst_transaxial;
520 else if ((image_orientation_patient[0] == -1) &&
521 (image_orientation_patient[4] == -1))
522 patient_slice_orientation = prone_headfirst_transaxial;
523 else if ((image_orientation_patient[0] == +1) &&
524 (image_orientation_patient[4] == -1))
525 patient_slice_orientation = prone_feetfirst_transaxial;
527 else if ((image_orientation_patient[1] == +1) &&
528 (image_orientation_patient[5] == -1))
529 patient_slice_orientation = supine_headfirst_sagittal;
530 else if ((image_orientation_patient[1] == +1) &&
531 (image_orientation_patient[5] == +1))
532 patient_slice_orientation = supine_feetfirst_sagittal;
533 else if ((image_orientation_patient[1] == -1) &&
534 (image_orientation_patient[5] == -1))
535 patient_slice_orientation = prone_headfirst_sagittal;
536 else if ((image_orientation_patient[1] == -1) &&
537 (image_orientation_patient[5] == +1))
538 patient_slice_orientation = prone_feetfirst_sagittal;
540 else if ((image_orientation_patient[0] == +1) &&
541 (image_orientation_patient[5] == -1))
542 patient_slice_orientation = supine_headfirst_coronal;
543 else if ((image_orientation_patient[0] == -1) &&
544 (image_orientation_patient[5] == +1))
545 patient_slice_orientation = supine_feetfirst_coronal;
546 else if ((image_orientation_patient[0] == -1) &&
547 (image_orientation_patient[5] == -1))
548 patient_slice_orientation = prone_headfirst_coronal;
549 else if ((image_orientation_patient[0] == +1) &&
550 (image_orientation_patient[5] == +1))
551 patient_slice_orientation = prone_feetfirst_coronal;
552 return patient_slice_orientation;
556 -------------------------------------------------------------------------*/
558 } // end namespace gdcm