from pylab import *
from scipy.ndimage import filters,measurements
from scipy import stats

def trivedges(m,threshold=1e-3):
    return abs(m-roll(m,1,0))+abs(m-roll(m,1,1))>threshold

def mondrian(dims=(200,200),sr=(20,50),n=100,noise=0.5):
    h,w = dims
    image = zeros(dims)
    for i in range(n):
        y = randint(h)
        x = randint(w)
        ry = randint(*sr)
        rx = randint(*sr)
        y0 = clip(y-ry,0,h)
        y1 = clip(y+ry,0,h)
        x0 = clip(x-rx,0,w)
        x1 = clip(x+rx,0,w)
        image[y0:y1,x0:x1] = rand()
    edges = trivedges(image)
    return image+noise*randn(*image.shape),edges

image,edges = mondrian(noise=0.1)
subplot(121); imshow(image,cmap=cm.gray)
subplot(122); imshow(edges,cmap=cm.gray)

imshow(trivedges(m,threshold=0.5),cmap=cm.gray)

def evaluate_edges(truth,pred,b=5):
    cases = array(2*(truth!=0)+(pred!=0),'i')[b:-b,b:-b]
    c00,c01,c10,c11 = measurements.sum(ones(cases.shape),cases,arange(4))
    precision = c11 * 1.0 / (c01 + c11)
    recall = c11 * 1.0 / (c10 + c11)
    fmeasure = precision * recall / (precision + recall)
    return fmeasure,precision,recall



evaluate_edges(edges,trivedges(m,threshold=0.2))

es = []
ts = linspace(0.0,1.5,100)
for threshold in ts:
    e = evaluate_edges(edges,trivedges(image,threshold=threshold))
    es.append(e)
es = array(es)
plot(ts,es[:,0])
plot(ts,es[:,1])
plot(ts,es[:,2])

def blobworld(dims=(200,200),sigma=10.0,frac=0.5,noise=0.5):
    h,w = dims
    image = randn(*dims)
    image = filters.gaussian_filter(image,sigma)
    threshold = stats.scoreatpercentile(image,frac*100)
    image = (image>threshold)
    image,n = measurements.label(image)
    image = rand(n+1)[image]
    edges = trivedges(image)
    return image+noise*randn(*image.shape),edges

bw,be = blobworld()
subplot(121); imshow(bw,cmap=cm.gray)
subplot(122); imshow(be,cmap=cm.gray)