diff --git a/src/smoothing.cpp b/src/smoothing.cpp index ec03dec..4d25f9e 100644 --- a/src/smoothing.cpp +++ b/src/smoothing.cpp @@ -42,8 +42,9 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) { area += triangle_area(p, pi, pb); count++; } - area /= 3.; - return sum / (2.*count); + // area /= 3.; + // return sum / (2.*area); + return sum / count; } @@ -73,46 +74,48 @@ Eigen::SparseMatrix laplacian_matrix(const MyMesh &mesh) { area += triangle_area(p, pi, pb); count++; } - area /= 3.; + // area /= 3.; cotangent.insert(vert.idx(), vert.idx()) = -sum; - mass.insert(vert.idx(), vert.idx()) = 1. / area; + mass.insert(vert.idx(), vert.idx()) = 1. / (4. * area); + // mass.insert(vert.idx(), vert.idx()) = 1. / count; } 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]); - // } - // } + 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 - 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); - } + // // Approche matricielle + // Eigen::SparseMatrix laplacian = laplacian_matrix(mesh); + // // laplacian = laplacian * laplacian; + // 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); + // } }