from pylab import *
from scipy.ndimage import filters
def F(a): return array(a,'f')
figsize(12,6)

temp = F(mean(imread("page.png"),2))
roi = (slice(500,1000),slice(500,1000))
target = F(1-(temp[roi]>0.5))
imshow(target)

image = filters.gaussian_filter(target+0.2*randn(*target.shape),3.0)
image -= amin(image); image /= amax(image)
imshow(image)

def sigmoid(x): return 1/(1+exp(-x))

r=10
filter = F(0.01*randn(2*r+1,2*r+1))

pred = sigmoid(filters.convolve(image,filter))
subplot(121); imshow(pred)
subplot(122); imshow(pred,vmin=0,vmax=1)
amin(pred),amax(pred)

delta0 = (target-pred)
subplot(121); imshow(delta0)
delta = (target-pred)*pred*(1-pred)
subplot(122); imshow(delta)

dw = array([[sum(delta*roll(roll(image,i,0),j,1)) for j in range(-r,r+1)] for i in range(-r,r+1)])
dw /= prod(image.shape)
imshow(dw)
amin(dw),amax(dw),dw.dtype

theta = 0.0
for iter in range(1000):
    pred = sigmoid(filters.convolve(image,filter)+theta)
    err = sum((pred-target)**2)
    delta = (pred-target)*pred*(1-pred)
    delta /= prod(image.shape)
    dw = array([[sum(delta*roll(roll(image,i,0),j,1)) for j in range(-r,r+1)] for i in range(-r,r+1)])
    if iter%50==0:
        print iter,err,":",(amin(pred),amax(pred)),sum(abs(delta)),(amin(dw),amax(dw)),(amin(filter),amax(filter)),theta
    filter -= dw
    theta -= sum(delta)

figsize(20,6)
subplot(131); imshow(target)
subplot(132); imshow(image)
subplot(133); imshow(sigmoid(filters.convolve(image,filter)+theta))

imshow(filter)