refactor quadratic patch code

CMakeLists.txt

@@ -26,7 +26,9 @@ target_sources(${PROJECT_NAME} PRIVATE
   src/curvature.cpp
   src/curvature.h
   src/quad_patch.cpp
-  src/quad_patch.h)
+  src/quad_patch.h
+  src/quad_patch_tesselator.cpp
+  src/quad_patch_tesselator.h)
 target_link_libraries(${PROJECT_NAME} PRIVATE
   Qt5::Core Qt5::Gui Qt5::Widgets
   OpenMeshCore

src/curvature.cpp

@@ -95,10 +95,8 @@ std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::
 // }
 
 
-MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh)
+QuadPatch Courbures::fit_quad(MyMesh::VertexHandle vh)
 {
-    MyQuad q ;
-
     std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh) ;
 
 //    std::cout << "-- Neighborhood" << std::endl ;
@@ -163,7 +161,7 @@ MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh)
         // Résolution aux moindres carrés par SVD
         Eigen::VectorXd coef(5) ;
         coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B) ;
-        MyQuad q(coef[0], coef[1], coef[2], coef[3], coef[4], r, t) ;
+        QuadPatch q(coef, ch_base) ;
 //        std::cout << "-- quad : " << std::endl ;
 //        std::cout << q[0] << ", " << q[1] << ", " << q[2] << ", " << q[3] << ", " << q[4] << ", " << q[5] << ", " << std::endl ;
         return q ;
@@ -221,7 +219,7 @@ void Courbures::compute_KH() {
 	_mesh.add_property(vprop_quad, "vprop_quad");
 
     for (VertexHandle vh : _mesh.vertices()) {
-        MyQuad q = fit_quad(vh);
+        QuadPatch q = fit_quad(vh);
 		_mesh.property(vprop_quad, vh) = q;
 
 		// K = det (K_P) = 4 a_0 a_2 - a_1 ^ 2

src/curvature.h

@@ -2,69 +2,7 @@
 #define CURVATURE_H
 
 #include "my_mesh.h"
-
+#include "quad_patch.h"
-
-class MyQuad {
-    double _coefs[5] ; // a_0 x^2 + a1 xy + a2 y^2 + a3 x + a4 y + a5
-public:
-    Eigen::AngleAxisd _r ;
-    Eigen::Translation3d _t ;
-
-	MyQuad(double *data,
-		   const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
-		   const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
-		: _r(r), _t(t)
-		{
-			for (size_t i=0; i<5; i++)
-				_coefs[i] = data[i] ;
-		}
-	MyQuad(const Eigen::VectorXd &v,
-		   const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
-		   const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
-		: _r(r), _t(t)
-		{
-			for (size_t i=0; i<5; i++)
-				_coefs[i] = v[i] ;
-		}
-	MyQuad(double a0=0, double a1=0, double a2=0, double a3=0, double a4=0,
-		   const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
-		   const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
-		: _r(r), _t(t)
-		{
-			_coefs[0] = a0 ;
-			_coefs[1] = a1 ;
-			_coefs[2] = a2 ;
-			_coefs[3] = a3 ;
-			_coefs[4] = a4 ;
-		}
-	MyQuad(const MyQuad &q) : _r(q._r), _t(q._t)
-		{
-			for (size_t i=0; i<5; i++)
-				_coefs[i] = q._coefs[i] ;
-		}
-
-	MyQuad &operator=(const MyQuad &q) {
-		_r = q._r;
-		_t = q._t;
-		for (size_t i=0; i<5; i++)
-			_coefs[i] = q._coefs[i] ;
-		return *this;
-	}
-
-    double & operator[] (size_t i) {return _coefs[i] ; }
-
-    double quad_fun (double x, double y)
-		{
-			return _coefs[0] * x*x + _coefs[1] * x*y + _coefs[2] * y*y + _coefs[3] * x + _coefs[4] * y  ;
-		}
-
-    double quad_fun (const OpenMesh::Vec2d &v)
-		{
-
-			return quad_fun(v[0], v[1]) ;
-		}
-
-};
 
 
 template <typename T>
@@ -147,14 +85,14 @@ private:
 public:
     OpenMesh::VPropHandleT<double> vprop_K;
     OpenMesh::VPropHandleT<double> vprop_H;
-	OpenMesh::VPropHandleT<MyQuad>    vprop_quad;
+	OpenMesh::VPropHandleT<QuadPatch> vprop_quad;
 
 	Courbures(MyMesh &mesh) : _mesh(mesh) {}
 
     void set_fixed_colors() ;
     void normales_locales() ;
     std::vector<MyMesh::VertexHandle> get_two_neighborhood(MyMesh::VertexHandle vh);
-    MyQuad fit_quad(MyMesh::VertexHandle vh) ;
+    QuadPatch fit_quad(MyMesh::VertexHandle vh);
     void compute_KH() ;
     void set_K_colors() ;
 };

src/mesh_processor.cpp

@@ -1,4 +1,4 @@
-#include "quad_patch.h"
+#include "quad_patch_tesselator.h"
 #include "mesh_processor.h"
 #include "util.h"
 #include "smoothing.h"
@@ -38,8 +38,8 @@ MeshProcessor::MeshProcessor(const QString &path) {
 	courb.compute_KH();
     courb.set_K_colors();
 	VertexHandle vh = mesh.vertex_handle(16);
-	MyQuad q = mesh.property(courb.vprop_quad, vh);
+	QuadPatch q = mesh.property(courb.vprop_quad, vh);
-	patch = quad_patch(q, mesh, vh);
+	patch = tesselate_quad_patch(q, mesh, vh);
 	// mesh.holes = findHoles(mesh);
 	// fillHoles();
 	// smooth(mesh);

src/my_quad.h

@@ -0,0 +1,33 @@
+#ifndef MY_QUAD_H
+#define MY_QUAD_H
+
+#include <Eigen/Geometry>
+
+
+class MyQuad {
+	typedef Eigen::Vector<double, 5> Vector5d;
+	Vector5d _coefficients; // a₀² + a₁xy + a₂y² + a₃x + a₄y + a₅
+    Eigen::AngleAxisd _rotation;
+    Eigen::Translation3d _translation;
+
+public:
+	MyQuad(Vector5d coefficients,
+		   const Eigen::AngleAxisd &rotation,
+		   const Eigen::Translation3d &translation)
+		:_coefficients(coefficients),
+		_rotation(rotation),
+		 _translation(translation) {}
+	const Vector5d &coefficients() { return _coefficients; }
+	const Eigen::AngleAxisd &rotation() { return _rotation; }
+	const Eigen::Translation3d &translation() { return _translation; }
+    double operator()(double x, double y) {
+		return _coefficients[0] * x*x
+			+ _coefficients[1] * x*y
+			+ _coefficients[2] * y*y
+			+ _coefficients[3] * x
+			+ _coefficients[4] * y;
+	}
+};
+
+
+#endif

src/quad_patch.cpp

@@ -1,46 +1,25 @@
 #include "quad_patch.h"
 
 
-MyMesh quad_patch(MyQuad q, MyMesh mesh, VertexHandle vh) {
+QuadPatch::QuadPatch() {}
-	MyMesh patch;
+
-	Eigen::Vector3d point
+
-		(mesh.point(vh)[0], mesh.point(vh)[1], mesh.point(vh)[2]);
+QuadPatch::QuadPatch(const Vector5d coefficients, const Transform transform)
-    Eigen::Transform<double, 3, Eigen::Affine> ch_base = q._r * q._t;
+		:_coefficients(coefficients),
-    Eigen::Transform<double, 3, Eigen::Affine> inv_ch_base = ch_base.inverse();
+		 _transform(transform),
-	point = ch_base * point;
+		 _inverse_transform(transform.inverse()) {
-	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 = ch_base * 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]);
 }
-	double xstep = (xmax-xmin)/64.;
+
-	double ystep = (ymax-ymin)/64.;
+
-	for (size_t y = 0; y < 64; y++) {
+double QuadPatch::operator()(double x, double y) {
-		for (size_t x = 0; x < 64; x++) {
+	return _coefficients[0] * x*x
-			double dx = xmin + x * xstep;
+		+ _coefficients[1] * x*y
-			double dy = ymin + y * ystep;
+		+ _coefficients[2] * y*y
-			Eigen::Vector3d point(dx, dy, -q.quad_fun(dx, dy));
+		+ _coefficients[3] * x
-			point = inv_ch_base * point;
+		+ _coefficients[4] * y;
-			patch.new_vertex(Point(point[0], point[1], point[2]));
 }
-	}
+
-	for (VertexHandle vhi : patch.vertices()) {
+
-		patch.set_color(vhi, MyMesh::Color(0, 1, .2));
+double QuadPatch::operator[](size_t i) {
-		size_t i = vhi.idx();
+	return _coefficients[i];
-		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),
-					   patch.vertex_handle(i + 1));
-		patch.add_face(patch.vertex_handle(i + 1),
-					   patch.vertex_handle(i + 64),
-					   patch.vertex_handle(i + 65));
-	}
-	return patch;
 }

src/quad_patch.h

@@ -1,11 +1,27 @@
-#ifndef QUAD_PATCH_H
+#ifndef MY_QUAD_H
-#define QUAD_PATCH_H
+#define MY_QUAD_H
 
-#include "curvature.h"
+#include <Eigen/Geometry>
-#include "my_mesh.h"
 
 
-MyMesh quad_patch(MyQuad q, MyMesh mesh, VertexHandle vh);
+class QuadPatch {
+public:
+	typedef Eigen::Vector<double, 5> Vector5d;
+	typedef Eigen::Transform<double, 3, Eigen::Affine> Transform;
+
+private:
+	Vector5d _coefficients; // a₀² + a₁xy + a₂y² + a₃x + a₄y + a₅
+	Transform _transform;
+	Transform _inverse_transform;
+
+public:
+	QuadPatch();
+	QuadPatch(const Vector5d coefficients, const Transform transform);
+    double operator()(double x, double y);
+	double operator[](size_t i);
+	constexpr Transform &transform() { return _transform; }
+	constexpr Transform &inverse_transform() { return _inverse_transform; }
+};
 
 
 #endif

src/quad_patch_tesselator.cpp

@@ -0,0 +1,44 @@
+#include "quad_patch_tesselator.h"
+
+
+MyMesh tesselate_quad_patch(QuadPatch q, MyMesh mesh, VertexHandle vh) {
+	MyMesh patch;
+	Eigen::Vector3d point
+		(mesh.point(vh)[0], mesh.point(vh)[1], mesh.point(vh)[2]);
+	point = q.transform() * point;
+	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]);
+	}
+	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++) {
+			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]));
+		}
+	}
+	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),
+					   patch.vertex_handle(i + 1));
+		patch.add_face(patch.vertex_handle(i + 1),
+					   patch.vertex_handle(i + 64),
+					   patch.vertex_handle(i + 65));
+	}
+	return patch;
+}

src/quad_patch_tesselator.h

@@ -0,0 +1,11 @@
+#ifndef QUAD_PATCH_H
+#define QUAD_PATCH_H
+
+#include "quad_patch.h"
+#include "my_mesh.h"
+
+
+MyMesh tesselate_quad_patch(QuadPatch q, MyMesh mesh, VertexHandle vh);
+
+
+#endif