add MeshReconstruction dep and some test meshes

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,14 @@ 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.h
+  src/hole_filling.cpp)
 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

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,12 @@
+add_library(MeshReconstruction MeshReconstruction.h MeshReconstruction.cpp Cube.h Cube.cpp DataStructs.h IO.h IO.cpp Triangulation.h Triangulation.cpp)
+
+target_link_libraries(MeshReconstruction m)
+
+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,149 @@
+#include "Cube.h"
+#include <cmath>
+
+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
+	{
+		int 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 (abs(isoLevel - v1) < Eps) return p1;
+	if (abs(isoLevel - v2) < Eps) return p2;
+	if (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);
+
+	for (auto e = 0; e<12; ++e)
+	{
+		if (signConfigToIntersectedEdges[intersect.signConfig] & 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,53 @@
+#include "IO.h"
+#include <stdexcept>
+using namespace std;
+
+void MeshReconstruction::WriteObjFile(Mesh const& mesh, string const& fileName)
+{
+	// FILE faster than streams.
+
+	FILE* file;
+	//auto status = fopen_s(&file, fileName.c_str(), "w");
+	//if (status != 0)
+	//{
+	//	throw runtime_error("Could not write obj 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 (auto 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 (auto 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 (auto 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,78 @@
+#include "MeshReconstruction.h"
+#include "Cube.h"
+#include "Triangulation.h"
+
+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 HalfCubeDiag = cubeSize.Norm() / 2.0;
+	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;
+				auto dist = abs(sdf(cubeCenter) - isoLevel);
+				if (dist > HalfCubeDiag) 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

@@ -54,7 +54,9 @@ QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh) {
     static std::vector<MyMesh::VertexHandle> neigh;
 	neigh.clear();
 	get_two_neighborhood(neigh, vh);
-    if (neigh.size() < 5) throw "Quad fitting: not enough neighbors";
+    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,7 +14,10 @@ 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();

src/hole_filling.cpp

@@ -0,0 +1,58 @@
+#include "hole_filling.h"
+#include "util.h"
+#include <MeshReconstruction.h>
+
+
+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);
+	}
+}
+
+
+void fillHolesImplicit(MyMesh &mesh) {
+}

src/hole_filling.h

@@ -0,0 +1,13 @@
+#ifndef HOLE_FILLING_H
+#define HOLE_FILLING_H
+
+#include "my_mesh.h"
+#include <vector>
+
+
+void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole);
+void fillHolesDumb(MyMesh &mesh);
+void fillHolesImplicit(MyMesh &mesh);
+
+
+#endif

src/main.cpp

@@ -5,11 +5,27 @@
 #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);
+	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::smoothClicked,
+					 mesh_processor, &MeshProcessor::smooth);
+	QObject::connect(&main_window, &MainWindow::patchViewToggled,
+					 mesh_processor, &MeshProcessor::setPatchView);
+	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,32 +33,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);
-						 }
-					 });
-
 	if (argc > 2) {
 		qWarning("Utilisation : %s [MAILLAGE]", argv[0]);
 		return 1;
 	} else if (argc == 2) {
-		mesh_processor = new MeshProcessor(argv[1], *mesh_viewer);
-		QObject::connect(mesh_viewer, &MeshViewer::clicked,
-						 mesh_processor, &MeshProcessor::click);
+		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);
-						 mesh_viewer->addMesh(mesh_processor->mesh);
-						 QObject::connect(mesh_viewer, &MeshViewer::clicked,
-										  mesh_processor, &MeshProcessor::click);
+						 if (mesh_processor) delete mesh_processor;
+						 mesh_processor = create_mesh_processor
+							 (path, main_window);
 					 });
 	main_window.show();
 	return app.exec();

src/main_window.cpp

@@ -3,20 +3,21 @@
 
 #include <QApplication>
 #include <QFileDialog>
+#include <QMenuBar>
+#include <QGroupBox>
+#include <QPushButton>
 
 
 MainWindow::MainWindow(QWidget *parent)
 	:QMainWindow(parent),
 	 toolbar(this),
 	 mesh_viewer() {
-	connect(&mesh_viewer, &MeshViewer::initialized, [&]() {
-		open_action->setEnabled(true);
-	});
 	setCentralWidget(&mesh_viewer);
-	addToolBar(Qt::RightToolBarArea, &toolbar);
-	open_action = toolbar.addAction("Ouvrir…", [&](){
-		emit open(QFileDialog::getOpenFileName(this, "Ouvrir un maillage"));
-	});
+	// addToolBar(Qt::RightToolBarArea, &toolbar);
+
+	// open_action = toolbar.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);
@@ -28,8 +29,62 @@ MainWindow::MainWindow(QWidget *parent)
 	// 		mesh_viewer.addOpenGLMeshFromOpenMesh(&fragment, mat);
 	// 	}
 	// }));
-	open_action->setEnabled(false);
-	for (QAction *a : toolbar_actions) {
-		a->setEnabled(false);
+
+	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"));
+	});
+	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");
+	QLayout *hole_vbox = new QVBoxLayout();
+	hole_box->setLayout(hole_vbox);
+	QPushButton *fill_holes_dumb = new QPushButton("Remplir bêtement");
+	connect(fill_holes_dumb, &QPushButton::clicked,
+			this, &MainWindow::fillHolesDumbClicked);
+	hole_vbox->addWidget(fill_holes_dumb);
+	QPushButton *fill_holes_implicit = new QPushButton("Remplir bêtement");
+	connect(fill_holes_implicit, &QPushButton::clicked,
+			this, &MainWindow::fillHolesImplicitClicked);
+	hole_vbox->addWidget(fill_holes_implicit);
+	toolbar.addWidget(hole_box);
+
+	// Smoothing tools
+	QGroupBox *smooth_box = new QGroupBox("Adoucissement");
+	QLayout *smooth_vbox = new QVBoxLayout();
+	smooth_box->setLayout(smooth_vbox);
+	QPushButton *smooth = new QPushButton("Adoucir");
+	connect(smooth, &QPushButton::clicked,
+			this, &MainWindow::smoothClicked);
+	smooth_vbox->addWidget(smooth);
+	toolbar.addWidget(smooth_box);
+
+	// Curvature tools
+	QGroupBox *curvature_box = new QGroupBox("Analyse de courbure");
+	QLayout *curvature_vbox = new QVBoxLayout();
+	curvature_box->setLayout(curvature_vbox);
+	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_vbox->addWidget(patch_mode);
+	toolbar.addWidget(curvature_box);
 }

src/main_window.h

@@ -3,6 +3,7 @@
 
 #include <QMainWindow>
 #include <QToolBar>
+#include <QVBoxLayout>
 
 #include "mesh_viewer.h"
 #include "my_mesh.h"
@@ -13,10 +14,15 @@ class MainWindow : public QMainWindow {
 
 	QToolBar toolbar;
 	QAction *open_action;
-	QList<QAction *> toolbar_actions;
+	QAction *save_action;
 
 signals:
 	void open(const QString &path);
+	void save(const QString &path);
+	void fillHolesDumbClicked();
+	void fillHolesImplicitClicked();
+	void smoothClicked();
+	void patchViewToggled(bool checked);
 
 public:
 	MeshViewer mesh_viewer;

src/mesh_processor.cpp

@@ -1,31 +1,13 @@
 #include "quad_patch_tesselator.h"
 #include "mesh_processor.h"
 #include "util.h"
+#include "hole_filling.h"
 #include "smoothing.h"
 #include "curvature.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, MeshViewer &mesh_viewer)
 	:mesh_viewer(mesh_viewer) {
 	OpenMesh::IO::Options options;
@@ -34,52 +16,68 @@ MeshProcessor::MeshProcessor(const QString &path, MeshViewer &mesh_viewer)
 		qWarning() << "Failed to read" << path;
 		return;
 	}
-	mesh.update_normals();
+	for (VertexHandle vh : mesh.vertices()) {
+		mesh.set_color(vh, MyMesh::Color(.5, .5, .5));
+	}
 	courbure = new Courbures(mesh);
-	courbure->compute_KH();
-    courbure->set_K_colors();
-	// mesh.holes = findHoles(mesh);
-	// fillHoles();
-	// smooth(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();
 }
 
 
 MeshProcessor::~MeshProcessor() {
+	if (mesh_viewer.isInitialized()) {
+		mesh_viewer.removeMesh(mesh);
+	}
 	if (courbure) delete courbure;
 }
 
 
-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);
-
-	for (HalfedgeHandle it : hole) {
-		mesh.add_face(mesh.from_vertex_handle(it),
-					  mesh.to_vertex_handle(it),
-					  center_handle);
+void MeshProcessor::updateView() const {
+	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() {
+	::fillHolesImplicit(mesh);
+	updateView();
+}
+
+
+void MeshProcessor::smooth() {
+	::smooth(mesh);
+	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;

src/mesh_processor.h

@@ -3,7 +3,6 @@
 
 #include "my_mesh.h"
 #include "curvature.h"
-#include "mesh_processor_painter.h"
 #include "mesh_viewer.h"
 #include <vector>
 #include <QObject>
@@ -14,6 +13,9 @@ class MeshProcessor : public QObject {
 
 	Courbures *courbure = nullptr;
 	MeshViewer &mesh_viewer;
+	bool view_patches = false;
+
+	void updateView() const;
 
 public:
 	MyMesh mesh;
@@ -23,7 +25,10 @@ public:
 	~MeshProcessor();
 
 public slots:
-	void fillHoles();
+	void fillHolesDumb();
+	void fillHolesImplicit();
+	void smooth();
+	void setPatchView(bool on);
 	void click(QVector3D position);
 };
 

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

@@ -64,6 +64,7 @@ void MeshViewer::initializeGL() {
 	glf->glEnable(GL_MULTISAMPLE);
 
 	qDebug("MeshViewer: initialization complete");
+	is_initialized = true;
 	emit initialized();
 }
 

src/mesh_viewer.h

@@ -35,6 +35,7 @@ class MeshViewer : public QOpenGLWidget {
 	QMatrix4x4 view = trans * rot;
 	QPoint mouse_pos;
 	void updateViewMatrix();
+	bool is_initialized = false;
 
 public:
     MeshViewer(QWidget *parent=nullptr);
@@ -42,6 +43,7 @@ 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);

src/my_mesh.h

@@ -17,14 +17,14 @@ struct MyTraits : public OpenMesh::DefaultTraits {
     VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
     HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
     FaceAttributes(OpenMesh::Attributes::Normal);
-    EdgeAttributes(OpenMesh::Attributes::Color);
+    // EdgeAttributes(OpenMesh::Attributes::Color);
 };
 
 class MyMesh : public OpenMesh::TriMesh_ArrayKernelT<MyTraits> {
 public:
+	Color default_color {.5, .5, .5};
 	QMatrix4x4 transform;
 	std::vector<std::vector<HalfedgeHandle>> holes;
-	// QColor color = {127, 127, 127};
 };
 
 typedef MyMesh::FaceHandle FaceHandle;

src/smoothing.cpp

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