Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

ShapeAnalysis.hpp

Go to the documentation of this file.
00001 /* ImLib3D
00002  * Copyright (c) 2001, ULP-IPB Strasbourg.
00003  *
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2 of the License, or (at
00007  * your option) any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017  */
00018 
00021 #ifndef _ShapeAnalysis_hpp
00022 #define _ShapeAnalysis_hpp
00023 
00024 #include<ImLib3D/ImageProcessor.hpp>
00025 #include<ImLib3D/Matrix.hpp>
00026 
00027 
00028 namespace IP3D
00029 {
00037     template <class ImageType>
00038     void CenterOfGravity(const ImageType& src,double &mass,Vect3Df &centerOfGravity)
00039     {
00040         mass = src.Sum();
00041         centerOfGravity.x=0;
00042         centerOfGravity.y=0;
00043         centerOfGravity.z=0;
00044         //SymmetricMatrix inertiaTensor(3); //! for object orientation
00045         typename ImageType::const_iteratorXYZ p;
00046         for (p=src.begin();p!=src.end();++p)
00047         {
00048             centerOfGravity.x += (*p)*p.x;
00049             centerOfGravity.y += (*p)*p.y;
00050             centerOfGravity.z += (*p)*p.z;
00051         }
00052         centerOfGravity.x /= mass;
00053         centerOfGravity.y /= mass;
00054         centerOfGravity.z /= mass;
00055     }
00056 };
00057 
00058 namespace IP3D
00059 {
00069     template <class ImageType>
00070     void PrincipalAxes(const ImageType& src,vector<Vect3Df>& pprincipalAxes,
00071                       double *_mass=NULL,
00072                       Vect3Df *_centerOfGravity=NULL)
00073     {
00074         double mass;
00075         Vect3Df centerOfGravity;
00076         if (_mass)
00077             mass=*_mass;
00078         else
00079             mass = src.Sum();
00080         if (!_centerOfGravity)
00081         {
00082             centerOfGravity.x=0;
00083             centerOfGravity.y=0;
00084             centerOfGravity.z=0;
00085         }
00086         else centerOfGravity=*_centerOfGravity;
00087         Matrix inertiaTensor; 
00088         typename ImageType::const_iteratorXYZ p;
00089         for (p=src.begin();p!=src.end();++p)
00090         {
00091             if (!_centerOfGravity)
00092             {
00093                 centerOfGravity.x += (*p)*p.x;
00094                 centerOfGravity.y += (*p)*p.y;
00095                 centerOfGravity.z += (*p)*p.z;
00096             }
00097             inertiaTensor(1,1)+= (*p)*p.x*p.x;
00098             inertiaTensor(2,2)+= (*p)*p.y*p.y;
00099             inertiaTensor(3,3)+= (*p)*p.z*p.z;
00100             inertiaTensor(2,1)+= (*p)*p.x*p.y;
00101             inertiaTensor(3,1)+= (*p)*p.x*p.z;
00102             inertiaTensor(3,2)+= (*p)*p.z*p.y;
00103         }
00104         if (!_centerOfGravity)
00105         {
00106             centerOfGravity.x /= mass;
00107             centerOfGravity.y /= mass;
00108             centerOfGravity.z /= mass;
00109         }
00110         inertiaTensor(1,1) -= centerOfGravity.x*centerOfGravity.x*mass;
00111         inertiaTensor(2,2) -= centerOfGravity.y*centerOfGravity.y*mass;
00112         inertiaTensor(3,3) -= centerOfGravity.z*centerOfGravity.z*mass;
00113             
00114         inertiaTensor(2,1) -= centerOfGravity.x*centerOfGravity.y*mass;
00115         inertiaTensor(3,1) -= centerOfGravity.x*centerOfGravity.z*mass;
00116         inertiaTensor(3,2) -= centerOfGravity.z*centerOfGravity.y*mass;
00117             
00118         inertiaTensor(2,1) *= (-1);
00119         inertiaTensor(3,1) *= (-1);
00120         inertiaTensor(3,2) *= (-1);
00121 
00122         vector<Vect3Df> &principalAxes=pprincipalAxes;
00123         principalAxes.clear();
00124         principalAxes.resize(3);
00125         Vect3Df &comp1=principalAxes[0];
00126         Vect3Df &comp2=principalAxes[1];
00127         Vect3Df &comp3=principalAxes[2];
00128         
00129         //      DiagonalMatrix eigenValues;
00130         Matrix eigenValues;
00131         Matrix eigenVectors;
00132         //more stable
00133         //Jacobi(inertiaTensor, eigenValues, eigenVectors);
00134         comp1.x = eigenVectors(1,1);
00135         comp1.y = eigenVectors(2,1);
00136         comp1.z = eigenVectors(3,1);
00137             
00138         comp2.x = eigenVectors(1,2);
00139         comp2.y = eigenVectors(2,2);
00140         comp2.z = eigenVectors(3,2);
00141             
00142         comp3.x = eigenVectors(1,3);
00143         comp3.y = eigenVectors(2,3);
00144         comp3.z = eigenVectors(3,3);
00145         
00146         // Sort in descending order
00147         if (eigenValues(1,1) < eigenValues(2,2)) // switch 1 and 2 ?
00148         {
00149             swap(comp1,comp2);
00150             
00151             float t=eigenValues(1,1);
00152             eigenValues(1,1)=eigenValues(2,2);
00153             eigenValues(2,2)=t;
00154         }
00155         // now there are only 3 possibilities left:
00156         if (eigenValues(1,1) < eigenValues(3,3)){swap(comp1,comp3);}
00157         else 
00158             if (eigenValues(2,2) < eigenValues(3,3)){swap(comp2,comp3);}
00159     }
00160 };
00161 
00162 namespace IP3D
00163 {
00170         void BoundingBox(const Mask3D& mask, RectZone3Di &bbox);
00171 };
00172 
00173 
00174 #endif //_ShapeAnalysis_hpp

Generated on Fri Jun 17 13:36:07 2005 for ImLib3D by  doxygen 1.4.2