CGAL 5.1 - 3D Generalized Barycentric Coordinates
|
Barycentric coordinates are an important tool in computer graphics and geometric modeling. Originally, these coordinates were used to represent a given point with respect to a simplex but have been later generalized to more complex shapes.
The package 3D Generalized Barycentric Coordinates offers an efficient and robust implementation of three-dimensional closed-form, generalized barycentric coordinates defined for convex simplicial polytopes.
In particular, this package includes an implementation of Wachspress, discrete harmonic, mean value, and one extra function to calculate barycentric coordinates with respect to tetrahedra. In this implementation, we restrict our polyhedra to convex ones with triangular faces, although some of the coordinates may accept more general shapes.
Wachspress, discrete harmonic, and mean value coordinates are all generalized barycentric coordinates that can be computed analytically. All of the three analytic coordinates can be computed either by instantiating a class or through a free function. Tetrahedron coordinates can be computed only through the free function.
Similarly to Barycentric_coordinates::Computation_policy_2
in the 2D package, we can specify a computation policy Barycentric_coordinates::Computation_policy_3 that can be FAST or FAST_WITH_EDGE_CASES for each of the three analytical coordinates. The difference between them is that FAST_WITH_EDGE_CASES treats points near the boundaries by projecting them into the face of the polyhedron and then calculating triangle coordinates. Note that, different from the 2D package, there is not yet an implementation for a PRECISE algorithm.
The output of a query point is the coordinate value with respect to each vertex, the number of coordinates values being the same as the number of vertices, and the ordering is also the same.
All class and function templates are parameterized by a traits class, which is a model of the concept BarycentricTraits_3. It provides all necessary geometric primitives, predicates, and constructions, which are required for the computation. All models of Kernel can be used. A polyhedron is provided as a model of the concept FaceListGraph with a property map that maps a vertex from the polyhedron to BarycentricTraits_3::Point_3.
To facilitate the process of learning this package, we provide various examples with basic usage of different barycentric components.
In this example, we use the global function tetrahedron_coordinates_in_array()
. We compute coordinates for the tetrahedron whose vertices are the points in the set {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}. We use points inside, outside and at the boundary of the tetrahedron.
File Barycentric_coordinates_3/tetrahedron_coordinates.cpp
In this example, we generate 250 random points inside a unit sphere, centered at the origin, then we take the convex hull of this set of points and use this as our polyhedron. Finally, we calculate Wachspress coordinates for all of these 250 points.
File Barycentric_coordinates_3/wachspress_coordinates.cpp
This example is very similar to the one used with tetrahedron coordinates. We start with a regular icosahedron and for points inside, outside and at the boundary, we calculate discrete harmonic coordinates. In this example, we use the fast with edge cases algorithm because it treats points very close to the boundaries.
File Barycentric_coordinates_3/discrete_harmonic_coordinates.cpp
This example shows how to compute mean value coordinates for a set of points in a star-shaped polyhedron. We note that this type of coordinate is well-defined for a concave polyhedron but it may yield negative coordinate values for points outside the polyhedron's kernel (shown in blue).
File Barycentric_coordinates_3/mean_value_coordinates.cpp
This is an example that shows how to deform a simple smooth sphere into another shape, topologically equivalent. To achieve this, we enclose the sphere inside a cube cage, then we calculate the barycentric coordinate of each vertex with respect to the cube cage. After deforming the cage, we calculate the linear combination of the points to get the deformed sphere vertices.
File Barycentric_coordinates_3/shape_deformation.cpp
The precision of each coordinate depends on the used Kernel
. If an inexact kernel is used and the user is not sure if the points are near the boundaries, FAST_WITH_EDGE_CASES
algorithm should be used. Implementation details are described in [1] for Wachspress and mean value coordinates, and in [3] for discrete harmonic coordinates.
For each coordinate, it is necessary to make divisions by the signed distance between the query point and each face. So, if one of these distances is zero or close to zero (query point at the boundary), it will cause a division by zero error or numerical instability, respectively. To avoid this, the main purpose of the FAST_WITH_EDGE_CASES algorithm is to extend the region where the analytical coordinates are well-defined. It adds the guarantee to calculate points near the boundaries. The way it works is very simple: before calculating any coordinate, the algorithm checks, for each face, if the distance between the query point and the plane is less than one predetermined tolerance. If so, instead of calculating the analytical form of the coordinates, it decomposes the query point with respect to this particular face and then calculates triangle coordinates. However, for Wachspress coordinates, the 2D version is used because the faces are not necessarily triangular. Note that for every vertex that does not belong to this face, the coordinate value will be zero.
To compute tetrahedron coordinates of the query point q
, we adopt the simple formula:
where \(V_i\) is the signed volume of the sub-tetrahedron opposite to the vertex \(i\) and \(V\) is the total volume of the tetrahedron, that is \(V = V_0 + V_1 + V_2 + V_3\).
These coordinates can be computed exactly if an exact number type is chosen, for any query point in the space and with respect to any non-degenerate tetrahedron. No special cases are handled. The computation always yields the correct result. The notion of correctness depends on the precision of the used number type. Note that for exterior points some coordinate values will be negative.
For each vertex \(v\), let \(f_1, f_2, ..., f_k\) be the \(k\) faces incident to \(v\). We are assuming that the faces are taken in counterclockwise order.
We can define \({p_f}({q}) = \frac{{n_f}}{h_f({q})}\), and \(h_f({q}) = ({v} - {q})\cdot {n_f}\) as the perpendicular distance of \(q\) to \(f\). For the face \(f\), let \({n_f}\) denote its unit outward normal.
So, to compute Wachspress coordinates of the query point q
, we adopt the simple formula:
In this implementation, Wachspress coordinates are well defined in the closure of any convex polyhedra. If an exact number type is chosen, they are computed exactly.
To compute discrete harmonic coordinates of the query point q
, we adopt the simple formula:
where, within a triangle face T = {v1, v2, v3}, \(\theta_i^T\) is the dihedral angle between T and triangle {x, \(v_{i+1}\), \(v_{iā1}\)}, and \(h_i^T\) is the edge length \(|v_{i+1} ā v_{iā1}|\).
Discrete harmonic coordinates cannot be computed exactly due to a square root operation. Although, if an exact number type is used, the default precision of the computation depends only on two CGAL functions: CGAL::to_double() and CGAL::sqrt(). In this implementation, discrete harmonic coordinates are well defined in the closure of any convex polyhedra with triangular faces. Unlike Wachspress coordinates, they are not necessarily positive.
To compute mean value coordinates of the query point q
, we adopt the simple formula:
where a vertex v is projected to the point (unit vector) \(e_v = \frac{({v} - {q})}{|{v - q}|}\), and \({m_i} = \frac{{e_i} \times {e_{i+1}}}{|{e_i} \times {e_{i+1}}|}\). \(\beta_j \in (0, \pi)\) is the angle between \(e_i\) and \(e_{i+1}\).
Like discrete harmonic, mean value coordinates cannot be computed exactly due to a square root operation. In this implementation, mean value coordinates are well defined everywhere in the space, but just for polyhedra with triangular faces. Also, they are non-negative in the kernel of a star-shaped polyhedron.
Efficiency is really important in this implementation. These coordinates are used in applications that require calculations for millions of points, thus developing metrics to evaluate performance is absolutely necessary. In this section, we present benchmark results for each algorithm.
To make the benchmark and evaluate runtimes, for each analytic coordinate, we regularly sample approximately \(n^3\) ( \(n\) varying from 1 to 100) strictly interior points with respect to a cube with unit length, and then calculate its coordinate values (see figure below). The results are represented in a log-log scale plot and are the mean value of 10 loop iterations (see plot below).
The performance strongly depends on the chosen kernel, for this test, we choose to use Simple_cartesian<double>
because is much faster than others. Also, we can see that time (WP) << time (DH) < time (MV). This happens because wp implementation has fewer instructions per loop than the other two, so naturally, the computation time tends to be faster.
Tetrahedron coordinates are not shown in the same plot because the test is slightly different. For this one, we simply show in the table below the results for some pre-defined quantity of points. The test is done by regularly sampling strictly interior points with respect to a tetrahedron with unit sides that lies on the coordinate axis.
Number of points | Total time (in seconds) |
---|---|
1000 | 0.002662682533 |
50000 | 0.09459688663 |
100000 | 0.1865084648 |
500000 | 0.8238497972 |
1000000 | 1.665773582 |
5000000 | 8.359417105 |
To benchmark each coordinate, we used a 2.6 GHz Intel Core i7 processor (6 cores) and 16 GB DDR4 2933MHz memory. The installed operating system was Ubuntu 20.04 LTS.
This package was introduced during GSoC 2021 and implemented by Antonio Gomes under the supervision of Dmitry Anisimov.