CGAL 5.1 - 3D Generalized Barycentric Coordinates
User Manual

Authors
Antonio Gomes, Dmitry Anisimov

Introduction

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.

Software Design

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.

Examples

To facilitate the process of learning this package, we provide various examples with basic usage of different barycentric components.

Tetrahedron Coordinates

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

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_3/tetrahedron_coordinates.h>
//Typedefs
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
int main(){
// Construct tetrahedron
const Point_3 p0(0.0, 0.0, 0.0);
const Point_3 p1(1.0, 0.0, 0.0);
const Point_3 p2(0.0, 1.0, 0.0);
const Point_3 p3(0.0, 0.0, 1.0);
// Instantiate some interior, boundary, and exterior query points for which we compute coordinates.
const std::vector<Point_3> queries = {
Point_3(0.25f , 0.25f, 0.25f), Point_3(0.3f, 0.2f, 0.3f), // interior query points
Point_3(0.1f, 0.1f, 0.1f), Point_3(0.2f, 0.5f, 0.3f), // interior query points
Point_3(0.0f , 0.0f, 0.5f), Point_3(0.4f, 0.4f, 0.0f), // boundary query points
Point_3(0.0f, 0.4f, 0.4f), Point_3(0.4f, 0.0f, 0.4f), // boundary query points
Point_3(0.5f, 0.5f, 0.5f), Point_3(2.0f, 0.0f, 0.0f), // exterior query points
Point_3(-1.0f, -1.0f, 1.0f), Point_3(0.5f, 0.5f, -2.0f)}; // exterior query point
std::cout << std::endl << "tetrahedra coordinates (all queries): " << std::endl
<< std::endl;
// Get an array of triangle coordinates for all query points
for(std::size_t i = 0; i < queries.size(); i++){
const auto coords_array =
std::cout << "tetrahedra coordinates (query " << i << "): " <<
coords_array[0] << " " << coords_array[1] << " " <<
coords_array[2] << " " << coords_array[3] << std::endl;
}
return EXIT_SUCCESS;
}

Wachspress Coordinates

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

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Wachspress_coordinates_3.h>
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
namespace PMP = CGAL::Polygon_mesh_processing;
int main()
{
std::vector<Point_3> points;
const std::size_t number_of_points = 250;
std::copy_n(gen, number_of_points, std::back_inserter(points));
Surface_mesh sm;
CGAL::convex_hull_3(points.begin(), points.end(), sm);
const std::size_t number_of_vertices = num_vertices(sm);
WP wp(sm, CP3::FAST_WITH_EDGE_CASES);
std::cout << "Computed Wachspress coordinates: " << std::endl << std::endl;
for(std::size_t i = 0; i < number_of_points; i++){
std::vector<FT> coordinates;
coordinates.reserve(number_of_vertices);
wp(points[i], std::back_inserter(coordinates));
std::cout << "Point " << i + 1 << ": " << std::endl;
for(std::size_t j = 0; j < number_of_vertices; j++)
std::cout << "Coordinate " << j + 1 << " = " << coordinates[j] << "; " << std::endl;
std::cout << std::endl;
}
return EXIT_SUCCESS;
}

Discrete Harmonic Coordinates

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

#define PHI 1.6180339887498948482
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Discrete_harmonic_coordinates_3.h>
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(){
Surface_mesh icosahedron;
CGAL::make_icosahedron(icosahedron, Point_3(0.0, 0.0, 0.0), 2.0);
PMP::triangulate_faces(faces(icosahedron), icosahedron);
std::vector<FT> coords;
std::vector<Point_3> queries{
Point_3(-1, 1 + PHI, PHI), Point_3(0.5, (1+3*PHI)/2, PHI/2), Point_3(1, 1+PHI, -PHI), //Boundary
Point_3(-1, 1, 1), Point_3(0, 0, 1), Point_3(0, 2, 1), //Interior
Point_3(0, 2*PHI, 4), Point_3(0, 3, 2*PHI), Point_3(4, 0, 0)}; //EXterior
std::cout << std::endl << "Discrete harmonic coordinates : " << std::endl << std::endl;
for (const auto& query : queries){
coords.clear();
icosahedron, query, std::back_inserter(coords), CP3::FAST_WITH_EDGE_CASES);
// Output discrete harmonics coordinates.
for (std::size_t i = 0; i < coords.size() -1; ++i) {
std::cout << coords[i] << ", ";
}
std::cout << coords[coords.size() -1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Mean Value Coordinates

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).

fig__mv_3_example Example's point pattern.



File Barycentric_coordinates_3/mean_value_coordinates.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Mean_value_coordinates_3.h>
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(){
Surface_mesh concave;
const Point_3 p0(0, 3, 0);
const Point_3 p1(1, 1, 0);
const Point_3 p2(3, 0, 0);
const Point_3 p3(0, 0, 0);
const Point_3 p4(0, 0, 3);
const Point_3 p5(0, 3, 3);
const Point_3 p6(1, 1, 3);
const Point_3 p7(3, 0, 3);
CGAL::make_hexahedron(p0, p1, p2, p3, p4, p5, p6, p7, concave);
PMP::triangulate_faces(faces(concave), concave);
std::vector<FT> coords;
std::vector<Point_3> queries{
Point_3(FT(1)/FT(2), FT(1)/FT(2), FT(1)), Point_3(FT(1)/FT(3), FT(1)/FT(3), FT(2)), // Only points in the kernel
Point_3(FT(4)/FT(3), FT(1)/FT(3), FT(1)), Point_3(FT(4)/FT(3), FT(1)/FT(3), FT(2)),
Point_3(FT(1)/FT(3), FT(4)/FT(3), FT(1)), Point_3(FT(1)/FT(3), FT(4)/FT(3), FT(2))};
std::cout << std::endl << "Mean value coordinates : " << std::endl << std::endl;
for (const auto& query : queries){
coords.clear();
concave, query, std::back_inserter(coords), CP3::FAST);
// Output mean value coordinates.
for (std::size_t i = 0; i < coords.size() -1; ++i) {
std::cout << coords[i] << ", ";
}
std::cout << coords[coords.size() -1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Shape deformation

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.

fig__shape_deform_example_3 The shape on the left is deformed into the shape on the right.



File Barycentric_coordinates_3/shape_deformation.cpp

#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Mean_value_coordinates_3.h>
#include <fstream>
// default triangulation for Surface_mesher
using Tr = CGAL::Surface_mesh_default_triangulation_3;
using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
using Sphere_3 = Kernel::Sphere_3;
using Point_3 = Kernel::Point_3;
using FT = Kernel::FT;
typedef FT (*Function)(Point_3);
using Surface_3 = CGAL::Implicit_surface_3<Kernel, Function>;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
namespace PMP = CGAL::Polygon_mesh_processing;
FT sphere_function (Point_3 p){
const FT x2=p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
return x2+y2+z2-1;
}
int main() {
Tr tr;
C2t3 c2t3(tr);
Surface_3 surface(sphere_function, Sphere_3(CGAL::ORIGIN, 2.));
CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30., 0.1, 0.1);
CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());
Surface_mesh sm;
Surface_mesh deformed;
CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, sm);
deformed = sm;
Surface_mesh quad_cage;
const Point_3 p0(2, -2, -2), p0_new(5, -5, -5);
const Point_3 p1(2, 2, -2), p1_new(3, 3, -3);
const Point_3 p2(-2, 2, -2), p2_new(-2, 2, -2);
const Point_3 p3(-2, -2, -2), p3_new(-3, -3, -3);
const Point_3 p4(-2, -2, 2), p4_new(-3, -3, 3);
const Point_3 p5(2, -2, 2), p5_new(4, -4, 4);
const Point_3 p6(2, 2, 2), p6_new(2, 2, 3);
const Point_3 p7(-2, 2, 2), p7_new(-3, 3, 3);
CGAL::make_hexahedron(p0, p1, p2, p3, p4, p5, p6, p7, quad_cage);
PMP::triangulate_faces(faces(quad_cage), quad_cage);
auto vertex_to_point_map = get_property_map(CGAL::vertex_point, deformed);
std::vector<FT> coords;
std::vector<Point_3> target_cube{p0_new, p1_new, p2_new, p3_new,
p4_new, p5_new, p6_new, p7_new};
for(auto& v : vertices(deformed)){
const Point_3 vertex_val = get(vertex_to_point_map, v);
coords.clear();
mv(vertex_val, std::back_inserter(coords));
FT x = FT(0), y = FT(0), z = FT(0);
for(std::size_t i = 0; i < 8; i++){
x += target_cube[i].x() * coords[i];
y += target_cube[i].y() * coords[i];
z += target_cube[i].z() * coords[i];
}
put(vertex_to_point_map, v, Point_3(x, y, z));
}
std::ofstream out_original("sphere.off");
out_original << sm << std::endl;
std::ofstream out_deformed("deformed_sphere.off");
out_deformed << deformed << std::endl;
return EXIT_SUCCESS;
}

Edge Cases

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.

Tetrahedron Coordinates

To compute tetrahedron coordinates of the query point q, we adopt the simple formula:

\(w_i = \frac{V_i}{V}\)

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.

Wachspress Coordinates

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:

\(w_v({q}) = \sum_{i=2}^{k-1}det({p_{f_1}}({q}), {p_{f_i}}({q}), {p_{f_{i+1}}}({q}))\)

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.

Discrete Harmonic Coordinates

To compute discrete harmonic coordinates of the query point q, we adopt the simple formula:

\(w_i = \sum_{T : v_i \in T} \frac{cot[\theta_i^T]h_i^T}{2}\)

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.

Mean Value Coordinates

To compute mean value coordinates of the query point q, we adopt the simple formula:

\(w_i = \frac{1}{2}\sum_{j=1}^3 \beta_j\frac{{m_j}\cdot {m_{i+1}}}{{e_i}\cdot {m_{i+1}}}\)

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.

Performance

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).

fig__bc_coords_bench_structure_3 The points shown in red are the sample points used to make the benchmark.


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.

fig__bc_coords_bench_3 Time in seconds to compute \(n^3\) coordinate values for a cube. Are represented in the graph: Wachspress (blue), discrete harmonic (orange), and mean value (green).


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.

History

This package was introduced during GSoC 2021 and implemented by Antonio Gomes under the supervision of Dmitry Anisimov.

Acknowledgments

CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Random_points_in_sphere_3
CGAL::Barycentric_coordinates::mean_value_coordinates_3
OutIterator mean_value_coordinates_3(const TriangleMesh &triangle_mesh, const Point_3 &query, OutIterator c_begin, const Computation_policy_3 policy=Computation_policy_3::FAST_WITH_EDGE_CASES)
computes 3D mean value coordinates.
Definition: Mean_value_coordinates_3.h:409
CGAL::Barycentric_coordinates::Computation_policy_3
Computation_policy_3
Computation_policy_3 provides a way to choose an asymptotic time complexity of the algorithm and its ...
Definition: barycentric_enum_3.h:36
copy_n
OutputIterator copy_n(InputIterator first, Size n, OutputIterator result)
CGAL::Barycentric_coordinates::Mean_value_coordinates_3
3D mean value coordinates.
Definition: Mean_value_coordinates_3.h:54
Kernel::Sphere_3
unspecified_type Sphere_3
Kernel::Point_3
unspecified_type Point_3
CGAL::Barycentric_coordinates::discrete_harmonic_coordinates_3
OutIterator discrete_harmonic_coordinates_3(const TriangleMesh &triangle_mesh, const Point_3 &query, OutIterator c_begin, const Computation_policy_3 policy=Computation_policy_3::FAST)
computes 3D discrete harmonic coordinates.
Definition: Discrete_harmonic_coordinates_3.h:400
ORIGIN
const CGAL::Origin ORIGIN
Kernel
CGAL::Barycentric_coordinates::Wachspress_coordinates_3
3D Wachspress coordinates.
Definition: Wachspress_coordinates_3.h:55
make_icosahedron
boost::graph_traits< Graph >::halfedge_descriptor make_icosahedron(Graph &g, const P &center=P(0, 0, 0), typename CGAL::Kernel_traits< P >::Kernel::FT radius=1.0)
CGAL::Barycentric_coordinates::tetrahedron_coordinates_in_array
std::array< typename GeomTraits::FT, 4 > tetrahedron_coordinates_in_array(const typename GeomTraits::Point_3 &p0, const typename GeomTraits::Point_3 &p1, const typename GeomTraits::Point_3 &p2, const typename GeomTraits::Point_3 &p3, const typename GeomTraits::Point_3 &query, const GeomTraits &traits)
computes tetrahedron coordinates.
Definition: tetrahedron_coordinates.h:145
make_hexahedron
boost::graph_traits< Graph >::halfedge_descriptor make_hexahedron(const P &p0, const P &p1, const P &p2, const P &p3, const P &p4, const P &p5, const P &p6, const P &p7, Graph &g)
Kernel::FT
unspecified_type FT
CGAL::Simple_cartesian