idk some changes, shit's still fucked

CMakeLists.txt

@@ -4,6 +4,7 @@ project(tp LANGUAGES C CXX)
 
 find_package(Qt5 REQUIRED COMPONENTS Core Gui Widgets)
 add_subdirectory(external/OpenMesh)
+add_subdirectory(external/Eigen)
 
 add_executable(${PROJECT_NAME})
 target_sources(${PROJECT_NAME} PRIVATE
@@ -24,7 +25,8 @@ target_sources(${PROJECT_NAME} PRIVATE
   src/util.h)
 target_link_libraries(${PROJECT_NAME} PRIVATE
   Qt5::Core Qt5::Gui Qt5::Widgets
-  OpenMeshCore)
+  OpenMeshCore
+  eigen)
 target_compile_features(${PROJECT_NAME} PRIVATE cxx_std_11)
 set_target_properties(${PROJECT_NAME} PROPERTIES
   AUTOMOC ON

data/bunnyLowPoly.obj

@@ -8128,7 +8128,6 @@ vn 0.8166 -0.5571 -0.1512
 vn 0.4879 -0.8093 -0.3271
 vn 0.4835 -0.8315 -0.2736
 vn 0.7111 -0.7029 -0.0163
-usemtl None
 s 1
 f 1//1 2//1 3//1
 f 4//2 5//2 6//2

external/Eigen/CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.15)
+
+include(FetchContent)
+
+FetchContent_Declare(
+  Eigen
+  URL https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
+  URL_HASH MD5=4c527a9171d71a72a9d4186e65bea559)
+
+FetchContent_MakeAvailable(Eigen)

src/my_mesh.h

@@ -7,14 +7,17 @@
 #include <QMatrix4x4>
 #include <QColor>
 #include <QVector3D>
+#include <OpenMesh/Core/Geometry/EigenVectorT.hh>
 
 
 struct MyTraits : public OpenMesh::DefaultTraits {
+    using Point = Eigen::Vector3<qreal>;
+    using Normal = Eigen::Vector3<qreal>;
+    using Color = Eigen::Vector3<qreal>;
     VertexAttributes(OpenMesh::Attributes::Normal);
     HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
     FaceAttributes(OpenMesh::Attributes::Normal);
     EdgeAttributes(OpenMesh::Attributes::Color);
-	typedef OpenMesh::Vec3f Color;
 };
 
 class MyMesh : public OpenMesh::TriMesh_ArrayKernelT<MyTraits> {

src/smoothing.cpp

@@ -3,6 +3,8 @@
 #include "OpenMesh/Core/Mesh/SmartHandles.hh"
 #include <OpenMesh/Core/Utils/PropertyManager.hh>
 #include "util.h"
+#include <Eigen/Sparse>
+#include <Eigen/Core>
 
 
 
@@ -25,32 +27,21 @@
 Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
 	Point sum {0, 0, 0};
 	qreal area = 0;
-	size_t count = 0;
-	Point v = mesh.point(vert);
+	Point p = mesh.point(vert);
 	for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
-		Point vi = mesh.point(mesh.to_vertex_handle(v_vi));
-
+		Point pi = mesh.point(mesh.to_vertex_handle(v_vi));
 		HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
 		HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
-		Point va = mesh.point(mesh.to_vertex_handle(v_va));
-
+		Point pa = mesh.point(mesh.to_vertex_handle(v_va));
 		HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
-		Point vb = mesh.point(mesh.to_vertex_handle(vi_vb));
-
-		qreal a = angle_between(vi, va, v);
-		qreal b = angle_between(v, vb, vi);
-		sum += (1/qTan(a) + 1/qTan(b)) * (vi - v);
-		// sum += (vi - v);
-		area += triangle_area(v, vi, vb);
-		// qDebug("a: %lf, b: %lf, sum: %f %f %f, area: %lf",
-		// 	   a, b, sum[0], sum[1], sum[2], area);
-		count++;
+		Point pb = mesh.point(mesh.to_vertex_handle(vi_vb));
+		qreal a = angle_between(pi, pa, p);
+		qreal b = angle_between(pi, pb, p);
+		sum += (cotan(a) + cotan(b)) * (pi - p);
+		area += triangle_area(p, pi, pb);
 	}
 	area /= 3.;
-	// Point tmp(sum / (2*area));
-	// qDebug() << tmp[0] << tmp[1] << tmp[2];
 	return sum / (2.*area);
-	// return sum / count;
 }
 
 
@@ -59,7 +50,6 @@ void smooth(MyMesh &mesh) {
 	for (VertexHandle v : mesh.vertices()) {
 		if (!mesh.is_boundary(v)) {
 			new_pos[v] = mesh.point(v) + laplace_beltrami(mesh, v);
-			// mesh.set_point(v, mesh.point(v) + laplace_beltrami(mesh, v));
 		}
 	}
 	for (VertexHandle v : mesh.vertices()) {
@@ -67,4 +57,54 @@ void smooth(MyMesh &mesh) {
 			mesh.set_point(v, new_pos[v]);
 		}
 	}
+
+	// // Approche matricielle
+	// size_t n_verts = mesh.n_vertices();
+	// Eigen::SparseMatrix<qreal> mass(n_verts, n_verts);
+	// Eigen::SparseMatrix<qreal> cotangent(n_verts, n_verts);
+	// for (VertexHandle vert : mesh.vertices()) {
+	// 	if (mesh.is_boundary(vert)) continue;
+	// 	qreal sum = 0;
+	// 	qreal area = 0;
+	// 	qreal count = 0;
+	// 	Point v = mesh.point(vert);
+	// 	for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
+	// 		VertexHandle vi_h = mesh.to_vertex_handle(v_vi);
+	// 		Point vi = mesh.point(vi_h);
+
+	// 		HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
+	// 		HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
+	// 		Point va = mesh.point(mesh.to_vertex_handle(v_va));
+
+	// 		HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
+	// 		Point vb = mesh.point(mesh.to_vertex_handle(vi_vb));
+	// 		qreal a = angle_between(vi, va, v);
+	// 		qreal b = angle_between(v, vb, vi);
+	// 		qreal w = -(1./qTan(a) + 1./qTan(b)) / 2.;
+	// 		sum += w;
+	// 		cotangent.insert(vert.idx(), vi_h.idx()) = w;
+	// 		area += triangle_area(v, vi, vb);
+	// 		count++;
+	// 	}
+	// 	cotangent.insert(vert.idx(), vert.idx()) = -sum;
+	// 	mass.insert(vert.idx(), vert.idx()) = 1. / area;
+	// }
+	// Eigen::SparseMatrix<qreal> laplacian = mass * cotangent;
+	// 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();
+	// 	mesh.set_point(vert, {X[id], Y[id], Z[id]});
+	// }
 }

src/util.cpp

@@ -4,3 +4,8 @@
 QDebug operator<<(QDebug dbg, const Point &p) {
 	return dbg << p[0] << p[1] << p[2];
 }
+
+
+qreal cotan(const qreal x) {
+	return 1. / qTan(x);
+}

src/util.h

@@ -47,27 +47,22 @@ QDebug operator<<(QDebug dbg, const Point &p);
 /* Returns the angle between the vector from p2 to p1 and the vector
  * from p2 to p3. */
 template <typename Point>
-qreal angle_between(Point p1,
-					Point p2,
-					Point p3) {
+qreal angle_between(Point p1, Point p2, Point p3) {
 	Point vec_a = p1 - p2;
 	Point vec_b = p3 - p2;
-	QVector3D a(vec_a[0], vec_a[1], vec_a[2]);
-	QVector3D b(vec_b[0], vec_b[1], vec_b[2]);
-	return qAcos(QVector3D::dotProduct(a.normalized(), b.normalized()));
+	return qAcos(vec_a.dot(vec_b) / (vec_a.norm() * vec_b.norm()));
 }
 
 
 template <typename Point>
-qreal triangle_area(Point p1,
-					Point p2,
-					Point p3) {
+qreal triangle_area(Point p1, Point p2, Point p3) {
 	Point vec_a = p1 - p2;
 	Point vec_b = p3 - p2;
-	QVector3D a(vec_a[0], vec_a[1], vec_a[2]);
-	QVector3D b(vec_b[0], vec_b[1], vec_b[2]);
-	return QVector3D::crossProduct(a, b).length() / 2.;
+	return vec_a.cross(vec_b).norm() / 2.;
 }
 
 
+qreal cotan(const qreal x);
+
+
 #endif