[SciPy-User] Gaussian Filter

[Sturla Molden]

Brian Thorne skrev:
> But I am curious as to what causes the edges to do that in the IIR 
> filter version?
>
> I made a pretty (IMHO) plot of the diff image seperating each channel:
> http://1.bp.blogspot.com/_lewp47C9PZI/SqJSUG7MJEI/AAAAAAAAAgs/s8aSE0FBCwI/s1600-h/lena_diff_ndfilt_iir.png

In your code you save to JPEG (lossy compression) before taking the 
difference. Don't do that, you incur some error from the compression.

Mind you that neither filter are 'correct'. The Gaussian filters in 
ndimage and OpenCV truncate the Gaussian, my version use a recursive 
approximation. The Gaussian filter in ndimage and OpenCV have larger 
truncation errors near the edges: As the kernel overlaps with an edge, 
the kernel is truncated further and truncation error increase. Thus, you 
were looking at the discreapancy between two approximations. You cannot 
attribute the edge difference to error in the IIR approximation from 
this result.

First, it seems it helps to pad with some zeros (I forgot to do that, 
sorry!), e.g. setting padding=3*sigma.

def gaussian2D(image, sigma, padding=None):
    if padding is None: padding = 3*sigma
    n,m = image.shape[0],image.shape[1]
    tmp = zeros((n + 2*padding, m + 2*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[padding:n+padding,padding:m+padding] = 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[padding:n+padding,padding:m+padding]

Now I get a maximum illumination difference of 9 (range 0-255) between 
ndimage and IIR. Another image (Forest.jpg, sample image in Windows 
Vista) gave a maximum of 5. The biggest discrepancy was near the edges 
here as well.


There is still discrepancies at the edges. Who are the culprit?

We can easily find 'facit' using an FFT. By truncating at 3 standard 
deviations, the truncation error from FFT convolution will be very small:

def gaussian2D_FFT(image, sigma, padding=None):
    if padding is None: padding = 3*sigma
    n,m = image.shape[0],image.shape[1]
    tmp = zeros((n + 2*padding, m + 2*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'
    tmp[:n,:m] = image
    x,y = meshgrid(range(tmp.shape[1]),range(tmp.shape[0]))
    B = exp(-(((x-padding)**2) + ((y-padding)**2))/(2.*sigma*sigma))
    B[:] *= 1./B.sum()
    retv = irfft2(rfft2(B)*rfft2(tmp))
    return retv[padding:n+padding,padding:m+padding]

This FFT-based filter also have loss of gain near the edges due to 
padding, so we adjust it accordingly:

img = np.fromstring(Image.open('lenna.png').tostring(), dtype=np.uint8)\
            .reshape((512,512,3)).astype(float)
gain = gaussian2D_FFT(np.ones(img.shape[:2]), 10)
for i in range(3):
    img[:,:,i] = gaussian2D_FFT(img[:,:,i], 10) / (gain + 1E-100)

Now, comparing what we get from scipy.ndimage.filters.gaussian_filter 
and IIR with the FFT 'facit' (very small truncation error) :

* IIR vs. FFT:     maximum illumination difference: 3
* ndimage vs. FFT: maximum illumination difference: 9 (close to edges)
* ndimage vs. IIR: maximum illumination difference: 9 (close to edges)

Which in fact means that the IIR approximation is more accurate than 
OpenCV and ndimage's Gaussian filter! The edge difference in fact comes 
from truncation error in ndimage and OpenVC, not from the IIR being 
inaccurate. It is due to increased truncation error near the edges.

I have prepared similar images like yours:

IIR vs. FFT: http://hostdump.com/images/iirvsfft.png

ndimage vs. FFT: http://hostdump.com/images/ndarrayvsf.png

ndimage vs. IIR: http://hostdump.com/images/ndarrayvsi.png

Thus, the IIR is the better approximation to the Gaussian filter.



Regards,
Sturla Molden