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

Modified PlaneFittingPoints to return eighevalues and eigenvectors and added backward compatibility wrapper.

This commit is contained in:
Marco Di Benedetto 2009-10-05 22:43:14 +00:00
parent ca16dcdf52
commit e54e0a7124
1 changed files with 45 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
vcg/space/fitting3.h
View File

@ -46,52 +46,61 @@ created
namespace vcg { namespace vcg {
template <class S> template <class S>
Point3<S> PlaneFittingPoints( std::vector< Point3<S> > & samples,Plane3<S> &p){ Point3<S> PlaneFittingPoints(const std::vector< Point3<S> > & samples, Plane3<S> & p, Point4<S> & eval, Matrix44<S> & evec)
{
int j;
Matrix44<S> m;m.SetZero(); Matrix44<S> m;m.SetZero();
typename std::vector< Point3<S> > ::iterator i; typename std::vector< Point3<S> >::const_iterator it;
Point3<S> c; c.SetZero();
for(i = samples.begin(); i != samples.end(); ++i)
c+=*i;
c/=samples.size();
for(i = samples.begin(); i != samples.end(); ++i){ Point3<S> c; c.SetZero();
Point3<S> p = (*i)-c; for(it = samples.begin(); it != samples.end(); ++it)
for(j = 0 ; j < 3;++j) {
*(Point3<S>*)&m[j][0] += p * p[j]; c += *it;
}
c /= samples.size();
for(it = samples.begin(); it != samples.end(); ++it)
{
Point3<S> p = (*it) - c;
for(int j = 0 ; j < 3;++j)
*(Point3<S>*)&m[j][0] += p * p[j];
} }
m[0][3]= m[1][3]=m[2][3]=0.0; m[0][3] = m[1][3] = m[2][3] = S(0);
m[3][3]= 1.0; m[3][3] = S(1);
m[3][0]= m[3][1]=m[3][2]=0.0; m[3][0] = m[3][1] = m[3][2] = S(0);
int n; int n;
Matrix44<S> res;
Point4<S> e;
Point3<S> d; Point3<S> d;
Jacobi(m,e,res,n); Jacobi(m, eval, evec, n);
//Sort eigenvalues (tarinisort) //Sort eigenvalues (tarinisort)
e[0] = fabs(e[0]); Point4<S> e;
e[1] = fabs(e[1]); e[0] = S(fabs(eval[0]));
e[2] = fabs(e[2]); e[1] = S(fabs(eval[1]));
Point3<S> eval; e[2] = S(fabs(eval[2]));
int maxi,mini,medi;
if (e[1] > e[0]) { maxi=1; mini=0; } else { maxi=0; mini=1;}
if (e[maxi] < e[2]) maxi=2; else if(e[mini] > e[2]) mini=2;
medi = 3 - maxi -mini;
eval = Point3<S>(e[mini], e[medi], e[maxi]);
d[0]=res[0][mini];
d[1]=res[1][mini];
d[2]=res[2][mini];
p.SetOffset(c.dot(d)/d.Norm()); int maxi, mini, medi;
p.SetDirection(d/d.Norm()); if (e[1] > e[0]) { maxi = 1; mini = 0; } else { maxi = 0; mini = 1;}
if (e[maxi] < e[2]) { maxi = 2; } else if (e[mini] > e[2]) { mini = 2; };
medi = 3 - maxi -mini;
return eval; d[0] = evec[0][mini];
d[1] = evec[1][mini];
d[2] = evec[2][mini];
const S norm = d.Norm();
p.SetOffset(c.dot(d)/norm);
p.SetDirection(d/norm);
return Point3<S>(e[mini], e[medi], e[maxi]);
}
template <class S>
Point3<S> PlaneFittingPoints(const std::vector< Point3<S> > & samples, Plane3<S> & p)
{
Point4<S> eval;
Matrix44<S> evec;
return PlaneFittingPoints(samples, p, eval, evec);
} }
template<class S> template<class S>