# 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*