manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

corrected wrong constant in Covariance (thanks F.Ponchio) 

and minor calculation simplifications
This commit is contained in:
ganovelli 2008-09-10 16:18:32 +00:00
parent b9ce07204e
commit 4971b69b13
1 changed files with 16 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
vcg/complex/trimesh/inertia.h
View File

@ -75,6 +75,7 @@ class Inertia
typedef typename MeshType::FaceType FaceType; typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::FacePointer FacePointer; typedef typename MeshType::FacePointer FacePointer;
typedef typename MeshType::FaceIterator FaceIterator; typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::ConstFaceIterator ConstFaceIterator;
typedef typename MeshType::FaceContainer FaceContainer; typedef typename MeshType::FaceContainer FaceContainer;
typedef typename MeshType::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
@ -307,10 +308,10 @@ void InertiaTensorEigen(Matrix44<ScalarType> &EV, Point4<ScalarType> &ev )
/** Compute covariance matrix of a mesh, i.e. the integral /** Compute covariance matrix of a mesh, i.e. the integral
int_{M} { (x-b)(x-b)^T }dx where b is the barycenter and x spans over the mesh M int_{M} { (x-b)(x-b)^T }dx where b is the barycenter and x spans over the mesh M
*/ */
static void Covariance(MeshType m, vcg::Point3<ScalarType> & bary, vcg::Matrix33<ScalarType> &C){ static void Covariance(const MeshType & m, vcg::Point3<ScalarType> & bary, vcg::Matrix33<ScalarType> &C){
// find the barycenter // find the barycenter
FaceIterator fi; ConstFaceIterator fi;
ScalarType area = 0.0; ScalarType area = 0.0;
bary.Zero(); bary.Zero();
for(fi = m.face.begin(); fi != m.face.end(); ++fi) for(fi = m.face.begin(); fi != m.face.end(); ++fi)
@ -320,7 +321,8 @@ static void Covariance(MeshType m, vcg::Point3<ScalarType> & bary, vcg::Matrix33
area+=vcg::DoubleArea(*fi); area+=vcg::DoubleArea(*fi);
} }
bary/=area; bary/=area;
C.SetZero(); C.SetZero();
// C as covariance of triangle (0,0,0)(1,0,0)(0,1,0) // C as covariance of triangle (0,0,0)(1,0,0)(0,1,0)
vcg::Matrix33<ScalarType> C0; vcg::Matrix33<ScalarType> C0;
@ -330,26 +332,26 @@ static void Covariance(MeshType m, vcg::Point3<ScalarType> & bary, vcg::Matrix33
C0*=1/24.0; C0*=1/24.0;
// integral of (x,y,0) in the same triangle // integral of (x,y,0) in the same triangle
CoordType X(2/3.0,2/3.0,0); CoordType X(1/6.0,1/6.0,0);
vcg::Matrix33<ScalarType> A,At,DC; // matrix that bring the vertices to (v1-v0,v2-v0,n) vcg::Matrix33<ScalarType> A, // matrix that bring the vertices to (v1-v0,v2-v0,n)
At,DC;
for(fi = m.face.begin(); fi != m.face.end(); ++fi) for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if(!(*fi).IsD()) if(!(*fi).IsD())
{ {
CoordType n = (*fi).N().Normalize();
const CoordType &P0 = (*fi).P(0); const CoordType &P0 = (*fi).P(0);
const CoordType &P1 = (*fi).P(1); const CoordType &P1 = (*fi).P(1);
const CoordType &P2 = (*fi).P(2); const CoordType &P2 = (*fi).P(2);
CoordType n = ((P1-P0)^(P2-P0));
const float da = n.Norm();
n/=da*da;
A.SetColumn(0,P1-P0); A.SetColumn(0,P1-P0);
A.SetColumn(1,P2-P0); A.SetColumn(1,P2-P0);
A.SetColumn(2,n); A.SetColumn(2,n);
CoordType delta = P0 - bary; CoordType delta = P0 - bary;
At= A; At= A;
At.Transpose(); At.Transpose();
/* DC is calculated as integral of (A*x+delta) * (A*x+delta)^T over the triangle, /* DC is calculated as integral of (A*x+delta) * (A*x+delta)^T over the triangle,
where delta = v0-bary where delta = v0-bary
*/ */
@ -357,13 +359,14 @@ static void Covariance(MeshType m, vcg::Point3<ScalarType> & bary, vcg::Matrix33
DC.SetZero(); DC.SetZero();
DC+= A*C0*At; DC+= A*C0*At;
vcg::Matrix33<ScalarType> tmp; vcg::Matrix33<ScalarType> tmp;
tmp.OuterProduct(X,delta); tmp.OuterProduct(A*X,delta);
DC+= A * tmp; DC+= tmp;
tmp.Transpose(); tmp.Transpose();
DC+= tmp * At; DC+= tmp;
tmp.OuterProduct(delta,delta); tmp.OuterProduct(delta,delta);
DC+=tmp*0.5; DC+=tmp*0.5;
DC*=fabs(A.Determinant()); // the determinant of A is the jacobian of the change of variables A*x+delta // DC*=fabs(A.Determinant()); // the determinant of A is the jacobian of the change of variables A*x+delta
DC*=da; // the determinant of A is also the double area of *fi
C+=DC; C+=DC;
} }