From 6905ed25b40991f88214c1e0691e36c6292ddacd Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Sat, 26 Mar 2011 07:39:56 -0700 Subject: [PATCH] Merge Cell Profiler Canny edge detector --- scikits/image/filter/canny.py | 192 +++++++++++++++++++++++ scikits/image/filter/smooth.py | 97 ++++++++++++ scikits/image/filter/tests/test_canny.py | 58 +++++++ 3 files changed, 347 insertions(+) create mode 100644 scikits/image/filter/canny.py create mode 100644 scikits/image/filter/smooth.py create mode 100644 scikits/image/filter/tests/test_canny.py diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py new file mode 100644 index 00000000..1cf3e63f --- /dev/null +++ b/scikits/image/filter/canny.py @@ -0,0 +1,192 @@ +'''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 +from smooth import smooth_with_function_and_mask +import scipy.ndimage as scind +from scipy.ndimage import gaussian_filter, convolve, generate_binary_structure, \ + binary_erosion, label + +def fix(whatever_it_returned): + if getattr(whatever_it_returned,"__getitem__",False): + return np.array(whatever_it_returned) + else: + return np.array([whatever_it_returned]) + + +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.ndarray((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 diff --git a/scikits/image/filter/smooth.py b/scikits/image/filter/smooth.py new file mode 100644 index 00000000..dad706af --- /dev/null +++ b/scikits/image/filter/smooth.py @@ -0,0 +1,97 @@ +"""smooth.py - smoothing of images + +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.linalg + +def smooth_with_noise(image, bits): + """Smooth the image with a per-pixel random multiplier + + image - the image to perturb + bits - the noise is this many bits below the pixel value + + The noise is random with normal distribution, so the individual pixels + get either multiplied or divided by a normally distributed # of bits + """ + + rr = np.random.RandomState() + rr.seed(0) + r = rr.normal(size=image.shape) + delta = pow(2.0,-bits) + image_copy = np.clip(image, delta, 1) + result = np.exp2(np.log2(image_copy + delta) * r + + (1-r) * np.log2(image_copy)) + result[result>1] = 1 + result[result<0] = 0 + return result + +def smooth_with_function_and_mask(image, function, mask): + """Smooth an image with a linear function, ignoring the contribution of masked pixels + + image - image to smooth + function - a function that takes an image and returns a smoothed image + mask - mask with 1's for significant pixels, 0 for masked pixels + + 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 circular_gaussian_kernel(sd,radius): + """Create a 2-d Gaussian convolution kernel + + sd - standard deviation of the gaussian in pixels + radius - build a circular kernel that convolves all points in the circle + bounded by this radius + """ + i,j = np.mgrid[-radius:radius+1,-radius:radius+1].astype(float) / radius + mask = i**2 + j**2 <= 1 + i = i * radius / sd + j = j * radius / sd + + kernel = np.zeros((2*radius+1,2*radius+1)) + kernel[mask] = np.e ** (-(i[mask]**2+j[mask]**2) / + (2 * sd **2)) + # + # Normalize the kernel so that there is no net effect on a uniform image + # + kernel = kernel / np.sum(kernel) + return kernel + +def fit_polynomial(pixel_data, mask): + '''Return an "image" which is a polynomial fit to the pixel data + + Fit the image to the polynomial Ax**2+By**2+Cxy+Dx+Ey+F + ''' + mask = np.logical_and(mask,pixel_data > 0) + if not np.any(mask): + return pixel_data + x,y = np.mgrid[0:pixel_data.shape[0],0:pixel_data.shape[1]] + x2 = x*x + y2 = y*y + xy = x*y + o = np.ones(pixel_data.shape) + a = np.array([x[mask],y[mask],x2[mask],y2[mask],xy[mask],o[mask]]) + coeffs = scipy.linalg.lstsq(a.transpose(),pixel_data[mask])[0] + output_pixels = np.sum([coeff * index for coeff, index in + zip(coeffs, [x,y,x2,y2,xy,o])],0) + output_pixels[output_pixels > 1] = 1 + output_pixels[output_pixels < 0] = 0 + return output_pixels diff --git a/scikits/image/filter/tests/test_canny.py b/scikits/image/filter/tests/test_canny.py new file mode 100644 index 00000000..00eaecc0 --- /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)),np.ones((20,20),bool), 4, 0, 0) + 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)),np.zeros((20,20),bool), + 4,0,0) + 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),np.ones(c.shape,bool), 4, 0, 0) + # + # 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,np.ones(c.shape,bool), 4, .1, .2) + # + # 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)