diff --git a/examples/event_handling/lasso_demo.py b/examples/event_handling/lasso_demo.py
index d441a203fa3f..c56d7c1c8e70 100644
--- a/examples/event_handling/lasso_demo.py
+++ b/examples/event_handling/lasso_demo.py
@@ -4,21 +4,20 @@
 selected points
 
 This is currently a proof-of-concept implementation (though it is
-usable as is).  There will be some refinement of the API and the
-inside polygon detection routine.
+usable as is).  There will be some refinement of the API.
 """
 from matplotlib.widgets import Lasso
-from matplotlib.nxutils import points_inside_poly
 from matplotlib.colors import colorConverter
 from matplotlib.collections import RegularPolyCollection
+from matplotlib import path
 
-from matplotlib.pyplot import figure, show
+import matplotlib.pyplot as plt
 from numpy import nonzero
 from numpy.random import rand
 
-class Datum:
+class Datum(object):
     colorin = colorConverter.to_rgba('red')
-    colorout = colorConverter.to_rgba('green')
+    colorout = colorConverter.to_rgba('blue')
     def __init__(self, x, y, include=False):
         self.x = x
         self.y = y
@@ -26,7 +25,7 @@ def __init__(self, x, y, include=False):
         else: self.color = self.colorout
 
 
-class LassoManager:
+class LassoManager(object):
     def __init__(self, ax, data):
         self.axes = ax
         self.canvas = ax.figure.canvas
@@ -46,13 +45,13 @@ def __init__(self, ax, data):
         ax.add_collection(self.collection)
 
         self.cid = self.canvas.mpl_connect('button_press_event', self.onpress)
-        self.ind = None
 
     def callback(self, verts):
         facecolors = self.collection.get_facecolors()
-        ind = nonzero(points_inside_poly(self.xys, verts))[0]
-        for i in range(self.Nxy):
-            if i in ind:
+        p = path.Path(verts)
+        ind = p.contains_points(self.xys)
+        for i in range(len(self.xys)):
+            if ind[i]:
                 facecolors[i] = Datum.colorin
             else:
                 facecolors[i] = Datum.colorout
@@ -60,7 +59,7 @@ def callback(self, verts):
         self.canvas.draw_idle()
         self.canvas.widgetlock.release(self.lasso)
         del self.lasso
-        self.ind = ind
+
     def onpress(self, event):
         if self.canvas.widgetlock.locked(): return
         if event.inaxes is None: return
@@ -72,8 +71,7 @@ def onpress(self, event):
 
     data = [Datum(*xy) for xy in rand(100, 2)]
 
-    fig = figure()
-    ax = fig.add_subplot(111, xlim=(0,1), ylim=(0,1), autoscale_on=False)
+    ax = plt.axes(xlim=(0,1), ylim=(0,1), autoscale_on=False)
     lman = LassoManager(ax, data)
 
-    show()
+    plt.show()
diff --git a/lib/matplotlib/nxutils.py b/lib/matplotlib/nxutils.py
new file mode 100644
index 000000000000..078cb6afa690
--- /dev/null
+++ b/lib/matplotlib/nxutils.py
@@ -0,0 +1,50 @@
+import warnings
+
+from matplotlib import path
+
+def pnpoly(x, y, xyverts):
+    """
+    inside = pnpoly(x, y, xyverts)
+
+    Return 1 if x,y is inside the polygon, 0 otherwise.
+
+    *xyverts*
+        a sequence of x,y vertices.
+
+    A point on the boundary may be treated as inside or outside.
+
+    .. deprecated:: 1.2.0
+        Use :meth:`~matplotlib.path.Path.contains_point` instead.
+    """
+    warings.warn(
+        "nxutils is deprecated.  Use matplotlib.path.Path.contains_point"
+        " instead.",
+        DeprecationWarning)
+
+    p = path.Path(xyverts)
+    return p.contains_point(x, y)
+
+def points_inside_poly(xypoints, xyverts):
+    """
+    mask = points_inside_poly(xypoints, xyverts)
+
+    Returns a boolean ndarray, True for points inside the polygon.
+
+    *xypoints*
+        a sequence of N x,y pairs.
+
+    *xyverts*
+        sequence of x,y vertices of the polygon.
+
+    A point on the boundary may be treated as inside or outside.
+
+    .. deprecated:: 1.2.0
+        Use :meth:`~matplotlib.path.Path.contains_points` instead.
+    """
+    warnings.warn(
+        "nxutils is deprecated.  Use matplotlib.path.Path.contains_points"
+        " instead.",
+        DeprecationWarning)
+
+    p = path.Path(xyverts)
+    return p.contains_points(xypoints)
diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py
index 8ebdfe99127f..a4a3d811e394 100644
--- a/lib/matplotlib/path.py
+++ b/lib/matplotlib/path.py
@@ -12,7 +12,7 @@
 from matplotlib._path import point_in_path, get_path_extents, \
     point_in_path_collection, get_path_collection_extents, \
     path_in_path, path_intersects_path, convert_path_to_polygons, \
-    cleanup_path
+    cleanup_path, points_in_path
 from matplotlib.cbook import simple_linear_interpolation, maxdict
 from matplotlib import rcParams
 
@@ -280,12 +280,31 @@ def contains_point(self, point, transform=None, radius=0.0):
 
         If *transform* is not *None*, the path will be transformed
         before performing the test.
+
+        *radius* allows the path to be made slightly larger or
+        smaller.
         """
         if transform is not None:
             transform = transform.frozen()
         result = point_in_path(point[0], point[1], radius, self, transform)
         return result
 
+    def contains_points(self, points, transform=None, radius=0.0):
+        """
+        Returns a bool array which is *True* if the path contains the
+        corresponding point.
+
+        If *transform* is not *None*, the path will be transformed
+        before performing the test.
+
+        *radius* allows the path to be made slightly larger or
+        smaller.
+        """
+        if transform is not None:
+            transform = transform.frozen()
+        result = points_in_path(points, radius, self, transform)
+        return result
+
     def contains_path(self, path, transform=None):
         """
         Returns *True* if this path completely contains the given path.
diff --git a/src/_path.cpp b/src/_path.cpp
index 2729b3b1042a..4a39788fd1eb 100644
--- a/src/_path.cpp
+++ b/src/_path.cpp
@@ -33,6 +33,8 @@ class _path_module : public Py::ExtensionModule<_path_module>
     {
         add_varargs_method("point_in_path", &_path_module::point_in_path,
                            "point_in_path(x, y, path, trans)");
+        add_varargs_method("points_in_path", &_path_module::points_in_path,
+                           "points_in_path(points, path, trans)");
         add_varargs_method("point_on_path", &_path_module::point_on_path,
                            "point_on_path(x, y, r, path, trans)");
         add_varargs_method("get_path_extents", &_path_module::get_path_extents,
@@ -66,6 +68,7 @@ class _path_module : public Py::ExtensionModule<_path_module>
 
 private:
     Py::Object point_in_path(const Py::Tuple& args);
+    Py::Object points_in_path(const Py::Tuple& args);
     Py::Object point_on_path(const Py::Tuple& args);
     Py::Object get_path_extents(const Py::Tuple& args);
     Py::Object update_path_extents(const Py::Tuple& args);
@@ -119,16 +122,28 @@ class _path_module : public Py::ExtensionModule<_path_module>
 // Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 // _point_, returns 1 if inside, 0 if outside.
 template<class T>
-bool
-point_in_path_impl(const double tx, const double ty, T& path)
+static void
+point_in_path_impl(const void* const points_, const size_t s0,
+                   const size_t s1, const size_t n, T& path,
+                   npy_bool* const inside_flag)
 {
-    int yflag0, yflag1, inside_flag;
-    double vtx0, vty0, vtx1, vty1, sx, sy;
+    int *yflag0;
+    int yflag1;
+    double vtx0, vty0, vtx1, vty1;
+    double tx, ty;
+    double sx, sy;
     double x, y;
+    size_t i;
+    int all_done;
+    const char *const points = (const char * const)points_;
+
+    yflag0 = (int *)malloc(n * sizeof(int));
 
     path.rewind(0);
 
-    inside_flag = 0;
+    for (i = 0; i < n; ++i) {
+        inside_flag[i] = 0;
+    }
 
     unsigned code = 0;
     do
@@ -138,23 +153,25 @@ point_in_path_impl(const double tx, const double ty, T& path)
             code = path.vertex(&x, &y);
         }
 
-        sx = vtx0 = x;
-        sy = vty0 = y;
+        sx = vtx0 = vtx1 = x;
+        sy = vty0 = vty1 = y;
 
-        // get test bit for above/below X axis
-        yflag0 = (vty0 >= ty);
+        for (i = 0; i < n; ++i) {
+            ty = *(double *)(points + s0 * i + s1);
 
-        vtx1 = x;
-        vty1 = y;
+            // get test bit for above/below X axis
+            yflag0[i] = (vty0 >= ty);
+
+            inside_flag[i] = 0;
+        }
 
-        inside_flag = 0;
         do
         {
             code = path.vertex(&x, &y);
 
             // The following cases denote the beginning on a new subpath
             if (code == agg::path_cmd_stop ||
-                    (code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly)
+                (code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly)
             {
                 x = sx;
                 y = sy;
@@ -164,38 +181,42 @@ point_in_path_impl(const double tx, const double ty, T& path)
                 break;
             }
 
-            yflag1 = (vty1 >= ty);
-            // Check if endpoints straddle (are on opposite sides) of
-            // X axis (i.e. the Y's differ); if so, +X ray could
-            // intersect this edge.  The old test also checked whether
-            // the endpoints are both to the right or to the left of
-            // the test point.  However, given the faster intersection
-            // point computation used below, this test was found to be
-            // a break-even proposition for most polygons and a loser
-            // for triangles (where 50% or more of the edges which
-            // survive this test will cross quadrants and so have to
-            // have the X intersection computed anyway).  I credit
-            // Joseph Samosky with inspiring me to try dropping the
-            // "both left or both right" part of my code.
-            if (yflag0 != yflag1)
-            {
-                // Check intersection of pgon segment with +X ray.
-                // Note if >= point's X; if so, the ray hits it.  The
-                // division operation is avoided for the ">=" test by
-                // checking the sign of the first vertex wrto the test
-                // point; idea inspired by Joseph Samosky's and Mark
-                // Haigh-Hutchinson's different polygon inclusion
-                // tests.
-                if (((vty1 - ty) * (vtx0 - vtx1) >=
-                        (vtx1 - tx) * (vty0 - vty1)) == yflag1)
-                {
-                    inside_flag ^= 1;
+            for (i = 0; i < n; ++i) {
+                tx = *(double *)(points + s0 * i);
+                ty = *(double *)(points + s0 * i + s1);
+
+                yflag1 = (vty1 >= ty);
+                // Check if endpoints straddle (are on opposite sides) of
+                // X axis (i.e. the Y's differ); if so, +X ray could
+                // intersect this edge.  The old test also checked whether
+                // the endpoints are both to the right or to the left of
+                // the test point.  However, given the faster intersection
+                // point computation used below, this test was found to be
+                // a break-even proposition for most polygons and a loser
+                // for triangles (where 50% or more of the edges which
+                // survive this test will cross quadrants and so have to
+                // have the X intersection computed anyway).  I credit
+                // Joseph Samosky with inspiring me to try dropping the
+                // "both left or both right" part of my code.
+                if (yflag0[i] != yflag1) {
+                    // Check intersection of pgon segment with +X ray.
+                    // Note if >= point's X; if so, the ray hits it.  The
+                    // division operation is avoided for the ">=" test by
+                    // checking the sign of the first vertex wrto the test
+                    // point; idea inspired by Joseph Samosky's and Mark
+                    // Haigh-Hutchinson's different polygon inclusion
+                    // tests.
+                    if (((vty1 - ty) * (vtx0 - vtx1) >=
+                         (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
+                        inside_flag[i] ^= 1;
+                    }
                 }
+
+                // Move to the next pair of vertices, retaining info as
+                // possible.
+                yflag0[i] = yflag1;
             }
 
-            // Move to the next pair of vertices, retaining info as
-            // possible.
-            yflag0 = yflag1;
             vtx0 = vtx1;
             vty0 = vty1;
 
@@ -205,38 +226,55 @@ point_in_path_impl(const double tx, const double ty, T& path)
         while (code != agg::path_cmd_stop &&
                (code & agg::path_cmd_end_poly) != agg::path_cmd_end_poly);
 
-        yflag1 = (vty1 >= ty);
-        if (yflag0 != yflag1)
-        {
-            if (((vty1 - ty) * (vtx0 - vtx1) >=
-                    (vtx1 - tx) * (vty0 - vty1)) == yflag1)
-            {
-                inside_flag ^= 1;
+        all_done = 1;
+        for (i = 0; i < n; ++i) {
+            tx = *(double *)(points + s0 * i);
+            ty = *(double *)(points + s0 * i + s1);
+
+            yflag1 = (vty1 >= ty);
+            if (yflag0[i] != yflag1) {
+                if (((vty1 - ty) * (vtx0 - vtx1) >=
+                     (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
+                    inside_flag[i] ^= 1;
+                }
+            }
+
+            if (inside_flag[i] == 0) {
+                all_done = 0;
             }
         }
 
-        if (inside_flag != 0)
-        {
-            return true;
+        if (all_done) {
+            goto exit;
         }
     }
     while (code != agg::path_cmd_stop);
 
-    return (inside_flag != 0);
+ exit:
+
+    free(yflag0);
 }
 
-inline bool
-point_in_path(double x, double y, double r, PathIterator& path,
-              const agg::trans_affine& trans)
+inline void
+points_in_path(const void* const points, const size_t s0,
+               const size_t s1, const size_t n,
+               const double r, PathIterator& path,
+               const agg::trans_affine& trans,
+               npy_bool* result)
 {
     typedef agg::conv_transform<PathIterator> transformed_path_t;
     typedef PathNanRemover<transformed_path_t> no_nans_t;
     typedef agg::conv_curve<no_nans_t> curve_t;
     typedef agg::conv_contour<curve_t> contour_t;
 
+    size_t i;
+    for (i = 0; i < n; ++i) {
+        result[i] = 0;
+    }
+
     if (path.total_vertices() < 3)
     {
-        return false;
+        return;
     }
 
     transformed_path_t trans_path(path, trans);
@@ -244,12 +282,29 @@ point_in_path(double x, double y, double r, PathIterator& path,
     curve_t curved_path(no_nans_path);
     contour_t contoured_path(curved_path);
     contoured_path.width(fabs(r));
-    return point_in_path_impl(x, y, contoured_path);
+    point_in_path_impl(points, s0, s1, n, contoured_path, result);
 }
 
 inline bool
-point_on_path(double x, double y, double r, PathIterator& path,
-              const agg::trans_affine& trans)
+point_in_path(const double x, const double y, const double r,
+              PathIterator& path, const agg::trans_affine& trans)
+{
+    double points[2];
+    npy_bool result;
+
+    points[0] = x;
+    points[1] = y;
+
+    points_in_path(points, 0, sizeof(double), 1, r, path, trans, &result);
+    return result;
+}
+
+inline void
+points_on_path(const void* const points, const size_t s0,
+               const size_t s1, const size_t n,
+               const double r, PathIterator& path,
+               const agg::trans_affine& trans,
+               npy_bool* result)
 {
     typedef agg::conv_transform<PathIterator> transformed_path_t;
     typedef PathNanRemover<transformed_path_t> no_nans_t;
@@ -261,32 +316,74 @@ point_on_path(double x, double y, double r, PathIterator& path,
     curve_t curved_path(nan_removed_path);
     stroke_t stroked_path(curved_path);
     stroked_path.width(r * 2.0);
-    return point_in_path_impl(x, y, stroked_path);
+    point_in_path_impl(points, s0, s1, n, stroked_path, result);
+}
+
+inline bool
+point_on_path(const double x, const double y, const double r,
+              PathIterator& path, const agg::trans_affine& trans)
+{
+    double points[2];
+    npy_bool result;
+
+    points[0] = x;
+    points[1] = y;
+
+    points_on_path(points, 0, sizeof(double), 1, r, path, trans, &result);
+    return result;
 }
 
 Py::Object
 _path_module::point_in_path(const Py::Tuple& args)
 {
-    args.verify_length(5);
-
     double x = Py::Float(args[0]);
     double y = Py::Float(args[1]);
     double r = Py::Float(args[2]);
     PathIterator path(args[3]);
     agg::trans_affine trans = py_to_agg_transformation_matrix(args[4].ptr(), false);
 
-    if (::point_in_path(x, y, r, path, trans))
-    {
+    if (::point_in_path(x, y, r, path, trans)) {
         return Py::Int(1);
     }
     return Py::Int(0);
 }
 
 Py::Object
-_path_module::point_on_path(const Py::Tuple& args)
+_path_module::points_in_path(const Py::Tuple& args)
 {
-    args.verify_length(5);
+    args.verify_length(4);
+
+    npy_intp n;
+    PyArrayObject* points_array;
+    points_array = (PyArrayObject*)PyArray_FromObject(args[0].ptr(), PyArray_DOUBLE, 2, 2);
+    if (points_array == NULL || PyArray_DIM(points_array, 1) != 2) {
+        throw Py::TypeError(
+            "Argument 0 to points_in_path must be an Nx2 numpy array");
+
+    }
+    double r = Py::Float(args[1]);
+    PathIterator path(args[2]);
+    agg::trans_affine trans = py_to_agg_transformation_matrix(args[3].ptr(), false);
+
+    n = PyArray_DIM(points_array, 0);
+    PyObject* result = PyArray_ZEROS(1, &n, PyArray_BOOL, 0);
+    if (result == NULL) {
+        throw Py::MemoryError("Could not allocate memory for result");
+    }
+
+    ::points_in_path(PyArray_DATA(points_array),
+                     PyArray_STRIDE(points_array, 0),
+                     PyArray_STRIDE(points_array, 1),
+                     n, r, path, trans,
+                     (npy_bool *)PyArray_DATA(result));
+    Py_DECREF(points_array);
+
+    return Py::Object(result, true);;
+}
 
+Py::Object
+_path_module::point_on_path(const Py::Tuple& args)
+{
     double x = Py::Float(args[0]);
     double y = Py::Float(args[1]);
     double r = Py::Float(args[2]);