diff --git a/skimage/_shared/__init__.py b/skimage/_shared/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/skimage/_shared/geometry.pxd b/skimage/_shared/geometry.pxd
new file mode 100644
index 00000000..b228b23a
--- /dev/null
+++ b/skimage/_shared/geometry.pxd
@@ -0,0 +1,7 @@
+
+cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp,
+                                           double x, double y)
+
+cdef void points_in_polygon(int nr_verts, double *xp, double *yp,
+                            int nr_points, double *x, double *y,
+                            unsigned char *result)
diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx
new file mode 100644
index 00000000..3f4850b0
--- /dev/null
+++ b/skimage/_shared/geometry.pyx
@@ -0,0 +1,54 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+
+cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp,
+                                           double x, double y):
+    """Test whether point lies inside a polygon.
+
+    Parameters
+    ----------
+    nr_verts : int
+        Number of vertices of polygon.
+    xp, yp : double array
+        Coordinates of polygon with length nr_verts.
+    x, y : double
+        Coordinates of point.
+    """
+    cdef int i
+    cdef unsigned char c = 0
+    cdef int j = nr_verts - 1
+    for i in range(nr_verts):
+        if (
+            (((yp[i] <= y) and (y < yp[j])) or
+            ((yp[j] <= y) and (y < yp[i])))
+            and (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])
+        ):
+            c = not c
+        j = i
+    return c
+
+
+cdef void points_in_polygon(int nr_verts, double *xp, double *yp,
+                            int nr_points, double *x, double *y,
+                            unsigned char *result):
+    """Test whether points lie inside a polygon.
+
+    Parameters
+    ----------
+    nr_verts : int
+        Number of vertices of polygon.
+    xp, yp : double array
+        Coordinates of polygon with length nr_verts.
+    nr_points : int
+        Number of points to test.
+    x, y : double array
+        Coordinates of points.
+    result : unsigned char array
+        Test results for each point.
+    """
+    cdef int n
+    for n in range(nr_points):
+        result[n] = point_in_polygon(nr_verts, xp, yp, x[n], y[n])
diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
new file mode 100644
index 00000000..c883d00d
--- /dev/null
+++ b/skimage/_shared/interpolation.pxd
@@ -0,0 +1,13 @@
+
+cdef inline double nearest_neighbour(double* image, int rows, int cols,
+                                     double r, double c, char mode,
+                                     double cval=*)
+
+cdef inline double bilinear_interpolation(double* image, int rows, int cols,
+                                          double r, double c, char mode,
+                                          double cval=*)
+
+cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
+                             char mode, double cval=*)
+
+cdef inline int coord_map(int dim, int coord, char mode)
diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
new file mode 100644
index 00000000..137004e5
--- /dev/null
+++ b/skimage/_shared/interpolation.pyx
@@ -0,0 +1,128 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+from libc.math cimport ceil, floor, round
+
+
+cdef inline double nearest_neighbour(double* image, int rows, int cols,
+                                     double r, double c, char mode,
+                                     double cval=0):
+    """Nearest neighbour interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols: int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'M'}
+        Wrapping mode. Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for constant mode.
+
+    """
+
+    return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
+                     mode, cval)
+
+
+cdef inline double bilinear_interpolation(double* image, int rows, int cols,
+                                          double r, double c, char mode,
+                                          double cval=0):
+    """Bilinear interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols: int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'M'}
+        Wrapping mode. Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for constant mode.
+
+    """
+    cdef double dr, dc
+    cdef int minr, minc, maxr, maxc
+
+    minr = <int>floor(r)
+    minc = <int>floor(c)
+    maxr = <int>ceil(r)
+    maxc = <int>ceil(c)
+    dr = r - minr
+    dc = c - minc
+    top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
+          + dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
+    bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
+             + dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
+    return (1 - dr) * top + dr * bottom
+
+
+cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
+                             char mode, double cval=0):
+    """Get a pixel from the image, taking wrapping mode into consideration.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols: int
+        Shape of image.
+    r, c : int
+        Position at which to get the pixel.
+    mode : {'C', 'W', 'M'}
+        Wrapping mode. Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for constant mode.
+
+    """
+    if mode == 'C':
+        if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
+            return cval
+        else:
+            return image[r * cols + c]
+    else:
+        return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
+
+
+cdef inline int coord_map(int dim, int coord, char mode):
+    """
+    Wrap a coordinate, according to a given mode.
+
+    Parameters
+    ----------
+    dim : int
+        Maximum coordinate.
+    coord : int
+        Coord provided by user.  May be < 0 or > dim.
+    mode : {'W', 'M'}
+        Whether to wrap or mirror the coordinate if it
+        falls outside [0, dim).
+
+    """
+    dim = dim - 1
+    if mode == 'M': # mirror
+        if (coord < 0):
+            # How many times times does the coordinate wrap?
+            if (<int>(-coord / dim) % 2 != 0):
+                return dim - <int>(-coord % dim)
+            else:
+                return <int>(-coord % dim)
+        elif (coord > dim):
+            if (<int>(coord / dim) % 2 != 0):
+                return <int>(dim - (coord % dim))
+            else:
+                return <int>(coord % dim)
+    elif mode == 'W': # wrap
+        if (coord < 0):
+            return <int>(dim - (-coord % dim))
+        elif (coord > dim):
+            return <int>(coord % dim)
+
+    return coord
diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py
new file mode 100644
index 00000000..81113656
--- /dev/null
+++ b/skimage/_shared/setup.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+import os
+
+from skimage._build import cython
+
+base_path = os.path.abspath(os.path.dirname(__file__))
+
+
+def configuration(parent_package='', top_path=None):
+    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+
+    config = Configuration('_shared', parent_package, top_path)
+    config.add_data_dir('tests')
+
+    cython(['geometry.pyx'], working_path=base_path)
+    cython(['interpolation.pyx'], working_path=base_path)
+    cython(['transform.pyx'], working_path=base_path)
+
+    config.add_extension('geometry', sources=['geometry.c'])
+    config.add_extension('interpolation', sources=['interpolation.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('transform', sources=['transform.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
+    return config
+
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(maintainer='Scikits-image Developers',
+          author='Scikits-image Developers',
+          maintainer_email='scikits-image@googlegroups.com',
+          description='Transforms',
+          url='https://github.com/scikits-image/scikits-image',
+          license='SciPy License (BSD Style)',
+          **(configuration(top_path='').todict())
+          )
diff --git a/skimage/_shared/transform.pxd b/skimage/_shared/transform.pxd
new file mode 100644
index 00000000..0edc22a4
--- /dev/null
+++ b/skimage/_shared/transform.pxd
@@ -0,0 +1,5 @@
+cimport numpy as cnp
+
+
+cdef float integrate(cnp.ndarray[float, ndim=2,  mode="c"] sat,
+                     int r0, int c0, int r1, int c1)
diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx
new file mode 100644
index 00000000..ba0efc71
--- /dev/null
+++ b/skimage/_shared/transform.pyx
@@ -0,0 +1,44 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+cimport numpy as cnp
+
+
+cdef float integrate(cnp.ndarray[float, ndim=2,  mode="c"] sat,
+                     int r0, int c0, int r1, int c1):
+    """
+    Using a summed area table / integral image, calculate the sum
+    over a given window.
+
+    This function is the same as the `integrate` function in
+    `skimage.transform.integrate`, but this Cython version significantly
+    speeds up the code.
+
+    Parameters
+    ----------
+    sat : ndarray of float
+        Summed area table / integral image.
+    r0, c0 : int
+        Top-left corner of block to be summed.
+    r1, c1 : int
+        Bottom-right corner of block to be summed.
+
+    Returns
+    -------
+    S : int
+        Sum over the given window.
+    """
+    cdef float S = 0
+
+    S += sat[r1, c1]
+
+    if (r0 - 1 >= 0) and (c0 - 1 >= 0):
+        S += sat[r0 - 1, c0 - 1]
+
+    if (r0 - 1 >= 0):
+        S -= sat[r0 - 1, c1]
+
+    if (c0 - 1 >= 0):
+        S -= sat[r1, c0 - 1]
+    return S
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index 9dc1d307..fbd525fa 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -3,11 +3,7 @@ import math
 from libc.math cimport sqrt
 cimport numpy as np
 cimport cython
-
-
-cdef extern from "../morphology/_pnpoly.h":
-     int pnpoly(int nr_verts, double *xp, double *yp,
-                double x, double y)
+from skimage._shared.geometry cimport point_in_polygon
 
 
 @cython.boundscheck(False)
@@ -119,7 +115,7 @@ def polygon(y, x, shape=None):
 
     for r in range(minr, maxr+1):
         for c in range(minc, maxc+1):
-            if pnpoly(nr_verts, cptr, rptr, c, r):
+            if point_in_polygon(nr_verts, cptr, rptr, c, r):
                 rr.append(r)
                 cc.append(c)
 
diff --git a/skimage/draw/setup.py b/skimage/draw/setup.py
index 4503d664..59067f9a 100644
--- a/skimage/draw/setup.py
+++ b/skimage/draw/setup.py
@@ -15,7 +15,7 @@ def configuration(parent_package='', top_path=None):
     cython(['_draw.pyx'], working_path=base_path)
 
     config.add_extension('_draw', sources=['_draw.c'],
-                         include_dirs=[get_numpy_include_dirs()])
+                         include_dirs=[get_numpy_include_dirs(), '../shared'])
 
     return config
 
diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx
index b83761a8..58d48524 100644
--- a/skimage/feature/_template.pyx
+++ b/skimage/feature/_template.pyx
@@ -35,51 +35,8 @@ cimport numpy as np
 import numpy as np
 from scipy.signal import fftconvolve
 from skimage.transform import integral
-
-
-cdef extern from "math.h":
-    float sqrt(float x)
-    float fabs(float x)
-
-
-@cython.boundscheck(False)
-cdef float integrate(np.ndarray[float, ndim=2,  mode="c"] sat,
-        int r0, int c0, int r1, int c1):
-    """
-    Using a summed area table / integral image, calculate the sum
-    over a given window.
-
-    This function is the same as the `integrate` function in
-    `skimage.transform.integrate`, but this Cython version significantly
-    speeds up the code.
-
-    Parameters
-    ----------
-    sat : ndarray of float
-        Summed area table / integral image.
-    r0, c0 : int
-        Top-left corner of block to be summed.
-    r1, c1 : int
-        Bottom-right corner of block to be summed.
-
-    Returns
-    -------
-    S : int
-        Sum over the given window.
-    """
-    cdef float S = 0
-
-    S += sat[r1, c1]
-
-    if (r0 - 1 >= 0) and (c0 - 1 >= 0):
-        S += sat[r0 - 1, c0 - 1]
-
-    if (r0 - 1 >= 0):
-        S -= sat[r0 - 1, c1]
-
-    if (c0 - 1 >= 0):
-        S -= sat[r1, c0 - 1]
-    return S
+from libc.math cimport sqrt, fabs
+from skimage._shared.transform cimport integrate
 
 
 @cython.boundscheck(False)
diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx
index ca51e177..0e82c306 100644
--- a/skimage/feature/_texture.pyx
+++ b/skimage/feature/_texture.pyx
@@ -1,11 +1,11 @@
-#cython: cdivison=True
+#cython: cdivision=True
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
 import numpy as np
 cimport numpy as np
 from libc.math cimport sin, cos, abs
-from skimage.transform._project cimport bilinear_interpolation
+from skimage._shared.interpolation cimport bilinear_interpolation
 
 
 def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py
index 0b0b80bd..9c9074be 100644
--- a/skimage/feature/setup.py
+++ b/skimage/feature/setup.py
@@ -16,10 +16,9 @@ def configuration(parent_package='', top_path=None):
     cython(['_template.pyx'], working_path=base_path)
 
     config.add_extension('_texture', sources=['_texture.c'],
-                         include_dirs=[get_numpy_include_dirs(),
-                                       '../transform'])
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
     config.add_extension('_template', sources=['_template.c'],
-                         include_dirs=[get_numpy_include_dirs()])
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config
 
diff --git a/skimage/graph/_spath.pyx b/skimage/graph/_spath.pyx
index f342021d..1624bc7d 100644
--- a/skimage/graph/_spath.pyx
+++ b/skimage/graph/_spath.pyx
@@ -2,9 +2,8 @@
 
 import _mcp
 cimport _mcp
+from libc.math cimport fabs
 
-cdef extern from "math.h":
-    double fabs(double f)
 
 cdef class MCP_Diff(_mcp.MCP):
     """MCP_Diff(costs, offsets=None, fully_connected=True)
diff --git a/skimage/io/_plugins/_colormixer.pyx b/skimage/io/_plugins/_colormixer.pyx
index 5d0dd1b9..ace45c94 100644
--- a/skimage/io/_plugins/_colormixer.pyx
+++ b/skimage/io/_plugins/_colormixer.pyx
@@ -8,15 +8,10 @@ integers, so currently the only way to clip results efficiently
 one.
 
 """
-
+import cython
 import numpy as np
 cimport numpy as np
-
-import cython
-
-cdef extern from "math.h":
-    float exp(float) nogil
-    float pow(float, float) nogil
+from libc.math cimport exp, pow
 
 
 @cython.boundscheck(False)
@@ -189,7 +184,6 @@ def sigmoid_gamma(np.ndarray[np.uint8_t, ndim=3] img,
                 img[i,j,2] = lut[stateimg[i,j,2]]
 
 
-
 @cython.boundscheck(False)
 def gamma(np.ndarray[np.uint8_t, ndim=3] img,
           np.ndarray[np.uint8_t, ndim=3] stateimg,
@@ -219,7 +213,6 @@ def gamma(np.ndarray[np.uint8_t, ndim=3] img,
                 img[i,j,2] = lut[stateimg[i,j,2]]
 
 
-
 @cython.cdivision(True)
 cdef void rgb_2_hsv(float* RGB, float* HSV) nogil:
     cdef float R, G, B, H, S, V, MAX, MIN
@@ -283,6 +276,7 @@ cdef void rgb_2_hsv(float* RGB, float* HSV) nogil:
     HSV[1] = S
     HSV[2] = V
 
+
 @cython.cdivision(True)
 cdef void hsv_2_rgb(float* HSV, float* RGB) nogil:
     cdef float H, S, V
@@ -388,6 +382,7 @@ def py_hsv_2_rgb(H, S, V):
 
     return (R, G, B)
 
+
 def py_rgb_2_hsv(R, G, B):
     '''Convert an HSV value to RGB.
 
@@ -561,9 +556,3 @@ def hsv_multiply(np.ndarray[np.uint8_t, ndim=3] img,
                 img[i, j, 0] = <np.uint8_t>RGB[0]
                 img[i, j, 1] = <np.uint8_t>RGB[1]
                 img[i, j, 2] = <np.uint8_t>RGB[2]
-
-
-
-
-
-
diff --git a/skimage/morphology/_pnpoly.h b/skimage/morphology/_pnpoly.h
deleted file mode 100644
index 95c89bcb..00000000
--- a/skimage/morphology/_pnpoly.h
+++ /dev/null
@@ -1,72 +0,0 @@
-/* `pnpoly` is from
-
-   http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
-
-   Copyright (c) 1970-2003, Wm. Randolph Franklin
-
-   Permission is hereby granted, free of charge, to any person
-   obtaining a copy of this software and associated documentation
-   files (the "Software"), to deal in the Software without
-   restriction, including without limitation the rights to use, copy,
-   modify, merge, publish, distribute, sublicense, and/or sell copies
-   of the Software, and to permit persons to whom the Software is
-   furnished to do so, subject to the following conditions:
-
-    1. Redistributions of source code must retain the above copyright
-       notice, this list of conditions and the following disclaimers.
-    2. Redistributions in binary form must reproduce the above
-       copyright notice in the documentation and/or other materials
-       provided with the distribution.
-    3. The name of W. Randolph Franklin may not be used to endorse or
-       promote products derived from this Software without specific
-       prior written permission.
-
-   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
-   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
-   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
-   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-   SOFTWARE. */
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-unsigned char pnpoly(int nr_verts, double *xp, double *yp, double x, double y)
-{
-    int i, j;
-    unsigned char c = 0;
-    for (i = 0, j = nr_verts-1; i < nr_verts; j = i++) {
-        if ((((yp[i]<=y) && (y<yp[j])) ||
-             ((yp[j]<=y) && (y<yp[i]))) &&
-            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
-	    
-	    c = !c;
-    }
-    return c;
-}
-
-void npnpoly(int nr_verts, double *xp, double *yp,
-	     int nr_points, double *x, double *y,
-	     unsigned char *result)
-/*
- * For N provided points, calculate whether they are in
- * the polygon defined by vertices *xp, *yp.
- *
- * nr_verts : number of vertices
- * *xp, *yp : x and y coordinates of vertices
- * nr_points : number of data points provided
- * *x, *y : data points
- */
-{
-    unsigned char n = 0;
-    for (n = 0; n < nr_points; n++) {
-	result[n] = pnpoly(nr_verts, xp, yp, x[n], y[n]);
-    }
-}
-
-#ifdef __cplusplus
-}
-#endif
diff --git a/skimage/morphology/_pnpoly.pyx b/skimage/morphology/_pnpoly.pyx
index c492eb0d..7deb6a50 100644
--- a/skimage/morphology/_pnpoly.pyx
+++ b/skimage/morphology/_pnpoly.pyx
@@ -2,14 +2,7 @@
 
 cimport numpy as np
 import numpy as np
-
-cdef extern from "_pnpoly.h":
-     int pnpoly(int nr_verts, double *xp, double *yp,
-                double x, double y)
-
-     void npnpoly(int nr_verts, double *xp, double *yp,
-                  int nr_points, double *x, double *y,
-                  unsigned char *result)
+from skimage._shared.geometry cimport point_in_polygon, points_in_polygon
 
 
 def grid_points_inside_poly(shape, verts):
@@ -49,45 +42,45 @@ def grid_points_inside_poly(shape, verts):
 
     for m in range(M):
         for n in range(N):
-            out[m, n] = pnpoly(V, <double*>vx.data, <double*>vy.data, m, n)
+            out[m, n] = point_in_polygon(V, <double*>vx.data, <double*>vy.data, m, n)
 
     return out.view(bool)
-    
+
 
 def points_inside_poly(points, verts):
-     """Test whether points lie inside a polygon.
+    """Test whether points lie inside a polygon.
 
-     Parameters
-     ----------
-     points : (N, 2) array
-         Input points, ``(x, y)``.
-     verts : (M, 2) array
-         Vertices of the polygon, sorted either clockwise or anti-clockwise.
-         The first point may (but does not need to be) duplicated.
+    Parameters
+    ----------
+    points : (N, 2) array
+     Input points, ``(x, y)``.
+    verts : (M, 2) array
+     Vertices of the polygon, sorted either clockwise or anti-clockwise.
+     The first point may (but does not need to be) duplicated.
 
-     Returns
-     -------
-     mask : (N,) array of bool
-         True if corresponding point is inside the polygon.
+    Returns
+    -------
+    mask : (N,) array of bool
+     True if corresponding point is inside the polygon.
 
-     """
-     cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy
+    """
+    cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy
 
-     points = np.asarray(points)
-     verts = np.asarray(verts)
+    points = np.asarray(points)
+    verts = np.asarray(verts)
 
-     x = points[:, 0].astype(np.double)
-     y = points[:, 1].astype(np.double)
+    x = points[:, 0].astype(np.double)
+    y = points[:, 1].astype(np.double)
 
-     vx = verts[:, 0].astype(np.double)
-     vy = verts[:, 1].astype(np.double)
+    vx = verts[:, 0].astype(np.double)
+    vy = verts[:, 1].astype(np.double)
 
-     cdef np.ndarray[np.uint8_t, ndim=1] out = \
-          np.zeros(x.shape[0], dtype=np.uint8)
-     
-     npnpoly(vx.shape[0], <double*>vx.data, <double*>vy.data,
-             x.shape[0], <double*>x.data, <double*>y.data,
-             <unsigned char*>out.data)
+    cdef np.ndarray[np.uint8_t, ndim=1] out = \
+      np.zeros(x.shape[0], dtype=np.uint8)
 
-     return out.astype(bool)
+    points_in_polygon(vx.shape[0], <double*>vx.data, <double*>vy.data,
+         x.shape[0], <double*>x.data, <double*>y.data,
+         <unsigned char*>out.data)
+
+    return out.astype(bool)
 
diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py
index f5dd7f19..02c39ae7 100644
--- a/skimage/morphology/setup.py
+++ b/skimage/morphology/setup.py
@@ -29,7 +29,7 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_skeletonize_cy', sources=['_skeletonize_cy.c'],
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_pnpoly', sources=['_pnpoly.c'],
-                         include_dirs=[get_numpy_include_dirs()])
+                         include_dirs=[get_numpy_include_dirs(), '../shared'])
     config.add_extension('_convex_hull', sources=['_convex_hull.c'],
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'],
diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx
index 57009ad1..57be1040 100644
--- a/skimage/segmentation/_quickshift.pyx
+++ b/skimage/segmentation/_quickshift.pyx
@@ -1,6 +1,7 @@
 import numpy as np
 cimport numpy as np
 cimport cython
+from libc.math cimport exp, sqrt
 
 from itertools import product
 from scipy import ndimage
@@ -9,11 +10,6 @@ from ..util import img_as_float
 from ..color import rgb2lab
 
 
-cdef extern from "math.h":
-    double exp(double)
-    double sqrt(double)
-
-
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
diff --git a/skimage/setup.py b/skimage/setup.py
index 02c0b52d..1082ba07 100644
--- a/skimage/setup.py
+++ b/skimage/setup.py
@@ -6,6 +6,7 @@ def configuration(parent_package='', top_path=None):
 
     config = Configuration('skimage', parent_package, top_path)
 
+    config.add_subpackage('_shared')
     config.add_subpackage('color')
     config.add_subpackage('data')
     config.add_subpackage('draw')
diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx
index b34b28e2..906e4464 100644
--- a/skimage/transform/_hough_transform.pyx
+++ b/skimage/transform/_hough_transform.pyx
@@ -2,27 +2,20 @@ cimport cython
 import numpy as np
 cimport numpy as np
 from random import randint
+from libc.math cimport abs, fabs, sqrt, ceil, floor, round
+from libc.stdlib cimport rand
+
+
 np.import_array()
 
-cdef extern from "stdlib.h":
-    int rand()
-        
-cdef extern from "math.h":
-    int abs(int)
-    double fabs(double)
-    double sqrt(double)
-    double ceil(double)
-    double floor(double)
-
-cdef double round(double val):
-    return floor(val + 0.5);
 
 cdef double PI_2 = 1.5707963267948966
 cdef double NEG_PI_2 = -PI_2
 
+
 @cython.boundscheck(False)
 def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
-    
+
     if img.ndim != 2:
         raise ValueError('The input image must be 2D.')
 
@@ -31,7 +24,7 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     cdef np.ndarray[ndim=1, dtype=np.double_t] stheta
 
     if theta is None:
-        theta = np.linspace(PI_2, NEG_PI_2, 180) 
+        theta = np.linspace(PI_2, NEG_PI_2, 180)
 
     ctheta = np.cos(theta)
     stheta = np.sin(theta)
@@ -39,14 +32,14 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     # compute the bins and allocate the accumulator array
     cdef np.ndarray[ndim=2, dtype=np.uint64_t] accum
     cdef np.ndarray[ndim=1, dtype=np.double_t] bins
-    cdef int max_distance, offset 
+    cdef int max_distance, offset
 
-    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
+    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
                                        img.shape[1] * img.shape[1])))
     accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
     bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
     offset = max_distance / 2
-    
+
     # compute the nonzero indexes
     cdef np.ndarray[ndim=1, dtype=np.npy_intp] x_idxs, y_idxs
     y_idxs, x_idxs = np.PyArray_Nonzero(img)
@@ -58,7 +51,7 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
     nthetas = theta.shape[0]
     for i in range(nidxs):
         x = x_idxs[i]
-        y = y_idxs[i]   
+        y = y_idxs[i]
         for j in range(nthetas):
             accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
             accum[accum_idx, j] += 1
@@ -94,7 +87,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
     # maximum line number cutoff
     cdef int lines_max = 2 ** 15
     cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count
-    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] + 
+    max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
                                        img.shape[1] * img.shape[1])))
     accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64)
     offset = max_distance / 2
@@ -114,11 +107,11 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
         # select random non-zero point
         count = len(points)
         if count == 0:
-            break        
+            break
         index = rand() % (count)
         x = points[index][0]
         y = points[index][1]
-        del points[index]        
+        del points[index]
         # if previously eliminated, skip
         if not mask[y, x]:
             continue
@@ -147,7 +140,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                 dx0 = 1
             else:
                 dx0 = -1
-            dy0 = <int>round(b * (1 << shift) / fabs(a)) 
+            dy0 = <int>round(b * (1 << shift) / fabs(a))
             y0 = (y0 << shift) + (1 << (shift - 1))
         else:
             if b > 0:
@@ -156,7 +149,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                 dy0 = -1
             dx0 = <int>round(a * (1 << shift) / fabs(b))
             x0 = (x0 << shift) + (1 << (shift - 1))
-  
+
         # pass 1: walk the line, merging lines less than specified gap length
         for k in range(2):
             gap = 0
@@ -208,9 +201,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
                     x1 = px >> shift
                     y1 = py
                 # if non-zero point found, continue the line
-                if mask[y1, x1]:  
-                    if good_line:                  
-                        accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset                    
+                if mask[y1, x1]:
+                    if good_line:
+                        accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
                         accum[accum_idx, max_theta] -= 1
                         mask[y1, x1] = 0
                 # exit when the point is the line end
diff --git a/skimage/transform/_project.pxd b/skimage/transform/_project.pxd
deleted file mode 100644
index dee51026..00000000
--- a/skimage/transform/_project.pxd
+++ /dev/null
@@ -1,7 +0,0 @@
-cimport numpy as np
-import numpy as np
-
-
-cdef inline double bilinear_interpolation(double* image, int rows, int cols,
-                                          double r, double c, char mode,
-                                          double cval=*)
diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx
index 3f31f229..4df6e001 100644
--- a/skimage/transform/_project.pyx
+++ b/skimage/transform/_project.pyx
@@ -1,111 +1,12 @@
-#cython: cdivison=True
+#cython: cdivision=True
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
 
 cimport numpy as np
 import numpy as np
-from cython.operator import dereference
-from libc.math cimport ceil, floor
-
-
-cdef inline double bilinear_interpolation(double* image, int rows, int cols,
-                                          double r, double c, char mode,
-                                          double cval=0):
-    """Bilinear interpolation at a given position in the image.
-
-    Parameters
-    ----------
-    image : double array
-        Input image.
-    rows, cols: int
-        Shape of image.
-    r, c : int
-        Position at which to interpolate.
-    mode : {'C', 'W', 'M'}
-        Wrapping mode. Constant, Wrap or Mirror.
-    cval : double
-        Constant value to use for constant mode.
-
-    """
-    cdef double dr, dc
-    cdef int minr, minc, maxr, maxc
-
-    minr = <int>floor(r)
-    minc = <int>floor(c)
-    maxr = <int>ceil(r)
-    maxc = <int>ceil(c)
-    dr = r - minr
-    dc = c - minc
-    top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
-          + dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
-    bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
-             + dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
-    return (1 - dr) * top + dr * bottom
-
-
-cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
-                             char mode, double cval=0):
-    """Get a pixel from the image, taking wrapping mode into consideration.
-
-    Parameters
-    ----------
-    image : double array
-        Input image.
-    rows, cols: int
-        Shape of image.
-    r, c : int
-        Position at which to get the pixel.
-    mode : {'C', 'W', 'M'}
-        Wrapping mode. Constant, Wrap or Mirror.
-    cval : double
-        Constant value to use for constant mode.
-
-    """
-    if mode == 'C':
-        if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
-            return cval
-        else:
-            return image[r * cols + c]
-    else:
-        return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
-
-
-cdef inline int coord_map(int dim, int coord, char mode):
-    """
-    Wrap a coordinate, according to a given mode.
-
-    Parameters
-    ----------
-    dim : int
-        Maximum coordinate.
-    coord : int
-        Coord provided by user.  May be < 0 or > dim.
-    mode : {'W', 'M'}
-        Whether to wrap or mirror the coordinate if it
-        falls outside [0, dim).
-
-    """
-    dim = dim - 1
-    if mode == 'M': # mirror
-        if (coord < 0):
-            # How many times times does the coordinate wrap?
-            if (<int>(-coord / dim) % 2 != 0):
-                return dim - <int>(-coord % dim)
-            else:
-                return <int>(-coord % dim)
-        elif (coord > dim):
-            if (<int>(coord / dim) % 2 != 0):
-                return <int>(dim - (coord % dim))
-            else:
-                return <int>(coord % dim)
-    elif mode == 'W': # wrap
-        if (coord < 0):
-            return <int>(dim - (-coord % dim))
-        elif (coord > dim):
-            return <int>(coord % dim)
-
-    return coord
+from skimage._shared.interpolation cimport (nearest_neighbour,
+                                            bilinear_interpolation)
 
 
 cdef inline _matrix_transform(double x, double y, double* H, double *x_,
@@ -132,7 +33,7 @@ cdef inline _matrix_transform(double x, double y, double* H, double *x_,
     y_[0] = yy / zz
 
 
-def homography(np.ndarray image, np.ndarray H, output_shape=None,
+def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
                mode='constant', double cval=0):
     """
     Projective transformation (homography).
@@ -167,6 +68,10 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
         Transformation matrix H that defines the homography.
     output_shape : tuple (rows, cols)
         Shape of the output image generated.
+    order : {0, 1}
+        Order of interpolation::
+        * 0: Nearest-neighbour interpolation.
+        * 1: Bilinear interpolation (default).
     mode : {'constant', 'mirror', 'wrap'}
         How to handle values outside the image borders.
     cval : string
@@ -175,13 +80,15 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
 
     """
 
-    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image.astype(np.double)
+    cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \
+         np.ascontiguousarray(image, dtype=np.double)
     cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \
          np.ascontiguousarray(np.linalg.inv(H))
 
     if mode not in ('constant', 'wrap', 'mirror'):
         raise ValueError("Invalid mode specified.  Please use "
                          "`constant`, `wrap` or `mirror`.")
+    cdef char mode_c
     if mode == 'constant':
         mode_c = ord('C')
     elif mode == 'wrap':
@@ -189,6 +96,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     elif mode == 'mirror':
         mode_c = ord('M')
 
+    cdef int out_r, out_c
     if output_shape is None:
         out_r = img.shape[0]
         out_c = img.shape[1]
@@ -207,7 +115,11 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     for tfr in range(out_r):
         for tfc in range(out_c):
             _matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
-            out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
-                                                   cols, r, c, mode_c)
+            if order == 0:
+                out[tfr, tfc] = nearest_neighbour(<double*>img.data, rows,
+                                                  cols, r, c, mode_c)
+            elif order == 1:
+                out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
+                                                       cols, r, c, mode_c)
 
     return out
diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py
index 4ecb6ba4..75210b54 100644
--- a/skimage/transform/setup.py
+++ b/skimage/transform/setup.py
@@ -20,7 +20,7 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
 
     config.add_extension('_project', sources=['_project.c'],
-                         include_dirs=[get_numpy_include_dirs()])
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config
 
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index ac241722..7be3151b 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -54,20 +54,23 @@ def test_fast_homography():
     H[:2, :2] = [[C, -S], [S, C]]
     H[:2, 2] = [tx, ty]
 
-    for mode in ('constant', 'mirror', 'wrap'):
-        p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1)
-        p1 = fast_homography(img, H, mode=mode)
+    tform = ProjectiveTransform(H)
 
-        # import matplotlib.pyplot as plt
-        # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
-        # ax0.imshow(img)
-        # ax1.imshow(p0, cmap=plt.cm.gray)
-        # ax2.imshow(p1, cmap=plt.cm.gray)
-        # ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
-        # plt.show()
+    for order in range(2):
+        for mode in ('constant', 'mirror', 'wrap'):
+            p0 = warp(img, tform.inverse, mode=mode, order=order)
+            p1 = fast_homography(img, H, mode=mode, order=order)
 
-        d = np.mean(np.abs(p0 - p1))
-        assert d < 0.001
+            # import matplotlib.pyplot as plt
+            # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
+            # ax0.imshow(img)
+            # ax1.imshow(p0, cmap=plt.cm.gray)
+            # ax2.imshow(p1, cmap=plt.cm.gray)
+            # ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+            # plt.show()
+
+            d = np.mean(np.abs(p0 - p1))
+            assert d < 0.001
 
 
 def test_swirl():