merge curvature stuff into this

CMakeLists.txt

@@ -22,7 +22,11 @@ target_sources(${PROJECT_NAME} PRIVATE
   src/smoothing.cpp
   src/smoothing.h
   src/util.cpp
-  src/util.h)
+  src/util.h
+  src/curvature.cpp
+  src/curvature.h
+  src/quad_patch.cpp
+  src/quad_patch.h)
 target_link_libraries(${PROJECT_NAME} PRIVATE
   Qt5::Core Qt5::Gui Qt5::Widgets
   OpenMeshCore

data/bunnyLowPoly.obj

@@ -8128,6 +8128,7 @@ vn 0.8166 -0.5571 -0.1512
 vn 0.4879 -0.8093 -0.3271
 vn 0.4835 -0.8315 -0.2736
 vn 0.7111 -0.7029 -0.0163
+usemtl None
 s 1
 f 1//1 2//1 3//1
 f 4//2 5//2 6//2

data/cube.obj

@@ -1,30 +1,37 @@
-# Blender v2.90.1 OBJ File: ''
-# www.blender.org
 o Cube
-v -1.000000 -1.000000 1.000000
-v -1.000000 1.000000 1.000000
-v -1.000000 -1.000000 -1.000000
-v -1.000000 1.000000 -1.000000
-v 1.000000 -1.000000 1.000000
-v 1.000000 1.000000 1.000000
-v 1.000000 -1.000000 -1.000000
 v 1.000000 1.000000 -1.000000
-vn -1.0000 0.0000 0.0000
-vn 0.0000 0.0000 -1.0000
-vn 1.0000 0.0000 0.0000
-vn 0.0000 0.0000 1.0000
-vn 0.0000 -1.0000 0.0000
+v 1.000000 -1.000000 -1.000000
+v 1.000000 1.000000 1.000000
+v 1.000000 -1.000000 1.000000
+v -1.000000 1.000000 -1.000000
+v -1.000000 -1.000000 -1.000000
+v -1.000000 1.000000 1.000000
+v -1.000000 -1.000000 1.000000
+vt 0.625000 0.500000
+vt 0.875000 0.500000
+vt 0.875000 0.750000
+vt 0.625000 0.750000
+vt 0.375000 0.750000
+vt 0.625000 1.000000
+vt 0.375000 1.000000
+vt 0.375000 0.000000
+vt 0.625000 0.000000
+vt 0.625000 0.250000
+vt 0.375000 0.250000
+vt 0.125000 0.500000
+vt 0.375000 0.500000
+vt 0.125000 0.750000
 vn 0.0000 1.0000 0.0000
+vn 0.0000 0.0000 1.0000
+vn -1.0000 0.0000 0.0000
+vn 0.0000 -1.0000 0.0000
+vn 1.0000 0.0000 0.0000
+vn 0.0000 0.0000 -1.0000
+usemtl Material
 s off
-f 2//1 3//1 1//1
-f 4//2 7//2 3//2
-f 8//3 5//3 7//3
-f 6//4 1//4 5//4
-f 7//5 1//5 3//5
-f 4//6 6//6 8//6
-f 2//1 4//1 3//1
-f 4//2 8//2 7//2
-f 8//3 6//3 5//3
-f 6//4 2//4 1//4
-f 7//5 5//5 1//5
-f 4//6 2//6 6//6
+f 1/1/1 5/2/1 7/3/1 3/4/1
+f 4/5/2 3/4/2 7/6/2 8/7/2
+f 8/8/3 7/9/3 5/10/3 6/11/3
+f 6/12/4 2/13/4 4/5/4 8/14/4
+f 2/13/5 1/1/5 3/4/5 4/5/5
+f 6/11/6 5/10/6 1/1/6 2/13/6

src/curvature.cpp

@@ -0,0 +1,265 @@
+#include "curvature.h"
+#include "util.h"
+
+
+void Courbures::normales_locales() {
+    _mesh.request_vertex_normals();
+    size_t i;
+    for (VertexHandle vh : _mesh.vertices()) {
+		MyMesh::Normal normal(0,0,0);
+        i = 0;
+		for (MyMesh::FaceHandle vf : _mesh.vf_range(vh)) {
+            i++ ;
+            normal += _mesh.calc_face_normal(vf);
+        }
+        if (i != 0) {
+            normal /= i;
+            normal /= norm(normal);
+        }
+        _mesh.set_normal(vh, normal);
+    }
+}
+
+
+std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh)
+{
+    OpenMesh::VPropHandleT<bool>    vprop_flag;
+    _mesh.add_property(vprop_flag,  "vprop_flag");
+
+    // Initialisation
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end(); ++v_it)
+    {
+        _mesh.property(vprop_flag, *v_it) = false ;
+    }
+    // Circulateur sur le premier cercle
+    std::vector<MyMesh::VertexHandle> neigh, neigh2 ;
+
+	_mesh.property(vprop_flag, vh) = true ;
+    for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it)
+    {
+        neigh.push_back(*vv_it) ; // ajout du point à la liste
+        _mesh.property(vprop_flag, *vv_it) = true ;
+    }
+
+    // Parcours du premier cercle et ajout du second cercle par circulateurs
+    for (int i=0; i<neigh.size(); i++)
+    {
+        MyMesh::VertexHandle vh = neigh.at(i) ;
+        for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it)
+        {
+            if (!_mesh.property(vprop_flag, *vv_it)) // sommet non encore rencontré
+            neigh2.push_back(*vv_it) ; // ajout du point à la liste
+            _mesh.property(vprop_flag, *vv_it) = true ;
+        }
+    }
+
+    // Concaténation des deux cercles
+    neigh.insert(neigh.end(), neigh2.begin(), neigh2.end()) ;
+
+    _mesh.remove_property(vprop_flag);
+
+    return neigh ;
+}
+// std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) {
+//     OpenMesh::VPropHandleT<bool> vprop_flag;
+//     _mesh.add_property(vprop_flag, "vprop_flag");
+
+//     // 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;
+//     for(VertexHandle vv : _mesh.vv_range(vh)) {
+//         neigh.push_back(vv); // ajout du point à la liste
+//         _mesh.property(vprop_flag, 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);
+// 			}
+//             _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;
+// }
+
+
+MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh)
+{
+    MyQuad q ;
+
+    std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh) ;
+
+//    std::cout << "-- Neighborhood" << std::endl ;
+    for (int i = 0; i< neigh.size(); i++)
+    {
+        MyMesh::Point p = _mesh.point(neigh.at(i)) ;
+//        std::cout << p[0] << ", " << p[1] << ", " << p[2] << " ; " << std::endl ;
+    }
+
+    // Calcul de la matrice de changement de base
+    MyMesh::Normal n = _mesh.normal(vh) ;
+    MyMesh::Point p = _mesh.point(vh);
+//    std::cout << "-- normale" << std::endl ;
+//    std::cout << n[0] << ", " << n[1] << ", " << n[2] << std::endl ;
+//    std::cout << "-- point" << std::endl ;
+//    std::cout << p[0] << ", " << p[1] << ", " << p[2] << std::endl ;
+
+    Eigen::Vector3d ne(n[0], n[1], n[2]), Oz(0,0,1), axis;
+    Eigen::Vector3d p_e(p[0], p[1], p[2]), pi_e ;
+
+    axis = ne.cross(Oz) ;
+    double sina = axis.norm(), cosa = ne.dot(Oz), angle ;
+    if (sina >= 0)
+        angle = acos(cosa) ;
+    else
+        angle = -acos(cosa) ;
+    axis = axis.normalized() ;
+
+    Eigen::AngleAxisd r(angle, axis) ;
+//    std::cout << "-- rotation" << std:: endl ;
+//    std::cout << r.matrix()(0,0) << ", " << r.matrix()(0, 1) << ", " << r.matrix()(0,2) << "; " ;
+//    std::cout << r.matrix()(1,0) << ", " << r.matrix()(1, 1) << ", " << r.matrix()(1,2) << "; " ;
+//    std::cout << r.matrix()(2,0) << ", " << r.matrix()(2, 1) << ", " << r.matrix()(2,2) << std::endl ;
+    Eigen::Translation3d t(-p_e[0], -p_e[1], -p_e[2]) ;
+    Eigen::Transform<double, 3, Eigen::Affine> ch_base = r * t ;
+
+
+    // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen)
+
+    if (neigh.size() >= 5)
+    {
+        int n(neigh.size()) ;
+        Eigen::MatrixXd A(n,5);
+        Eigen::VectorXd B(n);
+
+//        std::cout << "-- Après changement de base" << std::endl ;
+        for(int i=0; i<neigh.size(); i++)
+        {
+            MyMesh::Point p = _mesh.point(neigh.at(i)) ;
+            pi_e << p[0], p[1], p[2] ;
+            pi_e = ch_base * pi_e ; // Application du changement de base
+//            std::cout << pi_e[0] << ", " << pi_e[1] << ", " << pi_e[2] << std::endl ;
+
+			B[i] = -pi_e[2]; // -zi
+			A(i, 0) = pi_e[0] * pi_e[0]; // xi^2
+			A(i, 1) = pi_e[0] * pi_e[1]; // xi*yi
+			A(i, 2) = pi_e[1] * pi_e[1]; // yi^2
+			A(i, 3) = pi_e[0]; // xi
+			A(i, 4) = pi_e[1]; // yi
+        }
+
+        // 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) ;
+//        std::cout << "-- quad : " << std::endl ;
+//        std::cout << q[0] << ", " << q[1] << ", " << q[2] << ", " << q[3] << ", " << q[4] << ", " << q[5] << ", " << std::endl ;
+        return q ;
+    }
+    else
+    {
+        std::cout << "Quad fitting : not enough neighbors" ;
+        throw "Quad fitting : not enough neighbors" ;
+    }
+}
+// MyQuad 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";
+
+//     // Calcul de la matrice de changement de base
+// 	Eigen::Vector3d Oz(0,0,1);
+// 	Eigen::Vector3d axis = _mesh.normal(vh).cross(Oz);
+//     double sina = axis.norm();
+// 	double cosa = _mesh.normal(vh).dot(Oz);
+// 	double angle;
+//     if (sina >= 0) {
+//         angle = acos(cosa);
+// 	} else {
+//         angle = -acos(cosa);
+// 	}
+//     axis = axis.normalized();
+//     Eigen::AngleAxisd rot(angle, axis);
+//     Eigen::Translation3d trans(-_mesh.point(vh));
+//     Eigen::Transform<double, 3, Eigen::Affine> ch_base = rot * trans;
+
+//     // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen)
+// 	Eigen::MatrixXd A(neigh.size(), 5);
+// 	Eigen::VectorXd B(neigh.size());
+// 	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 MyQuad(coef[0], coef[1], coef[2], coef[3], coef[4], rot, trans);
+// }
+
+
+void Courbures::compute_KH() {
+    _mesh.add_property(vprop_K,  "vprop_K");
+    _mesh.add_property(vprop_H,  "vprop_H");
+	_mesh.add_property(vprop_quad, "vprop_quad");
+
+    for (VertexHandle vh : _mesh.vertices()) {
+        MyQuad q = fit_quad(vh);
+		_mesh.property(vprop_quad, vh) = q;
+
+		// K = det (K_P) = 4 a_0 a_2 - a_1 ^ 2
+        _mesh.property(vprop_K, vh) = 4 * q[0] * q[2] - q[1]*q[1];
+
+		// H = Tr (K_P) = a_0 + a_2
+        _mesh.property(vprop_H, vh) = q[0] + q[2];
+    }
+}
+
+
+void Courbures::set_fixed_colors() {
+    size_t i = 0;
+    for (VertexHandle vh : _mesh.vertices()) {
+        if (i % 2)
+            _mesh.set_color(vh, MyMesh::Color(0, 1, 0));
+        else
+            _mesh.set_color(vh, MyMesh::Color(0, 0, 1));
+		i++;
+    }
+}
+
+
+void Courbures::set_K_colors() {
+    MyStats<double> dist;
+    for (VertexHandle vh : _mesh.vertices()) {
+        dist.push_back(_mesh.property(vprop_K, vh));
+    }
+    double m = dist.mean();
+	double sigma = dist.stdev(m);
+    std::cout << "K min : " << dist.min()
+			  << " - max : " << dist.max() << std::endl ;
+    std::cout << "K mean : " << m << " - sigma : " << sigma << std::endl;
+
+    for (VertexHandle vh : _mesh.vertices()) {
+        double tmp  = _mesh.property(vprop_K, vh) - m;
+		tmp = clamp<double>(tmp, -sigma, +sigma);
+		tmp = (tmp + sigma) / (2 * sigma);
+		_mesh.set_color(vh, MyMesh::Color(tmp, .23, .5));
+    }
+}

src/curvature.h

@@ -0,0 +1,163 @@
+#ifndef CURVATURE_H
+#define CURVATURE_H
+
+#include "my_mesh.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>
+class MyStats {
+private:
+    std::vector<T> _distrib ;
+public:
+    MyStats () {} ;
+
+    void push_back (T data)
+		{
+			_distrib.push_back(data) ;
+		}
+
+    T min ()
+		{
+			T tmp = _distrib.at(0) ;
+			for (size_t i=1 ; i<_distrib.size(); i++)
+				{
+					if (_distrib.at(i) < tmp)
+						tmp = _distrib.at(i) ;
+				}
+			return tmp ;
+		}
+
+    T max ()
+		{
+			T tmp = _distrib.at(0) ;
+			for (size_t i=1 ; i<_distrib.size(); i++)
+				{
+					if (_distrib.at(i) > tmp)
+						tmp = _distrib.at(i) ;
+				}
+			return tmp ;
+		}
+
+    T mean ()
+		{
+			T acc(0) ;
+			std::cout << "acc : " << acc << std::endl;
+			if (_distrib.size() > 0)
+				{
+					for(size_t i=0; i<_distrib.size(); i++)
+						{
+							acc += _distrib.at(i) ;
+						}
+					return acc/(_distrib.size()) ;
+				}
+			else
+				return acc ;
+		}
+
+    T stdev ()
+		{
+			return stdev(mean());
+		}
+
+    T stdev (T m)
+		{
+			T acc(0), tmp ;
+			if (_distrib.size() > 0)
+				{
+					for(size_t i=0; i<_distrib.size(); i++)
+						{
+							tmp = _distrib.at(i) - m ;
+							acc += tmp * tmp ;
+						}
+					return sqrt(acc/(_distrib.size())) ;
+				}
+			else
+				return acc ;
+		}
+};
+
+
+class Courbures {
+private:
+    MyMesh &_mesh ;
+
+public:
+    OpenMesh::VPropHandleT<double>    vprop_K;
+    OpenMesh::VPropHandleT<double>    vprop_H;
+	OpenMesh::VPropHandleT<MyQuad>    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) ;
+    void compute_KH() ;
+    void set_K_colors() ;
+};
+
+
+#endif

src/main.cpp

@@ -21,7 +21,8 @@ int main(int argc, char *argv[]) {
 	QObject::connect(mesh_viewer, &MeshViewer::initialized,
 					 [&]() {
 						 if (mesh_processor) {
-							 mesh_viewer->addMesh(*mesh_processor);
+							 mesh_viewer->addMesh(mesh_processor->mesh);
+							 mesh_viewer->addMesh(mesh_processor->patch);
 						 }
 					 });
 
@@ -34,11 +35,11 @@ int main(int argc, char *argv[]) {
 	QObject::connect(&main_window, &MainWindow::open,
 					 [&](const QString &path) {
 						 if (mesh_processor) {
-							 mesh_viewer->removeMesh(*mesh_processor);
+							 mesh_viewer->removeMesh(mesh_processor->mesh);
 							 delete mesh_processor;
 						 }
 						 mesh_processor = new MeshProcessor(path);
-						 mesh_viewer->addMesh(*mesh_processor);
+						 mesh_viewer->addMesh(mesh_processor->mesh);
 					 });
 	main_window.show();
 	return app.exec();

src/mesh_processor.cpp

@@ -1,6 +1,8 @@
+#include "quad_patch.h"
 #include "mesh_processor.h"
 #include "util.h"
 #include "smoothing.h"
+#include "curvature.h"
 #include <QGenericMatrix>
 #include <unordered_set>
 
@@ -31,9 +33,16 @@ MeshProcessor::MeshProcessor(const QString &path) {
 		qWarning() << "Failed to read" << path;
 		return;
 	}
-	// holes = findHoles(mesh);
+	mesh.update_normals();
+	Courbures courb(mesh);
+	courb.compute_KH();
+    courb.set_K_colors();
+	VertexHandle vh = mesh.vertex_handle(16);
+	MyQuad q = mesh.property(courb.vprop_quad, vh);
+	patch = quad_patch(q, mesh, vh);
+	// mesh.holes = findHoles(mesh);
 	// fillHoles();
-	smooth(mesh);
+	// smooth(mesh);
 }
 
 
@@ -60,7 +69,7 @@ void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole) {
 }
 
 void MeshProcessor::fillHoles() {
-	for (auto hole : holes) {
+	for (auto hole : mesh.holes) {
 		fillHoleDumb(mesh, hole);
 	}
 }

src/mesh_processor.h

@@ -12,7 +12,7 @@ class MeshProcessor : public QObject {
 
 public:
 	MyMesh mesh;
-	std::vector<std::vector<HalfedgeHandle>> holes;
+	MyMesh patch;
 
 	MeshProcessor(const QString &path);
 

src/mesh_view.cpp

@@ -4,21 +4,22 @@
 #include <QOpenGLFunctions_2_1>
 
 
-MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program)
+MeshView::MeshView(const MyMesh &mesh, QOpenGLShaderProgram &program)
 	:vertex_buffer(QOpenGLBuffer::VertexBuffer),
 	 index_buffer(QOpenGLBuffer::IndexBuffer),
 	 holes_index_buffer(QOpenGLBuffer::IndexBuffer),
-	 n_faces(mesh_processor.mesh.n_faces()),
-	 n_vertices(mesh_processor.mesh.n_vertices()),
-	 mesh_processor(mesh_processor) {
-	const MyMesh &mesh = mesh_processor.mesh;
-
-    QVector<GLfloat> vertices(n_vertices * 3);
+	 n_faces(mesh.n_faces()),
+	 n_vertices(mesh.n_vertices()),
+	 mesh(mesh) {
+    QVector<GLfloat> vertices(n_vertices * 6);
 	size_t i = 0;
 	for (const VertexHandle it : mesh.vertices()) {
-		vertices[3*i + 0] = mesh.point(it)[0];
-		vertices[3*i + 1] = mesh.point(it)[1];
-		vertices[3*i + 2] = mesh.point(it)[2];
+		vertices[6*i + 0] = mesh.point(it)[0];
+		vertices[6*i + 1] = mesh.point(it)[1];
+		vertices[6*i + 2] = mesh.point(it)[2];
+		vertices[6*i + 3] = mesh.color(it)[0];
+		vertices[6*i + 4] = mesh.color(it)[1];
+		vertices[6*i + 5] = mesh.color(it)[2];
 		i++;
 	}
 
@@ -36,7 +37,7 @@ MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &pr
 
 	QVector<GLuint> holes_indices;
 	n_boundary_halfedges = 0;
-	for (std::vector<HalfedgeHandle> hole : mesh_processor.holes) {
+	for (std::vector<HalfedgeHandle> hole : mesh.holes) {
 		for (HalfedgeHandle it : hole) {
 			holes_indices.push_back(mesh.from_vertex_handle(it).idx());
 			holes_indices.push_back(mesh.to_vertex_handle(it).idx());
@@ -48,18 +49,24 @@ MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &pr
 	{QOpenGLVertexArrayObject::Binder binder(&vao);
 		vertex_buffer.create();
 		vertex_buffer.bind();
-		vertex_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
-		vertex_buffer.allocate(vertices.constData(), n_vertices * 3 * sizeof (GLfloat));
+		vertex_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
+		vertex_buffer.allocate(vertices.constData(), n_vertices * 6 * sizeof (GLfloat));
 		index_buffer.create();
 		index_buffer.bind();
-		index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
+		index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
 		index_buffer.allocate(indices.constData(), n_faces * 3 * sizeof (GLuint));
 		holes_index_buffer.create();
 		holes_index_buffer.bind();
-		holes_index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
+		holes_index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
 		holes_index_buffer.allocate(holes_indices.constData(), n_boundary_halfedges * 2 * sizeof (GLuint));
-		program.setAttributeBuffer("pos", GL_FLOAT, 0, 3);
+
+		vertex_buffer.bind();
+		program.setAttributeBuffer("pos", GL_FLOAT, 0,
+								   3, 6 * sizeof (GLfloat));
 		program.enableAttributeArray("pos");
+		program.setAttributeBuffer("col", GL_FLOAT, 3 * sizeof (GLfloat),
+								   3, 6 * sizeof (GLfloat));
+		program.enableAttributeArray("col");
 	}
 }
 
@@ -72,7 +79,6 @@ MeshView::~MeshView() {
 
 
 void MeshView::paint(QOpenGLShaderProgram &program) {
-	const MyMesh &mesh = mesh_processor.mesh;
 	QOpenGLContext *ctx = QOpenGLContext::currentContext();
 	QOpenGLExtraFunctions *glf = ctx->extraFunctions();
 
@@ -82,7 +88,7 @@ void MeshView::paint(QOpenGLShaderProgram &program) {
 		qFatal("Failed to get OpenGL 2.1 functions");
 
 	program.setUniformValue("model", mesh.transform);
-	program.setUniformValue("col", mesh.color.redF(), mesh.color.greenF(), mesh.color.blueF());
+	// program.setUniformValue("col", mesh.color.redF(), mesh.color.greenF(), mesh.color.blueF());
 
 	{QOpenGLVertexArrayObject::Binder binder(&vao);
 		/* Mesh */

src/mesh_view.h

@@ -1,7 +1,7 @@
 #ifndef MESH_VIEW_H
 #define MESH_VIEW_H
 
-#include "mesh_processor.h"
+#include "my_mesh.h"
 #include <QOpenGLExtraFunctions>
 #include <QOpenGLShaderProgram>
 #include <QOpenGLVertexArrayObject>
@@ -18,9 +18,9 @@ class MeshView {
 	unsigned n_boundary_halfedges;
 
 public:
-	const MeshProcessor &mesh_processor;
+	const MyMesh &mesh;
 
-	MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program);
+	MeshView(const MyMesh &mesh_processor, QOpenGLShaderProgram &program);
 	~MeshView();
 	void paint(QOpenGLShaderProgram &program);
 };

src/mesh_viewer.cpp

@@ -83,19 +83,19 @@ void MeshViewer::paintGL() {
 }
 
 
-void MeshViewer::addMesh(const MeshProcessor &mesh_processor) {
+void MeshViewer::addMesh(const MyMesh &mesh) {
 	Q_ASSERT(isValid());
 	makeCurrent();
-	meshes.emplace_back(mesh_processor, program);
+	meshes.emplace_back(mesh, program);
 	doneCurrent();
 	update();
 }
 
 
-void MeshViewer::removeMesh(const MeshProcessor &mesh_processor) {
+void MeshViewer::removeMesh(const MyMesh &mesh) {
 	makeCurrent();
 	meshes.remove_if([&](const MeshView &mv) {
-		return &mv.mesh_processor == &mesh_processor;
+		return &mv.mesh == &mesh;
 	});
 	doneCurrent();
 	update();

src/mesh_viewer.h

@@ -2,7 +2,7 @@
 #define MESH_VIEWER_H
 
 #define GL_GLEXT_PROTOTYPES
-#include "mesh_processor.h"
+#include "my_mesh.h"
 #include "mesh_view.h"
 #include <list>
 #include <QMatrix4x4>
@@ -39,8 +39,8 @@ public:
     void paintGL() override;
 
 public slots:
-	void addMesh(const MeshProcessor &mesh);
-	void removeMesh(const MeshProcessor &mesh);
+	void addMesh(const MyMesh &mesh);
+	void removeMesh(const MyMesh &mesh);
 
 protected:
     virtual void mousePressEvent(QMouseEvent *e);

src/my_mesh.h

@@ -10,11 +10,18 @@
 #include <OpenMesh/Core/Geometry/EigenVectorT.hh>
 
 
+// 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);
+//     HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
+//     FaceAttributes(OpenMesh::Attributes::Normal);
+//     EdgeAttributes(OpenMesh::Attributes::Color);
+// };
 struct MyTraits : public OpenMesh::DefaultTraits {
-    using Point = Eigen::Vector3<qreal>;
-    using Normal = Eigen::Vector3<qreal>;
-    using Color = Eigen::Vector3<qreal>;
-    VertexAttributes(OpenMesh::Attributes::Normal);
+	typedef OpenMesh::Vec3f Color;
+    VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
     HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
     FaceAttributes(OpenMesh::Attributes::Normal);
     EdgeAttributes(OpenMesh::Attributes::Color);
@@ -23,7 +30,8 @@ struct MyTraits : public OpenMesh::DefaultTraits {
 class MyMesh : public OpenMesh::TriMesh_ArrayKernelT<MyTraits> {
 public:
 	QMatrix4x4 transform;
-	QColor color = {127, 127, 127};
+	std::vector<std::vector<HalfedgeHandle>> holes;
+	// QColor color = {127, 127, 127};
 };
 
 typedef MyMesh::FaceHandle FaceHandle;

src/quad_patch.cpp

@@ -0,0 +1,46 @@
+#include "quad_patch.h"
+
+
+MyMesh quad_patch(MyQuad q, MyMesh mesh, VertexHandle vh) {
+	MyMesh patch;
+	Eigen::Vector3d point
+		(mesh.point(vh)[0], mesh.point(vh)[1], mesh.point(vh)[2]);
+    Eigen::Transform<double, 3, Eigen::Affine> ch_base = q._r * q._t;
+    Eigen::Transform<double, 3, Eigen::Affine> inv_ch_base = ch_base.inverse();
+	point = ch_base * 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 = 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++) {
+		for (size_t x = 0; x < 64; x++) {
+			double dx = xmin + x * xstep;
+			double dy = ymin + y * ystep;
+			Eigen::Vector3d point(dx, dy, -q.quad_fun(dx, dy));
+			point = inv_ch_base * 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.h

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

src/shader.frag

@@ -1,3 +1,5 @@
+varying vec3 frag_col;
+
 uniform vec3 col;
 uniform vec3 wf_col;
 uniform bool wireframe;
@@ -7,5 +9,5 @@ void main() {
 	if (wireframe)
 		gl_FragColor = vec4(wf_col, alpha);
 	else
-		gl_FragColor = vec4(col, alpha);
+		gl_FragColor = vec4(frag_col, alpha);
 }

src/shader.vert

@@ -1,9 +1,13 @@
 attribute vec3 pos;
+attribute vec3 col;
+
+varying vec3 frag_col;
 
 uniform mat4 proj;
 uniform mat4 view;
 uniform mat4 model;
 
 void main() {
+	frag_col = col;
 	gl_Position = proj * view * model * vec4(pos, 1.0);
 }

src/smoothing.cpp

@@ -28,6 +28,7 @@ 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);
@@ -37,74 +38,81 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
 		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)) * (pi - p);
+		sum += (cotan(a) + cotan(b)) * (p - pi);
 		area += triangle_area(p, pi, pb);
+		count++;
 	}
 	area /= 3.;
-	return sum / (2.*area);
+	return sum / (2.*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;
+	}
+	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
-	// 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 v = mesh.point(vert);
-	// 	for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
-	// 		VertexHandle vi_h = mesh.to_vertex_handle(v_vi);
-	// 		Point vi = mesh.point(vi_h);
-
-	// 		HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
-	// 		HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
-	// 		Point va = mesh.point(mesh.to_vertex_handle(v_va));
-
-	// 		HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
-	// 		Point vb = mesh.point(mesh.to_vertex_handle(vi_vb));
-	// 		qreal a = angle_between(vi, va, v);
-	// 		qreal b = angle_between(v, vb, vi);
-	// 		qreal w = -(1./qTan(a) + 1./qTan(b)) / 2.;
-	// 		sum += w;
-	// 		cotangent.insert(vert.idx(), vi_h.idx()) = w;
-	// 		area += triangle_area(v, vi, vb);
-	// 		count++;
-	// 	}
-	// 	cotangent.insert(vert.idx(), vert.idx()) = -sum;
-	// 	mass.insert(vert.idx(), vert.idx()) = 1. / area;
-	// }
-	// Eigen::SparseMatrix<qreal> laplacian = mass * cotangent;
-	// 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();
-	// 	mesh.set_point(vert, {X[id], Y[id], Z[id]});
-	// }
+	// 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);
+	}
 }

src/util.h

@@ -65,4 +65,12 @@ qreal triangle_area(Point p1, Point p2, Point p3) {
 qreal cotan(const qreal x);
 
 
+template <typename T>
+constexpr const T& clamp(const T& v, const T& lo, const T& hi) {
+	if (v < lo) return lo;
+	if (v > hi) return hi;
+	return v;
+}
+
+
 #endif