TI/Programowanie dla Fizyków Medycznych/Morfologia matematyczna

import numpy as np
import pylab as py
a=np.zeros((100,100),dtype=np.bool)
a[30:50,30:50]=True
a[50:70,50:70]=True
for x in range(100):
    for y in range(100):
        if (np.random.random()<0.2): a[x,y]=False
        if (np.random.random()>0.8): a[x,y]=True
        


#brush=np.array([[0,1,1,1,0],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[0,1,1,1,0]],dtype=np.bool)
#brush=np.array([[0,0,1,0,0],[0,1,1,1,0],[1,1,1,1,1],[0,1,1,1,0],[0,0,1,0,0]],dtype=np.bool)
brush=np.ones((5,5))
print brush

def brush2list(brush):
    result=[]
    N=brush.shape[0]
    middle=N/2
    for x in range(N):
        for y in range(N):
            if brush[x,y]: result.append((x-middle,y-middle))
    return result

def medianowy(fig,brush=np.array([[0,1,0],[1,1,1],[0,1,0]],dtype=np.bool)):
    result=np.zeros(fig.shape)
    brush_list=brush2list(brush)
    for x in range(3,fig.shape[0]-3):
        for y in range (3,fig.shape[1]-3):
            result[x,y]=np.median([fig[x+x_shift,y+y_shift] for (x_shift,y_shift) in brush_list])
    return result


def dylacja(fig,brush=np.array([[0,1,0],[1,1,1],[0,1,0]],dtype=np.bool)):
    result=np.zeros(fig.shape)
    brush_list=brush2list(brush)
    for x in range(3,fig.shape[0]-3):
        for y in range (3,fig.shape[1]-3):
            result[x,y]=max([fig[x+x_shift,y+y_shift] for (x_shift,y_shift) in brush_list])
    return result

def erozja(fig,brush=np.array([[0,1,0],[1,1,1],[0,1,0]],dtype=np.bool)):
    result=np.zeros(fig.shape)
    brush_list=brush2list(brush)
    for x in range(3,fig.shape[0]-3):
        for y in range (3,fig.shape[1]-3):
            result[x,y]=min([fig[x+x_shift,y+y_shift] for (x_shift,y_shift) in brush_list])
    return result

def otwarcie(fig,brush=np.array([[0,1,0],[1,1,1],[0,1,0]],dtype=np.bool)):
    return dylacja(erozja(fig,brush),brush)

def zamkniecie(fig,brush=np.array([[0,1,0],[1,1,1],[0,1,0]],dtype=np.bool)):
    return erozja(dylacja(fig,brush),brush)


py.imshow(a, cmap='Greys',  interpolation='nearest')
py.show()
py.imshow(medianowy(a,brush), cmap='Greys',  interpolation='nearest')
py.show()
py.imshow(medianowy(medianowy(a,brush),brush), cmap='Greys',  interpolation='nearest')
py.show()
py.imshow(medianowy(medianowy(medianowy(a,brush),brush),brush), cmap='Greys',  interpolation='nearest')
py.show()