mod_geo-tp/tp.md

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.

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.

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

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.

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