# [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).

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):
# 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
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)) )
# Now, hysteresis thresholding: find seeds above a high threshold,
then
# expand out until we go below the low threshold

```