improve the implicit surface approximation

by scaling the normals based on the mesh size

CMakeLists.txt

@@ -36,7 +36,9 @@ target_sources(${PROJECT_NAME} PRIVATE
   src/hole_filling.cpp
   src/hole_filling.h
   src/double_input.cpp
-  src/double_input.h)
+  src/double_input.h
+  src/noise_removal.cpp
+  src/noise_removal.h)
 target_link_libraries(${PROJECT_NAME} PRIVATE
   Qt5::Core Qt5::Gui Qt5::Widgets
   OpenMeshCore

LISEZ-MOI.txt

@@ -6,14 +6,15 @@ Ceci comprend les fonctionalités suivantes :
  - analyse de courbure (dans ‘src/curvature.cpp’)
  - adoucissement d'un maillage par laplacien uniforme ou cotangent (dans
    ‘src/smoothing.cpp’)
+ - suppression de parasites (dans ‘src/noise_removal.cpp’)
 
 Des maillages exemples sont inclus, ‘data/gargoyle_trou.obj’ est un
 bon exemple pour le remplissage, ‘data/bunnyLowPoly-noisy.obj’ pour
-l'adoucissement.
+l'adoucissement et ‘data/bunny.obj’ pour la suppression de parasites.
 
 
 Compilation
-   ‘cmake -Bbuild -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS=-O2’
+   ‘cmake -Bbuild -DCMAKE_BUILD_TYPE=Release’
    puis
    ‘cmake --build build’
 

external/MeshReconstruction/lib/MeshReconstruction.cpp

@@ -1,6 +1,7 @@
 #include "MeshReconstruction.h"
 #include "Cube.h"
 #include "Triangulation.h"
+#include <iostream>
 
 using namespace MeshReconstruction;
 using namespace std;
@@ -44,8 +45,9 @@ Mesh MeshReconstruction::MarchCube(
 	auto const NumY = static_cast<int>(ceil(domain.size.y / cubeSize.y));
 	auto const NumZ = static_cast<int>(ceil(domain.size.z / cubeSize.z));
 
-	auto const HalfCubeDiag = cubeSize.Norm() / 2.0;
-	auto const HalfCubeSize = cubeSize * 0.5;
+	auto const CubeDiag = cubeSize.Norm();
+	// auto const HalfCubeDiag = cubeSize.Norm() * 0.5;
+	// auto const HalfCubeSize = cubeSize * 0.5;
 
 	Mesh mesh;
 
@@ -61,9 +63,10 @@ Mesh MeshReconstruction::MarchCube(
 				Vec3 min{ x, y, z };
 
 				// Process only if cube lies within narrow band around surface.
-				auto cubeCenter = min + HalfCubeSize;
-				double dist = std::fabs(sdf(cubeCenter) - isoLevel);
-				if (dist > HalfCubeDiag) continue;
+				// auto cubeCenter = min + HalfCubeSize;
+				// double dist = std::abs(sdf(cubeCenter) - isoLevel);
+				double dist = std::abs(sdf(min) - isoLevel);
+				// if (dist > CubeDiag) continue;
 
 				Cube cube({ min, cubeSize }, sdf);
 				auto intersect = cube.Intersect(isoLevel);

src/double_input.cpp

@@ -46,6 +46,7 @@ DoubleInput::DoubleInput(QObject *parent, double min, double max,
 	_spin_box->setValue(value);
 	_slider->setMaximum(_slider_resolution);
 	_slider->setValue(doubleToInt(value));
+	_slider->setTracking(false);
 	connect(_slider, &QSlider::valueChanged,
 			this, &DoubleInput::onSliderValueChanged);
 	connect(_spin_box, QOverload<double>::of(&QDoubleSpinBox::valueChanged),

src/hole_filling.cpp

@@ -2,6 +2,7 @@
 #include "IO.h"
 #include "util.h"
 #include <MeshReconstruction.h>
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
 
 
 static std::vector<std::vector<HalfedgeHandle>> findHoles(MyMesh &mesh) {
@@ -153,7 +154,7 @@ vector<MyMesh::VertexHandle> Hole_Filling::next_neighbors(const vector<MyMesh::V
 
 // ***** Computation of RBF
 
-pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist)
+pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist, double normal_scale)
 {
     const int n(vlist.size()), d(10) ;
     Eigen::MatrixXd & A = *(new Eigen::MatrixXd(3*n+d,3*n+d)) ;
@@ -170,12 +171,12 @@ pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_F
     //Append vertices+normals to pts_list
     for (int i=0; i<n; i++)
     {
-        pts_list.push_back(_mesh.point(vlist.at(i)) + _mesh.normal(vlist.at(i))) ;
+        pts_list.push_back(_mesh.point(vlist.at(i)) + _mesh.normal(vlist.at(i)) * normal_scale) ;
     }
     //Append vertices-normals to pts_list
     for (int i=0; i<n; i++)
     {
-        pts_list.push_back(_mesh.point(vlist.at(i)) - _mesh.normal(vlist.at(i))) ;
+        pts_list.push_back(_mesh.point(vlist.at(i)) - _mesh.normal(vlist.at(i)) * normal_scale) ;
     }
 
     int nn = pts_list.size() ;
@@ -302,18 +303,40 @@ Mesh Hole_Filling::poly_n_out(const Implicit_RBF &implicit, Rect3 domain)
 }
 
 
+/* Computes a mesh's bounding box and stores it in a mesh property
+ * named "bounding_box". */
+static void computeMeshBoundingBox(MyMesh &mesh, Hole_Filling &hf) {
+	try {
+		auto mesh_bb = OpenMesh::getProperty<void, Rect3>
+			(mesh, "bounding_box");
+	} catch (const std::runtime_error &e) {
+		auto mesh_bb = OpenMesh::getOrMakeProperty<void, Rect3>
+			(mesh, "bounding_box");
+		std::vector<VertexHandle> verts;
+		for (VertexHandle vh : mesh.vertices()) {
+			verts.push_back(vh);
+		}
+		*mesh_bb = hf.estimate_BB(verts);
+	}
+}
+
+
 MyMesh fillHoleImplicit(MyMesh &mesh, Hole_Filling &hf,
 						std::vector<HalfedgeHandle> &hole) {
+	computeMeshBoundingBox(mesh, hf);
+	Rect3 mesh_bb = *OpenMesh::getProperty<void, Rect3>(mesh, "bounding_box");
+	double diag = mesh_bb.size.Norm();
+
 	std::vector<VertexHandle> verts;
 	for (HalfedgeHandle hh : hole) {
 		verts.push_back(mesh.to_vertex_handle(hh));
 	}
+	auto bb = hf.estimate_BB(verts) ;
 	verts = hf.next_neighbors(verts);
-	auto [system, pts_list] = hf.compute_approx_mat(verts);
+	auto [system, pts_list] = hf.compute_approx_mat(verts, diag * .1);
 	auto [alpha, beta] = hf.solve_approx(system, pts_list.size(), 10);
 	Implicit_RBF rbf(alpha, beta, pts_list);
 
-	auto bb = hf.estimate_BB(verts) ;
 	Mesh filling = 	hf.poly_n_out(rbf, bb);
 	MyMesh ret;
 	for (const Vec3 &v : filling.vertices) {
@@ -331,27 +354,11 @@ MyMesh fillHoleImplicit(MyMesh &mesh, Hole_Filling &hf,
 
 std::vector<MyMesh> fillHolesImplicit(MyMesh &mesh,
 									  float scale, float discr) {
-	Hole_Filling hf(mesh, scale, discr);
 	mesh.holes = findHoles(mesh);
 	std::vector<MyMesh> fillings;
 	for (auto hole : mesh.holes) {
+		Hole_Filling hf(mesh, scale, discr);
 		fillings.push_back(fillHoleImplicit(mesh, hf, hole));
 	}
 	return fillings;
-
-	// auto sdf = [&](Vec3 const& v) { return v.Norm() - 10.; };
-	// Rect3 domain {{-10, -10, -10}, {20, 20, 20}};
-	// auto filling = MarchCube(sdf, domain);
-	// WriteObjFile(filling, "out.obj");
-	// MyMesh ret;
-	// for (const Vec3 &v : filling.vertices) {
-	// 	VertexHandle vh = ret.new_vertex({v.x, v.y, v.z});
-	// 	ret.set_color(vh, ret.default_color);
-	// }
-	// for (const Triangle &t : filling.triangles) {
-	// 	ret.add_face(ret.vertex_handle(t[0]),
-	// 				 ret.vertex_handle(t[1]),
-	// 				 ret.vertex_handle(t[2]));
-	// }
-	// return {ret};
 }

src/hole_filling.h

@@ -103,7 +103,7 @@ public:
     vector<MyMesh::VertexHandle> next_neighbors(const vector<MyMesh::VertexHandle> & bnd) ;
 
     // Computation of RBF
-    pair<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<MyMesh::Point> &> compute_approx_mat(vector<MyMesh::VertexHandle> vlist) ;
+    pair<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<MyMesh::Point> &> compute_approx_mat(vector<MyMesh::VertexHandle> vlist, double normal_scale=1) ;
     pair<vector<float>&, vector<float>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;
 
     // IO

src/main.cpp

@@ -25,6 +25,8 @@ static MeshProcessor *create_mesh_processor(const QString &path,
 					 mesh_processor, &MeshProcessor::setImplicitHoleFillingScale);
 	QObject::connect(&main_window, &MainWindow::fillHolesImplicitDiscrChanged,
 					 mesh_processor, &MeshProcessor::setImplicitHoleFillingDiscr);
+	QObject::connect(&main_window, &MainWindow::filterNoiseClicked,
+					 mesh_processor, &MeshProcessor::removeNoise);
 	return mesh_processor;
 }
 

src/main_window.cpp

@@ -9,6 +9,7 @@
 #include <QPushButton>
 #include <QSlider>
 #include <QLabel>
+#include <qnamespace.h>
 
 
 MainWindow::MainWindow(QWidget *parent)
@@ -113,4 +114,19 @@ MainWindow::MainWindow(QWidget *parent)
 			this, &MainWindow::patchViewToggled);
 	curvature_layout->addWidget(patch_mode);
 	toolbar.addWidget(curvature_box);
+
+
+	// Noise cleaning
+	QGroupBox *noise_box = new QGroupBox("Élagage");
+	QLayout *noise_layout = new QVBoxLayout();
+	noise_box->setLayout(noise_layout);
+	QSlider *noise_slider = new QSlider(Qt::Horizontal);
+	QPushButton *noise_button = new QPushButton("Élaguer");
+	connect(noise_button, &QPushButton::clicked,
+			[=]() {
+				emit filterNoiseClicked(noise_slider->value());
+			});
+	noise_layout->addWidget(noise_slider);
+	noise_layout->addWidget(noise_button);
+	toolbar.addWidget(noise_box);
 }

src/main_window.h

@@ -32,6 +32,7 @@ signals:
 	void smoothUniformClicked();
 	void smoothCotangentClicked(double factor);
 	void patchViewToggled(bool checked);
+	void filterNoiseClicked(int value);
 
 public:
 	MeshViewer mesh_viewer;

src/mesh_processor.cpp

@@ -4,6 +4,7 @@
 #include "hole_filling.h"
 #include "smoothing.h"
 #include "curvature.h"
+#include "noise_removal.h"
 #include <QGenericMatrix>
 #include <unordered_set>
 
@@ -31,6 +32,7 @@ MeshProcessor::MeshProcessor(const QString &path, MeshViewer &mesh_viewer,
 		qWarning() << "Curvature computation failed";
 	}
 	connect(&mesh_viewer, &MeshViewer::clicked, this, &MeshProcessor::click);
+
 	updateView();
 }
 
@@ -43,7 +45,7 @@ MeshProcessor::~MeshProcessor() {
 }
 
 
-void MeshProcessor::updateView() const {
+void MeshProcessor::updateView() {
 	if (mesh_viewer.isInitialized()) {
 		mesh_viewer.removeMesh(mesh);
 		mesh_viewer.addMesh(mesh);
@@ -131,3 +133,12 @@ void MeshProcessor::click(QVector3D position) {
 	mesh_viewer.addMesh(patch);
 	mesh_viewer.updateForReal();
 }
+
+
+void MeshProcessor::removeNoise(int value) {
+	std::cout << value << std::endl;
+	::remove_noise(mesh, value);
+	mesh_viewer.removeMesh(mesh);
+	mesh_viewer.addMesh(mesh);
+	mesh_viewer.updateForReal();
+}

src/mesh_processor.h

@@ -18,7 +18,7 @@ class MeshProcessor : public QObject {
 	double implicit_hole_filling_discr;
 	std::vector<MyMesh> fillings;
 
-	void updateView() const;
+	void updateView();
 
 public:
 	MyMesh mesh;
@@ -38,6 +38,7 @@ public slots:
 	void smoothUniform();
 	void setPatchView(bool on);
 	void click(QVector3D position);
+	void removeNoise(int value);
 };
 
 

src/mesh_viewer.cpp

@@ -89,9 +89,10 @@ void MeshViewer::paintGL() {
 }
 
 
-void MeshViewer::addMesh(const MyMesh &mesh) {
+void MeshViewer::addMesh(MyMesh &mesh) {
 	Q_ASSERT(isValid());
 	makeCurrent();
+	mesh.viewer_id = meshes.size();
 	meshes.emplace_back(mesh, program);
 	doneCurrent();
 	update();
@@ -100,9 +101,14 @@ void MeshViewer::addMesh(const MyMesh &mesh) {
 
 void MeshViewer::removeMesh(const MyMesh &mesh) {
 	makeCurrent();
-	meshes.remove_if([&](const MeshView &mv) {
-		return &mv.mesh == &mesh;
-	});
+	size_t i = 0;
+	for (auto it = meshes.begin(); it != meshes.end(); ++it) {
+		if (i == mesh.viewer_id) {
+			meshes.erase(it);
+			break;
+		}
+		i++;
+	}
 	doneCurrent();
 	update();
 }

src/mesh_viewer.h

@@ -15,7 +15,7 @@
 
 
 #define WIREFRAME_COLOR 0, 0, 0
-#define FOV 70
+#define FOV 40
 
 
 using namespace OpenMesh;
@@ -46,7 +46,7 @@ public:
 	constexpr bool isInitialized() { return is_initialized; }
 
 public slots:
-	void addMesh(const MyMesh &mesh);
+	void addMesh(MyMesh &mesh);
 	void removeMesh(const MyMesh &mesh);
 	void updateForReal();
 

src/my_mesh.h

@@ -14,9 +14,10 @@ struct MyTraits : public OpenMesh::DefaultTraits {
     using Point = Eigen::Vector3<qreal>;
     using Normal = Eigen::Vector3<qreal>;
     using Color = Eigen::Vector3<qreal>;
-    VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
+    VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color | OpenMesh::Attributes::Status);
+	EdgeAttributes(OpenMesh::Attributes::Status);
     HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
-    FaceAttributes(OpenMesh::Attributes::Normal);
+    FaceAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Status);
     // EdgeAttributes(OpenMesh::Attributes::Color);
 };
 
@@ -24,6 +25,7 @@ class MyMesh : public OpenMesh::TriMesh_ArrayKernelT<MyTraits> {
 public:
 	Color default_color {.5, .5, .5};
 	QMatrix4x4 transform;
+	size_t viewer_id;
 	std::vector<std::vector<HalfedgeHandle>> holes;
 };
 

src/noise_removal.cpp

@@ -0,0 +1,15 @@
+#include "noise_removal.h"
+#include "util.h"
+
+
+void remove_noise(MyMesh &mesh, unsigned threshold) {
+	auto components = find_connected_components(mesh);
+	for (auto component : components) {
+		if (component.size() < threshold) {
+			for (VertexHandle vh : component) {
+				mesh.delete_vertex(vh);
+			}
+		}
+	}
+	mesh.garbage_collection();
+}

src/noise_removal.h

@@ -0,0 +1,10 @@
+#ifndef NOISE_REMOVAL_H
+#define NOISE_REMOVAL_H
+
+#include "my_mesh.h"
+
+
+void remove_noise(MyMesh &mesh, unsigned threshold=20);
+
+
+#endif

src/smoothing.cpp

@@ -132,7 +132,7 @@ void smooth(MyMesh &mesh, SmoothingMethod method, double cotan_factor) {
 		laplacian = laplacian_matrix(mesh, uniform_weight, uniform_mass);
 	}
 
-	laplacian = laplacian * laplacian; // Mise au carré.
+	// laplacian = laplacian * laplacian; // Mise au carré.
 
 	// Transfert des coordonées de chaque point dans une matrice.
 	size_t n_verts = mesh.n_vertices();

src/util.cpp

@@ -1,5 +1,8 @@
 #include "util.h"
 
+#include <stack>
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
+
 
 QDebug operator<<(QDebug dbg, const Point &p) {
 	return dbg << p[0] << p[1] << p[2];
@@ -9,3 +12,44 @@ QDebug operator<<(QDebug dbg, const Point &p) {
 qreal cotan(const qreal x) {
 	return 1. / qTan(x);
 }
+
+
+static void traverse_connected(MyMesh &mesh, std::vector<VertexHandle> &out,
+							   VertexHandle start) {
+	auto prop =
+		OpenMesh::makeTemporaryProperty<MyMesh::VertexHandle, bool>(mesh);
+	for (VertexHandle vh : mesh.vertices()) {
+		prop[vh] = false;
+	}
+	std::stack<VertexHandle> stack;
+	stack.push(start);
+	prop[start] = true;
+	while (!stack.empty()) {
+		VertexHandle v = stack.top();
+		stack.pop();
+		out.emplace_back(v);
+		for (VertexHandle v2 : mesh.vv_range(v)) {
+			if (!prop[v2]) {
+				stack.push(v2);
+				prop[v2] = true;
+			}
+		}
+	}
+}
+
+
+std::vector<std::vector<VertexHandle>> find_connected_components(
+	MyMesh &mesh) {
+	auto prop =
+		OpenMesh::makeTemporaryProperty<MyMesh::VertexHandle, bool>(mesh);
+	for (VertexHandle vh : mesh.vertices()) {
+		prop[vh] = false;
+	}
+	std::vector<std::vector<VertexHandle>> connected_components;
+	for (VertexHandle vh : mesh.vertices()) {
+		if (prop[vh]) continue;
+		connected_components.push_back({});
+		traverse_connected(mesh, connected_components.back(), vh);
+	}
+	return connected_components;
+}

src/util.h

@@ -101,4 +101,8 @@ standard_deviation(const ForwardIt first,
 }
 
 
+std::vector<std::vector<VertexHandle>> find_connected_components(
+	MyMesh &mesh);
+
+
 #endif