beginning to believe

src/mesh_processing.cpp

@@ -1,55 +0,0 @@
-#include "mesh_processing.h"
-#include "util.h"
-#include <QGenericMatrix>
-#include <unordered_set>
-
-
-std::vector<HalfedgeHandle> findBorders(MyMesh &mesh) {
-	std::vector<HalfedgeHandle> borders;
-	std::set<HalfedgeHandle> ignore;
-	for (HalfedgeHandle it : mesh.halfedges()) {
-		if (mesh.is_boundary(it)
-			&& ignore.find(it) == ignore.end()) {
-			borders.push_back(it);
-			ignore.insert(it);
-			for (HalfedgeHandle it2 : HalfedgeLoopRange(mesh, it)) {
-				ignore.insert(it2);
-			}
-		}
-	}
-	return borders;
-}
-
-
-void setBorderColor(MyMesh &mesh, const HalfedgeHandle &border,
-					MyMesh::Color color) {
-	for (HalfedgeHandle it : HalfedgeLoopRange(mesh, border)) {
-		mesh.set_color(mesh.to_vertex_handle(it), color);
-	}
-}
-
-
-void fillHoleDumb(MyMesh &mesh, HalfedgeHandle &border) {
-	mesh.request_vertex_status();
-	mesh.request_edge_status();
-	mesh.request_face_status();
-
-	Point center(0, 0, 0);
-	size_t length = 0;
-	for (HalfedgeHandle it : HalfedgeLoopRange(mesh, border)) {
-		VertexHandle vert = mesh.to_vertex_handle(it);
-		center += mesh.point(vert);
-		length++;
-	}
-	center /= length;
-	VertexHandle center_handle = mesh.new_vertex_dirty(center);
-
-	for (HalfedgeHandle it = border;
-		 mesh.to_vertex_handle(it) != center_handle;) {
-		HalfedgeHandle next = mesh.next_halfedge_handle(it);
-		mesh.add_face(mesh.from_vertex_handle(it),
-					  mesh.to_vertex_handle(it),
-					  center_handle);
-		it = next;
-	}
-}

src/mesh_processor.cpp

@@ -1,5 +1,6 @@
 #include "mesh_processor.h"
 #include "util.h"
+#include "smoothing.h"
 #include <QGenericMatrix>
 #include <unordered_set>
 
@@ -13,7 +14,7 @@ static std::vector<std::vector<HalfedgeHandle>> findHoles(MyMesh &mesh) {
 			holes.emplace_back();
 			holes.back().push_back(it);
 			ignore.insert(it);
-			for (HalfedgeHandle it2 : HalfedgeLoopRange(mesh, it)) {
+			for (HalfedgeHandle it2 : HalfedgeLoopRange<MyMesh>(mesh, it)) {
 				holes.back().push_back(it2);
 				ignore.insert(it2);
 			}
@@ -30,8 +31,9 @@ MeshProcessor::MeshProcessor(const QString &path) {
 		qWarning() << "Failed to read" << path;
 		return;
 	}
-	holes = findHoles(mesh);
-	fillHoles();
+	// holes = findHoles(mesh);
+	// fillHoles();
+	smooth(mesh);
 }
 
 

src/mesh_view.cpp

@@ -95,7 +95,6 @@ void MeshView::paint(QOpenGLShaderProgram &program) {
 			/* Wireframe */
 			glf21->glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
 			glf->glDepthMask(GL_FALSE);
-			glf->glLineWidth(2);
 			glf->glDrawElements(GL_TRIANGLES, n_faces * 3, GL_UNSIGNED_INT, 0);
 			glf->glDepthMask(GL_TRUE);
 			glf21->glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

src/mesh_viewer.cpp

@@ -71,7 +71,7 @@ void MeshViewer::resizeGL(int w, int h) {
 
 
 void MeshViewer::paintGL() {
-	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+	// glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 	QMatrix4x4 trans;
 	trans.translate(0, 0, -cam_dist);
 	QMatrix4x4 view = trans * rot;

src/my_mesh.h

@@ -6,6 +6,7 @@
 #include <OpenMesh/Core/Geometry/VectorT.hh>
 #include <QMatrix4x4>
 #include <QColor>
+#include <QVector3D>
 
 
 struct MyTraits : public OpenMesh::DefaultTraits {

src/smoothing.cpp

@@ -0,0 +1,70 @@
+#include "smoothing.h"
+
+#include "OpenMesh/Core/Mesh/SmartHandles.hh"
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
+#include "util.h"
+
+
+
+
+/*
+   	  vα/--\
+   		/ 	 ------\   v
+   	   /  			---X
+   	  /			/---- /	\
+   	 /		/---	 /	 -\
+   	/  /----		/	   \
+ vi	---			   |	 	\
+  	  \--		   /	 	 -\
+  	     \-		  /		 	/--\
+  		   \--	 /	/-------
+   	   	      \-+---
+               vβ
+*/
+
+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);
+	for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
+		Point vi = 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));
+
+		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++;
+	}
+	area /= 3.;
+	// Point tmp(sum / (2*area));
+	// qDebug() << tmp[0] << tmp[1] << tmp[2];
+	return sum / (2.*area);
+	// return sum / count;
+}
+
+
+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);
+			// mesh.set_point(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]);
+		}
+	}
+}

src/smoothing.h

@@ -0,0 +1,11 @@
+#ifndef SMOOTHING_H
+#define SMOOTHING_H
+
+#include "my_mesh.h"
+
+
+Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert);
+void smooth(MyMesh &mesh);
+
+
+#endif

src/util.cpp

@@ -0,0 +1,6 @@
+#include "util.h"
+
+
+QDebug operator<<(QDebug dbg, const Point &p) {
+	return dbg << p[0] << p[1] << p[2];
+}

src/util.h

@@ -1,6 +1,11 @@
 #ifndef UTIL_H
 #define UTIL_H
 
+#include "my_mesh.h"
+
+#include <QVector3D>
+#include <QtMath>
+
 
 template <typename Mesh>
 class HalfedgeLoopRange {
@@ -36,4 +41,33 @@ class ConstHalfedgeLoopRange {
 };
 
 
+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) {
+	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()));
+}
+
+
+template <typename Point>
+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.;
+}
+
+
 #endif