CGAL 5.6 - Weights
Weights/weighted_laplacian.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Eigen_solver_traits.h>
#include <CGAL/Eigen_matrix.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Weights.h>
// Typedefs.
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Solver_traits = CGAL::Eigen_solver_traits<
Eigen::SparseLU<CGAL::Eigen_sparse_matrix<FT>::EigenType> >;
using Matrix = Solver_traits::Matrix;
using VD = boost::graph_traits<Mesh>::vertex_descriptor;
using HD = boost::graph_traits<Mesh>::halfedge_descriptor;
template<typename PointMap>
FT get_w_ij(const Mesh& mesh, const HD he, const PointMap pmap) {
const VD v0 = target(he, mesh);
const VD v1 = source(he, mesh);
const auto& q = get(pmap, v0); // query
const auto& p1 = get(pmap, v1); // neighbor j
if (is_border_edge(he, mesh)) {
const HD he_cw = opposite(next(he, mesh), mesh);
VD v2 = source(he_cw, mesh);
if (is_border_edge(he_cw, mesh)) {
const HD he_ccw = prev(opposite(he, mesh), mesh);
v2 = source(he_ccw, mesh);
const auto& p2 = get(pmap, v2); // neighbor jp
return CGAL::Weights::cotangent(p1, p2, q);
} else {
const auto& p0 = get(pmap, v2); // neighbor jm
return CGAL::Weights::cotangent(q, p0, p1);
}
}
const HD he_cw = opposite(next(he, mesh), mesh);
const VD v2 = source(he_cw, mesh);
const HD he_ccw = prev(opposite(he, mesh), mesh);
const VD v3 = source(he_ccw, mesh);
const auto& p0 = get(pmap, v2); // neighbor jm
const auto& p2 = get(pmap, v3); // neighbor jp
return CGAL::Weights::cotangent_weight(p0, p1, p2, q) / 2.0;
}
template<typename PointMap>
FT get_w_i(const Mesh& mesh, const VD v_i, const PointMap pmap) {
FT A_i = 0.0;
const VD v0 = v_i;
const HD init = halfedge(v_i, mesh);
for (const HD he : halfedges_around_target(init, mesh)) {
assert(v0 == target(he, mesh));
if (is_border(he, mesh)) { continue; }
const VD v1 = source(he, mesh);
const VD v2 = target(next(he, mesh), mesh);
const auto& p = get(pmap, v0);
const auto& q = get(pmap, v1);
const auto& r = get(pmap, v2);
}
assert(A_i != 0.0);
return 1.0 / (2.0 * A_i);
}
void set_laplacian_matrix(const Mesh& mesh, Matrix& L) {
const auto pmap = get(CGAL::vertex_point, mesh); // vertex to point map
const auto imap = get(CGAL::vertex_index, mesh); // vertex to index map
// Precompute Voronoi areas.
std::map<std::size_t, FT> w_i;
for (const VD v_i : vertices(mesh)) {
w_i[get(imap, v_i)] = get_w_i(mesh, v_i, pmap);
}
// Fill the matrix.
for (const HD he : halfedges(mesh)) {
const VD vi = source(he, mesh);
const VD vj = target(he, mesh);
const std::size_t i = get(imap, vi);
const std::size_t j = get(imap, vj);
if (i > j) { continue; }
const FT w_ij = w_i[j] * get_w_ij(mesh, he, pmap);
L.set_coef(i, j, w_ij, true);
L.set_coef(j, i, w_ij, true);
L.add_coef(i, i, -w_ij);
L.add_coef(j, j, -w_ij);
}
}
int main() {
// Create mesh.
Mesh mesh;
const VD v0 = mesh.add_vertex(Point_3(0, 2, 0));
const VD v1 = mesh.add_vertex(Point_3(2, 2, 0));
const VD v2 = mesh.add_vertex(Point_3(0, 0, 0));
const VD v3 = mesh.add_vertex(Point_3(2, 0, 0));
const VD v4 = mesh.add_vertex(Point_3(1, 1, 1));
mesh.add_face(v0, v2, v4);
mesh.add_face(v2, v3, v4);
mesh.add_face(v3, v1, v4);
mesh.add_face(v1, v0, v4);
mesh.add_face(v2, v3, v1);
mesh.add_face(v1, v0, v2);
assert(CGAL::is_triangle_mesh(mesh));
// Set discretized Laplacian.
const std::size_t n = 5; // we have 5 vertices
Matrix L(n, n);
set_laplacian_matrix(mesh, L);
std::cout << std::fixed; std::cout << std::showpos;
for (std::size_t i = 0; i < n; ++i) {
for (std::size_t j = 0; j < n; ++j) {
std::cout << L.get_coef(i, j) << " ";
}
std::cout << std::endl;
}
return EXIT_SUCCESS;
}