diff --git a/src/curvature.cpp b/src/curvature.cpp index 92c487d..936ca81 100644 --- a/src/curvature.cpp +++ b/src/curvature.cpp @@ -21,196 +21,80 @@ void Courbures::normales_locales() { } -std::vector Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) -{ - OpenMesh::VPropHandleT vprop_flag; - _mesh.add_property(vprop_flag, "vprop_flag"); +std::vector Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) { + OpenMesh::VPropHandleT vprop_flag; + _mesh.add_property(vprop_flag, "vprop_flag"); // Initialisation - for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end(); ++v_it) - { - _mesh.property(vprop_flag, *v_it) = false ; + for (VertexHandle vh : _mesh.vertices()) { + _mesh.property(vprop_flag, vh) = false; } // Circulateur sur le premier cercle - std::vector neigh, neigh2 ; + std::vector neigh, neigh2; - _mesh.property(vprop_flag, vh) = true ; - for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it) - { - neigh.push_back(*vv_it) ; // ajout du point à la liste - _mesh.property(vprop_flag, *vv_it) = true ; + _mesh.property(vprop_flag, vh) = true; + for(VertexHandle vv : _mesh.vv_range(vh)) { + neigh.push_back(vv); // ajout du point à la liste + _mesh.property(vprop_flag, vv) = true; } // Parcours du premier cercle et ajout du second cercle par circulateurs - for (int i=0; i Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) { -// OpenMesh::VPropHandleT vprop_flag; -// _mesh.add_property(vprop_flag, "vprop_flag"); - -// // Initialisation -// for (VertexHandle vh : _mesh.vertices()) { -// _mesh.property(vprop_flag, vh) = false; -// } -// // Circulateur sur le premier cercle -// std::vector neigh, neigh2; - -// _mesh.property(vprop_flag, vh) = true; -// for(VertexHandle vv : _mesh.vv_range(vh)) { -// neigh.push_back(vv); // ajout du point à la liste -// _mesh.property(vprop_flag, vv) = true; -// } - -// // Parcours du premier cercle et ajout du second cercle par circulateurs -// for (size_t i = 0; i < neigh.size(); i++) { -// MyMesh::VertexHandle vh = neigh.at(i) ; -// for (VertexHandle vv : _mesh.vv_range(vh)) { -// if (!_mesh.property(vprop_flag, vv)) { -// neigh2.push_back(vv); -// } -// _mesh.property(vprop_flag, vv) = true; -// } -// } - -// // Concaténation des deux cercles -// neigh.insert(neigh.end(), neigh2.begin(), neigh2.end()); -// _mesh.remove_property(vprop_flag); -// return neigh; -// } -QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) -{ - std::vector neigh = get_two_neighborhood(vh) ; - -// std::cout << "-- Neighborhood" << std::endl ; - for (int i = 0; i< neigh.size(); i++) - { - MyMesh::Point p = _mesh.point(neigh.at(i)) ; -// std::cout << p[0] << ", " << p[1] << ", " << p[2] << " ; " << std::endl ; - } +QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) { + std::vector neigh = get_two_neighborhood(vh); + if (neigh.size() < 5) throw "Quad fitting: not enough neighbors"; // Calcul de la matrice de changement de base - MyMesh::Normal n = _mesh.normal(vh) ; - MyMesh::Point p = _mesh.point(vh); -// std::cout << "-- normale" << std::endl ; -// std::cout << n[0] << ", " << n[1] << ", " << n[2] << std::endl ; -// std::cout << "-- point" << std::endl ; -// std::cout << p[0] << ", " << p[1] << ", " << p[2] << std::endl ; - - Eigen::Vector3d ne(n[0], n[1], n[2]), Oz(0,0,1), axis; - Eigen::Vector3d p_e(p[0], p[1], p[2]), pi_e ; - - axis = ne.cross(Oz) ; - double sina = axis.norm(), cosa = ne.dot(Oz), angle ; - if (sina >= 0) - angle = acos(cosa) ; - else - angle = -acos(cosa) ; - axis = axis.normalized() ; - - Eigen::AngleAxisd r(angle, axis) ; -// std::cout << "-- rotation" << std:: endl ; -// std::cout << r.matrix()(0,0) << ", " << r.matrix()(0, 1) << ", " << r.matrix()(0,2) << "; " ; -// std::cout << r.matrix()(1,0) << ", " << r.matrix()(1, 1) << ", " << r.matrix()(1,2) << "; " ; -// std::cout << r.matrix()(2,0) << ", " << r.matrix()(2, 1) << ", " << r.matrix()(2,2) << std::endl ; - Eigen::Translation3d t(-p_e[0], -p_e[1], -p_e[2]) ; - Eigen::Transform ch_base = r * t ; - + Eigen::Vector3d Oz(0,0,1); + Eigen::Vector3d axis = _mesh.normal(vh).cross(Oz); + double sina = axis.norm(); + double cosa = _mesh.normal(vh).dot(Oz); + double angle; + if (sina >= 0) { + angle = acos(cosa); + } else { + angle = -acos(cosa); + } + axis = axis.normalized(); + Eigen::AngleAxisd rot(angle, axis); + Eigen::Translation3d trans(-_mesh.point(vh)); + Eigen::Transform ch_base = rot * trans; // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen) + Eigen::MatrixXd A(neigh.size(), 5); + Eigen::VectorXd B(neigh.size()); + for(size_t i = 0; i < neigh.size(); i++) { + MyMesh::Point point = _mesh.point(neigh[i]); + point = ch_base * point; // Application du changement de base + B[i] = -point[2]; // -zi + A(i, 0) = point[0] * point[0]; // xi^2 + A(i, 1) = point[0] * point[1]; // xi*yi + A(i, 2) = point[1] * point[1]; // yi^2 + A(i, 3) = point[0]; // xi + A(i, 4) = point[1]; // yi + } - if (neigh.size() >= 5) - { - int n(neigh.size()) ; - Eigen::MatrixXd A(n,5); - Eigen::VectorXd B(n); - -// std::cout << "-- Après changement de base" << std::endl ; - for(int i=0; i neigh = get_two_neighborhood(vh); -// if (neigh.size() < 5) throw "Quad fitting: not enough neighbors"; - -// // Calcul de la matrice de changement de base -// Eigen::Vector3d Oz(0,0,1); -// Eigen::Vector3d axis = _mesh.normal(vh).cross(Oz); -// double sina = axis.norm(); -// double cosa = _mesh.normal(vh).dot(Oz); -// double angle; -// if (sina >= 0) { -// angle = acos(cosa); -// } else { -// angle = -acos(cosa); -// } -// axis = axis.normalized(); -// Eigen::AngleAxisd rot(angle, axis); -// Eigen::Translation3d trans(-_mesh.point(vh)); -// Eigen::Transform ch_base = rot * trans; - -// // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen) -// Eigen::MatrixXd A(neigh.size(), 5); -// Eigen::VectorXd B(neigh.size()); -// for(size_t i = 0; i < neigh.size(); i++) { -// MyMesh::Point point = _mesh.point(neigh[i]); -// point = ch_base * point; // Application du changement de base -// B[i] = -point[2]; // -zi -// A(i, 0) = point[0] * point[0]; // xi^2 -// A(i, 1) = point[0] * point[1]; // xi*yi -// A(i, 2) = point[1] * point[1]; // yi^2 -// A(i, 3) = point[0]; // xi -// A(i, 4) = point[1]; // yi -// } - -// // Résolution aux moindres carrés par SVD -// Eigen::VectorXd coef(5); -// coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B); -// return MyQuad(coef[0], coef[1], coef[2], coef[3], coef[4], rot, trans); -// } void Courbures::compute_KH() { diff --git a/src/my_mesh.h b/src/my_mesh.h index 48b35ba..79db8d0 100644 --- a/src/my_mesh.h +++ b/src/my_mesh.h @@ -10,17 +10,10 @@ #include -// struct MyTraits : public OpenMesh::DefaultTraits { -// using Point = Eigen::Vector3; -// using Normal = Eigen::Vector3; -// using Color = Eigen::Vector3; -// VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color); -// HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge); -// FaceAttributes(OpenMesh::Attributes::Normal); -// EdgeAttributes(OpenMesh::Attributes::Color); -// }; struct MyTraits : public OpenMesh::DefaultTraits { - typedef OpenMesh::Vec3f Color; + using Point = Eigen::Vector3; + using Normal = Eigen::Vector3; + using Color = Eigen::Vector3; VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color); HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge); FaceAttributes(OpenMesh::Attributes::Normal);