from scipy.ndimage import filters
gray()

blob = mean(imread("blob.png"),2)
imshow(blob)

f = filters.gaussian_gradient_magnitude(blob,2.0)
imshow(f)

fx = filters.gaussian_filter(f,2.0,(1,0))
fy = filters.gaussian_filter(f,2.0,(0,1))
imshow(fx)

u = fx.copy()
v = fy.copy()
mu = 0.1
for i in range(10000):
    m = (fx**2+fy**2)
    ut = mu*filters.laplace(u)-(u-fx)*m
    vt = mu*filters.laplace(v)-(v-fy)*m
    u += ut
    v += vt

# result of GVF gradient spreading
subplot(121); imshow(fx)
subplot(122); imshow(u)

t = linspace(0,2*pi,50)
path = array([150*cos(t)+200,200*sin(t)+300]).T
imshow(u); plot(path[:,1],path[:,0]); axis("off")

# iterative deformation of contour by GVF
t = linspace(0,2*pi,50)
path = array([150*cos(t)+200,200*sin(t)+300]).T
for t in range(10000):
    dx = array([u[int(y),int(x)] for y,x in path])
    dy = array([v[int(y),int(x)] for y,x in path])
    delta = array([dx,dy]).T
    path += 30.0*delta
imshow(u); plot(path[:,1],path[:,0]); axis("off")

# deformation with smoothness term
t = linspace(0,2*pi,50)
path = array([150*cos(t)+200,200*sin(t)+300]).T
for t in range(30000):
    x2 = filters.gaussian_filter(path,1.0,(2,0),mode='wrap')
    x4 = filters.gaussian_filter(x2,1.0,(2,0),mode='wrap')
    dx = array([u[int(y),int(x)] for y,x in path])
    dy = array([v[int(y),int(x)] for y,x in path])
    delta = array([dx,dy]).T
    path += 0.001*x2-0.001*x4+10.0*delta
imshow(u); plot(path[:,1],path[:,0]); axis("off")