diff --git a/scikits/image/filter/__init__.py b/scikits/image/filter/__init__.py index afcc9db5..29247baf 100644 --- a/scikits/image/filter/__init__.py +++ b/scikits/image/filter/__init__.py @@ -1,2 +1,3 @@ from lpi_filter import * from ctmf import median_filter +from canny import canny diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py new file mode 100644 index 00000000..49db7343 --- /dev/null +++ b/scikits/image/filter/canny.py @@ -0,0 +1,225 @@ +'''canny.py - Canny Edge detector + +Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans. + Pattern Analysis and Machine Intelligence, 8:679-714, 1986 + +Originally part of CellProfiler, code licensed under both GPL and BSD licenses. +Website: http://www.cellprofiler.org +Copyright (c) 2003-2009 Massachusetts Institute of Technology +Copyright (c) 2009-2011 Broad Institute +All rights reserved. +Original author: Lee Kamentsky + +''' + +import numpy as np +import scipy.ndimage as ndi +from scipy.ndimage import (gaussian_filter, convolve, + generate_binary_structure, binary_erosion, label) + + +def smooth_with_function_and_mask(image, function, mask): + """Smooth an image with a linear function, ignoring masked pixels + + Parameters + ---------- + image : array + The image to smooth + + function : callable + A function that takes an image and returns a smoothed image + + mask : array + Mask with 1's for significant pixels, 0 for masked pixels + + Notes + ------ + This function calculates the fractional contribution of masked pixels + by applying the function to the mask (which gets you the fraction of + the pixel data that's due to significant points). We then mask the image + and apply the function. The resulting values will be lower by the bleed-over + fraction, so you can recalibrate by dividing by the function on the mask + to recover the effect of smoothing from just the significant pixels. + """ + not_mask = np.logical_not(mask) + bleed_over = function(mask.astype(float)) + masked_image = np.zeros(image.shape, image.dtype) + masked_image[mask] = image[mask] + smoothed_image = function(masked_image) + output_image = smoothed_image / (bleed_over + np.finfo(float).eps) + return output_image + + +def canny(image, sigma, low_threshold, high_threshold, mask=None): + '''Edge filter an image using the Canny algorithm. + + Parameters + ----------- + image : array_like, dtype=float + The greyscale input image to detect edges on; should be normalized to 0.0 + to 1.0. + + sigma : float + The standard deviation of the Gaussian filter + + low_threshold : float + The lower bound for hysterisis thresholding (linking edges) + + high_threshold : float + The upper bound for hysterisis thresholding (linking edges) + + mask : array, dtype=bool, optional + An optional mask to limit the application of Canny to a certain area. + + Returns + ------- + output : array (image) + The binary edge map. + + References + ----------- + 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? + # + if mask is None: + mask = np.ones(image.shape, dtype=bool) + fsmooth = lambda x: gaussian_filter(x, sigma, mode='constant') + smoothed = smooth_with_function_and_mask(image, fsmooth, mask) + jsobel = ndi.sobel(smoothed, axis=1) + isobel = ndi.sobel(smoothed, axis=0) + abs_isobel = np.abs(isobel) + abs_jsobel = np.abs(jsobel) + magnitude = np.hypot(isobel, 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) + eroded_mask = binary_erosion(mask, s, border_value=0) + eroded_mask = eroded_mask & (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 = (isobel >= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel) + pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel) + pts = pts_plus | pts_minus + pts = eroded_mask & 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] = c_plus & c_minus + #----- 45 to 90 degrees ------ + # Mix diagonal and vertical + # + pts_plus = (isobel >= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel) + pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel) + pts = pts_plus | pts_minus + pts = eroded_mask & 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] = c_plus & c_minus + #----- 90 to 135 degrees ------ + # Mix anti-diagonal and vertical + # + pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel) + pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel) + pts = pts_plus | pts_minus + pts = eroded_mask & 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] = c_plus & c_minus + #----- 135 to 180 degrees ------ + # Mix anti-diagonal and anti-horizontal + # + pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel) + pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel) + pts = pts_plus | pts_minus + pts = eroded_mask & 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] = c_plus & c_minus + # + #---- Create two masks at the two thresholds. + # + high_mask = local_maxima & (magnitude >= high_threshold) + low_mask = 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.ndarray((3, 3),bool)) + if count == 0: + return low_mask + + sums = (np.array(ndi.sum(high_mask,labels, + np.arange(count,dtype=np.int32) + 1), + copy=False, ndmin=1)) + good_label = np.zeros((count + 1,),bool) + good_label[1:] = sums > 0 + output_mask = good_label[labels] + return output_mask diff --git a/scikits/image/filter/tests/test_canny.py b/scikits/image/filter/tests/test_canny.py new file mode 100644 index 00000000..a92a576e --- /dev/null +++ b/scikits/image/filter/tests/test_canny.py @@ -0,0 +1,58 @@ +import unittest +import numpy as np +from scipy.ndimage import binary_dilation, binary_erosion +import scikits.image.filter as F + +class TestCanny(unittest.TestCase): + def test_00_00_zeros(self): + '''Test that the Canny filter finds no points for a blank field''' + result = F.canny(np.zeros((20,20)), 4, 0, 0, np.ones((20,20),bool)) + self.assertFalse(np.any(result)) + + def test_00_01_zeros_mask(self): + '''Test that the Canny filter finds no points in a masked image''' + result = (F.canny(np.random.uniform(size=(20,20)), 4,0,0, + np.zeros((20,20),bool))) + self.assertFalse(np.any(result)) + + def test_01_01_circle(self): + '''Test that the Canny filter finds the outlines of a circle''' + i,j = np.mgrid[-200:200,-200:200].astype(float) / 200 + c = np.abs(np.sqrt(i*i+j*j) - .5) < .02 + result = F.canny(c.astype(float), 4, 0, 0,np.ones(c.shape,bool)) + # + # erode and dilate the circle to get rings that should contain the + # outlines + # + cd = binary_dilation(c, iterations=3) + ce = binary_erosion(c,iterations=3) + cde = np.logical_and(cd, np.logical_not(ce)) + self.assertTrue(np.all(cde[result])) + # + # The circle has a radius of 100. There are two rings here, one + # for the inside edge and one for the outside. So that's 100 * 2 * 2 * 3 + # for those places where pi is still 3. The edge contains both pixels + # if there's a tie, so we bump the count a little. + # + point_count = np.sum(result) + self.assertTrue(point_count > 1200) + self.assertTrue(point_count < 1600) + + def test_01_02_circle_with_noise(self): + '''Test that the Canny filter finds the circle outlines in a noisy image''' + np.random.seed(0) + i,j = np.mgrid[-200:200,-200:200].astype(float) / 200 + c = np.abs(np.sqrt(i*i+j*j) - .5) < .02 + cf = c.astype(float) * .5 + np.random.uniform(size=c.shape)*.5 + result = F.canny(cf, 4, .1, .2,np.ones(c.shape,bool)) + # + # erode and dilate the circle to get rings that should contain the + # outlines + # + cd = binary_dilation(c, iterations=4) + ce = binary_erosion(c,iterations=4) + cde = np.logical_and(cd, np.logical_not(ce)) + self.assertTrue(np.all(cde[result])) + point_count = np.sum(result) + self.assertTrue(point_count > 1200) + self.assertTrue(point_count < 1600)