From 37fbe47aa858644d017ae99e0373bd1762133aef Mon Sep 17 00:00:00 2001 From: papush! Date: Tue, 9 Nov 2021 22:18:01 +0100 Subject: [PATCH] idk some changes, shit's still fucked --- src/smoothing.cpp | 80 +++++++++++++++++++++++++++++++++++------------ src/util.cpp | 5 +++ src/util.h | 19 +++++------ 3 files changed, 72 insertions(+), 32 deletions(-) diff --git a/src/smoothing.cpp b/src/smoothing.cpp index f1eef94..d776026 100644 --- a/src/smoothing.cpp +++ b/src/smoothing.cpp @@ -3,6 +3,8 @@ #include "OpenMesh/Core/Mesh/SmartHandles.hh" #include #include "util.h" +#include +#include @@ -25,32 +27,21 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) { Point sum {0, 0, 0}; qreal area = 0; - size_t count = 0; - Point v = mesh.point(vert); + Point p = mesh.point(vert); for (HalfedgeHandle v_vi : mesh.voh_range(vert)) { - Point vi = mesh.point(mesh.to_vertex_handle(v_vi)); - + Point pi = mesh.point(mesh.to_vertex_handle(v_vi)); 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)); - + Point pa = 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); - sum += (1/qTan(a) + 1/qTan(b)) * (vi - v); - // sum += (vi - v); - area += triangle_area(v, vi, vb); - // qDebug("a: %lf, b: %lf, sum: %f %f %f, area: %lf", - // a, b, sum[0], sum[1], sum[2], area); - count++; + 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); + area += triangle_area(p, pi, pb); } area /= 3.; - // Point tmp(sum / (2*area)); - // qDebug() << tmp[0] << tmp[1] << tmp[2]; return sum / (2.*area); - // return sum / count; } @@ -59,7 +50,6 @@ void smooth(MyMesh &mesh) { for (VertexHandle v : mesh.vertices()) { if (!mesh.is_boundary(v)) { new_pos[v] = mesh.point(v) + laplace_beltrami(mesh, v); - // mesh.set_point(v, mesh.point(v) + laplace_beltrami(mesh, v)); } } for (VertexHandle v : mesh.vertices()) { @@ -67,4 +57,54 @@ void smooth(MyMesh &mesh) { 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++; + // } + // 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]}); + // } } diff --git a/src/util.cpp b/src/util.cpp index 0a51e1b..1600f9f 100644 --- a/src/util.cpp +++ b/src/util.cpp @@ -4,3 +4,8 @@ QDebug operator<<(QDebug dbg, const Point &p) { return dbg << p[0] << p[1] << p[2]; } + + +qreal cotan(const qreal x) { + return 1. / qTan(x); +} diff --git a/src/util.h b/src/util.h index e07a39a..0dceb6b 100644 --- a/src/util.h +++ b/src/util.h @@ -47,27 +47,22 @@ QDebug operator<<(QDebug dbg, const Point &p); /* Returns the angle between the vector from p2 to p1 and the vector * from p2 to p3. */ template -qreal angle_between(Point p1, - Point p2, - Point p3) { +qreal angle_between(Point p1, Point p2, Point p3) { Point vec_a = p1 - p2; Point vec_b = p3 - p2; - QVector3D a(vec_a[0], vec_a[1], vec_a[2]); - QVector3D b(vec_b[0], vec_b[1], vec_b[2]); - return qAcos(QVector3D::dotProduct(a.normalized(), b.normalized())); + return qAcos(vec_a.dot(vec_b) / (vec_a.norm() * vec_b.norm())); } template -qreal triangle_area(Point p1, - Point p2, - Point p3) { +qreal triangle_area(Point p1, Point p2, Point p3) { Point vec_a = p1 - p2; Point vec_b = p3 - p2; - QVector3D a(vec_a[0], vec_a[1], vec_a[2]); - QVector3D b(vec_b[0], vec_b[1], vec_b[2]); - return QVector3D::crossProduct(a, b).length() / 2.; + return vec_a.cross(vec_b).norm() / 2.; } +qreal cotan(const qreal x); + + #endif \ No newline at end of file