improve the implicit surface approximation

by scaling the normals based on the mesh size

CMakeLists.txt

@@ -6,6 +6,10 @@ find_package(Qt5 REQUIRED COMPONENTS Core Gui Widgets)
 add_subdirectory(external/OpenMesh)
 add_subdirectory(external/Eigen)
 
+add_subdirectory(external/MeshReconstruction)
+target_include_directories(MeshReconstruction INTERFACE
+  external/MeshReconstruction/lib)
+
 add_executable(${PROJECT_NAME})
 target_sources(${PROJECT_NAME} PRIVATE
   resources.qrc
@@ -28,11 +32,18 @@ target_sources(${PROJECT_NAME} PRIVATE
   src/quad_patch.cpp
   src/quad_patch.h
   src/quad_patch_tesselator.cpp
-  src/quad_patch_tesselator.h)
+  src/quad_patch_tesselator.h
+  src/hole_filling.cpp
+  src/hole_filling.h
+  src/double_input.cpp
+  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
-  eigen)
+  eigen
+  MeshReconstruction)
 target_compile_features(${PROJECT_NAME} PRIVATE cxx_std_17)
 set_target_properties(${PROJECT_NAME} PROPERTIES
   AUTOMOC ON

LISEZ-MOI.txt

@@ -0,0 +1,29 @@
+TP de modélisation géométrique réalisé dans le cadre du M2 GIG
+par Jérémy André et Cyril Colin
+
+Ceci comprend les fonctionalités suivantes :
+ - remplissage de trous (dans ‘src/hole_filling.cpp’)
+ - 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 et ‘data/bunny.obj’ pour la suppression de parasites.
+
+
+Compilation
+   ‘cmake -Bbuild -DCMAKE_BUILD_TYPE=Release’
+   puis
+   ‘cmake --build build’
+
+
+Exécution
+   ‘build/tp [FICHIER OBJ]’
+
+
+Dépendances
+   - Eigen (téléchargée automatiquement)
+   - Openmesh (téléchargée automatiquement)
+   - MeshReconstruction (inclue dans `external')

external/MeshReconstruction/.gitignore

@@ -0,0 +1,3 @@
+bin
+html/
+latex/

external/MeshReconstruction/CMakeLists.txt

@@ -0,0 +1,16 @@
+cmake_minimum_required (VERSION 3.0)
+project (MeshReconstruction)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Parallel compilation.
+if(WIN32)
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP")
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250")
+endif()
+
+set(CMAKE_DEBUG_POSTFIX "d")
+
+add_subdirectory(lib)
+add_subdirectory(demo)

external/MeshReconstruction/Readme.md

@@ -0,0 +1,13 @@
+# MeshReconstruction
+
+This is a small library that can reconstruct a triangle mesh from a <a href="http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm">signed distance function</a> using the <a href="https://en.wikipedia.org/wiki/Marching_cubes">Marching Cubes algorithm</a> and export it to a file in the <a href="https://de.wikipedia.org/wiki/Wavefront_OBJ">obj format</a>.
+
+The library is self-contained and has no dependencies. The library is fast due to precomputed lookup tables and a narrow-band approach which excludes a lot of marching cubes that are far away from the surface.
+
+The library requires C++14 and has been tested under Visual Studio 2017 and Windows 10 but should port to other systems without major problems.
+
+The library can be used under the terms of the MIT License.
+
+<p align="center">
+  <img src="overview.png" width="400" alt="Overview">
+</p>

external/MeshReconstruction/demo/CMakeLists.txt

@@ -0,0 +1,3 @@
+include_directories("." ${PROJECT_SOURCE_DIR}/lib)
+add_executable(Demo main.cpp)
+target_link_libraries (Demo MeshReconstruction)

external/MeshReconstruction/demo/main.cpp

@@ -0,0 +1,22 @@
+#include <MeshReconstruction.h>
+#include <IO.h>
+
+using namespace MeshReconstruction;
+
+int main()
+{
+	auto sphereSdf = [](Vec3 const& pos)
+	{
+		auto const Radius = 1.0;
+		return pos.Norm() - Radius;
+	};
+
+	Rect3 domain;
+	domain.min = { -1.1, -1.1, -1.1 };
+	domain.size = {2.2, 2.2, 2.2};
+
+	Vec3 cubeSize{ 0.05, 0.05, 0.05 };
+
+	auto mesh = MarchCube(sphereSdf, domain, cubeSize);
+	WriteObjFile(mesh, "sphere.obj");
+}

external/MeshReconstruction/lib/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_library(MeshReconstruction MeshReconstruction.h MeshReconstruction.cpp Cube.h Cube.cpp DataStructs.h IO.h IO.cpp Triangulation.h Triangulation.cpp)
+
+install(TARGETS MeshReconstruction
+        RUNTIME DESTINATION bin
+        LIBRARY DESTINATION lib
+        ARCHIVE DESTINATION lib)
+
+install(DIRECTORY .
+        DESTINATION .
+        FILES_MATCHING PATTERN "*.h")

external/MeshReconstruction/lib/Cube.cpp

@@ -0,0 +1,148 @@
+#include "Cube.h"
+
+using namespace MeshReconstruction;
+
+namespace
+{
+	// Cube has 8 vertices. Each vertex can have positive or negative sign.
+	// 2^8 = 256 possible configurations mapped to intersected edges in each case.
+	// The 12 edges are numbered as 1, 2, 4, ..., 2048 and are stored as a 12-bit bitstring for each configuration.
+	const int signConfigToIntersectedEdges[256] = {
+		0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
+		0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+		0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
+		0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+		0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
+		0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+		0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
+		0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+		0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
+		0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+		0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
+		0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+		0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
+		0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+		0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
+		0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+		0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
+		0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+		0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
+		0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+		0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
+		0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+		0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
+		0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
+		0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
+		0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
+		0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
+		0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
+		0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
+		0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
+		0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
+		0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 };
+
+	struct Edge
+	{
+		unsigned edgeFlag : 12; // flag: 1, 2, 4, ... 2048
+		int vert0; // 0-7
+		int vert1; // 0-7
+	};
+
+	const Edge edges[12] =
+	{
+		{ 1, 0, 1 }, // edge 0
+		{ 2, 1, 2 }, // edge 1
+		{ 4, 2, 3 }, // ...
+		{ 8, 3, 0 },
+		{ 16, 4, 5 },
+		{ 32, 5, 6 },
+		{ 64, 6, 7 },
+		{ 128, 7, 4 },
+		{ 256, 0, 4 },
+		{ 512, 1, 5 },
+		{ 1024, 2, 6 },
+		{ 2048, 3, 7 } // edge 11
+	};
+}
+
+Vec3 Cube::LerpVertex(double isoLevel, int i1, int i2) const
+{
+	auto const Eps = 1e-5;
+	auto const v1 = sdf[i1];
+	auto const v2 = sdf[i2];
+	auto const& p1 = pos[i1];
+	auto const& p2 = pos[i2];
+
+	if (std::abs(isoLevel - v1) < Eps) return p1;
+	if (std::abs(isoLevel - v2) < Eps) return p2;
+	if (std::abs(v1 - v2) < Eps) return p1;
+
+	auto mu = (isoLevel - v1) / (v2 - v1);
+	return p1 + (p2 - p1)*mu;
+}
+
+Cube::Cube(Rect3 const& space, Fun3s const& sdf)
+{
+	auto mx = space.min.x;
+	auto my = space.min.y;
+	auto mz = space.min.z;
+
+	auto sx = space.size.x;
+	auto sy = space.size.y;
+	auto sz = space.size.z;
+
+	pos[0] = space.min;
+	pos[1] = { mx + sx, my, mz };
+	pos[2] = { mx + sx, my, mz + sz };
+	pos[3] = { mx, my, mz + sz };
+	pos[4] = { mx, my + sy, mz };
+	pos[5] = { mx + sx, my + sy, mz };
+	pos[6] = { mx + sx, my + sy, mz + sz };
+	pos[7] = { mx, my + sy, mz + sz };
+
+	for (auto i = 0; i < 8; ++i)
+	{
+		auto sd = sdf(pos[i]);
+		if (sd == 0) sd += 1e-6;
+		this->sdf[i] = sd;
+	}
+}
+
+int Cube::SignConfig(double isoLevel) const
+{
+	auto edgeIndex = 0;
+
+	for (auto i = 0; i < 8; ++i)
+	{
+		if (sdf[i] < isoLevel)
+		{
+			edgeIndex |= (1 << i);
+		}
+	}
+
+	return edgeIndex;
+}
+
+IntersectInfo Cube::Intersect(double iso) const
+{
+	// idea:
+	// from signs at 8 corners of cube a sign configuration (256 possible ones) is computed
+	// this configuration can be used to index into a table that tells which of the 12 edges are intersected
+	// find vertices adjacent to edges and interpolate cut vertex and store it in IntersectionInfo object
+
+	IntersectInfo intersect;
+	intersect.signConfig = SignConfig(iso);
+	auto intersectedEdges = signConfigToIntersectedEdges[intersect.signConfig];
+	for (auto e = 0; e<12; ++e)
+	{
+		if (intersectedEdges & edges[e].edgeFlag)
+		{
+			auto v0 = edges[e].vert0;
+			auto v1 = edges[e].vert1;
+			auto vert = LerpVertex(iso, v0, v1);
+			intersect.edgeVertIndices[e] = vert;
+		}
+	}
+
+	return intersect;
+}

external/MeshReconstruction/lib/Cube.h

@@ -0,0 +1,30 @@
+#pragma once
+#include "DataStructs.h"
+
+namespace MeshReconstruction
+{
+	struct IntersectInfo
+	{
+		// 0 - 255
+		int signConfig;
+
+		// If it exists, vertex on edge i is stored at position i.
+		// For edge numbering and location see numberings.png.
+		std::array<Vec3, 12> edgeVertIndices;
+	};
+
+	class Cube
+	{
+		Vec3 pos[8];
+		double sdf[8];
+
+		Vec3 LerpVertex(double isoLevel, int i1, int i2) const;
+		int SignConfig(double isoLevel) const;
+
+	public:
+		Cube(Rect3 const& space, Fun3s const& sdf);
+
+		// Find the vertices where the surface intersects the cube.
+		IntersectInfo Intersect(double isoLevel = 0) const;
+	};
+}

external/MeshReconstruction/lib/DataStructs.h

@@ -0,0 +1,57 @@
+#pragma once
+#include <vector>
+#include <array>
+#include <functional>
+#include <cmath>
+
+namespace MeshReconstruction
+{
+	struct Vec3
+	{
+		double x, y, z;
+
+		Vec3 operator+(Vec3 const& o) const
+		{
+			return { x + o.x, y + o.y, z + o.z };
+		}
+
+		Vec3 operator-(Vec3 const& o) const
+		{
+			return { x - o.x, y - o.y, z - o.z };
+		}
+
+		Vec3 operator*(double c) const
+		{
+			return { c*x, c*y, c*z };
+		}
+
+		double Norm() const
+		{
+			return sqrt(x*x + y*y + z*z);
+		}
+
+		Vec3 Normalized() const
+		{
+			auto n = Norm();
+			return { x / n, y / n, z / n };
+		}
+	};
+
+	struct Rect3
+	{
+		Vec3 min;
+		Vec3 size;
+	};
+
+	using Triangle = std::array<int, 3>;
+
+	struct Mesh
+	{
+		std::vector<Vec3> vertices;
+		std::vector<Triangle> triangles;
+		std::vector<Vec3> vertexNormals;
+	};
+
+	using Fun3s = std::function<double(Vec3 const&)>;
+	using Fun3v = std::function<Vec3(Vec3 const&)>;
+}

external/MeshReconstruction/lib/IO.cpp

@@ -0,0 +1,47 @@
+#include "IO.h"
+#include <stdexcept>
+using namespace std;
+
+void MeshReconstruction::WriteObjFile(Mesh const& mesh, string const& fileName)
+{
+	// FILE faster than streams.
+
+	FILE* file = fopen(fileName.c_str(), "w");
+	if (file == NULL)
+	{
+		throw runtime_error("Could not write obj file.");
+	}
+
+	// write stats
+	fprintf(file, "# %d vertices, %d triangles\n\n",
+		static_cast<int>(mesh.vertices.size()),
+		static_cast<int>(mesh.triangles.size()));
+
+	// vertices
+	for (size_t vi = 0; vi < mesh.vertices.size(); ++vi)
+	{
+		auto const& v = mesh.vertices.at(vi);
+		fprintf(file, "v %f %f %f\n", v.x, v.y, v.z);
+	}
+
+	// vertex normals
+	fprintf(file, "\n");
+	for (size_t ni = 0; ni < mesh.vertices.size(); ++ni)
+	{
+		auto const& vn = mesh.vertexNormals.at(ni);
+		fprintf(file, "vn %f %f %f\n", vn.x, vn.y, vn.z);
+	}
+
+	// triangles (1-based)
+	fprintf(file, "\n");
+	for (size_t ti = 0; ti < mesh.triangles.size(); ++ti)
+	{
+		auto const& t = mesh.triangles.at(ti);
+		fprintf(file, "f %d//%d %d//%d %d//%d\n",
+			t[0] + 1, t[0] + 1,
+			t[1] + 1, t[1] + 1,
+			t[2] + 1, t[2] + 1);
+	}
+
+	fclose(file);
+}

external/MeshReconstruction/lib/IO.h

@@ -0,0 +1,9 @@
+#pragma once
+#include <string>
+#include "DataStructs.h"
+
+namespace MeshReconstruction
+{
+	/// Writes a mesh to a file in <a href="https://de.wikipedia.org/wiki/Wavefront_OBJ">obj format</a>.
+	void WriteObjFile(Mesh const& mesh, std::string const& file);
+}

external/MeshReconstruction/lib/MeshReconstruction.cpp

@@ -0,0 +1,79 @@
+#include "MeshReconstruction.h"
+#include "Cube.h"
+#include "Triangulation.h"
+#include <iostream>
+
+using namespace MeshReconstruction;
+using namespace std;
+
+// Adapted from here: http://paulbourke.net/geometry/polygonise/
+
+namespace
+{
+	Vec3 NumGrad(Fun3s const& f, Vec3 const& p)
+	{
+		auto const Eps = 1e-6;
+		Vec3 epsX{ Eps, 0, 0 }, epsY{ 0, Eps, 0 }, epsZ{ 0, 0, Eps };
+		auto gx = (f(p + epsX) - f(p - epsX)) / 2;
+		auto gy = (f(p + epsY) - f(p - epsY)) / 2;
+		auto gz = (f(p + epsZ) - f(p - epsZ)) / 2;
+		return { gx, gy, gz };
+	}
+}
+
+Mesh MeshReconstruction::MarchCube(Fun3s const& sdf, Rect3 const& domain)
+{
+	auto const NumCubes = 50;
+	auto cubeSize = domain.size * (1.0 / NumCubes);
+
+	return MarchCube(sdf, domain, cubeSize);
+}
+
+Mesh MeshReconstruction::MarchCube(
+	Fun3s const& sdf,
+	Rect3 const& domain,
+	Vec3 const& cubeSize,
+	double isoLevel,
+	Fun3v sdfGrad)
+{
+	// Default value.
+	sdfGrad = sdfGrad == nullptr
+		? [&sdf](Vec3 const& p) { return NumGrad(sdf, p); }
+		: sdfGrad;
+
+	auto const NumX = static_cast<int>(ceil(domain.size.x / cubeSize.x));
+	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 CubeDiag = cubeSize.Norm();
+	// auto const HalfCubeDiag = cubeSize.Norm() * 0.5;
+	// auto const HalfCubeSize = cubeSize * 0.5;
+
+	Mesh mesh;
+
+	for (auto ix = 0; ix < NumX; ++ix)
+	{
+		auto x = domain.min.x + ix * cubeSize.x;
+		for (auto iy = 0; iy < NumY; ++iy)
+		{
+			auto y = domain.min.y + iy * cubeSize.y;
+			for (auto iz = 0; iz < NumZ; ++iz)
+			{
+				auto z = domain.min.z + iz * cubeSize.z;
+				Vec3 min{ x, y, z };
+
+				// Process only if cube lies within narrow band around surface.
+				// 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);
+				Triangulate(intersect, sdfGrad, mesh);
+			}
+		}
+	}
+
+	return mesh;
+}

external/MeshReconstruction/lib/MeshReconstruction.h

@@ -0,0 +1,27 @@
+#pragma once
+#include "DataStructs.h"
+
+namespace MeshReconstruction
+{
+	/// Reconstructs a triangle mesh from a given signed distance function using <a href="https://en.wikipedia.org/wiki/Marching_cubes">Marching Cubes</a>.
+	/// @param sdf The <a href="http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm">Signed Distance Function</a>.
+	/// @param domain Domain of reconstruction.
+	/// @returns The reconstructed mesh.
+	Mesh MarchCube(
+		Fun3s const& sdf,
+		Rect3 const& domain);
+
+	/// Reconstructs a triangle mesh from a given signed distance function using <a href="https://en.wikipedia.org/wiki/Marching_cubes">Marching Cubes</a>.
+	/// @param sdf The <a href="http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm">Signed Distance Function</a>.
+	/// @param domain Domain of reconstruction.
+	/// @param cubeSize Size of marching cubes. Smaller cubes yields meshes of higher resolution.
+	/// @param isoLevel Level set of the SDF for which triangulation should be done. Changing this value moves the reconstructed surface.
+	/// @param sdfGrad Gradient of the SDF which yields the vertex normals of the reconstructed mesh. If none is provided a numerical approximation is used.
+	/// @returns The reconstructed mesh.
+	Mesh MarchCube(
+		Fun3s const& sdf,
+		Rect3 const& domain,
+		Vec3 const& cubeSize,
+		double isoLevel = 0,
+		Fun3v sdfGrad = nullptr);
+}

external/MeshReconstruction/lib/Triangulation.cpp

@@ -0,0 +1,301 @@
+#include "Triangulation.h"
+
+namespace
+{
+	// Indices into vertex buffer (0 - 11).
+	// Three successive entries make up one triangle.
+	// -1 means unused.
+	const int signConfigToTriangles[256][16] =
+		{ { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1 },
+		{ 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1 },
+		{ 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1 },
+		{ 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1 },
+		{ 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1 },
+		{ 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1 },
+		{ 10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1 },
+		{ 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1 },
+		{ 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1 },
+		{ 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1 },
+		{ 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1 },
+		{ 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1 },
+		{ 11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1 },
+		{ 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1 },
+		{ 11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1 },
+		{ 11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1 },
+		{ 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1 },
+		{ 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1 },
+		{ 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1 },
+		{ 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1 },
+		{ 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1 },
+		{ 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1 },
+		{ 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1 },
+		{ 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1 },
+		{ 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1 },
+		{ 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1 },
+		{ 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1 },
+		{ 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1 },
+		{ 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1 },
+		{ 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1 },
+		{ 10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1 },
+		{ 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1 },
+		{ 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1 },
+		{ 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1 },
+		{ 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1 },
+		{ 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1 },
+		{ 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1 },
+		{ 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1 },
+		{ 10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1 },
+		{ 10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1 },
+		{ 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1 },
+		{ 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1 },
+		{ 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1 },
+		{ 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1 },
+		{ 11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1 },
+		{ 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1 },
+		{ 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1 },
+		{ 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1 },
+		{ 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1 },
+		{ 10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1 },
+		{ 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1 },
+		{ 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1 },
+		{ 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1 },
+		{ 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1 },
+		{ 10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1 },
+		{ 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1 },
+		{ 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1 },
+		{ 10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1 },
+		{ 10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1 },
+		{ 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1 },
+		{ 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1 },
+		{ 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1 },
+		{ 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1 },
+		{ 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1 },
+		{ 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1 },
+		{ 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1 },
+		{ 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1 },
+		{ 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1 },
+		{ 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1 },
+		{ 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1 },
+		{ 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1 },
+		{ 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1 },
+		{ 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1 },
+		{ 11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1 },
+		{ 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1 },
+		{ 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1 },
+		{ 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1 },
+		{ 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1 },
+		{ 10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1 },
+		{ 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1 },
+		{ 11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1 },
+		{ 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1 },
+		{ 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1 },
+		{ 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1 },
+		{ 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1 },
+		{ 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1 },
+		{ 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1 },
+		{ 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1 },
+		{ 10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1 },
+		{ 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1 },
+		{ 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1 },
+		{ 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1 },
+		{ 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1 },
+		{ 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1 },
+		{ 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1 },
+		{ 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1 },
+		{ 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1 },
+		{ 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1 },
+		{ 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1 },
+		{ 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1 },
+		{ 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1 },
+		{ 11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1 },
+		{ 11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1 },
+		{ 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1 },
+		{ 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1 },
+		{ 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1 },
+		{ 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1 },
+		{ 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1 },
+		{ 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1 },
+		{ 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1 },
+		{ 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
+		{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } };
+}
+
+/// Given a grid cube and an isolevel the triangles (5 max)
+/// required to represent the isosurface in the cube are computed.
+void MeshReconstruction::Triangulate(
+	IntersectInfo const& intersect,
+	Fun3v const& grad,
+	Mesh& mesh)
+{
+	// Cube is entirely in/out of the surface. Generate no triangles.
+	if (intersect.signConfig == 0 || intersect.signConfig == 255) return;
+
+	auto const& tri = signConfigToTriangles[intersect.signConfig];
+
+	for (auto i = 0; tri[i] != -1; i += 3)
+	{
+		auto const& v0 = intersect.edgeVertIndices[tri[i]];
+		auto const& v1 = intersect.edgeVertIndices[tri[i + 1]];
+		auto const& v2 = intersect.edgeVertIndices[tri[i + 2]];
+
+		mesh.vertices.push_back(v0);
+		mesh.vertices.push_back(v1);
+		mesh.vertices.push_back(v2);
+
+		auto normal0 = grad(v0).Normalized();
+		auto normal1 = grad(v1).Normalized();
+		auto normal2 = grad(v2).Normalized();
+
+		mesh.vertexNormals.push_back(normal0);
+		mesh.vertexNormals.push_back(normal1);
+		mesh.vertexNormals.push_back(normal2);
+
+		auto last = static_cast<int>(mesh.vertices.size() - 1);
+
+		mesh.triangles.push_back({ last - 2, last - 1, last });
+	}
+}

external/MeshReconstruction/lib/Triangulation.h

@@ -0,0 +1,11 @@
+#pragma once
+#include "Cube.h"
+#include "DataStructs.h"
+
+namespace MeshReconstruction
+{
+	void Triangulate(
+		IntersectInfo const& intersect,
+		Fun3v const& grad,
+		Mesh& mesh);
+}

notes_réparation_bac.txt

@@ -0,0 +1,5 @@
+Notes du TP de réparation de maillage par surface implicite de Mme Bac
+
+Fonction implicite f(x) = Σ(1,n)[αi φ(||X-Xi||) + …]
+
+φ(r) = r²ln(r) : fonction radiale, une par point, =0 au point

src/curvature.cpp

@@ -1,5 +1,6 @@
 #include "curvature.h"
 #include "util.h"
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
 
 
 void Courbures::normales_locales() {
@@ -9,7 +10,7 @@ void Courbures::normales_locales() {
 		MyMesh::Normal normal(0,0,0);
         i = 0;
 		for (MyMesh::FaceHandle vf : _mesh.vf_range(vh)) {
-            i++ ;
+            i++;
             normal += _mesh.calc_face_normal(vf);
         }
         if (i != 0) {
@@ -21,44 +22,41 @@ void Courbures::normales_locales() {
 }
 
 
-std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) {
-    OpenMesh::VPropHandleT<bool> vprop_flag;
-    _mesh.add_property(vprop_flag, "vprop_flag");
+void Courbures::get_two_neighborhood(std::vector<MyMesh::VertexHandle> &out,
+									 const MyMesh::VertexHandle vh) {
+	auto flag_prop =
+		OpenMesh::getOrMakeProperty<VertexHandle, bool>(_mesh, "vprop_flag");
+	flag_prop.set_range(_mesh.vertices_begin(), _mesh.vertices_end(), false);
 
-    // Initialisation
-    for (VertexHandle vh : _mesh.vertices()) {
-        _mesh.property(vprop_flag, vh) = false;
-    }
-    // Circulateur sur le premier cercle
-    std::vector<MyMesh::VertexHandle> neigh, neigh2;
-
-	_mesh.property(vprop_flag, vh) = true;
+	// Parcours du 1-anneau
+	flag_prop[vh] = true;
     for(VertexHandle vv : _mesh.vv_range(vh)) {
-        neigh.push_back(vv); // ajout du point à la liste
-        _mesh.property(vprop_flag, vv) = true;
+        out.push_back(vv); // ajout du point à la liste
+        flag_prop[vv] = true;
     }
 
-    // Parcours du premier cercle et ajout du second cercle par circulateurs
-    for (size_t i = 0; i < neigh.size(); i++) {
-        MyMesh::VertexHandle vh = neigh.at(i) ;
-        for (VertexHandle vv : _mesh.vv_range(vh)) {
-            if (!_mesh.property(vprop_flag, vv)) {
-				neigh2.push_back(vv);
+	if (out.size() >= 5) return;
+
+    // Parcours du 1-anneau de chaque sommet du 1-anneau
+	size_t old_size = out.size();
+    for (size_t i = 0; i < old_size; i++) {
+        for (VertexHandle vv : _mesh.vv_range(out[i])) {
+            if (!flag_prop[vv]) {
+				out.push_back(vv);
+				flag_prop[vv] = true;
 			}
-            _mesh.property(vprop_flag, vv) = true;
         }
     }
-
-    // Concaténation des deux cercles
-    neigh.insert(neigh.end(), neigh2.begin(), neigh2.end());
-    _mesh.remove_property(vprop_flag);
-    return neigh;
 }
 
 
 QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) {
-    std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh);
-    if (neigh.size() < 5) throw "Quad fitting: not enough neighbors";
+    static std::vector<MyMesh::VertexHandle> neigh;
+	neigh.clear();
+	get_two_neighborhood(neigh, vh);
+    if (neigh.size() < 5) {
+		throw std::runtime_error("Quad fitting: not enough neighbors");
+	}
 
     // Calcul de la matrice de changement de base
 	Eigen::Vector3d Oz(0,0,1);

src/curvature.h

@@ -14,14 +14,18 @@ public:
     OpenMesh::VPropHandleT<double> vprop_H;
 	OpenMesh::VPropHandleT<QuadPatch> vprop_quad;
 
-	Courbures(MyMesh &mesh) : _mesh(mesh) {}
+	Courbures(MyMesh &mesh) : _mesh(mesh) {
+		_mesh.request_vertex_normals();
+		_mesh.update_normals();
+	}
 
-    void set_fixed_colors() ;
-    void normales_locales() ;
-    std::vector<MyMesh::VertexHandle> get_two_neighborhood(MyMesh::VertexHandle vh);
+    void set_fixed_colors();
+    void normales_locales();
+    void get_two_neighborhood(std::vector<MyMesh::VertexHandle> &out,
+							  MyMesh::VertexHandle vh);
     QuadPatch fit_quad(MyMesh::VertexHandle vh);
-    void compute_KH() ;
-    void set_K_colors() ;
+    void compute_KH();
+    void set_K_colors();
 };
 
 

src/double_input.cpp

@@ -0,0 +1,59 @@
+#include "double_input.h"
+
+
+double DoubleInput::intToDouble(int value) const {
+	return static_cast<double>(value) / _slider_resolution
+		* (_max - _min) + _min;
+}
+
+
+int DoubleInput::doubleToInt(double value) const {
+	return (value - _min) / (_max - _min) * _slider_resolution;
+}
+
+
+void DoubleInput::onSpinBoxValueChanged(double value) {
+	emit valueChanged(value);
+	if (_propagate) {
+		_propagate = false;
+		_slider->setValue(doubleToInt(value));
+	} else {
+		_propagate = true;
+	}
+}
+
+
+void DoubleInput::onSliderValueChanged(int value) {
+	if (_propagate) {
+		_propagate = false;
+		_spin_box->setValue(intToDouble(value));
+	} else {
+		_propagate = true;
+	}
+}
+
+
+DoubleInput::DoubleInput(QObject *parent, double min, double max,
+						 double value,
+						 int slider_resolution)
+	:QObject(parent),
+	 _min(min),
+	 _max(max),
+	 _slider_resolution(slider_resolution),
+	 _spin_box(new QDoubleSpinBox()),
+	 _slider(new QSlider(Qt::Horizontal)) {
+	_spin_box->setRange(_min, _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),
+			this, &DoubleInput::onSpinBoxValueChanged);
+}
+
+
+void DoubleInput::setValue(double value) {
+	_spin_box->setValue(value);
+}

src/double_input.h

@@ -0,0 +1,41 @@
+#ifndef DOUBLE_INPUT_H
+#define DOUBLE_INPUT_H
+
+#include <QWidget>
+#include <QSlider>
+#include <QDoubleSpinBox>
+
+
+class DoubleInput : public QObject {
+	Q_OBJECT
+
+	const double _min;
+	const double _max;
+	const int _slider_resolution;
+	QDoubleSpinBox *_spin_box;
+	QSlider *_slider;
+	bool _propagate = true;
+
+	double intToDouble(int value) const;
+	int doubleToInt(double value) const;
+
+private slots:
+	void onSpinBoxValueChanged(double value);
+	void onSliderValueChanged(int value);
+
+public:
+	DoubleInput(QObject *parent, double min, double max, double value,
+				int slider_resolution=100);
+	QWidget *spinBox() { return _spin_box; }
+	QWidget *slider() { return _slider; }
+	double value() const { return _spin_box->value(); }
+
+signals:
+	void valueChanged(double value);
+
+public slots:
+	void setValue(double value);
+};
+
+
+#endif

src/hole_filling.cpp

@@ -0,0 +1,364 @@
+#include "hole_filling.h"
+#include "IO.h"
+#include "util.h"
+#include <MeshReconstruction.h>
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
+
+
+static std::vector<std::vector<HalfedgeHandle>> findHoles(MyMesh &mesh) {
+	std::vector<std::vector<HalfedgeHandle>> holes;
+	std::set<HalfedgeHandle> ignore;
+	for (HalfedgeHandle it : mesh.halfedges()) {
+		if (mesh.is_boundary(it)
+			&& ignore.find(it) == ignore.end()) {
+			holes.emplace_back();
+			holes.back().push_back(it);
+			ignore.insert(it);
+			for (HalfedgeHandle it2 : HalfedgeLoopRange<MyMesh>(mesh, it)) {
+				holes.back().push_back(it2);
+				ignore.insert(it2);
+			}
+		}
+	}
+	return holes;
+}
+
+
+void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole) {
+	mesh.request_vertex_status();
+	mesh.request_edge_status();
+	mesh.request_face_status();
+
+	Point center(0, 0, 0);
+	size_t length = 0;
+	for (HalfedgeHandle it : hole) {
+		VertexHandle vert = mesh.to_vertex_handle(it);
+		center += mesh.point(vert);
+		length++;
+	}
+	center /= length;
+	VertexHandle center_handle = mesh.new_vertex_dirty(center);
+	mesh.set_color(center_handle, mesh.default_color);
+
+	for (HalfedgeHandle it : hole) {
+		mesh.add_face(mesh.from_vertex_handle(it),
+					  mesh.to_vertex_handle(it),
+					  center_handle);
+	}
+}
+
+
+void fillHolesDumb(MyMesh &mesh) {
+	mesh.holes = findHoles(mesh);
+	for (auto hole : mesh.holes) {
+		fillHoleDumb(mesh, hole);
+	}
+}
+
+
+using namespace std ;
+using namespace MeshReconstruction;
+
+/* ************** Implicit RBF ************** */
+
+Implicit_RBF::Implicit_RBF(vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) : _alpha(alpha), _beta(beta), _center(center)
+{
+    _n = _alpha.size() ;
+    _d = _beta.size() ;
+    if (_center.size() != _n)
+        throw runtime_error("Inconsistent size of alpha and centers in Implicit_RBF constructor.");
+}
+
+// Computation of the value of the implicit surface at point X
+float Implicit_RBF::val(MyMesh::Point X) const
+{
+    float res = 0 ;
+    // Computation of the sum of RBF at centers
+    for(int i=0; i<_n; i++)
+    {
+        res += _alpha.at(i) * myphi((X-_center.at(i)).norm()) ;
+    }
+    // Computation of the polynomial part
+    for(int j=0; j<_d; j++)
+    {
+        res += _beta.at(j) * myp(j, X) ;
+    }
+    return res ;
+}
+
+/* ************** Hole Filling ************** */
+
+Hole_Filling::Hole_Filling(MyMesh &mesh, float scale, float discr)
+	: _mesh(mesh),
+	  _scale(scale),
+	  _discr(discr)
+{
+    _mesh.add_property(_vprop,  "vprop_flag");
+    cout << "Starting hole filling ..." << endl ;
+}
+
+// ***** Computation of boundary and its neighborhood
+
+MyMesh::HalfedgeHandle Hole_Filling::find_boundary_edge()
+{
+    MyMesh::HalfedgeIter he_it = _mesh.halfedges_begin();
+    while ( (he_it != _mesh.halfedges_end()) && (!_mesh.is_boundary(*he_it)))
+    {
+        ++he_it ;
+    }
+    if (he_it != _mesh.halfedges_end())
+        return *he_it ;
+    else
+        throw std::runtime_error("Boundary HE does not exist") ;
+}
+
+vector<MyMesh::VertexHandle> Hole_Filling::find_boundary(MyMesh::HalfedgeHandle heh)
+{
+    MyMesh::HalfedgeHandle heh_ini = heh ;
+    vector<MyMesh::VertexHandle> boundary ;
+
+    // Follow (and memorize) boundary edges
+    do
+    {
+        boundary.push_back(_mesh.to_vertex_handle(heh));
+        heh = _mesh.next_halfedge_handle(heh);
+    } while (heh != heh_ini) ;
+    return boundary ;
+}
+
+void Hole_Filling::init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd)
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin() ; v_it != _mesh.vertices_end(); ++v_it)
+        _mesh.property(_vprop, *v_it) = false ;
+    for (int i=0; i<bnd.size(); ++i)
+        _mesh.property(_vprop, bnd.at(i)) = true ;
+}
+
+vector<MyMesh::VertexHandle> Hole_Filling::next_neighbors(const vector<MyMesh::VertexHandle> & bnd)
+{
+    // Visit bnd vertices to find and mark next circle
+    vector<MyMesh::VertexHandle> next_bnd ;
+    for (int i=0; i<bnd.size(); i++)
+    {
+        for (MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(bnd.at(i));vv_it.is_valid();++vv_it)
+        {
+            if (_mesh.property(_vprop, *vv_it) == false) // new vertex
+            {
+                _mesh.property(_vprop, *vv_it) = true ;
+                next_bnd.push_back(*vv_it);
+            }
+        }
+    }
+    return next_bnd ;
+}
+
+// ***** Computation of RBF
+
+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)) ;
+    Eigen::MatrixXd Phi(3*n,3*n);
+    Eigen::MatrixXd P(3*n,d);
+    Eigen::VectorXd & B = *(new Eigen::VectorXd(3*n+d));
+
+    vector<MyMesh::Point> & pts_list = *(new vector<MyMesh::Point>);
+    //Append vertices to pts_list
+    for (int i=0; i<n; i++)
+    {
+        pts_list.push_back(_mesh.point(vlist.at(i))) ;
+    }
+    //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)) * 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)) * normal_scale) ;
+    }
+
+    int nn = pts_list.size() ;
+    // Compute corresponding B vector (0 / 1 / -1 )
+    B << Eigen::VectorXd::Zero(n), Eigen::VectorXd::Ones(n), -Eigen::VectorXd::Ones(n), Eigen::VectorXd::Zero(d) ;
+
+    // Fill Phi matrix
+	//TODO
+	for (int i = 0; i < 3*n; i++) {
+		for (int j = 0; j <= i; j++) {
+			Phi(i, j) = myphi((pts_list[i] - pts_list[j]).norm());
+			Phi(j, i) = Phi(i, j);
+		}
+	}
+
+    // Fill P matrix
+    for (int i = 0; i < 3*n; i++) {
+		for (int j = 0; j < d; j++) {
+			P(i, j) = myp(j, pts_list[i]);
+		}
+	}
+
+    // Set final A matrix
+    /* A = Phi | P
+     *      P' | 0    */
+	A << Phi, P, P.transpose(), Eigen::MatrixXd::Zero(d, d);
+    // TODO
+
+    cout << "size of pts_list : " << nn << endl ;
+    cout << "End computation of matrices" << endl ;
+
+    return {{A,B},pts_list} ;
+}
+
+pair<vector<float>&, vector<float>&> Hole_Filling::solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d)
+{
+    Eigen::MatrixXd & A = p.first ;
+    Eigen::VectorXd & B = p.second ;
+
+    Eigen::VectorXd res = A.householderQr().solve(B) ;
+
+    cout << "End of solver" << endl ;
+    cout << "res : " << res.head(10) << endl ;
+    vector<float> & alpha = *(new vector<float>);
+    vector<float> & beta = *(new vector<float>);
+
+    if (res.size() != (n+d))
+    {
+        cout << "taille du res : " << res.size() << endl ;
+        throw std::runtime_error("Error in solve_approx") ;
+    }
+    for (int i=0; i<n; i++)
+    {
+        alpha.push_back(res(i)) ;
+    }
+
+    for (int j=0; j<d; j++)
+    {
+        beta.push_back(res(n+j)) ;
+    }
+
+    return {alpha, beta} ;
+}
+
+// ***** IO
+
+void Hole_Filling::colorize_prop()
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
+    {
+        if(_mesh.property(_vprop, *v_it) == true)
+            _mesh.set_color(*v_it, MyMesh::Color(255, 0, 0)) ;
+        else
+			_mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
+    }
+}
+
+void Hole_Filling::colorize_verts(const vector<MyMesh::VertexHandle> &vlist)
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
+    {
+           _mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
+    }
+
+    for(int i=0; i<vlist.size(); i++)
+    {
+        _mesh.set_color(vlist.at(i), MyMesh::Color(255, 0, 0)) ;
+    }
+}
+
+Rect3 Hole_Filling::estimate_BB(const vector<MyMesh::VertexHandle> &vlist)
+{
+    MyMesh::Point minp = _mesh.point(vlist.at(0)), maxp=minp ;
+    for (int i=0; i<vlist.size(); i++)
+    {
+        minp = min(minp, _mesh.point(vlist.at(i))) ;
+        maxp = max(maxp, _mesh.point(vlist.at(i))) ;
+    }
+    MyMesh::Point sizep(maxp-minp), centerp = (minp+maxp)/2 ;
+    minp = centerp - _scale*sizep/2 ;
+    sizep *= _scale ;
+
+    Rect3 domain ;
+    domain.min = {minp[0], minp[1], minp[2]} ;
+    domain.size = {sizep[0], sizep[1], sizep[2]} ;
+    return domain ;
+}
+
+Mesh Hole_Filling::poly_n_out(const Implicit_RBF &implicit, Rect3 domain)
+{
+    auto implicitEq = [&implicit](Vec3 const& pos)
+    {
+        MyMesh::Point X(pos.x, pos.y, pos.z) ;
+        return implicit.val(X) ;
+    };
+
+
+    Vec3 cubeSize(domain.size*(_discr)) ;
+    cout << "mind " << domain.min << endl ;
+    cout << "sized" << domain.size << endl ;
+    cout << "cubesize "<< cubeSize << endl ;
+
+    return MarchCube(implicitEq, domain, cubeSize);
+}
+
+
+/* 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, diag * .1);
+	auto [alpha, beta] = hf.solve_approx(system, pts_list.size(), 10);
+	Implicit_RBF rbf(alpha, beta, pts_list);
+
+	Mesh filling = 	hf.poly_n_out(rbf, bb);
+	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;
+}
+
+
+std::vector<MyMesh> fillHolesImplicit(MyMesh &mesh,
+									  float scale, float 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;
+}

src/hole_filling.h

@@ -0,0 +1,124 @@
+#ifndef HOLE_FILLING_H
+#define HOLE_FILLING_H
+
+#include "my_mesh.h"
+#include <vector>
+#include <iostream>
+#include <cmath>
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <Eigen/Dense>
+#include <MeshReconstruction.h>
+#include <IO.h>
+
+
+void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole);
+void fillHolesDumb(MyMesh &mesh);
+std::vector<MyMesh> fillHolesImplicit(MyMesh &mesh, float scale, float discr);
+
+using namespace std ;
+using namespace MeshReconstruction;
+
+// ******************************
+// Function computing polynomial P_i (degree 2 polynomials in 3 variables)
+inline float myp(int i, MyMesh::Point X) {
+    float x = X[0] ;
+    float y = X[1] ;
+    float z = X[2] ;
+    switch (i) {
+    case 0:
+        return x*x ;
+        break;
+    case 1:
+        return y*y ;
+        break;
+    case 2:
+        return z*z ;
+        break;
+    case 3:
+        return x*y ;
+        break;
+    case 4:
+        return y*z ;
+        break;
+    case 5:
+        return x*z ;
+        break;
+    case 6:
+        return x ;
+        break;
+    case 7:
+        return y ;
+        break;
+    case 8:
+        return z;
+        break;
+    case 9:
+        return 1;
+        break;
+    default:
+        throw std::runtime_error("Error on indice i : unknown polynomial P.");
+    }
+}
+
+// ******************************
+// Function computing thin spline RBF
+inline float myphi(float r) { if (r == 0) return 0 ; else return r*r*log(r) ; }
+
+// ******************************
+// Class encoding an implicit function built from the basis of RBF
+// f(X) = Sum_i=0^n-1 alpha_i * phi(|| X-c_i ||) + Sum_j=0^9 beta_j * P_j(X)
+class Implicit_RBF {
+private:
+    int _n ;
+    int _d ;
+    vector<float> _alpha, _beta ;
+    vector<MyMesh::Point> _center ;
+
+public:
+    Implicit_RBF (vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) ;
+    float val(MyMesh::Point) const ;
+};
+
+// ******************************
+// Class encoding all the functions required for implicit hole filling
+class Hole_Filling
+{
+private:
+    MyMesh &_mesh ;
+    OpenMesh::VPropHandleT<bool> _vprop ;
+
+    // Constantes servant au rendu
+    // TODO : à mettre dans des boutons ...
+    const float _scale; // Boîte englobante de rendu = _scale * BB des centres
+    const float _discr; // BB discrétisée en 1/_discr voxels
+public:
+    Hole_Filling(MyMesh &mesh, float scale, float discr);
+    inline ~Hole_Filling() { std::cout << "Ending hole filling" << std::endl ; }
+
+    // Computation of boundary and its neighborhood
+    MyMesh::HalfedgeHandle find_boundary_edge() ;
+    vector<MyMesh::VertexHandle> find_boundary(MyMesh::HalfedgeHandle heh) ;
+    void init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd) ;
+    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, double normal_scale=1) ;
+    pair<vector<float>&, vector<float>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;
+
+    // IO
+    void colorize_prop() ;
+    void colorize_verts(const vector<MyMesh::VertexHandle> &vlist) ;
+    Rect3 estimate_BB(const vector<MyMesh::VertexHandle> &vlist) ;
+    Mesh poly_n_out (const Implicit_RBF & implicit, Rect3 domain) ;
+};
+
+// Other ....
+
+inline MyMesh::Point min (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = min(X[i], Y[i]) ; return P ;}
+inline MyMesh::Point max (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = max(X[i], Y[i]) ; return P ;}
+
+inline ostream & operator<< (ostream & out, Vec3 v) { out << v.x << ", " << v.y << ", " << v.z ; return out ;}
+
+
+#endif

src/main.cpp

@@ -5,11 +5,37 @@
 #include <QApplication>
 
 
+static MeshProcessor *create_mesh_processor(const QString &path,
+											MainWindow &main_window) {
+	MeshViewer &mesh_viewer = main_window.mesh_viewer;
+	MeshProcessor *mesh_processor = new MeshProcessor(path, mesh_viewer,
+													  main_window.fillHolesImplicitScale(),
+													  main_window.fillHolesImplicitDiscr());
+	QObject::connect(&main_window, &MainWindow::fillHolesDumbClicked,
+					 mesh_processor, &MeshProcessor::fillHolesDumb);
+	QObject::connect(&main_window, &MainWindow::fillHolesImplicitClicked,
+					 mesh_processor, &MeshProcessor::fillHolesImplicit);
+	QObject::connect(&main_window, &MainWindow::smoothUniformClicked,
+					 mesh_processor, &MeshProcessor::smoothUniform);
+	QObject::connect(&main_window, &MainWindow::smoothCotangentClicked,
+					 mesh_processor, &MeshProcessor::smoothCotangent);
+	QObject::connect(&main_window, &MainWindow::patchViewToggled,
+					 mesh_processor, &MeshProcessor::setPatchView);
+	QObject::connect(&main_window, &MainWindow::fillHolesImplicitScaleChanged,
+					 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;
+}
+
+
 int main(int argc, char *argv[]) {
 	using namespace std;
 	QSurfaceFormat format;
 	format.setRenderableType(QSurfaceFormat::OpenGL);
-#ifndef QT_DEBUG
+#ifdef QT_DEBUG
 	qDebug("Debug build");
 	format.setOption(QSurfaceFormat::DebugContext);
 #endif
@@ -17,29 +43,17 @@ int main(int argc, char *argv[]) {
 	QApplication app(argc, argv);
 	MeshProcessor *mesh_processor = nullptr;
 	MainWindow main_window;
-	MeshViewer *mesh_viewer = &main_window.mesh_viewer;
-	QObject::connect(mesh_viewer, &MeshViewer::initialized,
-					 [&]() {
-						 if (mesh_processor) {
-							 mesh_viewer->addMesh(mesh_processor->mesh);
-							 mesh_viewer->addMesh(mesh_processor->patch);
-						 }
-					 });
-
 	if (argc > 2) {
 		qWarning("Utilisation : %s [MAILLAGE]", argv[0]);
 		return 1;
 	} else if (argc == 2) {
-		mesh_processor = new MeshProcessor(argv[1]);
+		mesh_processor = create_mesh_processor(argv[1], main_window);
 	}
 	QObject::connect(&main_window, &MainWindow::open,
 					 [&](const QString &path) {
-						 if (mesh_processor) {
-							 mesh_viewer->removeMesh(mesh_processor->mesh);
-							 delete mesh_processor;
-						 }
-						 mesh_processor = new MeshProcessor(path);
-						 mesh_viewer->addMesh(mesh_processor->mesh);
+						 if (mesh_processor) delete mesh_processor;
+						 mesh_processor = create_mesh_processor
+							 (path, main_window);
 					 });
 	main_window.show();
 	return app.exec();

src/main_window.cpp

@@ -1,35 +1,132 @@
 #include "main_window.h"
 #include "mesh_processor.h"
+#include "double_input.h"
 
 #include <QApplication>
 #include <QFileDialog>
+#include <QMenuBar>
+#include <QGroupBox>
+#include <QPushButton>
+#include <QSlider>
+#include <QLabel>
+#include <qnamespace.h>
 
 
 MainWindow::MainWindow(QWidget *parent)
 	:QMainWindow(parent),
 	 toolbar(this),
-	 mesh_viewer(this) {
-	connect(&mesh_viewer, &MeshViewer::initialized, [&]() {
-		open_action->setEnabled(true);
-	});
+	 mesh_viewer() {
 	setCentralWidget(&mesh_viewer);
-	addToolBar(Qt::RightToolBarArea, &toolbar);
-	open_action = toolbar.addAction("Ouvrir…", [&](){
+
+	QMenuBar *menu_bar = new QMenuBar();
+	setMenuBar(menu_bar);
+
+	// File menu
+	QMenu *file_menu = new QMenu("Fichier");
+	open_action = file_menu->addAction("Ouvrir…", [&](){
 		emit open(QFileDialog::getOpenFileName(this, "Ouvrir un maillage"));
 	});
-	// toolbar_actions.append(toolbar.addAction("Fractionner", [&](){
-	// 	QVector<QPair<MyMesh::Point, MyMesh>> fragments = shatter(mesh);
-	// 	mesh_viewer.removeOpenGLMesh(glm);
-	// 	for (auto &[pos, fragment] : fragments) {
-	// 		fragment.triangulate();
-	// 		QMatrix4x4 mat;
-	// 		float scale = 1.2;
-	// 		mat.translate(pos[0] * scale, pos[1] * scale, pos[2] * scale);
-	// 		mesh_viewer.addOpenGLMeshFromOpenMesh(&fragment, mat);
-	// 	}
-	// }));
-	open_action->setEnabled(false);
-	for (QAction *a : toolbar_actions) {
-		a->setEnabled(false);
+	save_action = file_menu->addAction("Enregistrer sous…", [&]() {
+		emit save(QFileDialog::getSaveFileName(this,
+											   "Enregistrer un maillage"));
+	});
+	menu_bar->addMenu(file_menu);
+	if (!mesh_viewer.isInitialized()) {
+		open_action->setEnabled(false);
+		connect(&mesh_viewer, &MeshViewer::initialized, [&]() {
+			open_action->setEnabled(true);
+		});
 	}
+
+	addToolBar(Qt::RightToolBarArea, &toolbar);
+
+	// Hole filling tools
+	QGroupBox *hole_box = new QGroupBox("Remplissage de trous");
+	QGridLayout *hole_layout = new QGridLayout();
+	hole_box->setLayout(hole_layout);
+
+	QPushButton *fill_holes_dumb = new QPushButton("Remplir bêtement");
+	connect(fill_holes_dumb, &QPushButton::clicked,
+			this, &MainWindow::fillHolesDumbClicked);
+	hole_layout->addWidget(fill_holes_dumb, 0, 0);
+
+	QPushButton *fill_holes_implicit =
+		new QPushButton("Remplir par une surface implicite");
+	connect(fill_holes_implicit, &QPushButton::clicked,
+			this, &MainWindow::fillHolesImplicitClicked);
+	hole_layout->addWidget(fill_holes_implicit, 1, 0);
+
+	QLabel *implicit_scale_text =
+		new QLabel("Échelle du remplissage implicite", this);
+	hole_layout->addWidget(implicit_scale_text, 2, 0);
+	fill_holes_implicit_scale = new DoubleInput(this, 0, 10, 4);
+	connect(fill_holes_implicit_scale, &DoubleInput::valueChanged,
+			this, &MainWindow::fillHolesImplicitScaleChanged);
+	hole_layout->addWidget(fill_holes_implicit_scale->slider(), 3, 0);
+	hole_layout->addWidget(fill_holes_implicit_scale->spinBox(), 3, 1);
+
+	QLabel *implicit_discr_text =
+		new QLabel("Taux de discrétisation du remplissage implicite", this);
+	hole_layout->addWidget(implicit_discr_text, 4, 0);
+	fill_holes_implicit_discr = new DoubleInput(this, .02, .1, 1/40.);
+	connect(fill_holes_implicit_discr, &DoubleInput::valueChanged,
+			this, &MainWindow::fillHolesImplicitDiscrChanged);
+	hole_layout->addWidget(fill_holes_implicit_discr->slider(), 5, 0);
+	hole_layout->addWidget(fill_holes_implicit_discr->spinBox(), 5, 1);
+
+	toolbar.addWidget(hole_box);
+
+
+	// Smoothing tools
+	QGroupBox *smooth_box = new QGroupBox("Adoucissement");
+	QGridLayout *smooth_layout = new QGridLayout();
+	smooth_box->setLayout(smooth_layout);
+	QPushButton *smooth = new QPushButton("Adoucir (uniforme)");
+	connect(smooth, &QPushButton::clicked,
+			this, &MainWindow::smoothUniformClicked);
+	smooth_layout->addWidget(smooth, 1, 0);
+	QPushButton *smooth_cotan = new QPushButton("Adoucir (cotangent)");
+	smooth_cotangent_factor_input =
+		new DoubleInput(this, .00001, .001, .0001);
+	connect(smooth_cotangent_factor_input, &DoubleInput::valueChanged,
+			[&](double value) { smooth_cotangent_factor = value; });
+	connect(smooth_cotan, &QPushButton::clicked,
+			[&]() { emit smoothCotangentClicked(smooth_cotangent_factor); });
+	smooth_layout->addWidget(smooth_cotan, 2, 0);
+	QLabel *smooth_cotan_text =
+		new QLabel("Facteur de l'adoucissement cotangentiel", this);
+	smooth_layout->addWidget(smooth_cotan_text, 3, 0);
+	smooth_layout->addWidget(smooth_cotangent_factor_input->slider(), 4, 0);
+	QDoubleSpinBox *sb = (QDoubleSpinBox *)(smooth_cotangent_factor_input->spinBox());
+	sb->setDecimals(5);
+	smooth_layout->addWidget(smooth_cotangent_factor_input->spinBox(), 4, 1);
+	toolbar.addWidget(smooth_box);
+
+
+	// Curvature tools
+	QGroupBox *curvature_box = new QGroupBox("Analyse de courbure");
+	QLayout *curvature_layout = new QVBoxLayout();
+	curvature_box->setLayout(curvature_layout);
+	QPushButton *patch_mode = new QPushButton(
+		"Afficher le patch de la sélection");
+	patch_mode->setCheckable(true);
+	connect(patch_mode, &QPushButton::toggled,
+			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

@@ -1,11 +1,14 @@
 #ifndef MAIN_WINDOW_H
 #define MAIN_WINDOW_H
 
-#include <QMainWindow>
-#include <QToolBar>
-
 #include "mesh_viewer.h"
 #include "my_mesh.h"
+#include "double_input.h"
+
+#include <QMainWindow>
+#include <QToolBar>
+#include <QVBoxLayout>
+#include <QDoubleSpinBox>
 
 
 class MainWindow : public QMainWindow {
@@ -13,15 +16,34 @@ class MainWindow : public QMainWindow {
 
 	QToolBar toolbar;
 	QAction *open_action;
-	QList<QAction *> toolbar_actions;
+	QAction *save_action;
+	DoubleInput *fill_holes_implicit_scale;
+	DoubleInput *fill_holes_implicit_discr;
+	DoubleInput *smooth_cotangent_factor_input;
+	double smooth_cotangent_factor;
 
 signals:
 	void open(const QString &path);
+	void save(const QString &path);
+	void fillHolesDumbClicked();
+	void fillHolesImplicitClicked();
+	void fillHolesImplicitScaleChanged(float value);
+	void fillHolesImplicitDiscrChanged(float value);
+	void smoothUniformClicked();
+	void smoothCotangentClicked(double factor);
+	void patchViewToggled(bool checked);
+	void filterNoiseClicked(int value);
 
 public:
 	MeshViewer mesh_viewer;
 
 	MainWindow(QWidget *parent=nullptr);
+	double fillHolesImplicitScale() const {
+		return fill_holes_implicit_scale->value();
+	}
+	double fillHolesImplicitDiscr() const {
+		return fill_holes_implicit_discr->value();
+	}
 };
 
 

src/mesh_processor.cpp

@@ -1,75 +1,144 @@
 #include "quad_patch_tesselator.h"
 #include "mesh_processor.h"
 #include "util.h"
+#include "hole_filling.h"
 #include "smoothing.h"
 #include "curvature.h"
+#include "noise_removal.h"
 #include <QGenericMatrix>
 #include <unordered_set>
 
 
-static std::vector<std::vector<HalfedgeHandle>> findHoles(MyMesh &mesh) {
-	std::vector<std::vector<HalfedgeHandle>> holes;
-	std::set<HalfedgeHandle> ignore;
-	for (HalfedgeHandle it : mesh.halfedges()) {
-		if (mesh.is_boundary(it)
-			&& ignore.find(it) == ignore.end()) {
-			holes.emplace_back();
-			holes.back().push_back(it);
-			ignore.insert(it);
-			for (HalfedgeHandle it2 : HalfedgeLoopRange<MyMesh>(mesh, it)) {
-				holes.back().push_back(it2);
-				ignore.insert(it2);
-			}
-		}
-	}
-	return holes;
-}
-
-
-MeshProcessor::MeshProcessor(const QString &path) {
+MeshProcessor::MeshProcessor(const QString &path, MeshViewer &mesh_viewer,
+							 double implicit_hole_filling_scale,
+							 double implicit_hole_filling_discr)
+	:mesh_viewer(mesh_viewer),
+	 implicit_hole_filling_scale(implicit_hole_filling_scale),
+	 implicit_hole_filling_discr(implicit_hole_filling_discr) {
 	OpenMesh::IO::Options options;
 	options.set(OpenMesh::IO::Options::VertexColor);
 	if (!OpenMesh::IO::read_mesh(mesh, path.toUtf8().constData(), options)) {
 		qWarning() << "Failed to read" << path;
 		return;
 	}
-	mesh.update_normals();
-	Courbures courb(mesh);
-	courb.compute_KH();
-    courb.set_K_colors();
-	VertexHandle vh = mesh.vertex_handle(16);
-	QuadPatch q = mesh.property(courb.vprop_quad, vh);
-	patch = tesselate_quad_patch(q, mesh, vh);
-	// mesh.holes = findHoles(mesh);
-	// fillHoles();
-	// smooth(mesh);
+	for (VertexHandle vh : mesh.vertices()) {
+		mesh.set_color(vh, MyMesh::Color(.5, .5, .5));
+	}
+	courbure = new Courbures(mesh);
+	try {
+		courbure->compute_KH();
+		courbure->set_K_colors();
+	} catch (std::runtime_error &e) {
+		qWarning() << "Curvature computation failed";
+	}
+	connect(&mesh_viewer, &MeshViewer::clicked, this, &MeshProcessor::click);
+
+	updateView();
 }
 
 
-void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole) {
-	mesh.request_vertex_status();
-	mesh.request_edge_status();
-	mesh.request_face_status();
-
-	Point center(0, 0, 0);
-	size_t length = 0;
-	for (HalfedgeHandle it : hole) {
-		VertexHandle vert = mesh.to_vertex_handle(it);
-		center += mesh.point(vert);
-		length++;
+MeshProcessor::~MeshProcessor() {
+	if (mesh_viewer.isInitialized()) {
+		mesh_viewer.removeMesh(mesh);
 	}
-	center /= length;
-	VertexHandle center_handle = mesh.new_vertex_dirty(center);
+	if (courbure) delete courbure;
+}
 
-	for (HalfedgeHandle it : hole) {
-		mesh.add_face(mesh.from_vertex_handle(it),
-					  mesh.to_vertex_handle(it),
-					  center_handle);
+
+void MeshProcessor::updateView() {
+	if (mesh_viewer.isInitialized()) {
+		mesh_viewer.removeMesh(mesh);
+		mesh_viewer.addMesh(mesh);
+		mesh_viewer.updateForReal();
+	} else {
+		connect(&mesh_viewer, &MeshViewer::initialized, [&]() {
+			mesh_viewer.addMesh(mesh);
+			mesh_viewer.updateForReal();
+		});
 	}
 }
 
-void MeshProcessor::fillHoles() {
-	for (auto hole : mesh.holes) {
-		fillHoleDumb(mesh, hole);
+
+void MeshProcessor::fillHolesDumb() {
+	::fillHolesDumb(mesh);
+	updateView();
+}
+
+
+void MeshProcessor::fillHolesImplicit() {
+	for (MyMesh &filling : fillings) {
+		mesh_viewer.removeMesh(filling);
+	}
+	std::vector<MyMesh> new_fillings =
+		::fillHolesImplicit(mesh,
+							implicit_hole_filling_scale,
+							implicit_hole_filling_discr);
+	for (MyMesh &filling : new_fillings) {
+		fillings.push_back(filling);
+		mesh_viewer.addMesh(fillings[fillings.size() - 1]);
+	}
+	updateView();
+}
+
+
+void MeshProcessor::setImplicitHoleFillingScale(float scale) {
+	implicit_hole_filling_scale = scale;
+	if (fillings.size() > 0) {
+		fillHolesImplicit();
 	}
 }
+
+
+void MeshProcessor::setImplicitHoleFillingDiscr(float discr) {
+	implicit_hole_filling_discr = discr;
+	if (fillings.size() > 0) {
+		fillHolesImplicit();
+	}
+}
+
+
+void MeshProcessor::smoothCotangent(double factor) {
+	::smooth(mesh, SmoothingMethod::COTANGENT, factor);
+	updateView();
+}
+
+
+void MeshProcessor::smoothUniform() {
+	::smooth(mesh, SmoothingMethod::UNIFORM);
+	updateView();
+}
+
+
+void MeshProcessor::setPatchView(bool on) {
+	view_patches = on;
+}
+
+
+void MeshProcessor::click(QVector3D position) {
+	if (!view_patches) return;
+	Eigen::Vector3d pos {position.x(), position.y(), position.z()};
+	MyMesh::VertexIter it = mesh.vertices_begin();
+	VertexHandle min_vert = *it;
+	double min_dist = (mesh.point(*it) - pos).squaredNorm();
+	for (; it != mesh.vertices_end(); ++it) {
+		double dist = (mesh.point(*it) - pos).squaredNorm();
+		if (dist < min_dist) {
+			min_dist = dist;
+			min_vert = *it;
+		}
+	}
+	QuadPatch q = mesh.property(courbure->vprop_quad, min_vert);
+	patch = tesselate_quad_patch(q, mesh, min_vert);
+	mesh_viewer.removeMesh(patch);
+	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

@@ -2,7 +2,8 @@
 #define MESH_PROCESSING_H
 
 #include "my_mesh.h"
-#include "mesh_processor_painter.h"
+#include "curvature.h"
+#include "mesh_viewer.h"
 #include <vector>
 #include <QObject>
 
@@ -10,14 +11,34 @@
 class MeshProcessor : public QObject {
 	Q_OBJECT
 
+	Courbures *courbure = nullptr;
+	MeshViewer &mesh_viewer;
+	bool view_patches = false;
+	double implicit_hole_filling_scale;
+	double implicit_hole_filling_discr;
+	std::vector<MyMesh> fillings;
+
+	void updateView();
+
 public:
 	MyMesh mesh;
 	MyMesh patch;
 
-	MeshProcessor(const QString &path);
+	MeshProcessor(const QString &path, MeshViewer &mesh_viewer,
+				  double implicit_hole_filling_scale,
+				  double implicit_hole_filling_discr);
+	~MeshProcessor();
 
 public slots:
-	void fillHoles();
+	void fillHolesDumb();
+	void fillHolesImplicit();
+	void setImplicitHoleFillingScale(float scale);
+	void setImplicitHoleFillingDiscr(float discr);
+	void smoothCotangent(double factor);
+	void smoothUniform();
+	void setPatchView(bool on);
+	void click(QVector3D position);
+	void removeNoise(int value);
 };
 
 

src/mesh_processor_painter.cpp

@@ -1,6 +0,0 @@
-#include "mesh_processor_painter.h"
-
-
-void MeshProcessorPainter::paint() {
-	
-}

src/mesh_processor_painter.h

@@ -1,15 +0,0 @@
-#ifndef MESH_PAINTER_H
-#define MESH_PAINTER_H
-
-#include <QObject>
-
-
-class MeshProcessorPainter {
-	Q_OBJECT
-
-public slots:
-	void paint();
-};
-
-
-#endif

src/mesh_viewer.cpp

@@ -9,6 +9,12 @@
 MeshViewer::MeshViewer(QWidget *parent) : QOpenGLWidget(parent) {
 	setMouseTracking(true);
 	setFocus();
+	updateViewMatrix();
+}
+
+
+void MeshViewer::updateViewMatrix() {
+	view = zoom * rot * trans;
 }
 
 
@@ -59,22 +65,22 @@ void MeshViewer::initializeGL() {
 	glf->glEnable(GL_MULTISAMPLE);
 
 	qDebug("MeshViewer: initialization complete");
+	is_initialized = true;
 	emit initialized();
 }
 
 
 void MeshViewer::resizeGL(int w, int h) {
-	QMatrix4x4 projection;
-	projection.perspective(FOV, (float) w/h, .01, 100);
-	program.setUniformValue("proj", projection);
+	proj.setToIdentity();
+	proj.perspective(FOV, (float) w/h, .01, 100);
+	program.setUniformValue("proj", proj);
+	update();
 }
 
 
 void MeshViewer::paintGL() {
 	// glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 	QMatrix4x4 trans;
-	trans.translate(0, 0, -cam_dist);
-	QMatrix4x4 view = trans * rot;
 	program.bind();
 	program.setUniformValue("view", view);
 	for (MeshView &m : meshes) {
@@ -83,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();
@@ -94,39 +101,81 @@ 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();
 }
 
 
+void MeshViewer::updateForReal() {
+	setEnabled(true);
+	setVisible(true);
+	update();
+	repaint();
+	makeCurrent();
+	paintGL();
+	doneCurrent();
+}
+
+
 void MeshViewer::mousePressEvent(QMouseEvent *e) {
-	if (e->button() == Qt::LeftButton) {
+	if (e->button() & Qt::MiddleButton || e->button() & Qt::RightButton) {
 		mouse_pos = e->pos();
+	} else if (e->button() == Qt::LeftButton) {
+		float x = e->x();
+		float y = height() - e->y();
+		float depth;
+		makeCurrent();
+		QOpenGLFunctions *glf = context()->functions();
+		glf->glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &depth);
+		doneCurrent();
+		QVector3D position {x, y, depth};
+		position = position.unproject(view, proj, rect());
+		emit clicked(QVector3D(position));
+	} else {
+		e->ignore();
 	}
 }
 
 
 void MeshViewer::mouseReleaseEvent(QMouseEvent *e) {
 	(void) e;
-	rot_start = rot;
 }
 
 
 void MeshViewer::mouseMoveEvent(QMouseEvent *e) {
-	if (e->buttons() & Qt::LeftButton) {
-		QPoint delta = e->pos() - mouse_pos;
-		rot = rot_start;
+	QPoint delta = e->pos() - mouse_pos;
+	mouse_pos = e->pos();
+	if (e->buttons() & Qt::MiddleButton) {
 		rot.rotate(delta.x() / 5., 0, 1, 0);
 		rot.rotate(delta.y() / 5., QVector3D(1, 0, 0) * rot);
+		updateViewMatrix();
+		update();
+	}
+	if (e->buttons() & Qt::RightButton) {
+		QMatrix4x4 screen_trans;
+		screen_trans.translate(delta.x() / 1000., -delta.y() / 1000., 0);
+		trans = rot.inverted() * screen_trans * rot * trans;
+		// trans(0, 0) = 1;
+		// trans(0, 1) = 0;
+		// trans(1, 0) = 0;
+		// trans(1, 1) = 1;
+		updateViewMatrix();
 		update();
 	}
 }
 
 
 void MeshViewer::wheelEvent(QWheelEvent *e) {
-	cam_dist -= e->angleDelta().y() / 1000. * cam_dist;
+	zoom(2, 3) -= e->angleDelta().y() / 1000. * zoom(2, 3);
+	zoom.optimize();
+	updateViewMatrix();
 	update();
 }

src/mesh_viewer.h

@@ -15,7 +15,7 @@
 
 
 #define WIREFRAME_COLOR 0, 0, 0
-#define FOV 70
+#define FOV 40
 
 
 using namespace OpenMesh;
@@ -27,9 +27,15 @@ class MeshViewer : public QOpenGLWidget {
 	std::list<MeshView> meshes;
 	QOpenGLShaderProgram program;
 	QMatrix4x4 proj;
-	QMatrix4x4 rot, rot_start;
-	GLfloat cam_dist = 1;
+	QMatrix4x4 zoom = QMatrix4x4(1, 0, 0, 0,
+								 0, 1, 0, 0,
+								 0, 0, 1, -1,
+								 0, 0, 0, 1);
+	QMatrix4x4 rot, trans;
+	QMatrix4x4 view;
 	QPoint mouse_pos;
+	void updateViewMatrix();
+	bool is_initialized = false;
 
 public:
     MeshViewer(QWidget *parent=nullptr);
@@ -37,10 +43,12 @@ public:
     void initializeGL() override;
     void resizeGL(int w, int h) override;
     void paintGL() override;
+	constexpr bool isInitialized() { return is_initialized; }
 
 public slots:
-	void addMesh(const MyMesh &mesh);
+	void addMesh(MyMesh &mesh);
 	void removeMesh(const MyMesh &mesh);
+	void updateForReal();
 
 protected:
     virtual void mousePressEvent(QMouseEvent *e);
@@ -50,6 +58,7 @@ protected:
 
 signals:
 	void initialized();
+	void clicked(QVector3D position);
 };
 
 

src/my_mesh.h

@@ -14,17 +14,19 @@ 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);
-    EdgeAttributes(OpenMesh::Attributes::Color);
+    FaceAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Status);
+    // EdgeAttributes(OpenMesh::Attributes::Color);
 };
 
 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;
-	// QColor color = {127, 127, 127};
 };
 
 typedef MyMesh::FaceHandle FaceHandle;

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/quad_patch_tesselator.cpp

@@ -2,10 +2,13 @@
 
 
 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)) {
@@ -17,10 +20,12 @@ MyMesh tesselate_quad_patch(QuadPatch q, MyMesh mesh, VertexHandle vh) {
 		ymin = std::min(ymin, point[1]);
 		ymax = std::max(ymax, point[1]);
 	}
-	double xstep = (xmax-xmin)/64.;
-	double ystep = (ymax-ymin)/64.;
-	for (size_t y = 0; y < 64; y++) {
-		for (size_t x = 0; x < 64; x++) {
+
+	// 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));
@@ -28,17 +33,22 @@ MyMesh tesselate_quad_patch(QuadPatch q, MyMesh mesh, VertexHandle vh) {
 			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 % 64;
-		size_t y = i / 64;
-		if (x == 63 || y == 63) continue;
-		patch.add_face(vhi, patch.vertex_handle(i + 64),
+		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 + 64),
-					   patch.vertex_handle(i + 65));
+					   patch.vertex_handle(i + patch_divs),
+					   patch.vertex_handle(i + patch_divs + 1));
 	}
 	return patch;
 }

src/smoothing.cpp

@@ -10,92 +10,131 @@
 
 
 /*
-   	  vα/--\
-   		/ 	 ------\   v
-   	   /  			---X
-   	  /			/---- /	\
-   	 /		/---	 /	 -\
-   	/  /----		/	   \
- vi	---			   |	 	\
-  	  \--		   /	 	 -\
-  	     \-		  /		 	/--\
-  		   \--	 /	/-------
-   	   	      \-+---
+      vα/--\
+        /    ------\   vi
+       /            ---X
+      /         /---- / \
+     /      /---     /   -\
+    /  /----        /      \
+ vj ---            |        \
+      \--          /         -\
+         \-       /         /--\
+           \--   /  /-------
+              \-+---
                vβ
 */
 
-Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
-	Point sum {0, 0, 0};
-	qreal area = 0;
-	Point p = mesh.point(vert);
-	qreal count = 0;
-	for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
-		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 pa = mesh.point(mesh.to_vertex_handle(v_va));
-		HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
-		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)) * (p - pi);
-		area += triangle_area(p, pi, pb);
+
+/* Poids pour le laplacien uniforme, constant à 1. */
+double uniform_weight(MyMesh &mesh, HalfedgeHandle vi_vj) {
+	(void) mesh;
+	(void) vi_vj;
+	return 1;
+}
+
+
+/* Masse pour le laplacien uniforme, 1 / taille(N1(vi)). */
+double uniform_mass(MyMesh &mesh, VertexHandle vi) {
+	double count = 0;
+	for (VertexHandle v : mesh.vv_range(vi)) {
+		(void) v;
 		count++;
 	}
-	area /= 3.;
-	return sum / (2.*count);
+	return 1. / count;
 }
 
 
-Eigen::SparseMatrix<qreal> laplacian_matrix(const MyMesh &mesh) {
-	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 p = mesh.point(vert);
-		for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
-			VertexHandle vi = mesh.to_vertex_handle(v_vi);
-			Point pi = mesh.point(vi);
-			HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
-			HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
-			Point pa = mesh.point(mesh.to_vertex_handle(v_va));
-			HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
-			Point pb = mesh.point(mesh.to_vertex_handle(vi_vb));
-			qreal a = angle_between(pi, pa, p);
-			qreal b = angle_between(pi, pb, p);
-			qreal w = -(cotan(a) + cotan(b)) / 2.;
-			sum += w;
-			cotangent.insert(vert.idx(), vi.idx()) = w;
-			area += triangle_area(p, pi, pb);
-			count++;
-		}
-		area /= 3.;
-		cotangent.insert(vert.idx(), vert.idx()) = -sum;
-		mass.insert(vert.idx(), vert.idx()) = 1. / area;
+/* Calcule le poids cotangent pour l'arête reliant vi à vj. */
+double cotangent_weight(MyMesh &mesh, HalfedgeHandle vi_vj) {
+	// Voir le graphique en haut de ce fichier.
+	Point pj = mesh.point(mesh.to_vertex_handle(vi_vj));
+	HalfedgeHandle vj_vb = mesh.next_halfedge_handle(vi_vj);
+	Point pb = mesh.point(mesh.to_vertex_handle(vj_vb));
+	HalfedgeHandle vj_vi = mesh.opposite_halfedge_handle(vi_vj);
+	Point pi = mesh.point(mesh.to_vertex_handle(vj_vi));
+	HalfedgeHandle vi_va = mesh.next_halfedge_handle(vj_vi);
+	Point pa = mesh.point(mesh.to_vertex_handle(vi_va));
+	double a = angle_between(pi, pa, pj);
+	double b = angle_between(pj, pb, pi);
+	return cotan(a) + cotan(b);
+}
+
+
+/* Calcule l'aire de chaque face et la stoque dans une propriété
+ * "face_area" du maillage. */
+void compute_face_areas(MyMesh &mesh) {
+	auto area = OpenMesh::getOrMakeProperty<FaceHandle, double>
+		(mesh, "face_area");
+	for (FaceHandle fh : mesh.faces()) {
+		MyMesh::FaceVertexIter fv_it = mesh.fv_iter(fh);
+		Point pi = mesh.point(*fv_it);
+		Point pj = mesh.point(*++fv_it);
+		Point pk = mesh.point(*++fv_it);
+		area[fh] = triangle_area(pi, pj, pk);
 	}
-	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]);
-	// 	}
-	// }
+/* Calcule l'aire barycentrique incidente à vert. */
+double barycentric_vertex_area(MyMesh &mesh, VertexHandle vert) {
+	auto area = OpenMesh::getProperty<FaceHandle, double>(mesh, "face_area");
+	double sum = 0;
+	for (FaceHandle fh : mesh.vf_range(vert)) {
+		sum += area[fh];
+	}
+	return sum / 3;
+}
 
-	// Approche matricielle
-	Eigen::SparseMatrix<qreal> laplacian = laplacian_matrix(mesh);
+
+/* Fonction de masse pour le laplacien cotangent. */
+double cotangent_mass(MyMesh &mesh, VertexHandle vi) {
+	double area = barycentric_vertex_area(mesh, vi);
+	return 1. / (2 * area);
+}
+
+
+/* Calcule la matrice laplacienne avec les fonctions de poids d'arête
+ * et de masse de sommet spécifiées par l'utilisatrice. */
+Eigen::SparseMatrix<qreal> laplacian_matrix(MyMesh &mesh, double (*edge_weight)(MyMesh &, HalfedgeHandle), double (*vertex_mass)(MyMesh &, VertexHandle)) {
+	compute_face_areas(mesh);
+	size_t n_verts = mesh.n_vertices();
+	Eigen::SparseMatrix<double> weight(n_verts, n_verts);
+	Eigen::SparseMatrix<double> mass(n_verts, n_verts);
+	for (VertexHandle vi : mesh.vertices()) {
+		if (mesh.is_boundary(vi)) continue;
+		double sum = 0;
+		for (HalfedgeHandle vi_vj : mesh.voh_range(vi)) {
+			// For each outgoing halfedge
+			VertexHandle vj = mesh.to_vertex_handle(vi_vj);
+			double halfedge_weight = edge_weight(mesh, vi_vj);
+			weight.insert(vi.idx(), vj.idx()) = halfedge_weight;
+			sum -= halfedge_weight;
+		}
+		// weight(i, i) = -Σ(j in N1(i)) weight(i, j)
+		weight.insert(vi.idx(), vi.idx()) = sum;
+
+		mass.insert(vi.idx(), vi.idx()) = vertex_mass(mesh, vi);
+	}
+	return mass * weight;
+}
+
+
+/* Adoucissement du maillage. */
+void smooth(MyMesh &mesh, SmoothingMethod method, double cotan_factor) {
+	// Calcul de la matrice laplacienne en fonction de la méthode choisie.
+	double factor;
+	Eigen::SparseMatrix<qreal> laplacian;
+	if (method == SmoothingMethod::COTANGENT) {
+		factor = cotan_factor;
+		laplacian = laplacian_matrix(mesh, cotangent_weight, cotangent_mass);
+	} else {
+		factor = 1;
+		laplacian = laplacian_matrix(mesh, uniform_weight, uniform_mass);
+	}
+
+	// laplacian = laplacian * laplacian; // Mise au carré.
+
+	// Transfert des coordonées de chaque point dans une matrice.
 	size_t n_verts = mesh.n_vertices();
 	Eigen::VectorX<qreal> X(n_verts), Y(n_verts), Z(n_verts);
 	for (VertexHandle vert : mesh.vertices()) {
@@ -106,13 +145,17 @@ void smooth(MyMesh &mesh) {
 		Y(id) = p[1];
 		Z(id) = p[2];
 	}
+
+	// Application du laplacien pour calculer les décalages.
 	X = laplacian * X;
 	Y = laplacian * Y;
 	Z = laplacian * Z;
+
+	// Application des décalages.
 	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);
+		mesh.set_point(vert, mesh.point(vert) + offset * factor);
 	}
 }

src/smoothing.h

@@ -4,8 +4,12 @@
 #include "my_mesh.h"
 
 
-Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert);
-void smooth(MyMesh &mesh);
+enum class SmoothingMethod {
+	UNIFORM,
+	COTANGENT,
+};
+
+void smooth(MyMesh &mesh, SmoothingMethod method, double cotan_factor=.0001);
 
 
 #endif

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

tp.md

@@ -0,0 +1,116 @@
+# 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.
+```cpp
+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.
+```cpp
+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](k.png)
+
+*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.
+
+```cpp
+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](patch.png)
+
+*Visualisation d'un patch*