143 lines
3.9 KiB
C++
143 lines
3.9 KiB
C++
#include "smoothing.h"
|
||
|
||
#include "OpenMesh/Core/Mesh/SmartHandles.hh"
|
||
#include <OpenMesh/Core/Utils/PropertyManager.hh>
|
||
#include "util.h"
|
||
#include <Eigen/Sparse>
|
||
#include <Eigen/Core>
|
||
|
||
|
||
|
||
|
||
/*
|
||
vα/--\
|
||
/ ------\ vi
|
||
/ ---X
|
||
/ /---- / \
|
||
/ /--- / -\
|
||
/ /---- / \
|
||
vj --- | \
|
||
\-- / -\
|
||
\- / /--\
|
||
\-- / /-------
|
||
\-+---
|
||
vβ
|
||
*/
|
||
|
||
|
||
double uniform_weight(MyMesh &mesh, HalfedgeHandle vi_vj) {
|
||
(void) mesh;
|
||
(void) vi_vj;
|
||
return 1;
|
||
}
|
||
|
||
double uniform_mass(MyMesh &mesh, VertexHandle vi) {
|
||
double count = 0;
|
||
for (VertexHandle v : mesh.vv_range(vi)) {
|
||
(void) v;
|
||
count++;
|
||
}
|
||
return 1. / count;
|
||
}
|
||
|
||
|
||
/* Calcule le poids cotangent pour l'arête reliant vi à vj. */
|
||
double cotangent_weight(MyMesh &mesh, HalfedgeHandle vi_vj) {
|
||
Point pj = mesh.point(mesh.to_vertex_handle(vi_vj));
|
||
HalfedgeHandle vj_vb = mesh.next_halfedge_handle(vi_vj);
|
||
Point pb = mesh.point(mesh.to_vertex_handle(vj_vb));
|
||
HalfedgeHandle vj_vi = mesh.opposite_halfedge_handle(vi_vj);
|
||
Point pi = mesh.point(mesh.to_vertex_handle(vj_vi));
|
||
HalfedgeHandle vi_va = mesh.next_halfedge_handle(vj_vi);
|
||
Point pa = mesh.point(mesh.to_vertex_handle(vi_va));
|
||
double a = angle_between(pi, pa, pj);
|
||
double b = angle_between(pj, pb, pi);
|
||
return cotan(a) + cotan(b);
|
||
}
|
||
|
||
|
||
/* Calcule l'aire de chaque face et la stoque dans une propriété
|
||
* "face_area" du maillage. */
|
||
void compute_face_areas(MyMesh &mesh) {
|
||
auto area = OpenMesh::getOrMakeProperty<FaceHandle, double>
|
||
(mesh, "face_area");
|
||
for (FaceHandle fh : mesh.faces()) {
|
||
MyMesh::FaceVertexIter fv_it = mesh.fv_iter(fh);
|
||
Point pi = mesh.point(*fv_it);
|
||
Point pj = mesh.point(*++fv_it);
|
||
Point pk = mesh.point(*++fv_it);
|
||
area[fh] = triangle_area(pi, pj, pk);
|
||
}
|
||
}
|
||
|
||
|
||
/* Calcule l'aire barycentrique incidente à vert. */
|
||
double barycentric_vertex_area(MyMesh &mesh, VertexHandle vert) {
|
||
auto area = OpenMesh::getProperty<FaceHandle, double>(mesh, "face_area");
|
||
double sum = 0;
|
||
for (FaceHandle fh : mesh.vf_range(vert)) {
|
||
sum += area[fh];
|
||
}
|
||
return sum / 3;
|
||
}
|
||
|
||
|
||
double cotangent_mass(MyMesh &mesh, VertexHandle vi) {
|
||
double area = barycentric_vertex_area(mesh, vi);
|
||
return 1. / (2 * area);
|
||
}
|
||
|
||
|
||
Eigen::SparseMatrix<qreal> laplacian_matrix(MyMesh &mesh, double (*edge_weight)(MyMesh &, HalfedgeHandle), double (*vertex_mass)(MyMesh &, VertexHandle)) {
|
||
compute_face_areas(mesh);
|
||
size_t n_verts = mesh.n_vertices();
|
||
Eigen::SparseMatrix<double> weight(n_verts, n_verts);
|
||
Eigen::SparseMatrix<double> mass(n_verts, n_verts);
|
||
for (VertexHandle vi : mesh.vertices()) {
|
||
if (mesh.is_boundary(vi)) continue;
|
||
double sum = 0;
|
||
for (HalfedgeHandle vi_vj : mesh.voh_range(vi)) {
|
||
VertexHandle vj = mesh.to_vertex_handle(vi_vj);
|
||
double halfedge_weight = edge_weight(mesh, vi_vj);
|
||
weight.insert(vi.idx(), vj.idx()) = halfedge_weight;
|
||
sum -= halfedge_weight;
|
||
}
|
||
weight.insert(vi.idx(), vi.idx()) = sum;
|
||
mass.insert(vi.idx(), vi.idx()) = vertex_mass(mesh, vi);
|
||
}
|
||
return mass * weight;
|
||
}
|
||
|
||
|
||
void smooth(MyMesh &mesh, SmoothingMethod method, double cotan_factor) {
|
||
double factor;
|
||
Eigen::SparseMatrix<qreal> laplacian;
|
||
if (method == SmoothingMethod::COTANGENT) {
|
||
factor = cotan_factor;
|
||
laplacian = laplacian_matrix(mesh, cotangent_weight, cotangent_mass);
|
||
} else {
|
||
factor = 1;
|
||
laplacian = laplacian_matrix(mesh, uniform_weight, uniform_mass);
|
||
}
|
||
// laplacian = laplacian * laplacian;
|
||
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 * factor);
|
||
}
|
||
}
|