[SciPy-user] Some more useful image filters: anisotropic diffusion and Canny edge-finding
Zachary Pincus
zachary.pincus@yale....
Fri Feb 27 15:47:55 CST 2009
Hi all,
As I've been quite happily using with Nadav's bilateral filtering
code, I figured I should reciprocate and post some other useful basic
image-processing algorithms that I've recently numpy-ified. Here are
first-passes at Canny edge finding and Perona and Malik-style
anisotropic diffusion; both are pretty simplistic and not particularly
fast, but they work as advertised.
Zach
import numpy
import scipy.ndimage as ndimage
# Anisotropic Diffusion, as per Perona and Malik's paper (see section
V).
def _exp(image_gradient, scale):
return numpy.exp(-(numpy.absolute(image_gradient)/scale)**2)
def _inv(image_gradient, scale):
return 1 / (1 + (numpy.absolute(image_gradient)/scale)**2)
def anisotropic_diffusion(image, num_iters=10, scale=10,
step_size=0.2, conduction_function=_inv):
# 'step_size' is Perona and Malik's lambda parameter; scale is
their 'K' parameter.
# The 'conduction_function' is the function 'g' in the original
formulation;
# if this function simply returns a constant, the result is
Gaussian blurring.
if step_size > 0.25:
raise ValueError('step_size parameter must be <= 0.25 for
numerical stability.')
image = image.copy()
# simplistic boundary conditions -- no diffusion at the boundary
central = image[1:-1, 1:-1]
n = image[:-2, 1:-1]
s = image[2:, 1:-1]
e = image[1:-1, :-2]
w = image[1:-1, 2:]
directions = [s,e,w]
for i in range(num_iters):
di = n - central
accumulator = conduction_function(di, scale)*di
for direction in directions:
di = direction - central
accumulator += conduction_function(di, scale)*di
accumulator *= step_size
central += accumulator
return image
# Canny edge-finding, implemented as per the Wikipedia article
# Note that this takes four passes through the image to do the
# non-maximal suppression, whereas a c or cython loop could do
# it in one.
# Filter kernels for calculating the value of neighbors in several
directions
_N = numpy.array([[0, 1, 0],
[0, 0, 0],
[0, 1, 0]], dtype=bool)
_NE = numpy.array([[0, 0, 1],
[0, 0, 0],
[1, 0, 0]], dtype=bool)
_W = numpy.array([[0, 0, 0],
[1, 0, 1],
[0, 0, 0]], dtype=bool)
_NW = numpy.array([[1, 0, 0],
[0, 0, 0],
[0, 0, 1]], dtype=bool)
# After quantizing the angles, vertical (north-south) edges get values
of 3,
# northwest-southeast edges get values of 2, and so on, as below:
_NE_d = 0
_W_d = 1
_NW_d = 2
_N_d = 3
def canny(image, high_threshold, low_threshold):
grad_x = ndimage.sobel(image, 0)
grad_y = ndimage.sobel(image, 1)
grad_mag = numpy.sqrt(grad_x**2+grad_y**2)
grad_angle = numpy.arctan2(grad_y, grad_x)
# next, scale the angles in the range [0, 3] and then round to
quantize
quantized_angle = numpy.around(3 * (grad_angle + numpy.pi) /
(numpy.pi * 2))
# Non-maximal suppression: an edge pixel is only good if its
magnitude is
# greater than its neighbors normal to the edge direction. We
quantize
# edge direction into four angles, so we only need to look at four
# sets of neighbors
NE = ndimage.maximum_filter(grad_mag, footprint=_NE)
W = ndimage.maximum_filter(grad_mag, footprint=_W)
NW = ndimage.maximum_filter(grad_mag, footprint=_NW)
N = ndimage.maximum_filter(grad_mag, footprint=_N)
thinned = (((g > W) & (quantized_angle == _N_d )) |
((g > N) & (quantized_angle == _W_d )) |
((g > NW) & (quantized_angle == _NE_d)) |
((g > NE) & (quantized_angle == _NW_d)) )
thinned_grad = thinned * grad_mag
# Now, hysteresis thresholding: find seeds above a high threshold,
then
# expand out until we go below the low threshold
high = thinned_grad > high_threshold
low = thinned_grad > low_threshold
canny_edges = ndimage.binary_dilation(high, iterations=-1, mask=low)
return grad_mag, thinned_grad, canny_edges
More information about the SciPy-user
mailing list