mod_geo-tp/tp.md

117 lines
4.5 KiB
Markdown
Raw Permalink Normal View History

2021-11-24 16:25:13 +01:00
# 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<MyMesh::VertexHandle> &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<double, 3, Eigen::Affine> 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<double>(patch_divs);
double ystep = (ymax-ymin)/static_cast<double>(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*