from scipy.ndimage import filters
def ims(image,**kw):
    size = kw.get("s",8)
    if "s" in kw: del kw["s"]
    subplots(1,1,figsize=(size,size))
    gray(); imshow(image,origin='lower',**kw)
def imp(image,**kw):
    subplots(1,1,figsize=(6,6))
    gray(); imshow(image,origin='lower',interpolation='nearest',**kw)
def imrow(*args,**kw):
    gray()
    size = kw.get("s",8)
    if "s" in kw: del kw["s"]
    n = len(args)
    subplots(1,n,figsize=(n*size,size))
    for i,im in enumerate(args):
        subplot(1,n,i+1); imshow(im,origin='lower',**kw)
def rescale(image):
    return (image-amin(image))/max(1e-4,amax(image)-amin(image))

# test image examples
image = mean(imread("bridge.jpg")/255.0,axis=2)[::-1]
test = mean(imread("testimage.jpg")/255.0,axis=2)[::-1]
imrow(image,test)

# 2D array representing an impulse
impulse = zeros((64,64))
impulse[32,32] = 1.0
imp(impulse)

# filtering with a box filter, aka uniform filter
from scipy.ndimage.filters import uniform_filter as bf
imrow(bf(image,10.0),bf(test,10.0),bf(impulse,10.0))

# gaussian filter
from scipy.ndimage.filters import gaussian_filter as gf
imrow(gf(image,5.0),gf(test,5.0),gf(impulse,5.0))

# reflect, wrap, constant boundary conditions
from scipy.ndimage.filters import uniform_filter as bf
imrow(bf(image,50.0,mode='reflect'),
      bf(image,50.0,mode='wrap'),
      rescale(bf(image,50.0,mode='constant',cval=0.0)))

# noise reduction by blurring
noise = imread("imagenoise.jpg")
imrow(noise,gf(noise,(5.0,5.0,0.0)))

# stack of noisy images
noisy = array([image+0.1*randn(*image.shape) 
                  for i in range(50)])
ims(noisy[10])

# noise reduction by averaging a stack of noisy images
avg = mean(noisy,axis=0)
imrow(image,clip(avg,0.0,1.0))

# demonstration of linearity
def F(image): return gf(image,5.0)
result1 = F(image)+2.0*F(test)
result2 = F(image+2.0*test)
imrow(result1/3.0,result2/3.0)
print amax(abs(result1-result2))

# impulse response of Gaussian filter
mask = gf(impulse,5.0)
imrow(impulse,mask)

# convolving with the impulse response
result1 = gf(image,5.0)
mask = gf(impulse,5.0)
result2 = filters.convolve(image,mask)
imrow(result1,result2)
print amax(abs(result1-result2))

# composition of linear filters
result1 = filters.convolve(filters.convolve(image,mask),mask)
result2 = filters.convolve(image,filters.convolve(mask,mask))
imrow(result1,result2)
print amax(abs(result1-result2))

# composition of Gaussian filters
result1 = gf(gf(image,5.0),5.0)
result2 = gf(image,sqrt(5**2+5**2))
imrow(result1,result2)
print amax(abs(result1-result2))

# horizontal and vertical Gaussians
vf = gf(impulse,(5.0,0.0))
hf = gf(impulse,(0.0,5.0))
imrow(impulse,vf,hf,s=4)

# convolution of horizontal and vertical Gaussians
imrow(vf,hf,filters.convolve(vf,hf),s=4)

# separable convolution
imrow(impulse,gf(impulse,(5.0,5.0)),gf(gf(impulse,(5.0,0.0)),(0.0,5.0)),s=4)

# separable convolution
imrow(gf(image,(5.0,5.0)),gf(gf(image,(5.0,0.0)),(0.0,5.0)),s=4)

# speed of separable vs full 2D convolution
import timeit
print "2D convolution",
print timeit.timeit(lambda: filters.convolve(image,mask),number=1)
print "separable convolution",
print timeit.timeit(lambda: gf(gf(image,(5.0,0.0)),(0.0,5.0)),number=100)/100.0
print "Gaussian library",
print timeit.timeit(lambda: gf(image,(5.0,5.0)),number=100)/100.0

# Gaussians
xs = arange(-5.0,5.0,0.25)
ys = exp(-xs**2/2.0)
plot(xs,ys)

mask = outer(ys,ys)
mask /= sum(mask)
imshow(mask)

# convolution with an explicitly constructed Gaussian mask
imrow(image,filters.convolve(image,mask))

# nested loop implementation of a box filter
def dumb_boxfilter(image,r):
    r /= 2
    w,h = image.shape
    out = zeros(image.shape)
    for i in range(r,w-r):
        for j in range(r,h-r):
            out[i,j] = mean(image[i-r:i+r,j-r:j+r].ravel())
    return out

# library vs nested loop for box filter
imrow(filters.uniform_filter(image,20),dumb_boxfilter(image,20))

# data parallel implementation of box filter
def myaverage(image,r):
    out = zeros(image.shape)
    for i in range(-r/2,r/2):
        for j in range(-r/2,r/2):
            out += roll(roll(image,i,axis=0),j,axis=1)
    return out*1.0/(r*r)

imrow(filters.uniform_filter(image,10),myaverage(image,10))

# data parallel implementation of convolution
def myconvolve(image,mask):
    mw,mh = mask.shape
    out = zeros(image.shape)
    for i in range(mw):
        for j in range(mh):
            out += mask[i,j]*roll(roll(image,mw/2-i,axis=0),mh/2-j,axis=1)
    return out

# library vs data parallel convolution
imrow(filters.convolve(image,mask),myconvolve(image,mask))