diff --git a/k.png b/k.png new file mode 100644 index 0000000..aeb95ff Binary files /dev/null and b/k.png differ diff --git a/patch.png b/patch.png new file mode 100644 index 0000000..364f06e Binary files /dev/null and b/patch.png differ diff --git a/tp.md b/tp.md new file mode 100644 index 0000000..832ff40 --- /dev/null +++ b/tp.md @@ -0,0 +1,116 @@ +# TP - Courbures discrètes +**Cyril Colin et Jérémy André** + +## Exercice 1 +Récupération des voisins d'un sommet : +Si le nombre de sommets dans le 1-voisinage est < 5, alors on ajoute les points du 2-voisinage pour calculer le patch. +```cpp +void Courbures::get_neighborhood(std::vector &out, const MyMesh::VertexHandle vh) { + for(...) {...} // ajout 1N à out + if(out.size() >= 5) return; + for(...) {...} // ajout 2N à out +} +``` +Après avoir récupéré les voisins, la matrice de changement de base calculé pour positionner les sommets dans un repère local (pour avoir Z représentant la distance au plan XY). Puis on construit la matrice et le vecteur permettant à Eigen de calculer les coefficients du patch quadratique avec la méthode des moindres carrés. Le résultat est retourné dans une structure `QuadPatch` comprenant la matrice de changement de repère, son inverse, ainsi que les coefficients. +```cpp +QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) { + [...] + // Récupération des voisins + get_neighborhood(neigh, vh); + [...] + + [...] // Calcul de la matrice de changement de base + Eigen::Transform ch_base = rot * trans; + + [...] + // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen) + 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 QuadPatch(coef, ch_base); +} +``` +![Visualisation des courbures K](k.png) + +*Visualisation des courbures K* + +## Exercice 2 +Pour la seconde partie, nous avons choisi d'afficher le patch quadratique. +Avec la fonction `fit_quad` chaque sommet possèdes un attribut de type "QuadPatch" qui contient tout ce qu'il faut pour visualiser le patch. + +Pour construire le maillage d'un patch sur un sommet : + - On utilise son 1-voisinage pour déterminer la taille du maillage du patch, en cherchant les coordonnées min et max. + - Puis selon un niveau de discrétisation `patch_divs`, on créé des sommets en calculant leurs positions Z avec les coefficients du patch et en les multipliant par l'inverse de la matrice de changement de base (pour retrouver sa position par rapport au maillage originel). + - Enfin, la dernière étape est d'itérer sur tous ces nouveaux sommets pour ainsi créer des triangles, permettant l'affiche du patch. + +```cpp +MyMesh tesselate_quad_patch(QuadPatch q, MyMesh mesh, VertexHandle vh) { + const size_t patch_divs = 8; + + MyMesh patch; + Eigen::Vector3d point(mesh.point(vh)[0], mesh.point(vh)[1], mesh.point(vh)[2]); + point = q.transform() * point; + + // Recherche des dimensions englobantes du 1-anneau de vh + double xmin = point[0], xmax = point[0]; + double ymin = point[1], ymax = point[1]; + + for (VertexHandle vi : mesh.vv_range(vh)) { + Eigen::Vector3d point(mesh.point(vi)[0], mesh.point(vi)[1], mesh.point(vi)[2]); + + point = q.transform() * point; + xmin = std::min(xmin, point[0]); + xmax = std::max(xmax, point[0]); + ymin = std::min(ymin, point[1]); + ymax = std::max(ymax, point[1]); + } + + // Générations des sommets + double xstep = (xmax-xmin)/static_cast(patch_divs); + double ystep = (ymax-ymin)/static_cast(patch_divs); + + for (size_t y = 0; y < patch_divs; y++) { + for (size_t x = 0; x < patch_divs; x++) { + double dx = xmin + x * xstep; + double dy = ymin + y * ystep; + + Eigen::Vector3d point(dx, dy, -q(dx, dy)); + point = q.inverse_transform() * point; + patch.new_vertex(Point(point[0], point[1], point[2])); + } + } + + // Générations des triangles + for (VertexHandle vhi : patch.vertices()) { + patch.set_color(vhi, MyMesh::Color(0, 1, .2)); + size_t i = vhi.idx(); + size_t x = i % patch_divs; + size_t y = i / patch_divs; + + // On ignore la dernière colonne et dernière ligne + if (x == patch_divs - 1 || y == patch_divs - 1) continue; + + patch.add_face(vhi, patch.vertex_handle(i + patch_divs), patch.vertex_handle(i + 1)); + patch.add_face(patch.vertex_handle(i + 1), patch.vertex_handle(i + patch_divs), + patch.vertex_handle(i + patch_divs + 1)); + } + + return patch; +} +``` + +![Visualisation d'un patch](patch.png) + +*Visualisation d'un patch*