From 1140ca5a32b3935f6f0b73baf1dd21c412884a42 Mon Sep 17 00:00:00 2001 From: cnr-isti-vclab Date: Fri, 23 Jul 2010 15:48:45 +0000 Subject: [PATCH] Added multiscale curvature computation (UpdateCurvatureLocal) --- vcg/complex/trimesh/nring.h | 4 +- .../trimesh/update/curvature_fitting.h | 390 ++++++++++++++++++ 2 files changed, 392 insertions(+), 2 deletions(-) diff --git a/vcg/complex/trimesh/nring.h b/vcg/complex/trimesh/nring.h index 4b5b857f..0397898a 100644 --- a/vcg/complex/trimesh/nring.h +++ b/vcg/complex/trimesh/nring.h @@ -61,14 +61,14 @@ public: MeshType* m; - RingWalker(VertexType* v, MeshType* m) : m(m) + Nring(VertexType* v, MeshType* m) : m(m) { assert((v - &*m->vert.begin()) < m->vert.size()); insertAndFlag(v); } - ~RingWalker() + ~Nring() { clear(); } diff --git a/vcg/complex/trimesh/update/curvature_fitting.h b/vcg/complex/trimesh/update/curvature_fitting.h index 1fdf08e1..f8cb9e9c 100644 --- a/vcg/complex/trimesh/update/curvature_fitting.h +++ b/vcg/complex/trimesh/update/curvature_fitting.h @@ -44,7 +44,11 @@ #include #include #include +#include +// GG include +#include +#include namespace vcg { @@ -310,6 +314,392 @@ class Quadric // ---- end Eigen stuff } } + + // GG LOCAL CURVATURE + + class QuadricLocal + { + public: + + QuadricLocal () + { + a() = b() = c() = d() = e() = 1.0; + } + + QuadricLocal (double av, double bv, double cv, double dv, double ev) + { + a() = av; + b() = bv; + c() = cv; + d() = dv; + e() = ev; + } + + double& a() { return data[0];} + double& b() { return data[1];} + double& c() { return data[2];} + double& d() { return data[3];} + double& e() { return data[4];} + + double data[5]; + + double evaluate(double u, double v) + { + return a()*u*u + b()*u*v + c()*v*v + d()*u + e()*v; + } + + double du(double u, double v) + { + return 2.0*a()*u + b()*v + d(); + } + + double dv(double u, double v) + { + return 2.0*c()*v + b()*u + e(); + } + + double duv(double u, double v) + { + return b(); + } + + double duu(double u, double v) + { + return 2.0*a(); + } + + double dvv(double u, double v) + { + return 2.0*c(); + } + + + static QuadricLocal fit(std::vector &VV, bool svdRes, bool detCheck) + { + assert(VV.size() >= 5); + Eigen::MatrixXd A(VV.size(),5); + Eigen::MatrixXd b(VV.size(),1); + Eigen::MatrixXd sol(5,1); + + for(unsigned int c=0; c < VV.size(); ++c) + { + double u = VV[c].X(); + double v = VV[c].Y(); + double n = VV[c].Z(); + + A(c,0) = u*u; + A(c,1) = u*v; + A(c,2) = v*v; + A(c,3) = u; + A(c,4) = v; + + b[c] = n; + } + + + static int count = 0, index = 0; + double min = 0.000000000001; //1.0e-12 + /* + if (!count) + printf("GNE %e\n", min); + */ + + if (detCheck && ((A.transpose()*A).determinant() < min && (A.transpose()*A).determinant() > -min)) + { + //A.svd().solve(b, &sol); A.svd().solve(b, &sol); + //cout << sol << endl; + printf("Quadric: unsolvable vertex %d %d\n", count, ++index); + //return Quadric (1, 1, 1, 1, 1); + A.svd().solve(b, &sol); + return QuadricLocal(sol[0],sol[1],sol[2],sol[3],sol[4]); + } + count++; + + //for (int i = 0; i < 100; i++) + { + if (svdRes) + A.svd().solve(b, &sol); + else + sol = ((A.transpose()*A).inverse()*A.transpose())*b; + + } + + return QuadricLocal(sol[0],sol[1],sol[2],sol[3],sol[4]); + } + }; + + static void expandMaxLocal (MeshType & mesh, VertexType *v, int max, std::vector *vv) + { + Nring rw = Nring (v, &mesh); + do rw.expand (); while (rw.allV.size() < max+1); + if (rw.allV[0] != v) + printf ("rw.allV[0] != *v\n"); + vv->reserve ((size_t)max); + for (int i = 1; i < max+1; i++) + vv->push_back(rw.allV[i]); + + rw.clear(); + } + + + static void expandSphereLocal (MeshType & mesh, VertexType *v, float r, int min, std::vector *vv) + { + Nring rw = Nring (v, &mesh); + + bool isInside = true; + while (isInside) + { + rw.expand(); + vv->reserve(rw.allV.size()); + + typename std::vector::iterator b = rw.lastV.begin(); + typename std::vector::iterator e = rw.lastV.end(); + isInside = false; + while(b != e) + { + if (((*b)->P() - v->P()).Norm() < r) + { + vv->push_back(*b);; + isInside = true; + } + ++b; + } + } + //printf ("%d\n", vv->size()); + rw.clear(); + + if (vv->size() < min) + { + vv->clear(); + expandMaxLocal (mesh, v, min, vv); + } + } + + + static void getAverageNormal (VertexType *vp, std::vector &vv, CoordType *ppn) + { + *ppn = CoordType (0,0,0); + for (typename std::vector::iterator vpi = vv.begin(); vpi != vv.end(); ++vpi) + *ppn += (*vpi)->N(); + *ppn += (*vp).N(); + *ppn /= vv.size() + 1; + ppn->Normalize(); + } + + + static void applyProjOnPlane (CoordType ppn, std::vector &vin, std::vector *vout) + { + for (typename std::vector::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi) + if ((*vpi)->N() * ppn > 0.0f) + vout->push_back (*vpi); + } + + static CoordType projectLocal(VertexType* v, VertexType* vp, CoordType ppn) + { + return vp->P() - (ppn * ((vp->P() - v->P()) * ppn)); + } + + + static void computeReferenceFramesLocal (VertexType *v, CoordType ppn, std::vector *ref) + { + vcg::face::VFIterator vfi (v); + + int i = (vfi.I() + 1) % 3; + VertexTypeP vp = vfi.F()->V(i); + + CoordType x = (projectLocal (v, vp, ppn) - v->P()).Normalize(); + + assert(fabs(x * ppn) < 0.1); + + *ref = std::vector(3); + + (*ref)[0] = x; + (*ref)[1] = (ppn ^ (*ref)[0]).Normalize(); + (*ref)[2] = ppn.Normalize(); //ppn / ppn.Norm(); + } + + + static void fitQuadricLocal (VertexType *v, std::vector ref, std::vector &vv, QuadricLocal *q) + { + bool svdResolution = false; + bool zeroDeterminantCheck = false; + + std::vector points; + points.reserve (vv.size()); + + typename std::vector::iterator b = vv.begin(); + typename std::vector::iterator e = vv.end(); + + while(b != e) + { + CoordType cp = (*b)->P(); + + // vtang non e` il v tangente!!! + CoordType vTang = cp - v->P(); + + double x = vTang * ref[0]; + double y = vTang * ref[1]; + double z = vTang * ref[2]; + points.push_back(CoordType(x,y,z)); + ++b; + } + + *q = QuadricLocal::fit (points, svdResolution, zeroDeterminantCheck); + } + + + static void finalEigenStuff (VertexType *v, std::vector ref, QuadricLocal q) + { + double a = q.a(); + double b = q.b(); + double c = q.c(); + double d = q.d(); + double e = q.e(); + + double E = 1.0 + d*d; + double F = d*e; + double G = 1.0 + e*e; + + CoordType n = CoordType(-d,-e,1.0).Normalize(); + + v->N() = ref[0] * n[0] + ref[1] * n[1] + ref[2] * n[2]; + + double L = 2.0 * a * n.Z(); + double M = b * n.Z(); + double N = 2 * c * n.Z(); + + // ----------------- Eigen stuff + Eigen::Matrix2d m; + m << L*G - M*F, M*E-L*F, M*E-L*F, N*E-M*F; + m = m / (E*G-F*F); + Eigen::SelfAdjointEigenSolver eig(m); + + Eigen::Vector2d c_val = eig.eigenvalues(); + Eigen::Matrix2d c_vec = eig.eigenvectors(); + + c_val = -c_val; + + CoordType v1, v2; + v1[0] = c_vec[0]; + v1[1] = c_vec[1]; + v1[2] = d * v1[0] + e * v1[1]; + + v2[0] = c_vec[2]; + v2[1] = c_vec[3]; + v2[2] = d * v2[0] + e * v2[1]; + + v1 = v1.Normalize(); + v2 = v2.Normalize(); + + CoordType v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2]; + CoordType v2global = ref[0] * v2[0] + ref[1] * v2[1] + ref[2] * v2[2]; + + v1global.Normalize(); + v2global.Normalize(); + + v1global *= c_val[0]; + v2global *= c_val[1]; + + if (c_val[0] > c_val[1]) + { + (*v).PD1() = v1global; + (*v).PD2() = v2global; + (*v).K1() = c_val[0]; + (*v).K2() = c_val[1]; + } + else + { + (*v).PD1() = v2global; + (*v).PD2() = v1global; + (*v).K1() = c_val[1]; + (*v).K2() = c_val[0]; + } + // ---- end Eigen stuff + } + + + + static void updateCurvatureLocal (MeshType & mesh, float radiusSphere) + { + bool verbose = false; + bool projectionPlaneCheck = true; + int vertexesPerFit = 0; + + int i = 0; + VertexIterator vi; + for(vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi, i++) + { + std::vector vv; + std::vector vvtmp; + + int count; + if (verbose && !((count = (vi - mesh.vert.begin())) % 1000)) + printf ("vertex %d of %d\n",count,mesh.vert.size()); + + // if (kRing != 0) + // expandRing (&*vi, kRing, 5, &vv); + // else + expandSphereLocal (mesh, &*vi, radiusSphere, 5, &vv); + + assert (vv.size() >= 5); + + CoordType ppn; + // if (averageNormalMode) + // //ppn = (*vi).N(); + getAverageNormal (&*vi, vv, &ppn); + // else + // getProjPlaneNormal (&*vi, vv, &ppn); + + if (projectionPlaneCheck) + { + vvtmp.reserve (vv.size ()); + applyProjOnPlane (ppn, vv, &vvtmp); + if (vvtmp.size() >= 5) + vv = vvtmp; + } + + vvtmp.clear(); + + // if (montecarloMaxVertexNum) + // { + // //printf ("P: %d\n", vv.size()); + // vvtmp.reserve (vv.size ()); + // //printf ("TP: %d\n", vvtmp.size()); + // applyMontecarlo (montecarloMaxVertexNum, vv, &vvtmp); + // //printf ("TD: %d\n", vvtmp.size()); + // vv = vvtmp; + // //printf ("D: %d\n", vv.size()); + // //printf ("\n"); + // } + + assert (vv.size() >= 5); + + std::vector ref; + computeReferenceFramesLocal (&*vi, ppn, &ref); + + /* + printf ("%lf %lf %lf - %lf %lf %lf - %lf %lf %lf\n", + ref[0][0], ref[0][1], ref[0][2], + ref[1][0], ref[1][1], ref[1][2], + ref[2][0], ref[2][1], ref[2][2]); + */ + + vertexesPerFit += vv.size(); + //printf ("size: %d\n", vv.size()); + + QuadricLocal q; + fitQuadricLocal (&*vi, ref, vv, &q); + + finalEigenStuff (&*vi, ref, q); + + } + + //if (verbose) + //printf ("average vertex num in each fit: %f, total %d, vn %d\n", ((float) vertexesPerFit) / mesh.vn, vertexesPerFit, mesh.vn); + if (verbose) + printf ("average vertex num in each fit: %f\n", ((float) vertexesPerFit) / mesh.vn); + } + }; }