From d271736ef2afa3329562f177bbea08677b4bde40 Mon Sep 17 00:00:00 2001 From: Ian Thomas Date: Tue, 28 Oct 2014 08:06:56 +0000 Subject: [PATCH] Final decxx corrections to PR 3723 --- lib/matplotlib/tri/_tri.cpp | 278 ++++++++++++---------------- lib/matplotlib/tri/_tri.h | 121 +++++------- lib/matplotlib/tri/_tri_wrapper.cpp | 257 +++++++++++-------------- lib/matplotlib/tri/_tri_wrapper.h | 24 --- lib/matplotlib/tri/trifinder.py | 12 +- src/numpy_cpp.h | 39 +++- 6 files changed, 322 insertions(+), 409 deletions(-) delete mode 100644 lib/matplotlib/tri/_tri_wrapper.h diff --git a/lib/matplotlib/tri/_tri.cpp b/lib/matplotlib/tri/_tri.cpp index a4bfb77bcb7f..74aab7e079e7 100644 --- a/lib/matplotlib/tri/_tri.cpp +++ b/lib/matplotlib/tri/_tri.cpp @@ -8,7 +8,6 @@ #define NO_IMPORT_ARRAY #include "_tri.h" -#include "_tri_wrapper.h" #include #include @@ -218,15 +217,13 @@ void write_contour(const Contour& contour) -Triangulation::Triangulation(PyArrayObject* x, - PyArrayObject* y, - PyArrayObject* triangles, - PyArrayObject* mask, - PyArrayObject* edges, - PyArrayObject* neighbors) - : _npoints(PyArray_DIM(x,0)), - _ntri(PyArray_DIM(triangles,0)), - _x(x), +Triangulation::Triangulation(const CoordinateArray& x, + const CoordinateArray& y, + const TriangleArray& triangles, + const MaskArray& mask, + const EdgeArray& edges, + const NeighborArray& neighbors) + : _x(x), _y(y), _triangles(triangles), _mask(mask), @@ -236,16 +233,6 @@ Triangulation::Triangulation(PyArrayObject* x, correct_triangles(); } -Triangulation::~Triangulation() -{ - Py_XDECREF(_x); - Py_XDECREF(_y); - Py_XDECREF(_triangles); - Py_XDECREF(_mask); - Py_XDECREF(_edges); - Py_XDECREF(_neighbors); -} - void Triangulation::calculate_boundaries() { get_neighbors(); // Ensure _neighbors has been created. @@ -254,7 +241,7 @@ void Triangulation::calculate_boundaries() // have a neighbor triangle. typedef std::set BoundaryEdges; BoundaryEdges boundary_edges; - for (int tri = 0; tri < _ntri; ++tri) { + for (int tri = 0; tri < get_ntri(); ++tri) { if (!is_masked(tri)) { for (int edge = 0; edge < 3; ++edge) { if (get_neighbor(tri, edge) == -1) { @@ -304,13 +291,13 @@ void Triangulation::calculate_boundaries() void Triangulation::calculate_edges() { - Py_XDECREF(_edges); + assert(_edges.empty() && "Expected empty edges array"); // Create set of all edges, storing them with start point index less than // end point index. typedef std::set EdgeSet; EdgeSet edge_set; - for (int tri = 0; tri < _ntri; ++tri) { + for (int tri = 0; tri < get_ntri(); ++tri) { if (!is_masked(tri)) { for (int edge = 0; edge < 3; edge++) { int start = get_triangle_point(tri, edge); @@ -322,23 +309,28 @@ void Triangulation::calculate_edges() // Convert to python _edges array. npy_intp dims[2] = {static_cast(edge_set.size()), 2}; - _edges = (PyArrayObject*)PyArray_SimpleNew(2, dims, NPY_INT); - int* edges_ptr = (int*)PyArray_DATA(_edges); + _edges = EdgeArray(dims); + + int i = 0; for (EdgeSet::const_iterator it = edge_set.begin(); it != edge_set.end(); ++it) { - *edges_ptr++ = it->start; - *edges_ptr++ = it->end; + _edges(i, 0) = it->start; + _edges(i++, 1) = it->end; } } void Triangulation::calculate_neighbors() { - Py_XDECREF(_neighbors); + assert(_neighbors.empty() && "Expected empty neighbors array"); // Create _neighbors array with shape (ntri,3) and initialise all to -1. - npy_intp dims[2] = {_ntri,3}; - _neighbors = (PyArrayObject*)PyArray_SimpleNew(2, dims, NPY_INT); - int* neighbors_ptr = (int*)PyArray_DATA(_neighbors); - std::fill(neighbors_ptr, neighbors_ptr + 3*_ntri, -1); + npy_intp dims[2] = {get_ntri(), 3}; + _neighbors = NeighborArray(dims); + + int tri, edge; + for (tri = 0; tri < get_ntri(); ++tri) { + for (edge = 0; edge < 3; ++edge) + _neighbors(tri, edge) = -1; + } // For each triangle edge (start to end point), find corresponding neighbor // edge from end to start point. Do this by traversing all edges and @@ -347,9 +339,9 @@ void Triangulation::calculate_neighbors() // already found. typedef std::map EdgeToTriEdgeMap; EdgeToTriEdgeMap edge_to_tri_edge_map; - for (int tri = 0; tri < _ntri; ++tri) { + for (tri = 0; tri < get_ntri(); ++tri) { if (!is_masked(tri)) { - for (int edge = 0; edge < 3; ++edge) { + for (edge = 0; edge < 3; ++edge) { int start = get_triangle_point(tri, edge); int end = get_triangle_point(tri, (edge+1)%3); EdgeToTriEdgeMap::iterator it = @@ -361,8 +353,8 @@ void Triangulation::calculate_neighbors() } else { // Neighbor edge found, set the two elements of _neighbors // and remove edge from edge_to_tri_edge_map. - neighbors_ptr[3*tri + edge] = it->second.tri; - neighbors_ptr[3*it->second.tri + it->second.edge] = tri; + _neighbors(tri, edge)= it->second.tri; + _neighbors(it->second.tri, it->second.edge) = tri; edge_to_tri_edge_map.erase(it); } } @@ -373,20 +365,18 @@ void Triangulation::calculate_neighbors() // boundary edges, but the boundaries are calculated separately elsewhere. } -void Triangulation::calculate_plane_coefficients( - PyArrayObject* z, PyArrayObject* plane_coefficients) +Triangulation::TwoCoordinateArray Triangulation::calculate_plane_coefficients( + const CoordinateArray& z) { - const double* zs = (const double*)PyArray_DATA(z); - double* planes = (double*)PyArray_DATA(plane_coefficients); - 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) { + npy_intp dims[2] = {get_ntri(), 3}; + Triangulation::TwoCoordinateArray planes(dims); + + int point; + for (int tri = 0; tri < get_ntri(); ++tri) { if (is_masked(tri)) { - *planes++ = 0.0; - *planes++ = 0.0; - *planes++ = 0.0; - tris += 3; + planes(tri, 0) = 0.0; + planes(tri, 1) = 0.0; + planes(tri, 2) = 0.0; } else { // Equation of plane for all points r on plane is r.normal = p @@ -396,12 +386,12 @@ void Triangulation::calculate_plane_coefficients( // 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 side01 = XYZ(xs[*tris], ys[*tris], zs[*tris]) - point0; - tris++; - XYZ side02 = XYZ(xs[*tris], ys[*tris], zs[*tris]) - point0; - tris++; + point = _triangles(tri, 0); + XYZ point0(_x(point), _y(point), z(point)); + point = _triangles(tri, 1); + XYZ side01 = XYZ(_x(point), _y(point), z(point)) - point0; + point = _triangles(tri, 2); + XYZ side02 = XYZ(_x(point), _y(point), z(point)) - point0; XYZ normal = side01.cross(side02); @@ -413,32 +403,32 @@ void Triangulation::calculate_plane_coefficients( side02.x*side02.x + side02.y*side02.y); double a = (side01.x*side01.z + side02.x*side02.z) / sum2; double b = (side01.y*side01.z + side02.y*side02.z) / sum2; - *planes++ = a; - *planes++ = b; - *planes++ = point0.z - a*point0.x - b*point0.y; + planes(tri, 0) = a; + planes(tri, 1) = b; + planes(tri, 2) = point0.z - a*point0.x - b*point0.y; } else { - *planes++ = -normal.x / normal.z; // x - *planes++ = -normal.y / normal.z; // y - *planes++ = normal.dot(point0) / normal.z; // constant + planes(tri, 0) = -normal.x / normal.z; // x + planes(tri, 1) = -normal.y / normal.z; // y + planes(tri, 2) = normal.dot(point0) / normal.z; // constant } } } + + return planes; } void Triangulation::correct_triangles() { - int* triangles_ptr = (int*)PyArray_DATA(_triangles); - int* neighbors_ptr = _neighbors != 0 ? (int*)PyArray_DATA(_neighbors) : 0; - for (int tri = 0; tri < _ntri; ++tri) { - XY point0 = get_point_coords(*triangles_ptr++); - XY point1 = get_point_coords(*triangles_ptr++); - XY point2 = get_point_coords(*triangles_ptr++); + for (int tri = 0; tri < get_ntri(); ++tri) { + XY point0 = get_point_coords(_triangles(tri, 0)); + XY point1 = get_point_coords(_triangles(tri, 1)); + XY point2 = get_point_coords(_triangles(tri, 2)); if ( (point1 - point0).cross_z(point2 - point0) < 0.0) { // Triangle points are clockwise, so change them to anticlockwise. - std::swap(*(triangles_ptr-2), *(triangles_ptr-1)); - if (neighbors_ptr) - std::swap(*(neighbors_ptr+3*tri+1), *(neighbors_ptr+3*tri+2)); + std::swap(_triangles(tri, 1), _triangles(tri, 2)); + if (!_neighbors.empty()) + std::swap(_neighbors(tri, 1), _neighbors(tri, 2)); } } } @@ -465,27 +455,29 @@ void Triangulation::get_boundary_edge(const TriEdge& triEdge, int Triangulation::get_edge_in_triangle(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."); - const int* triangles_ptr = get_triangles_ptr() + 3*tri; + assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds"); + assert(point >= 0 && point < get_npoints() && "Point index out of bounds."); for (int edge = 0; edge < 3; ++edge) { - if (*triangles_ptr++ == point) return edge; + if (_triangles(tri, edge) == point) + return edge; } return -1; // point is not in triangle. } -PyArrayObject* Triangulation::get_edges() +Triangulation::EdgeArray& Triangulation::get_edges() { - if (_edges == 0) + if (_edges.empty()) calculate_edges(); return _edges; } int Triangulation::get_neighbor(int tri, int edge) const { - assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds"); + assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds"); assert(edge >= 0 && edge < 3 && "Edge index out of bounds"); - return get_neighbors_ptr()[3*tri + edge]; + if (_neighbors.empty()) + const_cast(*this).calculate_neighbors(); + return _neighbors(tri, edge); } TriEdge Triangulation::get_neighbor_edge(int tri, int edge) const @@ -500,42 +492,34 @@ TriEdge Triangulation::get_neighbor_edge(int tri, int edge) const (edge+1)%3))); } -PyArrayObject* Triangulation::get_neighbors() +Triangulation::NeighborArray& Triangulation::get_neighbors() { - if (_neighbors == 0) + if (_neighbors.empty()) calculate_neighbors(); return _neighbors; } -const int* Triangulation::get_neighbors_ptr() const -{ - if (_neighbors == 0) - const_cast(this)->calculate_neighbors(); - return (const int*)PyArray_DATA(_neighbors); -} - int Triangulation::get_npoints() const { - return _npoints; + return _x.size(); } int Triangulation::get_ntri() const { - return _ntri; + return _triangles.size(); } XY Triangulation::get_point_coords(int point) const { - assert(point >= 0 && point < _npoints && "Point index out of bounds."); - return XY( ((const double*)PyArray_DATA(_x))[point], - ((const double*)PyArray_DATA(_y))[point] ); + assert(point >= 0 && point < get_npoints() && "Point index out of bounds."); + return XY(_x(point), _y(point)); } int Triangulation::get_triangle_point(int tri, int edge) const { - assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds"); + assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds"); assert(edge >= 0 && edge < 3 && "Edge index out of bounds"); - return get_triangles_ptr()[3*tri + edge]; + return _triangles(tri, edge); } int Triangulation::get_triangle_point(const TriEdge& tri_edge) const @@ -543,29 +527,19 @@ int Triangulation::get_triangle_point(const TriEdge& tri_edge) const return get_triangle_point(tri_edge.tri, tri_edge.edge); } -const int* Triangulation::get_triangles_ptr() const -{ - return (const int*)PyArray_DATA(_triangles); -} - bool Triangulation::is_masked(int tri) const { - assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds."); - return _mask && *((bool*)PyArray_DATA(_mask) + tri); + assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds."); + return !_mask.empty() && _mask(tri); } -void Triangulation::set_mask(PyArrayObject* mask) +void Triangulation::set_mask(const MaskArray& mask) { - // No longer need old mask, whether it was set or not. - Py_XDECREF(_mask); - _mask = mask; // Clear derived fields so they are recalculated when needed. - Py_XDECREF(_edges); - _edges = 0; - Py_XDECREF(_neighbors); - _neighbors = 0; + _edges = EdgeArray(); + _neighbors = NeighborArray(); _boundaries.clear(); } @@ -585,21 +559,15 @@ void Triangulation::write_boundaries() const -TriContourGenerator::TriContourGenerator(PyObject* triangulation, - PyArrayObject* z) +TriContourGenerator::TriContourGenerator(Triangulation& triangulation, + const CoordinateArray& z) : _triangulation(triangulation), _z(z), - _interior_visited(2*get_triangulation().get_ntri()), + _interior_visited(2*_triangulation.get_ntri()), _boundaries_visited(0), _boundaries_used(0) {} -TriContourGenerator::~TriContourGenerator() -{ - Py_XDECREF(_triangulation); - Py_XDECREF(_z); -} - void TriContourGenerator::clear_visited_flags(bool include_boundaries) { // Clear _interiorVisited. @@ -643,6 +611,7 @@ PyObject* TriContourGenerator::contour_to_segs(const Contour& contour) *p++ = it->y; } if (PyList_SetItem(segs, i, (PyObject*)py_line)) { + Py_XDECREF(segs); PyErr_SetString(PyExc_RuntimeError, "Unable to set contour segments"); return NULL; @@ -684,6 +653,7 @@ PyObject* TriContourGenerator::contour_to_segs_and_kinds(const Contour& contour) PyObject* result = PyTuple_New(2); if (PyTuple_SetItem(result, 0, (PyObject*)segs) || PyTuple_SetItem(result, 1, (PyObject*)kinds)) { + Py_XDECREF(result); PyErr_SetString(PyExc_RuntimeError, "Unable to set contour segments and kinds"); return NULL; @@ -717,8 +687,8 @@ PyObject* TriContourGenerator::create_filled_contour(const double& lower_level, XY TriContourGenerator::edge_interp(int tri, int edge, const double& level) { - return interp(get_triangulation().get_triangle_point(tri, edge), - get_triangulation().get_triangle_point(tri, (edge+1)%3), + return interp(_triangulation.get_triangle_point(tri, edge), + _triangulation.get_triangle_point(tri, (edge+1)%3), level); } @@ -728,7 +698,7 @@ void TriContourGenerator::find_boundary_lines(Contour& contour, // Traverse boundaries to find starting points for all contour lines that // intersect the boundaries. For each starting point found, follow the // line to its end before continuing. - const Triangulation& triang = get_triangulation(); + const Triangulation& triang = _triangulation; const Boundaries& boundaries = get_boundaries(); for (Boundaries::const_iterator it = boundaries.begin(); it != boundaries.end(); ++it) { @@ -761,7 +731,7 @@ void TriContourGenerator::find_boundary_lines_filled(Contour& contour, // Traverse boundaries to find starting points for all contour lines that // intersect the boundaries. For each starting point found, follow the // line to its end before continuing. - const Triangulation& triang = get_triangulation(); + const Triangulation& triang = _triangulation; const Boundaries& boundaries = get_boundaries(); for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) { const Boundary& boundary = boundaries[i]; @@ -826,7 +796,7 @@ void TriContourGenerator::find_interior_lines(Contour& contour, bool on_upper, bool filled) { - const Triangulation& triang = get_triangulation(); + const Triangulation& triang = _triangulation; int ntri = triang.get_ntri(); for (int tri = 0; tri < ntri; ++tri) { int visited_index = (on_upper ? tri+ntri : tri); @@ -864,7 +834,7 @@ bool TriContourGenerator::follow_boundary(ContourLine& contour_line, const double& upper_level, bool on_upper) { - const Triangulation& triang = get_triangulation(); + const Triangulation& triang = _triangulation; const Boundaries& boundaries = get_boundaries(); // Have TriEdge to start at, need equivalent boundary edge. @@ -937,7 +907,7 @@ void TriContourGenerator::follow_interior(ContourLine& contour_line, while (true) { int visited_index = tri; if (on_upper) - visited_index += get_triangulation().get_ntri(); + visited_index += _triangulation.get_ntri(); // Check for end not on boundary. if (!end_on_boundary && _interior_visited[visited_index]) @@ -954,7 +924,7 @@ void TriContourGenerator::follow_interior(ContourLine& contour_line, contour_line.push_back(edge_interp(tri, edge, level)); // Move to next triangle. - TriEdge next_tri_edge = get_triangulation().get_neighbor_edge(tri,edge); + TriEdge next_tri_edge = _triangulation.get_neighbor_edge(tri,edge); // Check if ending on a boundary. if (end_on_boundary && next_tri_edge.tri == -1) @@ -967,20 +937,20 @@ void TriContourGenerator::follow_interior(ContourLine& contour_line, const TriContourGenerator::Boundaries& TriContourGenerator::get_boundaries() const { - return get_triangulation().get_boundaries(); + return _triangulation.get_boundaries(); } int TriContourGenerator::get_exit_edge(int tri, const double& level, bool on_upper) const { - assert(tri >= 0 && tri < get_triangulation().get_ntri() && + assert(tri >= 0 && tri < _triangulation.get_ntri() && "Triangle index out of bounds."); unsigned int config = - (get_z(get_triangulation().get_triangle_point(tri, 0)) >= level) | - (get_z(get_triangulation().get_triangle_point(tri, 1)) >= level) << 1 | - (get_z(get_triangulation().get_triangle_point(tri, 2)) >= level) << 2; + (get_z(_triangulation.get_triangle_point(tri, 0)) >= level) | + (get_z(_triangulation.get_triangle_point(tri, 1)) >= level) << 1 | + (get_z(_triangulation.get_triangle_point(tri, 2)) >= level) << 2; if (on_upper) config = 7-config; @@ -997,35 +967,30 @@ int TriContourGenerator::get_exit_edge(int tri, } } -const Triangulation& TriContourGenerator::get_triangulation() const -{ - return *((PyTriangulation*)_triangulation)->ptr; -} - const double& TriContourGenerator::get_z(int point) const { - assert(point >= 0 && point < get_triangulation().get_npoints() && + assert(point >= 0 && point < _triangulation.get_npoints() && "Point index out of bounds."); - return ((const double*)PyArray_DATA(_z))[point]; + return _z(point); } XY TriContourGenerator::interp(int point1, int point2, const double& level) const { - assert(point1 >= 0 && point1 < get_triangulation().get_npoints() && + assert(point1 >= 0 && point1 < _triangulation.get_npoints() && "Point index 1 out of bounds."); - assert(point2 >= 0 && point2 < get_triangulation().get_npoints() && + assert(point2 >= 0 && point2 < _triangulation.get_npoints() && "Point index 2 out of bounds."); assert(point1 != point2 && "Identical points"); double fraction = (get_z(point2) - level) / (get_z(point2) - get_z(point1)); - return get_triangulation().get_point_coords(point1)*fraction + - get_triangulation().get_point_coords(point2)*(1.0 - fraction); + return _triangulation.get_point_coords(point1)*fraction + + _triangulation.get_point_coords(point2)*(1.0 - fraction); } -TrapezoidMapTriFinder::TrapezoidMapTriFinder(PyObject* triangulation) +TrapezoidMapTriFinder::TrapezoidMapTriFinder(Triangulation& triangulation) : _triangulation(triangulation), _points(0), _tree(0) @@ -1034,7 +999,6 @@ TrapezoidMapTriFinder::TrapezoidMapTriFinder(PyObject* triangulation) TrapezoidMapTriFinder::~TrapezoidMapTriFinder() { clear(); - Py_XDECREF(_triangulation); } bool @@ -1264,22 +1228,20 @@ TrapezoidMapTriFinder::clear() _tree = 0; } -PyArrayObject* -TrapezoidMapTriFinder::find_many(PyArrayObject* x, PyArrayObject* y) +TrapezoidMapTriFinder::TriIndexArray +TrapezoidMapTriFinder::find_many(const CoordinateArray& x, + const CoordinateArray& y) { // Create integer array to return. - PyArrayObject* tri = (PyArrayObject*)PyArray_SimpleNew( - PyArray_NDIM(x), PyArray_DIMS(x), NPY_INT); + npy_intp n = x.dim(0); + npy_intp dims[1] = {n}; + TriIndexArray tri_indices(dims); // 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++)); + for (npy_intp i = 0; i < n; ++i) + tri_indices(i) = find_one(XY(x(i), y(i))); - return tri; + return tri_indices; } int @@ -1349,17 +1311,11 @@ TrapezoidMapTriFinder::get_tree_stats() stats.sum_trapezoid_depth / stats.trapezoid_count); } -const Triangulation& -TrapezoidMapTriFinder::get_triangulation() const -{ - return *((PyTriangulation*)_triangulation)->ptr; -} - void TrapezoidMapTriFinder::initialize() { clear(); - const Triangulation& triang = get_triangulation(); + const Triangulation& triang = _triangulation; // Set up points array, which contains all of the points in the // triangulation plus the 4 corners of the enclosing rectangle. diff --git a/lib/matplotlib/tri/_tri.h b/lib/matplotlib/tri/_tri.h index a40c0fa673fc..7485831bbce6 100644 --- a/lib/matplotlib/tri/_tri.h +++ b/lib/matplotlib/tri/_tri.h @@ -63,8 +63,7 @@ #ifndef _TRI_H #define _TRI_H -#include "Python.h" -#include "numpy/arrayobject.h" +#include "src/numpy_cpp.h" #include #include @@ -162,6 +161,13 @@ void write_contour(const Contour& contour); class Triangulation { public: + typedef numpy::array_view CoordinateArray; + typedef numpy::array_view TwoCoordinateArray; + typedef numpy::array_view TriangleArray; + typedef numpy::array_view MaskArray; + typedef numpy::array_view EdgeArray; + typedef numpy::array_view NeighborArray; + /* A single boundary is a vector of the TriEdges that make up that boundary * following it around with unmasked triangles on the left. */ typedef std::vector Boundary; @@ -180,26 +186,20 @@ class Triangulation * once. * neighbors: Optional int array of shape (ntri,3) indicating which * triangles are the neighbors of which TriEdges, or -1 if - * there is no such neighbor. - * Argument reference counts are not incremented in the constructor, but - * are decremented when no longer needed. */ - Triangulation(PyArrayObject* x, - PyArrayObject* y, - PyArrayObject* triangles, - PyArrayObject* mask, - PyArrayObject* edges, - PyArrayObject* neighbors); - - ~Triangulation(); + * there is no such neighbor. */ + Triangulation(const CoordinateArray& x, + const CoordinateArray& y, + const TriangleArray& triangles, + const MaskArray& mask, + const EdgeArray& edges, + const NeighborArray& neighbors); /* 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 must have shape (npoints,3) and allows + * 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]. - * Does not alter reference counts of z or plane_coefficients. */ - void calculate_plane_coefficients(PyArrayObject* z, - PyArrayObject* plane_coefficients); + * z = array[tri,0]*x + array[tri,1]*y + array[tri,2]. */ + TwoCoordinateArray calculate_plane_coefficients(const CoordinateArray& z); // Return the boundaries collection, creating it if necessary. const Boundaries& get_boundaries() const; @@ -209,9 +209,8 @@ class Triangulation int& boundary, int& edge) const; - /* Return a borrowed reference to the edges array, creating it if - * necessary. */ - PyArrayObject* get_edges(); + /* Return the edges array, creating it if necessary. */ + EdgeArray& get_edges(); /* Return the triangle index of the neighbor of the specified triangle * edge. */ @@ -221,9 +220,8 @@ class Triangulation * or TriEdge(-1,-1) if there is no such neighbor. */ TriEdge get_neighbor_edge(int tri, int edge) const; - /* Return a borrowed reference to the neighbors array, creating it if - * necessary. */ - PyArrayObject* get_neighbors(); + /* Return the neighbors array, creating it if necessary. */ + NeighborArray& get_neighbors(); // Return the number of points in this triangulation. int get_npoints() const; @@ -244,9 +242,9 @@ class Triangulation /* Set or clear the mask array. Clears various derived fields so they are * recalculated when next needed. - * mask: New reference to bool array of shape (ntri) indicating which - * triangles are masked, or 0 to clear mask. */ - void set_mask(PyArrayObject* mask); + * mask: bool array of shape (ntri) indicating which triangles are + * masked, or an empty array to clear mask. */ + void set_mask(const MaskArray& mask); // Debug function to write boundaries. void write_boundaries() const; @@ -294,29 +292,21 @@ class Triangulation * the specified triangle, or -1 if the point is not in the triangle. */ int get_edge_in_triangle(int tri, int point) const; - // Return pointer to contents of neighbors array. - const int* get_neighbors_ptr() const; - - // Return pointer to contents of triangles array. - const int* get_triangles_ptr() const; - int _npoints, _ntri; - // Variables shared with python, always set. - PyArrayObject* _x; // double array (npoints). - PyArrayObject* _y; // double array (npoints). - PyArrayObject* _triangles; // int array (ntri,3) of triangle point indices, + CoordinateArray _x, _y; // double array (npoints). + TriangleArray _triangles; // int array (ntri,3) of triangle point indices, // ordered anticlockwise. // Variables shared with python, may be zero. - PyArrayObject* _mask; // bool array (ntri). + MaskArray _mask; // bool array (ntri). // Derived variables shared with python, may be zero. If zero, are // recalculated when needed. - PyArrayObject* _edges; // int array (?,2) of start & end point indices. - PyArrayObject* _neighbors; // int array (ntri,3), neighbor triangle indices + EdgeArray _edges; // int array (?,2) of start & end point indices. + NeighborArray _neighbors; // int array (ntri,3), neighbor triangle indices // or -1 if no neighbor. // Variables internal to C++ only. @@ -334,20 +324,18 @@ class Triangulation class TriContourGenerator { public: + typedef Triangulation::CoordinateArray CoordinateArray; + /* Constructor. * triangulation: Triangulation to generate contours for. * z: Double array of shape (npoints) of z-values at triangulation - * points. - * Argument reference counts are not incremented in the constructor, but - * are decremented when no longer needed. */ - TriContourGenerator(PyObject* triangulation, - PyArrayObject* z); - - ~TriContourGenerator(); + * points. */ + TriContourGenerator(Triangulation& triangulation, + const CoordinateArray& z); /* Create and return a non-filled contour. * level: Contour level. - * Returns reference to new python list [segs0, segs1, ...] where + * Returns new python list [segs0, segs1, ...] where * segs0: double array of shape (?,2) of point coordinates of first * contour line, etc. */ PyObject* create_contour(const double& level); @@ -355,7 +343,7 @@ class TriContourGenerator /* Create and return a filled contour. * lower_level: Lower contour level. * upper_level: Upper contour level. - * Returns reference to new python tuple (segs, kinds) where + * Returns new python tuple (segs, kinds) where * segs: double array of shape (n_points,2) of all point coordinates, * kinds: ubyte array of shape (n_points) of all point code types. */ PyObject* create_filled_contour(const double& lower_level, @@ -371,13 +359,13 @@ class TriContourGenerator void clear_visited_flags(bool include_boundaries); /* Convert a non-filled Contour from C++ to Python. - * Returns reference to new python list [segs0, segs1, ...] where + * Returns new python list [segs0, segs1, ...] where * segs0: double array of shape (?,2) of point coordinates of first * contour line, etc. */ PyObject* contour_to_segs(const Contour& contour); /* Convert a filled Contour from C++ to Python. - * Returns reference to new python tuple (segs, kinds) where + * Returns new python tuple (segs, kinds) where * segs: double array of shape (n_points,2) of all point coordinates, * kinds: ubyte array of shape (n_points) of all point code types. */ PyObject* contour_to_segs_and_kinds(const Contour& contour); @@ -455,9 +443,6 @@ class TriContourGenerator * on_upper: Whether following upper or lower contour level. */ int get_exit_edge(int tri, const double& level, bool on_upper) const; - // Return the Triangulation object. - const Triangulation& get_triangulation() const; - // Return the z-value at the specified point index. const double& get_z(int point) const; @@ -468,8 +453,8 @@ class TriContourGenerator // Variables shared with python, always set. - PyObject* _triangulation; - PyArrayObject* _z; // double array (npoints). + Triangulation& _triangulation; + CoordinateArray _z; // double array (npoints). // Variables internal to C++ only. typedef std::vector InteriorVisited; // Size 2*ntri @@ -514,21 +499,20 @@ class TriContourGenerator class TrapezoidMapTriFinder { public: + typedef Triangulation::CoordinateArray CoordinateArray; + typedef numpy::array_view TriIndexArray; + /* Constructor. A separate call to initialize() is required to initialize * the object before use. - * triangulation: Triangulation to find triangles in. - * Triangulation reference count is not incremented in the constructor, but - * is decremented when no longer needed. */ - TrapezoidMapTriFinder(PyObject* triangulation); + * triangulation: Triangulation to find triangles in. */ + TrapezoidMapTriFinder(Triangulation& triangulation); ~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. - * Argument reference counts are not modified, return is a reference to a - * new array. */ - PyArrayObject* find_many(PyArrayObject* x, PyArrayObject* y); + /* Return an array of triangle indices. Takes 1D arrays x and y of + * point coordinates, and returns an array of the same size containing the + * indices of the triangles at those points. */ + TriIndexArray find_many(const CoordinateArray& x, const CoordinateArray& y); /* Return a reference to a new python list containing the following * statistics about the tree: @@ -784,13 +768,10 @@ class TrapezoidMapTriFinder 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. - PyObject* _triangulation; + Triangulation& _triangulation; // Variables internal to C++ only. Point* _points; // Array of all points in triangulation plus corners of diff --git a/lib/matplotlib/tri/_tri_wrapper.cpp b/lib/matplotlib/tri/_tri_wrapper.cpp index ca20960db5c2..b06c95d6a12e 100644 --- a/lib/matplotlib/tri/_tri_wrapper.cpp +++ b/lib/matplotlib/tri/_tri_wrapper.cpp @@ -1,10 +1,16 @@ -#include "_tri_wrapper.h" +#include "_tri.h" #include "src/mplutils.h" #include "src/py_exceptions.h" /* Triangulation */ +typedef struct +{ + PyObject_HEAD; + Triangulation* ptr; +} PyTriangulation; + static PyTypeObject PyTriangulationType; static PyObject* PyTriangulation_new(PyTypeObject* type, PyObject* args, PyObject* kwds) @@ -24,103 +30,64 @@ const char* PyTriangulation_init__doc__ = static int PyTriangulation_init(PyTriangulation* self, PyObject* args, PyObject* kwds) { - PyObject* x_arg; - PyObject* y_arg ; - PyObject* triangles_arg; - PyObject* mask_arg; - PyObject* edges_arg; - PyObject* neighbors_arg; - if (!PyArg_ParseTuple(args, "OOOOOO", &x_arg, &y_arg, &triangles_arg, - &mask_arg, &edges_arg, &neighbors_arg)) { + Triangulation::CoordinateArray x, y; + Triangulation::TriangleArray triangles; + Triangulation::MaskArray mask; + Triangulation::EdgeArray edges; + Triangulation::NeighborArray neighbors; + + if (!PyArg_ParseTuple(args, + "O&O&O&O&O&O&", + &x.converter, &x, + &y.converter, &y, + &triangles.converter, &triangles, + &mask.converter, &mask, + &edges.converter, &edges, + &neighbors.converter, &neighbors)) { return -1; } // x and y. - PyArrayObject* x = (PyArrayObject*)PyArray_ContiguousFromObject( - x_arg, NPY_DOUBLE, 1, 1); - PyArrayObject* y = (PyArrayObject*)PyArray_ContiguousFromObject( - y_arg, NPY_DOUBLE, 1, 1); - if (x == 0 || y == 0 || PyArray_DIM(x,0) != PyArray_DIM(y,0)) { - Py_XDECREF(x); - Py_XDECREF(y); + if (x.empty() || y.empty() || x.dim(0) != y.dim(0)) { PyErr_SetString(PyExc_ValueError, "x and y must be 1D arrays of the same length"); } // triangles. - PyArrayObject* triangles = (PyArrayObject*)PyArray_ContiguousFromObject( - triangles_arg, NPY_INT, 2, 2); - if (triangles == 0 || PyArray_DIM(triangles,1) != 3) { - Py_XDECREF(x); - Py_XDECREF(y); - Py_XDECREF(triangles); + if (triangles.empty() || triangles.dim(1) != 3) { PyErr_SetString(PyExc_ValueError, "triangles must be a 2D array of shape (?,3)"); } // Optional mask. - PyArrayObject* mask = 0; - if (mask_arg != 0 && mask_arg != Py_None) - { - mask = (PyArrayObject*)PyArray_ContiguousFromObject( - mask_arg, NPY_BOOL, 1, 1); - if (mask == 0 || PyArray_DIM(mask,0) != PyArray_DIM(triangles,0)) { - Py_XDECREF(x); - Py_XDECREF(y); - Py_XDECREF(triangles); - Py_XDECREF(mask); - PyErr_SetString(PyExc_ValueError, - "mask must be a 1D array with the same length as the triangles array"); - } + if (!mask.empty() && mask.dim(0) != triangles.dim(0)) { + PyErr_SetString(PyExc_ValueError, + "mask must be a 1D array with the same length as the triangles array"); } // Optional edges. - PyArrayObject* edges = 0; - if (edges_arg != 0 && edges_arg != Py_None) - { - edges = (PyArrayObject*)PyArray_ContiguousFromObject( - edges_arg, NPY_INT, 2, 2); - if (edges == 0 || PyArray_DIM(edges,1) != 2) { - Py_XDECREF(x); - Py_XDECREF(y); - Py_XDECREF(triangles); - Py_XDECREF(mask); - Py_XDECREF(edges); - PyErr_SetString(PyExc_ValueError, - "edges must be a 2D array with shape (?,2)"); - } + if (!edges.empty() && edges.dim(1) != 2) { + PyErr_SetString(PyExc_ValueError, + "edges must be a 2D array with shape (?,2)"); } // Optional neighbors. - PyArrayObject* neighbors = 0; - if (neighbors_arg != 0 && neighbors_arg != Py_None) - { - neighbors = (PyArrayObject*)PyArray_ContiguousFromObject( - neighbors_arg, NPY_INT, 2, 2); - if (neighbors == 0 || - PyArray_DIM(neighbors,0) != PyArray_DIM(triangles,0) || - PyArray_DIM(neighbors,1) != PyArray_DIM(triangles,1)) { - Py_XDECREF(x); - Py_XDECREF(y); - Py_XDECREF(triangles); - Py_XDECREF(mask); - Py_XDECREF(edges); - Py_XDECREF(neighbors); - PyErr_SetString(PyExc_ValueError, - "neighbors must be a 2D array with the same shape as the triangles array"); - } + if (!neighbors.empty() && (neighbors.dim(0) != triangles.dim(0) || + neighbors.dim(1) != triangles.dim(1))) { + PyErr_SetString(PyExc_ValueError, + "neighbors must be a 2D array with the same shape as the triangles array"); } CALL_CPP_INIT("Triangulation", - (self->ptr = new Triangulation(x, y, triangles, mask, edges, - neighbors))); + (self->ptr = new Triangulation(x, y, triangles, mask, + edges, neighbors))); return 0; } static void PyTriangulation_dealloc(PyTriangulation* self) { delete self->ptr; - Py_TYPE(self)->tp_free((PyObject *)self); + Py_TYPE(self)->tp_free((PyObject*)self); } const char* PyTriangulation_calculate_plane_coefficients__doc__ = @@ -130,29 +97,21 @@ const char* PyTriangulation_calculate_plane_coefficients__doc__ = static PyObject* PyTriangulation_calculate_plane_coefficients(PyTriangulation* self, PyObject* args, PyObject* kwds) { - PyObject* z_arg; - if (!PyArg_ParseTuple(args, "O:calculate_plane_coefficients", &z_arg)) { + Triangulation::CoordinateArray z; + if (!PyArg_ParseTuple(args, "O&:calculate_plane_coefficients", + &z.converter, &z)) { return NULL; } - PyArrayObject* z = (PyArrayObject*)PyArray_ContiguousFromObject( - z_arg, NPY_DOUBLE, 1, 1); - if (z == 0 || PyArray_DIM(z,0) != self->ptr->get_npoints()) { - Py_XDECREF(z); + if (z.empty() || z.dim(0) != self->ptr->get_npoints()) { PyErr_SetString(PyExc_ValueError, "z array must have same length as triangulation x and y arrays"); } - npy_intp dims[2] = {self->ptr->get_ntri(), 3}; - PyArrayObject* result = (PyArrayObject*)PyArray_SimpleNew( - 2, dims, NPY_DOUBLE); - - CALL_CPP_CLEANUP("calculate_plane_coefficients", - (self->ptr->calculate_plane_coefficients(z, result)), - Py_XDECREF(z); Py_XDECREF(result)); - - Py_XDECREF(z); - return (PyObject*)result; + Triangulation::TwoCoordinateArray result; + CALL_CPP("calculate_plane_coefficients", + (result = self->ptr->calculate_plane_coefficients(z))); + return result.pyobj(); } const char* PyTriangulation_get_edges__doc__ = @@ -162,10 +121,14 @@ const char* PyTriangulation_get_edges__doc__ = static PyObject* PyTriangulation_get_edges(PyTriangulation* self, PyObject* args, PyObject* kwds) { - PyArrayObject* result; - CALL_CPP("get_edges", (result = self->ptr->get_edges())); - Py_XINCREF(result); - return (PyObject*)result; + Triangulation::EdgeArray* result; + CALL_CPP("get_edges", (result = &self->ptr->get_edges())); + + if (result->empty()) { + Py_RETURN_NONE; + } + else + return result->pyobj(); } const char* PyTriangulation_get_neighbors__doc__ = @@ -175,10 +138,14 @@ const char* PyTriangulation_get_neighbors__doc__ = static PyObject* PyTriangulation_get_neighbors(PyTriangulation* self, PyObject* args, PyObject* kwds) { - PyArrayObject* result; - CALL_CPP("get_neighbors", (result = self->ptr->get_neighbors())); - Py_XINCREF(result); - return (PyObject*)result; + Triangulation::NeighborArray* result; + CALL_CPP("get_neighbors", (result = &self->ptr->get_neighbors())); + + if (result->empty()) { + Py_RETURN_NONE; + } + else + return result->pyobj(); } const char* PyTriangulation_set_mask__doc__ = @@ -188,26 +155,18 @@ const char* PyTriangulation_set_mask__doc__ = static PyObject* PyTriangulation_set_mask(PyTriangulation* self, PyObject* args, PyObject* kwds) { - PyObject* mask_arg; - if (!PyArg_ParseTuple(args, "O:set_mask", &mask_arg)) { + Triangulation::MaskArray mask; + + if (!PyArg_ParseTuple(args, "O&:set_mask", &mask.converter, &mask)) { return NULL; } - // Optional mask. - PyArrayObject* mask = 0; - if (mask_arg != 0 && mask_arg != Py_None) - { - mask = (PyArrayObject*)PyArray_ContiguousFromObject( - mask_arg, NPY_BOOL, 1, 1); - if (mask == 0 || PyArray_DIM(mask,0) != self->ptr->get_ntri()) { - Py_XDECREF(mask); - PyErr_SetString(PyExc_ValueError, - "mask must be a 1D array with the same length as the triangles array"); - } + if (!mask.empty() && mask.dim(0) != self->ptr->get_ntri()) { + PyErr_SetString(PyExc_ValueError, + "mask must be a 1D array with the same length as the triangles array"); } CALL_CPP("set_mask", (self->ptr->set_mask(mask))); - Py_RETURN_NONE; } @@ -245,6 +204,13 @@ static PyTypeObject* PyTriangulation_init_type(PyObject* m, PyTypeObject* type) /* TriContourGenerator */ +typedef struct +{ + PyObject_HEAD; + TriContourGenerator* ptr; + PyTriangulation* py_triangulation; +} PyTriContourGenerator; + static PyTypeObject PyTriContourGeneratorType; static PyObject* PyTriContourGenerator_new(PyTypeObject* type, PyObject* args, PyObject* kwds) @@ -252,6 +218,7 @@ static PyObject* PyTriContourGenerator_new(PyTypeObject* type, PyObject* args, P PyTriContourGenerator* self; self = (PyTriContourGenerator*)type->tp_alloc(type, 0); self->ptr = NULL; + self->py_triangulation = NULL; return (PyObject*)self; } @@ -264,34 +231,34 @@ const char* PyTriContourGenerator_init__doc__ = static int PyTriContourGenerator_init(PyTriContourGenerator* self, PyObject* args, PyObject* kwds) { - PyObject* triangulation; - PyObject* z_arg; - if (!PyArg_ParseTuple(args, "O!O", &PyTriangulationType, &triangulation, - &z_arg)) { + PyObject* triangulation_arg; + TriContourGenerator::CoordinateArray z; + + if (!PyArg_ParseTuple(args, "O!O&", + &PyTriangulationType, &triangulation_arg, + &z.converter, &z)) { return -1; } - Py_INCREF(triangulation); - int npoints = ((PyTriangulation*)triangulation)->ptr->get_npoints(); + PyTriangulation* py_triangulation = (PyTriangulation*)triangulation_arg; + Py_INCREF(py_triangulation); + self->py_triangulation = py_triangulation; + Triangulation& triangulation = *(py_triangulation->ptr); - PyArrayObject* z = (PyArrayObject*)PyArray_ContiguousFromObject( - z_arg, NPY_DOUBLE, 1, 1); - if (z == 0 && PyArray_DIM(z,0) != npoints) { - Py_DECREF(triangulation); - Py_XDECREF(z); + if (z.empty() || z.dim(0) != triangulation.get_npoints()) { PyErr_SetString(PyExc_ValueError, "z must be a 1D array with the same length as the x and y arrays"); } CALL_CPP_INIT("TriContourGenerator", (self->ptr = new TriContourGenerator(triangulation, z))); - return 0; } static void PyTriContourGenerator_dealloc(PyTriContourGenerator* self) { delete self->ptr; + Py_XDECREF(self->py_triangulation); Py_TYPE(self)->tp_free((PyObject *)self); } @@ -363,6 +330,13 @@ static PyTypeObject* PyTriContourGenerator_init_type(PyObject* m, PyTypeObject* /* TrapezoidMapTriFinder */ +typedef struct +{ + PyObject_HEAD; + TrapezoidMapTriFinder* ptr; + PyTriangulation* py_triangulation; +} PyTrapezoidMapTriFinder; + static PyTypeObject PyTrapezoidMapTriFinderType; static PyObject* PyTrapezoidMapTriFinder_new(PyTypeObject* type, PyObject* args, PyObject* kwds) @@ -370,6 +344,7 @@ static PyObject* PyTrapezoidMapTriFinder_new(PyTypeObject* type, PyObject* args, PyTrapezoidMapTriFinder* self; self = (PyTrapezoidMapTriFinder*)type->tp_alloc(type, 0); self->ptr = NULL; + self->py_triangulation = NULL; return (PyObject*)self; } @@ -382,11 +357,16 @@ const char* PyTrapezoidMapTriFinder_init__doc__ = static int PyTrapezoidMapTriFinder_init(PyTrapezoidMapTriFinder* self, PyObject* args, PyObject* kwds) { - PyObject* triangulation; - if (!PyArg_ParseTuple(args, "O!", &PyTriangulationType, &triangulation)) { + PyObject* triangulation_arg; + if (!PyArg_ParseTuple(args, "O!", + &PyTriangulationType, &triangulation_arg)) { return -1; } - Py_INCREF(triangulation); + + PyTriangulation* py_triangulation = (PyTriangulation*)triangulation_arg; + Py_INCREF(py_triangulation); + self->py_triangulation = py_triangulation; + Triangulation& triangulation = *(py_triangulation->ptr); CALL_CPP_INIT("TrapezoidMapTriFinder", (self->ptr = new TrapezoidMapTriFinder(triangulation))); @@ -396,6 +376,7 @@ static int PyTrapezoidMapTriFinder_init(PyTrapezoidMapTriFinder* self, PyObject* static void PyTrapezoidMapTriFinder_dealloc(PyTrapezoidMapTriFinder* self) { delete self->ptr; + Py_XDECREF(self->py_triangulation); Py_TYPE(self)->tp_free((PyObject *)self); } @@ -406,37 +387,21 @@ const char* PyTrapezoidMapTriFinder_find_many__doc__ = static PyObject* PyTrapezoidMapTriFinder_find_many(PyTrapezoidMapTriFinder* self, PyObject* args, PyObject* kwds) { - PyObject* x_arg; - PyObject* y_arg; - if (!PyArg_ParseTuple(args, "OO:find_many", &x_arg, &y_arg)) { + TrapezoidMapTriFinder::CoordinateArray x, y; + if (!PyArg_ParseTuple(args, "O&O&:find_many", + &x.converter, &x, + &y.converter, &y)) { return NULL; } - PyArrayObject* x = (PyArrayObject*)PyArray_ContiguousFromObject( - x_arg, NPY_DOUBLE, 0, 0); - PyArrayObject* y = (PyArrayObject*)PyArray_ContiguousFromObject( - y_arg, NPY_DOUBLE, 0, 0); - - bool ok = (x != 0 && y != 0 && PyArray_NDIM(x) == PyArray_NDIM(y)); - int ndim = (x == 0 ? 0 : 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); + if (x.empty() || y.empty() || x.dim(0) != y.dim(0)) { PyErr_SetString(PyExc_ValueError, "x and y must be array_like with same shape"); } - PyArrayObject* result; - CALL_CPP_CLEANUP("find_many", - (result = self->ptr->find_many(x, y)), - Py_XDECREF(x); Py_XDECREF(y)); - - Py_XDECREF(x); - Py_XDECREF(y); - return (PyObject*)result; + TrapezoidMapTriFinder::TriIndexArray result; + CALL_CPP("find_many", (result = self->ptr->find_many(x, y))); + return result.pyobj(); } const char* PyTrapezoidMapTriFinder_get_tree_stats__doc__ = diff --git a/lib/matplotlib/tri/_tri_wrapper.h b/lib/matplotlib/tri/_tri_wrapper.h deleted file mode 100644 index 9cd9d8b56c9a..000000000000 --- a/lib/matplotlib/tri/_tri_wrapper.h +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef _TRI_WRAPPER_H -#define _TRI_WRAPPER_H - -#include "_tri.h" - -typedef struct -{ - PyObject_HEAD; - Triangulation* ptr; -} PyTriangulation; - -typedef struct -{ - PyObject_HEAD; - TriContourGenerator* ptr; -} PyTriContourGenerator; - -typedef struct -{ - PyObject_HEAD; - TrapezoidMapTriFinder* ptr; -} PyTrapezoidMapTriFinder; - -#endif // _TRI_WRAPPER_H diff --git a/lib/matplotlib/tri/trifinder.py b/lib/matplotlib/tri/trifinder.py index 3918e99d3b5f..247d061ea1d0 100644 --- a/lib/matplotlib/tri/trifinder.py +++ b/lib/matplotlib/tri/trifinder.py @@ -5,6 +5,7 @@ from matplotlib.tri import Triangulation import matplotlib._tri as _tri +import numpy as np class TriFinder(object): @@ -54,8 +55,15 @@ def __call__(self, x, y): Returns integer array with the same shape and *x* and *y*. """ - # C++ checks arguments are OK. - return self._cpp_trifinder.find_many(x, y) + 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 array-like with the same shape") + + # C++ does the heavy lifting, and expects 1D arrays. + indices = self._cpp_trifinder.find_many(x.ravel(), y.ravel()) + indices.shape = x.shape + return indices def _get_tree_stats(self): """ diff --git a/src/numpy_cpp.h b/src/numpy_cpp.h index 966ae599f824..8d74d26b0a6c 100644 --- a/src/numpy_cpp.h +++ b/src/numpy_cpp.h @@ -35,10 +35,17 @@ namespace numpy template struct type_num_of; -// This is dodgy - need sizeof(bool) == 1 consistently for this to be valid... -// template <> struct type_num_of { -// enum {value = NPY_BOOL}; -//}; +/* Be careful with bool arrays as python has sizeof(npy_bool) == 1, but it is + * not always the case that sizeof(bool) == 1. Using the array_view_accessors + * is always fine regardless of sizeof(bool), so do this rather than using + * array.data() and pointer arithmetic which will not work correctly if + * sizeof(bool) != 1. */ +template <> struct type_num_of +{ + enum { + value = NPY_BOOL + }; +}; template <> struct type_num_of { @@ -394,6 +401,20 @@ class array_view : public detail::array_view_accessors Py_XDECREF(m_arr); } + const array_view& operator=(const array_view &other) + { + if (this != &other) + { + Py_XDECREF(m_arr); + m_arr = other.m_arr; + Py_XINCREF(m_arr); + m_data = other.m_data; + m_shape = other.m_shape; + m_strides = other.m_strides; + } + return *this; + } + int set(PyObject *arr, bool contiguous = false) { PyArrayObject *tmp; @@ -440,7 +461,7 @@ class array_view : public detail::array_view_accessors return 1; } - npy_intp dim(size_t i) + npy_intp dim(size_t i) const { if (i > ND) { return 0; @@ -448,11 +469,17 @@ class array_view : public detail::array_view_accessors return m_shape[i]; } - size_t size() + size_t size() const { return (size_t)dim(0); } + bool empty() const + { + return size() == 0; + } + + // Do not use this for array_view. See comment near top of file. T *data() { return (T *)m_data;