Perona Malik diffusion in Cython

# cython: profile=False
# file: nld_aux.pyx

from __future__ import division
import numpy as np
import scipy
cimport numpy as np
#DTYPE = np.int
DTYPE = np.float32
ctypedef np.float32_t DTYPE_t

cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def nldiffusion(np.ndarray[DTYPE_t, ndim=3] img not None,
                np.ndarray[DTYPE_t, ndim=3] gradx not None,
                np.ndarray[DTYPE_t, ndim=3] grady not None,
                np.ndarray[DTYPE_t, ndim=2] cgradx not None,
                np.ndarray[DTYPE_t, ndim=2] cgrady not None,
                double p1, double p2):

    cdef int l,k

    for l in range(0,img.shape[0]):
        for k in range(0,img.shape[1]-1):
            gradx[l,k,0] = img[l,k+1,0] - img[l,k,0]
            gradx[l,k,1] = img[l,k+1,1] - img[l,k,1]
            gradx[l,k,2] = img[l,k+1,2] - img[l,k,2]
            cgradx[l,k] = p2*(fabs(gradx[l,k,0])+fabs(gradx[l,k,1])+fabs(gradx[l,k,2]))
            cgradx[l,k] = 1. / (1.+cgradx[l,k]*cgradx[l,k])
            # cgradx[l,k] = p2

    for l in range(0,img.shape[0]-1):
        for k in range(0,img.shape[1]):
            grady[l,k,0] = img[l+1,k,0] - img[l,k,0]
            grady[l,k,1] = img[l+1,k,1] - img[l,k,1]
            grady[l,k,2] = img[l+1,k,2] - img[l,k,2]
            cgrady[l,k] = p2*(fabs(grady[l,k,0])+fabs(grady[l,k,1])+fabs(grady[l,k,2]))
            cgrady[l,k] = 1. / (1.+cgrady[l,k]*cgrady[l,k])
            # cgrady[l,k] = p2

    for l in range(1,img.shape[0]-1):
        for k in range(1,img.shape[1]-1):
            img[l  ,k  ,0] += p1*(\
                -cgradx[l  ,k-1]*gradx[l  ,k-1,0] +\
                 cgradx[l  ,k  ]*gradx[l  ,k  ,0] +\
                -cgrady[l-1,k  ]*grady[l-1,k  ,0] +\
                 cgrady[l  ,k  ]*grady[l  ,k  ,0]     )
            img[l  ,k  ,1] += p1*(\
                -cgradx[l  ,k-1]*gradx[l  ,k-1,1] +\
                 cgradx[l  ,k  ]*gradx[l  ,k  ,1] +\
                -cgrady[l-1,k  ]*grady[l-1,k  ,1] +\
                 cgrady[l  ,k  ]*grady[l  ,k  ,1]     )
            img[l  ,k  ,2] += p1*(\
                -cgradx[l  ,k-1]*gradx[l  ,k-1,2] +\
                 cgradx[l  ,k  ]*gradx[l  ,k  ,2] +\
                -cgrady[l-1,k  ]*grady[l-1,k  ,2] +\
                 cgrady[l  ,k  ]*grady[l  ,k  ,2]     )


cdef extern double exp(double)
cdef extern double fabs(double)