Remove deps on pandas

<https://github.com/ydb-platform/ydb/pull/14418>

<https://github.com/ydb-platform/ydb/pull/14419>
\-- аналогичные правки в gh

Хочу залить в обход синка, чтобы посмотреть удалится ли pandas в нашей gh репе через piglet
commit_hash:abca127aa37d4dbb94b07e1e18cdb8eb5b711860

3 files changed, 0 insertions, 2931 deletions

diff --git a/contrib/python/matplotlib/py3/src/tri/_tri.cpp b/contrib/python/matplotlib/py3/src/tri/_tri.cpp
deleted file mode 100644
index 2674a3140b..0000000000
--- a/contrib/python/matplotlib/py3/src/tri/_tri.cpp
+++ /dev/null
@@ -1,2074 +0,0 @@
-/* This file contains liberal use of asserts to assist code development and
- * debugging.  Standard matplotlib builds disable asserts so they cause no
- * performance reduction.  To enable the asserts, you need to undefine the
- * NDEBUG macro, which is achieved by adding the following
- *     undef_macros=['NDEBUG']
- * to the appropriate make_extension call in setupext.py, and then rebuilding.
- */
-#include "../mplutils.h"
-#include "_tri.h"
-
-#include <algorithm>
-#include <random>
-#include <set>
-
-
-TriEdge::TriEdge()
-    : tri(-1), edge(-1)
-{}
-
-TriEdge::TriEdge(int tri_, int edge_)
-    : tri(tri_), edge(edge_)
-{}
-
-bool TriEdge::operator<(const TriEdge& other) const
-{
-    if (tri != other.tri)
-        return tri < other.tri;
-    else
-        return edge < other.edge;
-}
-
-bool TriEdge::operator==(const TriEdge& other) const
-{
-    return tri == other.tri && edge == other.edge;
-}
-
-bool TriEdge::operator!=(const TriEdge& other) const
-{
-    return !operator==(other);
-}
-
-std::ostream& operator<<(std::ostream& os, const TriEdge& tri_edge)
-{
-    return os << tri_edge.tri << ' ' << tri_edge.edge;
-}
-
-
-
-XY::XY()
-{}
-
-XY::XY(const double& x_, const double& y_)
-    : x(x_), y(y_)
-{}
-
-double XY::angle() const
-{
-    return atan2(y, x);
-}
-
-double XY::cross_z(const XY& other) const
-{
-    return x*other.y - y*other.x;
-}
-
-bool XY::is_right_of(const XY& other) const
-{
-    if (x == other.x)
-        return y > other.y;
-    else
-        return x > other.x;
-}
-
-bool XY::operator==(const XY& other) const
-{
-    return x == other.x && y == other.y;
-}
-
-bool XY::operator!=(const XY& other) const
-{
-    return x != other.x || y != other.y;
-}
-
-XY XY::operator*(const double& multiplier) const
-{
-    return XY(x*multiplier, y*multiplier);
-}
-
-const XY& XY::operator+=(const XY& other)
-{
-    x += other.x;
-    y += other.y;
-    return *this;
-}
-
-const XY& XY::operator-=(const XY& other)
-{
-    x -= other.x;
-    y -= other.y;
-    return *this;
-}
-
-XY XY::operator+(const XY& other) const
-{
-    return XY(x + other.x, y + other.y);
-}
-
-XY XY::operator-(const XY& other) const
-{
-    return XY(x - other.x, y - other.y);
-}
-
-std::ostream& operator<<(std::ostream& os, const XY& xy)
-{
-    return os << '(' << xy.x << ' ' << xy.y << ')';
-}
-
-
-
-XYZ::XYZ(const double& x_, const double& y_, const double& z_)
-    : x(x_), y(y_), z(z_)
-{}
-
-XYZ XYZ::cross(const XYZ& other) const
-{
-    return XYZ(y*other.z - z*other.y,
-               z*other.x - x*other.z,
-               x*other.y - y*other.x);
-}
-
-double XYZ::dot(const XYZ& other) const
-{
-    return x*other.x + y*other.y + z*other.z;
-}
-
-XYZ XYZ::operator-(const XYZ& other) const
-{
-    return XYZ(x - other.x, y - other.y, z - other.z);
-}
-
-std::ostream& operator<<(std::ostream& os, const XYZ& xyz)
-{
-    return os << '(' << xyz.x << ' ' << xyz.y << ' ' << xyz.z << ')';
-}
-
-
-
-BoundingBox::BoundingBox()
-    : empty(true), lower(0.0, 0.0), upper(0.0, 0.0)
-{}
-
-void BoundingBox::add(const XY& point)
-{
-    if (empty) {
-        empty = false;
-        lower = upper = point;
-    } else {
-        if      (point.x < lower.x) lower.x = point.x;
-        else if (point.x > upper.x) upper.x = point.x;
-
-        if      (point.y < lower.y) lower.y = point.y;
-        else if (point.y > upper.y) upper.y = point.y;
-    }
-}
-
-void BoundingBox::expand(const XY& delta)
-{
-    if (!empty) {
-        lower -= delta;
-        upper += delta;
-    }
-}
-
-
-
-ContourLine::ContourLine()
-    : std::vector<XY>()
-{}
-
-void ContourLine::push_back(const XY& point)
-{
-    if (empty() || point != back())
-        std::vector<XY>::push_back(point);
-}
-
-void ContourLine::write() const
-{
-    std::cout << "ContourLine of " << size() << " points:";
-    for (const_iterator it = begin(); it != end(); ++it)
-        std::cout << ' ' << *it;
-    std::cout << std::endl;
-}
-
-
-
-void write_contour(const Contour& contour)
-{
-    std::cout << "Contour of " << contour.size() << " lines." << std::endl;
-    for (Contour::const_iterator it = contour.begin(); it != contour.end(); ++it)
-        it->write();
-}
-
-
-
-Triangulation::Triangulation(const CoordinateArray& x,
-                             const CoordinateArray& y,
-                             const TriangleArray& triangles,
-                             const MaskArray& mask,
-                             const EdgeArray& edges,
-                             const NeighborArray& neighbors,
-                             bool correct_triangle_orientations)
-    : _x(x),
-      _y(y),
-      _triangles(triangles),
-      _mask(mask),
-      _edges(edges),
-      _neighbors(neighbors)
-{
-    if (_x.ndim() != 1 || _y.ndim() != 1 || _x.shape(0) != _y.shape(0))
-        throw std::invalid_argument("x and y must be 1D arrays of the same length");
-
-    if (_triangles.ndim() != 2 || _triangles.shape(1) != 3)
-        throw std::invalid_argument("triangles must be a 2D array of shape (?,3)");
-
-    // Optional mask.
-    if (_mask.size() > 0 &&
-        (_mask.ndim() != 1 || _mask.shape(0) != _triangles.shape(0)))
-        throw std::invalid_argument(
-            "mask must be a 1D array with the same length as the triangles array");
-
-    // Optional edges.
-    if (_edges.size() > 0 &&
-        (_edges.ndim() != 2 || _edges.shape(1) != 2))
-        throw std::invalid_argument("edges must be a 2D array with shape (?,2)");
-
-    // Optional neighbors.
-    if (_neighbors.size() > 0 &&
-        (_neighbors.ndim() != 2 || _neighbors.shape() != _triangles.shape()))
-        throw std::invalid_argument(
-            "neighbors must be a 2D array with the same shape as the triangles array");
-
-    if (correct_triangle_orientations)
-        correct_triangles();
-}
-
-void Triangulation::calculate_boundaries()
-{
-    get_neighbors();  // Ensure _neighbors has been created.
-
-    // Create set of all boundary TriEdges, which are those which do not
-    // have a neighbor triangle.
-    typedef std::set<TriEdge> BoundaryEdges;
-    BoundaryEdges boundary_edges;
-    for (int tri = 0; tri < get_ntri(); ++tri) {
-        if (!is_masked(tri)) {
-            for (int edge = 0; edge < 3; ++edge) {
-                if (get_neighbor(tri, edge) == -1) {
-                    boundary_edges.insert(TriEdge(tri, edge));
-                }
-            }
-        }
-    }
-
-    // Take any boundary edge and follow the boundary until return to start
-    // point, removing edges from boundary_edges as they are used.  At the same
-    // time, initialise the _tri_edge_to_boundary_map.
-    while (!boundary_edges.empty()) {
-        // Start of new boundary.
-        BoundaryEdges::iterator it = boundary_edges.begin();
-        int tri = it->tri;
-        int edge = it->edge;
-        _boundaries.push_back(Boundary());
-        Boundary& boundary = _boundaries.back();
-
-        while (true) {
-            boundary.push_back(TriEdge(tri, edge));
-            boundary_edges.erase(it);
-            _tri_edge_to_boundary_map[TriEdge(tri, edge)] =
-                BoundaryEdge(_boundaries.size()-1, boundary.size()-1);
-
-            // Move to next edge of current triangle.
-            edge = (edge+1) % 3;
-
-            // Find start point index of boundary edge.
-            int point = get_triangle_point(tri, edge);
-
-            // Find next TriEdge by traversing neighbors until find one
-            // without a neighbor.
-            while (get_neighbor(tri, edge) != -1) {
-                tri = get_neighbor(tri, edge);
-                edge = get_edge_in_triangle(tri, point);
-            }
-
-            if (TriEdge(tri,edge) == boundary.front())
-                break;  // Reached beginning of this boundary, so finished it.
-            else
-                it = boundary_edges.find(TriEdge(tri, edge));
-        }
-    }
-}
-
-void Triangulation::calculate_edges()
-{
-    assert(!has_edges() && "Expected empty edges array");
-
-    // Create set of all edges, storing them with start point index less than
-    // end point index.
-    typedef std::set<Edge> EdgeSet;
-    EdgeSet edge_set;
-    for (int tri = 0; tri < get_ntri(); ++tri) {
-        if (!is_masked(tri)) {
-            for (int edge = 0; edge < 3; edge++) {
-                int start = get_triangle_point(tri, edge);
-                int end   = get_triangle_point(tri, (edge+1)%3);
-                edge_set.insert(start > end ? Edge(start,end) : Edge(end,start));
-            }
-        }
-    }
-
-    // Convert to python _edges array.
-    py::ssize_t dims[2] = {static_cast<py::ssize_t>(edge_set.size()), 2};
-    _edges = EdgeArray(dims);
-    auto edges = _edges.mutable_data();
-
-    int i = 0;
-    for (EdgeSet::const_iterator it = edge_set.begin(); it != edge_set.end(); ++it) {
-        edges[i++] = it->start;
-        edges[i++] = it->end;
-    }
-}
-
-void Triangulation::calculate_neighbors()
-{
-    assert(!has_neighbors() && "Expected empty neighbors array");
-
-    // Create _neighbors array with shape (ntri,3) and initialise all to -1.
-    py::ssize_t dims[2] = {get_ntri(), 3};
-    _neighbors = NeighborArray(dims);
-    auto* neighbors = _neighbors.mutable_data();
-
-    int tri, edge;
-    std::fill(neighbors, neighbors+3*get_ntri(), -1);
-
-    // For each triangle edge (start to end point), find corresponding neighbor
-    // edge from end to start point.  Do this by traversing all edges and
-    // storing them in a map from edge to TriEdge.  If corresponding neighbor
-    // edge is already in the map, don't need to store new edge as neighbor
-    // already found.
-    typedef std::map<Edge, TriEdge> EdgeToTriEdgeMap;
-    EdgeToTriEdgeMap edge_to_tri_edge_map;
-    for (tri = 0; tri < get_ntri(); ++tri) {
-        if (!is_masked(tri)) {
-            for (edge = 0; edge < 3; ++edge) {
-                int start = get_triangle_point(tri, edge);
-                int end   = get_triangle_point(tri, (edge+1)%3);
-                EdgeToTriEdgeMap::iterator it =
-                    edge_to_tri_edge_map.find(Edge(end,start));
-                if (it == edge_to_tri_edge_map.end()) {
-                    // No neighbor edge exists in the edge_to_tri_edge_map, so
-                    // add this edge to it.
-                    edge_to_tri_edge_map[Edge(start,end)] = TriEdge(tri,edge);
-                } else {
-                    // Neighbor edge found, set the two elements of _neighbors
-                    // and remove edge from edge_to_tri_edge_map.
-                    neighbors[3*tri + edge] = it->second.tri;
-                    neighbors[3*it->second.tri + it->second.edge] = tri;
-                    edge_to_tri_edge_map.erase(it);
-                }
-            }
-        }
-    }
-
-    // Note that remaining edges in the edge_to_tri_edge_map correspond to
-    // boundary edges, but the boundaries are calculated separately elsewhere.
-}
-
-Triangulation::TwoCoordinateArray Triangulation::calculate_plane_coefficients(
-    const CoordinateArray& z)
-{
-    if (z.ndim() != 1 || z.shape(0) != _x.shape(0))
-        throw std::invalid_argument(
-            "z must be a 1D array with the same length as the triangulation x and y arrays");
-
-    int dims[2] = {get_ntri(), 3};
-    Triangulation::TwoCoordinateArray planes_array(dims);
-    auto planes = planes_array.mutable_unchecked<2>();
-    auto triangles = _triangles.unchecked<2>();
-    auto x = _x.unchecked<1>();
-    auto y = _y.unchecked<1>();
-    auto z_ptr = z.unchecked<1>();
-
-    int point;
-    for (int tri = 0; tri < get_ntri(); ++tri) {
-        if (is_masked(tri)) {
-            planes(tri, 0) = 0.0;
-            planes(tri, 1) = 0.0;
-            planes(tri, 2) = 0.0;
-        }
-        else {
-            // Equation of plane for all points r on plane is r.normal = p
-            // where normal is vector normal to the plane, and p is a
-            // constant.  Rewrite as
-            // r_x*normal_x + r_y*normal_y + r_z*normal_z = p
-            // and rearrange to give
-            // r_z = (-normal_x/normal_z)*r_x + (-normal_y/normal_z)*r_y +
-            //       p/normal_z
-            point = triangles(tri, 0);
-            XYZ point0(x(point), y(point), z_ptr(point));
-            point = triangles(tri, 1);
-            XYZ side01 = XYZ(x(point), y(point), z_ptr(point)) - point0;
-            point = triangles(tri, 2);
-            XYZ side02 = XYZ(x(point), y(point), z_ptr(point)) - point0;
-
-            XYZ normal = side01.cross(side02);
-
-            if (normal.z == 0.0) {
-                // Normal is in x-y plane which means triangle consists of
-                // colinear points. To avoid dividing by zero, we use the
-                // Moore-Penrose pseudo-inverse.
-                double sum2 = (side01.x*side01.x + side01.y*side01.y +
-                               side02.x*side02.x + side02.y*side02.y);
-                double a = (side01.x*side01.z + side02.x*side02.z) / sum2;
-                double b = (side01.y*side01.z + side02.y*side02.z) / sum2;
-                planes(tri, 0) = a;
-                planes(tri, 1) = b;
-                planes(tri, 2) = point0.z - a*point0.x - b*point0.y;
-            }
-            else {
-                planes(tri, 0) = -normal.x / normal.z;           // x
-                planes(tri, 1) = -normal.y / normal.z;           // y
-                planes(tri, 2) = normal.dot(point0) / normal.z;  // constant
-            }
-        }
-    }
-
-    return planes_array;
-}
-
-void Triangulation::correct_triangles()
-{
-    auto triangles = _triangles.mutable_data();
-    auto neighbors = _neighbors.mutable_data();
-
-    for (int tri = 0; tri < get_ntri(); ++tri) {
-        XY point0 = get_point_coords(triangles[3*tri]);
-        XY point1 = get_point_coords(triangles[3*tri+1]);
-        XY point2 = get_point_coords(triangles[3*tri+2]);
-        if ( (point1 - point0).cross_z(point2 - point0) < 0.0) {
-            // Triangle points are clockwise, so change them to anticlockwise.
-            std::swap(triangles[3*tri+1], triangles[3*tri+2]);
-            if (has_neighbors())
-                std::swap(neighbors[3*tri+1], neighbors[3*tri+2]);
-        }
-    }
-}
-
-const Triangulation::Boundaries& Triangulation::get_boundaries() const
-{
-    if (_boundaries.empty())
-        const_cast<Triangulation*>(this)->calculate_boundaries();
-    return _boundaries;
-}
-
-void Triangulation::get_boundary_edge(const TriEdge& triEdge,
-                                      int& boundary,
-                                      int& edge) const
-{
-    get_boundaries();  // Ensure _tri_edge_to_boundary_map has been created.
-    TriEdgeToBoundaryMap::const_iterator it =
-        _tri_edge_to_boundary_map.find(triEdge);
-    assert(it != _tri_edge_to_boundary_map.end() &&
-           "TriEdge is not on a boundary");
-    boundary = it->second.boundary;
-    edge = it->second.edge;
-}
-
-int Triangulation::get_edge_in_triangle(int tri, int point) const
-{
-    assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");
-    assert(point >= 0 && point < get_npoints() && "Point index out of bounds.");
-
-    auto triangles = _triangles.data();
-
-    for (int edge = 0; edge < 3; ++edge) {
-        if (triangles[3*tri + edge] == point)
-            return edge;
-    }
-    return -1;  // point is not in triangle.
-}
-
-Triangulation::EdgeArray& Triangulation::get_edges()
-{
-    if (!has_edges())
-        calculate_edges();
-    return _edges;
-}
-
-int Triangulation::get_neighbor(int tri, int edge) const
-{
-    assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");
-    assert(edge >= 0 && edge < 3 && "Edge index out of bounds");
-    if (!has_neighbors())
-        const_cast<Triangulation&>(*this).calculate_neighbors();
-    return _neighbors.data()[3*tri + edge];
-}
-
-TriEdge Triangulation::get_neighbor_edge(int tri, int edge) const
-{
-    int neighbor_tri = get_neighbor(tri, edge);
-    if (neighbor_tri == -1)
-        return TriEdge(-1,-1);
-    else
-        return TriEdge(neighbor_tri,
-                       get_edge_in_triangle(neighbor_tri,
-                                            get_triangle_point(tri,
-                                                               (edge+1)%3)));
-}
-
-Triangulation::NeighborArray& Triangulation::get_neighbors()
-{
-    if (!has_neighbors())
-        calculate_neighbors();
-    return _neighbors;
-}
-
-int Triangulation::get_npoints() const
-{
-    return _x.shape(0);
-}
-
-int Triangulation::get_ntri() const
-{
-    return _triangles.shape(0);
-}
-
-XY Triangulation::get_point_coords(int point) const
-{
-    assert(point >= 0 && point < get_npoints() && "Point index out of bounds.");
-    return XY(_x.data()[point], _y.data()[point]);
-}
-
-int Triangulation::get_triangle_point(int tri, int edge) const
-{
-    assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");
-    assert(edge >= 0 && edge < 3 && "Edge index out of bounds");
-    return _triangles.data()[3*tri + edge];
-}
-
-int Triangulation::get_triangle_point(const TriEdge& tri_edge) const
-{
-    return get_triangle_point(tri_edge.tri, tri_edge.edge);
-}
-
-bool Triangulation::has_edges() const
-{
-    return _edges.size() > 0;
-}
-
-bool Triangulation::has_mask() const
-{
-    return _mask.size() > 0;
-}
-
-bool Triangulation::has_neighbors() const
-{
-    return _neighbors.size() > 0;
-}
-
-bool Triangulation::is_masked(int tri) const
-{
-    assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds.");
-    return has_mask() && _mask.data()[tri];
-}
-
-void Triangulation::set_mask(const MaskArray& mask)
-{
-    if (mask.size() > 0 &&
-        (mask.ndim() != 1 || mask.shape(0) != _triangles.shape(0)))
-        throw std::invalid_argument(
-            "mask must be a 1D array with the same length as the triangles array");
-
-    _mask = mask;
-
-    // Clear derived fields so they are recalculated when needed.
-    _edges = EdgeArray();
-    _neighbors = NeighborArray();
-    _boundaries.clear();
-}
-
-void Triangulation::write_boundaries() const
-{
-    const Boundaries& bs = get_boundaries();
-    std::cout << "Number of boundaries: " << bs.size() << std::endl;
-    for (Boundaries::const_iterator it = bs.begin(); it != bs.end(); ++it) {
-        const Boundary& b = *it;
-        std::cout << "  Boundary of " << b.size() << " points: ";
-        for (Boundary::const_iterator itb = b.begin(); itb != b.end(); ++itb) {
-            std::cout << *itb << ", ";
-        }
-        std::cout << std::endl;
-    }
-}
-
-
-
-TriContourGenerator::TriContourGenerator(Triangulation& triangulation,
-                                         const CoordinateArray& z)
-    : _triangulation(triangulation),
-      _z(z),
-      _interior_visited(2*_triangulation.get_ntri()),
-      _boundaries_visited(0),
-      _boundaries_used(0)
-{
-    if (_z.ndim() != 1 || _z.shape(0) != _triangulation.get_npoints())
-        throw std::invalid_argument(
-            "z must be a 1D array with the same length as the x and y arrays");
-}
-
-void TriContourGenerator::clear_visited_flags(bool include_boundaries)
-{
-    // Clear _interiorVisited.
-    std::fill(_interior_visited.begin(), _interior_visited.end(), false);
-
-    if (include_boundaries) {
-        if (_boundaries_visited.empty()) {
-            const Boundaries& boundaries = get_boundaries();
-
-            // Initialise _boundaries_visited.
-            _boundaries_visited.reserve(boundaries.size());
-            for (Boundaries::const_iterator it = boundaries.begin();
-                    it != boundaries.end(); ++it)
-                _boundaries_visited.push_back(BoundaryVisited(it->size()));
-
-            // Initialise _boundaries_used.
-            _boundaries_used = BoundariesUsed(boundaries.size());
-        }
-
-        // Clear _boundaries_visited.
-        for (BoundariesVisited::iterator it = _boundaries_visited.begin();
-                it != _boundaries_visited.end(); ++it)
-            std::fill(it->begin(), it->end(), false);
-
-        // Clear _boundaries_used.
-        std::fill(_boundaries_used.begin(), _boundaries_used.end(), false);
-    }
-}
-
-py::tuple TriContourGenerator::contour_line_to_segs_and_kinds(const Contour& contour)
-{
-    // Convert all of the lines generated by a call to create_contour() into
-    // their Python equivalents for return to the calling function.
-    // A line is either a closed line loop (in which case the last point is
-    // identical to the first) or an open line strip.  Two NumPy arrays are
-    // created for each line:
-    //   vertices is a double array of shape (npoints, 2) containing the (x, y)
-    //     coordinates of the points in the line
-    //   codes is a uint8 array of shape (npoints,) containing the 'kind codes'
-    //     which are defined in the Path class
-    // and they are appended to the Python lists vertices_list and codes_list
-    // respectively for return to the Python calling function.
-
-    py::list vertices_list(contour.size());
-    py::list codes_list(contour.size());
-
-    for (Contour::size_type i = 0; i < contour.size(); ++i) {
-        const ContourLine& contour_line = contour[i];
-        py::ssize_t npoints = static_cast<py::ssize_t>(contour_line.size());
-
-        py::ssize_t segs_dims[2] = {npoints, 2};
-        CoordinateArray segs(segs_dims);
-        double* segs_ptr = segs.mutable_data();
-
-        py::ssize_t codes_dims[1] = {npoints};
-        CodeArray codes(codes_dims);
-        unsigned char* codes_ptr = codes.mutable_data();
-
-        for (ContourLine::const_iterator it = contour_line.begin();
-             it != contour_line.end(); ++it) {
-            *segs_ptr++ = it->x;
-            *segs_ptr++ = it->y;
-            *codes_ptr++ = (it == contour_line.begin() ? MOVETO : LINETO);
-        }
-
-        // Closed line loop has identical first and last (x, y) points.
-        if (contour_line.size() > 1 &&
-            contour_line.front() == contour_line.back())
-            *(codes_ptr-1) = CLOSEPOLY;
-
-        vertices_list[i] = segs;
-        codes_list[i] = codes;
-    }
-
-    return py::make_tuple(vertices_list, codes_list);
-}
-
-py::tuple TriContourGenerator::contour_to_segs_and_kinds(const Contour& contour)
-{
-    // Convert all of the polygons generated by a call to
-    // create_filled_contour() into their Python equivalents for return to the
-    // calling function.  All of the polygons' points and kinds codes are
-    // combined into single NumPy arrays for each; this avoids having
-    // to determine which polygons are holes as this will be determined by the
-    // renderer.  If there are ntotal points in all of the polygons, the two
-    // NumPy arrays created are:
-    //   vertices is a double array of shape (ntotal, 2) containing the (x, y)
-    //     coordinates of the points in the polygons
-    //   codes is a uint8 array of shape (ntotal,) containing the 'kind codes'
-    //     which are defined in the Path class
-    // and they are returned in the Python lists vertices_list and codes_list
-    // respectively.
-
-    Contour::const_iterator line;
-    ContourLine::const_iterator point;
-
-    // Find total number of points in all contour lines.
-    py::ssize_t n_points = 0;
-    for (line = contour.begin(); line != contour.end(); ++line)
-        n_points += static_cast<py::ssize_t>(line->size());
-
-    // Create segs array for point coordinates.
-    py::ssize_t segs_dims[2] = {n_points, 2};
-    TwoCoordinateArray segs(segs_dims);
-    double* segs_ptr = segs.mutable_data();
-
-    // Create kinds array for code types.
-    py::ssize_t codes_dims[1] = {n_points};
-    CodeArray codes(codes_dims);
-    unsigned char* codes_ptr = codes.mutable_data();
-
-    for (line = contour.begin(); line != contour.end(); ++line) {
-        for (point = line->begin(); point != line->end(); point++) {
-            *segs_ptr++ = point->x;
-            *segs_ptr++ = point->y;
-            *codes_ptr++ = (point == line->begin() ? MOVETO : LINETO);
-        }
-
-        if (line->size() > 1)
-            *(codes_ptr-1) = CLOSEPOLY;
-    }
-
-    py::list vertices_list(1);
-    vertices_list[0] = segs;
-
-    py::list codes_list(1);
-    codes_list[0] = codes;
-
-    return py::make_tuple(vertices_list, codes_list);
-}
-
-py::tuple TriContourGenerator::create_contour(const double& level)
-{
-    clear_visited_flags(false);
-    Contour contour;
-
-    find_boundary_lines(contour, level);
-    find_interior_lines(contour, level, false, false);
-
-    return contour_line_to_segs_and_kinds(contour);
-}
-
-py::tuple TriContourGenerator::create_filled_contour(const double& lower_level,
-                                                     const double& upper_level)
-{
-    if (lower_level >= upper_level)
-        throw std::invalid_argument("filled contour levels must be increasing");
-
-    clear_visited_flags(true);
-    Contour contour;
-
-    find_boundary_lines_filled(contour, lower_level, upper_level);
-    find_interior_lines(contour, lower_level, false, true);
-    find_interior_lines(contour, upper_level, true,  true);
-
-    return contour_to_segs_and_kinds(contour);
-}
-
-XY TriContourGenerator::edge_interp(int tri, int edge, const double& level)
-{
-    return interp(_triangulation.get_triangle_point(tri, edge),
-                  _triangulation.get_triangle_point(tri, (edge+1)%3),
-                  level);
-}
-
-void TriContourGenerator::find_boundary_lines(Contour& contour,
-                                              const double& level)
-{
-    // Traverse boundaries to find starting points for all contour lines that
-    // intersect the boundaries.  For each starting point found, follow the
-    // line to its end before continuing.
-    const Triangulation& triang = _triangulation;
-    const Boundaries& boundaries = get_boundaries();
-    for (Boundaries::const_iterator it = boundaries.begin();
-            it != boundaries.end(); ++it) {
-        const Boundary& boundary = *it;
-        bool startAbove, endAbove = false;
-        for (Boundary::const_iterator itb = boundary.begin();
-                itb != boundary.end(); ++itb) {
-            if (itb == boundary.begin())
-                startAbove = get_z(triang.get_triangle_point(*itb)) >= level;
-            else
-                startAbove = endAbove;
-            endAbove = get_z(triang.get_triangle_point(itb->tri,
-                                                       (itb->edge+1)%3)) >= level;
-            if (startAbove && !endAbove) {
-                // This boundary edge is the start point for a contour line,
-                // so follow the line.
-                contour.push_back(ContourLine());
-                ContourLine& contour_line = contour.back();
-                TriEdge tri_edge = *itb;
-                follow_interior(contour_line, tri_edge, true, level, false);
-            }
-        }
-    }
-}
-
-void TriContourGenerator::find_boundary_lines_filled(Contour& contour,
-                                                     const double& lower_level,
-                                                     const double& upper_level)
-{
-    // Traverse boundaries to find starting points for all contour lines that
-    // intersect the boundaries.  For each starting point found, follow the
-    // line to its end before continuing.
-    const Triangulation& triang = _triangulation;
-    const Boundaries& boundaries = get_boundaries();
-    for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
-        const Boundary& boundary = boundaries[i];
-        for (Boundary::size_type j = 0; j < boundary.size(); ++j) {
-            if (!_boundaries_visited[i][j]) {
-                // z values of start and end of this boundary edge.
-                double z_start = get_z(triang.get_triangle_point(boundary[j]));
-                double z_end = get_z(triang.get_triangle_point(
-                                   boundary[j].tri, (boundary[j].edge+1)%3));
-
-                // Does this boundary edge's z increase through upper level
-                // and/or decrease through lower level?
-                bool incr_upper = (z_start < upper_level && z_end >= upper_level);
-                bool decr_lower = (z_start >= lower_level && z_end < lower_level);
-
-                if (decr_lower || incr_upper) {
-                    // Start point for contour line, so follow it.
-                    contour.push_back(ContourLine());
-                    ContourLine& contour_line = contour.back();
-                    TriEdge start_tri_edge = boundary[j];
-                    TriEdge tri_edge = start_tri_edge;
-
-                    // Traverse interior and boundaries until return to start.
-                    bool on_upper = incr_upper;
-                    do {
-                        follow_interior(contour_line, tri_edge, true,
-                            on_upper ? upper_level : lower_level, on_upper);
-                        on_upper = follow_boundary(contour_line, tri_edge,
-                                       lower_level, upper_level, on_upper);
-                    } while (tri_edge != start_tri_edge);
-
-                    // Close polygon.
-                    contour_line.push_back(contour_line.front());
-                }
-            }
-        }
-    }
-
-    // Add full boundaries that lie between the lower and upper levels.  These
-    // are boundaries that have not been touched by an internal contour line
-    // which are stored in _boundaries_used.
-    for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
-        if (!_boundaries_used[i]) {
-            const Boundary& boundary = boundaries[i];
-            double z = get_z(triang.get_triangle_point(boundary[0]));
-            if (z >= lower_level && z < upper_level) {
-                contour.push_back(ContourLine());
-                ContourLine& contour_line = contour.back();
-                for (Boundary::size_type j = 0; j < boundary.size(); ++j)
-                    contour_line.push_back(triang.get_point_coords(
-                                      triang.get_triangle_point(boundary[j])));
-
-                // Close polygon.
-                contour_line.push_back(contour_line.front());
-            }
-        }
-    }
-}
-
-void TriContourGenerator::find_interior_lines(Contour& contour,
-                                              const double& level,
-                                              bool on_upper,
-                                              bool filled)
-{
-    const Triangulation& triang = _triangulation;
-    int ntri = triang.get_ntri();
-    for (int tri = 0; tri < ntri; ++tri) {
-        int visited_index = (on_upper ? tri+ntri : tri);
-
-        if (_interior_visited[visited_index] || triang.is_masked(tri))
-            continue;  // Triangle has already been visited or is masked.
-
-        _interior_visited[visited_index] = true;
-
-        // Determine edge via which to leave this triangle.
-        int edge = get_exit_edge(tri, level, on_upper);
-        assert(edge >= -1 && edge < 3 && "Invalid exit edge");
-        if (edge == -1)
-            continue;  // Contour does not pass through this triangle.
-
-        // Found start of new contour line loop.
-        contour.push_back(ContourLine());
-        ContourLine& contour_line = contour.back();
-        TriEdge tri_edge = triang.get_neighbor_edge(tri, edge);
-        follow_interior(contour_line, tri_edge, false, level, on_upper);
-
-        // Close line loop
-        contour_line.push_back(contour_line.front());
-    }
-}
-
-bool TriContourGenerator::follow_boundary(ContourLine& contour_line,
-                                          TriEdge& tri_edge,
-                                          const double& lower_level,
-                                          const double& upper_level,
-                                          bool on_upper)
-{
-    const Triangulation& triang = _triangulation;
-    const Boundaries& boundaries = get_boundaries();
-
-    // Have TriEdge to start at, need equivalent boundary edge.
-    int boundary, edge;
-    triang.get_boundary_edge(tri_edge, boundary, edge);
-    _boundaries_used[boundary] = true;
-
-    bool stop = false;
-    bool first_edge = true;
-    double z_start, z_end = 0;
-    while (!stop)
-    {
-        assert(!_boundaries_visited[boundary][edge] && "Boundary already visited");
-        _boundaries_visited[boundary][edge] = true;
-
-        // z values of start and end points of boundary edge.
-        if (first_edge)
-            z_start = get_z(triang.get_triangle_point(tri_edge));
-        else
-            z_start = z_end;
-        z_end = get_z(triang.get_triangle_point(tri_edge.tri,
-                                                (tri_edge.edge+1)%3));
-
-        if (z_end > z_start) {  // z increasing.
-            if (!(!on_upper && first_edge) &&
-                z_end >= lower_level && z_start < lower_level) {
-                stop = true;
-                on_upper = false;
-            } else if (z_end >= upper_level && z_start < upper_level) {
-                stop = true;
-                on_upper = true;
-            }
-        } else {  // z decreasing.
-            if (!(on_upper && first_edge) &&
-                z_start >= upper_level && z_end < upper_level) {
-                stop = true;
-                on_upper = true;
-            } else if (z_start >= lower_level && z_end < lower_level) {
-                stop = true;
-                on_upper = false;
-            }
-        }
-
-        first_edge = false;
-
-        if (!stop) {
-            // Move to next boundary edge, adding point to contour line.
-            edge = (edge+1) % (int)boundaries[boundary].size();
-            tri_edge = boundaries[boundary][edge];
-            contour_line.push_back(triang.get_point_coords(
-                                       triang.get_triangle_point(tri_edge)));
-        }
-    }
-
-    return on_upper;
-}
-
-void TriContourGenerator::follow_interior(ContourLine& contour_line,
-                                          TriEdge& tri_edge,
-                                          bool end_on_boundary,
-                                          const double& level,
-                                          bool on_upper)
-{
-    int& tri = tri_edge.tri;
-    int& edge = tri_edge.edge;
-
-    // Initial point.
-    contour_line.push_back(edge_interp(tri, edge, level));
-
-    while (true) {
-        int visited_index = tri;
-        if (on_upper)
-            visited_index += _triangulation.get_ntri();
-
-        // Check for end not on boundary.
-        if (!end_on_boundary && _interior_visited[visited_index])
-            break;  // Reached start point, so return.
-
-        // Determine edge by which to leave this triangle.
-        edge = get_exit_edge(tri, level, on_upper);
-        assert(edge >= 0 && edge < 3 && "Invalid exit edge");
-
-        _interior_visited[visited_index] = true;
-
-        // Append new point to point set.
-        assert(edge >= 0 && edge < 3 && "Invalid triangle edge");
-        contour_line.push_back(edge_interp(tri, edge, level));
-
-        // Move to next triangle.
-        TriEdge next_tri_edge = _triangulation.get_neighbor_edge(tri,edge);
-
-        // Check if ending on a boundary.
-        if (end_on_boundary && next_tri_edge.tri == -1)
-            break;
-
-        tri_edge = next_tri_edge;
-        assert(tri_edge.tri != -1 && "Invalid triangle for internal loop");
-    }
-}
-
-const TriContourGenerator::Boundaries& TriContourGenerator::get_boundaries() const
-{
-    return _triangulation.get_boundaries();
-}
-
-int TriContourGenerator::get_exit_edge(int tri,
-                                       const double& level,
-                                       bool on_upper) const
-{
-    assert(tri >= 0 && tri < _triangulation.get_ntri() &&
-           "Triangle index out of bounds.");
-
-    unsigned int config =
-        (get_z(_triangulation.get_triangle_point(tri, 0)) >= level) |
-        (get_z(_triangulation.get_triangle_point(tri, 1)) >= level) << 1 |
-        (get_z(_triangulation.get_triangle_point(tri, 2)) >= level) << 2;
-
-    if (on_upper) config = 7-config;
-
-    switch (config) {
-        case 0: return -1;
-        case 1: return  2;
-        case 2: return  0;
-        case 3: return  2;
-        case 4: return  1;
-        case 5: return  1;
-        case 6: return  0;
-        case 7: return -1;
-        default: assert(0 && "Invalid config value"); return -1;
-    }
-}
-
-const double& TriContourGenerator::get_z(int point) const
-{
-    assert(point >= 0 && point < _triangulation.get_npoints() &&
-           "Point index out of bounds.");
-    return _z.data()[point];
-}
-
-XY TriContourGenerator::interp(int point1,
-                               int point2,
-                               const double& level) const
-{
-    assert(point1 >= 0 && point1 < _triangulation.get_npoints() &&
-           "Point index 1 out of bounds.");
-    assert(point2 >= 0 && point2 < _triangulation.get_npoints() &&
-           "Point index 2 out of bounds.");
-    assert(point1 != point2 && "Identical points");
-    double fraction = (get_z(point2) - level) / (get_z(point2) - get_z(point1));
-    return _triangulation.get_point_coords(point1)*fraction +
-           _triangulation.get_point_coords(point2)*(1.0 - fraction);
-}
-
-
-
-TrapezoidMapTriFinder::TrapezoidMapTriFinder(Triangulation& triangulation)
-    : _triangulation(triangulation),
-      _points(0),
-      _tree(0)
-{}
-
-TrapezoidMapTriFinder::~TrapezoidMapTriFinder()
-{
-    clear();
-}
-
-bool
-TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge)
-{
-    std::vector<Trapezoid*> trapezoids;
-    if (!find_trapezoids_intersecting_edge(edge, trapezoids))
-        return false;
-    assert(!trapezoids.empty() && "No trapezoids intersect edge");
-
-    const Point* p = edge.left;
-    const Point* q = edge.right;
-    Trapezoid* left_old = 0;    // old trapezoid to the left.
-    Trapezoid* left_below = 0;  // below trapezoid to the left.
-    Trapezoid* left_above = 0;  // above trapezoid to the left.
-
-    // Iterate through trapezoids intersecting edge from left to right.
-    // Replace each old trapezoid with 2+ new trapezoids, and replace its
-    // corresponding nodes in the search tree with new nodes.
-    size_t ntraps = trapezoids.size();
-    for (size_t i = 0; i < ntraps; ++i) {
-        Trapezoid* old = trapezoids[i];  // old trapezoid to replace.
-        bool start_trap = (i == 0);
-        bool end_trap = (i == ntraps-1);
-        bool have_left = (start_trap && edge.left != old->left);
-        bool have_right = (end_trap && edge.right != old->right);
-
-        // Old trapezoid is replaced by up to 4 new trapezoids: left is to the
-        // left of the start point p, below/above are below/above the edge
-        // inserted, and right is to the right of the end point q.
-        Trapezoid* left = 0;
-        Trapezoid* below = 0;
-        Trapezoid* above = 0;
-        Trapezoid* right = 0;
-
-        // There are 4 different cases here depending on whether the old
-        // trapezoid in question is the start and/or end trapezoid of those
-        // that intersect the edge inserted.  There is some code duplication
-        // here but it is much easier to understand this way rather than
-        // interleave the 4 different cases with many more if-statements.
-        if (start_trap && end_trap) {
-            // Edge intersects a single trapezoid.
-            if (have_left)
-                left = new Trapezoid(old->left, p, old->below, old->above);
-            below = new Trapezoid(p, q, old->below, edge);
-            above = new Trapezoid(p, q, edge, old->above);
-            if (have_right)
-                right = new Trapezoid(q, old->right, old->below, old->above);
-
-            // Set pairs of trapezoid neighbours.
-            if (have_left) {
-                left->set_lower_left(old->lower_left);
-                left->set_upper_left(old->upper_left);
-                left->set_lower_right(below);
-                left->set_upper_right(above);
-            }
-            else {
-                below->set_lower_left(old->lower_left);
-                above->set_upper_left(old->upper_left);
-            }
-
-            if (have_right) {
-                right->set_lower_right(old->lower_right);
-                right->set_upper_right(old->upper_right);
-                below->set_lower_right(right);
-                above->set_upper_right(right);
-            }
-            else {
-                below->set_lower_right(old->lower_right);
-                above->set_upper_right(old->upper_right);
-            }
-        }
-        else if (start_trap) {
-            // Old trapezoid is the first of 2+ trapezoids that the edge
-            // intersects.
-            if (have_left)
-                left = new Trapezoid(old->left, p, old->below, old->above);
-            below = new Trapezoid(p, old->right, old->below, edge);
-            above = new Trapezoid(p, old->right, edge, old->above);
-
-            // Set pairs of trapezoid neighbours.
-            if (have_left) {
-                left->set_lower_left(old->lower_left);
-                left->set_upper_left(old->upper_left);
-                left->set_lower_right(below);
-                left->set_upper_right(above);
-            }
-            else {
-                below->set_lower_left(old->lower_left);
-                above->set_upper_left(old->upper_left);
-            }
-
-            below->set_lower_right(old->lower_right);
-            above->set_upper_right(old->upper_right);
-        }
-        else if (end_trap) {
-            // Old trapezoid is the last of 2+ trapezoids that the edge
-            // intersects.
-            if (left_below->below == old->below) {
-                below = left_below;
-                below->right = q;
-            }
-            else
-                below = new Trapezoid(old->left, q, old->below, edge);
-
-            if (left_above->above == old->above) {
-                above = left_above;
-                above->right = q;
-            }
-            else
-                above = new Trapezoid(old->left, q, edge, old->above);
-
-            if (have_right)
-                right = new Trapezoid(q, old->right, old->below, old->above);
-
-            // Set pairs of trapezoid neighbours.
-            if (have_right) {
-                right->set_lower_right(old->lower_right);
-                right->set_upper_right(old->upper_right);
-                below->set_lower_right(right);
-                above->set_upper_right(right);
-            }
-            else {
-                below->set_lower_right(old->lower_right);
-                above->set_upper_right(old->upper_right);
-            }
-
-            // Connect to new trapezoids replacing prevOld.
-            if (below != left_below) {
-                below->set_upper_left(left_below);
-                if (old->lower_left == left_old)
-                    below->set_lower_left(left_below);
-                else
-                    below->set_lower_left(old->lower_left);
-            }
-
-            if (above != left_above) {
-                above->set_lower_left(left_above);
-                if (old->upper_left == left_old)
-                    above->set_upper_left(left_above);
-                else
-                    above->set_upper_left(old->upper_left);
-            }
-        }
-        else {  // Middle trapezoid.
-            // Old trapezoid is neither the first nor last of the 3+ trapezoids
-            // that the edge intersects.
-            if (left_below->below == old->below) {
-                below = left_below;
-                below->right = old->right;
-            }
-            else
-                below = new Trapezoid(old->left, old->right, old->below, edge);
-
-            if (left_above->above == old->above) {
-                above = left_above;
-                above->right = old->right;
-            }
-            else
-                above = new Trapezoid(old->left, old->right, edge, old->above);
-
-            // Connect to new trapezoids replacing prevOld.
-            if (below != left_below) {  // below is new.
-                below->set_upper_left(left_below);
-                if (old->lower_left == left_old)
-                    below->set_lower_left(left_below);
-                else
-                    below->set_lower_left(old->lower_left);
-            }
-
-            if (above != left_above) {  // above is new.
-                above->set_lower_left(left_above);
-                if (old->upper_left == left_old)
-                    above->set_upper_left(left_above);
-                else
-                    above->set_upper_left(old->upper_left);
-            }
-
-            below->set_lower_right(old->lower_right);
-            above->set_upper_right(old->upper_right);
-        }
-
-        // Create new nodes to add to search tree.  Below and above trapezoids
-        // may already have owning trapezoid nodes, in which case reuse them.
-        Node* new_top_node = new Node(
-            &edge,
-            below == left_below ? below->trapezoid_node : new Node(below),
-            above == left_above ? above->trapezoid_node : new Node(above));
-        if (have_right)
-            new_top_node = new Node(q, new_top_node, new Node(right));
-        if (have_left)
-            new_top_node = new Node(p, new Node(left), new_top_node);
-
-        // Insert new_top_node in correct position or positions in search tree.
-        Node* old_node = old->trapezoid_node;
-        if (old_node == _tree)
-            _tree = new_top_node;
-        else
-            old_node->replace_with(new_top_node);
-
-        // old_node has been removed from all of its parents and is no longer
-        // needed.
-        assert(old_node->has_no_parents() && "Node should have no parents");
-        delete old_node;
-
-        // Clearing up.
-        if (!end_trap) {
-            // Prepare for next loop.
-            left_old = old;
-            left_above = above;
-            left_below = below;
-        }
-    }
-
-    return true;
-}
-
-void
-TrapezoidMapTriFinder::clear()
-{
-    delete [] _points;
-    _points = 0;
-
-    _edges.clear();
-
-    delete _tree;
-    _tree = 0;
-}
-
-TrapezoidMapTriFinder::TriIndexArray
-TrapezoidMapTriFinder::find_many(const CoordinateArray& x,
-                                 const CoordinateArray& y)
-{
-    if (x.ndim() != 1 || x.shape(0) != y.shape(0))
-        throw std::invalid_argument(
-            "x and y must be array-like with same shape");
-
-    // Create integer array to return.
-    auto n = x.shape(0);
-    TriIndexArray tri_indices_array(n);
-    auto tri_indices = tri_indices_array.mutable_unchecked<1>();
-    auto x_data = x.data();
-    auto y_data = y.data();
-
-    // Fill returned array.
-    for (py::ssize_t i = 0; i < n; ++i)
-        tri_indices(i) = find_one(XY(x_data[i], y_data[i]));
-
-    return tri_indices_array;
-}
-
-int
-TrapezoidMapTriFinder::find_one(const XY& xy)
-{
-    const Node* node = _tree->search(xy);
-    assert(node != 0 && "Search tree for point returned null node");
-    return node->get_tri();
-}
-
-bool
-TrapezoidMapTriFinder::find_trapezoids_intersecting_edge(
-    const Edge& edge,
-    std::vector<Trapezoid*>& trapezoids)
-{
-    // This is the FollowSegment algorithm of de Berg et al, with some extra
-    // checks to deal with simple colinear (i.e. invalid) triangles.
-    trapezoids.clear();
-    Trapezoid* trapezoid = _tree->search(edge);
-    if (trapezoid == 0) {
-        assert(trapezoid != 0 && "search(edge) returns null trapezoid");
-        return false;
-    }
-
-    trapezoids.push_back(trapezoid);
-    while (edge.right->is_right_of(*trapezoid->right)) {
-        int orient = edge.get_point_orientation(*trapezoid->right);
-        if (orient == 0) {
-            if (edge.point_below == trapezoid->right)
-                orient = +1;
-            else if (edge.point_above == trapezoid->right)
-                orient = -1;
-            else {
-                assert(0 && "Unable to deal with point on edge");
-                return false;
-            }
-        }
-
-        if (orient == -1)
-            trapezoid = trapezoid->lower_right;
-        else if (orient == +1)
-            trapezoid = trapezoid->upper_right;
-
-        if (trapezoid == 0) {
-            assert(0 && "Expected trapezoid neighbor");
-            return false;
-        }
-        trapezoids.push_back(trapezoid);
-    }
-
-    return true;
-}
-
-py::list
-TrapezoidMapTriFinder::get_tree_stats()
-{
-    NodeStats stats;
-    _tree->get_stats(0, stats);
-
-    py::list ret(7);
-    ret[0] = stats.node_count;
-    ret[1] = stats.unique_nodes.size(),
-    ret[2] = stats.trapezoid_count,
-    ret[3] = stats.unique_trapezoid_nodes.size(),
-    ret[4] = stats.max_parent_count,
-    ret[5] = stats.max_depth,
-    ret[6] = stats.sum_trapezoid_depth / stats.trapezoid_count;
-    return ret;
-}
-
-void
-TrapezoidMapTriFinder::initialize()
-{
-    clear();
-    const Triangulation& triang = _triangulation;
-
-    // Set up points array, which contains all of the points in the
-    // triangulation plus the 4 corners of the enclosing rectangle.
-    int npoints = triang.get_npoints();
-    _points = new Point[npoints + 4];
-    BoundingBox bbox;
-    for (int i = 0; i < npoints; ++i) {
-        XY xy = triang.get_point_coords(i);
-        // Avoid problems with -0.0 values different from 0.0
-        if (xy.x == -0.0)
-            xy.x = 0.0;
-        if (xy.y == -0.0)
-            xy.y = 0.0;
-        _points[i] = Point(xy);
-        bbox.add(xy);
-    }
-
-    // Last 4 points are corner points of enclosing rectangle.  Enclosing
-    // rectangle made slightly larger in case corner points are already in the
-    // triangulation.
-    if (bbox.empty) {
-        bbox.add(XY(0.0, 0.0));
-        bbox.add(XY(1.0, 1.0));
-    }
-    else {
-        const double small = 0.1;  // Any value > 0.0
-        bbox.expand( (bbox.upper - bbox.lower)*small );
-    }
-    _points[npoints  ] = Point(bbox.lower);                  // SW point.
-    _points[npoints+1] = Point(bbox.upper.x, bbox.lower.y);  // SE point.
-    _points[npoints+2] = Point(bbox.lower.x, bbox.upper.y);  // NW point.
-    _points[npoints+3] = Point(bbox.upper);                  // NE point.
-
-    // Set up edges array.
-    // First the bottom and top edges of the enclosing rectangle.
-    _edges.push_back(Edge(&_points[npoints],  &_points[npoints+1],-1,-1,0,0));
-    _edges.push_back(Edge(&_points[npoints+2],&_points[npoints+3],-1,-1,0,0));
-
-    // Add all edges in the triangulation that point to the right.  Do not
-    // explicitly include edges that point to the left as the neighboring
-    // triangle will supply that, unless there is no such neighbor.
-    int ntri = triang.get_ntri();
-    for (int tri = 0; tri < ntri; ++tri) {
-        if (!triang.is_masked(tri)) {
-            for (int edge = 0; edge < 3; ++edge) {
-                Point* start = _points + triang.get_triangle_point(tri,edge);
-                Point* end   = _points +
-                               triang.get_triangle_point(tri,(edge+1)%3);
-                Point* other = _points +
-                               triang.get_triangle_point(tri,(edge+2)%3);
-                TriEdge neighbor = triang.get_neighbor_edge(tri,edge);
-                if (end->is_right_of(*start)) {
-                    const Point* neighbor_point_below = (neighbor.tri == -1) ?
-                        0 : _points + triang.get_triangle_point(
-                                          neighbor.tri, (neighbor.edge+2)%3);
-                    _edges.push_back(Edge(start, end, neighbor.tri, tri,
-                                          neighbor_point_below, other));
-                }
-                else if (neighbor.tri == -1)
-                    _edges.push_back(Edge(end, start, tri, -1, other, 0));
-
-                // Set triangle associated with start point if not already set.
-                if (start->tri == -1)
-                    start->tri = tri;
-            }
-        }
-    }
-
-    // Initial trapezoid is enclosing rectangle.
-    _tree = new Node(new Trapezoid(&_points[npoints], &_points[npoints+1],
-                                   _edges[0], _edges[1]));
-    _tree->assert_valid(false);
-
-    // Randomly shuffle all edges other than first 2.
-    std::mt19937 rng(1234);
-    std::shuffle(_edges.begin()+2, _edges.end(), rng);
-
-    // Add edges, one at a time, to tree.
-    size_t nedges = _edges.size();
-    for (size_t index = 2; index < nedges; ++index) {
-        if (!add_edge_to_tree(_edges[index]))
-            throw std::runtime_error("Triangulation is invalid");
-        _tree->assert_valid(index == nedges-1);
-    }
-}
-
-void
-TrapezoidMapTriFinder::print_tree()
-{
-    assert(_tree != 0 && "Null Node tree");
-    _tree->print();
-}
-
-TrapezoidMapTriFinder::Edge::Edge(const Point* left_,
-                                  const Point* right_,
-                                  int triangle_below_,
-                                  int triangle_above_,
-                                  const Point* point_below_,
-                                  const Point* point_above_)
-    : left(left_),
-      right(right_),
-      triangle_below(triangle_below_),
-      triangle_above(triangle_above_),
-      point_below(point_below_),
-      point_above(point_above_)
-{
-    assert(left != 0 && "Null left point");
-    assert(right != 0 && "Null right point");
-    assert(right->is_right_of(*left) && "Incorrect point order");
-    assert(triangle_below >= -1 && "Invalid triangle below index");
-    assert(triangle_above >= -1 && "Invalid triangle above index");
-}
-
-int
-TrapezoidMapTriFinder::Edge::get_point_orientation(const XY& xy) const
-{
-    double cross_z = (xy - *left).cross_z(*right - *left);
-    return (cross_z > 0.0) ? +1 : ((cross_z < 0.0) ? -1 : 0);
-}
-
-double
-TrapezoidMapTriFinder::Edge::get_slope() const
-{
-    // Divide by zero is acceptable here.
-    XY diff = *right - *left;
-    return diff.y / diff.x;
-}
-
-double
-TrapezoidMapTriFinder::Edge::get_y_at_x(const double& x) const
-{
-    if (left->x == right->x) {
-        // If edge is vertical, return lowest y from left point.
-        assert(x == left->x && "x outside of edge");
-        return left->y;
-    }
-    else {
-        // Equation of line: left + lambda*(right - left) = xy.
-        // i.e. left.x + lambda(right.x - left.x) = x and similar for y.
-        double lambda = (x - left->x) / (right->x - left->x);
-        assert(lambda >= 0 && lambda <= 1.0 && "Lambda out of bounds");
-        return left->y + lambda*(right->y - left->y);
-    }
-}
-
-bool
-TrapezoidMapTriFinder::Edge::has_point(const Point* point) const
-{
-    assert(point != 0 && "Null point");
-    return (left == point || right == point);
-}
-
-bool
-TrapezoidMapTriFinder::Edge::operator==(const Edge& other) const
-{
-    return this == &other;
-}
-
-void
-TrapezoidMapTriFinder::Edge::print_debug() const
-{
-    std::cout << "Edge " << *this << " tri_below=" << triangle_below
-        << " tri_above=" << triangle_above << std::endl;
-}
-
-TrapezoidMapTriFinder::Node::Node(const Point* point, Node* left, Node* right)
-    : _type(Type_XNode)
-{
-    assert(point != 0 && "Invalid point");
-    assert(left != 0 && "Invalid left node");
-    assert(right != 0 && "Invalid right node");
-    _union.xnode.point = point;
-    _union.xnode.left = left;
-    _union.xnode.right = right;
-    left->add_parent(this);
-    right->add_parent(this);
-}
-
-TrapezoidMapTriFinder::Node::Node(const Edge* edge, Node* below, Node* above)
-    : _type(Type_YNode)
-{
-    assert(edge != 0 && "Invalid edge");
-    assert(below != 0 && "Invalid below node");
-    assert(above != 0 && "Invalid above node");
-    _union.ynode.edge = edge;
-    _union.ynode.below = below;
-    _union.ynode.above = above;
-    below->add_parent(this);
-    above->add_parent(this);
-}
-
-TrapezoidMapTriFinder::Node::Node(Trapezoid* trapezoid)
-    : _type(Type_TrapezoidNode)
-{
-    assert(trapezoid != 0 && "Null Trapezoid");
-    _union.trapezoid = trapezoid;
-    trapezoid->trapezoid_node = this;
-}
-
-TrapezoidMapTriFinder::Node::~Node()
-{
-    switch (_type) {
-        case Type_XNode:
-            if (_union.xnode.left->remove_parent(this))
-                delete _union.xnode.left;
-            if (_union.xnode.right->remove_parent(this))
-                delete _union.xnode.right;
-            break;
-        case Type_YNode:
-            if (_union.ynode.below->remove_parent(this))
-                delete _union.ynode.below;
-            if (_union.ynode.above->remove_parent(this))
-                delete _union.ynode.above;
-            break;
-        case Type_TrapezoidNode:
-            delete _union.trapezoid;
-            break;
-    }
-}
-
-void
-TrapezoidMapTriFinder::Node::add_parent(Node* parent)
-{
-    assert(parent != 0 && "Null parent");
-    assert(parent != this && "Cannot be parent of self");
-    assert(!has_parent(parent) && "Parent already in collection");
-    _parents.push_back(parent);
-}
-
-void
-TrapezoidMapTriFinder::Node::assert_valid(bool tree_complete) const
-{
-#ifndef NDEBUG
-    // Check parents.
-    for (Parents::const_iterator it = _parents.begin();
-         it != _parents.end(); ++it) {
-        Node* parent = *it;
-        assert(parent != this && "Cannot be parent of self");
-        assert(parent->has_child(this) && "Parent missing child");
-    }
-
-    // Check children, and recurse.
-    switch (_type) {
-        case Type_XNode:
-            assert(_union.xnode.left != 0 && "Null left child");
-            assert(_union.xnode.left->has_parent(this) && "Incorrect parent");
-            assert(_union.xnode.right != 0 && "Null right child");
-            assert(_union.xnode.right->has_parent(this) && "Incorrect parent");
-            _union.xnode.left->assert_valid(tree_complete);
-            _union.xnode.right->assert_valid(tree_complete);
-            break;
-        case Type_YNode:
-            assert(_union.ynode.below != 0 && "Null below child");
-            assert(_union.ynode.below->has_parent(this) && "Incorrect parent");
-            assert(_union.ynode.above != 0 && "Null above child");
-            assert(_union.ynode.above->has_parent(this) && "Incorrect parent");
-            _union.ynode.below->assert_valid(tree_complete);
-            _union.ynode.above->assert_valid(tree_complete);
-            break;
-        case Type_TrapezoidNode:
-            assert(_union.trapezoid != 0 && "Null trapezoid");
-            assert(_union.trapezoid->trapezoid_node == this &&
-                   "Incorrect trapezoid node");
-            _union.trapezoid->assert_valid(tree_complete);
-            break;
-    }
-#endif
-}
-
-void
-TrapezoidMapTriFinder::Node::get_stats(int depth,
-                                       NodeStats& stats) const
-{
-    stats.node_count++;
-    if (depth > stats.max_depth)
-        stats.max_depth = depth;
-    bool new_node = stats.unique_nodes.insert(this).second;
-    if (new_node)
-        stats.max_parent_count = std::max(stats.max_parent_count,
-                                          static_cast<long>(_parents.size()));
-
-    switch (_type) {
-        case Type_XNode:
-            _union.xnode.left->get_stats(depth+1, stats);
-            _union.xnode.right->get_stats(depth+1, stats);
-            break;
-        case Type_YNode:
-            _union.ynode.below->get_stats(depth+1, stats);
-            _union.ynode.above->get_stats(depth+1, stats);
-            break;
-        default:  // Type_TrapezoidNode:
-            stats.unique_trapezoid_nodes.insert(this);
-            stats.trapezoid_count++;
-            stats.sum_trapezoid_depth += depth;
-            break;
-    }
-}
-
-int
-TrapezoidMapTriFinder::Node::get_tri() const
-{
-    switch (_type) {
-        case Type_XNode:
-            return _union.xnode.point->tri;
-        case Type_YNode:
-            if (_union.ynode.edge->triangle_above != -1)
-                return _union.ynode.edge->triangle_above;
-            else
-                return _union.ynode.edge->triangle_below;
-        default:  // Type_TrapezoidNode:
-            assert(_union.trapezoid->below.triangle_above ==
-                   _union.trapezoid->above.triangle_below &&
-                   "Inconsistent triangle indices from trapezoid edges");
-            return _union.trapezoid->below.triangle_above;
-    }
-}
-
-bool
-TrapezoidMapTriFinder::Node::has_child(const Node* child) const
-{
-    assert(child != 0 && "Null child node");
-    switch (_type) {
-        case Type_XNode:
-            return (_union.xnode.left == child || _union.xnode.right == child);
-        case Type_YNode:
-            return (_union.ynode.below == child ||
-                    _union.ynode.above == child);
-        default:  // Type_TrapezoidNode:
-            return false;
-    }
-}
-
-bool
-TrapezoidMapTriFinder::Node::has_no_parents() const
-{
-    return _parents.empty();
-}
-
-bool
-TrapezoidMapTriFinder::Node::has_parent(const Node* parent) const
-{
-    return (std::find(_parents.begin(), _parents.end(), parent) !=
-            _parents.end());
-}
-
-void
-TrapezoidMapTriFinder::Node::print(int depth /* = 0 */) const
-{
-    for (int i = 0; i < depth; ++i) std::cout << "  ";
-    switch (_type) {
-        case Type_XNode:
-            std::cout << "XNode " << *_union.xnode.point << std::endl;
-            _union.xnode.left->print(depth + 1);
-            _union.xnode.right->print(depth + 1);
-            break;
-        case Type_YNode:
-            std::cout << "YNode " << *_union.ynode.edge << std::endl;
-            _union.ynode.below->print(depth + 1);
-            _union.ynode.above->print(depth + 1);
-            break;
-        case Type_TrapezoidNode:
-            std::cout << "Trapezoid ll="
-                << _union.trapezoid->get_lower_left_point()  << " lr="
-                << _union.trapezoid->get_lower_right_point() << " ul="
-                << _union.trapezoid->get_upper_left_point()  << " ur="
-                << _union.trapezoid->get_upper_right_point() << std::endl;
-            break;
-    }
-}
-
-bool
-TrapezoidMapTriFinder::Node::remove_parent(Node* parent)
-{
-    assert(parent != 0 && "Null parent");
-    assert(parent != this && "Cannot be parent of self");
-    Parents::iterator it = std::find(_parents.begin(), _parents.end(), parent);
-    assert(it != _parents.end() && "Parent not in collection");
-    _parents.erase(it);
-    return _parents.empty();
-}
-
-void
-TrapezoidMapTriFinder::Node::replace_child(Node* old_child, Node* new_child)
-{
-    switch (_type) {
-        case Type_XNode:
-            assert((_union.xnode.left == old_child ||
-                    _union.xnode.right == old_child) && "Not a child Node");
-            assert(new_child != 0 && "Null child node");
-            if (_union.xnode.left == old_child)
-                _union.xnode.left = new_child;
-            else
-                _union.xnode.right = new_child;
-            break;
-        case Type_YNode:
-            assert((_union.ynode.below == old_child ||
-                    _union.ynode.above == old_child) && "Not a child node");
-            assert(new_child != 0 && "Null child node");
-            if (_union.ynode.below == old_child)
-                _union.ynode.below = new_child;
-            else
-                _union.ynode.above = new_child;
-            break;
-        case Type_TrapezoidNode:
-            assert(0 && "Invalid type for this operation");
-            break;
-    }
-    old_child->remove_parent(this);
-    new_child->add_parent(this);
-}
-
-void
-TrapezoidMapTriFinder::Node::replace_with(Node* new_node)
-{
-    assert(new_node != 0 && "Null replacement node");
-    // Replace child of each parent with new_node.  As each has parent has its
-    // child replaced it is removed from the _parents collection.
-    while (!_parents.empty())
-        _parents.front()->replace_child(this, new_node);
-}
-
-const TrapezoidMapTriFinder::Node*
-TrapezoidMapTriFinder::Node::search(const XY& xy)
-{
-    switch (_type) {
-        case Type_XNode:
-            if (xy == *_union.xnode.point)
-                return this;
-            else if (xy.is_right_of(*_union.xnode.point))
-                return _union.xnode.right->search(xy);
-            else
-                return _union.xnode.left->search(xy);
-        case Type_YNode: {
-            int orient = _union.ynode.edge->get_point_orientation(xy);
-            if (orient == 0)
-                return this;
-            else if (orient < 0)
-                return _union.ynode.above->search(xy);
-            else
-                return _union.ynode.below->search(xy);
-        }
-        default:  // Type_TrapezoidNode:
-            return this;
-    }
-}
-
-TrapezoidMapTriFinder::Trapezoid*
-TrapezoidMapTriFinder::Node::search(const Edge& edge)
-{
-    switch (_type) {
-        case Type_XNode:
-            if (edge.left == _union.xnode.point)
-                return _union.xnode.right->search(edge);
-            else {
-                if (edge.left->is_right_of(*_union.xnode.point))
-                    return _union.xnode.right->search(edge);
-                else
-                    return _union.xnode.left->search(edge);
-            }
-        case Type_YNode:
-            if (edge.left == _union.ynode.edge->left) {
-                // Coinciding left edge points.
-                if (edge.get_slope() == _union.ynode.edge->get_slope()) {
-                    if (_union.ynode.edge->triangle_above ==
-                        edge.triangle_below)
-                        return _union.ynode.above->search(edge);
-                    else if (_union.ynode.edge->triangle_below ==
-                             edge.triangle_above)
-                        return _union.ynode.below->search(edge);
-                    else {
-                        assert(0 &&
-                               "Invalid triangulation, common left points");
-                        return 0;
-                    }
-                }
-                if (edge.get_slope() > _union.ynode.edge->get_slope())
-                    return _union.ynode.above->search(edge);
-                else
-                    return _union.ynode.below->search(edge);
-            }
-            else if (edge.right == _union.ynode.edge->right) {
-                // Coinciding right edge points.
-                if (edge.get_slope() == _union.ynode.edge->get_slope()) {
-                    if (_union.ynode.edge->triangle_above ==
-                        edge.triangle_below)
-                        return _union.ynode.above->search(edge);
-                    else if (_union.ynode.edge->triangle_below ==
-                             edge.triangle_above)
-                        return _union.ynode.below->search(edge);
-                    else {
-                        assert(0 &&
-                               "Invalid triangulation, common right points");
-                        return 0;
-                    }
-                }
-                if (edge.get_slope() > _union.ynode.edge->get_slope())
-                    return _union.ynode.below->search(edge);
-                else
-                    return _union.ynode.above->search(edge);
-            }
-            else {
-                int orient =
-                    _union.ynode.edge->get_point_orientation(*edge.left);
-                if (orient == 0) {
-                    // edge.left lies on _union.ynode.edge
-                    if (_union.ynode.edge->point_above != 0 &&
-                        edge.has_point(_union.ynode.edge->point_above))
-                        orient = -1;
-                    else if (_union.ynode.edge->point_below != 0 &&
-                             edge.has_point(_union.ynode.edge->point_below))
-                        orient = +1;
-                    else {
-                        assert(0 && "Invalid triangulation, point on edge");
-                        return 0;
-                    }
-                }
-                if (orient < 0)
-                    return _union.ynode.above->search(edge);
-                else
-                    return _union.ynode.below->search(edge);
-            }
-        default:  // Type_TrapezoidNode:
-            return _union.trapezoid;
-    }
-}
-
-TrapezoidMapTriFinder::Trapezoid::Trapezoid(const Point* left_,
-                                            const Point* right_,
-                                            const Edge& below_,
-                                            const Edge& above_)
-    : left(left_), right(right_), below(below_), above(above_),
-      lower_left(0), lower_right(0), upper_left(0), upper_right(0),
-      trapezoid_node(0)
-{
-    assert(left != 0 && "Null left point");
-    assert(right != 0 && "Null right point");
-    assert(right->is_right_of(*left) && "Incorrect point order");
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::assert_valid(bool tree_complete) const
-{
-#ifndef NDEBUG
-    assert(left != 0 && "Null left point");
-    assert(right != 0 && "Null right point");
-
-    if (lower_left != 0) {
-        assert(lower_left->below == below &&
-               lower_left->lower_right == this &&
-               "Incorrect lower_left trapezoid");
-        assert(get_lower_left_point() == lower_left->get_lower_right_point() &&
-               "Incorrect lower left point");
-    }
-
-    if (lower_right != 0) {
-        assert(lower_right->below == below &&
-               lower_right->lower_left == this &&
-               "Incorrect lower_right trapezoid");
-        assert(get_lower_right_point() == lower_right->get_lower_left_point() &&
-               "Incorrect lower right point");
-    }
-
-    if (upper_left != 0) {
-        assert(upper_left->above == above &&
-               upper_left->upper_right == this &&
-               "Incorrect upper_left trapezoid");
-        assert(get_upper_left_point() == upper_left->get_upper_right_point() &&
-               "Incorrect upper left point");
-    }
-
-    if (upper_right != 0) {
-        assert(upper_right->above == above &&
-               upper_right->upper_left == this &&
-               "Incorrect upper_right trapezoid");
-        assert(get_upper_right_point() == upper_right->get_upper_left_point() &&
-               "Incorrect upper right point");
-    }
-
-    assert(trapezoid_node != 0 && "Null trapezoid_node");
-
-    if (tree_complete) {
-        assert(below.triangle_above == above.triangle_below &&
-               "Inconsistent triangle indices from trapezoid edges");
-    }
-#endif
-}
-
-XY
-TrapezoidMapTriFinder::Trapezoid::get_lower_left_point() const
-{
-    double x = left->x;
-    return XY(x, below.get_y_at_x(x));
-}
-
-XY
-TrapezoidMapTriFinder::Trapezoid::get_lower_right_point() const
-{
-    double x = right->x;
-    return XY(x, below.get_y_at_x(x));
-}
-
-XY
-TrapezoidMapTriFinder::Trapezoid::get_upper_left_point() const
-{
-    double x = left->x;
-    return XY(x, above.get_y_at_x(x));
-}
-
-XY
-TrapezoidMapTriFinder::Trapezoid::get_upper_right_point() const
-{
-    double x = right->x;
-    return XY(x, above.get_y_at_x(x));
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::print_debug() const
-{
-    std::cout << "Trapezoid " << this
-        << " left=" << *left
-        << " right=" << *right
-        << " below=" << below
-        << " above=" << above
-        << " ll=" << lower_left
-        << " lr=" << lower_right
-        << " ul=" << upper_left
-        << " ur=" << upper_right
-        << " node=" << trapezoid_node
-        << " llp=" << get_lower_left_point()
-        << " lrp=" << get_lower_right_point()
-        << " ulp=" << get_upper_left_point()
-        << " urp=" << get_upper_right_point() << std::endl;
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::set_lower_left(Trapezoid* lower_left_)
-{
-    lower_left = lower_left_;
-    if (lower_left != 0)
-        lower_left->lower_right = this;
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::set_lower_right(Trapezoid* lower_right_)
-{
-    lower_right = lower_right_;
-    if (lower_right != 0)
-        lower_right->lower_left = this;
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::set_upper_left(Trapezoid* upper_left_)
-{
-    upper_left = upper_left_;
-    if (upper_left != 0)
-        upper_left->upper_right = this;
-}
-
-void
-TrapezoidMapTriFinder::Trapezoid::set_upper_right(Trapezoid* upper_right_)
-{
-    upper_right = upper_right_;
-    if (upper_right != 0)
-        upper_right->upper_left = this;
-}
diff --git a/contrib/python/matplotlib/py3/src/tri/_tri.h b/contrib/python/matplotlib/py3/src/tri/_tri.h
deleted file mode 100644
index c176b4c0e8..0000000000
--- a/contrib/python/matplotlib/py3/src/tri/_tri.h
+++ /dev/null
@@ -1,799 +0,0 @@
-/*
- * Unstructured triangular grid functions, particularly contouring.
- *
- * There are two main classes: Triangulation and TriContourGenerator.
- *
- * Triangulation
- * -------------
- * Triangulation is an unstructured triangular grid with npoints and ntri
- * triangles.  It consists of point x and y coordinates, and information about
- * the triangulation stored in an integer array of shape (ntri,3) called
- * triangles.  Each triangle is represented by three point indices (in the
- * range 0 to npoints-1) that comprise the triangle, ordered anticlockwise.
- * There is an optional mask of length ntri which can be used to mask out
- * triangles and has the same result as removing those triangles from the
- * 'triangles' array.
- *
- * A particular edge of a triangulation is termed a TriEdge, which is a
- * triangle index and an edge index in the range 0 to 2.  TriEdge(tri,edge)
- * refers to the edge that starts at point index triangles(tri,edge) and ends
- * at point index triangles(tri,(edge+1)%3).
- *
- * Various derived fields are calculated when they are first needed.  The
- * triangle connectivity is stored in a neighbors array of shape (ntri,3) such
- * that neighbors(tri,edge) is the index of the triangle that adjoins the
- * TriEdge(tri,edge), or -1 if there is no such neighbor.
- *
- * A triangulation has one or more boundaries, each of which is a 1D array of
- * the TriEdges that comprise the boundary, in order following the boundary
- * with non-masked triangles on the left.
- *
- * TriContourGenerator
- * -------------------
- * A TriContourGenerator generates contours for a particular Triangulation.
- * The process followed is different for non-filled and filled contours, with
- * one and two contour levels respectively.  In both cases boundary contour
- * lines are found first, then interior lines.
- *
- * Boundary lines start and end on a boundary.  They are found by traversing
- * the triangulation boundary edges until a suitable start point is found, and
- * then the contour line is followed across the interior of the triangulation
- * until it ends on another boundary edge.  For a non-filled contour this
- * completes a line, whereas a filled contour continues by following the
- * boundary around until either another boundary start point is found or the
- * start of the contour line is reached.  Filled contour generation stores
- * boolean flags to indicate which boundary edges have already been traversed
- * so that they are not dealt with twice.  Similar flags are used to indicate
- * which triangles have been used when following interior lines.
- *
- * Interior lines do not intersect any boundaries.  They are found by
- * traversing all triangles that have not yet been visited until a suitable
- * starting point is found, and then the contour line is followed across the
- * interior of the triangulation until it returns to the start point.  For
- * filled contours this process is repeated for both lower and upper contour
- * levels, and the direction of traversal is reversed for upper contours.
- *
- * Working out in which direction a contour line leaves a triangle uses the
- * a lookup table.  A triangle has three points, each of which has a z-value
- * which is either less than the contour level or not.  Hence there are 8
- * configurations to deal with, 2 of which do not have a contour line (all
- * points below or above (including the same as) the contour level) and 6 that
- * do.  See the function get_exit_edge for details.
- */
-#ifndef MPL_TRI_H
-#define MPL_TRI_H
-
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-
-#include <iostream>
-#include <list>
-#include <map>
-#include <set>
-#include <vector>
-
-namespace py = pybind11;
-
-
-/* An edge of a triangle consisting of an triangle index in the range 0 to
- * ntri-1 and an edge index in the range 0 to 2.  Edge i goes from the
- * triangle's point i to point (i+1)%3. */
-struct TriEdge
-{
-    TriEdge();
-    TriEdge(int tri_, int edge_);
-    bool operator<(const TriEdge& other) const;
-    bool operator==(const TriEdge& other) const;
-    bool operator!=(const TriEdge& other) const;
-    friend std::ostream& operator<<(std::ostream& os, const TriEdge& tri_edge);
-
-    int tri, edge;
-};
-
-// 2D point with x,y coordinates.
-struct XY
-{
-    XY();
-    XY(const double& x_, const double& y_);
-    double angle() const;           // Angle in radians with respect to x-axis.
-    double cross_z(const XY& other) const;     // z-component of cross product.
-    bool is_right_of(const XY& other) const;   // Compares x then y.
-    bool operator==(const XY& other) const;
-    bool operator!=(const XY& other) const;
-    XY operator*(const double& multiplier) const;
-    const XY& operator+=(const XY& other);
-    const XY& operator-=(const XY& other);
-    XY operator+(const XY& other) const;
-    XY operator-(const XY& other) const;
-    friend std::ostream& operator<<(std::ostream& os, const XY& xy);
-
-    double x, y;
-};
-
-// 3D point with x,y,z coordinates.
-struct XYZ
-{
-    XYZ(const double& x_, const double& y_, const double& z_);
-    XYZ cross(const XYZ& other) const;
-    double dot(const XYZ& other) const;
-    XYZ operator-(const XYZ& other) const;
-    friend std::ostream& operator<<(std::ostream& os, const XYZ& xyz);
-
-    double x, y, z;
-};
-
-// 2D bounding box, which may be empty.
-class BoundingBox
-{
-public:
-    BoundingBox();
-    void add(const XY& point);
-    void expand(const XY& delta);
-
-    // Consider these member variables read-only.
-    bool empty;
-    XY lower, upper;
-};
-
-/* A single line of a contour, which may be a closed line loop or an open line
- * strip.  Identical adjacent points are avoided using push_back(), and a closed
- * line loop should also not have identical first and last points. */
-class ContourLine : public std::vector<XY>
-{
-public:
-    ContourLine();
-    void push_back(const XY& point);
-    void write() const;
-};
-
-// A Contour is a collection of zero or more ContourLines.
-typedef std::vector<ContourLine> Contour;
-
-// Debug contour writing function.
-void write_contour(const Contour& contour);
-
-
-
-
-/* Triangulation with npoints points and ntri triangles.  Derived fields are
- * calculated when they are first needed. */
-class Triangulation
-{
-public:
-    typedef py::array_t<double, py::array::c_style | py::array::forcecast> CoordinateArray;
-    typedef py::array_t<double, py::array::c_style | py::array::forcecast> TwoCoordinateArray;
-    typedef py::array_t<int,    py::array::c_style | py::array::forcecast> TriangleArray;
-    typedef py::array_t<bool,   py::array::c_style | py::array::forcecast> MaskArray;
-    typedef py::array_t<int,    py::array::c_style | py::array::forcecast> EdgeArray;
-    typedef py::array_t<int,    py::array::c_style | py::array::forcecast> NeighborArray;
-
-    /* A single boundary is a vector of the TriEdges that make up that boundary
-     * following it around with unmasked triangles on the left. */
-    typedef std::vector<TriEdge> Boundary;
-    typedef std::vector<Boundary> Boundaries;
-
-    /* Constructor with optional mask, edges and neighbors.  The latter two
-     * are calculated when first needed.
-     *   x: double array of shape (npoints) of points' x-coordinates.
-     *   y: double array of shape (npoints) of points' y-coordinates.
-     *   triangles: int array of shape (ntri,3) of triangle point indices.
-     *              Those ordered clockwise are changed to be anticlockwise.
-     *   mask: Optional bool array of shape (ntri) indicating which triangles
-     *         are masked.
-     *   edges: Optional int array of shape (?,2) of start and end point
-     *          indices, each edge (start,end and end,start) appearing only
-     *          once.
-     *   neighbors: Optional int array of shape (ntri,3) indicating which
-     *              triangles are the neighbors of which TriEdges, or -1 if
-     *              there is no such neighbor.
-     *   correct_triangle_orientations: Whether or not should correct triangle
-     *                                  orientations so that vertices are
-     *                                  ordered anticlockwise. */
-    Triangulation(const CoordinateArray& x,
-                  const CoordinateArray& y,
-                  const TriangleArray& triangles,
-                  const MaskArray& mask,
-                  const EdgeArray& edges,
-                  const NeighborArray& neighbors,
-                  bool correct_triangle_orientations);
-
-    /* Calculate plane equation coefficients for all unmasked triangles from
-     * the point (x,y) coordinates and point z-array of shape (npoints) passed
-     * in via the args.  Returned array has shape (npoints,3) and allows
-     * z-value at (x,y) coordinates in triangle tri to be calculated using
-     *      z = array[tri,0]*x + array[tri,1]*y + array[tri,2]. */
-    TwoCoordinateArray calculate_plane_coefficients(const CoordinateArray& z);
-
-    // Return the boundaries collection, creating it if necessary.
-    const Boundaries& get_boundaries() const;
-
-    // Return which boundary and boundary edge the specified TriEdge is.
-    void get_boundary_edge(const TriEdge& triEdge,
-                           int& boundary,
-                           int& edge) const;
-
-    /* Return the edges array, creating it if necessary. */
-    EdgeArray& get_edges();
-
-    /* Return the triangle index of the neighbor of the specified triangle
-     * edge. */
-    int get_neighbor(int tri, int edge) const;
-
-    /* Return the TriEdge that is the neighbor of the specified triangle edge,
-     * or TriEdge(-1,-1) if there is no such neighbor. */
-    TriEdge get_neighbor_edge(int tri, int edge) const;
-
-    /* Return the neighbors array, creating it if necessary. */
-    NeighborArray& get_neighbors();
-
-    // Return the number of points in this triangulation.
-    int get_npoints() const;
-
-    // Return the number of triangles in this triangulation.
-    int get_ntri() const;
-
-    /* Return the index of the point that is at the start of the specified
-     * triangle edge. */
-    int get_triangle_point(int tri, int edge) const;
-    int get_triangle_point(const TriEdge& tri_edge) const;
-
-    // Return the coordinates of the specified point index.
-    XY get_point_coords(int point) const;
-
-    // Indicates if the specified triangle is masked or not.
-    bool is_masked(int tri) const;
-
-    /* Set or clear the mask array.  Clears various derived fields so they are
-     * recalculated when next needed.
-     *   mask: bool array of shape (ntri) indicating which triangles are
-     *         masked, or an empty array to clear mask. */
-    void set_mask(const MaskArray& mask);
-
-    // Debug function to write boundaries.
-    void write_boundaries() const;
-
-private:
-    // An edge of a triangulation, composed of start and end point indices.
-    struct Edge
-    {
-        Edge() : start(-1), end(-1) {}
-        Edge(int start_, int end_) : start(start_), end(end_) {}
-        bool operator<(const Edge& other) const {
-            return start != other.start ? start < other.start : end < other.end;
-        }
-        int start, end;
-    };
-
-    /* An edge of a boundary of a triangulation, composed of a boundary index
-     * and an edge index within that boundary.  Used to index into the
-     * boundaries collection to obtain the corresponding TriEdge. */
-    struct BoundaryEdge
-    {
-        BoundaryEdge() : boundary(-1), edge(-1) {}
-        BoundaryEdge(int boundary_, int edge_)
-            : boundary(boundary_), edge(edge_) {}
-        int boundary, edge;
-    };
-
-    /* Calculate the boundaries collection.  Should normally be accessed via
-     * get_boundaries(), which will call this function if necessary. */
-    void calculate_boundaries();
-
-    /* Calculate the edges array.  Should normally be accessed via
-     * get_edges(), which will call this function if necessary. */
-    void calculate_edges();
-
-    /* Calculate the neighbors array. Should normally be accessed via
-     * get_neighbors(), which will call this function if necessary. */
-    void calculate_neighbors();
-
-    /* Correct each triangle so that the vertices are ordered in an
-     * anticlockwise manner. */
-    void correct_triangles();
-
-    /* Determine which edge index (0,1 or 2) the specified point index is in
-     * the specified triangle, or -1 if the point is not in the triangle. */
-    int get_edge_in_triangle(int tri, int point) const;
-
-    bool has_edges() const;
-
-    bool has_mask() const;
-
-    bool has_neighbors() const;
-
-
-    // Variables shared with python, always set.
-    CoordinateArray _x, _y;    // double array (npoints).
-    TriangleArray _triangles;  // int array (ntri,3) of triangle point indices,
-                               //     ordered anticlockwise.
-
-    // Variables shared with python, may be unset (size == 0).
-    MaskArray _mask;           // bool array (ntri).
-
-    // Derived variables shared with python, may be unset (size == 0).
-    // If unset, are recalculated when needed.
-    EdgeArray _edges;          // int array (?,2) of start & end point indices.
-    NeighborArray _neighbors;  // int array (ntri,3), neighbor triangle indices
-                               //     or -1 if no neighbor.
-
-    // Variables internal to C++ only.
-    Boundaries _boundaries;
-
-    // Map used to look up BoundaryEdges from TriEdges.  Normally accessed via
-    // get_boundary_edge().
-    typedef std::map<TriEdge, BoundaryEdge> TriEdgeToBoundaryMap;
-    TriEdgeToBoundaryMap _tri_edge_to_boundary_map;
-};
-
-
-
-// Contour generator for a triangulation.
-class TriContourGenerator
-{
-public:
-    typedef Triangulation::CoordinateArray CoordinateArray;
-    typedef Triangulation::TwoCoordinateArray TwoCoordinateArray;
-    typedef py::array_t<unsigned char> CodeArray;
-
-    /* Constructor.
-     *   triangulation: Triangulation to generate contours for.
-     *   z: Double array of shape (npoints) of z-values at triangulation
-     *      points. */
-    TriContourGenerator(Triangulation& triangulation,
-                        const CoordinateArray& z);
-
-    /* Create and return a non-filled contour.
-     *   level: Contour level.
-     * Returns new python list [segs0, segs1, ...] where
-     *   segs0: double array of shape (?,2) of point coordinates of first
-     *   contour line, etc. */
-    py::tuple create_contour(const double& level);
-
-    /* Create and return a filled contour.
-     *   lower_level: Lower contour level.
-     *   upper_level: Upper contour level.
-     * Returns new python tuple (segs, kinds) where
-     *   segs: double array of shape (n_points,2) of all point coordinates,
-     *   kinds: ubyte array of shape (n_points) of all point code types. */
-    py::tuple create_filled_contour(const double& lower_level,
-                                    const double& upper_level);
-
-private:
-    typedef Triangulation::Boundary Boundary;
-    typedef Triangulation::Boundaries Boundaries;
-
-    /* Clear visited flags.
-     *   include_boundaries: Whether to clear boundary flags or not, which are
-     *                       only used for filled contours. */
-    void clear_visited_flags(bool include_boundaries);
-
-    /* Convert a non-filled Contour from C++ to Python.
-     * Returns new python tuple ([segs0, segs1, ...], [kinds0, kinds1...])
-     * where
-     *   segs0: double array of shape (n_points,2) of point coordinates of first
-     *   contour line, etc.
-     *   kinds0: ubyte array of shape (n_points) of kinds codes of first contour
-     *   line, etc. */
-    py::tuple contour_line_to_segs_and_kinds(const Contour& contour);
-
-    /* Convert a filled Contour from C++ to Python.
-     * Returns new python tuple ([segs], [kinds]) where
-     *   segs: double array of shape (n_points,2) of all point coordinates,
-     *   kinds: ubyte array of shape (n_points) of all point code types. */
-    py::tuple contour_to_segs_and_kinds(const Contour& contour);
-
-    /* Return the point on the specified TriEdge that intersects the specified
-     * level. */
-    XY edge_interp(int tri, int edge, const double& level);
-
-    /* Find and follow non-filled contour lines that start and end on a
-     * boundary of the Triangulation.
-     *   contour: Contour to add new lines to.
-     *   level: Contour level. */
-    void find_boundary_lines(Contour& contour,
-                             const double& level);
-
-    /* Find and follow filled contour lines at either of the specified contour
-     * levels that start and end of a boundary of the Triangulation.
-     *   contour: Contour to add new lines to.
-     *   lower_level: Lower contour level.
-     *   upper_level: Upper contour level. */
-    void find_boundary_lines_filled(Contour& contour,
-                                    const double& lower_level,
-                                    const double& upper_level);
-
-    /* Find and follow lines at the specified contour level that are
-     * completely in the interior of the Triangulation and hence do not
-     * intersect any boundary.
-     *   contour: Contour to add new lines to.
-     *   level: Contour level.
-     *   on_upper: Whether on upper or lower contour level.
-     *   filled: Whether contours are filled or not. */
-    void find_interior_lines(Contour& contour,
-                             const double& level,
-                             bool on_upper,
-                             bool filled);
-
-    /* Follow contour line around boundary of the Triangulation from the
-     * specified TriEdge to its end which can be on either the lower or upper
-     * levels.  Only used for filled contours.
-     *   contour_line: Contour line to append new points to.
-     *   tri_edge: On entry, TriEdge to start from.  On exit, TriEdge that is
-     *             finished on.
-     *   lower_level: Lower contour level.
-     *   upper_level: Upper contour level.
-     *   on_upper: Whether starts on upper level or not.
-     * Return true if finishes on upper level, false if lower. */
-    bool follow_boundary(ContourLine& contour_line,
-                         TriEdge& tri_edge,
-                         const double& lower_level,
-                         const double& upper_level,
-                         bool on_upper);
-
-    /* Follow contour line across interior of Triangulation.
-     *   contour_line: Contour line to append new points to.
-     *   tri_edge: On entry, TriEdge to start from.  On exit, TriEdge that is
-     *             finished on.
-     *   end_on_boundary: Whether this line ends on a boundary, or loops back
-     *                    upon itself.
-     *   level: Contour level to follow.
-     *   on_upper: Whether following upper or lower contour level. */
-    void follow_interior(ContourLine& contour_line,
-                         TriEdge& tri_edge,
-                         bool end_on_boundary,
-                         const double& level,
-                         bool on_upper);
-
-    // Return the Triangulation boundaries.
-    const Boundaries& get_boundaries() const;
-
-    /* Return the edge by which the a level leaves a particular triangle,
-     * which is 0, 1 or 2 if the contour passes through the triangle or -1
-     * otherwise.
-     *   tri: Triangle index.
-     *   level: Contour level to follow.
-     *   on_upper: Whether following upper or lower contour level. */
-    int get_exit_edge(int tri, const double& level, bool on_upper) const;
-
-    // Return the z-value at the specified point index.
-    const double& get_z(int point) const;
-
-    /* Return the point at which the a level intersects the line connecting the
-     * two specified point indices. */
-    XY interp(int point1, int point2, const double& level) const;
-
-
-
-    // Variables shared with python, always set.
-    Triangulation _triangulation;
-    CoordinateArray _z;        // double array (npoints).
-
-    // Variables internal to C++ only.
-    typedef std::vector<bool> InteriorVisited;    // Size 2*ntri
-    typedef std::vector<bool> BoundaryVisited;
-    typedef std::vector<BoundaryVisited> BoundariesVisited;
-    typedef std::vector<bool> BoundariesUsed;
-
-    InteriorVisited _interior_visited;
-    BoundariesVisited _boundaries_visited;  // Only used for filled contours.
-    BoundariesUsed _boundaries_used;        // Only used for filled contours.
-};
-
-
-
-/* TriFinder class implemented using the trapezoid map algorithm from the book
- * "Computational Geometry, Algorithms and Applications", second edition, by
- * M. de Berg, M. van Kreveld, M. Overmars and O. Schwarzkopf.
- *
- * The domain of interest is composed of vertical-sided trapezoids that are
- * bounded to the left and right by points of the triangulation, and below and
- * above by edges of the triangulation.  Each triangle is represented by 1 or
- * more of these trapezoids.  Edges are inserted one a time in a random order.
- *
- * As the trapezoid map is created, a search tree is also created which allows
- * fast lookup O(log N) of the trapezoid containing the point of interest.
- * There are 3 types of node in the search tree: all leaf nodes represent
- * trapezoids and all branch nodes have 2 child nodes and are either x-nodes or
- * y-nodes.  X-nodes represent points in the triangulation, and their 2 children
- * refer to those parts of the search tree to the left and right of the point.
- * Y-nodes represent edges in the triangulation, and their 2 children refer to
- * those parts of the search tree below and above the edge.
- *
- * Nodes can be repeated throughout the search tree, and each is reference
- * counted through the multiple parent nodes it is a child of.
- *
- * The algorithm is only intended to work with valid triangulations, i.e. it
- * must not contain duplicate points, triangles formed from colinear points, or
- * overlapping triangles.  It does have some tolerance to triangles formed from
- * colinear points but only in the simplest of cases.  No explicit testing of
- * the validity of the triangulation is performed as this is a computationally
- * more complex task than the trifinding itself. */
-class TrapezoidMapTriFinder
-{
-public:
-    typedef Triangulation::CoordinateArray CoordinateArray;
-    typedef py::array_t<int, py::array::c_style | py::array::forcecast> TriIndexArray;
-
-    /* Constructor.  A separate call to initialize() is required to initialize
-     * the object before use.
-     *   triangulation: Triangulation to find triangles in. */
-    TrapezoidMapTriFinder(Triangulation& triangulation);
-
-    ~TrapezoidMapTriFinder();
-
-    /* Return an array of triangle indices.  Takes 1D arrays x and y of
-     * point coordinates, and returns an array of the same size containing the
-     * indices of the triangles at those points. */
-    TriIndexArray find_many(const CoordinateArray& x, const CoordinateArray& y);
-
-    /* Return a reference to a new python list containing the following
-     * statistics about the tree:
-     *   0: number of nodes (tree size)
-     *   1: number of unique nodes (number of unique Node objects in tree)
-     *   2: number of trapezoids (tree leaf nodes)
-     *   3: number of unique trapezoids
-     *   4: maximum parent count (max number of times a node is repeated in
-     *          tree)
-     *   5: maximum depth of tree (one more than the maximum number of
-     *          comparisons needed to search through the tree)
-     *   6: mean of all trapezoid depths (one more than the average number of
-     *          comparisons needed to search through the tree) */
-    py::list get_tree_stats();
-
-    /* Initialize this object before use.  May be called multiple times, if,
-     * for example, the triangulation is changed by setting the mask. */
-    void initialize();
-
-    // Print the search tree as text to stdout; useful for debug purposes.
-    void print_tree();
-
-private:
-    /* A Point consists of x,y coordinates as well as the index of a triangle
-     * associated with the point, so that a search at this point's coordinates
-     * can return a valid triangle index. */
-    struct Point : XY
-    {
-        Point() : XY(), tri(-1) {}
-        Point(const double& x, const double& y) : XY(x,y), tri(-1) {}
-        explicit Point(const XY& xy) : XY(xy), tri(-1) {}
-
-        int tri;
-    };
-
-    /* An Edge connects two Points, left and right.  It is always true that
-     * right->is_right_of(*left).  Stores indices of triangles below and above
-     * the Edge which are used to map from trapezoid to triangle index.  Also
-     * stores pointers to the 3rd points of the below and above triangles,
-     * which are only used to disambiguate triangles with colinear points. */
-    struct Edge
-    {
-        Edge(const Point* left_,
-             const Point* right_,
-             int triangle_below_,
-             int triangle_above_,
-             const Point* point_below_,
-             const Point* point_above_);
-
-        // Return -1 if point to left of edge, 0 if on edge, +1 if to right.
-        int get_point_orientation(const XY& xy) const;
-
-        // Return slope of edge, even if vertical (divide by zero is OK here).
-        double get_slope() const;
-
-        /* Return y-coordinate of point on edge with specified x-coordinate.
-         * x must be within the x-limits of this edge. */
-        double get_y_at_x(const double& x) const;
-
-        // Return true if the specified point is either of the edge end points.
-        bool has_point(const Point* point) const;
-
-        bool operator==(const Edge& other) const;
-
-        friend std::ostream& operator<<(std::ostream& os, const Edge& edge)
-        {
-            return os << *edge.left << "->" << *edge.right;
-        }
-
-        void print_debug() const;
-
-
-        const Point* left;        // Not owned.
-        const Point* right;       // Not owned.
-        int triangle_below;       // Index of triangle below (to right of) Edge.
-        int triangle_above;       // Index of triangle above (to left of) Edge.
-        const Point* point_below; // Used only for resolving ambiguous cases;
-        const Point* point_above; //     is 0 if corresponding triangle is -1
-    };
-
-    class Node;  // Forward declaration.
-
-    // Helper structure used by TrapezoidMapTriFinder::get_tree_stats.
-    struct NodeStats
-    {
-        NodeStats()
-            : node_count(0), trapezoid_count(0), max_parent_count(0),
-              max_depth(0), sum_trapezoid_depth(0.0)
-        {}
-
-        long node_count, trapezoid_count, max_parent_count, max_depth;
-        double sum_trapezoid_depth;
-        std::set<const Node*> unique_nodes, unique_trapezoid_nodes;
-    };
-
-    struct Trapezoid;  // Forward declaration.
-
-    /* Node of the trapezoid map search tree.  There are 3 possible types:
-     * Type_XNode, Type_YNode and Type_TrapezoidNode.  Data members are
-     * represented using a union: an XNode has a Point and 2 child nodes
-     * (left and right of the point), a YNode has an Edge and 2 child nodes
-     * (below and above the edge), and a TrapezoidNode has a Trapezoid.
-     * Each Node has multiple parents so it can appear in the search tree
-     * multiple times without having to create duplicate identical Nodes.
-     * The parent collection acts as a reference count to the number of times
-     * a Node occurs in the search tree.  When the parent count is reduced to
-     * zero a Node can be safely deleted. */
-    class Node
-    {
-    public:
-        Node(const Point* point, Node* left, Node* right);// Type_XNode.
-        Node(const Edge* edge, Node* below, Node* above); // Type_YNode.
-        Node(Trapezoid* trapezoid);                       // Type_TrapezoidNode.
-
-        ~Node();
-
-        void add_parent(Node* parent);
-
-        /* Recurse through the search tree and assert that everything is valid.
-         * Reduces to a no-op if NDEBUG is defined. */
-        void assert_valid(bool tree_complete) const;
-
-        // Recurse through the tree to return statistics about it.
-        void get_stats(int depth, NodeStats& stats) const;
-
-        // Return the index of the triangle corresponding to this node.
-        int get_tri() const;
-
-        bool has_child(const Node* child) const;
-        bool has_no_parents() const;
-        bool has_parent(const Node* parent) const;
-
-        /* Recurse through the tree and print a textual representation to
-         * stdout.  Argument depth used to indent for readability. */
-        void print(int depth = 0) const;
-
-        /* Remove a parent from this Node.  Return true if no parents remain
-         * so that this Node can be deleted. */
-        bool remove_parent(Node* parent);
-
-        void replace_child(Node* old_child, Node* new_child);
-
-        // Replace this node with the specified new_node in all parents.
-        void replace_with(Node* new_node);
-
-        /* Recursive search through the tree to find the Node containing the
-         * specified XY point. */
-        const Node* search(const XY& xy);
-
-        /* Recursive search through the tree to find the Trapezoid containing
-         * the left endpoint of the specified Edge.  Return 0 if fails, which
-         * can only happen if the triangulation is invalid. */
-        Trapezoid* search(const Edge& edge);
-
-        /* Copy constructor and assignment operator defined but not implemented
-         * to prevent objects being copied. */
-        Node(const Node& other);
-        Node& operator=(const Node& other);
-
-    private:
-        typedef enum {
-            Type_XNode,
-            Type_YNode,
-            Type_TrapezoidNode
-        } Type;
-        Type _type;
-
-        union {
-            struct {
-                const Point* point;  // Not owned.
-                Node* left;          // Owned.
-                Node* right;         // Owned.
-            } xnode;
-            struct {
-                const Edge* edge;    // Not owned.
-                Node* below;         // Owned.
-                Node* above;         // Owned.
-            } ynode;
-            Trapezoid* trapezoid;    // Owned.
-        } _union;
-
-        typedef std::list<Node*> Parents;
-        Parents _parents;            // Not owned.
-    };
-
-    /* A Trapezoid is bounded by Points to left and right, and Edges below and
-     * above.  Has up to 4 neighboring Trapezoids to lower/upper left/right.
-     * Lower left neighbor is Trapezoid to left that shares the below Edge, or
-     * is 0 if there is no such Trapezoid (and similar for other neighbors).
-     * To obtain the index of the triangle corresponding to a particular
-     * Trapezoid, use the Edge member variables below.triangle_above or
-     * above.triangle_below. */
-    struct Trapezoid
-    {
-        Trapezoid(const Point* left_,
-                  const Point* right_,
-                  const Edge& below_,
-                  const Edge& above_);
-
-        /* Assert that this Trapezoid is valid.  Reduces to a no-op if NDEBUG
-         * is defined. */
-        void assert_valid(bool tree_complete) const;
-
-        /* Return one of the 4 corner points of this Trapezoid.  Only used for
-         * debugging purposes. */
-        XY get_lower_left_point() const;
-        XY get_lower_right_point() const;
-        XY get_upper_left_point() const;
-        XY get_upper_right_point() const;
-
-        void print_debug() const;
-
-        /* Set one of the 4 neighbor trapezoids and the corresponding reverse
-         * Trapezoid of the new neighbor (if it is not 0), so that they are
-         * consistent. */
-        void set_lower_left(Trapezoid* lower_left_);
-        void set_lower_right(Trapezoid* lower_right_);
-        void set_upper_left(Trapezoid* upper_left_);
-        void set_upper_right(Trapezoid* upper_right_);
-
-        /* Copy constructor and assignment operator defined but not implemented
-         * to prevent objects being copied. */
-        Trapezoid(const Trapezoid& other);
-        Trapezoid& operator=(const Trapezoid& other);
-
-
-        const Point* left;     // Not owned.
-        const Point* right;    // Not owned.
-        const Edge& below;
-        const Edge& above;
-
-        // 4 neighboring trapezoids, can be 0, not owned.
-        Trapezoid* lower_left;   // Trapezoid to left  that shares below
-        Trapezoid* lower_right;  // Trapezoid to right that shares below
-        Trapezoid* upper_left;   // Trapezoid to left  that shares above
-        Trapezoid* upper_right;  // Trapezoid to right that shares above
-
-        Node* trapezoid_node;    // Node that owns this Trapezoid.
-    };
-
-
-    // Add the specified Edge to the search tree, returning true if successful.
-    bool add_edge_to_tree(const Edge& edge);
-
-    // Clear all memory allocated by this object.
-    void clear();
-
-    // Return the triangle index at the specified point, or -1 if no triangle.
-    int find_one(const XY& xy);
-
-    /* Determine the trapezoids that the specified Edge intersects, returning
-     * true if successful. */
-    bool find_trapezoids_intersecting_edge(const Edge& edge,
-                                           std::vector<Trapezoid*>& trapezoids);
-
-
-
-    // Variables shared with python, always set.
-    Triangulation& _triangulation;
-
-    // Variables internal to C++ only.
-    Point* _points;    // Array of all points in triangulation plus corners of
-                       // enclosing rectangle.  Owned.
-
-    typedef std::vector<Edge> Edges;
-    Edges _edges;   // All Edges in triangulation plus bottom and top Edges of
-                    // enclosing rectangle.
-
-    Node* _tree;    // Root node of the trapezoid map search tree.  Owned.
-};
-
-#endif
diff --git a/contrib/python/matplotlib/py3/src/tri/_tri_wrapper.cpp b/contrib/python/matplotlib/py3/src/tri/_tri_wrapper.cpp
deleted file mode 100644
index 1b0c3d7555..0000000000
--- a/contrib/python/matplotlib/py3/src/tri/_tri_wrapper.cpp
+++ /dev/null
@@ -1,58 +0,0 @@
-#include "_tri.h"
-
-PYBIND11_MODULE(_tri, m) {
-    py::class_<Triangulation>(m, "Triangulation")
-        .def(py::init<const Triangulation::CoordinateArray&,
-                      const Triangulation::CoordinateArray&,
-                      const Triangulation::TriangleArray&,
-                      const Triangulation::MaskArray&,
-                      const Triangulation::EdgeArray&,
-                      const Triangulation::NeighborArray&,
-                      bool>(),
-            py::arg("x"),
-            py::arg("y"),
-            py::arg("triangles"),
-            py::arg("mask"),
-            py::arg("edges"),
-            py::arg("neighbors"),
-            py::arg("correct_triangle_orientations"),
-            "Create a new C++ Triangulation object.\n"
-            "This should not be called directly, use the python class\n"
-            "matplotlib.tri.Triangulation instead.\n")
-        .def("calculate_plane_coefficients", &Triangulation::calculate_plane_coefficients,
-            "Calculate plane equation coefficients for all unmasked triangles.")
-        .def("get_edges", &Triangulation::get_edges,
-            "Return edges array.")
-        .def("get_neighbors", &Triangulation::get_neighbors,
-            "Return neighbors array.")
-        .def("set_mask", &Triangulation::set_mask,
-            "Set or clear the mask array.");
-
-    py::class_<TriContourGenerator>(m, "TriContourGenerator")
-        .def(py::init<Triangulation&,
-                      const TriContourGenerator::CoordinateArray&>(),
-            py::arg("triangulation"),
-            py::arg("z"),
-            "Create a new C++ TriContourGenerator object.\n"
-            "This should not be called directly, use the functions\n"
-            "matplotlib.axes.tricontour and tricontourf instead.\n")
-        .def("create_contour", &TriContourGenerator::create_contour,
-            "Create and return a non-filled contour.")
-        .def("create_filled_contour", &TriContourGenerator::create_filled_contour,
-            "Create and return a filled contour.");
-
-    py::class_<TrapezoidMapTriFinder>(m, "TrapezoidMapTriFinder")
-        .def(py::init<Triangulation&>(),
-            py::arg("triangulation"),
-            "Create a new C++ TrapezoidMapTriFinder object.\n"
-            "This should not be called directly, use the python class\n"
-            "matplotlib.tri.TrapezoidMapTriFinder instead.\n")
-        .def("find_many", &TrapezoidMapTriFinder::find_many,
-            "Find indices of triangles containing the point coordinates (x, y).")
-        .def("get_tree_stats", &TrapezoidMapTriFinder::get_tree_stats,
-            "Return statistics about the tree used by the trapezoid map.")
-        .def("initialize", &TrapezoidMapTriFinder::initialize,
-            "Initialize this object, creating the trapezoid map from the triangulation.")
-        .def("print_tree", &TrapezoidMapTriFinder::print_tree,
-            "Print the search tree as text to stdout; useful for debug purposes.");
-}