From 068dede83280bd7e750675fefc35d84c7a78949d Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Thu, 25 Nov 2021 14:59:52 +0100
Subject: [PATCH] smoothing tweaks

---
 src/smoothing.cpp | 75 ++++++++++++++++++++++++-----------------------
 1 file changed, 39 insertions(+), 36 deletions(-)

diff --git a/src/smoothing.cpp b/src/smoothing.cpp
index ec03dec..4d25f9e 100644
--- a/src/smoothing.cpp
+++ b/src/smoothing.cpp
@@ -42,8 +42,9 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
 		area += triangle_area(p, pi, pb);
 		count++;
 	}
-	area /= 3.;
-	return sum / (2.*count);
+	// area /= 3.;
+	// return sum / (2.*area);
+	return sum / count;
 }
 
 
@@ -73,46 +74,48 @@ Eigen::SparseMatrix<qreal> laplacian_matrix(const MyMesh &mesh) {
 			area += triangle_area(p, pi, pb);
 			count++;
 		}
-		area /= 3.;
+		// area /= 3.;
 		cotangent.insert(vert.idx(), vert.idx()) = -sum;
-		mass.insert(vert.idx(), vert.idx()) = 1. / area;
+		mass.insert(vert.idx(), vert.idx()) = 1. / (4. * area);
+		// mass.insert(vert.idx(), vert.idx()) = 1. / count;
 	}
 	return mass * cotangent;
 }
 
 
 void smooth(MyMesh &mesh) {
-	// auto new_pos = OpenMesh::makeTemporaryProperty<VertexHandle, Point>(mesh);
-	// for (VertexHandle v : mesh.vertices()) {
-	// 	if (!mesh.is_boundary(v)) {
-	// 		new_pos[v] = mesh.point(v) - laplace_beltrami(mesh, v);
-	// 	}
-	// }
-	// for (VertexHandle v : mesh.vertices()) {
-	// 	if (!mesh.is_boundary(v)) {
-	// 		mesh.set_point(v, new_pos[v]);
-	// 	}
-	// }
+	auto new_pos = OpenMesh::makeTemporaryProperty<VertexHandle, Point>(mesh);
+	for (VertexHandle v : mesh.vertices()) {
+		if (!mesh.is_boundary(v)) {
+			new_pos[v] = mesh.point(v) - laplace_beltrami(mesh, v);
+		}
+	}
+	for (VertexHandle v : mesh.vertices()) {
+		if (!mesh.is_boundary(v)) {
+			mesh.set_point(v, new_pos[v]);
+		}
+	}
 
-	// Approche matricielle
-	Eigen::SparseMatrix<qreal> laplacian = laplacian_matrix(mesh);
-	size_t n_verts = mesh.n_vertices();
-	Eigen::VectorX<qreal> X(n_verts), Y(n_verts), Z(n_verts);
-	for (VertexHandle vert : mesh.vertices()) {
-		if (mesh.is_boundary(vert)) continue;
-		size_t id = vert.idx();
-		Point p = mesh.point(vert);
-		X(id) = p[0];
-		Y(id) = p[1];
-		Z(id) = p[2];
-	}
-	X = laplacian * X;
-	Y = laplacian * Y;
-	Z = laplacian * Z;
-	for (VertexHandle vert : mesh.vertices()) {
-		if (mesh.is_boundary(vert)) continue;
-		size_t id = vert.idx();
-		Point offset {X(id), Y(id), Z(id)};
-		mesh.set_point(vert, mesh.point(vert) - offset);
-	}
+	// // Approche matricielle
+	// Eigen::SparseMatrix<qreal> laplacian = laplacian_matrix(mesh);
+	// // laplacian = laplacian * laplacian;
+	// size_t n_verts = mesh.n_vertices();
+	// Eigen::VectorX<qreal> X(n_verts), Y(n_verts), Z(n_verts);
+	// for (VertexHandle vert : mesh.vertices()) {
+	// 	if (mesh.is_boundary(vert)) continue;
+	// 	size_t id = vert.idx();
+	// 	Point p = mesh.point(vert);
+	// 	X(id) = p[0];
+	// 	Y(id) = p[1];
+	// 	Z(id) = p[2];
+	// }
+	// X = laplacian * X;
+	// Y = laplacian * Y;
+	// Z = laplacian * Z;
+	// for (VertexHandle vert : mesh.vertices()) {
+	// 	if (mesh.is_boundary(vert)) continue;
+	// 	size_t id = vert.idx();
+	// 	Point offset {X(id), Y(id), Z(id)};
+	// 	mesh.set_point(vert, mesh.point(vert) - offset);
+	// }
 }