[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