00001 /*========================================================================= 00002 00003 Program: gdcm 00004 Module: $RCSfile: gdcmOrientation.cxx,v $ 00005 Language: C++ 00006 Date: $Date: 2007/05/23 14:18:10 $ 00007 Version: $Revision: 1.25 $ 00008 00009 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de 00010 l'Image). All rights reserved. See Doc/License.txt or 00011 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details. 00012 00013 This software is distributed WITHOUT ANY WARRANTY; without even 00014 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00015 PURPOSE. See the above copyright notices for more information. 00016 00017 =========================================================================*/ 00018 00019 #include "gdcmOrientation.h" 00020 #include "gdcmFile.h" 00021 #include "gdcmDebug.h" 00022 #include <math.h> // for sqrt 00023 00024 namespace GDCM_NAME_SPACE 00025 { 00026 //-------------------------------------------------------------------- 00027 // THERALYS Algorithm to determine the most similar basic orientation 00028 // 00029 // Transliterated from original Python code. 00030 // Kept as close as possible to the original code 00031 // in order to speed up any further modif of Python code :-( 00032 //----------------------------------------------------------------------- 00033 00055 static const char *OrientationTypeStrings[] = { 00056 "Not Applicable", 00057 "Axial", 00058 "Coronal", 00059 "Sagital", 00060 "Heart Axial", 00061 "Heart Coronal", 00062 "Heart Sagital", 00063 "Axial invert", 00064 "Coronal invert", 00065 "Sagital invert", 00066 "Heart Axial invert", 00067 "Heart Coronal invert", 00068 "Heart Sagital invert", 00069 NULL 00070 }; 00071 00074 const char* Orientation::GetOrientationTypeString(OrientationType const o) 00075 { 00076 int k = (int)o; 00077 if (k < 0) 00078 k = -k + 6; 00079 00080 return OrientationTypeStrings[k]; 00081 } 00082 00085 OrientationType Orientation::GetOrientationType( File *f ) 00086 { 00087 float iop[6]; 00088 bool succ = f->GetImageOrientationPatient( iop ); 00089 if ( !succ ) 00090 { 00091 gdcmWarningMacro( "No Image Orientation (0020,0037)/(0020,0032) found in the file, cannot proceed." ) 00092 return NotApplicable; 00093 } 00094 vector3D ori1; 00095 vector3D ori2; 00096 00097 ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2]; 00098 ori2.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5]; 00099 00100 // two perpendicular vectors describe one plane 00101 double dicPlane[6][2][3] = 00102 { { { 1, 0, 0 },{ 0, 1, 0 } }, // Axial 00103 { { 1, 0, 0 },{ 0, 0, -1 } }, // Coronal 00104 { { 0, 1, 0 },{ 0, 0, -1 } }, // Sagittal 00105 { { 0.8, 0.5, 0.0 },{-0.1, 0.1 , -0.95 } }, // Axial - HEART 00106 { { 0.8, 0.5, 0.0 },{-0.6674, 0.687, 0.1794} }, // Coronal - HEART 00107 { {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794} } // Sagittal - HEART 00108 }; 00109 00110 vector3D refA; 00111 vector3D refB; 00112 int i = 0; 00113 Res res; // [ <result> , <memory of the last succes calcule> ] 00114 res.first = 0; 00115 res.second = 99999; 00116 00117 for (int numDicPlane=0; numDicPlane<6; numDicPlane++) 00118 { 00119 ++i; 00120 // refA=plane[0] 00121 refA.x = dicPlane[numDicPlane][0][0]; 00122 refA.y = dicPlane[numDicPlane][0][1]; 00123 refA.z = dicPlane[numDicPlane][0][2]; 00124 // refB=plane[1] 00125 refB.x = dicPlane[numDicPlane][1][0]; 00126 refB.y = dicPlane[numDicPlane][1][1]; 00127 refB.z = dicPlane[numDicPlane][1][2]; 00128 res=VerfCriterion( i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res ); 00129 res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res ); 00130 } 00131 // res thought looks like is a float value, but is indeed an int 00132 // try casting it to int first then enum value to please VS7: 00133 int int_res = (int)res.first; 00134 gdcmAssertMacro( int_res <= 6 && int_res >= -6); 00135 return (OrientationType)int_res; 00136 } 00137 00138 Res 00139 Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const &in) 00140 { 00141 Res res; 00142 double type = in.first; 00143 double criterion = in.second; 00144 if (/*criterionNew < 0.1 && */criterionNew < criterion) 00145 { 00146 type = typeCriterion; 00147 criterion = criterionNew; 00148 } 00149 res.first = type; 00150 res.second = criterion; 00151 return res; 00152 } 00153 00154 inline double square_dist(vector3D const &v1, vector3D const &v2) 00155 { 00156 double res; 00157 res = (v1.x - v2.x)*(v1.x - v2.x) + 00158 (v1.y - v2.y)*(v1.y - v2.y) + 00159 (v1.z - v2.z)*(v1.z - v2.z); 00160 return res; 00161 } 00162 00163 //------------------------- Purpose : ----------------------------------- 00164 //- This function determines the orientation similarity of two planes. 00165 // Each plane is described by two vectors. 00166 //------------------------- Parameters : -------------------------------- 00167 //- <refA> : - type : vector 3D (double) 00168 //- <refB> : - type : vector 3D (double) 00169 // - Description of the first plane 00170 //- <ori1> : - type : vector 3D (double) 00171 //- <ori2> : - type : vector 3D (double) 00172 // - Description of the second plane 00173 //------------------------- Return : ------------------------------------ 00174 // double : 0 if the planes are perpendicular. While the difference of 00175 // the orientation between the planes are big more enlarge is 00176 // the criterion. 00177 //------------------------- Other : ------------------------------------- 00178 // The calculus is based with vectors normalice 00179 double 00180 Orientation::CalculLikelyhood2Vec(vector3D const &refA, vector3D const &refB, 00181 vector3D const &ori1, vector3D const &ori2 ) 00182 { 00183 00184 vector3D ori3 = ProductVectorial(ori1,ori2); 00185 vector3D refC = ProductVectorial(refA,refB); 00186 double res = square_dist(refC, ori3); 00187 00188 return sqrt(res); 00189 } 00190 00191 //------------------------- Purpose : ----------------------------------- 00192 //- Calculus of the poduct vectorial between two vectors 3D 00193 //------------------------- Parameters : -------------------------------- 00194 //- <vec1> : - type : vector 3D (double) 00195 //- <vec2> : - type : vector 3D (double) 00196 //------------------------- Return : ------------------------------------ 00197 // (vec) : - Vector 3D 00198 //------------------------- Other : ------------------------------------- 00199 vector3D 00200 Orientation::ProductVectorial(vector3D const &vec1, vector3D const &vec2) 00201 { 00202 vector3D vec3; 00203 vec3.x = vec1.y*vec2.z - vec1.z*vec2.y; 00204 vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x); 00205 vec3.z = vec1.x*vec2.y - vec1.y*vec2.x; 00206 00207 return vec3; 00208 } 00209 00210 00211 // ------------------------------------------------------------------------ 00212 /* 00213 2.2.2 Orientation of DICOM images 00214 00215 00216 http://www.dclunie.com/medical-image-faq/html/part2.html#DICOMOrientation 00217 says : 00218 00219 A question that is frequently asked in comp.protocols.dicom is how to determine 00220 which side of an image is which (e.g. left, right) and so on. 00221 The short answer is that for projection radiographs this is specified 00222 explicitly using the "Patient Orientation" attribute, and for cross-sectional 00223 images it needs to be derived from the "Image Orientation (Patient)" direction 00224 cosines. In the standard these are explained as follows: 00225 00226 * "C.7.6.1.1.1 Patient Orientation. 00227 The Patient Orientation (0020,0020) relative to the image 00228 plane shall be specified by two values that designate the 00229 anatomical direction of the positive row axis (left to right) 00230 and the positive column axis (top to bottom). 00231 The first entry is the direction of the rows, given by the 00232 direction of the last pixel in the first row from the first 00233 pixel in that row. 00234 The second entry is the direction of the columns, given by 00235 the direction of the last pixel in the first column from the 00236 first pixel in that column. 00237 Anatomical direction shall be designated by the capital 00238 letters: A (anterior), P (posterior), R (right),L (left), 00239 H (head), F (foot). 00240 Each value of the orientation attribute shall contain at 00241 least one of these characters. 00242 If refinements in the orientation descriptions are to be 00243 specified, then they shall be designated by one or two 00244 additional letters in each value. 00245 Within each value, the letters shall be ordered with the 00246 principal orientation designated in the first character." 00247 00248 * "C.7.6.2.1.1 Image Position And Image Orientation. 00249 The "Image Position (Patient)" (0020,0032) specifies the x, y, 00250 and z coordinates of the upper left hand corner of the image; 00251 it is the center of the first voxel transmitted. 00252 The "Image Orientation (Patient)" (0020,0037) specifies the 00253 direction cosines of the first row and the first column 00254 with respect to the patient. 00255 These Attributes shall be provided as a pair. 00256 Row value for the x, y, and z axes respectively followed by 00257 the Column value for the x, y, and z axes respectively. 00258 The direction of the axes is defined fully by the patient's 00259 orientation. 00260 The x-axis is increasing to the left hand side of the patient. 00261 The y-axis is increasing to the posterior side of the patient. 00262 The z-axis is increasing toward the head of the patient. 00263 The patient based coordinate system is a right handed system, 00264 i.e. the vector cross product of a unit vector along the 00265 positive x-axis and a unit vector along the positive y-axis 00266 is equal to a unit vector along the positive z-axis." 00267 00268 Some simple code to take one of the direction cosines (vectors) from the 00269 Image Orientation (Patient) attribute and generate strings equivalent to one 00270 of the values of Patient Orientation looks like this (noting that if the vector 00271 is not aligned exactly with one of the major axes, the resulting string will 00272 have multiple letters in as described under "refinements" in C.7.6.1.1.1): 00273 00274 */ 00275 00294 std::string Orientation::GetOrientation ( File *f ) 00295 { 00296 float iop[6]; 00297 if ( !f->GetImageOrientationPatient( iop ) ) 00298 return GDCM_UNFOUND; 00299 00300 std::string orientation; 00301 orientation = GetSingleOrientation ( iop ) 00302 + "\\" 00303 + GetSingleOrientation ( iop + 3 ); 00304 return orientation; 00305 } 00306 00307 00308 std::string Orientation::GetSingleOrientation ( float *iop) 00309 { 00310 std::string orientation; 00311 00312 char orientationX = iop[0] < 0 ? 'R' : 'L'; 00313 char orientationY = iop[1] < 0 ? 'A' : 'P'; 00314 char orientationZ = iop[2] < 0 ? 'F' : 'H'; 00315 00316 double absX = iop[0]; 00317 if (absX < 0) absX = -absX; 00318 double absY = iop[1]; 00319 if (absY < 0) absY = -absY; 00320 double absZ = iop[2]; 00321 if (absZ < 0) absZ = -absZ; 00322 00323 for (int i=0; i<3; ++i) 00324 { 00325 if (absX>.0001 && absX>absY && absX>absZ) 00326 { 00327 orientation = orientation + orientationX; 00328 absX=0; 00329 } 00330 else if (absY>.0001 && absY>absX && absY>absZ) 00331 { 00332 orientation = orientation + orientationY; 00333 absY=0; 00334 } 00335 else if (absZ>.0001 && absZ>absX && absZ>absY) 00336 { 00337 orientation = orientation + orientationZ; 00338 absZ=0; 00339 } 00340 else 00341 break; 00342 } 00343 return orientation; 00344 } 00345 00346 00347 /*------------------------------------------------------------------- 00348 00349 Some more stuff, from XMedcon 00350 00351 ---> Personal remark from JPRx : 00352 --> patient_position (0x0018,0x5100) can be "HFS ", "FFS ", "HFP ", "FFP " 00353 --> or, not so common, 00354 // HFDR = Head First-Decubitus Right 00355 // HFDL = Head First-Decubitus Left 00356 // FFDR = Feet First-Decubitus Right 00357 // FFDL = Feet First-Decubitus Left 00358 --> the cosines may have any value -1.< <+1., for MR images ! 00359 00360 enum patient_slice_orientation_type 00361 { 00362 patient_slice_orientation_unknown = 0, 00363 supine_headfirst_transaxial, 00364 supine_headfirst_sagittal, 00365 supine_headfirst_coronal, 00366 supine_feetfirst_transaxial, 00367 supine_feetfirst_sagittal, 00368 supine_feetfirst_coronal, 00369 prone_headfirst_transaxial, 00370 prone_headfirst_sagittal, 00371 prone_headfirst_coronal, 00372 prone_feetfirst_transaxial, 00373 prone_feetfirst_sagittal, 00374 prone_feetfirst_coronal 00375 }; 00376 00377 void GetImageOrientationPatient(gdcm::File &h,F32 image_orientation_patient[6]) 00378 { 00379 h.GetImageOrientationPatient(image_orientation_patient); 00380 } 00381 00382 #if 0 00383 // 00384 // this is all completely cribbed from the xmedcon library, since 00385 // we're trying to do what it does, mostly. 00386 patient_slice_orientation_type 00387 GetPatSliceOrient(gdcm::File &h) 00388 { 00389 F32 image_orientation_patient[6]; 00390 00391 // protected, do it the hard way 00392 // h.GetImageOrientationPatient(image_orientation_patient); 00393 GetImageOrientationPatient(h,image_orientation_patient); 00394 00395 enum { headfirst, feetfirst } patient_orientation; 00396 enum { supine, prone } patient_rotation; 00397 enum { transaxial, sagittal, coronal } slice_orientation; 00398 00399 std::string patient_position = h.GetEntryByNumber(0x0018,0x5100); 00400 if(patient_position == "gdcm::Unfound") 00401 { 00402 patient_position = "HF"; 00403 } 00404 if(patient_position.find("HF") != std::string::npos) 00405 { 00406 patient_orientation = headfirst; 00407 } 00408 else if(patient_position.find("FF") != std::string::npos) 00409 { 00410 patient_orientation = feetfirst; 00411 } 00412 00413 if(patient_position.find("S") != std::string::npos) 00414 { 00415 patient_rotation = supine; 00416 } 00417 else if(patient_position.find("P") != std::string::npos) 00418 { 00419 patient_rotation = prone; 00420 } 00421 00422 if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) && 00423 (image_orientation_patient[4] == +1 || image_orientation_patient[4] == -1)) 00424 { 00425 slice_orientation = transaxial; 00426 } 00427 else if((image_orientation_patient[1] == 1 || image_orientation_patient[1] == -1) && 00428 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1)) 00429 { 00430 slice_orientation = sagittal; 00431 } 00432 else if((image_orientation_patient[0] == 1 || image_orientation_patient[0] == -1) && 00433 (image_orientation_patient[5] == +1 || image_orientation_patient[5] == -1)) 00434 { 00435 slice_orientation = coronal; 00436 } 00437 // 00438 // combine results 00439 patient_slice_orientation_type patient_slice_orientation = 00440 patient_slice_orientation_unknown; 00441 switch (patient_rotation) 00442 { 00443 case supine: 00444 switch (patient_orientation) 00445 { 00446 case headfirst: 00447 switch (slice_orientation) 00448 { 00449 case transaxial: 00450 patient_slice_orientation = supine_headfirst_transaxial; 00451 break; 00452 case sagittal: 00453 patient_slice_orientation = supine_headfirst_sagittal; 00454 break; 00455 case coronal: 00456 patient_slice_orientation = supine_headfirst_coronal; 00457 break; 00458 } 00459 break; 00460 case feetfirst: 00461 switch (slice_orientation) 00462 { 00463 case transaxial: 00464 patient_slice_orientation = supine_feetfirst_transaxial; 00465 break; 00466 case sagittal: 00467 patient_slice_orientation = supine_feetfirst_sagittal; 00468 break; 00469 case coronal: 00470 patient_slice_orientation = supine_feetfirst_coronal; 00471 break; 00472 } 00473 break; 00474 } 00475 break; 00476 case prone: 00477 switch (patient_orientation) 00478 { 00479 case headfirst: 00480 switch (slice_orientation) 00481 { 00482 case transaxial: 00483 patient_slice_orientation = prone_headfirst_transaxial; 00484 break; 00485 case sagittal: 00486 patient_slice_orientation = prone_headfirst_sagittal; 00487 break; 00488 case coronal: 00489 patient_slice_orientation = prone_headfirst_coronal; 00490 break; 00491 } 00492 break; 00493 case feetfirst: 00494 switch (slice_orientation) 00495 { 00496 case transaxial: 00497 patient_slice_orientation = prone_feetfirst_transaxial; 00498 break; 00499 case sagittal: 00500 patient_slice_orientation = prone_feetfirst_sagittal; 00501 break; 00502 case coronal: 00503 patient_slice_orientation = prone_feetfirst_coronal; 00504 break; 00505 } 00506 break; 00507 } 00508 break; 00509 } 00510 if(patient_slice_orientation != patient_slice_orientation_unknown) 00511 return patient_slice_orientation; 00512 // 00513 // this is what xmedcon does 00514 if ((image_orientation_patient[0] == +1) && 00515 (image_orientation_patient[4] == +1)) 00516 patient_slice_orientation = supine_headfirst_transaxial; 00517 else if ((image_orientation_patient[0] == -1) && 00518 (image_orientation_patient[4] == +1)) 00519 patient_slice_orientation = supine_feetfirst_transaxial; 00520 else if ((image_orientation_patient[0] == -1) && 00521 (image_orientation_patient[4] == -1)) 00522 patient_slice_orientation = prone_headfirst_transaxial; 00523 else if ((image_orientation_patient[0] == +1) && 00524 (image_orientation_patient[4] == -1)) 00525 patient_slice_orientation = prone_feetfirst_transaxial; 00526 00527 else if ((image_orientation_patient[1] == +1) && 00528 (image_orientation_patient[5] == -1)) 00529 patient_slice_orientation = supine_headfirst_sagittal; 00530 else if ((image_orientation_patient[1] == +1) && 00531 (image_orientation_patient[5] == +1)) 00532 patient_slice_orientation = supine_feetfirst_sagittal; 00533 else if ((image_orientation_patient[1] == -1) && 00534 (image_orientation_patient[5] == -1)) 00535 patient_slice_orientation = prone_headfirst_sagittal; 00536 else if ((image_orientation_patient[1] == -1) && 00537 (image_orientation_patient[5] == +1)) 00538 patient_slice_orientation = prone_feetfirst_sagittal; 00539 00540 else if ((image_orientation_patient[0] == +1) && 00541 (image_orientation_patient[5] == -1)) 00542 patient_slice_orientation = supine_headfirst_coronal; 00543 else if ((image_orientation_patient[0] == -1) && 00544 (image_orientation_patient[5] == +1)) 00545 patient_slice_orientation = supine_feetfirst_coronal; 00546 else if ((image_orientation_patient[0] == -1) && 00547 (image_orientation_patient[5] == -1)) 00548 patient_slice_orientation = prone_headfirst_coronal; 00549 else if ((image_orientation_patient[0] == +1) && 00550 (image_orientation_patient[5] == +1)) 00551 patient_slice_orientation = prone_feetfirst_coronal; 00552 return patient_slice_orientation; 00553 } 00554 #else 00555 00556 -------------------------------------------------------------------------*/ 00557 00558 } // end namespace gdcm