Index: lib/matplotlib/tri/tricontour.py
===================================================================
--- lib/matplotlib/tri/tricontour.py	(revision 0)
+++ lib/matplotlib/tri/tricontour.py	(revision 0)
@@ -0,0 +1,288 @@
+from matplotlib.contour import ContourSet
+from matplotlib.tri.triangulation import Triangulation
+import matplotlib._tri as _tri
+import numpy as np
+
+
+class TriContourSet(ContourSet):
+    """
+    Create and store a set of contour lines or filled regions for
+    a triangular grid.
+
+    User-callable method: clabel
+
+    Useful attributes:
+      ax:
+        the axes object in which the contours are drawn
+      collections:
+        a silent_list of LineCollections or PolyCollections
+      levels:
+        contour levels
+      layers:
+        same as levels for line contours; half-way between
+        levels for filled contours.  See _process_colors method.
+    """
+    def __init__(self, ax, *args, **kwargs):
+        """
+        Draw triangular grid contour lines or filled regions,
+        depending on whether keyword arg 'filled' is False
+        (default) or True.
+
+        The first argument of the initializer must be an axes
+        object.  The remaining arguments and keyword arguments
+        are described in TriContourSet.tricontour_doc.
+        """
+        ContourSet.__init__(self, ax, *args, **kwargs)
+
+    def _process_args(self, *args, **kwargs):
+        """
+        Process args and kwargs.
+        """
+        if isinstance(args[0], TriContourSet):
+            C = args[0].cppContourGenerator
+            if self.levels is None:
+                self.levels = args[0].levels
+        else:
+            tri, z = self._contour_args(args, kwargs)
+            C = _tri.TriContourGenerator(tri.get_cpp_triangulation(), z)
+            x0 = tri.x.min()
+            x1 = tri.x.max()
+            y0 = tri.y.min()
+            y1 = tri.y.max()
+            self.ax.update_datalim([(x0,y0), (x1,y1)])
+            self.ax.autoscale_view()
+
+        self.cppContourGenerator = C
+
+    def _get_allsegs_and_allkinds(self):
+        """
+        Create and return allsegs and allkinds by calling underlying C code.
+        """
+        allsegs = []
+        if self.filled:
+            lowers, uppers = self._get_lowers_and_uppers()
+            allkinds = []
+            for lower, upper in zip(lowers, uppers):
+                segs, kinds = self.cppContourGenerator.create_filled_contour(
+                                                                 lower, upper)
+                allsegs.append([segs])
+                allkinds.append([kinds])
+        else:
+            allkinds = None
+            for level in self.levels:
+                segs = self.cppContourGenerator.create_contour(level)
+                allsegs.append(segs)
+        return allsegs, allkinds
+
+    def _contour_args(self, args, kwargs):
+        if self.filled: fn = 'contourf'
+        else:           fn = 'contour'
+        tri, args, kwargs = \
+            Triangulation.get_from_args_and_kwargs(*args, **kwargs)
+        z = np.asarray(args[0])
+        if z.shape != tri.x.shape:
+            raise ValueError('z array must have same length as triangulation x'
+                             'and y arrays')
+        self.zmax = z.max()
+        self.zmin = z.min()
+        if self.logscale and self.zmin <= 0:
+            raise ValueError('Cannot %s log of negative values.' % fn)
+        self._contour_level_args(z, args[1:])
+        return (tri, z)
+
+    tricontour_doc = """
+        :func:`~matplotlib.pyplot.tricontour` and
+        :func:`~matplotlib.pyplot.tricontourf` draw contour lines and
+        filled contours, respectively.  Except as noted, function
+        signatures and return values are the same for both versions.
+
+        The triangulation can be specified in one of two ways; either::
+
+          tricontour(triangulation, ...)
+
+        where triangulation is a :class:`~matplotlib.tri.Triangulation`
+        object, or
+
+        ::
+
+          tricontour(x, y, ...)
+          tricontour(x, y, triangles, ...)
+          tricontour(x, y, triangles=triangles, ...)
+          tricontour(x, y, mask, ...)
+          tricontour(x, y, mask=mask, ...)
+          tricontour(x, y, triangles, mask, ...)
+          tricontour(x, y, triangles, mask=mask, ...)
+
+        in which case a Triangulation object will be created.  See
+        :class:`~matplotlib.tri.Triangulation` for a explanation of
+        these possibilities.
+
+        The remaining arguments may be::
+
+          tricontour(..., Z)
+
+        where *Z* is the array of values to contour, one per point
+        in the triangulation.  The level values are chosen
+        automatically.
+
+        ::
+
+          tricontour(..., Z, N)
+
+        contour *N* automatically-chosen levels.
+
+        ::
+
+          tricontour(..., Z, V)
+
+        draw contour lines at the values specified in sequence *V*
+
+        ::
+
+          tricontourf(..., Z, V)
+
+        fill the (len(*V*)-1) regions between the values in *V*
+
+        ::
+
+          tricontour(Z, **kwargs)
+
+        Use keyword args to control colors, linewidth, origin, cmap ... see
+        below for more details.
+
+        ``C = tricontour(...)`` returns a
+        :class:`~matplotlib.contour.TriContourSet` object.
+
+        Optional keyword arguments:
+
+          *colors*: [ None | string | (mpl_colors) ]
+            If *None*, the colormap specified by cmap will be used.
+
+            If a string, like 'r' or 'red', all levels will be plotted in this
+            color.
+
+            If a tuple of matplotlib color args (string, float, rgb, etc),
+            different levels will be plotted in different colors in the order
+            specified.
+
+          *alpha*: float
+            The alpha blending value
+
+          *cmap*: [ None | Colormap ]
+            A cm :class:`~matplotlib.cm.Colormap` instance or
+            *None*. If *cmap* is *None* and *colors* is *None*, a
+            default Colormap is used.
+
+          *norm*: [ None | Normalize ]
+            A :class:`matplotlib.colors.Normalize` instance for
+            scaling data values to colors. If *norm* is *None* and
+            *colors* is *None*, the default linear scaling is used.
+
+          *levels* [level0, level1, ..., leveln]
+            A list of floating point numbers indicating the level
+            curves to draw; eg to draw just the zero contour pass
+            ``levels=[0]``
+
+          *origin*: [ None | 'upper' | 'lower' | 'image' ]
+            If *None*, the first value of *Z* will correspond to the
+            lower left corner, location (0,0). If 'image', the rc
+            value for ``image.origin`` will be used.
+
+            This keyword is not active if *X* and *Y* are specified in
+            the call to contour.
+
+          *extent*: [ None | (x0,x1,y0,y1) ]
+
+            If *origin* is not *None*, then *extent* is interpreted as
+            in :func:`matplotlib.pyplot.imshow`: it gives the outer
+            pixel boundaries. In this case, the position of Z[0,0]
+            is the center of the pixel, not a corner. If *origin* is
+            *None*, then (*x0*, *y0*) is the position of Z[0,0], and
+            (*x1*, *y1*) is the position of Z[-1,-1].
+
+            This keyword is not active if *X* and *Y* are specified in
+            the call to contour.
+
+          *locator*: [ None | ticker.Locator subclass ]
+            If *locator* is None, the default
+            :class:`~matplotlib.ticker.MaxNLocator` is used. The
+            locator is used to determine the contour levels if they
+            are not given explicitly via the *V* argument.
+
+          *extend*: [ 'neither' | 'both' | 'min' | 'max' ]
+            Unless this is 'neither', contour levels are automatically
+            added to one or both ends of the range so that all data
+            are included. These added ranges are then mapped to the
+            special colormap values which default to the ends of the
+            colormap range, but can be set via
+            :meth:`matplotlib.colors.Colormap.set_under` and
+            :meth:`matplotlib.colors.Colormap.set_over` methods.
+
+          *xunits*, *yunits*: [ None | registered units ]
+            Override axis units by specifying an instance of a
+            :class:`matplotlib.units.ConversionInterface`.
+
+
+        tricontour-only keyword arguments:
+
+          *linewidths*: [ None | number | tuple of numbers ]
+            If *linewidths* is *None*, the default width in
+            ``lines.linewidth`` in ``matplotlibrc`` is used.
+
+            If a number, all levels will be plotted with this linewidth.
+
+            If a tuple, different levels will be plotted with different
+            linewidths in the order specified
+
+          *linestyles*: [None | 'solid' | 'dashed' | 'dashdot' | 'dotted' ]
+            If *linestyles* is *None*, the 'solid' is used.
+
+            *linestyles* can also be an iterable of the above strings
+            specifying a set of linestyles to be used. If this
+            iterable is shorter than the number of contour levels
+            it will be repeated as necessary.
+
+            If contour is using a monochrome colormap and the contour
+            level is less than 0, then the linestyle specified
+            in ``contour.negative_linestyle`` in ``matplotlibrc``
+            will be used.
+
+        tricontourf-only keyword arguments:
+
+          *antialiased*: [ True | False ]
+            enable antialiasing
+
+          *nchunk*: [ 0 | integer ]
+            If 0, no subdivision of the domain. Specify a positive integer to
+            divide the domain into subdomains of roughly *nchunk* by *nchunk*
+            points. This may never actually be advantageous, so this option may
+            be removed. Chunking introduces artifacts at the chunk boundaries
+            unless *antialiased* is *False*.
+
+        Note: tricontourf fills intervals that are closed at the top; that
+        is, for boundaries *z1* and *z2*, the filled region is::
+
+            z1 < z <= z2
+
+        There is one exception: if the lowest boundary coincides with
+        the minimum value of the *z* array, then that minimum value
+        will be included in the lowest interval.
+
+        **Examples:**
+
+        .. plot:: mpl_examples/pylab_examples/tricontour_demo.py
+        """
+
+
+def tricontour(ax, *args, **kwargs):
+    if not ax._hold: ax.cla()
+    kwargs['filled'] = False
+    return TriContourSet(ax, *args, **kwargs)
+tricontour.__doc__ = TriContourSet.tricontour_doc
+
+
+def tricontourf(ax, *args, **kwargs):
+    if not ax._hold: ax.cla()
+    kwargs['filled'] = True
+    return TriContourSet(ax, *args, **kwargs)
+tricontourf.__doc__ = TriContourSet.tricontour_doc
Index: lib/matplotlib/tri/_tri.h
===================================================================
--- lib/matplotlib/tri/_tri.h	(revision 0)
+++ lib/matplotlib/tri/_tri.h	(revision 0)
@@ -0,0 +1,473 @@
+/*
+ * Unstructured triangular grid functions, particularly contouring.
+ *
+ * There are two main classes: Triangulation and TriContourGenerator.
+ *
+ * Triangulation
+ * -------------
+ * Triangulation is an unstructured triangular grid with npoints and ntri
+ * triangles.  It consists of point x and y coordinates, and information about
+ * the triangulation stored in an integer array of shape (ntri,3) called
+ * triangles.  Each triangle is represented by three point indices (in the
+ * range 0 to npoints-1) that comprise the triangle, ordered anticlockwise.
+ * There is an optional mask of length ntri which can be used to mask out
+ * triangles and has the same result as removing those triangles from the
+ * 'triangles' array.
+ *
+ * A particular edge of a triangulation is termed a TriEdge, which is a
+ * triangle index and an edge index in the range 0 to 2.  TriEdge(tri,edge)
+ * refers to the edge that starts at point index triangles(tri,edge) and ends
+ * at point index triangles(tri,(edge+1)%3).
+ *
+ * Various derived fields are calculated when they are first needed.  The
+ * triangle connectivity is stored in a neighbors array of shape (ntri,3) such
+ * that neighbors(tri,edge) is the index of the triangle that adjoins the
+ * TriEdge(tri,edge), or -1 if there is no such neighbor.
+ *
+ * A triangulation has one or more boundaries, each of which is a 1D array of
+ * the TriEdges that comprise the boundary, in order following the boundary
+ * with non-masked triangles on the left.
+ *
+ * TriContourGenerator
+ * -------------------
+ * A TriContourGenerator generates contours for a particular Triangulation.
+ * The process followed is different for non-filled and filled contours, with
+ * one and two contour levels respectively.  In both cases boundary contour
+ * lines are found first, then interior lines.
+ *
+ * Boundary lines start and end on a boundary.  They are found by traversing
+ * the triangulation boundary edges until a suitable start point is found, and
+ * then the contour line is followed across the interior of the triangulation
+ * until it ends on another boundary edge.  For a non-filled contour this
+ * completes a line, whereas a filled contour continues by following the
+ * boundary around until either another boundary start point is found or the
+ * start of the contour line is reached.  Filled contour generation stores
+ * boolean flags to indicate which boundary edges have already been traversed
+ * so that they are not dealt with twice.  Similar flags are used to indicate
+ * which triangles have been used when following interior lines.
+ *
+ * Interior lines do not intersect any boundaries.  They are found by
+ * traversing all triangles that have not yet been visited until a suitable
+ * starting point is found, and then the contour line is followed across the
+ * interior of the triangulation until it returns to the start point.  For
+ * filled contours this process is repeated for both lower and upper contour
+ * levels, and the direction of traversal is reversed for upper contours.
+ *
+ * Working out in which direction a contour line leaves a triangle uses the
+ * a lookup table.  A triangle has three points, each of which has a z-value
+ * which is either less than the contour level or not.  Hence there are 8
+ * configurations to deal with, 2 of which do not have a contour line (all
+ * points below or above (including the same as) the contour level) and 6 that
+ * do.  See the function get_exit_edge for details.
+ */
+#ifndef _TRI_H
+#define _TRI_H
+
+#include "CXX/Extensions.hxx"
+#include "CXX/Objects.hxx"
+#include "numpy/arrayobject.h"
+
+#include <map>
+#include <vector>
+
+
+
+/* An edge of a triangle consisting of an triangle index in the range 0 to
+ * ntri-1 and an edge index in the range 0 to 2.  Edge i goes from the
+ * triangle's point i to point (i+1)%3. */
+struct TriEdge
+{
+  TriEdge();
+  TriEdge(int tri_, int edge_);
+  bool operator<(const TriEdge& other) const;
+  bool operator==(const TriEdge& other) const;
+  bool operator!=(const TriEdge& other) const;
+  friend std::ostream& operator<<(std::ostream& os, const TriEdge& tri_edge);
+
+  int tri, edge;
+};
+
+// 2D point with x,y coordinates.
+struct XY
+{
+  XY();
+  XY(const double& x_, const double& y_);
+  double angle() const;           // Angle in radians with respect to x-axis.
+  double cross_z(const XY& other) const;     // z-component of cross product.
+  bool operator==(const XY& other) const;
+  bool operator!=(const XY& other) const;
+  bool operator<(const XY& other) const;
+  XY operator*(const double& multiplier) const;
+  XY operator+(const XY& other) const;
+  XY operator-(const XY& other) const;
+  friend std::ostream& operator<<(std::ostream& os, const XY& xy);
+
+  double x, y;
+};
+
+// 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);
+
+  // Consider these member variables read-only.
+  bool empty;
+  XY lower, upper;
+};
+
+/* A single line of a contour, which may be a closed line loop or an open line
+ * strip.  Identical adjacent points are avoided using insert_unique() and
+ * push_back(), and a closed line loop should also not have identical first and
+ * last points. */
+class ContourLine : public std::vector<XY>
+{
+public:
+    ContourLine();
+    void insert_unique(iterator pos, const XY& point);
+    void push_back(const XY& point);
+    void write() const;
+};
+
+// A Contour is a collection of zero or more ContourLines.
+typedef std::vector<ContourLine> Contour;
+
+// Debug contour writing function.
+void write_contour(const Contour& contour);
+
+
+
+
+/* Triangulation with npoints points and ntri triangles.  Derived fields are
+ * calculated when they are first needed. */
+class Triangulation : public Py::PythonExtension<Triangulation>
+{
+public:
+    /* A single boundary is a vector of the TriEdges that make up that boundary
+     * following it around with unmasked triangles on the left. */
+    typedef std::vector<TriEdge> Boundary;
+    typedef std::vector<Boundary> Boundaries;
+
+    /* Constructor with optional mask, edges and neighbors.  The latter two
+     * are calculated when first needed.
+     *   x: double array of shape (npoints) of points' x-coordinates.
+     *   y: double array of shape (npoints) of points' y-coordinates.
+     *   triangles: int array of shape (ntri,3) of triangle point indices.
+     *              Those ordered clockwise are changed to be anticlockwise.
+     *   mask: Optional bool array of shape (ntri) indicating which triangles
+     *         are masked.
+     *   edges: Optional int array of shape (?,2) of start and end point
+     *          indices, each edge (start,end and end,start) appearing only
+     *          once.
+     *   neighbors: Optional int array of shape (ntri,3) indicating which
+     *              triangles are the neighbors of which TriEdges, or -1 if
+     *              there is no such neighbor. */
+    Triangulation(PyArrayObject* x,
+                  PyArrayObject* y,
+                  PyArrayObject* triangles,
+                  PyArrayObject* mask,
+                  PyArrayObject* edges,
+                  PyArrayObject* neighbors);
+
+    virtual ~Triangulation();
+
+    // Return the boundaries collection, creating it if necessary.
+    const Boundaries& get_boundaries() const;
+
+    // Return which boundary and boundary edge the specified TriEdge is.
+    void get_boundary_edge(const TriEdge& triEdge,
+                           int& boundary,
+                           int& edge) const;
+
+    // Return the edges array, creating it if necessary.
+    Py::Object get_edges();
+
+    /* Return the triangle index of the neighbor of the specified triangle
+     * edge. */
+    int get_neighbor(int tri, int edge) const;
+
+    /* Return the TriEdge that is the neighbor of the specified triangle edge,
+     * or TriEdge(-1,-1) if there is no such neighbor. */
+    TriEdge get_neighbor_edge(int tri, int edge) const;
+
+    // Return the neighbors array, creating it if necessary.
+    Py::Object get_neighbors();
+
+    // Return the number of points in this triangulation.
+    int get_npoints() const;
+
+    // Return the number of triangles in this triangulation.
+    int get_ntri() const;
+
+    /* Return the index of the point that is at the start of the specified
+     * triangle edge. */
+    int get_triangle_point(int tri, int edge) const;
+    int get_triangle_point(const TriEdge& tri_edge) const;
+
+    // Return the coordinates of the specified point index.
+    XY get_point_coords(int point) const;
+
+    // CXX initialisation function.
+    static void init_type();
+
+    // Indicates if the specified triangle is masked or not.
+    bool is_masked(int tri) const;
+
+    /* Set or clear the mask array.  Clears various derived fields so they are
+     * recalculated when next needed.
+     *   args[0]: mask, either Py::None or boolean array of shape (ntri). */
+    Py::Object set_mask(const Py::Tuple &args);
+
+    // Debug function to write boundaries.
+    void write_boundaries() const;
+
+private:
+    // An edge of a triangulation, composed of start and end point indices.
+    struct Edge
+    {
+        Edge() : start(-1), end(-1) {}
+        Edge(int start_, int end_) : start(start_), end(end_) {}
+        bool operator<(const Edge& other) const {
+            return start != other.start ? start < other.start : end < other.end;
+        }
+        int start, end;
+    };
+
+    /* An edge of a boundary of a triangulation, composed of a boundary index
+     * and an edge index within that boundary.  Used to index into the
+     * boundaries collection to obtain the corresponding TriEdge. */
+    struct BoundaryEdge
+    {
+        BoundaryEdge() : boundary(-1), edge(-1) {}
+        BoundaryEdge(int boundary_, int edge_)
+            : boundary(boundary_), edge(edge_) {}
+        int boundary, edge;
+    };
+
+    /* Calculate the boundaries collection.  Should normally be accessed via
+     * get_boundaries(), which will call this function if necessary. */
+    void calculate_boundaries();
+
+    /* Calculate the edges array.  Should normally be accessed via
+     * get_edges(), which will call this function if necessary. */
+    void calculate_edges();
+
+    /* Calculate the neighbors array. Should normally be accessed via
+     * get_neighbors(), which will call this function if necessary. */
+    void calculate_neighbors();
+
+    /* Correct each triangle so that the vertices are ordered in an
+     * anticlockwise manner. */
+    void correct_triangles();
+
+    /* Determine which edge index (0,1 or 2) the specified point index is in
+     * the specified triangle, or -1 if the point is not in the triangle. */
+    int get_edge_in_triangle(int tri, int point) const;
+
+    // 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,
+                               //     ordered anticlockwise.
+
+    // Variables shared with python, may be zero.
+    PyArrayObject* _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
+                               //     or -1 if no neighbor.
+
+    // Variables internal to C++ only.
+    Boundaries _boundaries;
+
+    // Map used to look up BoundaryEdges from TriEdges.  Normally accessed via
+    // get_boundary_edge().
+    typedef std::map<TriEdge, BoundaryEdge> TriEdgeToBoundaryMap;
+    TriEdgeToBoundaryMap _tri_edge_to_boundary_map;
+};
+
+
+
+
+// Contour generator for a triangulation.
+class TriContourGenerator : public Py::PythonExtension<TriContourGenerator>
+{
+public:
+    /* Constructor.
+     *   triangulation: Triangulation to generate contours for.
+     *   z: Double array of shape (npoints) of z-values at triangulation
+     *      points. */
+    TriContourGenerator(Py::Object triangulation,
+                        PyArrayObject* z);
+
+    virtual ~TriContourGenerator();
+
+    /* Create and return a non-filled contour.
+     *   level: Contour level.
+     * Returns python list [segs0, segs1, ...] where
+     *   segs0: double array of shape (?,2) of point coordinates of first
+     *   contour line, etc. */
+    Py::Object create_contour(const Py::Tuple &args);
+
+    /* Create and return a filled contour.
+     *   lower_level: Lower contour level.
+     *   upper_level: Upper contour level.
+     * Returns python tuple (segs, kinds) where
+     *   segs: double array of shape (n_points,2) of all point coordinates,
+     *   kinds: ubyte array of shape (n_points) of all point code types. */
+    Py::Object create_filled_contour(const Py::Tuple &args);
+
+    // CXX initialisation function.
+    static void init_type();
+
+private:
+    typedef Triangulation::Boundary Boundary;
+    typedef Triangulation::Boundaries Boundaries;
+
+    /* Clear visited flags.
+     *   include_boundaries: Whether to clear boundary flags or not, which are
+     *                       only used for filled contours. */
+    void clear_visited_flags(bool include_boundaries);
+
+    /* Convert a non-filled Contour from C++ to Python.
+     * Returns python list [segs0, segs1, ...] where
+     *   segs0: double array of shape (?,2) of point coordinates of first
+     *   contour line, etc. */
+    Py::Object contour_to_segs(const Contour& contour);
+
+    /* Convert a filled Contour from C++ to Python.
+     * Returns python tuple (segs, kinds) where
+     *   segs: double array of shape (n_points,2) of all point coordinates,
+     *   kinds: ubyte array of shape (n_points) of all point code types. */
+    Py::Object contour_to_segs_and_kinds(const Contour& contour);
+
+    /* Return the point on the specified TriEdge that intersects the specified
+     * level. */
+    XY edge_interp(int tri, int edge, const double& level);
+
+    /* Find and follow non-filled contour lines that start and end on a
+     * boundary of the Triangulation.
+     *   contour: Contour to add new lines to.
+     *   level: Contour level. */
+    void find_boundary_lines(Contour& contour,
+                             const double& level);
+
+    /* Find and follow filled contour lines at either of the specified contour
+     * levels that start and end of a boundary of the Triangulation.
+     *   contour: Contour to add new lines to.
+     *   lower_level: Lower contour level.
+     *   upper_level: Upper contour level. */
+    void find_boundary_lines_filled(Contour& contour,
+                                    const double& lower_level,
+                                    const double& upper_level);
+
+    /* Find and follow lines at the specified contour level that are
+     * completely in the interior of the Triangulation and hence do not
+     * intersect any boundary.
+     *   contour: Contour to add new lines to.
+     *   level: Contour level.
+     *   on_upper: Whether on upper or lower contour level.
+     *   filled: Whether contours are filled or not. */
+    void find_interior_lines(Contour& contour,
+                             const double& level,
+                             bool on_upper,
+                             bool filled);
+
+    /* Follow contour line around boundary of the Triangulation from the
+     * specified TriEdge to its end which can be on either the lower or upper
+     * levels.  Only used for filled contours.
+     *   contour_line: Contour line to append new points to.
+     *   tri_edge: On entry, TriEdge to start from.  On exit, TriEdge that is
+     *             finished on.
+     *   lower_level: Lower contour level.
+     *   upper_level: Upper contour level.
+     *   on_upper: Whether starts on upper level or not.
+     * Return true if finishes on upper level, false if lower. */
+    bool follow_boundary(ContourLine& contour_line,
+                         TriEdge& tri_edge,
+                         const double& lower_level,
+                         const double& upper_level,
+                         bool on_upper);
+
+    /* Follow contour line across interior of Triangulation.
+     *   contour_line: Contour line to append new points to.
+     *   tri_edge: On entry, TriEdge to start from.  On exit, TriEdge that is
+     *             finished on.
+     *   end_on_boundary: Whether this line ends on a boundary, or loops back
+     *                    upon itself.
+     *   level: Contour level to follow.
+     *   on_upper: Whether following upper or lower contour level. */
+    void follow_interior(ContourLine& contour_line,
+                         TriEdge& tri_edge,
+                         bool end_on_boundary,
+                         const double& level,
+                         bool on_upper);
+
+    // Return the Triangulation boundaries.
+    const Boundaries& get_boundaries() const;
+
+    /* Return the edge by which the a level leaves a particular triangle,
+     * which is 0, 1 or 2 if the contour passes through the triangle or -1
+     * otherwise.
+     *   tri: Triangle index.
+     *   level: Contour level to follow.
+     *   on_upper: Whether following upper or lower contour level. */
+    int get_exit_edge(int tri, const double& level, bool on_upper) const;
+
+    // Return the Triangulation object.
+    const Triangulation& get_triangulation() const;
+
+    // Return the z-value at the specified point index.
+    const double& get_z(int point) const;
+
+    /* Return the point at which the a level intersects the line connecting the
+     * two specified point indices. */
+    XY interp(int point1, int point2, const double& level) const;
+
+
+
+    // Variables shared with python, always set.
+    Py::Object _triangulation;
+    PyArrayObject* _z;        // double array (npoints).
+
+    // Variables internal to C++ only.
+    typedef std::vector<bool> InteriorVisited;    // Size 2*ntri
+    typedef std::vector<bool> BoundaryVisited;
+    typedef std::vector<BoundaryVisited> BoundariesVisited;
+    typedef std::vector<bool> BoundariesUsed;
+
+    InteriorVisited _interior_visited;
+    BoundariesVisited _boundaries_visited;  // Only used for filled contours.
+    BoundariesUsed _boundaries_used;        // Only used for filled contours.
+};
+
+
+
+
+// The extension module.
+class TriModule : public Py::ExtensionModule<TriModule>
+{
+public:
+    TriModule();
+
+private:
+    Py::Object new_triangulation(const Py::Tuple &args);
+    Py::Object new_tricontourgenerator(const Py::Tuple &args);
+};
+
+#endif
Index: lib/matplotlib/tri/tripcolor.py
===================================================================
--- lib/matplotlib/tri/tripcolor.py	(revision 0)
+++ lib/matplotlib/tri/tripcolor.py	(revision 0)
@@ -0,0 +1,100 @@
+from matplotlib.collections import PolyCollection
+from matplotlib.colors import Normalize
+from matplotlib.tri.triangulation import Triangulation
+import numpy as np
+
+def tripcolor(ax, *args, **kwargs):
+    """
+    Create a pseudocolor plot of an unstructured triangular grid to
+    the :class:`~matplotlib.axes.Axes`.
+
+    The triangulation can be specified in one of two ways; either::
+
+      tripcolor(triangulation, ...)
+
+    where triangulation is a :class:`~matplotlib.tri.Triangulation`
+    object, or
+
+    ::
+
+      tripcolor(x, y, ...)
+      tripcolor(x, y, triangles, ...)
+      tripcolor(x, y, triangles=triangles, ...)
+      tripcolor(x, y, mask, ...)
+      tripcolor(x, y, mask=mask, ...)
+      tripcolor(x, y, triangles, mask, ...)
+      tripcolor(x, y, triangles, mask=mask, ...)
+
+    in which case a Triangulation object will be created.  See
+    :class:`~matplotlib.tri.Triangulation` for a explanation of these
+    possibilities.
+
+    The next argument must be *C*, the array of color values, one per
+    point in the triangulation.  The colors used for each triangle
+    are from the mean C of the triangle's three points.
+
+    The remaining kwargs are the same as for
+    :meth:`~matplotlib.axes.Axes.pcolor`.
+
+    **Example:**
+
+        .. plot:: mpl_examples/pylab_examples/tripcolor_demo.py
+    """
+    if not ax._hold: ax.cla()
+
+    alpha = kwargs.pop('alpha', 1.0)
+    norm = kwargs.pop('norm', None)
+    cmap = kwargs.pop('cmap', None)
+    vmin = kwargs.pop('vmin', None)
+    vmax = kwargs.pop('vmax', None)
+    shading = kwargs.pop('shading', 'flat')
+
+    tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args, **kwargs)
+    x = tri.x
+    y = tri.y
+    triangles = tri.get_masked_triangles()
+
+    # Vertices of triangles.
+    verts = np.concatenate((x[triangles][...,np.newaxis],
+                            y[triangles][...,np.newaxis]), axis=2)
+
+    C = np.asarray(args[0])
+    if C.shape != x.shape:
+        raise ValueError('C array must have same length as triangulation x and'
+                         ' y arrays')
+
+    # Color values, one per triangle, mean of the 3 vertex color values.
+    C = C[triangles].mean(axis=1)
+
+    if shading == 'faceted':
+        edgecolors = (0,0,0,1),
+        linewidths = (0.25,)
+    else:
+        edgecolors = 'face'
+        linewidths = (1.0,)
+    kwargs.setdefault('edgecolors', edgecolors)
+    kwargs.setdefault('antialiaseds', (0,))
+    kwargs.setdefault('linewidths', linewidths)
+
+    collection = PolyCollection(verts, **kwargs)
+
+    collection.set_alpha(alpha)
+    collection.set_array(C)
+    if norm is not None: assert(isinstance(norm, Normalize))
+    collection.set_cmap(cmap)
+    collection.set_norm(norm)
+    if vmin is not None or vmax is not None:
+        collection.set_clim(vmin, vmax)
+    else:
+        collection.autoscale_None()
+    ax.grid(False)
+
+    minx = tri.x.min()
+    maxx = tri.x.max()
+    miny = tri.y.min()
+    maxy = tri.y.max()
+    corners = (minx, miny), (maxx, maxy)
+    ax.update_datalim( corners)
+    ax.autoscale_view()
+    ax.add_collection(collection)
+    return collection
Index: lib/matplotlib/tri/triplot.py
===================================================================
--- lib/matplotlib/tri/triplot.py	(revision 0)
+++ lib/matplotlib/tri/triplot.py	(revision 0)
@@ -0,0 +1,60 @@
+from matplotlib.tri.triangulation import Triangulation
+
+def triplot(ax, *args, **kwargs):
+    """
+    Draw a unstructured triangular grid as lines and/or markers to
+    the :class:`~matplotlib.axes.Axes`.
+
+    The triangulation to plot can be specified in one of two ways;
+    either::
+
+      triplot(triangulation, ...)
+
+    where triangulation is a :class:`~matplotlib.tri.Triangulation`
+    object, or
+
+    ::
+
+      triplot(x, y, ...)
+      triplot(x, y, triangles, ...)
+      triplot(x, y, triangles=triangles, ...)
+      triplot(x, y, mask, ...)
+      triplot(x, y, mask=mask, ...)
+      triplot(x, y, triangles, mask, ...)
+      triplot(x, y, triangles, mask=mask, ...)
+
+    in which case a Triangulation object will be created.  See
+    :class:`~matplotlib.tri.Triangulation` for a explanation of these
+    possibilities.
+
+    The remaining args and kwargs are the same as for
+    :meth:`~matplotlib.axes.Axes.plot`.
+
+    **Example:**
+
+        .. plot:: mpl_examples/pylab_examples/triplot_demo.py
+    """
+    tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args, **kwargs)
+
+    x = tri.x
+    y = tri.y
+    edges = tri.edges
+
+    # If draw both lines and markers at the same time, e.g.
+    #     ax.plot(x[edges].T, y[edges].T, *args, **kwargs)
+    # then the markers are drawn more than once which is incorrect if alpha<1.
+    # Hence draw of lines and markers separately.
+
+    # Draw lines without markers.
+    marker = kwargs.pop('marker', None)
+    kwargs['marker'] = ''
+    ax.plot(x[edges].T, y[edges].T, *args, **kwargs)
+    if marker is None:
+        kwargs.pop('marker')
+    else:
+        kwargs['marker'] = marker
+
+    # Draw markers without lines.
+    # Should avoid drawing markers for points that are not in any triangle?
+    kwargs['linestyle'] = ''
+    ax.plot(x, y, *args, **kwargs)
Index: lib/matplotlib/tri/__init__.py
===================================================================
--- lib/matplotlib/tri/__init__.py	(revision 0)
+++ lib/matplotlib/tri/__init__.py	(revision 0)
@@ -0,0 +1,8 @@
+"""
+Unstructured triangular grid functions.
+"""
+
+from triangulation import *
+from tricontour import *
+from tripcolor import *
+from triplot import *
Index: lib/matplotlib/tri/triangulation.py
===================================================================
--- lib/matplotlib/tri/triangulation.py	(revision 0)
+++ lib/matplotlib/tri/triangulation.py	(revision 0)
@@ -0,0 +1,219 @@
+import matplotlib.delaunay as delaunay
+import matplotlib._tri as _tri
+import numpy as np
+
+class Triangulation(object):
+    """
+    An unstructured triangular grid consisting of npoints points and
+    ntri triangles.  The triangles can either be specified by the user
+    or automatically generated using a Delaunay triangulation.
+
+    Read-only attributes:
+
+      *x*: array of shape (npoints).
+        x-coordinates of grid points.
+
+      *y*: array of shape (npoints).
+        y-coordinates of grid points.
+
+      *triangles*: integer array of shape (ntri,3).
+        For each triangle, the indices of the three points that make
+        up the triangle, ordered in an anticlockwise manner.
+
+      *mask*: optional boolean array of shape (ntri).
+        Which triangles are masked out.
+
+      *edges*: integer array of shape (?,2).
+        All edges of non-masked triangles.  Each edge is the start
+        point index and end point index.  Each edge (start,end and
+        end,start) appears only once.
+
+      *neighbors*: integer array of shape (ntri,3).
+        For each triangle, the indices of the three triangles that
+        share the same edges, or -1 if there is no such neighboring
+        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].
+    """
+    def __init__(self, x, y, triangles=None, mask=None):
+        """
+        Create a Triangulation object.
+
+        The first two arguments must be:
+
+        *x*, *y*: arrays of shape (npoints).
+          Point coordinates.
+
+        Optional arguments (args or keyword args):
+
+        *triangles*: integer array of shape (ntri,3).
+          For each triangle, the indices of the three points that make
+          up the triangle.  If the points are ordered in a clockwise
+          manner, they are converted to anticlockwise.
+
+          If not specified, matplotlib.delaunay is used to create a
+          Delaunay triangulation of the points.
+
+        *mask*: optional boolean array of shape (ntri).
+          Which triangles are masked out.
+        """
+        self.x = np.asarray(x, dtype=np.float64)
+        self.y = np.asarray(y, dtype=np.float64)
+        if self.x.shape != self.y.shape or len(self.x.shape) != 1:
+            raise ValueError("x and y must be equal-length 1-D arrays")
+
+        self.mask = None
+        self._edges = None
+        self._neighbors = None
+
+        if triangles is None:
+            # No triangulation specified, so use matplotlib.delaunay.
+            dt = delaunay.Triangulation(self.x, self.y)
+            self.triangles = np.asarray(dt.triangle_nodes, dtype=np.int32)
+            if mask is None:
+                self._edges = np.asarray(dt.edge_db, dtype=np.int32)
+                # Delaunay triangle_neighbors uses different edge indexing,
+                # so convert.
+                neighbors = np.asarray(dt.triangle_neighbors, dtype=np.int32)
+                self._neighbors = np.roll(neighbors, 1, axis=1)
+        else:
+            # Triangulation specified.
+            self.triangles = np.asarray(triangles, dtype=np.int32)
+            if self.triangles.ndim != 2 or self.triangles.shape[1] != 3:
+                raise ValueError('triangles must be a (?,3) array')
+            if self.triangles.max() >= len(self.x):
+                raise ValueError('triangles max element is out of bounds')
+            if self.triangles.min() < 0:
+                raise ValueError('triangles min element is out of bounds')
+
+        if mask is not None:
+            self.mask = np.asarray(mask, dtype=np.bool)
+            if len(self.mask.shape) != 1 or \
+                    self.mask.shape[0] != self.triangles.shape[0]:
+                raise ValueError('mask array must have same length as '
+                                 'triangles array')
+
+        # Underlying C++ object is not created until first needed.
+        self._cpp_triangulation = None
+
+    @property
+    def edges(self):
+        if self._edges is None:
+            self._edges = self.get_cpp_triangulation().get_edges()
+        return self._edges
+
+    def get_cpp_triangulation(self):
+        """
+        Return the underlying C++ Triangulation object, creating it
+        if necessary.
+        """
+        if self._cpp_triangulation is None:
+            self._cpp_triangulation = _tri.Triangulation(
+                self.x, self.y, self.triangles, self.mask, self._edges,
+                self._neighbors)
+        return self._cpp_triangulation
+
+    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)
+        else:
+            return self.triangles
+
+    @staticmethod
+    def get_from_args_and_kwargs(*args, **kwargs):
+        """
+        Return a Triangulation object from the args and kwargs, and
+        the remaining args and kwargs with the consumed values removed.
+
+        There are two alternatives: either the first argument is a
+        Triangulation object, in which case it is returned, or the args
+        and kwargs are sufficient to create a new Triangulation to
+        return.  In the latter case, see Triangulation.__init__ for
+        the possible args and kwargs.
+        """
+        if isinstance(args[0], Triangulation):
+            triangulation = args[0]
+            args = args[1:]
+        else:
+            x = args[0]
+            y = args[1]
+            args = args[2:]  # Consumed first two args.
+            ignore_remaining_args = True
+
+            # Check triangles in kwargs then args.
+            triangles = kwargs.pop('triangles', None)
+            from_args = False
+            if triangles is None and len(args) > 0:
+                triangles = args[0]
+                from_args = True
+
+            if triangles is not None:
+                try:
+                    triangles = np.asarray(triangles, dtype=np.int32)
+                except ValueError:
+                    triangles = None
+
+            if triangles is not None and (triangles.ndim != 2 or
+                                          triangles.shape[1] != 3):
+                triangles = None
+
+            if triangles is not None and from_args:
+                    args = args[1:]  # Consumed first item in args.
+                    ignore_remaining_args = False
+
+            # Check for mask in kwargs then args.
+            mask = kwargs.pop('mask', None)
+            from_args = False
+            if mask is None and not ignore_remaining_args and len(args) > 0:
+                mask = args[0]
+                from_args = True
+
+            if mask is not None:
+                try:
+                    mask = np.asarray(mask, dtype=np.bool)
+                except ValueError:
+                    mask = None
+
+            if mask is not None and mask.ndim != 1:
+                mask = None
+
+            if mask is not None and triangles is not None and \
+                    len(mask) != triangles.shape[0]:
+                mask = None
+
+            if mask is not None and from_args:
+                args = args[1:]  # Consumed first item in args.
+
+            triangulation = Triangulation(x, y, triangles, mask)
+        return triangulation, args, kwargs
+
+    @property
+    def neighbors(self):
+        if self._neighbors is None:
+            self._neighbors = self._get_cpp_triangulation().get_neighbors()
+        return self._neighbors
+
+    def set_mask(self, mask):
+        """
+        Set or clear the mask array.  This is either None, or a boolean
+        array of shape (ntri).
+        """
+        if mask is None:
+            self.mask = None
+        else:
+            self.mask = np.asarray(mask, dtype=np.bool)
+            if len(self.mask.shape) != 1 or \
+                    self.mask.shape[0] != self.triangles.shape[0]:
+                raise ValueError('mask array must have same length as '
+                                 'triangles array')
+
+        # Set mask in C++ Triangulation.
+        if self._cpp_triangulation is not None:
+            self._cpp_triangulation.set_mask(self.mask)
+
+        # Clear derived fields so they are recalculated when needed.
+        self._edges = None
+        self._neighbors = None
Index: lib/matplotlib/tri/_tri.cpp
===================================================================
--- lib/matplotlib/tri/_tri.cpp	(revision 0)
+++ lib/matplotlib/tri/_tri.cpp	(revision 0)
@@ -0,0 +1,1112 @@
+#include "_tri.h"
+#include "src/mplutils.h"
+
+#include <algorithm>
+#include <iostream>
+#include <set>
+
+#define MOVETO 1
+#define LINETO 2
+
+
+
+TriEdge::TriEdge()
+    : tri(-1), edge(-1)
+{}
+
+TriEdge::TriEdge(int tri_, int edge_)
+    : tri(tri_), edge(edge_)
+{}
+
+bool TriEdge::operator<(const TriEdge& other) const
+{
+    if (tri != other.tri)
+        return tri < other.tri;
+    else
+        return edge < other.edge;
+}
+
+bool TriEdge::operator==(const TriEdge& other) const
+{
+    return tri == other.tri && edge == other.edge;
+}
+
+bool TriEdge::operator!=(const TriEdge& other) const
+{
+    return !operator==(other);
+}
+
+std::ostream& operator<<(std::ostream& os, const TriEdge& tri_edge)
+{
+    return os << tri_edge.tri << ' ' << tri_edge.edge;
+}
+
+
+
+XY::XY()
+{}
+
+XY::XY(const double& x_, const double& y_)
+    : x(x_), y(y_)
+{}
+
+double XY::angle() const
+{
+    return atan2(y, x);
+}
+
+double XY::cross_z(const XY& other) const
+{
+    return x*other.y - y*other.x;
+}
+
+bool XY::operator==(const XY& other) const
+{
+    return x == other.x && y == other.y;
+}
+
+bool XY::operator!=(const XY& other) const
+{
+    return x != other.x || y != other.y;
+}
+
+bool XY::operator<(const XY& other) const
+{
+    if (y == other.y)
+        return x < other.x;
+    else
+        return y < other.y;
+}
+
+XY XY::operator*(const double& multiplier) const
+{
+    return XY(x*multiplier, y*multiplier);
+}
+
+XY XY::operator+(const XY& other) const
+{
+    return XY(x + other.x, y + other.y);
+}
+
+XY XY::operator-(const XY& other) const
+{
+    return XY(x - other.x, y - other.y);
+}
+
+std::ostream& operator<<(std::ostream& os, const XY& xy)
+{
+    return os << '(' << xy.x << ' ' << xy.y << ')';
+}
+
+
+
+BoundingBox::BoundingBox()
+    : empty(true)
+{}
+
+void BoundingBox::add(const XY& point)
+{
+    if (empty) {
+        empty = false;
+        lower = upper = point;
+    } else {
+        if      (point.x < lower.x) lower.x = point.x;
+        else if (point.x > upper.x) upper.x = point.x;
+
+        if      (point.y < lower.y) lower.y = point.y;
+        else if (point.y > upper.y) upper.y = point.y;
+    }
+}
+
+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)
+{
+    if (box.empty)
+        return os << "<empty>";
+    else
+        return os << box.lower << " -> " << box.upper;
+}
+
+
+
+ContourLine::ContourLine()
+    : std::vector<XY>()
+{}
+
+void ContourLine::insert_unique(iterator pos, const XY& point)
+{
+    if (empty() || pos == end() || point != *pos)
+        std::vector<XY>::insert(pos, point);
+}
+
+void ContourLine::push_back(const XY& point)
+{
+    if (empty() || point != back())
+        std::vector<XY>::push_back(point);
+}
+
+void ContourLine::write() const
+{
+    std::cout << "ContourLine of " << size() << " points:";
+    for (const_iterator it = begin(); it != end(); ++it)
+        std::cout << ' ' << *it;
+    std::cout << std::endl;
+}
+
+
+
+void write_contour(const Contour& contour)
+{
+    std::cout << "Contour of " << contour.size() << " lines." << std::endl;
+    for (Contour::const_iterator it = contour.begin(); it != contour.end(); ++it)
+        it->write();
+}
+
+
+
+
+Triangulation::Triangulation(PyArrayObject* x,
+                             PyArrayObject* y,
+                             PyArrayObject* triangles,
+                             PyArrayObject* mask,
+                             PyArrayObject* edges,
+                             PyArrayObject* neighbors)
+    : _npoints(PyArray_DIM(x,0)),
+      _ntri(PyArray_DIM(triangles,0)),
+      _x(x),
+      _y(y),
+      _triangles(triangles),
+      _mask(mask),
+      _edges(edges),
+      _neighbors(neighbors)
+{
+    _VERBOSE("Triangulation::Triangulation");
+    correct_triangles();
+}
+
+Triangulation::~Triangulation()
+{
+    _VERBOSE("Triangulation::~Triangulation");
+    Py_XDECREF(_x);
+    Py_XDECREF(_y);
+    Py_XDECREF(_triangles);
+    Py_XDECREF(_mask);
+    Py_XDECREF(_edges);
+    Py_XDECREF(_neighbors);
+}
+
+void Triangulation::calculate_boundaries()
+{
+    _VERBOSE("Triangulation::calculate_boundaries");
+
+    get_neighbors();  // Ensure _neighbors has been created.
+
+    // Create set of all boundary TriEdges, which are those which do not
+    // have a neighbor triangle.
+    typedef std::set<TriEdge> BoundaryEdges;
+    BoundaryEdges boundary_edges;
+    for (int tri = 0; tri < _ntri; ++tri) {
+        if (!is_masked(tri)) {
+            for (int edge = 0; edge < 3; ++edge) {
+                if (get_neighbor(tri, edge) == -1) {
+                    boundary_edges.insert(TriEdge(tri, edge));
+                }
+            }
+        }
+    }
+
+    // Take any boundary edge and follow the boundary until return to start
+    // point, removing edges from boundary_edges as they are used.  At the same
+    // time, initialise the _tri_edge_to_boundary_map.
+    while (!boundary_edges.empty()) {
+        // Start of new boundary.
+        BoundaryEdges::iterator it = boundary_edges.begin();
+        int tri = it->tri;
+        int edge = it->edge;
+        _boundaries.push_back(Boundary());
+        Boundary& boundary = _boundaries.back();
+
+        while (true) {
+            boundary.push_back(TriEdge(tri, edge));
+            boundary_edges.erase(it);
+            _tri_edge_to_boundary_map[TriEdge(tri, edge)] =
+                BoundaryEdge(_boundaries.size()-1, boundary.size()-1);
+
+            // Move to next edge of current triangle.
+            edge = (edge+1) % 3;
+
+            // Find start point index of boundary edge.
+            int point = get_triangle_point(tri, edge);
+
+            // Find next TriEdge by traversing neighbors until find one
+            // without a neighbor.
+            while (get_neighbor(tri, edge) != -1) {
+                tri = get_neighbor(tri, edge);
+                edge = get_edge_in_triangle(tri, point);
+            }
+
+            if (TriEdge(tri,edge) == boundary.front())
+                break;  // Reached beginning of this boundary, so finished it.
+            else
+                it = boundary_edges.find(TriEdge(tri, edge));
+        }
+    }
+}
+
+void Triangulation::calculate_edges()
+{
+    _VERBOSE("Triangulation::calculate_edges");
+    Py_XDECREF(_edges);
+
+    // Create set of all edges, storing them with start point index less than
+    // end point index.
+    typedef std::set<Edge> EdgeSet;
+    EdgeSet edge_set;
+    for (int tri = 0; tri < _ntri; ++tri) {
+        if (!is_masked(tri)) {
+            for (int edge = 0; edge < 3; edge++) {
+                int start = get_triangle_point(tri, edge);
+                int end   = get_triangle_point(tri, (edge+1)%3);
+                edge_set.insert(start > end ? Edge(start,end) : Edge(end,start));
+            }
+        }
+    }
+
+    // Convert to python _edges array.
+    npy_intp dims[2] = {edge_set.size(), 2};
+    _edges = (PyArrayObject*)PyArray_SimpleNew(2, dims, PyArray_INT);
+    int* edges_ptr = (int*)PyArray_DATA(_edges);
+    for (EdgeSet::const_iterator it = edge_set.begin(); it != edge_set.end(); ++it) {
+        *edges_ptr++ = it->start;
+        *edges_ptr++ = it->end;
+    }
+}
+
+void Triangulation::calculate_neighbors()
+{
+    _VERBOSE("Triangulation::calculate_neighbors");
+    Py_XDECREF(_neighbors);
+
+    // Create _neighbors array with shape (ntri,3) and initialise all to -1.
+    npy_intp dims[2] = {_ntri,3};
+    _neighbors = (PyArrayObject*)PyArray_SimpleNew(2, dims, PyArray_INT);
+    int* neighbors_ptr = (int*)PyArray_DATA(_neighbors);
+    std::fill(neighbors_ptr, neighbors_ptr + 3*_ntri, -1);
+
+    // For each triangle edge (start to end point), find corresponding neighbor
+    // edge from end to start point.  Do this by traversing all edges and
+    // storing them in a map from edge to TriEdge.  If corresponding neighbor
+    // edge is already in the map, don't need to store new edge as neighbor
+    // already found.
+    typedef std::map<Edge, TriEdge> EdgeToTriEdgeMap;
+    EdgeToTriEdgeMap edge_to_tri_edge_map;
+    for (int tri = 0; tri < _ntri; ++tri) {
+        if (!is_masked(tri)) {
+            for (int edge = 0; edge < 3; ++edge) {
+                int start = get_triangle_point(tri, edge);
+                int end   = get_triangle_point(tri, (edge+1)%3);
+                EdgeToTriEdgeMap::iterator it =
+                    edge_to_tri_edge_map.find(Edge(end,start));
+                if (it == edge_to_tri_edge_map.end()) {
+                    // No neighbor edge exists in the edge_to_tri_edge_map, so
+                    // add this edge to it.
+                    edge_to_tri_edge_map[Edge(start,end)] = TriEdge(tri,edge);
+                } else {
+                    // Neighbor edge found, set the two elements of _neighbors
+                    // and remove edge from edge_to_tri_edge_map.
+                    neighbors_ptr[3*tri + edge] = it->second.tri;
+                    neighbors_ptr[3*it->second.tri + it->second.edge] = tri;
+                    edge_to_tri_edge_map.erase(it);
+                }
+            }
+        }
+    }
+
+    // Note that remaining edges in the edge_to_tri_edge_map correspond to
+    // boundary edges, but the boundaries are calculated separately elsewhere.
+}
+
+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++);
+        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));
+        }
+    }
+}
+
+const Triangulation::Boundaries& Triangulation::get_boundaries() const
+{
+    _VERBOSE("Triangulation::get_boundaries");
+    if (_boundaries.empty())
+        const_cast<Triangulation*>(this)->calculate_boundaries();
+    return _boundaries;
+}
+
+void Triangulation::get_boundary_edge(const TriEdge& triEdge,
+                                      int& boundary,
+                                      int& edge) const
+{
+    get_boundaries();  // Ensure _tri_edge_to_boundary_map has been created.
+    TriEdgeToBoundaryMap::const_iterator it =
+        _tri_edge_to_boundary_map.find(triEdge);
+    assert(it != _tri_edge_to_boundary_map.end() &&
+           "TriEdge is not on a boundary");
+    boundary = it->second.boundary;
+    edge = it->second.edge;
+}
+
+int Triangulation::get_edge_in_triangle(int tri, int point) const
+{
+    assert(tri >= 0 && tri < _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;
+    for (int edge = 0; edge < 3; ++edge) {
+        if (*triangles_ptr++ == point) return edge;
+    }
+    return -1;  // point is not in triangle.
+}
+
+Py::Object Triangulation::get_edges()
+{
+    _VERBOSE("Triangulation::get_edges");
+    if (_edges == 0)
+        calculate_edges();
+    return Py::asObject(Py::new_reference_to((PyObject*)_edges));
+}
+
+int Triangulation::get_neighbor(int tri, int edge) const
+{
+    assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds");
+    assert(edge >= 0 && edge < 3 && "Edge index out of bounds");
+    return get_neighbors_ptr()[3*tri + edge];
+}
+
+TriEdge Triangulation::get_neighbor_edge(int tri, int edge) const
+{
+    int neighbor_tri = get_neighbor(tri, edge);
+    if (neighbor_tri == -1)
+        return TriEdge(-1,-1);
+    else
+        return TriEdge(neighbor_tri,
+                       get_edge_in_triangle(neighbor_tri,
+                                            get_triangle_point(tri,
+                                                               (edge+1)%3)));
+}
+
+Py::Object Triangulation::get_neighbors()
+{
+    _VERBOSE("Triangulation::get_neighbors");
+    if (_neighbors == 0) calculate_neighbors();
+    return Py::asObject(Py::new_reference_to((PyObject*)_neighbors));
+}
+
+const int* Triangulation::get_neighbors_ptr() const
+{
+    if (_neighbors == 0)
+        const_cast<Triangulation*>(this)->calculate_neighbors();
+    return (const int*)PyArray_DATA(_neighbors);
+}
+
+int Triangulation::get_npoints() const
+{
+    return _npoints;
+}
+
+int Triangulation::get_ntri() const
+{
+    return _ntri;
+}
+
+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] );
+}
+
+int Triangulation::get_triangle_point(int tri, int edge) const
+{
+    assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds");
+    assert(edge >= 0 && edge < 3 && "Edge index out of bounds");
+    return get_triangles_ptr()[3*tri + edge];
+}
+
+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);
+}
+
+void Triangulation::init_type()
+{
+    _VERBOSE("Triangulation::init_type");
+
+    behaviors().name("Triangulation");
+    behaviors().doc("Triangulation");
+
+    add_noargs_method("get_edges", &Triangulation::get_edges,
+                      "get_edges()");
+    add_noargs_method("get_neighbors", &Triangulation::get_neighbors,
+                      "get_neighbors()");
+    add_varargs_method("set_mask", &Triangulation::set_mask,
+                       "set_mask(mask)");
+}
+
+bool Triangulation::is_masked(int tri) const
+{
+    assert(tri >= 0 && tri < _ntri && "Triangle index out of bounds.");
+    return _mask && *((bool*)PyArray_DATA(_mask) + tri);
+}
+
+Py::Object Triangulation::set_mask(const Py::Tuple &args)
+{
+    _VERBOSE("Triangulation::set_mask");
+    args.verify_length(1);
+
+    Py_XDECREF(_mask);
+    _mask = 0;
+    if (args[0] != Py::None())
+    {
+        _mask = (PyArrayObject*)PyArray_ContiguousFromObject(
+                    args[0].ptr(), PyArray_BOOL, 1, 1);
+        if (_mask == 0 || PyArray_DIM(_mask,0) != PyArray_DIM(_triangles,0)) {
+            Py_XDECREF(_mask);
+            throw Py::ValueError(
+                "mask must be a 1D array with the same length as the triangles array");
+        }
+    }
+
+    // Clear derived fields so they are recalculated when needed.
+     Py_XDECREF(_edges);
+    _edges = 0;
+    Py_XDECREF(_neighbors);
+    _neighbors = 0;
+    _boundaries.clear();
+
+    return Py::None();
+}
+
+void Triangulation::write_boundaries() const
+{
+    const Boundaries& bs = get_boundaries();
+    std::cout << "Number of boundaries: " << bs.size() << std::endl;
+    for (Boundaries::const_iterator it = bs.begin(); it != bs.end(); ++it) {
+        const Boundary& b = *it;
+        std::cout << "  Boundary of " << b.size() << " points: ";
+        for (Boundary::const_iterator itb = b.begin(); itb != b.end(); ++itb) {
+            std::cout << *itb << ", ";
+        }
+        std::cout << std::endl;
+    }
+}
+
+
+
+
+TriContourGenerator::TriContourGenerator(Py::Object triangulation,
+                                         PyArrayObject* z)
+    : _triangulation(triangulation),
+      _z(z),
+      _interior_visited(2*get_triangulation().get_ntri()),
+      _boundaries_visited(0),
+      _boundaries_used(0)
+{
+    _VERBOSE("TriContourGenerator::TriContourGenerator");
+}
+
+TriContourGenerator::~TriContourGenerator()
+{
+    _VERBOSE("TriContourGenerator::~TriContourGenerator");
+    Py_XDECREF(_z);
+}
+
+void TriContourGenerator::clear_visited_flags(bool include_boundaries)
+{
+    // Clear _interiorVisited.
+    std::fill(_interior_visited.begin(), _interior_visited.end(), false);
+
+    if (include_boundaries) {
+        if (_boundaries_visited.empty()) {
+            const Boundaries& boundaries = get_boundaries();
+
+            // Initialise _boundaries_visited.
+            _boundaries_visited.reserve(boundaries.size());
+            for (Boundaries::const_iterator it = boundaries.begin();
+                    it != boundaries.end(); ++it)
+                _boundaries_visited.push_back(BoundaryVisited(it->size()));
+
+            // Initialise _boundaries_used.
+            _boundaries_used = BoundariesUsed(boundaries.size());
+        }
+
+        // Clear _boundaries_visited.
+        for (BoundariesVisited::iterator it = _boundaries_visited.begin();
+                it != _boundaries_visited.end(); ++it)
+            fill(it->begin(), it->end(), false);
+
+        // Clear _boundaries_used.
+        fill(_boundaries_used.begin(), _boundaries_used.end(), false);
+    }
+}
+
+Py::Object TriContourGenerator::contour_to_segs(const Contour& contour)
+{
+    Py::List segs(contour.size());
+    for (Contour::size_type i = 0; i < contour.size(); ++i) {
+        const ContourLine& line = contour[i];
+        npy_intp dims[2] = {line.size(),2};
+        PyArrayObject* py_line = (PyArrayObject*)PyArray_SimpleNew(
+                                                     2, dims, PyArray_DOUBLE);
+        double* p = (double*)PyArray_DATA(py_line);
+        for (ContourLine::const_iterator it = line.begin(); it != line.end(); ++it) {
+            *p++ = it->x;
+            *p++ = it->y;
+        }
+        segs[i] = Py::asObject((PyObject*)py_line);
+    }
+    return segs;
+}
+
+Py::Object TriContourGenerator::contour_to_segs_and_kinds(const Contour& contour)
+{
+    Contour::const_iterator line;
+    ContourLine::const_iterator point;
+
+    // Find total number of points in all contour lines.
+    int n_points = 0;
+    for (line = contour.begin(); line != contour.end(); ++line)
+        n_points += line->size();
+
+    // Create segs array for point coordinates.
+    npy_intp segs_dims[2] = {n_points, 2};
+    PyArrayObject* segs = (PyArrayObject*)PyArray_SimpleNew(
+                                                2, segs_dims, PyArray_DOUBLE);
+    double* segs_ptr = (double*)PyArray_DATA(segs);
+
+    // Create kinds array for code types.
+    npy_intp kinds_dims[1] = {n_points};
+    PyArrayObject* kinds = (PyArrayObject*)PyArray_SimpleNew(
+                                                1, kinds_dims, PyArray_UBYTE);
+    unsigned char* kinds_ptr = (unsigned char*)PyArray_DATA(kinds);
+
+    for (line = contour.begin(); line != contour.end(); ++line) {
+        for (point = line->begin(); point != line->end(); point++) {
+            *segs_ptr++ = point->x;
+            *segs_ptr++ = point->y;
+            *kinds_ptr++ = (point == line->begin() ? MOVETO : LINETO);
+        }
+    }
+
+    Py::Tuple result(2);
+    result[0] = Py::asObject((PyObject*)segs);
+    result[1] = Py::asObject((PyObject*)kinds);
+    return result;
+}
+
+Py::Object TriContourGenerator::create_contour(const Py::Tuple &args)
+{
+    _VERBOSE("TriContourGenerator::create_contour");
+    args.verify_length(1);
+
+    double level = (Py::Float)args[0];
+
+    clear_visited_flags(false);
+    Contour contour;
+
+    find_boundary_lines(contour, level);
+    find_interior_lines(contour, level, false, false);
+
+    return contour_to_segs(contour);
+}
+
+Py::Object TriContourGenerator::create_filled_contour(const Py::Tuple &args)
+{
+    _VERBOSE("TriContourGenerator::create_filled_contour");
+    args.verify_length(2);
+
+    double lower_level = (Py::Float)args[0];
+    double upper_level = (Py::Float)args[1];
+
+    clear_visited_flags(true);
+    Contour contour;
+
+    find_boundary_lines_filled(contour, lower_level, upper_level);
+    find_interior_lines(contour, lower_level, false, true);
+    find_interior_lines(contour, upper_level, true,  true);
+
+    return contour_to_segs_and_kinds(contour);
+}
+
+XY TriContourGenerator::edge_interp(int tri, int edge, const double& level)
+{
+    return interp(get_triangulation().get_triangle_point(tri, edge),
+                  get_triangulation().get_triangle_point(tri, (edge+1)%3),
+                  level);
+}
+
+void TriContourGenerator::find_boundary_lines(Contour& contour,
+                                              const double& level)
+{
+    // Traverse boundaries to find starting points for all contour lines that
+    // intersect the boundaries.  For each starting point found, follow the
+    // line to its end before continuing.
+    const Triangulation& triang = get_triangulation();
+    const Boundaries& boundaries = get_boundaries();
+    for (Boundaries::const_iterator it = boundaries.begin();
+            it != boundaries.end(); ++it) {
+        const Boundary& boundary = *it;
+        bool startAbove, endAbove = false;
+        for (Boundary::const_iterator itb = boundary.begin();
+                itb != boundary.end(); ++itb) {
+            if (itb == boundary.begin())
+                startAbove = get_z(triang.get_triangle_point(*itb)) >= level;
+            else
+                startAbove = endAbove;
+            endAbove = get_z(triang.get_triangle_point(itb->tri,
+                                                       (itb->edge+1)%3)) >= level;
+            if (startAbove && !endAbove) {
+                // This boundary edge is the start point for a contour line,
+                // so follow the line.
+                contour.push_back(ContourLine());
+                ContourLine& contour_line = contour.back();
+                TriEdge tri_edge = *itb;
+                follow_interior(contour_line, tri_edge, true, level, false);
+            }
+        }
+    }
+}
+
+void TriContourGenerator::find_boundary_lines_filled(Contour& contour,
+                                                     const double& lower_level,
+                                                     const double& upper_level)
+{
+    // Traverse boundaries to find starting points for all contour lines that
+    // intersect the boundaries.  For each starting point found, follow the
+    // line to its end before continuing.
+    const Triangulation& triang = get_triangulation();
+    const Boundaries& boundaries = get_boundaries();
+    for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
+        const Boundary& boundary = boundaries[i];
+        for (Boundary::size_type j = 0; j < boundary.size(); ++j) {
+            if (!_boundaries_visited[i][j]) {
+                // z values of start and end of this boundary edge.
+                double z_start = get_z(triang.get_triangle_point(boundary[j]));
+                double z_end = get_z(triang.get_triangle_point(
+                                   boundary[j].tri, (boundary[j].edge+1)%3));
+
+                // Does this boundary edge's z increase through upper level
+                // and/or decrease through lower level?
+                bool incr_upper = (z_start < upper_level && z_end >= upper_level);
+                bool decr_lower = (z_start >= lower_level && z_end < lower_level);
+
+                if (decr_lower || incr_upper) {
+                    // Start point for contour line, so follow it.
+                    contour.push_back(ContourLine());
+                    ContourLine& contour_line = contour.back();
+                    TriEdge start_tri_edge = boundary[j];
+                    TriEdge tri_edge = start_tri_edge;
+
+                    // Traverse interior and boundaries until return to start.
+                    bool on_upper = incr_upper;
+                    do {
+                        follow_interior(contour_line, tri_edge, true,
+                            on_upper ? upper_level : lower_level, on_upper);
+                        on_upper = follow_boundary(contour_line, tri_edge,
+                                       lower_level, upper_level, on_upper);
+                    } while (tri_edge != start_tri_edge);
+
+                    // Filled contour lines must not have same first and last
+                    // points.
+                    if (contour_line.size() > 1 &&
+                            contour_line.front() == contour_line.back())
+                        contour_line.pop_back();
+                }
+            }
+        }
+    }
+
+    // Add full boundaries that lie between the lower and upper levels.  These
+    // are boundaries that have not been touched by an internal contour line
+    // which are stored in _boundaries_used.
+    for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {
+        if (!_boundaries_used[i]) {
+            const Boundary& boundary = boundaries[i];
+            double z = get_z(triang.get_triangle_point(boundary[0]));
+            if (z >= lower_level && z < upper_level) {
+                contour.push_back(ContourLine());
+                ContourLine& contour_line = contour.back();
+                for (Boundary::size_type j = 0; j < boundary.size(); ++j)
+                    contour_line.push_back(triang.get_point_coords(
+                                      triang.get_triangle_point(boundary[j])));
+            }
+        }
+    }
+}
+
+void TriContourGenerator::find_interior_lines(Contour& contour,
+                                              const double& level,
+                                              bool on_upper,
+                                              bool filled)
+{
+    const Triangulation& triang = get_triangulation();
+    int ntri = triang.get_ntri();
+    for (int tri = 0; tri < ntri; ++tri) {
+        int visited_index = (on_upper ? tri+ntri : tri);
+
+        if (_interior_visited[visited_index] || triang.is_masked(tri))
+            continue;  // Triangle has already been visited or is masked.
+
+        _interior_visited[visited_index] = true;
+
+        // Determine edge via which to leave this triangle.
+        int edge = get_exit_edge(tri, level, on_upper);
+        assert(edge >= -1 && edge < 3 && "Invalid exit edge");
+        if (edge == -1)
+            continue;  // Contour does not pass through this triangle.
+
+        // Found start of new contour line loop.
+        contour.push_back(ContourLine());
+        ContourLine& contour_line = contour.back();
+        TriEdge tri_edge = triang.get_neighbor_edge(tri, edge);
+        follow_interior(contour_line, tri_edge, false, level, on_upper);
+
+        if (!filled)
+            // Non-filled contour lines must be closed.
+            contour_line.push_back(contour_line.front());
+        else if (contour_line.size() > 1 &&
+                 contour_line.front() == contour_line.back())
+            // Filled contour lines must not have same first and last points.
+            contour_line.pop_back();
+    }
+}
+
+bool TriContourGenerator::follow_boundary(ContourLine& contour_line,
+                                          TriEdge& tri_edge,
+                                          const double& lower_level,
+                                          const double& upper_level,
+                                          bool on_upper)
+{
+    const Triangulation& triang = get_triangulation();
+    const Boundaries& boundaries = get_boundaries();
+
+    // Have TriEdge to start at, need equivalent boundary edge.
+    int boundary, edge;
+    triang.get_boundary_edge(tri_edge, boundary, edge);
+    _boundaries_used[boundary] = true;
+
+    bool stop = false;
+    bool first_edge = true;
+    double z_start, z_end = 0;
+    while (!stop)
+    {
+        assert(!_boundaries_visited[boundary][edge] && "Boundary already visited");
+        _boundaries_visited[boundary][edge] = true;
+
+        // z values of start and end points of boundary edge.
+        if (first_edge)
+            z_start = get_z(triang.get_triangle_point(tri_edge));
+        else
+            z_start = z_end;
+        z_end = get_z(triang.get_triangle_point(tri_edge.tri,
+                                                (tri_edge.edge+1)%3));
+
+        if (z_end > z_start) {  // z increasing.
+            if (!(!on_upper && first_edge) &&
+                z_end >= lower_level && z_start < lower_level) {
+                stop = true;
+                on_upper = false;
+            } else if (z_end >= upper_level && z_start < upper_level) {
+                stop = true;
+                on_upper = true;
+            }
+        } else {  // z decreasing.
+            if (!(on_upper && first_edge) &&
+                z_start >= upper_level && z_end < upper_level) {
+                stop = true;
+                on_upper = true;
+            } else if (z_start >= lower_level && z_end < lower_level) {
+                stop = true;
+                on_upper = false;
+            }
+        }
+
+        first_edge = false;
+
+        if (!stop) {
+            // Move to next boundary edge, adding point to contour line.
+            edge = (edge+1) % (int)boundaries[boundary].size();
+            tri_edge = boundaries[boundary][edge];
+            contour_line.push_back(triang.get_point_coords(
+                                       triang.get_triangle_point(tri_edge)));
+        }
+    }
+
+    return on_upper;
+}
+
+void TriContourGenerator::follow_interior(ContourLine& contour_line,
+                                          TriEdge& tri_edge,
+                                          bool end_on_boundary,
+                                          const double& level,
+                                          bool on_upper)
+{
+    int& tri = tri_edge.tri;
+    int& edge = tri_edge.edge;
+
+    // Initial point.
+    contour_line.push_back(edge_interp(tri, edge, level));
+
+    while (true) {
+        int visited_index = tri;
+        if (on_upper)
+            visited_index += get_triangulation().get_ntri();
+
+        // Check for end not on boundary.
+        if (!end_on_boundary && _interior_visited[visited_index])
+            break;  // Reached start point, so return.
+
+        // Determine edge by which to leave this triangle.
+        edge = get_exit_edge(tri, level, on_upper);
+        assert(edge >= 0 && edge < 3 && "Invalid exit edge");
+
+        _interior_visited[visited_index] = true;
+
+        // Append new point to point set.
+        assert(edge >= 0 && edge < 3 && "Invalid triangle edge");
+        contour_line.push_back(edge_interp(tri, edge, level));
+
+        // Move to next triangle.
+        TriEdge next_tri_edge = get_triangulation().get_neighbor_edge(tri,edge);
+
+        // Check if ending on a boundary.
+        if (end_on_boundary && next_tri_edge.tri == -1)
+            break;
+
+        tri_edge = next_tri_edge;
+        assert(tri_edge.tri != -1 && "Invalid triangle for internal loop");
+    }
+}
+
+const TriContourGenerator::Boundaries& TriContourGenerator::get_boundaries() const
+{
+    return get_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() &&
+           "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;
+
+    if (on_upper) config = 7-config;
+
+    switch (config) {
+        case 0: return -1;
+        case 1: return  2;
+        case 2: return  0;
+        case 3: return  2;
+        case 4: return  1;
+        case 5: return  1;
+        case 6: return  0;
+        case 7: return -1;
+        default: assert(0 && "Invalid config value"); return -1;
+    }
+}
+
+const Triangulation& TriContourGenerator::get_triangulation() const
+{
+    return *(Triangulation*)_triangulation.ptr();
+}
+
+const double& TriContourGenerator::get_z(int point) const
+{
+    assert(point >= 0 && point < get_triangulation().get_npoints() &&
+           "Point index out of bounds.");
+    return ((const double*)PyArray_DATA(_z))[point];
+}
+
+void TriContourGenerator::init_type()
+{
+    _VERBOSE("TriContourGenerator::init_type");
+
+    behaviors().name("TriContourGenerator");
+    behaviors().doc("TriContourGenerator");
+
+    add_varargs_method("create_contour",
+                       &TriContourGenerator::create_contour,
+                       "create_contour(level)");
+    add_varargs_method("create_filled_contour",
+                       &TriContourGenerator::create_filled_contour,
+                       "create_filled_contour(lower_level, upper_level)");
+}
+
+XY TriContourGenerator::interp(int point1,
+                               int point2,
+                               const double& level) const
+{
+    assert(point1 >= 0 && point1 < get_triangulation().get_npoints() &&
+           "Point index 1 out of bounds.");
+    assert(point2 >= 0 && point2 < get_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);
+}
+
+
+
+
+
+#if defined(_MSC_VER)
+DL_EXPORT(void)
+#elif defined(__cplusplus)
+extern "C" void
+#else
+void
+#endif
+init_tri()
+{
+    static TriModule* triModule = new TriModule;
+    import_array();
+}
+
+TriModule::TriModule()
+    : Py::ExtensionModule<TriModule>("tri")
+{
+    Triangulation::init_type();
+    TriContourGenerator::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");
+
+    initialize("Module for unstructured triangular grids");
+}
+
+Py::Object TriModule::new_triangulation(const Py::Tuple &args)
+{
+    _VERBOSE("TriModule::new_triangulation");
+    args.verify_length(6);
+
+    // x and y.
+    PyArrayObject* x = (PyArrayObject*)PyArray_ContiguousFromObject(
+                           args[0].ptr(), PyArray_DOUBLE, 1, 1);
+    PyArrayObject* y = (PyArrayObject*)PyArray_ContiguousFromObject(
+                           args[1].ptr(), PyArray_DOUBLE, 1, 1);
+    if (x == 0 || y == 0 || PyArray_DIM(x,0) != PyArray_DIM(y,0)) {
+        Py_XDECREF(x);
+        Py_XDECREF(y);
+        throw Py::ValueError("x and y must be 1D arrays of the same length");
+    }
+
+    // triangles.
+    PyArrayObject* triangles = (PyArrayObject*)PyArray_ContiguousFromObject(
+                                   args[2].ptr(), PyArray_INT, 2, 2);
+    if (triangles == 0 || PyArray_DIM(triangles,1) != 3) {
+        Py_XDECREF(x);
+        Py_XDECREF(y);
+        Py_XDECREF(triangles);
+        throw Py::ValueError("triangles must be a 2D array of shape (?,3)");
+    }
+
+    // Optional mask.
+    PyArrayObject* mask = 0;
+    if (args[3].ptr() != 0 && args[3] != Py::None())
+    {
+        mask = (PyArrayObject*)PyArray_ContiguousFromObject(
+                   args[3].ptr(), PyArray_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);
+            throw Py::ValueError(
+                "mask must be a 1D array with the same length as the triangles array");
+        }
+    }
+
+    // Optional edges.
+    PyArrayObject* edges = 0;
+    if (args[4].ptr() != 0 && args[4] != Py::None())
+    {
+        edges = (PyArrayObject*)PyArray_ContiguousFromObject(
+                    args[4].ptr(), PyArray_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);
+            throw Py::ValueError("edges must be a 2D array with shape (?,2)");
+        }
+    }
+
+    // Optional neighbors.
+    PyArrayObject* neighbors = 0;
+    if (args[5].ptr() != 0 && args[5] != Py::None())
+    {
+        neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(
+                        args[5].ptr(), PyArray_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);
+            throw Py::ValueError(
+                "neighbors must be a 2D array with the same shape as the triangles array");
+        }
+    }
+
+    return Py::asObject(new Triangulation(x, y, triangles, mask, edges, neighbors));
+}
+
+Py::Object TriModule::new_tricontourgenerator(const Py::Tuple &args)
+{
+    _VERBOSE("TriModule::new_tricontourgenerator");
+    args.verify_length(2);
+
+    Py::Object tri = args[0];
+    if (!Triangulation::check(tri))
+        throw Py::ValueError("Expecting a C++ Triangulation object");
+
+    PyArrayObject* z = (PyArrayObject*)PyArray_ContiguousFromObject(
+                           args[1].ptr(), PyArray_DOUBLE, 1, 1);
+    if (z == 0 ||
+        PyArray_DIM(z,0) != ((Triangulation*)tri.ptr())->get_npoints()) {
+        Py_XDECREF(z);
+        throw Py::ValueError(
+            "z must be a 1D array with the same length as the x and y arrays");
+    }
+
+    return Py::asObject(new TriContourGenerator(tri, z));
+}
+
Index: lib/matplotlib/pyplot.py
===================================================================
--- lib/matplotlib/pyplot.py	(revision 8232)
+++ lib/matplotlib/pyplot.py	(working copy)
@@ -2427,6 +2427,78 @@
 
 # This function was autogenerated by boilerplate.py.  Do not edit as
 # changes will be lost
+@autogen_docstring(Axes.tricontour)
+def tricontour(*args, **kwargs):
+    ax = gca()
+    # allow callers to override the hold state by passing hold=True|False
+    washold = ax.ishold()
+    hold = kwargs.pop('hold', None)
+    if hold is not None:
+        ax.hold(hold)
+    try:
+        ret = ax.tricontour(*args, **kwargs)
+        draw_if_interactive()
+    finally:
+        ax.hold(washold)
+    if ret._A is not None: sci(ret)
+    return ret
+
+# This function was autogenerated by boilerplate.py.  Do not edit as
+# changes will be lost
+@autogen_docstring(Axes.tricontourf)
+def tricontourf(*args, **kwargs):
+    ax = gca()
+    # allow callers to override the hold state by passing hold=True|False
+    washold = ax.ishold()
+    hold = kwargs.pop('hold', None)
+    if hold is not None:
+        ax.hold(hold)
+    try:
+        ret = ax.tricontourf(*args, **kwargs)
+        draw_if_interactive()
+    finally:
+        ax.hold(washold)
+    if ret._A is not None: sci(ret)
+    return ret
+
+# This function was autogenerated by boilerplate.py.  Do not edit as
+# changes will be lost
+@autogen_docstring(Axes.tripcolor)
+def tripcolor(*args, **kwargs):
+    ax = gca()
+    # allow callers to override the hold state by passing hold=True|False
+    washold = ax.ishold()
+    hold = kwargs.pop('hold', None)
+    if hold is not None:
+        ax.hold(hold)
+    try:
+        ret = ax.tripcolor(*args, **kwargs)
+        draw_if_interactive()
+    finally:
+        ax.hold(washold)
+    sci(ret)
+    return ret
+
+# This function was autogenerated by boilerplate.py.  Do not edit as
+# changes will be lost
+@autogen_docstring(Axes.triplot)
+def triplot(*args, **kwargs):
+    ax = gca()
+    # allow callers to override the hold state by passing hold=True|False
+    washold = ax.ishold()
+    hold = kwargs.pop('hold', None)
+    if hold is not None:
+        ax.hold(hold)
+    try:
+        ret = ax.triplot(*args, **kwargs)
+        draw_if_interactive()
+    finally:
+        ax.hold(washold)
+
+    return ret
+
+# This function was autogenerated by boilerplate.py.  Do not edit as
+# changes will be lost
 @autogen_docstring(Axes.vlines)
 def vlines(x, ymin, ymax, colors='k', linestyles='solid', label='', hold=None, **kwargs):
     ax = gca()
Index: lib/matplotlib/axes.py
===================================================================
--- lib/matplotlib/axes.py	(revision 8232)
+++ lib/matplotlib/axes.py	(working copy)
@@ -32,6 +32,7 @@
 import matplotlib.text as mtext
 import matplotlib.ticker as mticker
 import matplotlib.transforms as mtransforms
+import matplotlib.tri as mtri
 
 
 iterable = cbook.iterable
@@ -7829,7 +7830,23 @@
         self.xaxis.set_minor_locator(mticker.NullLocator())
         self.yaxis.set_minor_locator(mticker.NullLocator())
 
+    def tricontour(self, *args, **kwargs):
+        return mtri.tricontour(self, *args, **kwargs)
+    tricontour.__doc__ = mtri.TriContourSet.tricontour_doc
 
+    def tricontourf(self, *args, **kwargs):
+        return mtri.tricontourf(self, *args, **kwargs)
+    tricontourf.__doc__ = mtri.TriContourSet.tricontour_doc
+
+    def tripcolor(self, *args, **kwargs):
+        return mtri.tripcolor(self, *args, **kwargs)
+    tripcolor.__doc__ = mtri.tripcolor.__doc__
+
+    def triplot(self, *args, **kwargs):
+        mtri.triplot(self, *args, **kwargs)
+    triplot.__doc__ = mtri.triplot.__doc__
+
+
 class SubplotBase:
     """
     Base class for subplots, which are :class:`Axes` instances with
Index: setup.py
===================================================================
--- setup.py	(revision 8232)
+++ setup.py	(working copy)
@@ -42,7 +42,7 @@
      check_for_qt, check_for_qt4, check_for_cairo, \
      check_provide_pytz, check_provide_dateutil,\
      check_for_dvipng, check_for_ghostscript, check_for_latex, \
-     check_for_pdftops, check_for_datetime, options, build_png
+     check_for_pdftops, check_for_datetime, options, build_png, build_tri
 #import distutils.sysconfig
 
 # jdh
@@ -68,6 +68,7 @@
     'matplotlib.numerix.linear_algebra',
     'matplotlib.numerix.random_array',
     'matplotlib.numerix.fft',
+    'matplotlib.tri',
 
     ]
 
@@ -129,6 +130,7 @@
 build_delaunay(ext_modules, packages)
 build_nxutils(ext_modules, packages)
 build_path(ext_modules, packages)
+build_tri(ext_modules, packages)
 
 print_raw("")
 print_raw("OPTIONAL BACKEND DEPENDENCIES")
Index: boilerplate.py
===================================================================
--- boilerplate.py	(revision 8232)
+++ boilerplate.py	(working copy)
@@ -89,6 +89,10 @@
     #'spy',
     'stem',
     'step',
+    'tricontour',
+    'tricontourf',
+    'tripcolor',
+    'triplot',
     'vlines',
     'xcorr',
     'barbs',
@@ -114,6 +118,9 @@
     #'spy'    :  'sci(%(ret)s)',  ### may return image or Line2D
     'quiver' :  'sci(%(ret)s)',
     'specgram'  : 'sci(%(ret)s[-1])',
+    'tricontour' : 'if %(ret)s._A is not None: sci(%(ret)s)',
+    'tricontourf': 'if %(ret)s._A is not None: sci(%(ret)s)',
+    'tripcolor'  : 'sci(%(ret)s)',
 
 }
 
Index: setupext.py
===================================================================
--- setupext.py	(revision 8232)
+++ setupext.py	(working copy)
@@ -102,6 +102,7 @@
 BUILT_CONTOUR   = False
 BUILT_GDK       = False
 BUILT_PATH      = False
+BUILT_TRI       = False
 
 AGG_VERSION = 'agg24'
 TCL_TK_CACHE = None
@@ -1414,3 +1415,17 @@
 
     BUILT_GDK = True
 
+def build_tri(ext_modules, packages):
+    global BUILT_TRI
+    if BUILT_TRI: return # only build it if you you haven't already
+
+    deps = ['lib/matplotlib/tri/_tri.cpp', 'src/mplutils.cpp']
+    deps.extend(glob.glob('CXX/*.cxx'))
+    deps.extend(glob.glob('CXX/*.c'))
+
+    module = Extension('matplotlib._tri', deps,
+                       define_macros=defines)
+    add_numpy_flags(module)
+    add_base_flags(module)
+    ext_modules.append(module)
+    BUILT_TRI = True
Index: examples/pylab_examples/tricontour_demo.py
===================================================================
--- examples/pylab_examples/tricontour_demo.py	(revision 0)
+++ examples/pylab_examples/tricontour_demo.py	(revision 0)
@@ -0,0 +1,98 @@
+"""
+Contour plots of unstructured triangular grids.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+import math
+
+# Creating a Triangulation without specifying the triangles results in the
+# Delaunay triangulation of the points.
+
+# First create the x and y coordinates of the points.
+n_angles = 48
+n_radii = 8
+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()
+z = (np.cos(radii)*np.cos(angles*3.0)).flatten()
+
+# Create the Triangulation; no triangles so Delaunay triangulation created.
+triang = tri.Triangulation(x, y)
+
+# Mask off unwanted triangles.
+xmid = x[triang.triangles].mean(axis=1)
+ymid = y[triang.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triang.set_mask(mask)
+
+# pcolor plot.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.tricontourf(triang, z)
+plt.colorbar()
+plt.tricontour(triang, z, colors='k')
+plt.title('Contour plot of Delaunay triangulation')
+
+# You can specify your own triangulation rather than perform a Delaunay
+# triangulation of the points, where each triangle is given by the indices of
+# the three points that make up the triangle, ordered in either a clockwise or
+# anticlockwise manner.
+
+xy = np.asarray([
+    [-0.101,0.872],[-0.080,0.883],[-0.069,0.888],[-0.054,0.890],[-0.045,0.897],
+    [-0.057,0.895],[-0.073,0.900],[-0.087,0.898],[-0.090,0.904],[-0.069,0.907],
+    [-0.069,0.921],[-0.080,0.919],[-0.073,0.928],[-0.052,0.930],[-0.048,0.942],
+    [-0.062,0.949],[-0.054,0.958],[-0.069,0.954],[-0.087,0.952],[-0.087,0.959],
+    [-0.080,0.966],[-0.085,0.973],[-0.087,0.965],[-0.097,0.965],[-0.097,0.975],
+    [-0.092,0.984],[-0.101,0.980],[-0.108,0.980],[-0.104,0.987],[-0.102,0.993],
+    [-0.115,1.001],[-0.099,0.996],[-0.101,1.007],[-0.090,1.010],[-0.087,1.021],
+    [-0.069,1.021],[-0.052,1.022],[-0.052,1.017],[-0.069,1.010],[-0.064,1.005],
+    [-0.048,1.005],[-0.031,1.005],[-0.031,0.996],[-0.040,0.987],[-0.045,0.980],
+    [-0.052,0.975],[-0.040,0.973],[-0.026,0.968],[-0.020,0.954],[-0.006,0.947],
+    [ 0.003,0.935],[ 0.006,0.926],[ 0.005,0.921],[ 0.022,0.923],[ 0.033,0.912],
+    [ 0.029,0.905],[ 0.017,0.900],[ 0.012,0.895],[ 0.027,0.893],[ 0.019,0.886],
+    [ 0.001,0.883],[-0.012,0.884],[-0.029,0.883],[-0.038,0.879],[-0.057,0.881],
+    [-0.062,0.876],[-0.078,0.876],[-0.087,0.872],[-0.030,0.907],[-0.007,0.905],
+    [-0.057,0.916],[-0.025,0.933],[-0.077,0.990],[-0.059,0.993] ])
+x = xy[:,0]*180/3.14159
+y = xy[:,1]*180/3.14159
+x0 = -5
+y0 = 52
+z = np.exp(-0.01*( (x-x0)*(x-x0) + (y-y0)*(y-y0) ))
+
+triangles = np.asarray([
+    [67,66, 1],[65, 2,66],[ 1,66, 2],[64, 2,65],[63, 3,64],[60,59,57],
+    [ 2,64, 3],[ 3,63, 4],[ 0,67, 1],[62, 4,63],[57,59,56],[59,58,56],
+    [61,60,69],[57,69,60],[ 4,62,68],[ 6, 5, 9],[61,68,62],[69,68,61],
+    [ 9, 5,70],[ 6, 8, 7],[ 4,70, 5],[ 8, 6, 9],[56,69,57],[69,56,52],
+    [70,10, 9],[54,53,55],[56,55,53],[68,70, 4],[52,56,53],[11,10,12],
+    [69,71,68],[68,13,70],[10,70,13],[51,50,52],[13,68,71],[52,71,69],
+    [12,10,13],[71,52,50],[71,14,13],[50,49,71],[49,48,71],[14,16,15],
+    [14,71,48],[17,19,18],[17,20,19],[48,16,14],[48,47,16],[47,46,16],
+    [16,46,45],[23,22,24],[21,24,22],[17,16,45],[20,17,45],[21,25,24],
+    [27,26,28],[20,72,21],[25,21,72],[45,72,20],[25,28,26],[44,73,45],
+    [72,45,73],[28,25,29],[29,25,31],[43,73,44],[73,43,40],[72,73,39],
+    [72,31,25],[42,40,43],[31,30,29],[39,73,40],[42,41,40],[72,33,31],
+    [32,31,33],[39,38,72],[33,72,38],[33,38,34],[37,35,38],[34,38,35],
+    [35,37,36] ])
+
+# Rather than create a Triangulation object, can simply pass x, y and triangles
+# arrays to tripcolor directly.  It would be better to use a Triangulation object
+# if the same triangulation was to be used more than once to save duplicated
+# calculations.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.tricontourf(x, y, triangles, z)
+plt.colorbar()
+plt.title('Contour plot of user-specified triangulation')
+plt.xlabel('Longitude (degrees)')
+plt.ylabel('Latitude (degrees)')
+
+plt.show()
Index: examples/pylab_examples/tricontour_vs_griddata.py
===================================================================
--- examples/pylab_examples/tricontour_vs_griddata.py	(revision 0)
+++ examples/pylab_examples/tricontour_vs_griddata.py	(revision 0)
@@ -0,0 +1,47 @@
+"""
+Comparison of griddata and tricontour for an unstructured triangular grid.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+from numpy.random import uniform, seed
+from matplotlib.mlab import griddata
+import time
+
+seed(0)
+npts = 200
+ngridx = 100
+ngridy = 200
+x = uniform(-2,2,npts)
+y = uniform(-2,2,npts)
+z = x*np.exp(-x**2-y**2)
+
+# griddata and contour.
+start = time.clock()
+plt.subplot(211)
+xi = np.linspace(-2.1,2.1,ngridx)
+yi = np.linspace(-2.1,2.1,ngridy)
+zi = griddata(x,y,z,xi,yi,interp='linear')
+plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
+plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
+plt.colorbar() # draw colorbar
+plt.plot(x, y, 'ko', ms=3)
+plt.xlim(-2,2)
+plt.ylim(-2,2)
+plt.title('griddata and contour (%d points, %d grid points)' % (npts, ngridx*ngridy))
+print 'griddata and contour seconds:', time.clock() - start
+
+# tricontour.
+start = time.clock()
+plt.subplot(212)
+triang = tri.Triangulation(x, y)
+plt.tricontour(x, y, z, 15, linewidths=0.5, colors='k')
+plt.tricontourf(x, y, z, 15, cmap=plt.cm.jet)
+plt.colorbar()
+plt.plot(x, y, 'ko', ms=3)
+plt.xlim(-2,2)
+plt.ylim(-2,2)
+plt.title('tricontour (%d points)' % npts)
+print 'tricontour seconds:', time.clock() - start
+
+plt.show()
Index: examples/pylab_examples/tripcolor_demo.py
===================================================================
--- examples/pylab_examples/tripcolor_demo.py	(revision 0)
+++ examples/pylab_examples/tripcolor_demo.py	(revision 0)
@@ -0,0 +1,98 @@
+"""
+Pseudocolor plots of unstructured triangular grids.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+import math
+
+# Creating a Triangulation without specifying the triangles results in the
+# Delaunay triangulation of the points.
+
+# First create the x and y coordinates of the points.
+n_angles = 36
+n_radii = 8
+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()
+z = (np.cos(radii)*np.cos(angles*3.0)).flatten()
+
+# Create the Triangulation; no triangles so Delaunay triangulation created.
+triang = tri.Triangulation(x, y)
+
+# Mask off unwanted triangles.
+xmid = x[triang.triangles].mean(axis=1)
+ymid = y[triang.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triang.set_mask(mask)
+
+# pcolor plot.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.tripcolor(triang, z, shading='faceted')
+plt.colorbar()
+plt.title('tripcolor of Delaunay triangulation')
+
+
+# You can specify your own triangulation rather than perform a Delaunay
+# triangulation of the points, where each triangle is given by the indices of
+# the three points that make up the triangle, ordered in either a clockwise or
+# anticlockwise manner.
+
+xy = np.asarray([
+    [-0.101,0.872],[-0.080,0.883],[-0.069,0.888],[-0.054,0.890],[-0.045,0.897],
+    [-0.057,0.895],[-0.073,0.900],[-0.087,0.898],[-0.090,0.904],[-0.069,0.907],
+    [-0.069,0.921],[-0.080,0.919],[-0.073,0.928],[-0.052,0.930],[-0.048,0.942],
+    [-0.062,0.949],[-0.054,0.958],[-0.069,0.954],[-0.087,0.952],[-0.087,0.959],
+    [-0.080,0.966],[-0.085,0.973],[-0.087,0.965],[-0.097,0.965],[-0.097,0.975],
+    [-0.092,0.984],[-0.101,0.980],[-0.108,0.980],[-0.104,0.987],[-0.102,0.993],
+    [-0.115,1.001],[-0.099,0.996],[-0.101,1.007],[-0.090,1.010],[-0.087,1.021],
+    [-0.069,1.021],[-0.052,1.022],[-0.052,1.017],[-0.069,1.010],[-0.064,1.005],
+    [-0.048,1.005],[-0.031,1.005],[-0.031,0.996],[-0.040,0.987],[-0.045,0.980],
+    [-0.052,0.975],[-0.040,0.973],[-0.026,0.968],[-0.020,0.954],[-0.006,0.947],
+    [ 0.003,0.935],[ 0.006,0.926],[ 0.005,0.921],[ 0.022,0.923],[ 0.033,0.912],
+    [ 0.029,0.905],[ 0.017,0.900],[ 0.012,0.895],[ 0.027,0.893],[ 0.019,0.886],
+    [ 0.001,0.883],[-0.012,0.884],[-0.029,0.883],[-0.038,0.879],[-0.057,0.881],
+    [-0.062,0.876],[-0.078,0.876],[-0.087,0.872],[-0.030,0.907],[-0.007,0.905],
+    [-0.057,0.916],[-0.025,0.933],[-0.077,0.990],[-0.059,0.993] ])
+x = xy[:,0]*180/3.14159
+y = xy[:,1]*180/3.14159
+x0 = -5
+y0 = 52
+z = np.exp(-0.01*( (x-x0)*(x-x0) + (y-y0)*(y-y0) ))
+
+triangles = np.asarray([
+    [67,66, 1],[65, 2,66],[ 1,66, 2],[64, 2,65],[63, 3,64],[60,59,57],
+    [ 2,64, 3],[ 3,63, 4],[ 0,67, 1],[62, 4,63],[57,59,56],[59,58,56],
+    [61,60,69],[57,69,60],[ 4,62,68],[ 6, 5, 9],[61,68,62],[69,68,61],
+    [ 9, 5,70],[ 6, 8, 7],[ 4,70, 5],[ 8, 6, 9],[56,69,57],[69,56,52],
+    [70,10, 9],[54,53,55],[56,55,53],[68,70, 4],[52,56,53],[11,10,12],
+    [69,71,68],[68,13,70],[10,70,13],[51,50,52],[13,68,71],[52,71,69],
+    [12,10,13],[71,52,50],[71,14,13],[50,49,71],[49,48,71],[14,16,15],
+    [14,71,48],[17,19,18],[17,20,19],[48,16,14],[48,47,16],[47,46,16],
+    [16,46,45],[23,22,24],[21,24,22],[17,16,45],[20,17,45],[21,25,24],
+    [27,26,28],[20,72,21],[25,21,72],[45,72,20],[25,28,26],[44,73,45],
+    [72,45,73],[28,25,29],[29,25,31],[43,73,44],[73,43,40],[72,73,39],
+    [72,31,25],[42,40,43],[31,30,29],[39,73,40],[42,41,40],[72,33,31],
+    [32,31,33],[39,38,72],[33,72,38],[33,38,34],[37,35,38],[34,38,35],
+    [35,37,36] ])
+
+# Rather than create a Triangulation object, can simply pass x, y and triangles
+# arrays to tripcolor directly.  It would be better to use a Triangulation object
+# if the same triangulation was to be used more than once to save duplicated
+# calculations.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.tripcolor(x, y, triangles, z, shading='faceted')
+plt.colorbar()
+plt.title('tripcolor of user-specified triangulation')
+plt.xlabel('Longitude (degrees)')
+plt.ylabel('Latitude (degrees)')
+
+plt.show()
Index: examples/pylab_examples/triplot_demo.py
===================================================================
--- examples/pylab_examples/triplot_demo.py	(revision 0)
+++ examples/pylab_examples/triplot_demo.py	(revision 0)
@@ -0,0 +1,92 @@
+"""
+Creating and plotting unstructured triangular grids.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+import math
+
+# Creating a Triangulation without specifying the triangles results in the
+# Delaunay triangulation of the points.
+
+# First create the x and y coordinates of the points.
+n_angles = 36
+n_radii = 8
+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()
+
+# Create the Triangulation; no triangles so Delaunay triangulation created.
+triang = tri.Triangulation(x, y)
+
+# Mask off unwanted triangles.
+xmid = x[triang.triangles].mean(axis=1)
+ymid = y[triang.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triang.set_mask(mask)
+
+# Plot the triangulation.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.triplot(triang, 'bo-')
+plt.title('triplot of Delaunay triangulation')
+
+
+# You can specify your own triangulation rather than perform a Delaunay
+# triangulation of the points, where each triangle is given by the indices of
+# the three points that make up the triangle, ordered in either a clockwise or
+# anticlockwise manner.
+
+xy = np.asarray([
+    [-0.101,0.872],[-0.080,0.883],[-0.069,0.888],[-0.054,0.890],[-0.045,0.897],
+    [-0.057,0.895],[-0.073,0.900],[-0.087,0.898],[-0.090,0.904],[-0.069,0.907],
+    [-0.069,0.921],[-0.080,0.919],[-0.073,0.928],[-0.052,0.930],[-0.048,0.942],
+    [-0.062,0.949],[-0.054,0.958],[-0.069,0.954],[-0.087,0.952],[-0.087,0.959],
+    [-0.080,0.966],[-0.085,0.973],[-0.087,0.965],[-0.097,0.965],[-0.097,0.975],
+    [-0.092,0.984],[-0.101,0.980],[-0.108,0.980],[-0.104,0.987],[-0.102,0.993],
+    [-0.115,1.001],[-0.099,0.996],[-0.101,1.007],[-0.090,1.010],[-0.087,1.021],
+    [-0.069,1.021],[-0.052,1.022],[-0.052,1.017],[-0.069,1.010],[-0.064,1.005],
+    [-0.048,1.005],[-0.031,1.005],[-0.031,0.996],[-0.040,0.987],[-0.045,0.980],
+    [-0.052,0.975],[-0.040,0.973],[-0.026,0.968],[-0.020,0.954],[-0.006,0.947],
+    [ 0.003,0.935],[ 0.006,0.926],[ 0.005,0.921],[ 0.022,0.923],[ 0.033,0.912],
+    [ 0.029,0.905],[ 0.017,0.900],[ 0.012,0.895],[ 0.027,0.893],[ 0.019,0.886],
+    [ 0.001,0.883],[-0.012,0.884],[-0.029,0.883],[-0.038,0.879],[-0.057,0.881],
+    [-0.062,0.876],[-0.078,0.876],[-0.087,0.872],[-0.030,0.907],[-0.007,0.905],
+    [-0.057,0.916],[-0.025,0.933],[-0.077,0.990],[-0.059,0.993] ])
+x = xy[:,0]*180/3.14159
+y = xy[:,1]*180/3.14159
+
+triangles = np.asarray([
+    [67,66, 1],[65, 2,66],[ 1,66, 2],[64, 2,65],[63, 3,64],[60,59,57],
+    [ 2,64, 3],[ 3,63, 4],[ 0,67, 1],[62, 4,63],[57,59,56],[59,58,56],
+    [61,60,69],[57,69,60],[ 4,62,68],[ 6, 5, 9],[61,68,62],[69,68,61],
+    [ 9, 5,70],[ 6, 8, 7],[ 4,70, 5],[ 8, 6, 9],[56,69,57],[69,56,52],
+    [70,10, 9],[54,53,55],[56,55,53],[68,70, 4],[52,56,53],[11,10,12],
+    [69,71,68],[68,13,70],[10,70,13],[51,50,52],[13,68,71],[52,71,69],
+    [12,10,13],[71,52,50],[71,14,13],[50,49,71],[49,48,71],[14,16,15],
+    [14,71,48],[17,19,18],[17,20,19],[48,16,14],[48,47,16],[47,46,16],
+    [16,46,45],[23,22,24],[21,24,22],[17,16,45],[20,17,45],[21,25,24],
+    [27,26,28],[20,72,21],[25,21,72],[45,72,20],[25,28,26],[44,73,45],
+    [72,45,73],[28,25,29],[29,25,31],[43,73,44],[73,43,40],[72,73,39],
+    [72,31,25],[42,40,43],[31,30,29],[39,73,40],[42,41,40],[72,33,31],
+    [32,31,33],[39,38,72],[33,72,38],[33,38,34],[37,35,38],[34,38,35],
+    [35,37,36] ])
+
+# Rather than create a Triangulation object, can simply pass x, y and triangles
+# arrays to triplot directly.  It would be better to use a Triangulation object
+# if the same triangulation was to be used more than once to save duplicated
+# calculations.
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.triplot(x, y, triangles, 'go-')
+plt.title('triplot of user-specified triangulation')
+plt.xlabel('Longitude (degrees)')
+plt.ylabel('Latitude (degrees)')
+
+plt.show()
Index: examples/misc/contour_manual.py
===================================================================
--- examples/misc/contour_manual.py	(revision 0)
+++ examples/misc/contour_manual.py	(revision 0)
@@ -0,0 +1,50 @@
+"""
+Example of displaying your own contour lines and polygons using ContourSet.
+"""
+import matplotlib.pyplot as plt
+from matplotlib.contour import ContourSet
+import matplotlib.cm as cm
+
+# Contour lines for each level are a list/tuple of polygons.
+lines0 = [ [[0,0],[0,4]] ]
+lines1 = [ [[2,0],[1,2],[1,3]] ]
+lines2 = [ [[3,0],[3,2]], [[3,3],[3,4]] ]  # Note two lines.
+
+# Filled contours between two levels are also a list/tuple of polygons.
+# Points can be ordered clockwise or anticlockwise.
+filled01 = [ [[0,0],[0,4],[1,3],[1,2],[2,0]] ]
+filled12 = [ [[2,0],[3,0],[3,2],[1,3],[1,2]],   # Note two polygons.
+             [[1,4],[3,4],[3,3]] ]
+
+
+plt.figure()
+
+# Filled contours using filled=True.
+cs = ContourSet(plt.gca(), [0,1,2], [filled01, filled12], filled=True, cmap=cm.bone)
+cbar = plt.colorbar(cs)
+
+# Contour lines (non-filled).
+lines = ContourSet(plt.gca(), [0,1,2], [lines0, lines1, lines2], cmap=cm.cool,
+                   linewidths=3)
+cbar.add_lines(lines)
+
+plt.axis([-0.5, 3.5, -0.5, 4.5])
+plt.title('User-specified contours')
+
+
+
+# Multiple filled contour lines can be specified in a single list of polygon
+# vertices along with a list of vertex kinds (code types) as described in the
+# Path class.  This is particularly useful for polygons with holes.
+# Here a code type of 1 is a MOVETO, and 2 is a LINETO.
+
+plt.figure()
+filled01 = [ [[0,0],[3,0],[3,3],[0,3],[1,1],[1,2],[2,2],[2,1]] ]
+kinds01  = [ [1,2,2,2,1,2,2,2] ]
+cs = ContourSet(plt.gca(), [0,1], [filled01], [kinds01], filled=True)
+cbar = plt.colorbar(cs)
+
+plt.axis([-0.5, 3.5, -0.5, 3.5])
+plt.title('User specified filled contours with holes')
+
+plt.show()
\ No newline at end of file
