[SciPy-User] Gaussian Filter

[Sturla Molden]

Brian Thorne skrev:
 > Similar question, but now a bit harder. I have this code (pieced 
together from a few files) that does a gaussian filter on a single image 
in both OpenCV and in SciPy.
 > It is now at a point where I cannot tell them apart with a visual 
inspection, but a imshow(image1 - image2) begs to differ. Is it going to 
be possible to get the exact same output?
 >
 >
 > from opencv import cv
 > from opencv import adaptors
 > from __future__ import division
 > import numpy as np
 > from numpy import array, uint8
 > from scipy import signal, ndimage
 >
 > @scipyFromOpenCV
 > def gaussianBlur(np_image):
 >     """Blur an image with scipy"""
 >     sigma = opencvFilt2sigma(43.0)
 >    
 >     result = ndimage.filters.gaussian_filter(np_image,
 >                             sigma=(sigma, sigma, 0),
 >                             order=0,
 >                             mode='reflect'
 >                             )
 >     return result
 >
 > def gaussianBlur(image, filterSize=43, sigma=opencvFilt2sigma(43)):
 >     """Blur an image with a particular strength filter.
 >     Default is 43, 139 gives a very strong blur, but takes a while
 >     """
 >    
 >     # Carry out the filter operation
 >     cv.cvSmooth(image, image, cv.CV_GAUSSIAN, filterSize, 0, sigma)
 >     return image


For Gaussian filtering (and Gaussian blur in particular) one can also 
use a fast IIR approximation. It will not be faster if you use small 
truncated kernels like 3 x 3 pixels, but is quite efficient and run-time 
does not depend on sigma. (I think I've posted this code before, though.)

Regards,
Sturla Molden



from numpy import array, zeros, ones, flipud, fliplr
from scipy.signal import lfilter
from math import sqrt

def __gausscoeff(s):

    # Young, I.T. and van Vliet, L.J. (1995). Recursive implementation
    # of the Gaussian filter, Signal Processing, 44: 139-151.

    if s < .5: raise ValueError, \
        'Sigma for Gaussian filter must be >0.5 samples'
    q = 0.98711*s - 0.96330 if s > 0.5 else 3.97156 \
        - 4.14554*sqrt(1.0 - 0.26891*s)
    b = zeros(4)
    b[0] = 1.5785 + 2.44413*q + 1.4281*q**2 + 0.422205*q**3
    b[1] = 2.44413*q + 2.85619*q**2 + 1.26661*q**3
    b[2] = -(1.4281*q**2 + 1.26661*q**3)
    b[3] = 0.422205*q**3
    B = 1.0 - ((b[1] + b[2] + b[3])/b[0])

    # convert to a format compatible with lfilter's
    # difference equation

    B = array([B])
    A = ones(4)
    A[1:] = -b[1:]/b[0]
    return B,A

def gaussian1D(signal, sigma, padding=0):
    n = signal.shape[0]
    tmp = zeros(n + padding)
    if tmp.shape[0] < 4: raise ValueError, \
        'Signal and padding too short'
    tmp[:n] = signal
    B,A = __gausscoeff(sigma)
    tmp = lfilter(B, A, tmp)
    tmp = tmp[::-1]
    tmp = lfilter(B, A, tmp)
    tmp = tmp[::-1]
    return tmp[:n]

def gaussian2D(image, sigma, padding=0):
    n,m = image.shape[0],image.shape[1]
    tmp = zeros((n + padding, m + padding))
    if tmp.shape[0] < 4: raise ValueError, \
        'Image and padding too small'
    if tmp.shape[1] < 4: raise ValueError, \
        'Image and padding too small'
    B,A = __gausscoeff(sigma)
    tmp[:n,:m] = image
    tmp = lfilter(B, A, tmp, axis=0)
    tmp = flipud(tmp)
    tmp = lfilter(B, A, tmp, axis=0)
    tmp = flipud(tmp)
    tmp = lfilter(B, A, tmp, axis=1)
    tmp = fliplr(tmp)
    tmp = lfilter(B, A, tmp, axis=1)
    tmp = fliplr(tmp)
    return tmp[:n,:m]