diff --git a/src/meshlabplugins/editsegment/curvaturetensor.h b/src/meshlabplugins/editsegment/curvaturetensor.h index 699a5cafa..c44b6195d 100644 --- a/src/meshlabplugins/editsegment/curvaturetensor.h +++ b/src/meshlabplugins/editsegment/curvaturetensor.h @@ -1,5 +1,7 @@ #include #include +#include + #include #include #include @@ -12,10 +14,10 @@ namespace vcg { class CurvData { public: - Point3 T1; - Point3 T2; - double k1; - double k2; + Point3 T1; + Point3 T2; + float k1; + float k2; }; template class CurvatureTensor { @@ -32,9 +34,9 @@ namespace vcg { MESH_TYPE * mesh; SimpleTempData *TDCurvPtr; - void givens(double a, double b, double *c, double *s) + void givens(float a, float b, float *c, float *s) { - double tao; + float tao; if (b == 0) { *c = 1; *s = 0; @@ -53,16 +55,16 @@ namespace vcg { } } - void ComputePerVertexMatrix(VertexType & i, VertexType & j, Matrix33& Mvi, Matrix33 & n_nMatrix, double wij) { - Point3 ViVj; - Point3 Tij; - double kij = 0.0f; - Matrix33 tempMatrix; + void ComputePerVertexMatrix(VertexType & i, VertexType & j, Matrix33& Mvi, Matrix33 & n_nMatrix, float wij) { + Point3 ViVj; + Point3 Tij; + float kij = 0.0f; + Matrix33 tempMatrix; - Point3 tempViVj = i.P() - j.P(); - Point3 n = Point3((double)i.N()[0], (double)i.N()[1], (double)i.N()[2]); + + Point3 n = i.N(); n = n.Normalize(); - ViVj = Point3((double)tempViVj[0], (double)tempViVj[1], (double)tempViVj[2]); + ViVj = i.P() - j.P(); Tij = (n_nMatrix * ViVj) / Norm(n_nMatrix * ViVj); kij = (2 * (n * ViVj)) / Norm(ViVj); @@ -88,25 +90,26 @@ namespace vcg { for (vi = mesh->vert.begin(); vi != mesh->vert.end(); ++vi) { if (!vi->IsD() ) { - Matrix33 Mvi; + Matrix33 Mvi; Mvi.SetZero(); - Matrix33 I; + Matrix33 I; I.SetIdentity(); - Matrix33 n_nM; - Point3 normal_d = Point3((double)vi->N()[0], (double)vi->N()[1],(double)vi->N()[2]); - normal_d = normal_d.Normalize(); - n_nM.ExternalProduct(normal_d, normal_d); + Matrix33 n_nM; + //Point3 normal_d = Point3((double)vi->N()[0], (double)vi->N()[1],(double)vi->N()[2]); + Point3 normal = vi->N(); + normal = normal.Normalize(); + n_nM.ExternalProduct(normal, normal); n_nM = I - n_nM; //Inizio Calcolo dei wij per tutti i vertici adiacenti FaceType* first_face = vi->VFp(); vcg::face::Pos pos2(first_face, vi->VFi(), &(*vi)); - vector wij_container; - double totalDoubleAreaSize = 0; + vector wij_container; + float totalDoubleAreaSize = 0; do { - double doubleArea = vcg::DoubleArea(*pos2.F()); + float doubleArea = vcg::DoubleArea(*pos2.F()); totalDoubleAreaSize += doubleArea; wij_container.push_back(doubleArea); @@ -152,17 +155,17 @@ namespace vcg { //Mvi matrix ready for the vertex vi //calculate principal directions and curvature - Point3d Wvi; - Matrix33d Qvi; - Matrix33d QviT; - Matrix33d tempMatrix; + Point3f Wvi; + Matrix33f Qvi; + Matrix33f QviT; + Matrix33f tempMatrix; - Matrix33d A; - Point3d E1(1.0,0.0,0.0); - Point3d Nvi = Point3d::Construct((*vi).N()); + Matrix33f A; + Point3f E1(1.0,0.0,0.0); + Point3f Nvi = (*vi).N(); Nvi.Normalize(); - Point3d E1nNvi = E1 - Nvi; - Point3d E1pNvi = E1 + Nvi; + Point3f E1nNvi = E1 - Nvi; + Point3f E1pNvi = E1 + Nvi; if (E1nNvi.Norm() > E1pNvi.Norm() ) Wvi = E1nNvi / E1nNvi.Norm(); else Wvi = E1pNvi / E1pNvi.Norm(); @@ -176,25 +179,16 @@ namespace vcg { A = QviT * Mvi * Qvi; - //loop to set at 0 extremely small values - /*for (int i = 0; i<3; ++i) { - for (int j = 0; j<3; ++j) { - if (abs(A[i][j]) < 0.000001) { - A[i][j] = 0.0; - } - } - } */ - - double c; - double s; + float c; + float s; givens(A[1][1], A[2][1], &c, &s); - Point3 T1 = (A.GetColumn(1) * c) - (A.GetColumn(2) * s); - Point3 T2 = (A.GetColumn(1) * s) + (A.GetColumn(2) * c); + Point3 T1 = (A.GetColumn(1) * c) - (A.GetColumn(2) * s); + Point3 T2 = (A.GetColumn(1) * s) + (A.GetColumn(2) * c); - double k1 = 3*T1[1] - T2[2]; - double k2 = 3*T2[2] - T1[1]; + float k1 = 3*T1[1] - T2[2]; + float k2 = 3*T2[2] - T1[1]; (*TDCurvPtr)[*vi].T1 = T1; (*TDCurvPtr)[*vi].T2 = T2; @@ -203,7 +197,56 @@ namespace vcg { } } } +/* + //re-write, to be continued + void newComputeCurvatureTensor() { + vcg::tri::UpdateNormals::PerVertex(*mesh); + VertexIterator vi; + for (vi = mesh->vert.begin(); vi != mesh->vert.end(); ++vi) { + + if ( ! (*vi).IsD()) { + std::vector faces; + std::vector vertices; + std::vector weights; + std::vector vertices; + + vcg::face::JumpingPos pos((*vi).VFp(), (*vi)); + + VertexType* firstV = pos.VFlip(); + VertexType* tempV; + float totalDoubleAreaSize = 0.0f; + do + { + pos.NextE(); + tempV = pos.VFlip(); + + AdjVertex v; + + v.isBorder = pos.IsBorder() + v.vertex = tempV; + v.doubleArea = pos.f->DoubleArea(); + totalDoubleAreaSize += v.doubleArea; + + vertices.push_back(v); + } + while(tempV != firstV); + + std::vector Weights; + + for (int i = 0; i