From 6905ed25b40991f88214c1e0691e36c6292ddacd Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Sat, 26 Mar 2011 07:39:56 -0700 Subject: [PATCH 1/7] 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) From 22f983ce17c881866294bed3154106f8f0f134c4 Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Wed, 30 Mar 2011 14:20:21 -0700 Subject: [PATCH 2/7] Add init entry so importing works right --- scikits/image/filter/__init__.py | 1 + 1 file changed, 1 insertion(+) 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 From 2961147e2767a7f48b78eecad4bd0a91df4b0702 Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Thu, 31 Mar 2011 22:42:13 -0700 Subject: [PATCH 3/7] Format for pep8 compliance, improve documentation, and make mask optional. --- scikits/image/filter/canny.py | 129 +++++++++++++---------- scikits/image/filter/tests/test_canny.py | 10 +- 2 files changed, 78 insertions(+), 61 deletions(-) diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py index 1cf3e63f..aeaba27d 100644 --- a/scikits/image/filter/canny.py +++ b/scikits/image/filter/canny.py @@ -1,6 +1,6 @@ '''canny.py - Canny Edge detector -Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans. +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. @@ -14,26 +14,39 @@ 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]) +import scipy.ndimage as ndi +from scipy.ndimage import (gaussian_filter, convolve, + generate_binary_structure, binary_erosion, label) -def canny(image, mask, sigma, low_threshold, high_threshold): +def canny(image, sigma, low_threshold, high_threshold, mask=None): '''Edge filter an image using the Canny algorithm. + + Parameters + ----------- + image : array_like + The input image to detect edges on. - 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 + sigma : float + The standard deviation of the Gaussian filter - Canny, J., A Computational Approach To Edge Detection, IEEE Trans. + low_threshold : float + The lower bound for hysterisis thresholding (linking edges) + + high_threshold : float + The upper bound for hysterisis thresholding (linking edges) + + mask : array of booleans, 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 @@ -68,21 +81,23 @@ def canny(image, mask, sigma, low_threshold, high_threshold): # 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 = convolve(smoothed, [[-1,0,1],[-2,0,2],[-1,0,1]]) + 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) + 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) + eroded_mask = binary_erosion(mask, s, border_value=0) + eroded_mask = np.logical_and(eroded_mask, magnitude > 0) # #--------- Find local maxima -------------- # @@ -91,102 +106,104 @@ def canny(image, mask, sigma, low_threshold, high_threshold): # local_maxima = np.zeros(image.shape,bool) #----- 0 to 45 degrees ------ - pts_plus = np.logical_and(isobel >= 0, - np.logical_and(jsobel >= 0, + 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) + pts = np.logical_and(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 + 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 + 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, + 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, + np.logical_and(jsobel <= 0, abs_isobel <= abs_jsobel)) pts = np.logical_or(pts_plus, pts_minus) - pts = np.logical_and(emask, pts) + pts = np.logical_and(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 + 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 + 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, + 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, + np.logical_and(jsobel <= 0, abs_isobel <= abs_jsobel)) pts = np.logical_or(pts_plus, pts_minus) - pts = np.logical_and(emask, pts) + pts = np.logical_and(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 + 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 + 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, + 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, + np.logical_and(jsobel <= 0, abs_isobel >= abs_jsobel)) pts = np.logical_or(pts_plus, pts_minus) - pts = np.logical_and(emask, pts) + pts = np.logical_and(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 + 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 + 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) + 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 + # 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) + 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 + return output_mask diff --git a/scikits/image/filter/tests/test_canny.py b/scikits/image/filter/tests/test_canny.py index 00eaecc0..a92a576e 100644 --- a/scikits/image/filter/tests/test_canny.py +++ b/scikits/image/filter/tests/test_canny.py @@ -6,20 +6,20 @@ 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) + 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)),np.zeros((20,20),bool), - 4,0,0) + 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),np.ones(c.shape,bool), 4, 0, 0) + 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 @@ -44,7 +44,7 @@ class TestCanny(unittest.TestCase): 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) + 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 From 58b632a536b5bf072b461e60dd744d499c773145 Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Fri, 1 Apr 2011 21:35:21 -0700 Subject: [PATCH 4/7] Remove smooth.py for now --- scikits/image/filter/canny.py | 33 +++++++++++- scikits/image/filter/smooth.py | 97 ---------------------------------- 2 files changed, 32 insertions(+), 98 deletions(-) delete mode 100644 scikits/image/filter/smooth.py diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py index aeaba27d..58d179f9 100644 --- a/scikits/image/filter/canny.py +++ b/scikits/image/filter/canny.py @@ -13,12 +13,43 @@ Original author: Lee Kamentsky ''' import numpy as np -from smooth import smooth_with_function_and_mask 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. diff --git a/scikits/image/filter/smooth.py b/scikits/image/filter/smooth.py deleted file mode 100644 index dad706af..00000000 --- a/scikits/image/filter/smooth.py +++ /dev/null @@ -1,97 +0,0 @@ -"""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 From 1ce0936e3db90c3b83e63e7c1726a7592f6749b2 Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Mon, 4 Apr 2011 21:33:44 -0700 Subject: [PATCH 5/7] Update docs and convert convolution to ndimage sobel function --- scikits/image/filter/canny.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py index 58d179f9..51d4f200 100644 --- a/scikits/image/filter/canny.py +++ b/scikits/image/filter/canny.py @@ -55,8 +55,9 @@ def canny(image, sigma, low_threshold, high_threshold, mask=None): Parameters ----------- - image : array_like - The input image to detect edges on. + 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 @@ -67,7 +68,7 @@ def canny(image, sigma, low_threshold, high_threshold, mask=None): high_threshold : float The upper bound for hysterisis thresholding (linking edges) - mask : array of booleans, optional + mask : array, dtype=bool, optional An optional mask to limit the application of Canny to a certain area. Returns @@ -116,8 +117,8 @@ def canny(image, sigma, low_threshold, high_threshold, mask=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 = convolve(smoothed, [[-1,0,1], [-2,0,2], [-1,0,1]]) - isobel = convolve(smoothed, [[-1,-2,-1],[0,0,0],[1,2,1]]) + jsobel = ndi.sobel(smoothed, axis=1) + isobel = ndi.sobel(smoothed, axis=0) abs_isobel = np.abs(isobel) abs_jsobel = np.abs(jsobel) magnitude = np.sqrt(isobel * isobel + jsobel * jsobel) From e60d38892e3fc4e524a9dc02bf0a2fa39501f29c Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Mon, 4 Apr 2011 22:48:40 -0700 Subject: [PATCH 6/7] Converted logical ops to bitwise ops and made cosmetic changes wrt pep8 compliance --- scikits/image/filter/canny.py | 98 +++++++++++++++-------------------- 1 file changed, 41 insertions(+), 57 deletions(-) diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py index 51d4f200..39c71b5b 100644 --- a/scikits/image/filter/canny.py +++ b/scikits/image/filter/canny.py @@ -121,15 +121,15 @@ def canny(image, sigma, low_threshold, high_threshold, mask=None): isobel = ndi.sobel(smoothed, axis=0) abs_isobel = np.abs(isobel) abs_jsobel = np.abs(jsobel) - magnitude = np.sqrt(isobel * isobel + jsobel * 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) + s = generate_binary_structure(2, 2) eroded_mask = binary_erosion(mask, s, border_value=0) - eroded_mask = np.logical_and(eroded_mask, magnitude > 0) + eroded_mask = eroded_mask & (magnitude > 0) # #--------- Find local maxima -------------- # @@ -138,97 +138,81 @@ def canny(image, sigma, low_threshold, high_threshold, mask=None): # 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(eroded_mask, pts) + 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]] + 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:]] + 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) + local_maxima[pts] = 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(eroded_mask, pts) - c1 = magnitude[:,1:][pts[:,:-1]] - c2 = magnitude[1:,1:][pts[:-1,:-1]] + 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:]] + 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) + local_maxima[pts] = 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(eroded_mask, pts) - c1a = magnitude[:,1:][pts[:,:-1]] - c2a = magnitude[:-1,1:][pts[1:,:-1]] + 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:]] + 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) + local_maxima[pts] = 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(eroded_mask, pts) - c1 = magnitude[:-1,:][pts[1:,:]] - c2 = magnitude[:-1,1:][pts[1:,:-1]] + 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,:]] + 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) + local_maxima[pts] = 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) + 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)) + labels,count = label(low_mask, np.ndarray((3, 3),bool)) if count == 0: return low_mask From 5264e9d1d473047bf3e751e1200e1609740e4e8f Mon Sep 17 00:00:00 2001 From: Dan Farmer Date: Wed, 6 Apr 2011 22:11:17 -0700 Subject: [PATCH 7/7] Fix docstrings and one more pep8 issue --- scikits/image/filter/canny.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/scikits/image/filter/canny.py b/scikits/image/filter/canny.py index 39c71b5b..49db7343 100644 --- a/scikits/image/filter/canny.py +++ b/scikits/image/filter/canny.py @@ -24,13 +24,13 @@ def smooth_with_function_and_mask(image, function, mask): Parameters ---------- image : array - The image to smooth + The image to smooth function : callable - A function that takes an image and returns a smoothed image + A function that takes an image and returns a smoothed image mask : array - Mask with 1's for significant pixels, 0 for masked pixels + Mask with 1's for significant pixels, 0 for masked pixels Notes ------ @@ -41,12 +41,12 @@ def smooth_with_function_and_mask(image, function, mask): 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) + 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 @@ -56,25 +56,25 @@ def canny(image, sigma, low_threshold, high_threshold, mask=None): Parameters ----------- image : array_like, dtype=float - The greyscale input image to detect edges on; should be normalized to 0.0 - to 1.0. + 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 + The standard deviation of the Gaussian filter low_threshold : float - The lower bound for hysterisis thresholding (linking edges) + The lower bound for hysterisis thresholding (linking edges) high_threshold : float - The upper bound for hysterisis thresholding (linking edges) + 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. + An optional mask to limit the application of Canny to a certain area. Returns ------- output : array (image) - The binary edge map. + The binary edge map. References -----------