From b903eecc73b7d7fa08338a106b663ac9af2adbfb Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Wed, 5 Dec 2012 17:15:08 +0000 Subject: [PATCH 1/5] Added TriFinder and TrapezoidMapTriFinder classes. --- doc/api/tri_api.rst | 9 +- .../event_handling/trifinder_event_demo.py | 59 + lib/matplotlib/tests/test_triangulation.py | 88 ++ lib/matplotlib/tri/__init__.py | 1 + lib/matplotlib/tri/_tri.cpp | 1189 ++++++++++++++++- lib/matplotlib/tri/_tri.h | 368 ++++- lib/matplotlib/tri/triangulation.py | 35 +- lib/matplotlib/tri/trifinder.py | 84 ++ 8 files changed, 1806 insertions(+), 27 deletions(-) create mode 100644 examples/event_handling/trifinder_event_demo.py create mode 100644 lib/matplotlib/tri/trifinder.py diff --git a/doc/api/tri_api.rst b/doc/api/tri_api.rst index 1ecc8446c11e..9e0167aec856 100644 --- a/doc/api/tri_api.rst +++ b/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__ diff --git a/examples/event_handling/trifinder_event_demo.py b/examples/event_handling/trifinder_event_demo.py new file mode 100644 index 000000000000..fd7ccc77d400 --- /dev/null +++ b/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() diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py index 74fb37bdb62b..0239bd379417 100644 --- a/lib/matplotlib/tests/test_triangulation.py +++ b/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]) diff --git a/lib/matplotlib/tri/__init__.py b/lib/matplotlib/tri/__init__.py index f7fd7ca8d0f5..039ca110091e 100644 --- a/lib/matplotlib/tri/__init__.py +++ b/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 * diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp index 40fdb97380e0..f9deef5b8d2c 100644 --- a/lib/matplotlib/tri/_tri.cpp +++ b/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 << ""; - 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 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& 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(stats.unique_nodes.size())); + list[2] = Py::Int(stats.trapezoid_count); + list[3] = Py::Int(static_cast(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(_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)); +} diff --git a/lib/matplotlib/tri/_tri.h b/lib/matplotlib/tri/_tri.h index 842e6d39bf1f..1e6a95556998 100644 --- a/lib/matplotlib/tri/_tri.h +++ b/lib/matplotlib/tri/_tri.h @@ -67,7 +67,9 @@ #include "CXX/Objects.hxx" #include "numpy/arrayobject.h" +#include #include +#include #include @@ -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 // 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 +/* 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 +{ +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 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 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& 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 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 { @@ -468,6 +829,7 @@ class TriModule : public Py::ExtensionModule 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 diff --git a/lib/matplotlib/tri/triangulation.py b/lib/matplotlib/tri/triangulation.py index 122054b930ab..fb56d4d9bd76 100644 --- a/lib/matplotlib/tri/triangulation.py +++ b/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() diff --git a/lib/matplotlib/tri/trifinder.py b/lib/matplotlib/tri/trifinder.py new file mode 100644 index 000000000000..3d3427507267 --- /dev/null +++ b/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() From 413e18e371bfd910899d9269e8d8ed22b728ef6f Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Fri, 7 Dec 2012 07:53:23 +0000 Subject: [PATCH 2/5] Added LinearTriInterpolator class. --- CHANGELOG | 4 + doc/api/tri_api.rst | 6 ++ doc/users/whats_new.rst | 11 ++- examples/pylab_examples/triinterp_demo.py | 36 ++++++++ lib/matplotlib/tests/test_triangulation.py | 31 +++++++ lib/matplotlib/tri/__init__.py | 1 + lib/matplotlib/tri/_tri.cpp | 97 ++++++++++++++++++++ lib/matplotlib/tri/_tri.h | 8 ++ lib/matplotlib/tri/triangulation.py | 10 ++ lib/matplotlib/tri/triinterpolate.py | 101 +++++++++++++++++++++ 10 files changed, 303 insertions(+), 2 deletions(-) create mode 100644 examples/pylab_examples/triinterp_demo.py create mode 100644 lib/matplotlib/tri/triinterpolate.py diff --git a/CHANGELOG b/CHANGELOG index 352c176fa11b..df9bb8e4e2b6 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,7 @@ +2012-12-22 Added classes for interpolation within triangular grids + (LinearTriInterpolator) and to find the triangles in which points + lie (TrapezoidMapTriFinder) to matplotlib.tri module. - IMT + 2012-12-05 Added MatplotlibDeprecationWarning class for signaling deprecation. Matplotlib developers can use this class as follows: diff --git a/doc/api/tri_api.rst b/doc/api/tri_api.rst index 9e0167aec856..0a91c2549d8d 100644 --- a/doc/api/tri_api.rst +++ b/doc/api/tri_api.rst @@ -13,3 +13,9 @@ triangular grids .. autoclass:: matplotlib.tri.TrapezoidMapTriFinder :members: __call__ + +.. autoclass:: matplotlib.tri.TriInterpolator + :members: + +.. autoclass:: matplotlib.tri.LinearTriInterpolator + :members: __call__ diff --git a/doc/users/whats_new.rst b/doc/users/whats_new.rst index 5ba32223fae4..7171325099a1 100644 --- a/doc/users/whats_new.rst +++ b/doc/users/whats_new.rst @@ -49,6 +49,13 @@ using Inkscape for example, while preserving their intended position. For `svg` please note that you'll have to disable the default text-to-path conversion (`mpl.rc('svg', fonttype='none')`). +Triangular grid interpolation +----------------------------- +Ian Thomas added classes to perform interpolation within triangular grids +(:class:`~matplotlib.tri.LinearTriInterpolator`) and a utility class to find +the triangles in which points lie ( +:class:`~matplotlib.tri.TrapezoidMapTriFinder`). + .. _whats-new-1-2: new in matplotlib-1.2 @@ -298,7 +305,7 @@ to address the most common layout issues. fig, axes_list = plt.subplots(2, 1) for ax in axes_list.flat: ax.set(xlabel="x-label", ylabel="y-label", title="before tight_layout") - ax.locator_params(nbins=3) + ax.locator_params(nbins=3) plt.show() @@ -308,7 +315,7 @@ to address the most common layout issues. fig, axes_list = plt.subplots(2, 1) for ax in axes_list.flat: ax.set(xlabel="x-label", ylabel="y-label", title="after tight_layout") - ax.locator_params(nbins=3) + ax.locator_params(nbins=3) plt.tight_layout() plt.show() diff --git a/examples/pylab_examples/triinterp_demo.py b/examples/pylab_examples/triinterp_demo.py new file mode 100644 index 000000000000..6b5d0efbdf52 --- /dev/null +++ b/examples/pylab_examples/triinterp_demo.py @@ -0,0 +1,36 @@ +""" +Interpolation from triangular grid to quad grid. +""" +import matplotlib.pyplot as plt +import matplotlib.tri as mtri +import numpy as np +import math + + +# Create triangulation. +x = np.asarray([0, 1, 2, 3, 0.5, 1.5, 2.5, 1, 2, 1.5]) +y = np.asarray([0, 0, 0, 0, 1, 1, 1, 2, 2, 3]) +triangles = [[0,1,4], [1,2,5], [2,3,6], [1,5,4], [2,6,5], [4,5,7], [5,6,8], + [5,8,7], [7,8,9]] +triang = mtri.Triangulation(x, y, triangles) + +# Interpolate to regularly-spaced quad grid. +z = np.cos(1.5*x)*np.cos(1.5*y) +interp = mtri.LinearTriInterpolator(triang, z) +xi, yi = np.meshgrid(np.linspace(0, 3, 20), np.linspace(0, 3, 20)) +zi = interp(xi, yi) + +# Plot the triangulation. +plt.subplot(121) +plt.tricontourf(triang, z) +plt.triplot(triang, 'ko-') +plt.title('Triangular grid') + +# Plot interpolation to quad grid. +plt.subplot(122) +plt.contourf(xi, yi, zi) +plt.plot(xi, yi, 'k-', alpha=0.5) +plt.plot(xi.T, yi.T, 'k-', alpha=0.5) +plt.title('Linear interpolation') + +plt.show() diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py index 0239bd379417..254038c3ccd7 100644 --- a/lib/matplotlib/tests/test_triangulation.py +++ b/lib/matplotlib/tests/test_triangulation.py @@ -219,3 +219,34 @@ def test_trifinder(): assert_equal(trifinder, triang.get_trifinder()) tris = trifinder(xs, ys) assert_array_equal(tris, [-1, -1, 1, -1]) + +def test_triinterp(): + # Test points within triangles of masked triangulation. + x,y = np.meshgrid(np.arange(4), np.arange(4)) + x = x.ravel() + y = y.ravel() + z = 1.23*x - 4.79*y + 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) + linear_interp = mtri.LinearTriInterpolator(triang, z) + + xs = np.linspace(0.25, 2.75, 6) + ys = [0.25, 0.75, 2.25, 2.75] + xs,ys = np.meshgrid(xs,ys) + xs = xs.ravel() + ys = ys.ravel() + zs = linear_interp(xs, ys) + assert_array_almost_equal(zs, (1.23*xs - 4.79*ys)) + + # Test points outside triangulation. + xs = [-0.25, 1.25, 1.75, 3.25] + ys = xs + xs, ys = np.meshgrid(xs,ys) + xs = xs.ravel() + ys = ys.ravel() + zs = linear_interp(xs, ys) + assert_array_equal(zs.mask, [True]*16) diff --git a/lib/matplotlib/tri/__init__.py b/lib/matplotlib/tri/__init__.py index 039ca110091e..5624c30932de 100644 --- a/lib/matplotlib/tri/__init__.py +++ b/lib/matplotlib/tri/__init__.py @@ -6,5 +6,6 @@ from triangulation import * from tricontour import * from trifinder import * +from triinterpolate import * from tripcolor import * from triplot import * diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp index f9deef5b8d2c..11767423c42b 100644 --- a/lib/matplotlib/tri/_tri.cpp +++ b/lib/matplotlib/tri/_tri.cpp @@ -137,6 +137,11 @@ double XYZ::dot(const XYZ& other) const return x*other.x + y*other.y + z*other.z; } +double XYZ::length_squared() const +{ + return x*x + y*y + z*z; +} + XYZ XYZ::operator-(const XYZ& other) const { return XYZ(x - other.x, y - other.y, z - other.z); @@ -375,6 +380,95 @@ void Triangulation::calculate_neighbors() // boundary edges, but the boundaries are calculated separately elsewhere. } +Py::Object Triangulation::calculate_plane_coefficients(const Py::Tuple &args) +{ + _VERBOSE("Triangulation::calculate_plane_coefficients"); + args.verify_length(1); + + PyArrayObject* z = (PyArrayObject*)PyArray_ContiguousFromObject( + args[0].ptr(), PyArray_DOUBLE, 1, 1); + if (z == 0 || PyArray_DIM(z,0) != PyArray_DIM(_x,0)) { + Py_XDECREF(z); + throw Py::ValueError( + "z array must have same length as triangulation x and y arrays"); + } + const double* zs = (const double*)PyArray_DATA(z); + + npy_intp dims[2] = {_ntri, 3}; + PyArrayObject* planes_array = (PyArrayObject*)PyArray_SimpleNew( + 2, dims, PyArray_DOUBLE); + double* planes = (double*)PyArray_DATA(planes_array); + const int* tris = get_triangles_ptr(); + const double* xs = (const double*)PyArray_DATA(_x); + const double* ys = (const double*)PyArray_DATA(_y); + for (int tri = 0; tri < _ntri; ++tri) + { + if (is_masked(tri)) + { + *planes++ = 0.0; + *planes++ = 0.0; + *planes++ = 0.0; + tris += 3; + } + 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 + XYZ point0(xs[*tris], ys[*tris], zs[*tris]); + tris++; + XYZ point1(xs[*tris], ys[*tris], zs[*tris]); + tris++; + XYZ point2(xs[*tris], ys[*tris], zs[*tris]); + tris++; + + XYZ normal = (point1 - point0).cross(point2 - point0); + + if (normal.z == 0.0) + { + // Normal is in x-y plane which means triangle consists of + // colinear points. Try to do the best we can by taking plane + // through longest side of triangle. + double length_sqr_01 = (point1 - point0).length_squared(); + double length_sqr_12 = (point2 - point1).length_squared(); + double length_sqr_20 = (point0 - point2).length_squared(); + if (length_sqr_01 > length_sqr_12) + { + if (length_sqr_01 > length_sqr_20) + normal = normal.cross(point1 - point0); + else + normal = normal.cross(point0 - point2); + } + else + { + if (length_sqr_12 > length_sqr_20) + normal = normal.cross(point2 - point1); + else + normal = normal.cross(point0 - point2); + } + + if (normal.z == 0.0) + { + // The 3 triangle points have identical x and y! The best + // we can do here is take normal = (0,0,1) and for the + // constant p take the mean of the 3 points' z-values. + normal = XYZ(0.0, 0.0, 1.0); + point0.z = (point0.z + point1.z + point2.z) / 3.0; + } + } + + *planes++ = -normal.x / normal.z; // x + *planes++ = -normal.y / normal.z; // y + *planes++ = normal.dot(point0) / normal.z; // constant + } + } + + return Py::asObject((PyObject*)planes_array); +} + void Triangulation::correct_triangles() { int* triangles_ptr = (int*)PyArray_DATA(_triangles); @@ -506,6 +600,9 @@ void Triangulation::init_type() behaviors().name("Triangulation"); behaviors().doc("Triangulation"); + add_varargs_method("calculate_plane_coefficients", + &Triangulation::calculate_plane_coefficients, + "calculate_plane_coefficients(z)"); add_noargs_method("get_edges", &Triangulation::get_edges, "get_edges()"); add_noargs_method("get_neighbors", &Triangulation::get_neighbors, diff --git a/lib/matplotlib/tri/_tri.h b/lib/matplotlib/tri/_tri.h index 1e6a95556998..aa4fa145d277 100644 --- a/lib/matplotlib/tri/_tri.h +++ b/lib/matplotlib/tri/_tri.h @@ -115,6 +115,7 @@ struct XYZ XYZ(const double& x_, const double& y_, const double& z_); XYZ cross(const XYZ& other) const; double dot(const XYZ& other) const; + double length_squared() const; XYZ operator-(const XYZ& other) const; friend std::ostream& operator<<(std::ostream& os, const XYZ& xyz); @@ -189,6 +190,13 @@ class Triangulation : public Py::PythonExtension virtual ~Triangulation(); + /* 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]. */ + Py::Object calculate_plane_coefficients(const Py::Tuple &args); + // Return the boundaries collection, creating it if necessary. const Boundaries& get_boundaries() const; diff --git a/lib/matplotlib/tri/triangulation.py b/lib/matplotlib/tri/triangulation.py index fb56d4d9bd76..d87cddd8a878 100644 --- a/lib/matplotlib/tri/triangulation.py +++ b/lib/matplotlib/tri/triangulation.py @@ -88,6 +88,16 @@ def __init__(self, x, y, triangles=None, mask=None): # Default TriFinder not created until needed. self._trifinder = None + def calculate_plane_coefficients(self, z): + """ + Calculate plane equation coefficients for all unmasked triangles from + the point (x,y) coordinates and specified z-array of shape (npoints). + Returned array has shape (npoints,3) and allows z-value at (x,y) + position in triangle tri to be calculated using + z = array[tri,0]*x + array[tri,1]*y + array[tri,2]. + """ + return self.get_cpp_triangulation().calculate_plane_coefficients(z) + @property def edges(self): if self._edges is None: diff --git a/lib/matplotlib/tri/triinterpolate.py b/lib/matplotlib/tri/triinterpolate.py new file mode 100644 index 000000000000..b14c91ab5894 --- /dev/null +++ b/lib/matplotlib/tri/triinterpolate.py @@ -0,0 +1,101 @@ +from __future__ import print_function +from matplotlib.tri import Triangulation +from matplotlib.tri.trifinder import TriFinder +import numpy as np + + +class TriInterpolator(object): + """ + Abstract base class for classes used to perform interpolation on + triangular grids. + + Derived classes implement __call__(x,y) where x,y are array_like point + coordinates of the same shape, and that returns a masked array of the same + shape containing the interpolated z-values. + """ + def __init__(self, triangulation, z, trifinder=None): + if not isinstance(triangulation, Triangulation): + raise ValueError('Expected a Triangulation object') + self._triangulation = triangulation + + self._z = np.asarray(z) + if self._z.shape != self._triangulation.x.shape: + raise ValueError('z array must have same length as triangulation x' + ' and y arrays') + + if trifinder is not None and not isinstance(trifinder, TriFinder): + raise ValueError('Expected a TriFinder object') + self._trifinder = trifinder or self._triangulation.get_trifinder() + + +class LinearTriInterpolator(TriInterpolator): + """ + A LinearTriInterpolator performs linear interpolation on a triangular grid. + + Each triangle is represented by a plane so that an interpolated value at + point (x,y) lies on the plane of the triangle containing (x,y). + Interpolated values are therefore continuous across the triangulation, but + their first derivatives are discontinuous at edges between triangles. + """ + def __init__(self, triangulation, z, trifinder=None): + """ + *triangulation*: the :class:`~matplotlib.tri.Triangulation` to + interpolate over. + + *z*: array_like of shape (npoints). + Array of values, defined at grid points, to interpolate between. + + *trifinder*: optional :class:`~matplotlib.tri.TriFinder` object. + If this is not specified, the Triangulation's default TriFinder will + be used by calling :func:`matplotlib.tri.Triangulation.get_trifinder`. + """ + TriInterpolator.__init__(self, triangulation, z, trifinder) + + # Store plane coefficients for fast interpolation calculations. + self._plane_coefficients = \ + self._triangulation.calculate_plane_coefficients(self._z) + + # Store vectorized interpolation function, so can pass in arbitrarily + # shape arrays of x, y and tri and the _single_interp function is + # called in turn with scalar x, y and tri. + self._multi_interp = np.vectorize(self._single_interp, + otypes=[np.float]) + + def __call__(self, x, y): + """ + Return a masked array containing linearly interpolated values at the + specified x,y points. + + *x*, *y* are array_like x and y coordinates of the same shape and any + number of dimensions. + + Returned masked array has the same shape as *x* and *y*; values + corresponding to (x,y) points outside of the triangulation are masked + out. + """ + # Check arguments. + x = np.asarray(x, dtype=np.float64) + y = np.asarray(y, dtype=np.float64) + if x.shape != y.shape: + raise ValueError("x and y must be equal-shaped arrays") + + # Indices of triangles containing x, y points, or -1 for no triangles. + tris = self._trifinder(x, y) + + # Perform interpolation. + z = self._multi_interp(x, y, tris) + + # Return masked array. + return np.ma.masked_invalid(z, copy=False) + + def _single_interp(self, x, y, tri): + """ + Return single interpolated value at specified (*x*, *y*) coordinates + within triangle with index *tri*. Returns np.nan if tri == -1. + """ + if tri == -1: + return np.nan + else: + return (self._plane_coefficients[tri,0] * x + + self._plane_coefficients[tri,1] * y + + self._plane_coefficients[tri,2]) From 0a9f09c6d55337c521c0679ae2b6cda41b704640 Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Wed, 19 Dec 2012 09:36:21 +0000 Subject: [PATCH 3/5] Improved PEP7 compliance. --- lib/matplotlib/tri/_tri.cpp | 375 +++++++++++++++++------------------- lib/matplotlib/tri/_tri.h | 4 - 2 files changed, 177 insertions(+), 202 deletions(-) diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp index 11767423c42b..8c50c2449305 100644 --- a/lib/matplotlib/tri/_tri.cpp +++ b/lib/matplotlib/tri/_tri.cpp @@ -617,15 +617,6 @@ 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"); @@ -1145,7 +1136,8 @@ TrapezoidMapTriFinder::~TrapezoidMapTriFinder() clear(); } -bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) +bool +TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) { std::vector trapezoids; if (!find_trapezoids_intersecting_edge(edge, trapezoids)) @@ -1162,8 +1154,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) // 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) - { + 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); @@ -1183,8 +1174,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) // 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) - { + if (start_trap && end_trap) { // Edge intersects a single trapezoid. if (have_left) left = new Trapezoid(old->left, p, old->below, old->above); @@ -1194,34 +1184,29 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) right = new Trapezoid(q, old->right, old->below, old->above); // Set pairs of trapezoid neighbours. - if (have_left) - { + 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 - { + else { below->set_lower_left(old->lower_left); above->set_upper_left(old->upper_left); } - if (have_right) - { + 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 - { + else { below->set_lower_right(old->lower_right); above->set_upper_right(old->upper_right); } } - else if (start_trap) - { + else if (start_trap) { // Old trapezoid is the first of 2+ trapezoids that the edge // intersects. if (have_left) @@ -1230,15 +1215,13 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) above = new Trapezoid(p, old->right, edge, old->above); // Set pairs of trapezoid neighbours. - if (have_left) - { + 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 - { + else { below->set_lower_left(old->lower_left); above->set_upper_left(old->upper_left); } @@ -1246,20 +1229,17 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) below->set_lower_right(old->lower_right); above->set_upper_right(old->upper_right); } - else if (end_trap) - { + else if (end_trap) { // Old trapezoid is the last of 2+ trapezoids that the edge // intersects. - if (left_below->below == old->below) - { + 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) - { + if (left_above->above == old->above) { above = left_above; above->right = q; } @@ -1270,22 +1250,19 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) right = new Trapezoid(q, old->right, old->below, old->above); // Set pairs of trapezoid neighbours. - if (have_right) - { + 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 - { + 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) - { + if (below != left_below) { below->set_upper_left(left_below); if (old->lower_left == left_old) below->set_lower_left(left_below); @@ -1293,8 +1270,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) below->set_lower_left(old->lower_left); } - if (above != left_above) - { + if (above != left_above) { above->set_lower_left(left_above); if (old->upper_left == left_old) above->set_upper_left(left_above); @@ -1302,20 +1278,17 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) above->set_upper_left(old->upper_left); } } - else // Middle trapezoid. - { + 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) - { + 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) - { + if (left_above->above == old->above) { above = left_above; above->right = old->right; } @@ -1323,8 +1296,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) above = new Trapezoid(old->left, old->right, edge, old->above); // Connect to new trapezoids replacing prevOld. - if (below != left_below) // below is new. - { + 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); @@ -1332,8 +1304,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) below->set_lower_left(old->lower_left); } - if (above != left_above) // above is new. - { + 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); @@ -1369,8 +1340,7 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) delete old_node; // Clearing up. - if (!end_trap) - { + if (!end_trap) { // Prepare for next loop. left_old = old; left_above = above; @@ -1381,7 +1351,8 @@ bool TrapezoidMapTriFinder::add_edge_to_tree(const Edge& edge) return true; } -void TrapezoidMapTriFinder::clear() +void +TrapezoidMapTriFinder::clear() { delete [] _points; _points = 0; @@ -1392,7 +1363,8 @@ void TrapezoidMapTriFinder::clear() _tree = 0; } -Py::Object TrapezoidMapTriFinder::find_many(const Py::Tuple& args) +Py::Object +TrapezoidMapTriFinder::find_many(const Py::Tuple& args) { args.verify_length(2); @@ -1406,8 +1378,7 @@ Py::Object TrapezoidMapTriFinder::find_many(const Py::Tuple& args) for (int i = 0; ok && i < ndim; ++i) ok = (PyArray_DIM(x,i) == PyArray_DIM(y,i)); - if (!ok) - { + if (!ok) { Py_XDECREF(x); Py_XDECREF(y); throw Py::ValueError("x and y must be array_like with same shape"); @@ -1430,14 +1401,16 @@ Py::Object TrapezoidMapTriFinder::find_many(const Py::Tuple& args) return Py::asObject((PyObject*)tri); } -int TrapezoidMapTriFinder::find_one(const XY& xy) +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( +bool +TrapezoidMapTriFinder::find_trapezoids_intersecting_edge( const Edge& edge, std::vector& trapezoids) { @@ -1445,24 +1418,20 @@ bool TrapezoidMapTriFinder::find_trapezoids_intersecting_edge( // checks to deal with simple colinear (i.e. invalid) triangles. trapezoids.clear(); Trapezoid* trapezoid = _tree->search(edge); - if (trapezoid == 0) - { + 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)) - { + while (edge.right->is_right_of(*trapezoid->right)) { int orient = edge.get_point_orientation(*trapezoid->right); - if (orient == 0) - { + if (orient == 0) { if (edge.point_below == trapezoid->right) orient = +1; else if (edge.point_above == trapezoid->right) orient = -1; - else - { + else { assert(0 && "Unable to deal with point on edge"); return false; } @@ -1473,8 +1442,7 @@ bool TrapezoidMapTriFinder::find_trapezoids_intersecting_edge( else if (orient == +1) trapezoid = trapezoid->upper_right; - if (trapezoid == 0) - { + if (trapezoid == 0) { assert(0 && "Expected trapezoid neighbor"); return false; } @@ -1484,7 +1452,8 @@ bool TrapezoidMapTriFinder::find_trapezoids_intersecting_edge( return true; } -Py::Object TrapezoidMapTriFinder::get_tree_stats() +Py::Object +TrapezoidMapTriFinder::get_tree_stats() { _VERBOSE("TrapezoidMapTriFinder::get_tree_stats"); @@ -1502,12 +1471,14 @@ Py::Object TrapezoidMapTriFinder::get_tree_stats() return list; } -const Triangulation& TrapezoidMapTriFinder::get_triangulation() const +const Triangulation& +TrapezoidMapTriFinder::get_triangulation() const { return *(Triangulation*)_triangulation.ptr(); } -void TrapezoidMapTriFinder::init_type() +void +TrapezoidMapTriFinder::init_type() { _VERBOSE("TrapezoidMapTriFinder::init_type"); @@ -1528,7 +1499,8 @@ void TrapezoidMapTriFinder::init_type() "print_tree()"); } -Py::Object TrapezoidMapTriFinder::initialize() +Py::Object +TrapezoidMapTriFinder::initialize() { _VERBOSE("TrapezoidMapTriFinder::initialize"); @@ -1540,8 +1512,7 @@ Py::Object TrapezoidMapTriFinder::initialize() int npoints = triang.get_npoints(); _points = new Point[npoints + 4]; BoundingBox bbox; - for (int i = 0; i < npoints; ++i) - { + 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) @@ -1555,13 +1526,11 @@ Py::Object TrapezoidMapTriFinder::initialize() // 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) - { + if (bbox.empty) { bbox.add(XY(0.0, 0.0)); bbox.add(XY(1.0, 1.0)); } - else - { + else { const double small = 0.1; // Any value > 0.0 bbox.expand( (bbox.upper - bbox.lower)*small ); } @@ -1572,25 +1541,23 @@ Py::Object TrapezoidMapTriFinder::initialize() // 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)); + _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) - { + 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); + 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)) - { + 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); @@ -1618,8 +1585,7 @@ Py::Object TrapezoidMapTriFinder::initialize() // Add edges, one at a time, to tree. unsigned int nedges = _edges.size(); - for (unsigned int index = 2; index < nedges; ++index) - { + 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); @@ -1628,7 +1594,8 @@ Py::Object TrapezoidMapTriFinder::initialize() return Py::None(); } -Py::Object TrapezoidMapTriFinder::print_tree() +Py::Object +TrapezoidMapTriFinder::print_tree() { _VERBOSE("TrapezoidMapTriFinder::print_tree"); @@ -1658,29 +1625,30 @@ TrapezoidMapTriFinder::Edge::Edge(const Point* left_, assert(triangle_above >= -1 && "Invalid triangle above index"); } -int TrapezoidMapTriFinder::Edge::get_point_orientation(const XY& xy) const +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 +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 +double +TrapezoidMapTriFinder::Edge::get_y_at_x(const double& x) const { - if (left->x == right->x) - { + 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 - { + 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); @@ -1689,18 +1657,21 @@ double TrapezoidMapTriFinder::Edge::get_y_at_x(const double& x) const } } -bool TrapezoidMapTriFinder::Edge::has_point(const Point* point) const +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 +bool +TrapezoidMapTriFinder::Edge::operator==(const Edge& other) const { return this == &other; } -void TrapezoidMapTriFinder::Edge::print_debug() const +void +TrapezoidMapTriFinder::Edge::print_debug() const { std::cout << "Edge " << *this << " tri_below=" << triangle_below << " tri_above=" << triangle_above << std::endl; @@ -1742,8 +1713,7 @@ TrapezoidMapTriFinder::Node::Node(Trapezoid* trapezoid) TrapezoidMapTriFinder::Node::~Node() { - switch (_type) - { + switch (_type) { case Type_XNode: if (_union.xnode.left->remove_parent(this)) delete _union.xnode.left; @@ -1762,7 +1732,8 @@ TrapezoidMapTriFinder::Node::~Node() } } -void TrapezoidMapTriFinder::Node::add_parent(Node* parent) +void +TrapezoidMapTriFinder::Node::add_parent(Node* parent) { assert(parent != 0 && "Null parent"); assert(parent != this && "Cannot be parent of self"); @@ -1770,21 +1741,20 @@ void TrapezoidMapTriFinder::Node::add_parent(Node* parent) _parents.push_back(parent); } -void TrapezoidMapTriFinder::Node::assert_valid(bool tree_complete) const +void +TrapezoidMapTriFinder::Node::assert_valid(bool tree_complete) const { #ifndef NDEBUG // Check parents. for (Parents::const_iterator it = _parents.begin(); - it != _parents.end(); ++it) - { + 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) - { + switch (_type) { case Type_XNode: assert(_union.xnode.left != 0 && "Null left child"); assert(_union.xnode.left->has_parent(this) && "Incorrect parent"); @@ -1811,8 +1781,9 @@ void TrapezoidMapTriFinder::Node::assert_valid(bool tree_complete) const #endif } -void TrapezoidMapTriFinder::Node::get_stats(int depth, - NodeStats& stats) const +void +TrapezoidMapTriFinder::Node::get_stats(int depth, + NodeStats& stats) const { stats.node_count++; if (depth > stats.max_depth) @@ -1822,8 +1793,7 @@ void TrapezoidMapTriFinder::Node::get_stats(int depth, stats.max_parent_count = std::max(stats.max_parent_count, static_cast(_parents.size())); - switch (_type) - { + switch (_type) { case Type_XNode: _union.xnode.left->get_stats(depth+1, stats); _union.xnode.right->get_stats(depth+1, stats); @@ -1840,10 +1810,10 @@ void TrapezoidMapTriFinder::Node::get_stats(int depth, } } -int TrapezoidMapTriFinder::Node::get_tri() const +int +TrapezoidMapTriFinder::Node::get_tri() const { - switch (_type) - { + switch (_type) { case Type_XNode: return _union.xnode.point->tri; case Type_YNode: @@ -1859,36 +1829,39 @@ int TrapezoidMapTriFinder::Node::get_tri() const } } -bool TrapezoidMapTriFinder::Node::has_child(const Node* child) const +bool +TrapezoidMapTriFinder::Node::has_child(const Node* child) const { assert(child != 0 && "Null child node"); - switch (_type) - { + 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); + return (_union.ynode.below == child || + _union.ynode.above == child); default: // Type_TrapezoidNode: return false; } } -bool TrapezoidMapTriFinder::Node::has_no_parents() const +bool +TrapezoidMapTriFinder::Node::has_no_parents() const { return _parents.empty(); } -bool TrapezoidMapTriFinder::Node::has_parent(const Node* parent) const +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 +void +TrapezoidMapTriFinder::Node::print(int depth /* = 0 */) const { for (int i = 0; i < depth; ++i) std::cout << " "; - switch (_type) - { + switch (_type) { case Type_XNode: std::cout << "XNode " << *_union.xnode.point << std::endl; _union.xnode.left->print(depth + 1); @@ -1909,7 +1882,8 @@ void TrapezoidMapTriFinder::Node::print(int depth /* = 0 */) const } } -bool TrapezoidMapTriFinder::Node::remove_parent(Node* parent) +bool +TrapezoidMapTriFinder::Node::remove_parent(Node* parent) { assert(parent != 0 && "Null parent"); assert(parent != this && "Cannot be parent of self"); @@ -1919,11 +1893,10 @@ bool TrapezoidMapTriFinder::Node::remove_parent(Node* parent) return _parents.empty(); } -void TrapezoidMapTriFinder::Node::replace_child(Node* old_child, - Node* new_child) +void +TrapezoidMapTriFinder::Node::replace_child(Node* old_child, Node* new_child) { - switch (_type) - { + switch (_type) { case Type_XNode: assert((_union.xnode.left == old_child || _union.xnode.right == old_child) && "Not a child Node"); @@ -1950,7 +1923,8 @@ void TrapezoidMapTriFinder::Node::replace_child(Node* old_child, new_child->add_parent(this); } -void TrapezoidMapTriFinder::Node::replace_with(Node* new_node) +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 @@ -1959,11 +1933,10 @@ void TrapezoidMapTriFinder::Node::replace_with(Node* new_node) _parents.front()->replace_child(this, new_node); } -const TrapezoidMapTriFinder::Node* TrapezoidMapTriFinder::Node::search( - const XY& xy) +const TrapezoidMapTriFinder::Node* +TrapezoidMapTriFinder::Node::search(const XY& xy) { - switch (_type) - { + switch (_type) { case Type_XNode: if (xy == *_union.xnode.point) return this; @@ -1971,8 +1944,7 @@ const TrapezoidMapTriFinder::Node* TrapezoidMapTriFinder::Node::search( return _union.xnode.right->search(xy); else return _union.xnode.left->search(xy); - case Type_YNode: - { + case Type_YNode: { int orient = _union.ynode.edge->get_point_orientation(xy); if (orient == 0) return this; @@ -1986,34 +1958,32 @@ const TrapezoidMapTriFinder::Node* TrapezoidMapTriFinder::Node::search( } } -TrapezoidMapTriFinder::Trapezoid* TrapezoidMapTriFinder::Node::search( - const Edge& edge) +TrapezoidMapTriFinder::Trapezoid* +TrapezoidMapTriFinder::Node::search(const Edge& edge) { - switch (_type) - { + switch (_type) { case Type_XNode: if (edge.left == _union.xnode.point) return _union.xnode.right->search(edge); - else - { + 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) - { + 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) + 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) + else if (_union.ynode.edge->triangle_below == + edge.triangle_above) return _union.ynode.below->search(edge); - else - { - assert(0 && "Invalid triangulation, common left points"); + else { + assert(0 && + "Invalid triangulation, common left points"); return 0; } } @@ -2022,18 +1992,18 @@ TrapezoidMapTriFinder::Trapezoid* TrapezoidMapTriFinder::Node::search( else return _union.ynode.below->search(edge); } - else if (edge.right == _union.ynode.edge->right) - { + 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) + 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) + else if (_union.ynode.edge->triangle_below == + edge.triangle_above) return _union.ynode.below->search(edge); - else - { - assert(0 && "Invalid triangulation, common right points"); + else { + assert(0 && + "Invalid triangulation, common right points"); return 0; } } @@ -2042,11 +2012,10 @@ TrapezoidMapTriFinder::Trapezoid* TrapezoidMapTriFinder::Node::search( else return _union.ynode.above->search(edge); } - else - { - int orient = _union.ynode.edge->get_point_orientation(*edge.left); - if (orient == 0) - { + 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)) @@ -2054,8 +2023,7 @@ TrapezoidMapTriFinder::Trapezoid* TrapezoidMapTriFinder::Node::search( else if (_union.ynode.edge->point_below != 0 && edge.has_point(_union.ynode.edge->point_below)) orient = +1; - else - { + else { assert(0 && "Invalid triangulation, point on edge"); return 0; } @@ -2083,39 +2051,40 @@ TrapezoidMapTriFinder::Trapezoid::Trapezoid(const Point* left_, assert(right->is_right_of(*left) && "Incorrect point order"); } -void TrapezoidMapTriFinder::Trapezoid::assert_valid(bool tree_complete) const +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 && + 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 && + 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 && + 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 && + 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"); @@ -2123,39 +2092,43 @@ void TrapezoidMapTriFinder::Trapezoid::assert_valid(bool tree_complete) const assert(trapezoid_node != 0 && "Null trapezoid_node"); - if (tree_complete) - { + 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 +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 +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 +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 +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 +void +TrapezoidMapTriFinder::Trapezoid::print_debug() const { std::cout << "Trapezoid " << this << " left=" << *left @@ -2173,28 +2146,32 @@ void TrapezoidMapTriFinder::Trapezoid::print_debug() const << " urp=" << get_upper_right_point() << std::endl; } -void TrapezoidMapTriFinder::Trapezoid::set_lower_left(Trapezoid* lower_left_) +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_) +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_) +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_) +void +TrapezoidMapTriFinder::Trapezoid::set_upper_right(Trapezoid* upper_right_) { upper_right = upper_right_; if (upper_right != 0) @@ -2206,7 +2183,8 @@ RandomNumberGenerator::RandomNumberGenerator(unsigned long seed) : _M(21870), _A(1291), _C(4621), _seed(seed % _M) {} -unsigned long RandomNumberGenerator::operator()(unsigned long max_value) +unsigned long +RandomNumberGenerator::operator()(unsigned long max_value) { _seed = (_seed*_A + _C) % _M; return (_seed*max_value) / _M; @@ -2354,7 +2332,8 @@ 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) +Py::Object +TriModule::new_TrapezoidMapTriFinder(const Py::Tuple &args) { _VERBOSE("TriModule::new_TrapezoidMapTriFinder"); args.verify_length(1); diff --git a/lib/matplotlib/tri/_tri.h b/lib/matplotlib/tri/_tri.h index aa4fa145d277..3662678fd8f1 100644 --- a/lib/matplotlib/tri/_tri.h +++ b/lib/matplotlib/tri/_tri.h @@ -239,10 +239,6 @@ class Triangulation : public Py::PythonExtension // 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). */ From 01ddc8aa220947149583bae11b454dc9466952f3 Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Tue, 8 Jan 2013 08:57:38 +0000 Subject: [PATCH 4/5] Added try..catch block to Triangulation::calculate_plane_coefficients. --- lib/matplotlib/tri/_tri.cpp | 141 ++++++++++++++++++++---------------- 1 file changed, 78 insertions(+), 63 deletions(-) diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp index 8c50c2449305..379e72870af0 100644 --- a/lib/matplotlib/tri/_tri.cpp +++ b/lib/matplotlib/tri/_tri.cpp @@ -392,80 +392,95 @@ Py::Object Triangulation::calculate_plane_coefficients(const Py::Tuple &args) throw Py::ValueError( "z array must have same length as triangulation x and y arrays"); } - const double* zs = (const double*)PyArray_DATA(z); - - npy_intp dims[2] = {_ntri, 3}; - PyArrayObject* planes_array = (PyArrayObject*)PyArray_SimpleNew( - 2, dims, PyArray_DOUBLE); - double* planes = (double*)PyArray_DATA(planes_array); - const int* tris = get_triangles_ptr(); - const double* xs = (const double*)PyArray_DATA(_x); - const double* ys = (const double*)PyArray_DATA(_y); - for (int tri = 0; tri < _ntri; ++tri) + + PyArrayObject* planes_array = 0; // Array to return. + + try { - if (is_masked(tri)) - { - *planes++ = 0.0; - *planes++ = 0.0; - *planes++ = 0.0; - tris += 3; - } - else + const double* zs = (const double*)PyArray_DATA(z); + + npy_intp dims[2] = {_ntri, 3}; + planes_array = (PyArrayObject*)PyArray_SimpleNew(2, dims, + PyArray_DOUBLE); + double* planes = (double*)PyArray_DATA(planes_array); + const int* tris = get_triangles_ptr(); + const double* xs = (const double*)PyArray_DATA(_x); + const double* ys = (const double*)PyArray_DATA(_y); + for (int tri = 0; tri < _ntri; ++tri) { - // 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 - XYZ point0(xs[*tris], ys[*tris], zs[*tris]); - tris++; - XYZ point1(xs[*tris], ys[*tris], zs[*tris]); - tris++; - XYZ point2(xs[*tris], ys[*tris], zs[*tris]); - tris++; - - XYZ normal = (point1 - point0).cross(point2 - point0); - - if (normal.z == 0.0) + if (is_masked(tri)) { - // Normal is in x-y plane which means triangle consists of - // colinear points. Try to do the best we can by taking plane - // through longest side of triangle. - double length_sqr_01 = (point1 - point0).length_squared(); - double length_sqr_12 = (point2 - point1).length_squared(); - double length_sqr_20 = (point0 - point2).length_squared(); - if (length_sqr_01 > length_sqr_12) - { - if (length_sqr_01 > length_sqr_20) - normal = normal.cross(point1 - point0); - else - normal = normal.cross(point0 - point2); - } - else - { - if (length_sqr_12 > length_sqr_20) - normal = normal.cross(point2 - point1); - else - normal = normal.cross(point0 - point2); - } + *planes++ = 0.0; + *planes++ = 0.0; + *planes++ = 0.0; + tris += 3; + } + 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 + XYZ point0(xs[*tris], ys[*tris], zs[*tris]); + tris++; + XYZ point1(xs[*tris], ys[*tris], zs[*tris]); + tris++; + XYZ point2(xs[*tris], ys[*tris], zs[*tris]); + tris++; + + XYZ normal = (point1 - point0).cross(point2 - point0); if (normal.z == 0.0) { - // The 3 triangle points have identical x and y! The best - // we can do here is take normal = (0,0,1) and for the - // constant p take the mean of the 3 points' z-values. - normal = XYZ(0.0, 0.0, 1.0); - point0.z = (point0.z + point1.z + point2.z) / 3.0; + // Normal is in x-y plane which means triangle consists of + // colinear points. Try to do the best we can by taking + // plane through longest side of triangle. + double length_sqr_01 = (point1 - point0).length_squared(); + double length_sqr_12 = (point2 - point1).length_squared(); + double length_sqr_20 = (point0 - point2).length_squared(); + if (length_sqr_01 > length_sqr_12) + { + if (length_sqr_01 > length_sqr_20) + normal = normal.cross(point1 - point0); + else + normal = normal.cross(point0 - point2); + } + else + { + if (length_sqr_12 > length_sqr_20) + normal = normal.cross(point2 - point1); + else + normal = normal.cross(point0 - point2); + } + + if (normal.z == 0.0) + { + // The 3 triangle points have identical x and y! The + // best we can do here is take normal = (0,0,1) and for + // the constant p take the mean of the 3 points' + // z-values. + normal = XYZ(0.0, 0.0, 1.0); + point0.z = (point0.z + point1.z + point2.z) / 3.0; + } } - } - *planes++ = -normal.x / normal.z; // x - *planes++ = -normal.y / normal.z; // y - *planes++ = normal.dot(point0) / normal.z; // constant + *planes++ = -normal.x / normal.z; // x + *planes++ = -normal.y / normal.z; // y + *planes++ = normal.dot(point0) / normal.z; // constant + } } } + catch (...) + { + Py_DECREF(z); + Py_XDECREF(planes_array); + throw; + } + Py_DECREF(z); return Py::asObject((PyObject*)planes_array); } From 88dc57ed2ed116622710e548d7db1da6d1d88f21 Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Tue, 8 Jan 2013 16:54:03 +0000 Subject: [PATCH 5/5] Improved PEP8 compliance. --- examples/pylab_examples/triinterp_demo.py | 5 +- lib/matplotlib/tests/test_triangulation.py | 63 ++++++++++++---------- 2 files changed, 36 insertions(+), 32 deletions(-) diff --git a/examples/pylab_examples/triinterp_demo.py b/examples/pylab_examples/triinterp_demo.py index 6b5d0efbdf52..7aefd844c23a 100644 --- a/examples/pylab_examples/triinterp_demo.py +++ b/examples/pylab_examples/triinterp_demo.py @@ -4,14 +4,13 @@ import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np -import math # Create triangulation. x = np.asarray([0, 1, 2, 3, 0.5, 1.5, 2.5, 1, 2, 1.5]) y = np.asarray([0, 0, 0, 0, 1, 1, 1, 2, 2, 3]) -triangles = [[0,1,4], [1,2,5], [2,3,6], [1,5,4], [2,6,5], [4,5,7], [5,6,8], - [5,8,7], [7,8,9]] +triangles = [[0, 1, 4], [1, 2, 5], [2, 3, 6], [1, 5, 4], [2, 6, 5], [4, 5, 7], + [5, 6, 8], [5, 8, 7], [7, 8, 9]] triang = mtri.Triangulation(x, y, triangles) # Interpolate to regularly-spaced quad grid. diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py index 254038c3ccd7..889b02bfd948 100644 --- a/lib/matplotlib/tests/test_triangulation.py +++ b/lib/matplotlib/tests/test_triangulation.py @@ -134,12 +134,13 @@ def test_no_modify(): def test_trifinder(): # Test points within triangles of masked triangulation. - x,y = np.meshgrid(np.arange(4), np.arange(4)) + 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]] + 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) @@ -147,25 +148,27 @@ def test_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, 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]) + 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]) + 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] + 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]) + 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]) + 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. @@ -173,16 +176,17 @@ def test_trifinder(): # 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]] + 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) + 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]]) + 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. @@ -190,23 +194,23 @@ def test_trifinder(): # 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]] + 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] + 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) + 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]]) + 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]] + triangles = [[0, 1, 2], [1, 3, 2]] triang = mtri.Triangulation(x, y, triangles) trifinder = triang.get_trifinder() @@ -222,13 +226,14 @@ def test_trifinder(): def test_triinterp(): # Test points within triangles of masked triangulation. - x,y = np.meshgrid(np.arange(4), np.arange(4)) + x, y = np.meshgrid(np.arange(4), np.arange(4)) x = x.ravel() y = y.ravel() z = 1.23*x - 4.79*y - 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]] + 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) @@ -236,7 +241,7 @@ def test_triinterp(): xs = np.linspace(0.25, 2.75, 6) ys = [0.25, 0.75, 2.25, 2.75] - xs,ys = np.meshgrid(xs,ys) + xs, ys = np.meshgrid(xs, ys) xs = xs.ravel() ys = ys.ravel() zs = linear_interp(xs, ys) @@ -245,7 +250,7 @@ def test_triinterp(): # Test points outside triangulation. xs = [-0.25, 1.25, 1.75, 3.25] ys = xs - xs, ys = np.meshgrid(xs,ys) + xs, ys = np.meshgrid(xs, ys) xs = xs.ravel() ys = ys.ravel() zs = linear_interp(xs, ys)