User Manual

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>
#include <CGAL/Weights.h>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
using FT      = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
using Point_3 = typename Kernel::Point_3;
 
int main() {
 
  // 2D configuration.
  const Point_2 t2(-1,  0); // p0
  const Point_2 r2( 0, -1); // p1
  const Point_2 p2( 1,  0); // p2
  const Point_2 q2( 0,  0); // query
 
  // 3D configuration.
  const Point_3 t3(-1,  0, 1); // p0
  const Point_3 r3( 0, -1, 1); // p1
  const Point_3 p3( 1,  0, 1); // p2
  const Point_3 q3( 0,  0, 1); // query
 
  std::cout << "2D/3D tangent weight: ";
  std::cout << CGAL::Weights::tangent_weight(t2, r2, p2, q2) << "/";
  std::cout << CGAL::Weights::tangent_weight(t3, r3, p3, q3) << std::endl;
 
  std::cout << "2D/3D Shepard weight: ";
  std::cout << CGAL::Weights::shepard_weight(r2, q2, 2.0) << "/";
  std::cout << CGAL::Weights::shepard_weight(r3, q3, 2.0) << std::endl;
 
  std::cout << "2D/3D barycentric area: ";
  std::cout << CGAL::Weights::barycentric_area(p2, q2, r2) << "/";
  std::cout << CGAL::Weights::barycentric_area(p3, q3, r3) << std::endl;
  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>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
using FT      = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
 
int main() {
 
  // Create a polygon and a query point.
  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);
 
  // Allocate memory for weights and coordinates.
  std::vector<FT> weights;
  weights.reserve(polygon.size());
  std::vector<FT> coordinates;
  coordinates.reserve(polygon.size());
 
  // Compute barycentric weights.
  CGAL::Weights::discrete_harmonic_weights_2(polygon, query, std::back_inserter(weights));
  assert(weights.size() == polygon.size());
 
  std::cout << "2D weights: ";
  for (const FT weight : weights) {
    std::cout << weight << " ";
  }
  std::cout << std::endl;
 
  // Normalize weights in order to get barycentric coordinates.
  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>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
using FT      = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
 
using PointRange = std::vector<Point_2>;
using Creator    = CGAL::Creator_uniform_2<FT, Point_2>;
using Generator  = CGAL::Random_points_in_square_2<Point_2, Creator>;
using Wachspress = CGAL::Weights::Wachspress_weights_2<PointRange, Kernel>;
 
int main() {
 
  // Choose how many query points we want to generate.
  const std::size_t num_queries = 10;
 
  // Create a polygon.
  const std::vector<Point_2> polygon =
    { Point_2(0, 0), Point_2(1, 0), Point_2(1, 1), Point_2(0, 1) };
 
  // Generate a set of query points.
  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);
 
  // Allocate memory for weights and coordinates.
  std::vector<FT> weights;
  weights.reserve(num_queries * polygon.size());
  std::vector<FT> coordinates;
  coordinates.reserve(num_queries * polygon.size());
 
  // Compute barycentric weights.
  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;
  }
 
  // Normalize weights in order to get barycentric coordinates.
  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>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
using Point_2 = typename Kernel::Point_2;
using Point_3 = typename Kernel::Point_3;
 
// Custom traits class that has only two objects.
struct Custom_traits {
  using FT = typename Kernel::FT;
  using Point_2 = typename Kernel::Point_2;
  using Point_3 = typename Kernel::Point_3;
 
  decltype(auto) compute_squared_distance_2_object() const {
    return Kernel::Compute_squared_distance_2();
  }
  decltype(auto) compute_squared_distance_3_object() const {
    return Kernel::Compute_squared_distance_3();
  }
};
 
int main() {
 
  // 2D configuration.
  const Point_2 p2(0, 0);
  const Point_2 q2(0, 1);
 
  // 3D configuration.
  const Point_3 p3(0, 0, 1);
  const Point_3 q3(0, 1, 1);
 
  // Create custom traits.
  const Custom_traits ctraits;
 
  // Compute weights.
  std::cout << "2D/3D weight: ";
  std::cout << CGAL::Weights::inverse_distance_weight(p2, q2, ctraits) << "/";
  std::cout << CGAL::Weights::inverse_distance_weight(p3, q3, ctraits) << std::endl;
  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>
#include <CGAL/Weights.h>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
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 Mesh = CGAL::Surface_mesh<Point_3>;
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);
    A_i += CGAL::Weights::mixed_voronoi_area(p, q, r);
  }
  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;
}

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/three_point_family_weights.h>
 
// Typedefs.
using Kernel  = CGAL::Simple_cartesian<double>;
using FT      = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
 
int main() {
  std::cout << std::fixed;
 
  // 3D configuration.
  const Point_2 p0(0, 1);
  const Point_2 p1(2, 0);
  const Point_2 p2(7, 1);
  const Point_2 q0(3, 1);
 
  // Choose a type of the weight:
  // e.g. 0 - Wachspress (WP) weight; 1 - mean value (MV);
  const FT wp = 0.0;
  const FT mv = 1.0;
 
  // Compute WP and MV weights.
  std::cout << "3D Wachspress (WP, q0): ";
  std::cout << CGAL::Weights::three_point_family_weight(p0, p1, p2, q0, wp) << std::endl;
  std::cout << "3D mean value (MV, q0): ";
  std::cout << CGAL::Weights::three_point_family_weight(p0, p1, p2, q0, mv) << std::endl;
 
  // Converge WP towards MV.
  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: ";
    std::cout << CGAL::Weights::three_point_family_weight(p0, p1, p2, q0, x) << std::endl;
  }
  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.