CGAL 5.5.3 - 2D Generalized Barycentric Coordinates
User Manual

Authors
Dmitry Anisimov, David Bommes, Kai Hormann, and Pierre Alliez

Introduction

Barycentric coordinates are widely used in computer graphics and computational mechanics to determine a position of a point in the plane with respect to a triangle. These coordinates have been later generalized to support simple polygons in 2D and polyhedra in 3D.

This package offers an efficient and robust implementation of 2D generalized barycentric coordinates defined for simple polygons in the plane. If coordinates with respect to multivariate scattered points instead of a polygon are required, please refer to natural neighbor coordinates from the package 2D and Surface Function Interpolation.

In particular, this package includes an implementation of Wachspress, discrete harmonic, mean value, and harmonic coordinates, and provides some extra functions to compute barycentric coordinates with respect to segments and triangles.

overview.svg
Figure 100.1 Wachspress (WP), discrete harmonic (DH), mean value (MV), and harmonic (HM) coordinate functions for a convex polygon plotted with respect to the marked vertex.

Software Design

Mean value and harmonic coordinates are the most generic coordinates in this package, because they allow an arbitrary simple polygon as input. Wachspress and discrete harmonic coordinates are, by definition, limited to strictly convex polygons. Segment coordinates take as input any non-degenerate segment, and triangle coordinates allow an arbitrary non-degenerate triangle.

Wachspress, discrete harmonic, mean value, and harmonic coordinates are all generalized barycentric coordinates. However, while Wachspress, discrete harmonic, and mean value coordinates can be computed analytically, harmonic coordinates cannot. They first need to be approximated over a triangulation of the interior part of the polygon. Once approximated, they can be evaluated analytically at any point inside the polygon.

For all analytic coordinates, we provide two algorithms. One has a linear time complexity, but may suffer imprecisions near the polygon boundary, while the second one is precise but has a quadratic time complexity. The user can choose the preferred algorithm by specifying a computation policy Barycentric_coordinates::Computation_policy_2.

All analytic barycentric coordinates for polygons can be computed either by instantiating a class or through one of the free functions. Harmonic coordinates can be computed only by instantiating a class that must be parameterized by a model of the concept DiscretizedDomain_2. Segment and triangle coordinates can be computed only through the free functions. For more information see the Reference Manual.

Any point in the plane may be taken as a query point. However, we do not recommend using Wachspress and discrete harmonic coordinates with query points outside the closure of a polygon, because they are not well-defined for some of these points. The same holds for harmonic coordinates, which are not defined everywhere outside the polygon. For more information see Section Edge Cases.

The output of the computation is a set of coordinate values at the given query point with respect to the polygon vertices. That means that the number of returned coordinates per query point equates the number of polygon vertices. The ordering of the coordinates is the same as the ordering of polygon vertices.

All class and function templates are parameterized by a traits class, which is a model of the concept BarycentricTraits_2. It provides all necessary geometric primitives, predicates, and constructions, which are required for the computation. All models of Kernel can be used. A polygon is provided as a range of vertices with a property map that maps a vertex from the polygon to CGAL::Point_2.

If you do not know which coordinate function best fits your application, you can address the table below for some advise.

Coordinates Properties Valid domain Closed form Queries Speed
Segment All 2D non-degenerate segments Yes Everywhere on the supporting line +++
Triangle All 2D non-degenerate triangles Yes Everywhere in 2D +++
Discrete harmonic May be negative Strongly convex polygons Yes Everywhere inside the polygon ++
Wachspress All Strongly convex polygons Yes Everywhere inside the polygon ++
Mean value May be negative Simple polygons Yes Everywhere in 2D ++
Harmonic All Simple polygons No Everywhere inside the polygon +
Note
This is the second version of the package with the modified and improved API. The package still supports the old API. See more details here.

Examples

In order to facilitate the process of learning this package, we provide various examples with a basic usage of different barycentric components.

Segment Coordinates

This example illustrates the use of the global function segment_coordinates_2(). We compute coordinates at three green points along the segment \([v_0, v_1]\) and at two blue points outside this segment but along its supporting line. The symmetry of the query points helps recognizing errors that may have occurred during construction of the example. The used Kernel is exact.

seg_coord_example.svg
Figure 100.2 Example's point pattern.


File Barycentric_coordinates_2/segment_coordinates.cpp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/segment_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
int main() {
const FT y = FT(2) / FT(5);
// Construct a segment.
const Point_2 source(FT(0), y);
const Point_2 target(FT(2), y);
// Construct three interior and two exterior query points.
const std::vector<Point_2> queries = {
Point_2(FT(2) / FT(5), y), // interior query points
Point_2(FT(5) / FT(5), y),
Point_2(FT(8) / FT(5), y),
Point_2(-FT(1) / FT(5), y), // exterior query points
Point_2(FT(11) / FT(5), y) };
// Compute segment coordinates.
std::vector<FT> coordinates;
coordinates.reserve(queries.size() * 2);
for (const auto& query : queries) {
source, target, query, std::back_inserter(coordinates));
}
// Output all segment coordinates.
std::cout << std::endl << "segment coordinates (all queries): " << std::endl << std::endl;
for (std::size_t i = 0; i < coordinates.size(); i += 2) {
std::cout <<
coordinates[i + 0] << ", " <<
coordinates[i + 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Triangle Coordinates

In this example, we show how to use the global function triangle_coordinates_2(). We compute coordinates for three sets of points: interior (green), boundary (red), and exterior (blue). Note that some of the coordinate values for the exterior points are negative. The used Kernel is inexact.

tri_coord_example.svg
Figure 100.3 Example's point pattern.


File Barycentric_coordinates_2/triangle_coordinates.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_2/triangle_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
int main() {
// Construct a triangle.
const Point_2 p0(0.0, 0.0);
const Point_2 p1(2.0, 0.5);
const Point_2 p2(1.0, 2.0);
// Construct several interior, boundary, and exterior query points.
const std::vector<Point_2> queries = {
Point_2(0.50, 0.50), // interior query points
Point_2(1.00, 0.50), Point_2(1.0, 0.75), Point_2(1.00, 1.0),
Point_2(1.00, 1.25), Point_2(1.0, 1.50), Point_2(0.75, 1.0),
Point_2(1.25, 1.00), Point_2(1.5, 0.75),
Point_2(2.0, 0.50), Point_2(1.0, 2.00), // boundary query points
Point_2(1.0, 0.25), Point_2(1.5, 1.25), Point_2(0.5, 1.0),
Point_2(0.25, 1.00), Point_2(0.50, 1.75), // exterior query points
Point_2(1.50, 1.75), Point_2(1.75, 1.50) };
// Compute triangle coordinates.
std::vector<FT> coordinates;
coordinates.reserve(queries.size() * 3);
for (const auto& query : queries) {
p0, p1, p2, query, std::back_inserter(coordinates));
}
// Output all triangle coordinates.
std::cout << std::endl << "triangle coordinates (all queries): " << std::endl << std::endl;
for (std::size_t i = 0; i < coordinates.size(); i += 3) {
std::cout <<
coordinates[i + 0] << ", " <<
coordinates[i + 1] << ", " <<
coordinates[i + 2] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Wachspress Coordinates

In the following example, we generate 100 random points (green/red/black), then we take the convex hull (red/black) of this set of points as our polygon (black), and compute Wachspress coordinates at all the generated points. The used Kernel is inexact.

wp_coord_example.svg
Figure 100.4 Example's point pattern.


File Barycentric_coordinates_2/wachspress_coordinates.cpp

#include <CGAL/convex_hull_2.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Barycentric_coordinates_2/Wachspress_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
int main() {
// Choose how many query points we want to generate.
const std::size_t num_queries = 100;
// Create vectors to store query points and polygon vertices.
std::vector<Point_2> queries, convex;
// Generate a set of query points.
queries.reserve(num_queries);
Generator generator(1.0);
std::copy_n(generator, num_queries, std::back_inserter(queries));
// Find the convex hull of the generated query points.
// This convex hull gives the vertices of a convex polygon
// that contains all the generated points.
queries.begin(), queries.end(), std::back_inserter(convex));
// Compute Wachspress coordinates for all query points.
std::cout << std::endl << "Wachspress coordinates (interior + boundary): " << std::endl << std::endl;
std::vector<FT> coordinates;
coordinates.reserve(convex.size());
for (const auto& query : queries) {
coordinates.clear();
convex, query, std::back_inserter(coordinates));
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Discrete Harmonic Coordinates

In this example, we compute discrete harmonic coordinates for a set of green (interior), red (boundary), and blue (exterior) points with respect to a unit square. We also demonstrate the use of various containers, both random access and serial access, different property maps, and the ability to choose a computation policy. For points on the polygon boundary, we use the free function boundary_coordinates_2(). The used Kernel is exact.

dh_coord_example.svg
Figure 100.5 Example's point pattern.


File Barycentric_coordinates_2/discrete_harmonic_coordinates.cpp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/boundary_coordinates_2.h>
#include <CGAL/Barycentric_coordinates_2/Discrete_harmonic_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
struct Info {
Info(const std::string _name) :
name(_name) { }
std::string name;
};
using Vertex = std::pair<Point_2, Info>;
using Vertex_range = std::vector<Vertex>;
using Discrete_harmonic_coordinates_2 =
using Policy =
int main() {
Kernel kernel;
Point_map point_map;
// Construct a unit square.
const std::vector<Vertex> square = {
std::make_pair(Point_2(0, 0), Info("1")), std::make_pair(Point_2(1, 0), Info("2")),
std::make_pair(Point_2(1, 1), Info("3")), std::make_pair(Point_2(0, 1), Info("4"))
};
// Construct the class with discrete harmonic weights.
// We do not check for edge cases since we know the exact positions
// of all our points. We speed up the computation by using the O(n) algorithm.
const Policy policy = Policy::FAST;
Discrete_harmonic_coordinates_2 discrete_harmonic_2(square, policy, kernel, point_map);
// Construct the center point of the unit square.
const Point_2 center(FT(1) / FT(2), FT(1) / FT(2));
// Compute discrete harmonic weights for the center point.
std::list<FT> weights;
discrete_harmonic_2.weights(center, std::back_inserter(weights));
std::cout << std::endl << "discrete harmonic weights (center): ";
for (const FT& weight : weights) {
std::cout << weight << " ";
}
std::cout << std::endl;
// Compute discrete harmonic coordinates for the center point.
std::list<FT> coordinates;
discrete_harmonic_2(center, std::back_inserter(coordinates));
std::cout << std::endl << "discrete harmonic coordinates (center): ";
for (const FT& coordinate : coordinates) {
std::cout << coordinate << " ";
}
std::cout << std::endl;
// Construct several interior points.
const std::vector<Point_2> interior_points = {
Point_2(FT(1) / FT(5), FT(1) / FT(5)),
Point_2(FT(4) / FT(5), FT(1) / FT(5)),
Point_2(FT(4) / FT(5), FT(4) / FT(5)),
Point_2(FT(1) / FT(5), FT(4) / FT(5)) };
// Compute discrete harmonic weights for all interior points.
std::cout << std::endl << "discrete harmonic weights (interior): " << std::endl << std::endl;
std::vector<FT> ws;
for (const auto& query : interior_points) {
ws.clear();
discrete_harmonic_2.weights(query, std::back_inserter(ws));
for (std::size_t i = 0; i < ws.size() - 1; ++i) {
std::cout << ws[i] << ", ";
}
std::cout << ws[ws.size() - 1] << std::endl;
}
// Compute discrete harmonic coordinates for all interior point.
std::cout << std::endl << "discrete harmonic coordinates (interior): " << std::endl << std::endl;
std::vector<FT> bs;
for (const auto& query : interior_points) {
bs.clear();
discrete_harmonic_2(query, std::back_inserter(bs));
for (std::size_t i = 0; i < bs.size() - 1; ++i) {
std::cout << bs[i] << ", ";
}
std::cout << bs[bs.size() - 1] << std::endl;
}
// Construct 2 boundary points on the second and fourth edges.
const Point_2 e2(1, FT(4) / FT(5));
const Point_2 e4(0, FT(4) / FT(5));
// Compute discrete harmonic coordinates = boundary coordinates
// for these 2 points one by one.
coordinates.clear();
square, e2, std::back_inserter(coordinates), kernel, point_map);
square, e4, std::back_inserter(coordinates), kernel, point_map);
std::cout << std::endl << "boundary coordinates (edge 2 and edge 4): ";
for (const FT& coordinate : coordinates) {
std::cout << coordinate << " ";
}
std::cout << std::endl;
// Construct 6 other boundary points: 2 on the first and third edges respectively
// and 4 at the vertices.
const std::vector<Point_2> es13 = {
Point_2(FT(1) / FT(2), 0), // edges
Point_2(FT(1) / FT(2), 1),
// vertices
Point_2(0, 0), Point_2(1, 0),
Point_2(1, 1), Point_2(0, 1)
};
// Compute discrete harmonic coordinates = boundary coordinates for all 6 points.
std::cout << std::endl << "boundary coordinates (edge 1, edge 3, and vertices): " << std::endl << std::endl;
for (const auto& query : es13) {
bs.clear();
square, query, std::back_inserter(bs), point_map); // we can skip kernel here
for (std::size_t i = 0; i < bs.size() - 1; ++i) {
std::cout << bs[i] << ", ";
}
std::cout << bs[bs.size() - 1] << std::endl;
}
// Construct 2 points outside the unit square - one from the left and one from the right.
// Even if discrete harmonic coordinates may not be valid for some exterior points,
// we can still do it.
const Point_2 l(FT(-1) / FT(2), FT(1) / FT(2));
const Point_2 r(FT(3) / FT(2), FT(1) / FT(2));
// Compute discrete harmonic coordinates for all exterior points.
coordinates.clear();
discrete_harmonic_2(l, std::back_inserter(coordinates));
discrete_harmonic_2(r, std::back_inserter(coordinates));
std::cout << std::endl << "discrete harmonic coordinates (exterior): ";
for (const FT& coordinate : coordinates) {
std::cout << coordinate << " ";
}
std::cout << std::endl << std::endl;
return EXIT_SUCCESS;
}

Mean Value Coordinates

This is an example that illustrates how to compute mean value coordinates for a set of green points in a star-shaped polygon. We note that this type of coordinates is well-defined for such a concave polygon while Wachspress and discrete harmonic coordinates are not. However, it may yield negative coordinate values for points outside the polygon's kernel (shown in red). We speed up the computation using the linear time complexity algorithm by specifying a computation policy Barycentric_coordinates::Computation_policy_2. The used Kernel is inexact.

mv_coord_example.svg
Figure 100.6 Example's point pattern.


File Barycentric_coordinates_2/mean_value_coordinates.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Mean_value_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
int main() {
// Construct a star-shaped polygon.
const std::vector<Point_2> star_shaped = {
Point_2(0.0, 0.0), Point_2( 0.1, -0.8), Point_2(0.3, 0.0), Point_2(0.6, -0.5),
Point_2(0.6, 0.1), Point_2( 1.1, 0.6), Point_2(0.3, 0.2), Point_2(0.1, 0.8),
Point_2(0.1, 0.2), Point_2(-0.7, 0.0) };
// Construct some interior points in the polygon.
const std::vector<Point_2> interior_points = {
Point_2(0.12, -0.45), Point_2(0.55, -0.3), Point_2(0.9 , 0.45),
Point_2(0.15, 0.35), Point_2(-0.4, 0.04), Point_2(0.11, 0.11),
Point_2(0.28, 0.12), // the only point in the kernel of the star shaped polygon
Point_2(0.55, 0.11) };
// Choose a computation policy.
// We do not check for edge cases since we know
// that all our points are strictly interior.
const Policy policy = Policy::PRECISE;
// Create a vector `std::vector` to store coordinates.
std::vector<FT> coordinates;
coordinates.reserve(star_shaped.size());
// Compute mean value coordinates for all interior points.
std::cout << std::endl << "mean value coordinates (interior): " << std::endl << std::endl;
for (const auto& query : interior_points) {
coordinates.clear();
star_shaped, query, std::back_inserter(coordinates), policy);
// Output mean value coordinates.
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Harmonic Coordinates

This example illustrates how to discretize the interior part of the polygon and compute harmonic coordinates at the vertices of the discretized domain, which is represented by a 2D Delaunay triangulation. Once computed, harmonic coordinate functions can be evaluated at any point in the closure of the polygon. To illustrate such an evaluation, we compute the barycenter of each triangle and evaluate harmonic coordinates at this barycenter. Since harmonic coordinates can only be approximated, the used Kernel is inexact.


File Barycentric_coordinates_2/harmonic_coordinates.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_2/Delaunay_domain_2.h>
#include <CGAL/Barycentric_coordinates_2/Harmonic_coordinates_2.h>
// Typedefs.
using Point_2 = Kernel::Point_2;
using Point_range = std::vector<Point_2>;
using Domain =
using Harmonic_coordinates_2 =
int main() {
// Construct a simple polygon.
const std::vector<Point_2> polygon = {
Point_2(0.03, 0.05), Point_2(0.07, 0.04), Point_2(0.10, 0.04),
Point_2(0.14, 0.04), Point_2(0.17, 0.07), Point_2(0.20, 0.09),
Point_2(0.22, 0.11), Point_2(0.25, 0.11), Point_2(0.27, 0.10),
Point_2(0.30, 0.07), Point_2(0.31, 0.04), Point_2(0.34, 0.03),
Point_2(0.37, 0.02), Point_2(0.40, 0.03), Point_2(0.42, 0.04),
Point_2(0.44, 0.07), Point_2(0.45, 0.10), Point_2(0.46, 0.13),
Point_2(0.46, 0.19), Point_2(0.47, 0.26), Point_2(0.47, 0.31),
Point_2(0.47, 0.35), Point_2(0.45, 0.37), Point_2(0.41, 0.38),
Point_2(0.38, 0.37), Point_2(0.35, 0.36), Point_2(0.32, 0.35),
Point_2(0.30, 0.37), Point_2(0.28, 0.39), Point_2(0.25, 0.40),
Point_2(0.23, 0.39), Point_2(0.21, 0.37), Point_2(0.21, 0.34),
Point_2(0.23, 0.32), Point_2(0.24, 0.29), Point_2(0.27, 0.24),
Point_2(0.29, 0.21), Point_2(0.29, 0.18), Point_2(0.26, 0.16),
Point_2(0.24, 0.17), Point_2(0.23, 0.19), Point_2(0.24, 0.22),
Point_2(0.24, 0.25), Point_2(0.21, 0.26), Point_2(0.17, 0.26),
Point_2(0.12, 0.24), Point_2(0.07, 0.20), Point_2(0.03, 0.15),
Point_2(0.01, 0.10), Point_2(0.02, 0.07)
};
// Use seeds to mark the interior part of the polygon.
std::list<Point_2> seeds;
seeds.push_back(Point_2(0.1, 0.1));
// Construct a Delaunay domain.
const double max_edge_length = 0.01;
Domain domain(polygon);
domain.create(max_edge_length, seeds);
// Compute harmonic coordinates at the vertices of the domain.
Harmonic_coordinates_2 harmonic_coordinates_2(polygon, domain);
harmonic_coordinates_2.compute();
// Use it to store coordinates.
std::vector<double> coordinates;
coordinates.reserve(polygon.size());
// Output harmonic coordinates.
// We output only the first 20 results.
std::cout.precision(1);
std::cout << std::endl << "harmonic coordinates (computed): " << std::endl << std::endl;
for (std::size_t k = 0; k < 20; ++k) {
coordinates.clear();
harmonic_coordinates_2(k, std::back_inserter(coordinates));
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
// Evaluate harmonic coordinates at the barycenters of the domain triangles.
// We output only the first 20 results.
std::cout << std::endl << "harmonic coordinates (evaluated): " << std::endl << std::endl;
std::vector<Point_2> barycenters;
domain.barycenters(std::back_inserter(barycenters));
for (std::size_t k = 0; k < 20; ++k) {
coordinates.clear();
harmonic_coordinates_2(barycenters[k], std::back_inserter(coordinates));
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Terrain Modeling

This is an advanced example that illustrates how to use generalized barycentric coordinates for height interpolation with applications to terrain modeling. It also shows how to use a non-default traits class with our package instead of a Kernel traits class. Suppose we know the boundary of three-dimensional piece of terrain that can be represented as a polygon with several three-dimensional vertices, where the third dimension indicates the corresponding height. The task is to propagate the height from the known sample points on the boundary to the polygon's interior. This gives an approximate estimation of the terrain's surface in this region.

terrain.svg
Figure 100.7 A 2D polygon with 50 vertices representing a piece of terrain with convex and concave parts. The height is not shown.

In this example, we project a 3D polygon orthogonally onto the 2D plane using the class CGAL::Projection_traits_xy_3, triangulate its interior using the class Delaunay_domain_2, and compute mean value coordinates at the vertices of this triangulation with respect to the polygon vertices. Finally, we interpolate the height data from the polygon boundary to its interior using the computed coordinates and the global interpolation function from the package 2D and Surface Function Interpolation.


File Barycentric_coordinates_2/terrain_height_modeling.cpp

#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Delaunay_domain_2.h>
#include <CGAL/Barycentric_coordinates_2/Mean_value_coordinates_2.h>
// Typedefs.
using FT = typename Projection::FT;
using Point = typename Projection::Point_2;
using Point_range = std::vector<Point>;
using Domain =
using Mean_value_coordinates_2 =
using Vertex_function_value = std::map<Point, FT, typename Projection::Less_xy_2>;
using Function_value_access = CGAL::Data_access<Vertex_function_value>;
using Point_with_coordinate = std::pair<Point, FT>;
int main() {
// Construct a polygon that bounds a three-dimensional terrain.
// Note that the z-coordinate of each vertex represents the height function.
// Projection in 2D is performed automatically by the Projection traits class.
const std::vector<Point> polygon = {
Point(0.03, 0.05, 0.00), Point(0.07, 0.04, 0.02), Point(0.10, 0.04, 0.04),
Point(0.14, 0.04, 0.06), Point(0.17, 0.07, 0.08), Point(0.20, 0.09, 0.10),
Point(0.22, 0.11, 0.12), Point(0.25, 0.11, 0.14), Point(0.27, 0.10, 0.16),
Point(0.30, 0.07, 0.18), Point(0.31, 0.04, 0.20), Point(0.34, 0.03, 0.22),
Point(0.37, 0.02, 0.24), Point(0.40, 0.03, 0.26), Point(0.42, 0.04, 0.28),
Point(0.44, 0.07, 0.30), Point(0.45, 0.10, 0.32), Point(0.46, 0.13, 0.34),
Point(0.46, 0.19, 0.36), Point(0.47, 0.26, 0.38), Point(0.47, 0.31, 0.40),
Point(0.47, 0.35, 0.42), Point(0.45, 0.37, 0.44), Point(0.41, 0.38, 0.46),
Point(0.38, 0.37, 0.48), Point(0.35, 0.36, 0.50), Point(0.32, 0.35, 0.52),
Point(0.30, 0.37, 0.54), Point(0.28, 0.39, 0.56), Point(0.25, 0.40, 0.58),
Point(0.23, 0.39, 0.60), Point(0.21, 0.37, 0.62), Point(0.21, 0.34, 0.64),
Point(0.23, 0.32, 0.66), Point(0.24, 0.29, 0.68), Point(0.27, 0.24, 0.70),
Point(0.29, 0.21, 0.72), Point(0.29, 0.18, 0.74), Point(0.26, 0.16, 0.76),
Point(0.24, 0.17, 0.78), Point(0.23, 0.19, 0.80), Point(0.24, 0.22, 0.82),
Point(0.24, 0.25, 0.84), Point(0.21, 0.26, 0.86), Point(0.17, 0.26, 0.88),
Point(0.12, 0.24, 0.90), Point(0.07, 0.20, 0.92), Point(0.03, 0.15, 0.94),
Point(0.01, 0.10, 0.97), Point(0.02, 0.07, 1.00)
};
// Construct a Delaunay domain.
std::vector<Point> seeds;
seeds.push_back(Point(0.1, 0.1, 0.0));
Domain domain(polygon);
domain.create(0.05, seeds);
// Associate each polygon vertex with the corresponding function value.
Vertex_function_value vertex_function_value;
for (const auto& vertex : polygon) {
vertex_function_value.insert(
std::make_pair(vertex, vertex.z()));
}
// Construct the class with the mean value weights.
Mean_value_coordinates_2 mean_value_coordinates_2(polygon);
// Compute mean value coordinates and use them to interpolate data
// from the polygon boundary to its interior.
std::vector<FT> coordinates;
coordinates.reserve(polygon.size());
std::vector<Point_with_coordinate> boundary;
boundary.resize(polygon.size());
std::vector<Point> queries;
queries.reserve(domain.number_of_vertices());
for (std::size_t i = 0; i < domain.number_of_vertices(); ++i) {
const auto& query = domain.vertex(i);
coordinates.clear();
mean_value_coordinates_2(query, std::back_inserter(coordinates));
for (std::size_t i = 0; i < polygon.size(); ++i) {
boundary[i] = std::make_pair(polygon[i], coordinates[i]);
}
boundary.begin(), boundary.end(), FT(1),
Function_value_access(vertex_function_value));
queries.push_back(Point(query.x(), query.y(), f));
}
// Output interpolated heights.
std::cout << std::endl << "interpolated heights (all queries): " << std::endl << std::endl;
for (const auto& query : queries) {
std::cout << query.z() << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

As a result, we get a smooth function inside the polygon that approximates the underlying terrain surface.

terrain_interpolation.png
Figure 100.8 The interpolated data. The color bar represents the corresponding height.

Shape Deformation

This is another advanced example that shows how to use generalized barycentric coordinates in order to deform a given 2D shape into another shape as shown in the figure below. Harmonic coordinates satisfy all the properties of barycentric coordinates for complicated concave polygons and hence this is our choice to perform a shape deformation. Note that even though harmonic coordinates are guaranteed to be positive inside a polygon, they do not guarantee a bijective mapping between the source and target shapes that is the target mesh can fold-over the target polygon after the mapping (see the little fold-over in the left shoulder of the target shape).

shape_deformation.svg
Figure 100.9 The shape on the left is deformed into the shape on the right. The zoom shows a fold-over in the left shoulder of the target shape where the red triangle goes over the polygon boundary.


File Barycentric_coordinates_2/shape_deformation.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Delaunay_domain_2.h>
#include <CGAL/Barycentric_coordinates_2/Harmonic_coordinates_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
using Point_range = std::vector<Point_2>;
using Domain =
using Harmonic_coordinates_2 =
int main() {
// Construct the source and target shapes.
// The number of vertices in both shapes must be equal.
const std::vector<Point_2> source_shape = {
Point_2(1, 0), Point_2(2, 0), Point_2(3, 3), Point_2(4, 0), Point_2(5, 0),
Point_2(4, 3), Point_2(4, 5), Point_2(5, 4), Point_2(5, 5), Point_2(4, 6),
Point_2(2, 6), Point_2(1, 5), Point_2(1, 4), Point_2(2, 5), Point_2(2, 3)
};
const std::vector<Point_2> target_shape = {
Point_2(2, 0), Point_2(3, 0), Point_2(3, 3), Point_2(3, 0), Point_2(4, 0),
Point_2(4, 3), Point_2(4, 5), Point_2(5, 6), Point_2(5, 7), Point_2(4, 6),
Point_2(2, 6), Point_2(1, 7), Point_2(1, 6), Point_2(2, 5), Point_2(2, 3)
};
assert(target_shape.size() == source_shape.size());
// Use seeds to mark the interior part of the source shape.
const std::vector<Point_2> seeds = { Point_2(3, 5) };
// Construct a Delaunay domain.
const FT max_edge_length = FT(1) / FT(3);
Domain domain(source_shape);
domain.create(max_edge_length, seeds);
// Use it to store coordinates.
std::vector< std::vector<FT> > coordinates;
coordinates.reserve(domain.number_of_vertices());
// Compute harmonic coordinates at the vertices of the
// discretized interior domain of the source shape.
Harmonic_coordinates_2 harmonic_coordinates_2(source_shape, domain);
harmonic_coordinates_2.compute();
harmonic_coordinates_2(std::back_inserter(coordinates));
// Deform the source domain into the target domain.
// We output only the first 20 results.
std::cout << std::endl << "shape deformation: " << std::endl << std::endl;
for (std::size_t k = 0; k < 20; ++k) {
FT x = FT(0), y = FT(0);
for (std::size_t i = 0; i < coordinates[k].size(); ++i) {
x += coordinates[k][i] * target_shape[i].x();
y += coordinates[k][i] * target_shape[i].y();
}
std::cout << "deformed domain vertex: (" << x << ", " << y << ")" << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

But despite the possible fold-overs, a similar technique can be used for image warping in 2D and character articulation in 3D. For example in 2D, we first enclose an image, which we want to deform, into a simple polygon so-called cage, we then bound each image pixel to this cage using barycentric coordinates, and finally deform this cage into a new one, which also deforms the underlying image, as shown in the figure below for harmonic coordinates.

image_warping.png
Figure 100.10 An image on the left is deformed into a new image on the right using a 2D concave polygon (grey) and harmonic coordinates computed at each image pixel with respect to the vertices of this polygon.

Affine Coordinates

This is an example, where we show how a lambda expression can be used to define a model of generalized barycentric coordinates. To make this example useful, we implement affine generalized coordinates for a set of scattered points. The used Kernel is inexact.

aff_coord_example.svg
Figure 100.11 Example's point pattern.


File Barycentric_coordinates_2/affine_coordinates.cpp

#include <Eigen/Core>
#include <Eigen/Dense>
#include <CGAL/barycenter.h>
#include <CGAL/Simple_cartesian.h>
#include <boost/iterator/transform_iterator.hpp>
// Typedefs.
using Point_2 = Kernel::Point_2;
using VectorXd = Eigen::VectorXd;
using MatrixXd = Eigen::MatrixXd;
using Output_iterator =
std::back_insert_iterator< std::vector<double> >;
int main() {
// Create a set of vertices.
const std::vector<Point_2> vertices = {
Point_2(0.0, 0.0), Point_2(0.75, 0.25), Point_2(0.5, 0.5), Point_2(0.4, -0.2) };
// Create a set of query points.
const std::vector<Point_2> queries = {
Point_2(0.2, 0.2), Point_2(0.3, 0.3), Point_2(0.4, 0.4) };
// Create a lambda function with affine coordinates.
// This implementation is based on the following paper:
// S. Waldron. Affine generalized barycentric coordinates.
// Jaen Journal on Approximation, 3(2):209-226, 2011.
// This function is a model of the `AnalyticWeights_2` concept.
const auto affine = [&](
const Point_2& query,
Output_iterator coordinates) {
const std::size_t n = vertices.size();
const auto lambda = [](const Point_2& p){ return std::make_pair(p, 1.0); };
const Point_2 b = CGAL::barycenter(
boost::make_transform_iterator(vertices.begin(), lambda),
boost::make_transform_iterator(vertices.end() , lambda), Kernel());
MatrixXd V(2, n);
for (std::size_t i = 0; i < n; ++i) {
V(0, i) = vertices[i].x() - b.x();
V(1, i) = vertices[i].y() - b.y();
}
const auto A = V.adjoint();
const auto mat = V * A;
const auto inv = mat.inverse();
Point_2 diff; VectorXd vec(2);
for (std::size_t i = 0; i < n; ++i) {
const double x = query.x() - b.x();
const double y = query.y() - b.y();
diff = Point_2(x, y);
vec(0) = V(0, i);
vec(1) = V(1, i);
const auto res = inv * vec;
*(coordinates++) =
diff.x() * res(0) + diff.y() * res(1) + 1.0 / double(n);
}
};
// Compute affine coordinates for all query points.
std::cout << std::endl << "affine coordinates (all queries): " << std::endl << std::endl;
std::vector<double> coordinates;
coordinates.reserve(4);
for (const auto& query : queries) {
coordinates.clear();
affine(query, std::back_inserter(coordinates));
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}

Deprecated Coordinates

This example illustrates how to use the deprecated API of this package. The used Kernel is inexact and the used coordinates are mean value coordinates. The result is identical to the one from this example.


File Barycentric_coordinates_2/deprecated_coordinates.cpp

#include <CGAL/Installation/internal/disable_deprecation_warnings_and_errors.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
#include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
// Typedefs.
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
using Mean_value =
using Mean_value_coordinates =
int main() {
// Construct a star-shaped polygon.
const std::vector<Point_2> star_shaped = {
Point_2(0.0, 0.0), Point_2( 0.1, -0.8), Point_2(0.3, 0.0), Point_2(0.6, -0.5),
Point_2(0.6, 0.1), Point_2( 1.1, 0.6), Point_2(0.3, 0.2), Point_2(0.1, 0.8),
Point_2(0.1, 0.2), Point_2(-0.7, 0.0) };
// Construct the class with mean value coordinates
// for the star-shaped polygon defined above.
Mean_value_coordinates mean_value_coordinates(
star_shaped.begin(), star_shaped.end());
// Print some information about the polygon and coordinates.
mean_value_coordinates.print_information();
// Construct some interior points in the polygon.
const std::vector<Point_2> interior_points = {
Point_2(0.12, -0.45), Point_2(0.55, -0.3), Point_2(0.9 , 0.45),
Point_2(0.15, 0.35), Point_2(-0.4, 0.04), Point_2(0.11, 0.11),
Point_2(0.28, 0.12), // the only point in the kernel of the star shaped polygon
Point_2(0.55, 0.11) };
// We speed up the computation using the O(n) algorithm called with the
// parameter CGAL::Barycentric_coordinates::FAST.
// The default one is CGAL::Barycentric_coordinates::PRECISE.
const auto type_of_algorithm = CGAL::Barycentric_coordinates::FAST;
// We also speed up the computation by using the parameter
// query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE.
const auto query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE;
// Create a vector `std::vector` to store coordinates.
std::vector<FT> coordinates;
coordinates.reserve(star_shaped.size());
// Compute mean value coordinates for all interior points.
std::size_t count = 0;
for (const auto& query : interior_points) {
coordinates.clear();
const auto result = mean_value_coordinates(
query, std::back_inserter(coordinates), query_point_location, type_of_algorithm);
// Status of the computation.
const std::string status = (result ? "SUCCESS." : "FAILURE.");
std::cout << std::endl << "point: " << count << ", status of the computation: " << status << std::endl;
++count;
// Output the coordinates.
for (std::size_t i = 0; i < coordinates.size() - 1; ++i) {
std::cout << coordinates[i] << ", ";
}
std::cout << coordinates[coordinates.size() - 1] << std::endl;
}
std::cout << std::endl;
return EXIT_SUCCESS;
}
Note
The headers Segment_coordinates_2.h and Triangle_coordinates_2.h are not capitalized in the new version that is they are named segment_coordinates_2.h and triangle_coordinates_2.h.

Edge Cases

Not all presented coordinates are general enough to handle any query point in the plane, that is why we highly recommend reading this section in order to learn what can be expected from each coordinate function. If you want to get more mathematical details about each coordinate function as well as the complete history and theory behind barycentric coordinates, you should read [6]. You can also read an overview here (chapters 1 and 2).

Segment Coordinates

The segment coordinate function with respect to a given segment vertex is a linear function along the supporting line of this segment that grows from zero at the opposite vertex to one at the chosen vertex (see the figure below).

seg_coord.svg
Figure 100.12 The segment coordinate function with respect to the vertex \(v_0\).

Segment coordinates can be computed exactly if an exact number type is chosen. The segment itself, with respect to which we compute coordinates, must be non-degenerate. If both conditions are satisfied, then the computation never fails. However, to compute coordinates, the user must ensure that the query point lies exactly on the line \(L\) supporting the segment. Since in many applications this is not the case, and a query point may lie very close but not exactly on this line, we provide a solution to remedy this situation.

seg_coord_projection.svg
Figure 100.13 The orthogonal projection \(p'\) of the vector \(p\) (green) onto the vector \(q\) (red).

Suppose that some query point \(v\) does not lie exactly on the line \(L\), but is some distance \(d\) away as shown in the figure above. If we want to compute the segment coordinate \(b_1(v)\) with respect to the vertex \(v_1\), we first find the orthogonal projection \(p'\) of the vector \(p\) onto the vector \(q\) and then normalize it by the length of \(q\). This yields the segment coordinate \(b_1(v') = b_1(v)\) if \(v\) lies exactly on the line. The other segment coordinate \(b_0(v')\) that is equal to \(b_0(v)\) when \(v\) is on the line \(L\) is computed the same way but with the projection of the vector \(\vec{vv_1}\).

Warning: do not abuse the feature described above, because it does not yield correct segment coordinates for the point \(v\) but rather those for \(v'\). Moreover, segment coordinates for a point \(v\), which does not lie exactly on the line \(L\), do not exist. But if the non-zero distance \(d\) is due to some numerical instability when computing the location of the point \(v\) or any other problem, which causes the point to be not exactly on the line, the final segment coordinates will be, at least approximately, correct.

With inexact number types, the resulting coordinate values are correct up to the precision of the chosen type.

Triangle Coordinates

The triangle coordinate function with respect to a given triangle vertex is a linear function that grows from zero along the opposite edge to one at the chosen vertex (see the figure below).

tri_coord.svg
Figure 100.14 The triangle coordinate function with respect to the vertex \(v_0\).

To compute the triangle coordinates of the query point \(v\), we adopt the standard simple formula

\(b_i = \frac{A_i}{A}\) with \(i = 0\dots 2\)

where \(A_i\) is the signed area of the sub-triangle opposite to the vertex \(i\) and \(A\) is the total area of the triangle that is \(A = A_0 + A_1 + A_2\) (see the figure below).

tri_notations.svg
Figure 100.15 Notation for triangle coordinates.

These coordinates can be computed exactly if an exact number type is chosen, for any query point in the plane and with respect to any non-degenerate triangle. 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

Wachspress coordinates are well-defined in the closure of any strictly convex polygon. Therefore, when using an exact number type, for any query point from the polygon's closure, these coordinates are computed exactly and no false result is expected. For exterior query points, the coordinates can also be computed but not everywhere (see below for more details). For inexact number types, the resulting precision of the computation is due to the involved algorithm and a chosen number type. In the following paragraph, we discuss two available algorithms for computing Wachspress coordinates when an inexact number type is used. The chosen algorithm is specified by a computation policy Barycentric_coordinates::Computation_policy_2.

wp_notations.svg
Figure 100.16 Notation for Wachspress coordinates.

To compute Wachspress weights, we follow [2] and use the formula

\(w_i = \frac{C_i}{A_{i-1}A_i}\)

with \(i = 1\dots n\) where \(n\) is the number of polygon vertices. In order to compute the coordinates, we normalize these weights,

\(b_i = \frac{w_i}{W^{wp}}\qquad\) with \(\qquad W^{wp} = \sum_{j=1}^n w_j.\)

This formula becomes unstable when approaching the boundary of the polygon ( \(\approx 1.0e-10\) and closer). To fix the problem, we modify the weights \(w_i\) as

\(\bar{w}_i = C_i\prod_{j\not=i-1,i} A_j\).

After the above normalization, this gives us the precise algorithm to compute Wachspress coordinates but with \(O(n^2)\) performance only. The max speed \(O(n)\) algorithm uses the standard weights \(w_i\). Note that mathematically this modification does not change the coordinates. One should be cautious when using the unnormalized Wachspress weights. In that case, you must choose the \(O(n)\) type.

It is known that for strictly convex polygons the denominator's zero set of the Wachspress coordinates ( \(W^{wp} = 0~\)) is a curve, which (in many cases) lies quite far away from the polygon. More specifically, it interpolates the intersection points of the supporting lines of the polygon edges. Therefore, the computation of Wachspress coordinates outside the polygon is possible only at points that do not belong to this curve.

wp_zero_set.svg
Figure 100.17 Zero set (red) of the Wachspress coordinates' denominator \(W^{wp}\) for a non-regular hexagon.

Warning: we do not recommend using Wachspress coordinates for exterior points!

Discrete Harmonic Coordinates

Discrete harmonic coordinates have the same requirements as Wachspress coordinates. They are well-defined in the closure of any strictly convex polygon and, if an exact number type is chosen, they are computed exactly. However, and unlike Wachspress basis functions, these coordinates are not necessarily positive. In particular, the weight \(w_i\) is positive if and only if \(\alpha+\beta < \pi\) (see the figure below for notation). For inexact number types, the precision of the computation is due to the involved algorithm and a chosen number type. Again, we describe two algorithms to compute the coordinates when an inexact number type is used: one is of max precision and one is of max speed.

dh_notations.svg
Figure 100.18 Notation for discrete harmonic coordinates.

To compute discrete harmonic weights, we follow [2] and use the formula

\(w_i = \frac{r_{i+1}^2A_{i-1}-r_i^2B_i+r_{i-1}^2A_i}{A_{i-1}A_i}\)

with \(i = 1\dots n\) where \(n\) is the number of polygon vertices. In order to compute the coordinates, we normalize these weights,

\(b_i = \frac{w_i}{W^{dh}}\qquad\) with \(\qquad W^{dh} = \sum_{j=1}^n w_j.\)

This formula becomes unstable when approaching the boundary of the polygon ( \(\approx 1.0e-10\) and closer). To fix the problem, similarly to the previous subsection, we modify the weights \(w_i\) as

\(\bar{w}_i = (r_{i+1}^2A_{i-1}-r_i^2B_i+r_{i-1}^2A_i)\prod_{j\not=i-1,i} A_j\).

After the above normalization, this yields the precise algorithm to compute discrete harmonic coordinates but with \(O(n^2)\) performance only. The max speed \(O(n)\) algorithm uses the standard weights \(w_i\). Again, mathematically this modification does not change the coordinates, one should be cautious when using the unnormalized discrete harmonic weights. In that case, you must choose the \(O(n)\) type.

Warning: as for Wachspress coordinates, we do not recommend using discrete harmonic coordinates for exterior points, because the curve \(W^{dh} = 0\) may have several components, and one of them interpolates the polygon vertices. However, if you are sure that the query point does not belong to this curve, you can compute the coordinates as shown in this example.

Mean Value Coordinates

Unlike the previous coordinates, mean value coordinates cannot be computed exactly due to an inevitable 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(). It is worth saying that providing a number type that supports exact or nearly exact computation of the square root is possible, however since such types are usually impractical due to the large overhead, the conversion to a floating-point format above is always effective. On the other hand, mean value coordinates are well-defined everywhere in the plane for any simple polygon. In addition, if your traits class provides a more precise version of the square root function, the final precision of the computation with exact number types will depend only on the precision of that function.

mv_notations.svg
Figure 100.19 Notation for mean value coordinates.

For these coordinates, we provide two algorithms: one is of max precision and one is of max speed. The first one works everywhere in the plane, and the precision of the computation depends only on the chosen number type, including the remarks above. This algorithm is based on the following weight formula from [4]

\(w_i = \sigma_i\bar{w}_i\qquad\) with \(\qquad\bar{w}_i = (r_{i-1}r_{i+1}-d_{i-1}d_{i+1})^{1/2}\prod_{j\not= i-1,i}(r_jr_{j+1} + d_jd_{j+1})^{1/2}\qquad\) where \(\qquad r_i = \|d_i\|.\)

Since \(\bar{w}_i\) is always positive, we must append to it the proper sign \(\sigma_i\) of the signed mean value weight, which can be found efficiently (see the figures below). This weight is always positive to the left of the red piecewise linear curve, and it is negative to the right of this curve, moving in the counterclockwise direction.

mv_weight_signs.svg
Figure 100.20 Signs of the mean value weight \(w_i\) depending on the region with respect to a convex polygon \(P\) and a concave polygon \(P'\).

After the normalization of these weights as before

\(b_i = \frac{w_i}{W^{mv}}\qquad\) with \(\qquad W^{mv} = \sum_{j=1}^n w_j\)

we obtain the max precision \(O(n^2)\) algorithm. The max speed O(n) algorithm computes the weights \(w_i\) using the pseudocode from here. These weights

\(w_i = \frac{t_{i-1} + t_i}{r_i}\qquad\) with \(\qquad t_i = \frac{\text{det}(d_i, d_{i+1})}{r_ir_{i+1} + d_id_{i+1}}\)

are also normalized. Note that they are unstable if a query point is closer than \(\approx 1.0e-10\) to the polygon boundary, similarly to Wachspress and discrete harmonic coordinates and one should be cautious when using the unnormalized mean value weights. In that case, you must choose the \(O(n)\) type.

Harmonic Coordinates

The harmonic coordinates are computed by solving the Laplace equation

\(\Delta \boldsymbol{b} = \boldsymbol{0}\)

subject to suitable Dirichlet boundary conditions. Harmonic coordinates are the only coordinates in this package, which are guaranteed to be non-negative in the closure of any simple polygon and satisfy all properties of barycentric coordinates, however such desirable properties come with the fact that these coordinates are well-defined only inside a polygon. If an exterior query point is provided, its coordinates are set to zero.

Another disadvantage of these coordinates is that they cannot be computed exactly, because harmonic coordinates do not have a simple closed-form expression and must be approximated. The common way to approximate these coordinates is by discretizing over the space of piecewise linear functions with respect to a triangulation of the polygon. The denser triangulation of the interior part of the polygon, the better approximation of the coordinates. To get a high quality approximation of the coordinates, the user should provide a rather dense partition of the polygon's interior domain that in turn leads to larger running times when computing the coordinates.

terrain_triangulation.svg
Figure 100.21 Sparse triangulation of the polygon's interior domain (left): smaller running times, lower coordinates quality; dense triangulation (right): larger running times, higher coordinates quality.

From all this follows, that any exact Kernel will be rejected and it is not possible to compute analytic harmonic weights. However, once the coordinates are computed at the vertices of the triangulation, they can be evaluated analytically at any interior query point. For evaluation, we first locate a triangle that contains the query point and then linearly interpolate harmonic coordinates defined at the vertices of this triangle to the query point with the help of triangle coordinates as

\(b_i = b_0^{tr} b_i^0 + b_1^{tr} b_i^1 + b_2^{tr} b_i^2\)

with \(i = 1\dots n\), where \(n\) is the number of polygon vertices, \(b_{0}^{tr}\), \(b_{1}^{tr}\), and \(b_{2}^{tr}\) are the triangle coordinates of the query point with respect the three vertices of the located triangle, and \(b_i^{0}\), \(b_i^{1}\), and \(b_i^{2}\) are the harmonic coordinates pre-computed at the triangle vertices.

Performance

We strive for robustness and efficiency at the same time. Efficiency is especially important. These coordinates are used in many applications where they must be computed for millions of points and, thus, the real time computation of coordinates has been made possible. In this section, we present next the computation runtimes of the implemented algorithms.

The structure of the speed test that we use to evaluate the running times consists of computing coordinate values (or weights) at >= 1 million strictly interior points with respect to a polygon (or triangle, or segment). At each iteration of the loop, we create a query point and compute its coordinates. The time presented in the log-log scale plot at the end of the section is the arithmetic mean of all trials in the loop of 10 iterations. The time presented in the plot is for analytic coordinates only since harmonic coordinates of a reasonable (application-dependent) quality are substantially slower to compute and cannot be fairly compared to the analytic coordinate functions.

The time to compute coordinates depends on many factors such as memory allocation, input kernel, output container, number of points, etc. In our tests, we used the most standard C++ and CGAL features with minimum memory allocation. Therefore, the final time presented is the average time that can be expected without deep optimization but still with efficient memory allocation. It also means that it may vary depending on the usage of the package.

To benchmark analytic coordinates, we used a MacBook Pro 2011 with 2 GHz Intel Core i7 processor (2 cores) and 8 GB 1333 MHz DDR3 memory. The installed operating system was OS X 10.9 Maverick. The resulting timings for all closed-form coordinates can be found in the figure below.

analytic_timings.png
Figure 100.22 Time in seconds to compute \(n\) coordinate values for a polygon with \(n\) vertices at 1 million query points with the max speed \(O(n)\) algorithms (dashed) and the max precision \(0(n^2)\) algorithms (solid) for Wachspress (blue), discrete harmonic (red), and mean value (green) coordinates.

From the figure above we observe that the \(O(n^2)\) algorithm is as fast as the \(O(n)\) algorithm if we have a polygon with a small number of vertices. But as the number of vertices is increased, the linear algorithm outperforms the squared one, as expected. One of the reasons for this behavior is that for a small number of vertices the multiplications of \(n-2\) elements inside the \(O(n^2)\) algorithm take almost the same time as the corresponding divisions in the \(O(n)\) algorithm. For a polygon with many vertices, these multiplications are substantially slower.

To benchmark harmonic coordinates, we used a MacBook Pro 2018 with 2.2 GHz Intel Core i7 processor (6 cores) and 32 GB 2400 MHz DDR4 memory. The installed operating system was OS X 10.15 Catalina. The average time to compute harmonic coordinates in the loop of 10 iterations can be found in the tables below.

The first table shows how the time to compute the coordinates on a unit square depends on the number of triangulation vertices. We show separately the time to setup the matrix, factorize it, and solve it with respect to the four vertices of the unit square.

Number of queries (approx.) Setup (in seconds) Factorize (in seconds) Solve (in seconds) Total (in seconds)
100 0.000056 0.000099 0.000015 0.000170
500 0.000266 0.000574 0.000064 0.000904
1,000 0.000509 0.001194 0.000147 0.001850
25,000 0.014749 0.071152 0.008191 0.094092
50,000 0.034255 0.184237 0.018166 0.236658
100,000 0.065117 0.543177 0.044088 0.652382
500,000 0.576530 7.697143 0.310765 8.584438
1,000,000 1.295163 26.76945 0.737372 28.80199

The same results can be seen in the figure.

hm_4_bench.svg
Figure 100.23 Time in seconds to setup (red), factorize (green), and solve (blue) for harmonic coordinate values with respect to a unit square.

The second table shows how the time to compute the coordinates for 100k queries depends on the number of the polygon vertices. We show separately the time to setup the matrix, factorize it, and solve it with respect to the \(n\) vertices of the polygon. It can be seen that, unlike in the first table, the time to factorize the matrix here stays constant.

Number of vertices (approx.) Setup (in seconds) Factorize (in seconds) Solve (in seconds) Total (in seconds)
5 0.083444 0.631823 0.059827 0.775094
10 0.060294 0.450534 0.094583 0.605411
25 0.062760 0.478683 0.254953 0.796396
50 0.097359 0.492233 0.539654 1.129246
100 0.129487 0.450771 1.152544 1.732802
500 0.430694 0.460321 6.620061 7.511076
1000 0.812362 0.480052 16.14239 17.43480

The same results can be seen in the figure.

hm_n_bench.svg
Figure 100.24 Time in seconds to setup (red), factorize (green), and solve (blue) for harmonic coordinate values with respect to a polygon with \(n\) vertices at 100k query points.

While, in the first table, the most significant step is to factorize the matrix, in the second table, the slowest step is to solve for coordinates, as expected.

History

The package was first released in 2015 and included segment, triangle, Wachspress, discrete harmonic, and mean value coordinates. The API of that version is now deprecated but can still be used. An example of the old API can be found here. The docs of that API are also preserved and maintained here.

In 2018, this package was modified and improved by Keyu Chen and Dmitry Anisimov during the Google Summer of Code. The API was changed to the current version. In 2020, the new version was cleaned up and documented that includes:

  • the classes Segment_coordinates_2 and Triangle_coordinates_2 have been removed, only the free functions are preserved;
  • the entry class Generalized_barycentric_coordinates_2 was removed since it is not flexible enough to accommodate all types of 2D barycentric coordinates;
  • the classes Wachspress_2, Discrete harmonic_2, and Mean_value_2 have been renamed and modified so that they can be used now on their own without the class Generalized_barycentric_coordinates_2;
  • harmonic coordinates have been added;
  • the free functions for segment and triangle coordinates have been modified and improved;
  • the free functions for Wachspress, discrete harmonic, and mean value weights and coordinates have been added;
  • the free functions to compute barycentric coordinates for points on the polygon boundary have been added;
  • all functions and classes are now using ranges and property maps;
  • examples, tests, and benchmarks are modified/extended/improved;
  • the docs are refactored and simplified.

Acknowledgments

The authors wish to thank Teseo Schneider and Randolf Schaerfig for helpful comments and discussions. We also appreciate the great effort invested in this package by our reviewers Andreas Fabri, Sébastien Loriot, and Efi Fogel.