- Author
- Dmitry Anisimov
Introduction
Many geometric algorithms rely on the intermediate computation of scalars, so-called weights, which are then used for solving different linear systems or to favor one result over another, also known as weighting. This package provides a simple and unified interface to different types of weights.
A typical example of a geometric algorithm that requires weights is the Laplace smoothing of a triangle mesh:
\(v_i \leftarrow v_i + h \lambda\Delta v_i\),
where \(v_i\) is the position of the mesh vertex \(i\), \(h\) is a sufficiently small time step, \(\lambda\) is the scalar diffusion coefficient, and \(\Delta v_i\) is the discrete average of the Laplace-Beltrami operator at vertex \(v_i\) computed using the cotangent weights:
\(\Delta v_i = w_i\sum_{v_j \in N_1(v_i)} w_{ij} (v_j - v_i)\),
where \(w_i = \frac{1}{2A_i}\) and \(w_{ij} = \cot\beta_{ij} + \cot\gamma_{ij}\) and \(A_i\) is a local averaging domain.
Here, the weights \(w_{ij}\) can be computed using the cotangent weight and the local averaging domain can be computed using the mixed Voronoi region weight. The algorithm above smooths the mesh geometry, resulting in a higher quality version of the original mesh. The full example of the discretized Laplacian for all vertices of a triangle mesh can be found here.
There are many other scenarios where the weights from this package are used. In particular, the following CGAL packages make use of weights described in this package: 2D Generalized Barycentric Coordinates, Polygon Mesh Processing, Triangulated Surface Mesh Deformation, Triangulated Surface Mesh Parameterization, Triangulated Surface Mesh Skeletonization, and The Heat Method.
Weights
We call analytic weights all weights which can be computed with a simple analytic expression. All weights from this package can be computed analytically. However, for better navigation through all available weights and their applications, we distinguish three typical groups of weights:
- Analytic Weights include all basic weights which can be computed for a query point with respect to its local neighbors in 2D or 3D, however these neighbors are defined. Usually, the configuration is a query point and three other points. These weights return one unique value per query point.
- Barycentric Weights include all weights which can be computed for a query point with respect to the vertices of a planar polygon. These weights return \(n\) values per query point, where \(n\) is the number of polygon vertices. Barycentric weights are also used for computing 2D barycentric coordinates.
- Weighting Regions include all weights which are used to balance other weights but are rarely used on their own. Sometimes, such weights are also referred to as local averaging regions. These weights are usually lengths, areas, and volumes of 2D and 3D objects.
Implementation
All weight functions have a simple and unified interface. In particular, all analytic weight functions usually take a query point and three other points in 2D or 3D and return a unique scalar. They all have the same signature and are parameterized by a traits class that must be a model of AnalyticWeightTraits_2 for 2D computations or AnalyticWeightTraits_3 for 3D computations.
The barycentric weight functions are parameterized by a traits class of the concept AnalyticWeightTraits_2 and they are all models of the concept BarycentricWeights_2. They take an input polygon and a query point and compute the weights at this point with respect to all vertices of the polygon. The computed weights are then returned in a container providing the corresponding output iterator. These weight functions also provide a property map mechanism for mapping a user type of the polygon vertex to CGAL::Point_2.
All weighting regions have the same signature and are parameterized by a traits class of the concept AnalyticWeightTraits_2 or AnalyticWeightTraits_3. The returned weight is a unique scalar.
The traits parameter can be omitted for all functions and classes if it can be deduced from the input point type using CGAL::Kernel_traits.
Several weights in this package have different implementations. One reason to have it is explained in section about Coplanarity. Another reason is that the same weights are named and computed differently in different communities. If one searches for these weights, one needs to know all their alternative names which is problematic. We provide the most common names and implementations of these weights.
Coplanarity
When computing weights for a query point \(q\) with respect to its neighbors \(p_0\), \(p_1\), and \(p_2\), the local configuration is a quadrilateral [ \(p_0\), \(p_1\), \(p_2\), \(q\)] or two connected triangles [ \(q\), \(p_0\), \(p_1\)] and [ \(q\), \(p_1\), \(p_2\)]. When working in 3D, these triangles are not necessarily coplanar, in other words, they do not belong to the same common plane.
Certain weights in this package support only coplanar configurations, while other weights support both. The weights which support non-coplanar configurations, provide the corresponding overloads with 3D points, while other weights accept only 2D point types. For example, cotangent weights support both coplanar and non-coplanar configurations, while discrete harmonic weights support only coplanar configurations.
Edge Cases
None of the weights in this package are defined for query points which belong to
- end segments so-called edges, for example [ \(p_0\), \(p_1\)] or [ \(p_1\), \(p_2\)] or any polygon edge and to
- end points so-called corners, for example \(p_0\), \(p_1\), or \(p_2\) or any polygon corner.
For example, if \(q\) = \(p_0\) or \(q \in [p_0, p_1]\), a weight may be undefined due to invalid operations such as division by zero.
Several weights also require additional conditions to be satisfied. Consult the reference manual for more details.
Examples
In this section, you can find a few examples of how and when the provided weights can be used.
The First Example
This trivial example shows how to compute several analytic weights and weighting regions. Other weights have the same interface so they can be computed analogously.
File Weights/weights.cpp
#include <CGAL/Simple_cartesian.h>
int main() {
const Point_2 t2(-1, 0);
const Point_2 r2( 0, -1);
const Point_2 p2( 1, 0);
const Point_2 q2( 0, 0);
const Point_3 t3(-1, 0, 1);
const Point_3 r3( 0, -1, 1);
const Point_3 p3( 1, 0, 1);
const Point_3 q3( 0, 0, 1);
std::cout << "2D/3D tangent weight: ";
std::cout << "2D/3D Shepard weight: ";
std::cout << "2D/3D barycentric area: ";
return EXIT_SUCCESS;
}
Computing 2D Coordinates for One Query Point
This example shows how to compute barycentric weights and barycentric coordinates, which are normalized barycentric weights, for a query point with respect to a polygon in 2D. Since we have only one query point, we use a free function to show the simplified interface. For multiple query points though, calling a free function is not efficient (see the following example for more details). The used type of barycentric weights is Discrete_harmonic_weights_2.
File Weights/coordinates_one_query.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Weights/discrete_harmonic_weights.h>
int main() {
const std::vector<Point_2> polygon =
{ Point_2(0, 0), Point_2(1, 0), Point_2(1, 1), Point_2(0, 1) };
const Point_2 query(0.5, 0.5);
std::vector<FT> weights;
weights.reserve(polygon.size());
std::vector<FT> coordinates;
coordinates.reserve(polygon.size());
assert(weights.size() == polygon.size());
std::cout << "2D weights: ";
for (const FT weight : weights) {
std::cout << weight << " ";
}
std::cout << std::endl;
FT sum = 0.0;
for (const FT weight : weights) {
sum += weight;
}
assert(sum != 0.0);
for (const FT weight : weights) {
coordinates.push_back(weight / sum);
}
assert(coordinates.size() == weights.size());
std::cout << "2D coordinates: ";
for (const FT coordinate : coordinates) {
std::cout << coordinate << " ";
}
std::cout << std::endl;
return EXIT_SUCCESS;
}
Computing 2D Coordinates for Multiple Query Points
This example shows how to compute barycentric weights and barycentric coordinates, which are normalized barycentric weights, for a set of query points with respect to a polygon in 2D. Since we have multiple query points, we first create a class and then use it to compute the weights. Using a class for multiple query points is preferred, because in that case, the memory required for computing weights is allocated only once, while when using a free function as in the previous example, it is allocated for each query point. The used type of barycentric weights is Wachspress_weights_2.
File Weights/coordinates_multiple_queries.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Weights/wachspress_weights.h>
using PointRange = std::vector<Point_2>;
int main() {
const std::size_t num_queries = 10;
const std::vector<Point_2> polygon =
{ Point_2(0, 0), Point_2(1, 0), Point_2(1, 1), Point_2(0, 1) };
std::vector<Point_2> queries;
queries.reserve(num_queries);
Generator generator(1.0);
std::copy_n(generator, num_queries, std::back_inserter(queries));
assert(queries.size() == num_queries);
std::vector<FT> weights;
weights.reserve(num_queries * polygon.size());
std::vector<FT> coordinates;
coordinates.reserve(num_queries * polygon.size());
Wachspress wachspress(polygon);
for (const Point_2& query : queries) {
wachspress(query, std::back_inserter(weights));
}
assert(weights.size() == num_queries * polygon.size());
std::cout << "2D weights: " << std::endl;
for (std::size_t i = 0; i < weights.size(); i += polygon.size()) {
for (std::size_t j = 0; j < polygon.size(); ++j) {
std::cout << weights[i + j] << " ";
}
std::cout << std::endl;
}
for (std::size_t i = 0; i < weights.size(); i += polygon.size()) {
FT sum = 0.0;
for (std::size_t j = 0; j < polygon.size(); ++j) {
sum += weights[i + j];
}
assert(sum != 0.0);
for (std::size_t j = 0; j < polygon.size(); ++j) {
coordinates.push_back(weights[i + j] / sum);
}
}
assert(coordinates.size() == weights.size());
std::cout << "2D coordinates: " << std::endl;
for (std::size_t i = 0; i < coordinates.size(); i += polygon.size()) {
for (std::size_t j = 0; j < polygon.size(); ++j) {
std::cout << coordinates[i + j] << " ";
}
std::cout << std::endl;
}
return EXIT_SUCCESS;
}
Weights with Custom Traits
As you could see from the reference manual, it is possible to provide your own traits class with basic geometric objects, constructions, and predicates to the weight functions. All weights in this CGAL component are models of the AnalyticWeightTraits_2 and AnalyticWeightTraits_3 concepts. However, many weights do not require all objects from these concepts. This example shows that the inverse distance weight, for instance, requires only the squared distance object which is specified in the custom traits class.
File Weights/custom_traits.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Weights/inverse_distance_weights.h>
struct Custom_traits {
decltype(auto) compute_squared_distance_2_object() const {
}
decltype(auto) compute_squared_distance_3_object() const {
}
};
int main() {
const Point_2 p2(0, 0);
const Point_2 q2(0, 1);
const Point_3 p3(0, 0, 1);
const Point_3 q3(0, 1, 1);
const Custom_traits ctraits;
std::cout << "2D/3D weight: ";
return EXIT_SUCCESS;
}
Constructing Weighted Laplacian
A typical example of using weights is discretizing Poisson and Laplace equations which play an important role in various geometry processing applications such as, for example, Surface Mesh Deformation and Surface Mesh Parameterization (see also Introduction for more details). This example shows how to write the discretized Laplacian for all vertices of the given triangle mesh in matrix notation. We use the standard cotangent discretization weighted by the areas of the mixed Voronoi cells around each mesh vertex.
File 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>
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 auto v0 = target(he, mesh);
const auto v1 = source(he, mesh);
const auto& q = get(pmap, v0);
const auto& p1 = get(pmap, v1);
const auto he_cw =
opposite(next(he, mesh), mesh);
auto v2 = source(he_cw, mesh);
const auto he_ccw = prev(
opposite(he, mesh), mesh);
v2 = source(he_ccw, mesh);
const auto& p2 = get(pmap, v2);
return CGAL::Weights::cotangent(p1, p2, q);
} else {
const auto& p0 = get(pmap, v2);
return CGAL::Weights::cotangent(q, p0, p1);
}
}
const auto he_cw =
opposite(next(he, mesh), mesh);
const auto v2 = source(he_cw, mesh);
const auto he_ccw = prev(
opposite(he, mesh), mesh);
const auto v3 = source(he_ccw, mesh);
const auto& p0 = get(pmap, v2);
const auto& p2 = get(pmap, v3);
}
template<typename PointMap>
FT get_w_i(const Mesh& mesh, const VD v_i, const PointMap pmap) {
FT A_i = 0.0;
const auto v0 = v_i;
const auto init = halfedge(v_i, mesh);
assert(v0 == target(he, mesh));
const auto v1 = source(he, mesh);
const auto 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);
const auto imap = get(CGAL::vertex_index, mesh);
std::map<std::size_t, FT> w_i;
for (const auto& v_i : vertices(mesh)) {
w_i[get(imap, v_i)] = get_w_i(mesh, v_i, pmap);
}
for (const auto& he : halfedges(mesh)) {
const auto vi = source(he, mesh);
const auto 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() {
Mesh mesh;
const auto v0 = mesh.add_vertex(Point_3(0, 2, 0));
const auto v1 = mesh.add_vertex(Point_3(2, 2, 0));
const auto v2 = mesh.add_vertex(Point_3(0, 0, 0));
const auto v3 = mesh.add_vertex(Point_3(2, 0, 0));
const auto 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);
const std::size_t n = 5;
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;
}
Convergence
This little example shows how to use the family of weights which includes multiple types of weights in one function. In particular, we show how, by changing an input parameter, we converge from the Wachspress_weights_2 to the Mean_value_weights_2.
File Weights/convergence.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Weights/authalic_weights.h>
#include <CGAL/Weights/three_point_family_weights.h>
int main() {
std::cout << std::fixed;
const Point_2 p0(0, 1);
const Point_2 p1(2, 0);
const Point_2 p2(7, 1);
const Point_2 q0(3, 1);
const FT wp = 0.0;
const FT mv = 1.0;
std::cout << "3D Wachspress (WP, q0): ";
std::cout << "3D mean value (MV, q0): ";
std::cout << "Converge WP to MV on q0: " << std::endl;
const FT step = 0.1;
for (FT x = 0.0; x <= 1.0; x += step) {
std::cout << "3D x: ";
}
return EXIT_SUCCESS;
}
History
This package is a part of the weights unification effort inside CGAL that has been carried out by Dmitry Anisimov in 2020.
Acknowledgments
We wish to thank Guillaume Damiand, Andreas Fabri, and Mael Rouxel-Labbé for useful discussions and reviews.
1.8.13