From 69d2c65c131c9752b1d64aa9b21de29c3acfe62d Mon Sep 17 00:00:00 2001 From: malaterre Date: Sun, 24 Jul 2005 02:34:41 +0000 Subject: [PATCH] ENH: Moving Orientation stuff into its own class. gdcm::File is already too complicated... --- src/CMakeLists.txt | 1 + src/gdcmFile.cxx | 172 +---------------------------------- src/gdcmFile.h | 19 +--- src/gdcmOrientation.cxx | 192 ++++++++++++++++++++++++++++++++++++++++ src/gdcmOrientation.h | 56 ++++++++++++ 5 files changed, 253 insertions(+), 187 deletions(-) create mode 100644 src/gdcmOrientation.cxx create mode 100644 src/gdcmOrientation.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6a271f79..4898df50 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -63,6 +63,7 @@ SET(libgdcm_la_SOURCES gdcmJpeg16.cxx gdcmJpeg2000.cxx # gdcmMpeg.cxx + gdcmOrientation.cxx gdcmPixelReadConvert.cxx gdcmPixelWriteConvert.cxx gdcmRLEFrame.cxx diff --git a/src/gdcmFile.cxx b/src/gdcmFile.cxx index e21f6e0a..e1c6c35f 100644 --- a/src/gdcmFile.cxx +++ b/src/gdcmFile.cxx @@ -3,8 +3,8 @@ Program: gdcm Module: $RCSfile: gdcmFile.cxx,v $ Language: C++ - Date: $Date: 2005/07/24 02:10:48 $ - Version: $Revision: 1.262 $ + Date: $Date: 2005/07/24 02:34:41 $ + Version: $Revision: 1.263 $ Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). All rights reserved. See Doc/License.txt or @@ -78,7 +78,6 @@ #include #include //sscanf #include // for atoi -#include // for pow namespace gdcm { @@ -1788,173 +1787,6 @@ bool File::Load( std::string const &fileName ) } #endif -// ----------------------------------------------------------------------------------------- -// THERALYS Algorithm to determine the most similar basic orientation -// -// Transliterated from original Python code. -// Kept as close as possible to the original code -// in order to speed up any further modif of Python code :-( -// ------------------------------------------------------------------------------------------ - -/** - * \brief THERALYS' Algorithm to determine the most similar basic orientation - * (Axial, Coronal, Sagital) of the image - * \note Should be run on the first gdcm::File of a 'coherent' Serie - * @return orientation code - * @return orientation code - * # 0 : Not Applicable (neither 0020,0037 Image Orientation Patient - * # nor 0020,0032Image Position found ) - * # 1 : Axial - * # -1 : Axial invert - * # 2 : Coronal - * # -2 : Coronal invert - * # 3 : Sagital - * # -3 : Sagital invert - * # 4 : Heart Axial - * # -4 : Heart Axial invert - * # 5 : Heart Coronal - * # -5 : Heart Coronal invert - * # 6 : Heart Sagital - * # -6 : Heart Sagital invert - */ -double File::TypeOrientation( ) -{ - float iop[6]; - bool succ = GetImageOrientationPatient( iop ); - if ( !succ ) - { - return 0.; - } - - vector3D ori1; - vector3D ori2; - - ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2]; - ori1.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5]; - - // two perpendicular vectors describe one plane - double dicPlane[6][2][3] = - { { {1, 0, 0 },{0, 1, 0 } }, // Axial - { {1, 0, 0 },{0, 0, -1 } }, // Coronal - { {0, 1, 0 },{0, 0, -1 } }, // Sagittal - { { 0.8, 0.5, 0.0 },{-0.1, 0.1 , -0.95 } }, // Axial - HEART - { { 0.8, 0.5, 0.0 },{-0.6674, 0.687, 0.1794} }, // Coronal - HEART - { {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794} } // Sagittal - HEART - }; - - vector3D refA; - vector3D refB; - int i = 0; - Res res; // [ , ] - res.first = 0; - res.second = 99999; - for (int numDicPlane=0; numDicPlane<6; numDicPlane++) - { - ++i; - // refA=plane[0] - refA.x = dicPlane[numDicPlane][0][0]; - refA.y = dicPlane[numDicPlane][0][1]; - refA.z = dicPlane[numDicPlane][0][2]; - // refB=plane[1] - refB.x = dicPlane[numDicPlane][1][0]; - refB.y = dicPlane[numDicPlane][1][1]; - refB.z = dicPlane[numDicPlane][1][2]; - res=VerfCriterion( i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res ); - res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res ); - } - return res.first; -/* -// i=0 -// res=[0,99999] ## [ , ] -// for plane in dicPlane: -// i=i+1 -// refA=plane[0] -// refB=plane[1] -// res=self.VerfCriterion( i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res ) -// res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res ) -// return res[0] -*/ - -} - -Res -File::VerfCriterion(int typeCriterion, double criterionNew, Res const & in) -{ - Res res; - double criterion = in.second; - if (criterionNew < criterion) - { - res.first = criterionNew; - res.second = typeCriterion; - } -/* -// type = res[0] -// criterion = res[1] -// # if criterionNew<0.1 and criterionNew : - type : vector 3D (double) -//- : - type : vector 3D (double) -// - Description of the first plane -//- : - type : vector 3D (double) -//- : - type : vector 3D (double) -// - Description of the second plane -//------------------------- Return : ------------------------------------ -// double : 0 if the planes are perpendicular. While the difference of -// the orientation between the planes are big more enlarge is -// the criterion. -//------------------------- Other : ------------------------------------- -// The calculus is based with vectors normalice -double -File::CalculLikelyhood2Vec(vector3D const & refA, vector3D const & refB, - vector3D const & ori1, vector3D const & ori2 ) -{ - - vector3D ori3 = ProductVectorial(ori1,ori2); - vector3D refC = ProductVectorial(refA,refB); - double res = square_dist(refC, ori3); - - return sqrt(res); -} - -//------------------------- Purpose : ----------------------------------- -//- Calculus of the poduct vectorial between two vectors 3D -//------------------------- Parameters : -------------------------------- -//- : - type : vector 3D (double) -//- : - type : vector 3D (double) -//------------------------- Return : ------------------------------------ -// (vec) : - Vector 3D -//------------------------- Other : ------------------------------------- -vector3D -File::ProductVectorial(vector3D const & vec1, vector3D const & vec2) -{ - vector3D vec3; - vec3.x = vec1.y*vec2.z - vec1.z*vec2.y; - vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x); - vec3.z = vec1.x*vec2.y - vec1.y*vec2.x; - - return vec3; -} - //----------------------------------------------------------------------------- // Print diff --git a/src/gdcmFile.h b/src/gdcmFile.h index 2675a0d7..ef3fac78 100644 --- a/src/gdcmFile.h +++ b/src/gdcmFile.h @@ -3,8 +3,8 @@ Program: gdcm Module: $RCSfile: gdcmFile.h,v $ Language: C++ - Date: $Date: 2005/07/24 02:14:43 $ - Version: $Revision: 1.114 $ + Date: $Date: 2005/07/24 02:34:42 $ + Version: $Revision: 1.115 $ Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). All rights reserved. See Doc/License.txt or @@ -27,15 +27,6 @@ namespace gdcm { -typedef struct -{ - double x; - double y; - double z; -} vector3D; - -typedef std::pair Res; - class RLEFramesInfo; class JPEGFragmentsInfo; @@ -183,8 +174,6 @@ public: bool Write(std::string fileName, FileType filetype); - double TypeOrientation( ); - protected: /// Store the RLE frames info obtained during parsing of pixels. @@ -211,10 +200,6 @@ private: uint32_t ReadTagLength(uint16_t, uint16_t); void ReadAndSkipEncapsulatedBasicOffsetTable(); - Res VerfCriterion(int typeCriterion, double criterionNew, Res const & res); - double CalculLikelyhood2Vec(vector3D const & refA, vector3D const & refB, - vector3D const & ori1, vector3D const & ori2); - vector3D ProductVectorial(vector3D const & vec1, vector3D const & vec2); }; } // end namespace gdcm diff --git a/src/gdcmOrientation.cxx b/src/gdcmOrientation.cxx new file mode 100644 index 00000000..6632d4c8 --- /dev/null +++ b/src/gdcmOrientation.cxx @@ -0,0 +1,192 @@ +/*========================================================================= + + Program: gdcm + Module: $RCSfile: gdcmOrientation.cxx,v $ + Language: C++ + Date: $Date: 2005/07/24 02:34:42 $ + Version: $Revision: 1.1 $ + + Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de + l'Image). All rights reserved. See Doc/License.txt or + http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include "gdcmOrientation.h" +#include "gdcmFile.h" + +namespace gdcm +{ +//-------------------------------------------------------------------- +// THERALYS Algorithm to determine the most similar basic orientation +// +// Transliterated from original Python code. +// Kept as close as possible to the original code +// in order to speed up any further modif of Python code :-( +//----------------------------------------------------------------------- + +/** + * \brief THERALYS' Algorithm to determine the most similar basic orientation + * (Axial, Coronal, Sagital) of the image + * \note Should be run on the first gdcm::File of a 'coherent' Serie + * @return orientation code + * @return orientation code + * # 0 : Not Applicable (neither 0020,0037 Image Orientation Patient + * # nor 0020,0032Image Position found ) + * # 1 : Axial + * # -1 : Axial invert + * # 2 : Coronal + * # -2 : Coronal invert + * # 3 : Sagital + * # -3 : Sagital invert + * # 4 : Heart Axial + * # -4 : Heart Axial invert + * # 5 : Heart Coronal + * # -5 : Heart Coronal invert + * # 6 : Heart Sagital + * # -6 : Heart Sagital invert + */ +double Orientation::TypeOrientation( File *f ) +{ + float iop[6]; + bool succ = f->GetImageOrientationPatient( iop ); + if ( !succ ) + { + return 0.; + } + + vector3D ori1; + vector3D ori2; + + ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2]; + ori1.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5]; + + // two perpendicular vectors describe one plane + double dicPlane[6][2][3] = + { { {1, 0, 0 },{0, 1, 0 } }, // Axial + { {1, 0, 0 },{0, 0, -1 } }, // Coronal + { {0, 1, 0 },{0, 0, -1 } }, // Sagittal + { { 0.8, 0.5, 0.0 },{-0.1, 0.1 , -0.95 } }, // Axial - HEART + { { 0.8, 0.5, 0.0 },{-0.6674, 0.687, 0.1794} }, // Coronal - HEART + { {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794} } // Sagittal - HEART + }; + + vector3D refA; + vector3D refB; + int i = 0; + Res res; // [ , ] + res.first = 0; + res.second = 99999; + for (int numDicPlane=0; numDicPlane<6; numDicPlane++) + { + ++i; + // refA=plane[0] + refA.x = dicPlane[numDicPlane][0][0]; + refA.y = dicPlane[numDicPlane][0][1]; + refA.z = dicPlane[numDicPlane][0][2]; + // refB=plane[1] + refB.x = dicPlane[numDicPlane][1][0]; + refB.y = dicPlane[numDicPlane][1][1]; + refB.z = dicPlane[numDicPlane][1][2]; + res=VerfCriterion( i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res ); + res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res ); + } + return res.first; +/* +// i=0 +// res=[0,99999] ## [ , ] +// for plane in dicPlane: +// i=i+1 +// refA=plane[0] +// refB=plane[1] +// res=self.VerfCriterion( i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res ) +// res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res ) +// return res[0] +*/ + +} + +Res +Orientation::VerfCriterion(int typeCriterion, double criterionNew, Res const & in) +{ + Res res; + double criterion = in.second; + if (criterionNew < criterion) + { + res.first = criterionNew; + res.second = typeCriterion; + } +/* +// type = res[0] +// criterion = res[1] +// # if criterionNew<0.1 and criterionNew : - type : vector 3D (double) +//- : - type : vector 3D (double) +// - Description of the first plane +//- : - type : vector 3D (double) +//- : - type : vector 3D (double) +// - Description of the second plane +//------------------------- Return : ------------------------------------ +// double : 0 if the planes are perpendicular. While the difference of +// the orientation between the planes are big more enlarge is +// the criterion. +//------------------------- Other : ------------------------------------- +// The calculus is based with vectors normalice +double +Orientation::CalculLikelyhood2Vec(vector3D const & refA, vector3D const & refB, + vector3D const & ori1, vector3D const & ori2 ) +{ + + vector3D ori3 = ProductVectorial(ori1,ori2); + vector3D refC = ProductVectorial(refA,refB); + double res = square_dist(refC, ori3); + + return sqrt(res); +} + +//------------------------- Purpose : ----------------------------------- +//- Calculus of the poduct vectorial between two vectors 3D +//------------------------- Parameters : -------------------------------- +//- : - type : vector 3D (double) +//- : - type : vector 3D (double) +//------------------------- Return : ------------------------------------ +// (vec) : - Vector 3D +//------------------------- Other : ------------------------------------- +vector3D +Orientation::ProductVectorial(vector3D const & vec1, vector3D const & vec2) +{ + vector3D vec3; + vec3.x = vec1.y*vec2.z - vec1.z*vec2.y; + vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x); + vec3.z = vec1.x*vec2.y - vec1.y*vec2.x; + + return vec3; +} + +} // end namespace gdcm + diff --git a/src/gdcmOrientation.h b/src/gdcmOrientation.h new file mode 100644 index 00000000..a43c6e17 --- /dev/null +++ b/src/gdcmOrientation.h @@ -0,0 +1,56 @@ +/*========================================================================= + + Program: gdcm + Module: $RCSfile: gdcmOrientation.h,v $ + Language: C++ + Date: $Date: 2005/07/24 02:34:42 $ + Version: $Revision: 1.1 $ + + Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de + l'Image). All rights reserved. See Doc/License.txt or + http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#ifndef GDCMORIENTATION_H +#define GDCMORIENTATION_H + +#include "gdcmBase.h" + +namespace gdcm +{ +typedef struct +{ + double x; + double y; + double z; +} vector3D; + +typedef std::pair Res; +class File; + +//----------------------------------------------------------------------------- +/** + * \brief Orientation class for dealing with DICOM image orienation + */ +class GDCM_EXPORT Orientation : public Base +{ +public: + Orientation(); + ~Orientation(); + + double TypeOrientation( File* file ); + +private: + Res VerfCriterion(int typeCriterion, double criterionNew, Res const & res); + double CalculLikelyhood2Vec(vector3D const & refA, vector3D const & refB, + vector3D const & ori1, vector3D const & ori2); + vector3D ProductVectorial(vector3D const & vec1, vector3D const & vec2); +}; +} // end namespace gdcm +//----------------------------------------------------------------------------- +#endif -- 2.48.1