Perona Malik diffusion wrapper

#!/usr/bin/python
####################################################################
## This program executes a Perona Malik anisotropic diffusion to
## smooth images preserving edges. (I've heard this version is not
## actually anisotropic, but it's definitely non-linear!)
##
## The actual diffusion iteration algorithm is at
##   http://pastebin.com/Mi21mud3
##
## Coded by: nwerneck@gmail.com
####################################################################



import scipy
import scipy.ndimage

from scipy.linalg import eig
from scipy.misc import imsave as myimsave

import numpy as np
from numpy  import *
import Image
import sys



from nld_aux import *

if __name__=='__main__':
    from pylab import *

    ion()

    rc('image', interpolation='nearest')

    ion()

    im = Image.open (sys.argv[1])
    frame = np.array(im, dtype=np.float32)
    img = np.array(im, dtype=np.float32)
    img = img[:,:,:3] #remove alpha channel

    ## Parameters for the diffusion
    alpha=0.05
    glim=1./5


    ## Initialize the arrays used in the algorithm.
    gradx   = np.zeros((img.shape[0]  ,img.shape[1]-1,3), dtype=np.float32)
    grady   = np.zeros((img.shape[0]-1,img.shape[1]  ,3), dtype=np.float32)
    cgradx  = np.zeros((img.shape[0]  ,img.shape[1]-1), dtype=np.float32)
    cgrady  = np.zeros((img.shape[0]-1,img.shape[1]  ), dtype=np.float32)

    # PErform 2400 update iterations
    for k in range(2400):
        print k
        nldiffusion(img, gradx, grady, cgradx, cgrady, alpha, glim )

        # Save one image every 10 iterations to generate 10 seconds video.
        if not k%10:
            myimsave('aa%04d.png'%(k/10), img)

    figure(1)
    imshow(frame/256.)
    figure(2)
    imshow(img/256.)