mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-04 13:14:23 +08:00
Merge Cell Profiler Canny edge detector
This commit is contained in:
Dan Farmer
@@ -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
|
||||
@@ -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
|
||||
@@ -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)
|
||||
Reference in New Issue
Block a user