ccolin/mod_geo-tp
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

idk some changes, shit's still fucked

Browse Source
This commit is contained in:
papush! 2021-11-09 22:18:01 +01:00
parent 0c0cc8a907
commit 37fbe47aa8
3 changed files with 72 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/smoothing.cpp
View File

@ -3,6 +3,8 @@
#include "OpenMesh/Core/Mesh/SmartHandles.hh"
#include <OpenMesh/Core/Utils/PropertyManager.hh>
#include "util.h"
#include <Eigen/Sparse>
#include <Eigen/Core>
@ -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<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]});
// }
}

5
src/util.cpp
View File

@ -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);
}

19
src/util.h
View File

@ -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 <typename Point>
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 <typename Point>
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