def canny(image, mask, sigma, low_threshold, high_threshold):
    '''Edge filter an image using the Canny algorithm.

    sigma - the standard deviation of the Gaussian used
    low_threshold - threshold for edges that connect to high-threshold
                    edges
    high_threshold - threshold of a high-threshold edge

    Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
    Pattern Analysis and Machine Intelligence, 8:679-714, 1986

    William Green's Canny tutorial
    http://www.pages.drexel.edu/~weg22/can_tut.html
    '''
    #
    # The steps involved:
    #
    # * Smooth using the Gaussian with sigma above.
    #
    # * Apply the horizontal and vertical Sobel operators to get the gradients
    #   within the image. The edge strength is the sum of the magnitudes
    #   of the gradients in each direction.
    #
    # * Find the normal to the edge at each point using the arctangent of the
    #   ratio of the Y sobel over the X sobel - pragmatically, we can
    #   look at the signs of X and Y and the relative magnitude of X vs Y
    #   to sort the points into 4 categories: horizontal, vertical,
    #   diagonal and antidiagonal.
    #
    # * Look in the normal and reverse directions to see if the values
    #   in either of those directions are greater than the point in question.
    #   Use interpolation to get a mix of points instead of picking the one
    #   that's the closest to the normal.
    #
    # * Label all points above the high threshold as edges.
    # * Recursively label any point above the low threshold that is 8-connected
    #   to a labeled point as an edge.
    #
    # Regarding masks, any point touching a masked point will have a gradient
    # that is "infected" by the masked point, so it's enough to erode the
    # mask by one and then mask the output. We also mask out the border points
    # because who knows what lies beyond the edge of the image?
    #
    fsmooth = lambda x: gaussian_filter(x, sigma, mode='constant')
    smoothed = smooth_with_function_and_mask(image, fsmooth, mask)
    jsobel = convolve(smoothed, [[-1,0,1],[-2,0,2],[-1,0,1]])
    isobel = convolve(smoothed, [[-1,-2,-1],[0,0,0],[1,2,1]])
    abs_isobel = np.abs(isobel)
    abs_jsobel = np.abs(jsobel)
    magnitude = np.sqrt(isobel*isobel + jsobel*jsobel)
    #
    # Make the eroded mask. Setting the border value to zero will wipe
    # out the image edges for us.
    #
    s = generate_binary_structure(2,2)
    emask = binary_erosion(mask, s, border_value = 0)
    emask = np.logical_and(emask, magnitude > 0)
    #
    #--------- Find local maxima --------------
    #
    # Assign each point to have a normal of 0-45 degrees, 45-90 degrees,
    # 90-135 degrees and 135-180 degrees.
    #
    local_maxima = np.zeros(image.shape,bool)
    #----- 0 to 45 degrees ------
    pts_plus = np.logical_and(isobel >= 0,
                              np.logical_and(jsobel >= 0,
                                             abs_isobel >= abs_jsobel))
    pts_minus = np.logical_and(isobel <= 0,
                               np.logical_and(jsobel <= 0,
                                              abs_isobel >= abs_jsobel))
    pts = np.logical_or(pts_plus, pts_minus)
    pts = np.logical_and(emask, pts)
    # Get the magnitudes shifted left to make a matrix of the points to the
    # right of pts. Similarly, shift left and down to get the points to the
    # top right of pts.
    c1 = magnitude[1:,:][pts[:-1,:]]
    c2 = magnitude[1:,1:][pts[:-1,:-1]]
    m  = magnitude[pts]
    w  = abs_jsobel[pts] / abs_isobel[pts]
    c_plus  = c2 * w + c1 * (1-w) <= m
    c1 = magnitude[:-1,:][pts[1:,:]]
    c2 = magnitude[:-1,:-1][pts[1:,1:]]
    c_minus =  c2 * w + c1 * (1-w) <= m
    local_maxima[pts] = np.logical_and(c_plus, c_minus)
    #----- 45 to 90 degrees ------
    # Mix diagonal and vertical
    #
    pts_plus = np.logical_and(isobel >= 0,
                              np.logical_and(jsobel >= 0,
                                             abs_isobel <= abs_jsobel))
    pts_minus = np.logical_and(isobel <= 0,
                               np.logical_and(jsobel <= 0,
                                              abs_isobel <= abs_jsobel))
    pts = np.logical_or(pts_plus, pts_minus)
    pts = np.logical_and(emask, pts)
    c1 = magnitude[:,1:][pts[:,:-1]]
    c2 = magnitude[1:,1:][pts[:-1,:-1]]
    m  = magnitude[pts]
    w  = abs_isobel[pts] / abs_jsobel[pts]
    c_plus  = c2 * w + c1 * (1-w) <= m
    c1 = magnitude[:,:-1][pts[:,1:]]
    c2 = magnitude[:-1,:-1][pts[1:,1:]]
    c_minus =  c2 * w + c1 * (1-w) <= m
    local_maxima[pts] = np.logical_and(c_plus, c_minus)
    #----- 90 to 135 degrees ------
    # Mix anti-diagonal and vertical
    #
    pts_plus = np.logical_and(isobel <= 0,
                              np.logical_and(jsobel >= 0,
                                             abs_isobel <= abs_jsobel))
    pts_minus = np.logical_and(isobel >= 0,
                               np.logical_and(jsobel <= 0,
                                              abs_isobel <= abs_jsobel))
    pts = np.logical_or(pts_plus, pts_minus)
    pts = np.logical_and(emask, pts)
    c1a = magnitude[:,1:][pts[:,:-1]]
    c2a = magnitude[:-1,1:][pts[1:,:-1]]
    m  = magnitude[pts]
    w  = abs_isobel[pts] / abs_jsobel[pts]
    c_plus  = c2a * w + c1a * (1.0-w) <= m
    c1 = magnitude[:,:-1][pts[:,1:]]
    c2 = magnitude[1:,:-1][pts[:-1,1:]]
    c_minus =  c2 * w + c1 * (1.0-w) <= m
    cc = np.logical_and(c_plus,c_minus)
    local_maxima[pts] = np.logical_and(c_plus, c_minus)
    #----- 135 to 180 degrees ------
    # Mix anti-diagonal and anti-horizontal
    #
    pts_plus = np.logical_and(isobel <= 0,
                              np.logical_and(jsobel >= 0,
                                             abs_isobel >= abs_jsobel))
    pts_minus = np.logical_and(isobel >= 0,
                               np.logical_and(jsobel <= 0,
                                              abs_isobel >= abs_jsobel))
    pts = np.logical_or(pts_plus, pts_minus)
    pts = np.logical_and(emask, pts)
    c1 = magnitude[:-1,:][pts[1:,:]]
    c2 = magnitude[:-1,1:][pts[1:,:-1]]
    m  = magnitude[pts]
    w  = abs_jsobel[pts] / abs_isobel[pts]
    c_plus  = c2 * w + c1 * (1-w) <= m
    c1 = magnitude[1:,:][pts[:-1,:]]
    c2 = magnitude[1:,:-1][pts[:-1,1:]]
    c_minus =  c2 * w + c1 * (1-w) <= m
    local_maxima[pts] = np.logical_and(c_plus, c_minus)
    #
    #---- Create two masks at the two thresholds.
    #
    high_mask = np.logical_and(local_maxima, magnitude >= high_threshold)
    low_mask  = np.logical_and(local_maxima, magnitude >= low_threshold)
    #
    # Segment the low-mask, then only keep low-segments that have
    # some high_mask component in them
    #
    labels,count = label(low_mask, np.ones((3,3),bool))
    if count == 0:
        return low_mask

    sums = fix(scind.sum(high_mask, labels, np.arange(count,dtype=np.int32)+1))
    good_label = np.zeros((count+1,),bool)
    good_label[1:] = sums > 0
    output_mask = good_label[labels]
    return output_mask