CMakeLists.txt

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

data/bunnyLowPoly.obj

@@ -8128,7 +8128,6 @@ 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,37 +1,30 @@
+# 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
-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
+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 -1.0000 0.0000
-vn 1.0000 0.0000 0.0000
 vn 0.0000 0.0000 -1.0000
-usemtl Material
+vn 1.0000 0.0000 0.0000
+vn 0.0000 0.0000 1.0000
+vn 0.0000 -1.0000 0.0000
+vn 0.0000 1.0000 0.0000
 s off
-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
+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

src/curvature.cpp

@@ -1,265 +0,0 @@
-#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

@@ -1,163 +0,0 @@
-#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,8 +21,7 @@ int main(int argc, char *argv[]) {
 	QObject::connect(mesh_viewer, &MeshViewer::initialized,
 					 [&]() {
 						 if (mesh_processor) {
-							 mesh_viewer->addMesh(mesh_processor->mesh);
-							 mesh_viewer->addMesh(mesh_processor->patch);
+							 mesh_viewer->addMesh(*mesh_processor);
 						 }
 					 });
 
@@ -35,11 +34,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);
+							 mesh_viewer->removeMesh(*mesh_processor);
 							 delete mesh_processor;
 						 }
 						 mesh_processor = new MeshProcessor(path);
-						 mesh_viewer->addMesh(mesh_processor->mesh);
+						 mesh_viewer->addMesh(*mesh_processor);
 					 });
 	main_window.show();
 	return app.exec();

src/mesh_processor.cpp

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

src/mesh_processor.h

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

src/mesh_view.cpp

@@ -4,22 +4,21 @@
 #include <QOpenGLFunctions_2_1>
 
 
-MeshView::MeshView(const MyMesh &mesh, QOpenGLShaderProgram &program)
+MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program)
 	:vertex_buffer(QOpenGLBuffer::VertexBuffer),
 	 index_buffer(QOpenGLBuffer::IndexBuffer),
 	 holes_index_buffer(QOpenGLBuffer::IndexBuffer),
-	 n_faces(mesh.n_faces()),
-	 n_vertices(mesh.n_vertices()),
-	 mesh(mesh) {
-    QVector<GLfloat> vertices(n_vertices * 6);
+	 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);
 	size_t i = 0;
 	for (const VertexHandle it : mesh.vertices()) {
-		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];
+		vertices[3*i + 0] = mesh.point(it)[0];
+		vertices[3*i + 1] = mesh.point(it)[1];
+		vertices[3*i + 2] = mesh.point(it)[2];
 		i++;
 	}
 
@@ -37,7 +36,7 @@ MeshView::MeshView(const MyMesh &mesh, QOpenGLShaderProgram &program)
 
 	QVector<GLuint> holes_indices;
 	n_boundary_halfedges = 0;
-	for (std::vector<HalfedgeHandle> hole : mesh.holes) {
+	for (std::vector<HalfedgeHandle> hole : mesh_processor.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());
@@ -49,24 +48,18 @@ MeshView::MeshView(const MyMesh &mesh, QOpenGLShaderProgram &program)
 	{QOpenGLVertexArrayObject::Binder binder(&vao);
 		vertex_buffer.create();
 		vertex_buffer.bind();
-		vertex_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
-		vertex_buffer.allocate(vertices.constData(), n_vertices * 6 * sizeof (GLfloat));
+		vertex_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
+		vertex_buffer.allocate(vertices.constData(), n_vertices * 3 * sizeof (GLfloat));
 		index_buffer.create();
 		index_buffer.bind();
-		index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
+		index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
 		index_buffer.allocate(indices.constData(), n_faces * 3 * sizeof (GLuint));
 		holes_index_buffer.create();
 		holes_index_buffer.bind();
-		holes_index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
+		holes_index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
 		holes_index_buffer.allocate(holes_indices.constData(), n_boundary_halfedges * 2 * sizeof (GLuint));
-
-		vertex_buffer.bind();
-		program.setAttributeBuffer("pos", GL_FLOAT, 0,
-								   3, 6 * sizeof (GLfloat));
+		program.setAttributeBuffer("pos", GL_FLOAT, 0, 3);
 		program.enableAttributeArray("pos");
-		program.setAttributeBuffer("col", GL_FLOAT, 3 * sizeof (GLfloat),
-								   3, 6 * sizeof (GLfloat));
-		program.enableAttributeArray("col");
 	}
 }
 
@@ -79,6 +72,7 @@ MeshView::~MeshView() {
 
 
 void MeshView::paint(QOpenGLShaderProgram &program) {
+	const MyMesh &mesh = mesh_processor.mesh;
 	QOpenGLContext *ctx = QOpenGLContext::currentContext();
 	QOpenGLExtraFunctions *glf = ctx->extraFunctions();
 
@@ -88,7 +82,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 "my_mesh.h"
+#include "mesh_processor.h"
 #include <QOpenGLExtraFunctions>
 #include <QOpenGLShaderProgram>
 #include <QOpenGLVertexArrayObject>
@@ -18,9 +18,9 @@ class MeshView {
 	unsigned n_boundary_halfedges;
 
 public:
-	const MyMesh &mesh;
+	const MeshProcessor &mesh_processor;
 
-	MeshView(const MyMesh &mesh_processor, QOpenGLShaderProgram &program);
+	MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program);
 	~MeshView();
 	void paint(QOpenGLShaderProgram &program);
 };

src/mesh_viewer.cpp

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

src/mesh_viewer.h

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

src/my_mesh.h

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

src/quad_patch.cpp

@@ -1,46 +0,0 @@
-#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

@@ -1,11 +0,0 @@
-#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,5 +1,3 @@
-varying vec3 frag_col;
-
 uniform vec3 col;
 uniform vec3 wf_col;
 uniform bool wireframe;
@@ -9,5 +7,5 @@ void main() {
 	if (wireframe)
 		gl_FragColor = vec4(wf_col, alpha);
 	else
-		gl_FragColor = vec4(frag_col, alpha);
+		gl_FragColor = vec4(col, alpha);
 }

src/shader.vert

@@ -1,13 +1,9 @@
 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,7 +28,6 @@ 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);
@@ -38,81 +37,74 @@ 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)) * (p - pi);
+		sum += (cotan(a) + cotan(b)) * (pi - p);
 		area += triangle_area(p, pi, pb);
-		count++;
 	}
 	area /= 3.;
-	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;
+	return sum / (2.*area);
 }
 
 
 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
+	// 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]});
+	// }
 }

src/util.h

@@ -65,12 +65,4 @@ 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
+#endif