Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Added TriFinder and TrapezoidMapTriFinder classes.

  • Loading branch information...
commit b903eecc73b7d7fa08338a106b663ac9af2adbfb 1 parent 1787a1a
@ianthomas23 ianthomas23 authored
View
9 doc/api/tri_api.rst
@@ -5,6 +5,11 @@ triangular grids
:mod:`matplotlib.tri`
=====================
-.. automodule:: matplotlib.tri
- :members: Triangulation
+.. autoclass:: matplotlib.tri.Triangulation
+ :members:
+.. autoclass:: matplotlib.tri.TriFinder
+ :members:
+
+.. autoclass:: matplotlib.tri.TrapezoidMapTriFinder
+ :members: __call__
View
59 examples/event_handling/trifinder_event_demo.py
@@ -0,0 +1,59 @@
+"""
+Example showing the use of a TriFinder object. As the mouse is moved over the
+triangulation, the triangle under the cursor is highlighted and the index of
+the triangle is displayed in the plot title.
+"""
+import matplotlib.pyplot as plt
+from matplotlib.tri import Triangulation
+from matplotlib.patches import Polygon
+import numpy as np
+import math
+
+
+def update_polygon(tri):
+ if tri == -1:
+ points = [0, 0, 0]
+ else:
+ points = triangulation.triangles[tri]
+ xs = triangulation.x[points]
+ ys = triangulation.y[points]
+ polygon.set_xy(zip(xs, ys))
+
+
+def motion_notify(event):
+ if event.inaxes is None:
+ tri = -1
+ else:
+ tri = trifinder(event.xdata, event.ydata)
+ update_polygon(tri)
+ plt.title('In triangle %i' % tri)
+ event.canvas.draw()
+
+
+# Create a Triangulation.
+n_angles = 16
+n_radii = 5
+min_radius = 0.25
+radii = np.linspace(min_radius, 0.95, n_radii)
+angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
+angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
+angles[:, 1::2] += math.pi / n_angles
+x = (radii*np.cos(angles)).flatten()
+y = (radii*np.sin(angles)).flatten()
+triangulation = Triangulation(x, y)
+xmid = x[triangulation.triangles].mean(axis=1)
+ymid = y[triangulation.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triangulation.set_mask(mask)
+
+# Use the triangulation's default TriFinder object.
+trifinder = triangulation.get_trifinder()
+
+# Setup plot and callbacks.
+plt.subplot(111, aspect='equal')
+plt.triplot(triangulation, 'bo-')
+polygon = Polygon([[0, 0], [0, 0]], facecolor='y') # dummy data for xs,ys
+update_polygon(-1)
+plt.gca().add_patch(polygon)
+plt.gcf().canvas.mpl_connect('motion_notify_event', motion_notify)
+plt.show()
View
88 lib/matplotlib/tests/test_triangulation.py
@@ -131,3 +131,91 @@ def test_no_modify():
tri = mtri.Triangulation(points[:,0], points[:,1], triangles)
edges = tri.edges
assert_array_equal(old_triangles, triangles)
+
+def test_trifinder():
+ # Test points within triangles of masked triangulation.
+ x,y = np.meshgrid(np.arange(4), np.arange(4))
+ x = x.ravel()
+ y = y.ravel()
+ triangles = [[0,1,4], [1,5,4], [1,2,5], [2,6,5], [2,3,6], [3,7,6], [4,5,8],
+ [5,9,8], [5,6,9], [6,10,9], [6,7,10], [7,11,10], [8,9,12],
+ [9,13,12], [9,10,13], [10,14,13], [10,11,14], [11,15,14]]
+ mask = np.zeros(len(triangles))
+ mask[8:10] = 1
+ triang = mtri.Triangulation(x, y, triangles, mask)
+ trifinder = triang.get_trifinder()
+
+ xs = [0.25, 1.25, 2.25, 3.25]
+ ys = [0.25, 1.25, 2.25, 3.25]
+ xs,ys = np.meshgrid(xs,ys)
+ xs = xs.ravel()
+ ys = ys.ravel()
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [0,2,4,-1,6,-1,10,-1,12,14,16,-1,-1,-1,-1,-1])
+ tris = trifinder(xs-0.5, ys-0.5)
+ assert_array_equal(tris, [-1,-1,-1,-1,-1,1,3,5,-1,7,-1,11,-1,13,15,17])
+
+ # Test points exactly on boundary edges of masked triangulation.
+ xs = [0.5, 1.5, 2.5, 0.5, 1.5, 2.5, 1.5, 1.5, 0.0, 1.0, 2.0, 3.0]
+ ys = [0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 1.0, 2.0, 1.5, 1.5, 1.5, 1.5]
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [0,2,4,13,15,17,3,14,6,7,10,11])
+
+ # Test points exactly on boundary corners of masked triangulation.
+ xs = [0.0, 3.0]
+ ys = [0.0, 3.0]
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [0,17])
+
+ # Test triangles with horizontal colinear points. These are not valid
+ # triangulations, but we try to deal with the simplest violations.
+ delta = 0.0 # If +ve, triangulation is OK, if -ve triangulation invalid,
+ # if zero have colinear points but should pass tests anyway.
+ x = [1.5, 0, 1, 2, 3, 1.5, 1.5]
+ y = [-1, 0, 0, 0, 0, delta, 1]
+ triangles = [[0,2,1], [0,3,2], [0,4,3], [1,2,5], [2,3,5], [3,4,5], [1,5,6],
+ [4,6,5]]
+ triang = mtri.Triangulation(x, y, triangles)
+ trifinder = triang.get_trifinder()
+
+ xs = [-0.1, 0.4, 0.9, 1.4, 1.9, 2.4, 2.9]
+ ys = [-0.1,0.1]
+ xs,ys = np.meshgrid(xs, ys)
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [[-1,0,0,1,1,2,-1],[-1,6,6,6,7,7,-1]])
+
+ # Test triangles with vertical colinear points. These are not valid
+ # triangulations, but we try to deal with the simplest violations.
+ delta = 0.0 # If +ve, triangulation is OK, if -ve triangulation invalid,
+ # if zero have colinear points but should pass tests anyway.
+ x = [-1, -delta, 0, 0, 0, 0, 1]
+ y = [1.5, 1.5, 0, 1, 2, 3, 1.5]
+ triangles = [[0,1,2], [0,1,5], [1,2,3], [1,3,4], [1,4,5], [2,6,3], [3,6,4],
+ [4,6,5]]
+ triang = mtri.Triangulation(x, y, triangles)
+ trifinder = triang.get_trifinder()
+
+ xs = [-0.1,0.1]
+ ys = [-0.1, 0.4, 0.9, 1.4, 1.9, 2.4, 2.9]
+ xs,ys = np.meshgrid(xs, ys)
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [[-1,-1], [0,5], [0,5], [0,6], [1,6], [1,7],
+ [-1,-1]])
+
+ # Test that changing triangulation by setting a mask causes the trifinder
+ # to be reinitialised.
+ x = [0, 1, 0, 1]
+ y = [0, 0, 1, 1]
+ triangles = [[0,1,2], [1,3,2]]
+ triang = mtri.Triangulation(x, y, triangles)
+ trifinder = triang.get_trifinder()
+
+ xs = [-0.2, 0.2, 0.8, 1.2]
+ ys = [0.5, 0.5, 0.5, 0.5]
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [-1, 0, 1, -1])
+
+ triang.set_mask([1, 0])
+ assert_equal(trifinder, triang.get_trifinder())
+ tris = trifinder(xs, ys)
+ assert_array_equal(tris, [-1, -1, 1, -1])
View
1  lib/matplotlib/tri/__init__.py
@@ -5,5 +5,6 @@
from __future__ import print_function
from triangulation import *
from tricontour import *
+from trifinder import *
from tripcolor import *
from triplot import *
View
1,189 lib/matplotlib/tri/_tri.cpp
@@ -1,3 +1,10 @@
+/* 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 "_tri.h"
#include "src/mplutils.h"
@@ -60,6 +67,14 @@ 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;
@@ -70,17 +85,23 @@ bool XY::operator!=(const XY& other) const
return x != other.x || y != other.y;
}
-bool XY::operator<(const XY& other) const
+XY XY::operator*(const double& multiplier) const
{
- if (y == other.y)
- return x < other.x;
- else
- return y < other.y;
+ return XY(x*multiplier, y*multiplier);
}
-XY XY::operator*(const double& multiplier) const
+const XY& XY::operator+=(const XY& other)
{
- return XY(x*multiplier, y*multiplier);
+ 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
@@ -100,6 +121,34 @@ std::ostream& operator<<(std::ostream& os, const XY& xy)
+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)
{}
@@ -118,17 +167,13 @@ void BoundingBox::add(const XY& point)
}
}
-bool BoundingBox::contains_x(const double& x) const
-{
- return !empty && x >= lower.x && x <= upper.x;
-}
-
-std::ostream& operator<<(std::ostream& os, const BoundingBox& box)
+void BoundingBox::expand(const XY& delta)
{
- if (box.empty)
- return os << "<empty>";
- else
- return os << box.lower << " -> " << box.upper;
+ if (!empty)
+ {
+ lower -= delta;
+ upper += delta;
+ }
}
@@ -475,6 +520,15 @@ bool Triangulation::is_masked(int tri) const
return _mask && *((bool*)PyArray_DATA(_mask) + tri);
}
+bool Triangulation::is_triangle_point(int tri, int point) const
+{
+ assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds.");
+ assert(point >= 0 && point < _npoints && "Point index out of bounds.");
+ return (get_triangle_point(tri,0) == point ||
+ get_triangle_point(tri,1) == point ||
+ get_triangle_point(tri,2) == point);
+}
+
Py::Object Triangulation::set_mask(const Py::Tuple &args)
{
_VERBOSE("Triangulation::set_mask");
@@ -980,6 +1034,1091 @@ XY TriContourGenerator::interp(int point1,
+TrapezoidMapTriFinder::TrapezoidMapTriFinder(Py::Object triangulation)
+ : _triangulation(triangulation),
+ _points(0),
+ _tree(0)
+{
+ _VERBOSE("TrapezoidMapTriFinder::TrapezoidMapTriFinder");
+}
+
+TrapezoidMapTriFinder::~TrapezoidMapTriFinder()
+{
+ _VERBOSE("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.
+ unsigned int ntraps = trapezoids.size();
+ for (unsigned int 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;
+}
+
+Py::Object TrapezoidMapTriFinder::find_many(const Py::Tuple& args)
+{
+ args.verify_length(2);
+
+ // Check input arguments.
+ PyArrayObject* x = (PyArrayObject*)PyArray_ContiguousFromObject(
+ args[0].ptr(), PyArray_DOUBLE, 0, 0);
+ PyArrayObject* y = (PyArrayObject*)PyArray_ContiguousFromObject(
+ args[1].ptr(), PyArray_DOUBLE, 0, 0);
+ bool ok = (x != 0 && y != 0 && PyArray_NDIM(x) == PyArray_NDIM(y));
+ int ndim = PyArray_NDIM(x);
+ for (int i = 0; ok && i < ndim; ++i)
+ ok = (PyArray_DIM(x,i) == PyArray_DIM(y,i));
+
+ if (!ok)
+ {
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ throw Py::ValueError("x and y must be array_like with same shape");
+ }
+
+ // Create integer array to return.
+ PyArrayObject* tri = (PyArrayObject*)PyArray_SimpleNew(
+ ndim, PyArray_DIMS(x), PyArray_INT);
+
+ // Fill returned array.
+ double* x_ptr = (double*)PyArray_DATA(x);
+ double* y_ptr = (double*)PyArray_DATA(y);
+ int* tri_ptr = (int*)PyArray_DATA(tri);
+ int* tri_end = tri_ptr + PyArray_SIZE(tri);
+ while (tri_ptr < tri_end)
+ *tri_ptr++ = find_one(XY(*x_ptr++, *y_ptr++));
+
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ return Py::asObject((PyObject*)tri);
+}
+
+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::Object TrapezoidMapTriFinder::get_tree_stats()
+{
+ _VERBOSE("TrapezoidMapTriFinder::get_tree_stats");
+
+ NodeStats stats;
+ _tree->get_stats(0, stats);
+
+ Py::List list(7);
+ list[0] = Py::Int(stats.node_count);
+ list[1] = Py::Int(static_cast<long>(stats.unique_nodes.size()));
+ list[2] = Py::Int(stats.trapezoid_count);
+ list[3] = Py::Int(static_cast<long>(stats.unique_trapezoid_nodes.size()));
+ list[4] = Py::Int(stats.max_parent_count);
+ list[5] = Py::Int(stats.max_depth);
+ list[6] = Py::Float(stats.sum_trapezoid_depth / stats.trapezoid_count);
+ return list;
+}
+
+const Triangulation& TrapezoidMapTriFinder::get_triangulation() const
+{
+ return *(Triangulation*)_triangulation.ptr();
+}
+
+void TrapezoidMapTriFinder::init_type()
+{
+ _VERBOSE("TrapezoidMapTriFinder::init_type");
+
+ behaviors().name("TrapezoidMapTriFinder");
+ behaviors().doc("TrapezoidMapTriFinder");
+
+ add_varargs_method("find_many",
+ &TrapezoidMapTriFinder::find_many,
+ "find_many(x,y)");
+ add_noargs_method("get_tree_stats",
+ &TrapezoidMapTriFinder::get_tree_stats,
+ "get_tree_stats()");
+ add_noargs_method("initialize",
+ &TrapezoidMapTriFinder::initialize,
+ "initialize()");
+ add_noargs_method("print_tree",
+ &TrapezoidMapTriFinder::print_tree,
+ "print_tree()");
+}
+
+Py::Object TrapezoidMapTriFinder::initialize()
+{
+ _VERBOSE("TrapezoidMapTriFinder::initialize");
+
+ clear();
+ const Triangulation& triang = get_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.
+ RandomNumberGenerator rng(1234);
+ std::random_shuffle(_edges.begin()+2, _edges.end(), rng);
+
+ // Add edges, one at a time, to tree.
+ unsigned int nedges = _edges.size();
+ for (unsigned int index = 2; index < nedges; ++index)
+ {
+ if (!add_edge_to_tree(_edges[index]))
+ throw Py::RuntimeError("Triangulation is invalid");
+ _tree->assert_valid(index == nedges-1);
+ }
+
+ return Py::None();
+}
+
+Py::Object TrapezoidMapTriFinder::print_tree()
+{
+ _VERBOSE("TrapezoidMapTriFinder::print_tree");
+
+ assert(_tree != 0 && "Null Node tree");
+ _tree->print();
+
+ return Py::None();
+}
+
+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;
+}
+
+
+RandomNumberGenerator::RandomNumberGenerator(unsigned long seed)
+ : _M(21870), _A(1291), _C(4621), _seed(seed % _M)
+{}
+
+unsigned long RandomNumberGenerator::operator()(unsigned long max_value)
+{
+ _seed = (_seed*_A + _C) % _M;
+ return (_seed*max_value) / _M;
+}
+
+
+
+
+
#if PY_MAJOR_VERSION >= 3
PyMODINIT_FUNC
PyInit__tri(void)
@@ -1003,11 +2142,15 @@ TriModule::TriModule()
{
Triangulation::init_type();
TriContourGenerator::init_type();
+ TrapezoidMapTriFinder::init_type();
add_varargs_method("Triangulation", &TriModule::new_triangulation,
"Create and return new C++ Triangulation object");
add_varargs_method("TriContourGenerator", &TriModule::new_tricontourgenerator,
"Create and return new C++ TriContourGenerator object");
+ add_varargs_method("TrapezoidMapTriFinder",
+ &TriModule::new_TrapezoidMapTriFinder,
+ "Create and return new C++ TrapezoidMapTriFinder object");
initialize("Module for unstructured triangular grids");
}
@@ -1113,3 +2256,15 @@ Py::Object TriModule::new_tricontourgenerator(const Py::Tuple &args)
return Py::asObject(new TriContourGenerator(tri, z));
}
+
+Py::Object TriModule::new_TrapezoidMapTriFinder(const Py::Tuple &args)
+{
+ _VERBOSE("TriModule::new_TrapezoidMapTriFinder");
+ args.verify_length(1);
+
+ Py::Object triangulation = args[0];
+ if (!Triangulation::check(triangulation))
+ throw Py::ValueError("Expecting a C++ Triangulation object");
+
+ return Py::asObject(new TrapezoidMapTriFinder(triangulation));
+}
View
368 lib/matplotlib/tri/_tri.h
@@ -67,7 +67,9 @@
#include "CXX/Objects.hxx"
#include "numpy/arrayobject.h"
+#include <list>
#include <map>
+#include <set>
#include <vector>
@@ -94,10 +96,12 @@ struct 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;
- 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);
@@ -105,14 +109,25 @@ struct 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);
- bool contains_x(const double& x) const;
- friend std::ostream& operator<<(std::ostream& os, const BoundingBox& box);
+ void expand(const XY& delta);
// Consider these member variables read-only.
bool empty;
@@ -216,6 +231,10 @@ class Triangulation : public Py::PythonExtension<Triangulation>
// Indicates if the specified triangle is masked or not.
bool is_masked(int tri) const;
+ /* Indicates if the specified point is one of the vertices of the specified
+ * triangle. */
+ bool is_triangle_point(int tri, int point) const;
+
/* Set or clear the mask array. Clears various derived fields so they are
* recalculated when next needed.
* args[0]: mask, either Py::None or boolean array of shape (ntri). */
@@ -459,6 +478,348 @@ class TriContourGenerator : public Py::PythonExtension<TriContourGenerator>
+/* 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 Py::PythonExtension<TrapezoidMapTriFinder>
+{
+public:
+ /* Constructor. A separate call to initialize() is required to initialize
+ * the object before use.
+ * triangulation: Triangulation to find triangles in. */
+ TrapezoidMapTriFinder(Py::Object triangulation);
+
+ virtual ~TrapezoidMapTriFinder();
+
+ /* Return an array of triangle indices. Takes any-shaped arrays x and y of
+ * point coordinates, and returns an array of the same shape containing the
+ * indices of the triangles at those points. */
+ Py::Object find_many(const Py::Tuple& args);
+
+ /* Return a 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::Object get_tree_stats();
+
+ // CXX initialisation function.
+ static void init_type();
+
+ /* Initialize this object before use. May be called multiple times, if,
+ * for example, the triangulation is changed by setting the mask. */
+ Py::Object initialize();
+
+ // Print the search tree as text to stdout; useful for debug purposes.
+ Py::Object 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);
+
+ // Return the underlying C++ Triangulation object.
+ const Triangulation& get_triangulation() const;
+
+
+
+ // Variables shared with python, always set.
+ Py::Object _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.
+};
+
+
+/* Linear congruential random number generator. Edges in the triangulation are
+ * randomly shuffled before being added to the trapezoid map. Want the
+ * shuffling to be identical across different operating systems and the same
+ * regardless of previous random number use. Would prefer to use a STL or
+ * Boost random number generator, but support is not consistent across
+ * different operating systems so implementing own here.
+ *
+ * This is not particularly random, but is perfectly adequate for the use here.
+ * Coefficients taken from Numerical Recipes in C. */
+class RandomNumberGenerator
+{
+public:
+ RandomNumberGenerator(unsigned long seed);
+
+ // Return random integer in the range 0 to max_value-1.
+ unsigned long operator()(unsigned long max_value);
+
+private:
+ const unsigned long _M, _A, _C;
+ unsigned long _seed;
+};
+
+
+
// The extension module.
class TriModule : public Py::ExtensionModule<TriModule>
{
@@ -468,6 +829,7 @@ class TriModule : public Py::ExtensionModule<TriModule>
private:
Py::Object new_triangulation(const Py::Tuple &args);
Py::Object new_tricontourgenerator(const Py::Tuple &args);
+ Py::Object new_TrapezoidMapTriFinder(const Py::Tuple &args);
};
#endif
View
35 lib/matplotlib/tri/triangulation.py
@@ -3,6 +3,7 @@
import matplotlib._tri as _tri
import numpy as np
+
class Triangulation(object):
"""
An unstructured triangular grid consisting of npoints points and
@@ -35,6 +36,9 @@ class Triangulation(object):
triangle. neighbors[i,j] is the triangle that is the neighbor
to the edge from point index triangles[i,j] to point index
triangles[i,(j+1)%3].
+
+ For a Triangulation to be valid it must not have duplicate points,
+ triangles formed from colinear points, or overlapping triangles.
"""
def __init__(self, x, y, triangles=None, mask=None):
self.x = np.asarray(x, dtype=np.float64)
@@ -49,11 +53,13 @@ def __init__(self, x, y, triangles=None, mask=None):
if triangles is None:
# No triangulation specified, so use matplotlib.delaunay.
dt = delaunay.Triangulation(self.x, self.y)
- self.triangles = np.asarray(dt.to_client_point_indices(dt.triangle_nodes),
- dtype=np.int32)
+ self.triangles = np.asarray(
+ dt.to_client_point_indices(dt.triangle_nodes),
+ dtype=np.int32)
if mask is None:
- self._edges = np.asarray(dt.to_client_point_indices(dt.edge_db),
- dtype=np.int32)
+ self._edges = np.asarray(
+ dt.to_client_point_indices(dt.edge_db),
+ dtype=np.int32)
# Delaunay triangle_neighbors uses different edge indexing,
# so convert.
neighbors = np.asarray(dt.triangle_neighbors, dtype=np.int32)
@@ -79,6 +85,9 @@ def __init__(self, x, y, triangles=None, mask=None):
# Underlying C++ object is not created until first needed.
self._cpp_triangulation = None
+ # Default TriFinder not created until needed.
+ self._trifinder = None
+
@property
def edges(self):
if self._edges is None:
@@ -99,7 +108,7 @@ def get_masked_triangles(self):
Return an array of triangles that are not masked.
"""
if self.mask is not None:
- return self.triangles.compress(1-self.mask, axis=0)
+ return self.triangles.compress(1 - self.mask, axis=0)
else:
return self.triangles
@@ -149,6 +158,18 @@ def get_from_args_and_kwargs(*args, **kwargs):
triangulation = Triangulation(x, y, triangles, mask)
return triangulation, args, kwargs
+ def get_trifinder(self):
+ """
+ Return the default :class:`matplotlib.tri.TriFinder` of this
+ triangulation, creating it if necessary. This allows the same
+ TriFinder object to be easily shared.
+ """
+ if self._trifinder is None:
+ # Default TriFinder class.
+ from matplotlib.tri.trifinder import TrapezoidMapTriFinder
+ self._trifinder = TrapezoidMapTriFinder(self)
+ return self._trifinder
+
@property
def neighbors(self):
if self._neighbors is None:
@@ -176,3 +197,7 @@ def set_mask(self, mask):
# Clear derived fields so they are recalculated when needed.
self._edges = None
self._neighbors = None
+
+ # Recalculate TriFinder if it exists.
+ if self._trifinder is not None:
+ self._trifinder._initialize()
View
84 lib/matplotlib/tri/trifinder.py
@@ -0,0 +1,84 @@
+from __future__ import print_function
+from matplotlib.tri import Triangulation
+import matplotlib._tri as _tri
+
+
+class TriFinder(object):
+ """
+ Abstract base class for classes used to find the triangles of a
+ Triangulation in which (x,y) points lie.
+
+ Rather than instantiate an object of a class derived from TriFinder, it is
+ usually better to use the function
+ :func:`matplotlib.tri.Triangulation.get_trifinder`.
+
+ Derived classes implement __call__(x,y) where x,y are array_like point
+ coordinates of the same shape.
+ """
+ def __init__(self, triangulation):
+ if not isinstance(triangulation, Triangulation):
+ raise ValueError('Expected a Triangulation object')
+ self._triangulation = triangulation
+
+
+class TrapezoidMapTriFinder(TriFinder):
+ """
+ :class:`~matplotlib.tri.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 triangulation must be valid, i.e. it must not have duplicate points,
+ triangles formed from colinear points, or overlapping triangles. The
+ algorithm has some tolerance to triangles formed from colinear points, but
+ this should not be relied upon.
+ """
+ def __init__(self, triangulation):
+ TriFinder.__init__(self, triangulation)
+ self._cpp_trifinder = _tri.TrapezoidMapTriFinder(
+ triangulation.get_cpp_triangulation())
+ self._initialize()
+
+ def __call__(self, x, y):
+ """
+ Return an array containing the indices of the triangles in which the
+ specified x,y points lie, or -1 for points that do not lie within a
+ triangle.
+
+ *x*, *y* are array_like x and y coordinates of the same shape and any
+ number of dimensions.
+
+ Returns integer array with the same shape and *x* and *y*.
+ """
+ # C++ checks arguments are OK.
+ return self._cpp_trifinder.find_many(x, y)
+
+ def _get_tree_stats(self):
+ """
+ Return a python list containing the statistics about the node tree:
+ 0: number of nodes (tree size)
+ 1: number of unique nodes
+ 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)
+ """
+ return self._cpp_trifinder.get_tree_stats()
+
+ def _initialize(self):
+ """
+ Initialize the underlying C++ object. Can be called multiple times if,
+ for example, the triangulation is modified.
+ """
+ self._cpp_trifinder.initialize()
+
+ def _print_tree(self):
+ """
+ Print a text representation of the node tree, which is useful for
+ debugging purposes.
+ """
+ self._cpp_trifinder.print_tree()

0 comments on commit b903eec

Please sign in to comment.
Something went wrong with that request. Please try again.