ccolin/mod_geo-tp
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

switch back to an eigen-based MyMesh

Browse Source
This commit is contained in:
ccolin 2021-11-12 20:59:37 +01:00
parent 07093b6197
commit edcffa34c3
2 changed files with 54 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

218
src/curvature.cpp
View File

@ -21,196 +21,80 @@ void Courbures::normales_locales() {
} }
std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) {
{ OpenMesh::VPropHandleT<bool> vprop_flag;
OpenMesh::VPropHandleT<bool> vprop_flag; _mesh.add_property(vprop_flag, "vprop_flag");
_mesh.add_property(vprop_flag, "vprop_flag");
// Initialisation // Initialisation
for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end(); ++v_it) for (VertexHandle vh : _mesh.vertices()) {
{ _mesh.property(vprop_flag, vh) = false;
_mesh.property(vprop_flag, *v_it) = false ;
} }
// Circulateur sur le premier cercle // Circulateur sur le premier cercle
std::vector<MyMesh::VertexHandle> neigh, neigh2 ; std::vector<MyMesh::VertexHandle> neigh, neigh2;
_mesh.property(vprop_flag, vh) = true ; _mesh.property(vprop_flag, vh) = true;
for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it) for(VertexHandle vv : _mesh.vv_range(vh)) {
{ neigh.push_back(vv); // ajout du point à la liste
neigh.push_back(*vv_it) ; // ajout du point à la liste _mesh.property(vprop_flag, vv) = true;
_mesh.property(vprop_flag, *vv_it) = true ;
} }
// Parcours du premier cercle et ajout du second cercle par circulateurs // Parcours du premier cercle et ajout du second cercle par circulateurs
for (int i=0; i<neigh.size(); i++) for (size_t i = 0; i < neigh.size(); i++) {
{
MyMesh::VertexHandle vh = neigh.at(i) ; MyMesh::VertexHandle vh = neigh.at(i) ;
for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it) for (VertexHandle vv : _mesh.vv_range(vh)) {
{ if (!_mesh.property(vprop_flag, vv)) {
if (!_mesh.property(vprop_flag, *vv_it)) // sommet non encore rencontré neigh2.push_back(vv);
neigh2.push_back(*vv_it) ; // ajout du point à la liste }
_mesh.property(vprop_flag, *vv_it) = true ; _mesh.property(vprop_flag, vv) = true;
} }
} }
// Concaténation des deux cercles // Concaténation des deux cercles
neigh.insert(neigh.end(), neigh2.begin(), neigh2.end()) ; neigh.insert(neigh.end(), neigh2.begin(), neigh2.end());
_mesh.remove_property(vprop_flag); _mesh.remove_property(vprop_flag);
return neigh;
return neigh ;
} }
// std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) {
// OpenMesh::VPropHandleT<bool> 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<MyMesh::VertexHandle> 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) QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) {
{ std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh);
std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh) ; if (neigh.size() < 5) throw "Quad fitting: not enough neighbors";
// 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 ;
}
// Calcul de la matrice de changement de base // Calcul de la matrice de changement de base
MyMesh::Normal n = _mesh.normal(vh) ; Eigen::Vector3d Oz(0,0,1);
MyMesh::Point p = _mesh.point(vh); Eigen::Vector3d axis = _mesh.normal(vh).cross(Oz);
// std::cout << "-- normale" << std::endl ; double sina = axis.norm();
// std::cout << n[0] << ", " << n[1] << ", " << n[2] << std::endl ; double cosa = _mesh.normal(vh).dot(Oz);
// std::cout << "-- point" << std::endl ; double angle;
// std::cout << p[0] << ", " << p[1] << ", " << p[2] << std::endl ; if (sina >= 0) {
angle = acos(cosa);
Eigen::Vector3d ne(n[0], n[1], n[2]), Oz(0,0,1), axis; } else {
Eigen::Vector3d p_e(p[0], p[1], p[2]), pi_e ; angle = -acos(cosa);
}
axis = ne.cross(Oz) ; axis = axis.normalized();
double sina = axis.norm(), cosa = ne.dot(Oz), angle ; Eigen::AngleAxisd rot(angle, axis);
if (sina >= 0) Eigen::Translation3d trans(-_mesh.point(vh));
angle = acos(cosa) ; Eigen::Transform<double, 3, Eigen::Affine> ch_base = rot * trans;
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<double, 3, Eigen::Affine> ch_base = r * t ;
// Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen) // 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) // Résolution aux moindres carrés par SVD
{ Eigen::VectorXd coef(5);
int n(neigh.size()) ; coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B);
Eigen::MatrixXd A(n,5); return QuadPatch(coef, ch_base);
Eigen::VectorXd B(n);
// std::cout << "-- Après changement de base" << std::endl ;
for(int i=0; i<neigh.size(); i++)
{
MyMesh::Point p = _mesh.point(neigh.at(i)) ;
pi_e << p[0], p[1], p[2] ;
pi_e = ch_base * pi_e ; // Application du changement de base
// std::cout << pi_e[0] << ", " << pi_e[1] << ", " << pi_e[2] << std::endl ;
B[i] = -pi_e[2]; // -zi
A(i, 0) = pi_e[0] * pi_e[0]; // xi^2
A(i, 1) = pi_e[0] * pi_e[1]; // xi*yi
A(i, 2) = pi_e[1] * pi_e[1]; // yi^2
A(i, 3) = pi_e[0]; // xi
A(i, 4) = pi_e[1]; // yi
}
// Résolution aux moindres carrés par SVD
Eigen::VectorXd coef(5) ;
coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B) ;
QuadPatch q(coef, ch_base) ;
// std::cout << "-- quad : " << std::endl ;
// std::cout << q[0] << ", " << q[1] << ", " << q[2] << ", " << q[3] << ", " << q[4] << ", " << q[5] << ", " << std::endl ;
return q ;
}
else
{
std::cout << "Quad fitting : not enough neighbors" ;
throw "Quad fitting : not enough neighbors" ;
}
} }
// MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh) {
// std::vector<MyMesh::VertexHandle> 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<double, 3, Eigen::Affine> 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() { void Courbures::compute_KH() {

13
src/my_mesh.h
View File

@ -10,17 +10,10 @@
#include <OpenMesh/Core/Geometry/EigenVectorT.hh> #include <OpenMesh/Core/Geometry/EigenVectorT.hh>
// struct MyTraits : public OpenMesh::DefaultTraits {
// using Point = Eigen::Vector3<qreal>;
// using Normal = Eigen::Vector3<qreal>;
// using Color = Eigen::Vector3<qreal>;
// VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
// HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
// FaceAttributes(OpenMesh::Attributes::Normal);
// EdgeAttributes(OpenMesh::Attributes::Color);
// };
struct MyTraits : public OpenMesh::DefaultTraits { struct MyTraits : public OpenMesh::DefaultTraits {
typedef OpenMesh::Vec3f Color; using Point = Eigen::Vector3<qreal>;
using Normal = Eigen::Vector3<qreal>;
using Color = Eigen::Vector3<qreal>;
VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color); VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge); HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
FaceAttributes(OpenMesh::Attributes::Normal); FaceAttributes(OpenMesh::Attributes::Normal);