CGAL 6.2 - 3D Generalized Barycentric Coordinates

Loading...
Searching...
No Matches
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 non-convex polyhedra. The calculated coordinates are not normalized.

# Software Design

Wachspress, discrete harmonic, and mean value coordinates are all generalized barycentric coordinates that can be computed analytically. Each of the three analytic coordinates can be computed either by a free function or a function object (Functor). 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, which 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 onto the face of the polyhedron and then calculating triangle coordinates. Note that, unlike the 2D package, no implementation of a PRECISE algorithm is yet available.

The output of each method are the coordinates for a point with respect to each vertex of the mesh. The number of coordinates 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()`. We compute coordinates for the tetrahedron defined by the points (0,0,0), (1,0,0), (0,1,0), and (0,0,1). We use points inside, outside and at the boundary of the tetrahedron.

Example: Barycentric_coordinates_3/tetrahedron_coordinates_3.cpp

Show / Hide
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_3/tetrahedron_coordinates_3.h>
using Kernel = CGAL::Simple_cartesian<double>;
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.25 , 0.25, 0.25), Point_3(0.3, 0.2, 0.3), // interior query points
Point_3(0.1, 0.1, 0.1), Point_3(0.2, 0.5, 0.3), // interior query points
Point_3(0.0 , 0.0, 0.5), Point_3(0.4, 0.4, 0.0), // boundary query points
Point_3(0.0, 0.4, 0.4), Point_3(0.4, 0.0, 0.4), // boundary query points
Point_3(0.5, 0.5, 0.5), Point_3(2.0, 0.0, 0.0), // exterior query points
Point_3(-1.0, -1.0, 1.0), Point_3(0.5, 0.5, -2.0)}; // exterior query point
std::cout << std::endl << "tetrahedron 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 =
CGAL::Barycentric_coordinates::tetrahedron_coordinates_3(p0, p1, p2, p3, queries[i]);
std::cout << "tetrahedron coordinates (query " << i << "): " <<
coords_array[0] << " " << coords_array[1] << " " <<
coords_array[2] << " " << coords_array[3] << std::endl;
}
return EXIT_SUCCESS;
}
Kernel::Point_3
CGAL::Barycentric_coordinates::tetrahedron_coordinates_3
OutputIterator tetrahedron_coordinates_3(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, OutputIterator oi, const GeomTraits &traits)
computes barycentric coordinates with respect to a tetrahedron.
Definition: tetrahedron_coordinates_3.h:73
Kernel
CGAL::Simple_cartesian

## 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 each of these 250 points.

Example: Barycentric_coordinates_3/wachspress_coordinates_3.cpp

Show / Hide
#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/Barycentric_coordinates_3/Wachspress_coordinates_3.h>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Computation_policy_3 = CGAL::Barycentric_coordinates::Computation_policy_3;
using Wachspress = CGAL::Barycentric_coordinates::Wachspress_coordinates_3<Surface_mesh>;
int main() {
CGAL::Random_points_in_sphere_3<Point_3> gen(1.0);
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);
Wachspress wp(sm, Computation_policy_3::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<double> 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;
}
CGAL::Barycentric_coordinates::Wachspress_coordinates_3
computes 3D Wachspress coordinates with respect to a closed convex triangle mesh.
Definition: Wachspress_coordinates_3.h:56
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Random_points_in_sphere_3
CGAL::convex_hull_3
void convex_hull_3(InputIterator first, InputIterator last, PolygonMesh &pm, const Traits &ch_traits=Default_traits)
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
CGAL::Barycentric_coordinates::Computation_policy_3::FAST_WITH_EDGE_CASES
@ FAST_WITH_EDGE_CASES
Computation has a linear time complexity with respect to the number of face vertices,...

## 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 on the boundary, we calculate discrete harmonic coordinates. In this example, we use the fast with edge cases algorithm because it properly handles points very close to the boundaries.

Example: Barycentric_coordinates_3/discrete_harmonic_coordinates_3.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/generators.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Barycentric_coordinates_3/Discrete_harmonic_coordinates_3.h>
#define PHI 1.6180339887498948482
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Computation_policy_3 = CGAL::Barycentric_coordinates::Computation_policy_3;
int main() {
Surface_mesh icosahedron;
CGAL::make_icosahedron(icosahedron, Point_3(0.0, 0.0, 0.0), 2.0);
CGAL::Polygon_mesh_processing::triangulate_faces(faces(icosahedron), icosahedron);
std::vector<double> 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();
CGAL::Barycentric_coordinates::discrete_harmonic_coordinates_3(
icosahedron, query, std::back_inserter(coords));
// 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;
}
CGAL::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)
CGAL::Barycentric_coordinates::discrete_harmonic_coordinates_3
OutputIterator discrete_harmonic_coordinates_3(const TriangleMesh &tmesh, const Point &query, OutputIterator oi, const NamedParameters &np=parameters::default_values())
computes 3D discrete harmonic coordinates with respect to a closed convex triangle mesh.
Definition: Discrete_harmonic_coordinates_3.h:438

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

Figure 113.1 The star-shaped polyhedron used in mean_value_coordinates_3.cpp and the six query points, for which the mean value barycentric coordinates are computed.

Example: Barycentric_coordinates_3/mean_value_coordinates_3.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/generators.h>
#include <CGAL/Barycentric_coordinates_3/Mean_value_coordinates_3.h>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Computation_policy_3 = CGAL::Barycentric_coordinates::Computation_policy_3;
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,
CGAL::parameters::do_not_triangulate_faces(false));
std::vector<double> coords;
std::vector<Point_3> queries {
Point_3(0.5, 0.5, 1), Point_3(0.33, 0.33, 2), // Only points in the kernel
Point_3(1.33, 0.33, 1), Point_3(1.33, 0.33, 2),
Point_3(0.33, 1.33, 1), Point_3(0.33, 1.33, 2)};
std::cout << std::endl << "Mean value coordinates : " << std::endl << std::endl;
for (const auto& query : queries) {
coords.clear();
CGAL::Barycentric_coordinates::mean_value_coordinates_3(
concave, query, std::back_inserter(coords), CGAL::parameters::computation_policy(Computation_policy_3::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;
}
CGAL::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, const NamedParameters &np=parameters::default_values())
CGAL::Barycentric_coordinates::mean_value_coordinates_3
OutputIterator mean_value_coordinates_3(const TriangleMesh &tmesh, const Point &query, OutputIterator oi, const NamedParameters &np=parameters::default_values())
computes 3D mean value barycentric coordinates with respect to a closed triangle mesh.
Definition: Mean_value_coordinates_3.h:446
CGAL::Barycentric_coordinates::Computation_policy_3::FAST
@ FAST
Computation has a linear time complexity with respect to the number of face vertices,...

## Shape Deformation

This example shows how to deform a simple smooth sphere into another topologically equivalent shape. To achieve this, we create a box that encloses the sphere; then we calculate the barycentric coordinates of each vertex with respect to the box. After deforming the box, we use the barycentric coordinates to relocate the vertices of the sphere.

Figure 113.2 The shape on the left is deformed into the shape on the right.

Example: Barycentric_coordinates_3/shape_deformation_3.cpp

Show / Hide
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/generators.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <CGAL/Barycentric_coordinates_3.h>
#include <fstream>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv) {
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");
Surface_mesh sm;
if(!CGAL::IO::read_polygon_mesh(filename, sm)) {
std::cerr << "Invalid input: " << filename << std::endl;
return 1;
}
Surface_mesh deformed;
deformed = sm;
Surface_mesh quad_cage;
const Point_3 p0( 200, -200, -200), p0_new( 500, -500, -500);
const Point_3 p1( 200, 200, -200), p1_new( 300, 300, -300);
const Point_3 p2(-200, 200, -200), p2_new(-200, 200, -200);
const Point_3 p3(-200, -200, -200), p3_new(-300, -300, -300);
const Point_3 p4(-200, -200, 200), p4_new(-300, -300, 300);
const Point_3 p5( 200, -200, 200), p5_new( 400, -400, 400);
const Point_3 p6( 200, 200, 200), p6_new( 200, 200, 300);
const Point_3 p7(-200, 200, 200), p7_new(-300, 300, 300);
CGAL::make_hexahedron(p0, p1, p2, p3, p4, p5, p6, p7, quad_cage,
CGAL::parameters::do_not_triangulate_faces(false));
CGAL::Barycentric_coordinates::Mean_value_coordinates_3<Surface_mesh> mv(quad_cage);
auto vertex_point_map = get_property_map(CGAL::vertex_point, deformed);
std::vector<double> coords;
std::vector<Point_3> target_cube{p0_new, p1_new, p2_new, p3_new,
p4_new, p5_new, p6_new, p7_new};
for(Surface_mesh::Vertex_index v : vertices(deformed)) {
const Point_3 vertex_val = get(vertex_point_map, v);
coords.clear();
mv(vertex_val, std::back_inserter(coords));
Point_3 p = CGAL::Barycentric_coordinates::apply_barycentric_coordinates_3(target_cube, coords);
put(vertex_point_map, v, p);
}
std::ofstream out_original("armadillo.off");
out_original << sm << std::endl;
std::ofstream out_deformed("deformed_armadillo_mv.off");
out_deformed << deformed << std::endl;
return EXIT_SUCCESS;
}
Barycentric_coordinates_3.h
A convenience header that includes all free functions and classes for 3D barycentric coordinates in c...
CGAL::Barycentric_coordinates::Mean_value_coordinates_3
computes 3D mean value barycentric coordinates with respect to a closed triangle mesh.
Definition: Mean_value_coordinates_3.h:56
CGAL::IO::read_polygon_mesh
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
CGAL::Barycentric_coordinates::apply_barycentric_coordinates_3
boost::property_traits< VertexPointMap >::value_type apply_barycentric_coordinates_3(const TriangleMesh &tmesh, const CoordinateRange &coordinates, VertexPointMap vpm, GeomTraits geom_traits=GeomTraits())
computes a point location from barycentric coordinates with respect to a triangle mesh.
Definition: Barycentric_coordinates_3.h:76

# Edge Cases

The precision of each coordinate depends on the used `Kernel`. If an inexact kernel is used and the user is not sure whether the points are located 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 [2] 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 on the boundary), it will cause a division by zero error or numerical instability, respectively. To this end we introduce the FAST_WITH_EDGE_CASES algorithm. Its main purpose 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 simple: before calculating any coordinate, the algorithm checks, for each face, if the distance between the query point and the supporting plane of the face is less than a 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

We adopt the simple formula below to compute tetrahedron coordinates of the query point `q`:

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

where \(V_i\) is the signed volume of the sub-tetrahedron opposite to the vertex \(i\), i.e., the tetrahedron where the vertex \(i\) is replaced by the query point `q`. \(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 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_0, f_1, ..., f_{k-1}\) 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.

The following formula is adopted to compute Wachspress coordinates of the query point `q`:

\(w_v({q}) = \sum_{i=1}^{k-2}det({p_{f_0}}({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 in an exact manner.

## Discrete Harmonic Coordinates

We adopt the formula below to compute discrete harmonic coordinates of the query point `q`:

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

where, within a triangle face \(T = \{v0, v1, v2\}\), \(\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

The following formula is adopted to compute mean value coordinates of the query point `q`:

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

where a vertex v is projected onto 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 exactly computed 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 crucial in this implementation. These coordinates are used in applications that require calculations for millions of points; thus, developing metrics to evaluate performance is necessary. In this section, we present benchmark results for each algorithm.

The benchmark and runtimes are evaluated by regularly sampling \(n^3\) ( \(n\) varying from 1 to 100) points from the interior of the unit cube, and calculating their coordinate values (see figure below).

The results are averaged over 10 executions and represented in a log-log scale plot.

Figure 113.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 the Wachspress implementation has fewer instructions per loop than the other two. For 100k points those generalized barycentric coordinates are calculated within a convex hull with 500 vertices the timings are 17.5s (WP), 25.4s (DH) and 33.3s (MV).

Figure 113.4 Time in seconds to compute \(n^3\) coordinate values for a cube. 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 lie on the coordinate axis.

Number of points | Total time (in seconds)
--- | ---
50000 | 0.09
100000 | 0.19
500000 | 0.82
1000000 | 1.67
5000000 | 8.36

To benchmark each coordinate, we used a 2.5 GHz Intel Core i7 processor (8 cores) and 32 GB DDR4 2933MHz memory. The installed operating system was Windows 11.

# History

This package was introduced during GSoC 2021 and implemented by Antonio Gomes under the supervision of Dmitry Anisimov. Sven Oesau finished the final work on API and documentation.