User Manual

Authors

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

Introduction

The package 2D Generalized Barycentric Coordinates offers an efficient and robust implementation of two-dimensional closed-form generalized barycentric coordinates defined for simple two-dimensional polygons. 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, the package includes an implementation of Wachspress, mean value, and discrete harmonic coordinates and provides some extra functions to compute barycentric coordinates with respect to segments (segment coordinates) and triangles (triangle coordinates). The section Theory of 2D Generalized Barycentric Coordinates gives a short introduction to the topic of barycentric coordinates.

Interface

Each class that computes barycentric coordinates is parameterized by a traits class. This traits class specifies types and geometric primitives that are used in the computation and must be a model of the concept BarycentricTraits_2.

The main entry point to the component is an input iterator over the vertices of a polygon. The polygon's vertices must follow clockwise or anticlockwise ordering and can be of any type. However, internally the classes use the type CGAL::Point_2, that is why an appropriate traits class that converts the user's type to CGAL::Point_2 must be provided. The same argument holds for query points.

Mean value 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.

Segment and triangle coordinates can be computed by using either a global function or creating the corresponding class. All other generalized coordinates can be computed by creating an instance of the class CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2 parameterized by an appropriate coordinate type that must be a model of the concept BarycentricCoordinates_2.

Any point in the plane may be taken as a query point. However, we do not recommend to use Wachspress and discrete harmonic coordinates with query points outside the closure of a polygon because at some of those points these coordinates are not well-defined, as explained in the Section Degeneracies and Special Cases.

Once instantiated for some polygon, the coordinates can be computed multiple times for different query points with respect to all the vertices of the provided polygon. Use the Reference Manual for the detailed interface.

The output of the computation is a set of coordinate values at the current query point with respect to all the vertices of the polygon. This output can be stored in an arbitrary container providing an appropriate output iterator. In addition, all the classes return a pointer to the last stored element and a status of the computation (Boolean true or false).

Examples

Segment Coordinates

This is a simple example to show the use of the global function CGAL::Barycentric_coordinates::compute_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. We use the exact kernel and return coordinates as an array of two values. Again, the symmetry of the query points helps us to recognize errors that may have occurred during the computation.

Figure 96.1 Example's point pattern.

File Barycentric_coordinates_2/Segment_coordinates_example.cpp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Segment_coordinates_2.h>

// Namespace alias.
namespace BC = CGAL::Barycentric_coordinates;

// Some convenient typedefs.
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::array<Scalar,2> Pair;

using std::cout; using std::endl; using std::string;

int main()
{
    // Construct a segment.
    const Point  first_vertex(0, Scalar(2)/Scalar(5));
    const Point second_vertex(2, Scalar(2)/Scalar(5));

    // Instantiate three interior and two exterior query points.
    const Point query_points[5] = { Point(Scalar(2) /Scalar(5), Scalar(2)/Scalar(5)), // interior query points
                                    Point(1                   , Scalar(2)/Scalar(5)),
                                    Point(Scalar(8) /Scalar(5), Scalar(2)/Scalar(5)),
                                    Point(Scalar(-1)/Scalar(5), Scalar(2)/Scalar(5)), // exterior query points
                                    Point(Scalar(11)/Scalar(5), Scalar(2)/Scalar(5))
                                  };

    // Compute segment coordinates for all the defined points.
    // We use a global function and return the segment coordinates stored in an array of the type std::array<FT,2>.
    cout << endl << "Computed segment coordinates: " << endl << endl;
    for(int i = 0; i < 5; ++i) {
        const Pair pair = BC::compute_segment_coordinates_2(first_vertex, second_vertex, query_points[i], Kernel());

        // Output both coordinates for each point.
        cout << "Pair of coordinates # " << i + 1 << " = (" << pair[0] << ", " << pair[1] << ");" << endl;
    }
    cout << endl;

    return EXIT_SUCCESS;
}

Triangle Coordinates

In this example we show how to use the class CGAL::Barycentric_coordinates::Triangle_coordinates_2 with the Simple_cartesian kernel for double type. 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. We use a standard container of the type std::vector and std::insert_iterator to access and store the resulting coordinate values.

Figure 96.2 Example's point pattern.

File Barycentric_coordinates_2/Triangle_coordinates_example.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_2/Triangle_coordinates_2.h>

// Some convenient typedefs.
typedef CGAL::Simple_cartesian<double> Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::vector<Scalar> Scalar_vector;

typedef CGAL::Barycentric_coordinates::Triangle_coordinates_2<Kernel> Triangle_coordinates;

using std::cout; using std::endl; using std::string;

int main()
{
    // Construct a triangle.
    const Point  first_vertex(0.0f, 0.0f);
    const Point second_vertex(2.0f, 0.5f);
    const Point  third_vertex(1.0f, 2.0f);

    // Create an std::vector to store coordinates.
    Scalar_vector coordinates;

    // Instantiate the class Triangle_coordinates_2 for the triangle defined above.
    Triangle_coordinates triangle_coordinates(first_vertex, second_vertex, third_vertex);

    // Print some information about the triangle and coordinates.
    triangle_coordinates.print_information();

    // Instantiate some interior, boundary, and exterior query points for which we compute coordinates.
    const int number_of_query_points = 18;
    const Point query_points[] = { Point(0.5f , 0.5f ), Point(1.0f, 0.5f ), Point(1.0f , 0.75f), Point(1.0f , 1.0f),                     // interior query points
                                   Point(1.0f , 1.25f), Point(1.0f, 1.5f ), Point(0.75f, 1.0f ), Point(1.25f, 1.0f), Point(1.5f, 0.75f),
                                   Point(1.0f , 0.25f), Point(0.5f, 1.0f ), Point(1.5f , 1.25f), Point(1.0f , 2.0f), Point(2.0f, 0.5f ), // boundary query points
                                   Point(0.25f, 1.0f ), Point(0.5f, 1.75f), Point(1.5f , 1.75f), Point(1.75f, 1.5f)                      // exterior query points
                                 };

    // Reserve memory to store triangle coordinates for 18 query points.
    coordinates.reserve(number_of_query_points * 3);

    // Compute triangle coordinates for these points.
    cout << endl << "Computed triangle coordinates: " << endl << endl;
    for(int i = 0; i < number_of_query_points; ++i) {

        triangle_coordinates(query_points[i], std::inserter(coordinates, coordinates.end()));

        // Output the coordinates for each point.
        cout << "Point " << i + 1 << ": ";
        for(int j = 0; j < 3; ++j)
            cout << "coordinate " << j + 1 << " = " << coordinates[i * 3 + j] << "; ";
        cout << endl << endl;
    }

    return EXIT_SUCCESS;
}

Wachspress Coordinates

In the following example we create 1000 random points, then we take the convex hull of this set of points as our polygon, and compute Wachspress coordinates at all the defined points. We use the Simple_cartesian kernel with double type as a traits class and store obtained coordinate values in a container of the type std::vector. The output iterator is std::back_insert_iterator.

File Barycentric_coordinates_2/Wachspress_coordinates_example.cpp

#include <CGAL/convex_hull_2.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Barycentric_coordinates_2/Wachspress_2.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>

// Some convenient typedefs.
typedef CGAL::Simple_cartesian<double> Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point>  Point_vector;

typedef CGAL::Creator_uniform_2<double, Point> Creator;

typedef CGAL::Barycentric_coordinates::Wachspress_2<Kernel> Wachspress;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Wachspress, Kernel> Wachspress_coordinates;

using std::cout; using std::endl; using std::string;

int main()
{
    // Choose how many random points we want to generate.
    const int number_of_points = 1000;

    // Create vectors to store generated points and vertices of a convex polygon.
    Point_vector points, vertices;

    // Generate a set of random points.
    CGAL::Random_points_in_square_2<Point,Creator> point_generator(1.0);
    std::copy_n(point_generator, number_of_points, std::back_inserter(points));

    // Find the convex hull of the generated set of points.
    // This convex hull gives the vertices of a convex polygon that contains all the generated points.
    CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(vertices));

    const size_t number_of_vertices = vertices.size();

    // Instantiate the class with Wachspress coordinates for the convex polygon defined above.
    Wachspress_coordinates wachspress_coordinates(vertices.begin(), vertices.end());

    // Print some information about the polygon and coordinates.
    wachspress_coordinates.print_information();

    // Compute Wachspress coordinates for all the randomly defined points.
    cout << endl << "Computed Wachspress coordinates: " << endl << endl;
    for(int i = 0; i < number_of_points; ++i) {
        // Compute coordinates.
        Scalar_vector coordinates;
        coordinates.reserve(number_of_vertices);
        wachspress_coordinates(points[i], std::back_inserter(coordinates));

        // Output the computed coordinates.
        cout << "Point " << i + 1 << ": " << endl;
        for(int j = 0; j < int(number_of_vertices); ++j)  cout << "Coordinate " << j + 1 << " = " << coordinates[j] << "; " << endl;
        cout << 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 show how to specify the location of a query point using additional function parameters. The used kernel is exact, and we use an output container of the type std::vector. Since all the points are symmetric, it is easy to debug the correctness of the obtained coordinate values. The output iterator is std::back_insert_iterator.

Figure 96.3 Example's point pattern.

File Barycentric_coordinates_2/Discrete_harmonic_coordinates_example.cpp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Discrete_harmonic_2.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>

// Some convenient typedefs.
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point>  Point_vector;

typedef std::back_insert_iterator<Scalar_vector> Vector_insert_iterator;
typedef boost::optional<Vector_insert_iterator> Output_type;

typedef CGAL::Barycentric_coordinates::Discrete_harmonic_2<Kernel> Discrete_harmonic;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Discrete_harmonic, Kernel> Discrete_harmonic_coordinates;

using std::cout; using std::endl; using std::string;

int main()
{
    // Construct a unit square.
    const int number_of_vertices = 4;
    Point_vector vertices(number_of_vertices);

    vertices[0] = Point(0, 0); vertices[1] = Point(1, 0); vertices[2] = Point(1, 1); vertices[3] = Point(0, 1);

    // Create an std::vector to store coordinates.
    Scalar_vector coordinates;

    // Instantiate the class with discrete harmonic coordinates for the unit square defined above.
    Discrete_harmonic_coordinates discrete_harmonic_coordinates(vertices.begin(), vertices.end());

    // Print some information about the polygon and coordinates.
    discrete_harmonic_coordinates.print_information();

    // Instantiate the center point of the unit square.
    const Point center(Scalar(1)/Scalar(2), Scalar(1)/Scalar(2));

    // Compute discrete harmonic coordinates for the center point.
    // Use the parameter query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE.
    Output_type result = discrete_harmonic_coordinates(center, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE);

    // Instantiate other 4 interior points.
    const int number_of_interior_points = 4;
    const Point interior_points[number_of_interior_points] = { Point(Scalar(1)/Scalar(5), Scalar(1)/Scalar(5))  ,
                                                               Point(Scalar(4)/Scalar(5), Scalar(1)/Scalar(5))  ,
                                                               Point(Scalar(4)/Scalar(5), Scalar(4)/Scalar(5))  ,
                                                               Point(Scalar(1)/Scalar(5), Scalar(4)/Scalar(5)) };

    // Compute discrete harmonic coordinates for these points and store them at the same vector "coordinates" as before.
    for(int i = 0; i < number_of_interior_points; ++i)
        result = discrete_harmonic_coordinates(interior_points[i], std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE);

    // Instantiate 2 boundary points on the second and last edges.
    const Point second_edge(1, Scalar(4)/Scalar(5));
    const Point   last_edge(0, Scalar(4)/Scalar(5));

    // Compute discrete harmonic coordinates for these 2 points.
    // Use the parameter query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDARY.
    result = discrete_harmonic_coordinates(second_edge, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDARY);
    result = discrete_harmonic_coordinates(last_edge  , std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDARY);

    // Instantiate 2 other boundary points on the first and third edges.
    const Point first_edge(Scalar(1)/Scalar(2), 0);
    const Point third_edge(Scalar(1)/Scalar(2), 1);

    // Compute discrete harmonic coordinates using index of an appropriate edge.
    // Do not forget that index counting starts from zero.
    result = discrete_harmonic_coordinates.compute_on_edge(first_edge, 0, std::back_inserter(coordinates));
    result = discrete_harmonic_coordinates.compute_on_edge(third_edge, 2, std::back_inserter(coordinates));

    // Compute discrete harmonic coordinates for the points at the first and third vertex of the unit square.
    result = discrete_harmonic_coordinates.compute_on_vertex(0, std::back_inserter(coordinates));
    result = discrete_harmonic_coordinates.compute_on_vertex(2, std::back_inserter(coordinates));

    // Instantiate points at the second and fourth vertex of the unit square.
    const Point second_vertex(1, 0);
    const Point fourth_vertex(0, 1);

    // Compute discrete harmonic coordinates for these points.
    // Use the parameter query_point_location = CGAL::Barycentric_coordinates::ON_VERTEX.
    result = discrete_harmonic_coordinates(second_vertex, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_VERTEX);
    result = discrete_harmonic_coordinates(fourth_vertex, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_VERTEX);

    // Instantiate 2 points outside the unit square - one from the left and one from the right.
    const Point left_most(Scalar(-1)/Scalar(2), Scalar(1)/Scalar(2));
    const Point right_most(Scalar(3)/Scalar(2), Scalar(1)/Scalar(2));

    // Compute discrete harmonic coordinates for these 2 points.
    // Use the parameter query_point_location = CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE.
    result = discrete_harmonic_coordinates(left_most , std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE);
    result = discrete_harmonic_coordinates(right_most, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE);

    // Output the computed coordinate values.
    cout << endl << "Exact discrete harmonic coordinates for all the defined points: " << endl << endl;
    const size_t number_of_query_points = coordinates.size();
    for(int index = 0; index < int(number_of_query_points); ++index) {
        cout << "Coordinate " << index % number_of_vertices + 1 << " = " << coordinates[index] << " ";
        if((index +  1) %      number_of_vertices  == 0) cout << endl;
        if((index + 13) % (4 * number_of_vertices) == 0) cout << endl;
    }

    // Return status of the last computation.
    const string status = (result ? "SUCCESS." : "FAILURE.");
    cout << endl << "Status of the last computation: " << status << endl << endl;

    return EXIT_SUCCESS;
}

Mean Value Coordinates

This is an example that shows 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 give negative coordinate values for points outside the polygon's kernel (shown in red). We use an inexact data type, an output container of the type std::vector, and an output iterator of the type std::back_insert_iterator to compute, access, and store the resulting coordinate values. We also show how to choose different algorithms to compute generalized barycentric coordinates (one is more precise while the other is faster).

Figure 96.4 Example's point pattern.

File Barycentric_coordinates_2/Mean_value_coordinates_example.cpp

#include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>

// Some convenient typedefs.
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point>  Point_vector;

typedef std::back_insert_iterator<Scalar_vector> Vector_insert_iterator;
typedef boost::optional<Vector_insert_iterator> Output_type;

typedef CGAL::Barycentric_coordinates::Mean_value_2<Kernel> Mean_value;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Mean_value, Kernel> Mean_value_coordinates;

using std::cout; using std::endl; using std::string;

int main()
{
    // Construct a star-shaped polygon.
    const int number_of_vertices = 10;
    Point_vector vertices(number_of_vertices);

    vertices[0] = Point(0.0, 0.0); vertices[1] = Point(0.1, -0.8); vertices[2] = Point(0.3, 0.0); vertices[3] = Point(0.6, -0.5); vertices[4]  = Point(0.6 , 0.1);
    vertices[5] = Point(1.1, 0.6); vertices[6] = Point(0.3,  0.2); vertices[7] = Point(0.1, 0.8); vertices[8] = Point(0.1,  0.2); vertices[9] = Point(-0.7, 0.0);

    // Create an std::vector to store coordinates.
    Scalar_vector coordinates;

    // Instantiate the class with mean value coordinates for the polygon defined above.
    Mean_value_coordinates mean_value_coordinates(vertices.begin(), vertices.end());

    // Print some information about the polygon and coordinates.
    mean_value_coordinates.print_information();

    // Instantiate some interior points in the polygon.
    const int number_of_interior_points = 8;
    const Point interior_points[] = { Point(0.12, -0.45), Point(0.55, -0.3), Point(0.9 , 0.45),
                                      Point(0.15,  0.35), Point(-0.4, 0.04), Point(0.11, 0.11),
                                      Point(0.28,  0.12), // the only point in the kernel of the star shaped polygon
                                      Point(0.55,  0.11)
                                    };

    // Compute mean value coordinates for all the defined interior points.

    // 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 CGAL::Barycentric_coordinates::Type_of_algorithm 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 CGAL::Barycentric_coordinates::Query_point_location query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE;

    for(int i = 0; i < number_of_interior_points; ++i) {
        const Output_type result = mean_value_coordinates(interior_points[i], std::back_inserter(coordinates), query_point_location, type_of_algorithm);

        // Output the coordinates for each point.
        const string status = (result ? "SUCCESS." : "FAILURE.");
        cout << endl << "For the point " << i + 1 << " status of the computation: " << status << endl;

        for(int j = 0; j < number_of_vertices; ++j)
            cout << "Coordinate " << j + 1 << " = " << coordinates[i * number_of_vertices + j] << endl;
    }

    // If we need only the unnormalized weights for some point (lets take the last one), we can compute them as follows.

    // Instantiate an std::vector to store weights.
    Scalar_vector weights;

    // Compute mean value weights.
    const int last_point_index = number_of_interior_points - 1;
    const Output_type result = mean_value_coordinates.compute_weights(interior_points[last_point_index], std::back_inserter(weights));

    // Compute their sum.
    Scalar mv_denominator = Scalar(0);
    for(int j = 0; j < number_of_vertices; ++j) mv_denominator += weights[j];

    // Invert this sum.
    const Scalar mv_inverted_denominator = Scalar(1) / mv_denominator;

    // Output the mean value weights.
    const string status = (result ? "SUCCESS." : "FAILURE.");
    cout << endl << "Status of the weights' computation for the point " << last_point_index + 1 << ": " << status << endl;

    for(int j = 0; j < number_of_vertices; ++j)
        cout << "Weight " << j + 1 << " = " << weights[j] << endl;

    // Now, if we normalize the weights, we recover values of the mean value coordinates for the last point computed earlier.
    cout << endl << "After normalization, for the point " << last_point_index + 1 << " mean value coordinates are " << endl;
    for(int j = 0; j < number_of_vertices; ++j)
        cout << "Coordinate " << j + 1 << " = " << weights[j] * mv_inverted_denominator << endl;
    cout << endl;

    return EXIT_SUCCESS;
}

Height Interpolation for Terrain Modeling

This is an advanced example that shows 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 gives 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.

Figure 96.5 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 three-dimensional polygon orthogonally onto the two-dimensional plane using the class CGAL::Projection_traits_xy_3, triangulate its interior using the class CGAL::Delaunay_mesher_2, and compute mean value coordinates for all the obtained points with respect to all the polygon's vertices. Finally, we interpolate the height data from the polygon's 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/Delaunay_mesher_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>

// Some convenient typedefs.

// General.
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Projection_traits_xy_3<Kernel>                Projection;

typedef Projection::FT      Scalar;
typedef Projection::Point_2 Point;

typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point>  Point_vector;

// Coordinates related.
typedef CGAL::Barycentric_coordinates::Mean_value_2<Projection>                                      Mean_value;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Mean_value, Projection> Mean_value_coordinates;

// Triangulation related.
typedef CGAL::Delaunay_mesh_face_base_2<Projection>                  Face_base;
typedef CGAL::Triangulation_vertex_base_2<Projection>                Vertex_base;
typedef CGAL::Triangulation_data_structure_2<Vertex_base, Face_base> TDS;
typedef CGAL::Constrained_Delaunay_triangulation_2<Projection, TDS>  CDT;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>                     Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria>                       Mesher;
typedef CDT::Finite_vertices_iterator                                Vertex_iterator;
typedef CDT::Vertex_handle                                           Vertex_handle;

// Interpolation related.
typedef CGAL::Interpolation_traits_2<Projection>                             Interpolation_traits;
typedef CGAL::Data_access< std::map<Point, Scalar, Projection::Less_xy_2 > > Value_access;

// STD.
using std::cout; using std::endl; using std::string;

int main()
{
    // Construct a polygon bounding a piece of three-dimensional terrain.
    // Note that z-coordinate of each vertex represents the height function.
    // Projection in 2D is done automatically by the Projection traits class.
    const int number_of_vertices = 50;
    Point_vector vertices(number_of_vertices);

    vertices[0]  = Point(0.03, 0.05, 0.000); vertices[1]  = Point(0.07, 0.04, 10.00); vertices[2]  = Point(0.10, 0.04, 20.00);
    vertices[3]  = Point(0.14, 0.04, 30.00); vertices[4]  = Point(0.17, 0.07, 40.00); vertices[5]  = Point(0.19, 0.09, 50.00);
    vertices[6]  = Point(0.22, 0.11, 60.00); vertices[7]  = Point(0.25, 0.11, 70.00); vertices[8]  = Point(0.27, 0.10, 80.00);
    vertices[9]  = Point(0.30, 0.07, 90.00); vertices[10] = Point(0.31, 0.04, 100.0); vertices[11] = Point(0.34, 0.03, 110.0);
    vertices[12] = Point(0.37, 0.02, 120.0); vertices[13] = Point(0.40, 0.03, 130.0); vertices[14] = Point(0.42, 0.04, 140.0);
    vertices[15] = Point(0.44, 0.07, 150.0); vertices[16] = Point(0.45, 0.10, 160.0); vertices[17] = Point(0.46, 0.13, 170.0);
    vertices[18] = Point(0.46, 0.19, 180.0); vertices[19] = Point(0.47, 0.26, 190.0); vertices[20] = Point(0.47, 0.31, 200.0);
    vertices[21] = Point(0.47, 0.35, 210.0); vertices[22] = Point(0.45, 0.37, 220.0); vertices[23] = Point(0.41, 0.38, 230.0);
    vertices[24] = Point(0.38, 0.37, 240.0); vertices[25] = Point(0.35, 0.36, 250.0); vertices[26] = Point(0.32, 0.35, 260.0);
    vertices[27] = Point(0.30, 0.37, 270.0); vertices[28] = Point(0.28, 0.39, 280.0); vertices[29] = Point(0.25, 0.40, 290.0);
    vertices[30] = Point(0.23, 0.39, 300.0); vertices[31] = Point(0.21, 0.37, 310.0); vertices[32] = Point(0.21, 0.34, 320.0);
    vertices[33] = Point(0.23, 0.32, 330.0); vertices[34] = Point(0.24, 0.29, 340.0); vertices[35] = Point(0.27, 0.24, 350.0);
    vertices[36] = Point(0.29, 0.21, 360.0); vertices[37] = Point(0.29, 0.18, 370.0); vertices[38] = Point(0.26, 0.16, 380.0);
    vertices[39] = Point(0.24, 0.17, 390.0); vertices[40] = Point(0.23, 0.19, 400.0); vertices[41] = Point(0.24, 0.22, 410.0);
    vertices[42] = Point(0.24, 0.25, 420.0); vertices[43] = Point(0.21, 0.26, 430.0); vertices[44] = Point(0.17, 0.26, 440.0);
    vertices[45] = Point(0.12, 0.24, 450.0); vertices[46] = Point(0.07, 0.20, 460.0); vertices[47] = Point(0.03, 0.15, 470.0);
    vertices[48] = Point(0.01, 0.10, 480.0); vertices[49] = Point(0.02, 0.07, 490.0);

    // Mesh this polygon.

    // Create a constrained Delaunay triangulation.
    CDT cdt;

    std::vector<Vertex_handle> vertex_handle(number_of_vertices);

    // Insert vertices of the polygon as our initial point set.
    for(int i = 0; i < number_of_vertices; ++i) vertex_handle[i] = cdt.insert(vertices[i]);

    // Insert constraints - edges of the polygon - in order to mesh only the polygon's interior.
    for(int i = 0; i < number_of_vertices; ++i) cdt.insert_constraint(vertex_handle[i], vertex_handle[(i + 1) % number_of_vertices]);

    Mesher mesher(cdt);

    // Set a criteria on how to mesh.
    mesher.set_criteria(Criteria(0.01, 0.01));

    // Mesh the polygon.
    mesher.refine_mesh();

    // Compute mean value coordinates and use them to interpolate data from the polygon's boundary to its interior.

    // Associate each point with the corresponding function value and coordinates.
    std::map<Point, Scalar, Projection::Less_xy_2> point_function_value;
    std::vector< std::pair<Point, Scalar> >        point_coordinates(number_of_vertices);

    for(int i = 0; i < number_of_vertices; ++i)
            point_function_value.insert(std::make_pair(vertices[i], vertices[i].z()));

    // Create an instance of the class with mean value coordinates.
    Mean_value_coordinates mean_value_coordinates(vertices.begin(), vertices.end());

    // Store all new interior points with interpolated data here.
    std::vector<Point> points(cdt.number_of_vertices());

    cout << endl << "Result of the height interpolation: " << endl << endl;

    // Compute coordinates and interpolate the boundary data to the polygon's interior.
    int index = 0;
    for(Vertex_iterator vertex_iterator = cdt.finite_vertices_begin(); vertex_iterator != cdt.finite_vertices_end(); ++vertex_iterator) {
        Scalar_vector coordinates;

        const Point &point = vertex_iterator->point();
        mean_value_coordinates(point, std::back_inserter(coordinates));

        for(int j = 0; j < number_of_vertices; ++j)
            point_coordinates[j] = std::make_pair(vertices[j], coordinates[j]);

        Scalar f = CGAL::linear_interpolation(point_coordinates.begin(), point_coordinates.end(), Scalar(1), Value_access(point_function_value));
        points[index] = Point(point.x(), point.y(), f);
        cout << "The interpolated height with index " << index << " is " << f << ";" << endl;
        ++index;
    }
    cout << endl;

    return EXIT_SUCCESS;
}

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

Figure 96.6 The interpolated data. The color bar represents the corresponding height.

Degeneracies and Special Cases

Segment Coordinates

Segment coordinates can be computed exactly if an exact data 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 be sure about the query point being 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, the class is also able to handle this situation.

Figure 96.7 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 barycentric 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 gives the segment barycentric coordinate \(b_1(v') = b_1(v)\) if \(v\) lies exactly on the line.

Warning: do not abuse the feature described above because it does not give correct segment barycentric coordinates for the point \(v\) but rather those for \(v'\). Moreover, segment barycentric 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 data types, the resulting coordinate values are correct up to the precision of the chosen type.

Triangle Coordinates

These coordinates can be computed exactly if an exact data 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 gives the correct result. The notion of correctness depends on the precision of the used data 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, for any query point from the polygon's closure with an exact data type, these coordinates are computed exactly and no false result is expected. For inexact data types, the resulting precision of the computation is due to the involved algorithm and chosen data type. In the following paragraph we discuss two available algorithms for computing Wachspress coordinates. One of them is CGAL::Barycentric_coordinates::PRECISE, the other is CGAL::Barycentric_coordinates::FAST.

Figure 96.8 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 the polygon's 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 normalization as above, this gives us the precise algorithm to compute Wachspress coordinates but with \(O(n^2)\) performance only. The fast \(O(n)\) algorithm uses the standard weights \(w_i\). Note that mathematically this modification does not change the coordinates.

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. Speaking precisely, it interpolates the intersection points of the continuations of the polygon's edges. Therefore, the computation of Wachspress coordinates outside the polygon is possible only at points that do not belong to this curve.

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

Warning: we do not recommend to use 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 data type is chosen, they are computed exactly. But, 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 data types, the precision of the computation is due to the involved algorithm and chosen data type. Again, we describe two algorithms to compute the coordinates: one is precise and one is fast.

Figure 96.10 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 the polygon's 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 normalization as above, this gives the precise algorithm to compute discrete harmonic coordinates but with \(O(n^2)\) performance only. The fast \(O(n)\) algorithm uses the standard weights \(w_i\). Again, mathematically this modification does not change the coordinates.

Warning: as for Wachspress coordinates, we do not recommend to use discrete harmonic coordinates for exterior points because the curve \(W^{dh} = 0\) may have several components, and one of them interpolates the polygon's 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 data type is used, the default precision of the computation depends only on two CGAL functions: CGAL::to_double() and CGAL::sqrt(). 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 data types will depend only on the precision of that function.

Figure 96.11 Notation for mean value coordinates.

For these coordinates we also have two algorithms: one is precise and one is fast. The first one works everywhere in the plane, and the precision of the computation depends only on the chosen data 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 have to append to it the proper sign \(\sigma_i\) of the signed mean value weight, which can be found efficiently (see the figures below). Basically, 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 anticlockwise direction.

Figure 96.12 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 precise \(O(n^2)\) algorithm. The fast 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's boundary, similarly to Wachspress and discrete harmonic coordinates.

Performance

Apart from the most important requirement on barycentric coordinates to be as precise as possible, it is very important for them to be as fast as possible to evaluate. These coordinates are used in many applications where they must be computed for millions of points and, thus, the real time usage of coordinates is crucial. When writing the code, we tried to fulfill this important requirement, and in this section we present a few results about the computation times of the implemented coordinates.

The structure of the speed test that we ran for all functions consists of computing coordinate values (or weights) at >= 1 million strictly interior points with respect to some polygon (or triangle, or segment). At each iteration of the loop we create a query point, pass it to the function, and compute all the related coordinates. We run this loop 10 times in a row, and the time presented in the log-log scale plot at the end of the section is the arithmetic mean of all trials.

A typical example of this performance test for triangle coordinates with reduced number of query points can be found below. This example also illustrates how to construct an iterator and pass it to the class. In this example we create an iterator that writes coordinate values for each new query point over coordinate values of the previous point in the fixed-size standard C++ array, so that memory is allocated only once.

File Barycentric_coordinates_2/Triangle_coordinates_speed_test.cpp

#include <CGAL/Real_timer.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Triangle_coordinates_2.h>

// Construct an iterator that takes as input the current data type and pointer to the first element in the standard C++ array.
template<typename Scalar>
    class overwrite_iterator
{
private:
    Scalar* pointer;

public:
    explicit overwrite_iterator(Scalar* new_pointer) : pointer(new_pointer) { }

    // There are only two operations that we need to overload in order to use the class Triangle_coordinates_2.

    // This operation is intended to return the current coordinate value.
    inline Scalar& operator* () { return *pointer; }

    // This operation is intended to increase the index of the coordinate.
    inline void operator++ () { ++pointer; }
};

// Some convenient typedefs.
typedef CGAL::Real_timer Timer;

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef overwrite_iterator<Scalar> Overwrite_iterator;

typedef CGAL::Barycentric_coordinates::Triangle_coordinates_2<Kernel> Triangle_coordinates;

using std::cout; using std::endl; using std::string;

int main()
{
    // Number of x and y coordinates together gives the number of points.
    const int number_of_x_coordinates = 100000;
    const int number_of_y_coordinates = 1000;

    // Number of runs to compute the arithmetic mean of the time.
    const int number_of_runs = 10;

    // Compute the uniform step size along x and y directions to change coordinates.
    const Scalar zero = Scalar(0);
    const Scalar one  = Scalar(1);
    const Scalar two  = Scalar(2);
    const Scalar x_step = one / Scalar(number_of_x_coordinates);
    const Scalar y_step = one / Scalar(number_of_y_coordinates);

    // Create a right triangle with a slight offset from zero.
    const Point  first_vertex(zero - x_step, zero - x_step);
    const Point second_vertex(two  + y_step, zero - x_step);
    const Point  third_vertex(zero - x_step, two  + y_step);

    // Instantiate the class Triangle_coordinates_2 for the right triangle defined above.
    Triangle_coordinates triangle_coordinates(first_vertex, second_vertex, third_vertex);

    // Create an instance of the standard C++ array to store coordinates.
    // It has the fixed size = 3 = number of vertices.
    Scalar coordinates [3] = {0, 0, 0};

    // Pass pointer to the first element of the array with coordinates in order to overwrite them.
    Overwrite_iterator it( &(coordinates[0]) );

    // Create a timer.
    Timer time_to_compute;

    double time = 0.0;
    for(int i = 0; i < number_of_runs; ++i) { // Number of runs

        time_to_compute.start(); // Start clock
        for(Scalar x = zero; x <= one; x += x_step) {
            for(Scalar y = zero; y <= one; y += y_step)
                triangle_coordinates(Point(x, y), it); // Compute 3 coordinate values for each generated point
        }
        time_to_compute.stop(); // Stop clock

        time += time_to_compute.time(); // Add time of the current run to the whole time

        time_to_compute.reset(); // Reset time
    }

    // Compute the arithmetic mean of all the runs.
    const double mean_time = time / number_of_runs;

    // Output the resulting time.
    cout.precision(10);
    cout << endl << "CPU time to compute triangle coordinates for "
         << number_of_x_coordinates * number_of_y_coordinates << " points = " << mean_time << " seconds.";
    cout << endl << endl;

    return EXIT_SUCCESS;
}

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.

For all tests 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. In order to compile the speed test suite, we used the Clang 5.0 64bit compiler. The resulting timings can be found in the figure below.

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

From the figure above it is easy to see 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 behaviour is that for a small number of vertices the multiplications for \(n-2\) elements inside the \(O(n^2)\) algorithm with the fast \(O(n)\) algorithms (dashed) and the slow \(O(n^2)\) algorithms (solid) take almost the same time as the corresponding division in the \(O(n)\) algorithm. For a polygon with many vertices this multiplication is much slower.

Implementation Details

The generic design of the package was developed in 2013 by Dmitry Anisimov and David Bommes with many useful comments by Kai Hormann and Pierre Alliez. The package consists of 6 classes, 2 enumerations, and one namespace. Appropriate iterators are used to provide an efficient access to data and to pass them to one of the generic algorithms for computing coordinates. Once instantiated for a polygon (triangle, segment), the coordinates can be computed multiple times for different query points with respect to all the vertices of the provided polygon (triangle, segment). All the classes are fully templated and have a simple and similar design. In particular, we follow the same naming convention for all functions. Yet, the number of functions can differ from one class to another.

The implemented algorithms for computing coordinates do not depend on a particular kernel, and all the coordinates can be computed exactly, if an exact kernel is used, apart from mean value coordinates. The latter coordinates involve a square root operation, which results in a slightly worse precision with exact data types due to temporal conversion into a floating point type. The computed coordinates can be stored in an arbitrary container if an appropriate output iterator is provided.

It is worth noting that the class CGAL::Barycentric_coordinates::Segment_coordinates_2 is used to compute generalized barycentric coordinates along the polygon's boundary. Hence, one can use the trick for segment coordinates from Section Degeneracies and Special Cases if one is convinced that a point must lie exactly on the polygon's boundary but due to some numerical instabilities it does not.

The package is implemented in a way that later, if needed, other two-dimensional generalized barycentric coordinates can be easily added to this package.

Theory of 2D Generalized Barycentric Coordinates

In 1827, the German mathematician and theoretical astronomer August Ferdinand Möbius (1790–1868) proposed a method [7] to find coordinates of a point in the plane with respect to the vertices of a triangle. These coordinates are called triangle barycentric coordinates (sometimes area coordinates), and they are widely used in a variety of applications. Some of these applications are linear interpolation over a triangle and a triangle inclusion test. The first one is used for so-called shading, and the second one arises in the rasterization step when an image in vector graphics format needs to be converted into a raster image.

Triangle barycentric coordinates have many important properties, including constant and linear precision, the Lagrange property, and positivity inside a triangle. These properties make these coordinates a unique tool in many scientific fields. If we restrict triangle coordinates to one of the edges of a triangle and its supporting line, we get barycentric coordinates with respect to a segment and call them segment coordinates.

Let us show a couple of plots for the coordinates described above. To plot segment coordinates, we take a line \(y = 0.4\) and define a segment \([v_0, v_1]\) on this line. Then we sample this segment and compute segment coordinates for all the sample points. If we plot the segment coordinate function at all the defined points with respect to the vertex \(v_1\), we get the blue line depicted in the figure below. It grows from zero at the vertex \(v_0\) to one at the vertex \(v_1\).

Figure 96.14 Segment coordinates (blue) for all the segment points (green) with respect to the vertex \(v_1 = (2.0,\ 0.4)\).

If we want to plot triangle coordinates, we follow a similar approach. We take a triangle \([v_0, v_1, v_2]\) in the plane and sample its interior and boundary with a number of points. Once we have this sampling, we plot one of the triangle coordinate functions (here with respect to the third vertex of the triangle) at all the defined sample points. Likewise, we can plot the coordinate function with respect to the first or second vertex. The resulting function is linear (shown in the figure below) that grows from zero along the first edge \([v_0, v_1]\) to one at the chosen vertex \(v_2\).

Figure 96.15 Triangle coordinates with respect to \(v_2 = (1.0,\ 2.0)\). The color bar indicates the range of values for the chosen coordinate.

Since many applications require to work with more complex planar geometric shapes than segments and triangles, it seems natural to investigate a generalized version of triangle coordinates with respect to arbitrary polygons. The first attempt was taken in 1975 by E. L. Wachspress [9], and the resulting generalized barycentric coordinates are now called Wachspress coordinates [6]. These coordinates are well-defined for arbitrary strictly convex polygons and have all the properties of triangle coordinates [2]. Unfortunately, they are not well-defined for weakly convex and concave polygons.

Analogously to the previous cases, we want to plot the Wachspress coordinates and see how they look like. Let us choose a non-regular hexagon, slightly rotate it, and move one of its vertices towards the line through its two adjacent neighbors. We sample the interior and the boundary of this polygon as before and plot the coordinate function with respect to the vertex that we moved at all the sample points. We see that we get a smooth function, which is linear along all edges and grows from zero to one, as the color bar indicates.

Figure 96.16 The Wachspress coordinate function with respect to the indicated vertex with values from zero to one as the color bar indicates.

Another type of generalized barycentric coordinates goes back to Pinkall and Polthier in 1993 [8] and Eck et al. in 1995 [1] in the context of triangle mesh parameterization. They are called discrete harmonic coordinates. These coordinates are well-defined, similarly to Wachspress coordinates, for arbitrary strictly convex polygons and inherit all the properties of triangle coordinates apart from the positivity inside a polygon because they can take on negative values for some polygons. Another interesting property of these coordinate functions is that they coincide with Wachspress coordinates for any polygon whose vertices lie on a common circle.

To plot discrete harmonic coordinates we take the same polygon as for Wachspress coordinates and plot the coordinate function with respect to the same vertex. Again, we get a smooth function, which is linear along all edges and grows from zero to one. Isolines in the plot show the difference between discrete harmonic and Wachspress coordinates for the chosen polygon and vertex.

Figure 96.17 The discrete harmonic coordinate function with respect to the indicated vertex with values from zero to one as the color bar indicates.

The last type of generalized barycentric coordinates that we discuss are mean value coordinates [3] proposed by M. Floater in 2003. Based on the mean value theorem, these coordinates, unlike Wachspress and discrete harmonic coordinates, are well-defined for arbitrary simple polygons, inherit all the properties of triangle coordinates for any convex polygon, and lack only the positivity property for general concave polygons. Hormann and Floater prove in [5] that these coordinates are positive inside the kernel of a star-shaped polygon. They are also positive in the closure of any quadrilateral. Like discrete harmonic weights, mean value weights are often used in the context of triangle mesh parameterization.

In order to show the particular behaviour of mean value coordinates with an application to concave polygons, we take a star-shaped polygon with ten vertices \([v_0, \dots, v_9]\), sample its interior and boundary, and plot the coordinate function with respect to the fourth vertex \(v_3\). As the color bar indicates, the obtained function grows from a slightly negative value to one at the chosen vertex. It is also smooth inside the polygon and linear along all edges.

Figure 96.18 Mean value coordinates with respect to \(v_3\). The color bar indicates the range of values for the chosen coordinate function.

Interesting fact: all the coordinates discussed in this section and implemented in the package come from one and the same family of generalized barycentric coordinates named 3-point family of coordinates [2].

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 and Sébastien Loriot. Finally, to create pictures for this manual, we used two programs: Geogebra and Matlab.