From d01ba47b65d57e5220d1d7a88e3079bdfff75015 Mon Sep 17 00:00:00 2001 From: papush! Date: Fri, 12 Nov 2021 17:08:19 +0100 Subject: [PATCH] merge curvature stuff into this --- CMakeLists.txt | 6 +- data/bunnyLowPoly.obj | 1 + data/cube.obj | 59 +++++---- src/curvature.cpp | 265 +++++++++++++++++++++++++++++++++++++++++ src/curvature.h | 163 +++++++++++++++++++++++++ src/main.cpp | 7 +- src/mesh_processor.cpp | 15 ++- src/mesh_processor.h | 2 +- src/mesh_view.cpp | 42 ++++--- src/mesh_view.h | 6 +- src/mesh_viewer.cpp | 8 +- src/mesh_viewer.h | 6 +- src/my_mesh.h | 18 ++- src/quad_patch.cpp | 46 +++++++ src/quad_patch.h | 11 ++ src/shader.frag | 4 +- src/shader.vert | 4 + src/smoothing.cpp | 128 ++++++++++---------- src/util.h | 10 +- 19 files changed, 672 insertions(+), 129 deletions(-) create mode 100644 src/curvature.cpp create mode 100644 src/curvature.h create mode 100644 src/quad_patch.cpp create mode 100644 src/quad_patch.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 345b53c..69c69dc 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/data/bunnyLowPoly.obj b/data/bunnyLowPoly.obj index 826f74d..d5d81ee 100644 --- a/data/bunnyLowPoly.obj +++ b/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 diff --git a/data/cube.obj b/data/cube.obj index de4003f..17c9e47 100644 --- a/data/cube.obj +++ b/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 diff --git a/src/curvature.cpp b/src/curvature.cpp new file mode 100644 index 0000000..0405c3d --- /dev/null +++ b/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 Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) +{ + OpenMesh::VPropHandleT 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 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 Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) { +// OpenMesh::VPropHandleT 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 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 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 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 = 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 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 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(tmp, -sigma, +sigma); + tmp = (tmp + sigma) / (2 * sigma); + _mesh.set_color(vh, MyMesh::Color(tmp, .23, .5)); + } +} diff --git a/src/curvature.h b/src/curvature.h new file mode 100644 index 0000000..61632e9 --- /dev/null +++ b/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 +class MyStats { +private: + std::vector _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 vprop_K; + OpenMesh::VPropHandleT vprop_H; + OpenMesh::VPropHandleT vprop_quad; + + Courbures(MyMesh &mesh) : _mesh(mesh) {} + + void set_fixed_colors() ; + void normales_locales() ; + std::vector get_two_neighborhood(MyMesh::VertexHandle vh); + MyQuad fit_quad(MyMesh::VertexHandle vh) ; + void compute_KH() ; + void set_K_colors() ; +}; + + +#endif \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index d9eddf4..a0c9be1 100644 --- a/src/main.cpp +++ b/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(); diff --git a/src/mesh_processor.cpp b/src/mesh_processor.cpp index 889ae75..f6131ab 100644 --- a/src/mesh_processor.cpp +++ b/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 #include @@ -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 &hole) { } void MeshProcessor::fillHoles() { - for (auto hole : holes) { + for (auto hole : mesh.holes) { fillHoleDumb(mesh, hole); } } \ No newline at end of file diff --git a/src/mesh_processor.h b/src/mesh_processor.h index 489114d..dd045a8 100644 --- a/src/mesh_processor.h +++ b/src/mesh_processor.h @@ -12,7 +12,7 @@ class MeshProcessor : public QObject { public: MyMesh mesh; - std::vector> holes; + MyMesh patch; MeshProcessor(const QString &path); diff --git a/src/mesh_view.cpp b/src/mesh_view.cpp index 956874a..d857b6d 100644 --- a/src/mesh_view.cpp +++ b/src/mesh_view.cpp @@ -4,21 +4,22 @@ #include -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 vertices(n_vertices * 3); + n_faces(mesh.n_faces()), + n_vertices(mesh.n_vertices()), + mesh(mesh) { + QVector 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 holes_indices; n_boundary_halfedges = 0; - for (std::vector hole : mesh_processor.holes) { + for (std::vector 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 */ diff --git a/src/mesh_view.h b/src/mesh_view.h index 48f4bb8..7d8776a 100644 --- a/src/mesh_view.h +++ b/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 #include #include @@ -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); }; diff --git a/src/mesh_viewer.cpp b/src/mesh_viewer.cpp index ce45cfb..07105da 100644 --- a/src/mesh_viewer.cpp +++ b/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(); diff --git a/src/mesh_viewer.h b/src/mesh_viewer.h index 4aeef11..3c7a94f 100644 --- a/src/mesh_viewer.h +++ b/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 #include @@ -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); diff --git a/src/my_mesh.h b/src/my_mesh.h index b7512c5..48b35ba 100644 --- a/src/my_mesh.h +++ b/src/my_mesh.h @@ -10,11 +10,18 @@ #include +// struct MyTraits : public OpenMesh::DefaultTraits { +// using Point = Eigen::Vector3; +// using Normal = Eigen::Vector3; +// using Color = Eigen::Vector3; +// 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; - using Normal = Eigen::Vector3; - using Color = Eigen::Vector3; - 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 { public: QMatrix4x4 transform; - QColor color = {127, 127, 127}; + std::vector> holes; + // QColor color = {127, 127, 127}; }; typedef MyMesh::FaceHandle FaceHandle; diff --git a/src/quad_patch.cpp b/src/quad_patch.cpp new file mode 100644 index 0000000..ec2cc03 --- /dev/null +++ b/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 ch_base = q._r * q._t; + Eigen::Transform 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; +} diff --git a/src/quad_patch.h b/src/quad_patch.h new file mode 100644 index 0000000..097f487 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/shader.frag b/src/shader.frag index 0ccd914..a6c2a82 100644 --- a/src/shader.frag +++ b/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); } diff --git a/src/shader.vert b/src/shader.vert index f9c57af..164c6b7 100644 --- a/src/shader.vert +++ b/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); } diff --git a/src/smoothing.cpp b/src/smoothing.cpp index d776026..ec03dec 100644 --- a/src/smoothing.cpp +++ b/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 laplacian_matrix(const MyMesh &mesh) { + size_t n_verts = mesh.n_vertices(); + Eigen::SparseMatrix mass(n_verts, n_verts); + Eigen::SparseMatrix 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(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 mass(n_verts, n_verts); - // Eigen::SparseMatrix 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++; + // auto new_pos = OpenMesh::makeTemporaryProperty(mesh); + // for (VertexHandle v : mesh.vertices()) { + // if (!mesh.is_boundary(v)) { + // new_pos[v] = mesh.point(v) - laplace_beltrami(mesh, v); // } - // cotangent.insert(vert.idx(), vert.idx()) = -sum; - // mass.insert(vert.idx(), vert.idx()) = 1. / area; // } - // Eigen::SparseMatrix laplacian = mass * cotangent; - // Eigen::VectorX 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]}); + // for (VertexHandle v : mesh.vertices()) { + // if (!mesh.is_boundary(v)) { + // mesh.set_point(v, new_pos[v]); + // } // } + + // Approche matricielle + Eigen::SparseMatrix laplacian = laplacian_matrix(mesh); + size_t n_verts = mesh.n_vertices(); + Eigen::VectorX 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); + } } diff --git a/src/util.h b/src/util.h index 0dceb6b..068ebf7 100644 --- a/src/util.h +++ b/src/util.h @@ -65,4 +65,12 @@ qreal triangle_area(Point p1, Point p2, Point p3) { qreal cotan(const qreal x); -#endif \ No newline at end of file +template +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