mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-16 11:21:25 +08:00
Merge remote branch 'dfarmer/dfarmer-filters-canny'
This commit is contained in:
Stefan van der Walt
@@ -1,2 +1,3 @@
|
||||
from lpi_filter import *
|
||||
from ctmf import median_filter
|
||||
from canny import canny
|
||||
|
||||
@@ -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
|
||||
@@ -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)
|
||||
Reference in New Issue
Block a user