pybind11-wrapped CGAL

__init__.pyi

from typing import Any, List, Tuple

import numpy as np

ArrayNx2F64 = Any  # Actually, should be np.ndarray but mypy doesn't like it
ArrayNx3F64 = Any
ArrayNx3I32 = Any
ArrayNx2I32 = Any
ArrayNI32 = Any

Subtriangulation = Tuple[ArrayNx2F64, ArrayNx3I32, List[List[int]]]

def triangulate(barycentric_points: ArrayNx2F64, segment_indices: ArrayNx2I32) -> Subtriangulation:
    """
    Triangulate a set of barycentric coordinates, requiring edge-paths exist between points of specified indices.

    Returns a tuple of the form (new_point_coordinates, triangles, segments), where
    * `new_point_coordinates` contains the barycentric coordinates of *newly added points* (at intersections)
    * `triangles` contains the vertex indices of the triangulation, where the indices are obtained by stacking
      `barycentric_points` beneath `new_point_coordinates`
    * `segments` contains lists of vertex indices correcting the segments specified by `segment_indices`
    """
    ...

class Mesh:
    def __init__(self, xyz: ArrayNx3F64, faces: ArrayNx3I32): ...
    def closest_faces(self, xyz: ArrayNx3F64) -> ArrayNI32: ...
    def sector(self, vertex_halfedges: ArrayNx2I32) -> ArrayNI32: ...
    def __repr__(self) -> str: ...

arrays.cpp

#include "triangulation.h"

namespace cgal_tz {
    np_array_double double_array_2d(std::vector<double> data, size_t n_cols) {
        size_t n_rows = data.size() / n_cols;
        if (n_rows * n_cols != data.size()) {
            throw std::runtime_error("n_rows * n_cols != data.size()");
        }
        size_t n_dim{2};
        std::vector<size_t> shape{n_rows, n_cols};
        std::vector<size_t> strides{sizeof(double) * n_cols, sizeof(double)};
        return py::array(py::buffer_info(
                data.data(),
                sizeof(double),
                py::format_descriptor<double>::format(),
                n_dim,
                shape,
                strides));
    }

    np_array_int int_array_2d(std::vector<int> data, size_t n_cols) {
        size_t n_rows = data.size() / n_cols;
        if (n_rows * n_cols != data.size()) {
            throw std::runtime_error("n_rows * n_cols != data.size()");
        }
        size_t n_dim{2};
        std::vector<size_t> shape{n_rows, n_cols};
        std::vector<size_t> strides{sizeof(int) * n_cols, sizeof(int)};
        return py::array(py::buffer_info(
                data.data(),
                sizeof(int),
                py::format_descriptor<int>::format(),
                n_dim,
                shape,
                strides));
    }

    np_array_int int_array_1d(std::vector<int> data) {
        size_t n_dim{1};
        std::vector<size_t> shape{data.size()};
        std::vector<size_t> strides{sizeof(int)};
        return py::array(py::buffer_info(
                data.data(),
                sizeof(int),
                py::format_descriptor<int>::format(),
                n_dim,
                shape,
                strides));
    }
}

arrays.h

#ifndef CGAL_TZ_ARRAYS_H
#define CGAL_TZ_ARRAYS_H

#include "common.h"


namespace cgal_tz {
    np_array_int int_array_2d(std::vector<int> data, size_t n_cols);
    np_array_int int_array_1d(std::vector<int> data);
    np_array_double double_array_2d(std::vector<double> data, size_t n_cols);
}

#endif

cgal_tz_module.cpp

#include "cgal_tz_module.h"

#include "mesh.h"
#include "triangulation.h"

void pybind_mesh(py::module &m) {
    using cgal_tz::Mesh;
    py::class_<Mesh> mesh(m, "Mesh");
    mesh
            .def(py::init([](const np_array_double &xyz, const np_array_int &faces) {
                return std::make_unique<Mesh>(xyz, faces);
            }), "xyz"_a, "faces"_a)
            .def("closest_faces", &Mesh::closestFaces, "xyz"_a)
            .def("sector", &Mesh::sector, "vertex_halfedges"_a)
            .def("__repr__", &Mesh::toString);
}

PYBIND11_MODULE(compiled, m) {
    m.doc() = "Python bindings for CGAL functionality used for timezone location";

    pybind_mesh(m);
    m.def(
            "triangulate",
            &cgal_tz::triangulate,
            "barycentric_points"_a,
            "segment_indices"_a,
            "Generate a constrained delaunay triangulation. Returns new point coordinates and all triangles"
    );
}

common.h

#ifndef CGAL_TZ_COMMON_H
#define CGAL_TZ_COMMON_H

#include <cstdlib>
#include <iostream>
#include <sstream>
#include <fstream>
#include <memory>
#include <utility>
#include <vector>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Random.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_shortest_path.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>

#include <boost/lexical_cast.hpp>

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
#include <pybind11/numpy.h>
#include <pybind11/operators.h>
#include <pybind11/eigen.h>
#include <pybind11/functional.h>


// CGAL typedefs
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

// Surface_mesh typedefs
typedef CGAL::Surface_mesh<Kernel::Point_3> Triangle_mesh;
typedef boost::graph_traits<Triangle_mesh> Graph_traits;
typedef Graph_traits::vertex_descriptor vertex_descriptor;
typedef Graph_traits::face_descriptor face_descriptor;
typedef Graph_traits::halfedge_descriptor halfedge_descriptor;

// Surface_mesh_shortest_path typedefs
typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Triangle_mesh> TriangleMeshTraits;
typedef CGAL::Surface_mesh_shortest_path<TriangleMeshTraits> Surface_mesh_shortest_path;
typedef TriangleMeshTraits::Barycentric_coordinates Barycentric_coordinates;
typedef std::pair<face_descriptor, Barycentric_coordinates> face_location_descriptor;

// AABB_tree typedefs
typedef CGAL::AABB_face_graph_triangle_primitive<Triangle_mesh> AABBPrimitive;
typedef CGAL::AABB_traits<Kernel, AABBPrimitive> AABBTraits;
typedef CGAL::AABB_tree<AABBTraits> AABB_tree;


// pybind11 typedefs
namespace py = pybind11;
using namespace py::literals;
typedef py::array_t<double, py::array::c_style | py::array::forcecast> np_array_double;
typedef py::array_t<int, py::array::c_style | py::array::forcecast> np_array_int;

#endif

mesh.cpp

#include "arrays.h"
#include "mesh.h"
#include "triangulation.h"

namespace cgal_tz {
    using cgal_tz::Mesh;

    Mesh::Mesh(const np_array_double &xyz, np_array_int faces) {
        py::buffer_info xyz_buffer = xyz.request(), faces_buffer = faces.request();

        // Validate
        if (xyz_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `xyz` must be two");

        if (xyz_buffer.shape[1] != 3)
            throw std::runtime_error("Number of columns in `xyz` must be three");

        if (faces_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `faces` must be two");

        if (faces_buffer.shape[1] != 3)
            throw std::runtime_error("Number of columns in `faces` must be three");

        // Add vertices
        auto xyz_ptr = (double *) xyz_buffer.ptr;
        std::vector<vertex_descriptor> added_vertices{};
        for (size_t row_idx = 0; row_idx < xyz_buffer.shape[0]; ++row_idx) {
            vertex_descriptor v = tmesh_.add_vertex(
                    Kernel::Point_3(xyz_ptr[3 * row_idx],
                                    xyz_ptr[3 * row_idx + 1],
                                    xyz_ptr[3 * row_idx + 2])
            );
            added_vertices.push_back(v);
            ++n_vertices_;
        }

        // Add faces
        int *faces_ptr = (int *) faces_buffer.ptr;
        for (size_t row_idx = 0; row_idx < faces_buffer.shape[0]; ++row_idx) {
            for (int i = 0; i < 3; ++i) {
                if (faces_ptr[3 * row_idx + i] < 0) {
                    throw std::runtime_error("Negative face vertex index present");
                } else if (faces_ptr[3 * row_idx + i] >= n_vertices_) {
                    throw std::runtime_error("Face vertex index >= n_vertices present");
                }
            }
            tmesh_.add_face(added_vertices[faces_ptr[3 * row_idx]],
                            added_vertices[faces_ptr[3 * row_idx + 1]],
                            added_vertices[faces_ptr[3 * row_idx + 2]]);
            ++n_faces_;
        }

        // Build the AABB tree (and speed up queries)
        smsp_ = std::make_unique<Surface_mesh_shortest_path>(tmesh_);
        aabb_ = std::make_unique<AABB_tree>();
        smsp_->build_aabb_tree<AABBTraits>(*aabb_);
        aabb_->accelerate_distance_queries();
    }

    np_array_int Mesh::closestFaces(const np_array_double &xyz) {
        py::buffer_info xyz_buffer = xyz.request();
        // Validate
        if (xyz_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `xyz` must be two");

        if (xyz_buffer.shape[1] != 3)
            throw std::runtime_error("Number of columns in `xyz` must be three");

        // Get vertices
        std::vector<Kernel::Point_3> points{};

        auto xyz_ptr = (double *) xyz_buffer.ptr;
        for (size_t row_idx = 0; row_idx < xyz_buffer.shape[0]; ++row_idx) {
            points.emplace_back(
                    xyz_ptr[3 * row_idx],
                    xyz_ptr[3 * row_idx + 1],
                    xyz_ptr[3 * row_idx + 2]
            );
        }
        auto face_locations = closestFaceLocations(points);

        std::vector<int> output_faces{};
        for (auto const &face_location : face_locations) {
            output_faces.push_back(face_location.first);
        }

        return int_array_1d(output_faces);
    }


    np_array_int Mesh::intersections(np_array_double lonlat) {
        // Validate
        py::buffer_info lonlat_buffer = lonlat.request();
        if (lonlat_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `lonlat` must be two");

        if (lonlat_buffer.shape[1] != 2)
            throw std::runtime_error("Number of columns in `lonlat` must be two");

        Kernel::Vector_3 down{0.0, 0.0, -1.0};
        std::vector<Kernel::Ray_3> rays{};

        rays.reserve(lonlat_buffer.shape[0]);
        auto lonlat_ptr = (double *) lonlat_buffer.ptr;
        for (size_t row_idx = 0; row_idx < lonlat_buffer.shape[0]; ++row_idx) {
            Kernel::Point_3 point{lonlat_ptr[2 * row_idx], lonlat_ptr[2 * row_idx + 1], 1.0};
            rays.emplace_back(point, down);
        }

        std::vector<face_location_descriptor> intersections;
        intersections.reserve(rays.size());
        for (Kernel::Ray_3 ray : rays) {
            aabb_->all_intersected_primitives()
            auto intersection = smsp_->locate(ray, *aabb_);
            intersections.push_back(intersection);
        }

        std::vector<int> output_faces{};
        for (auto const &face_location : face_locations) {
            output_faces.push_back(face_location.first);
        }

        return int_array_1d(output_faces);

        std::vector<CGALMeshPoint> mesh_points{};
        mesh_points.reserve(intersections.size());
        for (face_location_descriptor face_location : intersections) {
            mesh_points.emplace_back(face_location, tmesh_);
        }
        return mesh_points;
    }

    np_array_int Mesh::sector(const np_array_int &vertex_halfedges) {
        py::buffer_info vh_buffer = vertex_halfedges.request();

        // Validate
        if (vh_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `vertex_halfedges` must be two");

        if (vh_buffer.shape[1] != 2)
            throw std::runtime_error("Number of columns in `vertex_halfedges` must be two");

        // Add vertices
        auto vh_ptr = (int *) vh_buffer.ptr;
        std::vector<halfedge_descriptor> halfedges_to_visit{};
        std::unordered_set<face_descriptor> visited_faces{};
        std::unordered_set<halfedge_descriptor> visited_halfedges{};
        std::unordered_set<halfedge_descriptor> border_halfedges{};
        for (size_t row_idx = 0; row_idx < vh_buffer.shape[0]; ++row_idx) {
            int source_index = vh_ptr[2 * row_idx];
            int target_index = vh_ptr[2 * row_idx + 1];
            if (source_index >= n_vertices_ || target_index >= n_vertices_) {
                throw std::runtime_error("An input vertex index was larger than the number of indices");
            }
            vertex_descriptor source(source_index);
            vertex_descriptor target(target_index);
            auto he = tmesh_.halfedge(source, target);
            halfedges_to_visit.push_back(he);
            border_halfedges.insert(he);
        }

        while (!halfedges_to_visit.empty()) {
            auto visited_halfedge = halfedges_to_visit.back();
            halfedges_to_visit.pop_back();
            if (visited_halfedges.find(visited_halfedge) != visited_halfedges.end()) {
                continue;
            }
            auto visited_face = tmesh_.face(visited_halfedge);
            if (visited_face == Triangle_mesh::null_face()) {
                continue;
            }
            visited_faces.insert(visited_face);
            auto face_halfedge = tmesh_.halfedge(visited_face);
            for (int i = 0; i < 3; ++i) {
                face_halfedge = tmesh_.next(face_halfedge);
                if (visited_halfedges.find(face_halfedge) != visited_halfedges.end())
                    continue;
                if (border_halfedges.find(face_halfedge) != border_halfedges.end())
                    continue;
                visited_halfedges.insert(face_halfedge);
                auto opposite_halfedge = tmesh_.opposite(face_halfedge);
                if (opposite_halfedge == Triangle_mesh::null_halfedge())
                    continue;
                if (visited_halfedges.find(opposite_halfedge) == visited_halfedges.end())
                    halfedges_to_visit.push_back(opposite_halfedge);
            }
        }
        std::vector<int> output_faces{};
        std::copy(visited_faces.begin(), visited_faces.end(), std::back_inserter(output_faces));
        return int_array_1d(output_faces);
    }


    std::string Mesh::toString() {
        std::stringstream ss;
        ss << "Mesh(#vertices=" << tmesh_.number_of_vertices()
           << ", #faces=" << tmesh_.number_of_faces() << ")";
        return ss.str();
    }

    std::vector<face_location_descriptor> Mesh::closestFaceLocations(
            std::vector<Kernel::Point_3> &points) {
        std::vector<face_location_descriptor> locations{};
        locations.reserve(points.size());
        for (Kernel::Point_3 point : points) {
            locations.push_back(closestFaceLocation(point));
        }
        return locations;
    }

    face_location_descriptor Mesh::closestFaceLocation(Kernel::Point_3 point) {
        return smsp_->locate(point, *aabb_);
    }
}

mesh.h

#ifndef CGAL_TZ_MESH_H
#define CGAL_TZ_MESH_H

#include "common.h"

namespace cgal_tz {
    class Mesh {
    private:
        Triangle_mesh tmesh_{};
        std::unique_ptr<Surface_mesh_shortest_path> smsp_{nullptr};
        std::unique_ptr<AABB_tree> aabb_{nullptr};
        int n_vertices_{0};
        int n_faces_{0};

    public:
        Mesh(const np_array_double& xyz, np_array_int faces);
        np_array_int closestFaces(const np_array_double& xyz);
        np_array_int sector(const np_array_int& vertex_halfedges);
        std::string toString();

    private:
        std::vector<face_location_descriptor> closestFaceLocations(std::vector<Kernel::Point_3> &points);
        face_location_descriptor closestFaceLocation(Kernel::Point_3 point);
    };
}
#endif

triangulation.cpp

#include "arrays.h"
#include "triangulation.h"

namespace cgal_tz {
    py::tuple triangulate(np_array_double barycentric_points, np_array_int segment_indices) {
        Triangulator triangulator{std::move(barycentric_points), segment_indices};
        return triangulator.getTriangulation();
    }

    Triangulator::Triangulator(np_array_double barycentric_points, np_array_int segment_indices) {
        py::buffer_info points_buffer = barycentric_points.request();
        py::buffer_info segment_indices_buffer = segment_indices.request();

        // Validate
        if (points_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `points` must be 2");

        if (points_buffer.shape[1] != 2)
            throw std::runtime_error("Number of columns in `points` must be 2");

        if (segment_indices_buffer.ndim != 2)
            throw std::runtime_error("Number of dimensions of `segment_indices_buffer` must be 2");

        if (segment_indices_buffer.shape[1] != 2)
            throw std::runtime_error("Number of columns in `segment_indices_buffer` must be 2");

        auto points_ptr = (double *) points_buffer.ptr;
        std::vector<Point> points{};
        points.reserve(points_buffer.shape[0]);
        for (size_t row_idx = 0; row_idx < points_buffer.shape[0]; ++row_idx) {
            Point p(points_ptr[2 * row_idx], points_ptr[2 * row_idx + 1]);
            points.push_back(p);
        }
        n_points_ = points.size();

        auto segment_indices_ptr = (int *) segment_indices_buffer.ptr;
        std::vector<std::pair<int, int>> segments{};
        for (size_t row_idx = 0; row_idx < segment_indices_buffer.shape[0]; ++row_idx) {
            std::pair<int, int> edge(segment_indices_ptr[2 * row_idx], segment_indices_ptr[2 * row_idx + 1]);
            segments.push_back(edge);
        }
        addPoints(points);
        addConstraints(points, segments);
    }

    void Triangulator::addPoints(std::vector<Point> points) {
        // Store vertex indices with the points
        unsigned idx = 0;
        for (auto point : points) {
            cdtp_.insert(point);
            // Update the vertex index of the newly added point; it will be at the end of the finite vertices
            auto vit = cdtp_.finite_vertices_end();
            --vit;
            vit->info() = idx;
            ++idx;
        }
    }

    void Triangulator::addConstraints(std::vector<Point> points, std::vector<std::pair<int, int>> segments) {
        for (std::pair<int, int> edge : segments) {
            if (edge.first >= points.size() or edge.second >= points.size()) {
                throw std::runtime_error("required_edges had an index that was too large");
            }
            cdtp_.insert_constraint(points[edge.first], points[edge.second]);
        }

        // Set the vertex indices for any vertices created due to constraints
        unsigned idx = 0;
        for (auto vit = cdtp_.finite_vertices_begin(); vit != cdtp_.finite_vertices_end(); ++vit) {
            if (idx >= n_points_) {
                vit->info() = idx;  // Set the (currently uninitialized) vertex index appropriately
            }
            ++idx;
        }
    }

    py::tuple Triangulator::getTriangulation() {
        auto new_point_coordinates = getNewPointCoordinates();
        auto triangles = getTriangles();
        auto segments = getSegments();

        return py::make_tuple(
                double_array_2d(new_point_coordinates, 2),
                int_array_2d(triangles, 3),
                segments
        );
    }

    std::vector<double> Triangulator::getNewPointCoordinates() {
        std::vector<double> new_point_coordinates{};

        // Set the vertex indices on new vertices
        unsigned idx = 0;
        for (auto vit = cdtp_.finite_vertices_begin(); vit != cdtp_.finite_vertices_end(); ++vit) {
            if (idx >= n_points_) {
                new_point_coordinates.push_back(vit->point().x());
                new_point_coordinates.push_back(vit->point().y());
            } else {
                ++idx;
            }
        }
        return new_point_coordinates;
    };

    std::vector<int> Triangulator::getTriangles() {
        std::vector<int> triangles{};
        for (auto fit = cdtp_.finite_faces_begin(); fit != cdtp_.finite_faces_end(); ++fit) {
            for (int i = 0; i < 3; ++i) {
                triangles.push_back(fit->vertex(i)->info());
            }
        }
        return triangles;
    }

    bool compareSegment(std::vector<int> i1, std::vector<int> i2)
    {
        if (i1[0] == i2[0]) {
            return i1.back() < i2.back();
        } else {
            return i1[0] < i2[0];
        }
    }

    std::vector<std::vector<int>> Triangulator::getSegments() {
        std::vector<std::vector<int>> segments{};
        for (auto cit = cdtp_.constraints_begin(); cit != cdtp_.constraints_end(); ++cit) {
            std::vector<int> segment{};
            for (CDTP::Vertices_in_constraint_iterator vicit = cdtp_.vertices_in_constraint_begin(*cit);
                 vicit != cdtp_.vertices_in_constraint_end(*cit); ++vicit) {
                segment.push_back((*vicit)->info());
            }
            segments.push_back(segment);
        }
        sort(segments.begin(), segments.end(), compareSegment);
        return segments;
    }
}

triangulation.h

#ifndef CGAL_TZ_TRIANGULATION_H
#define CGAL_TZ_TRIANGULATION_H

#include "common.h"
#include <CGAL/intersections.h>
#include <CGAL/Constrained_triangulation_face_base_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>


namespace cgal_tz {
    typedef CGAL::Exact_predicates_tag Itag;
    typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned, Kernel> Vb;
    typedef CGAL::Constrained_triangulation_face_base_2<Kernel> Fb;
    typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
    typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel, TDS, Itag> CDT;
    typedef CGAL::Constrained_triangulation_plus_2<CDT> CDTP;
    typedef CDTP::Point Point;

    py::tuple triangulate(np_array_double barycentric_points, np_array_int segment_indices);

    class Triangulator {
    private:
        CDTP cdtp_{};
        unsigned long n_points_{0};

    public:
        Triangulator(np_array_double barycentric_points, np_array_int segment_indices);
        py::tuple getTriangulation();

    private:
        void addPoints(std::vector<Point> points);
        void addConstraints(std::vector<Point> points, std::vector<std::pair<int, int>> segments);
        std::vector<double> getNewPointCoordinates();
        std::vector<int> getTriangles();
        std::vector<std::vector<int>> getSegments();
    };
}

#endif