In [71]:

    

%pylab inline
from skimage import io, exposure

plt.rcParams['image.cmap'] = 'gray_r'
plt.rcParams['image.interpolation'] = 'none'


    

    




Populating the interactive namespace from numpy and matplotlib

In [72]:

    

%qtconsole


    

In [ ]:

    

img = io.imread('unstained/view_u13v0ch0z0.png')

img = img[750:760, 950:960].astype(int)

def imshow_num(img, num, **kwargs):
    ax = subplot(num)
    im = imshow(img, **kwargs)
    colorbar(im)

    y,x = meshgrid(range(img.shape[0]),range(img.shape[1]))

    for i,j in zip(y.flatten(), x.flatten()):
        v = img[i,j]
        c = 'black'
        if v > 0.5*img.max()+0.5*img.min():
            c = 'white'
        ax.text(j,i,v,va='center',ha='center',color=c)


def imshow_arrow(angle, mag, num, **kwargs):
    ax = subplot(num)
    xlimit = (-0.5, angle.shape[1]-0.5)
    ylimit = (angle.shape[0]-0.5, -0.5)
    xlim(*xlimit)
    ylim(*ylimit)
    #imshow(img, **kwargs)
    #colorbar()

    y,x = meshgrid(range(angle.shape[0]),range(angle.shape[1]))

    for i,j in zip(y.flatten(), x.flatten()):
        v = angle[i,j]
        m = mag[i,j]
        dx = m*sin(v/180*pi)
        dy = m*cos(v/180*pi)
        
        ax.text(j,i,v,va='center',ha='center')
        ax.arrow(j,i,dx,dy,head_width=0.1,head_length=0.2)


def hog(angle, mag):
    h = np.zeros(180, dtype=float)
    y,x = meshgrid(range(angle.shape[0]),range(angle.shape[1]))
    for i,j in zip(y.flatten(), x.flatten()):
        a = angle[i,j]
        m = mag[i,j]
        h[a] += m
    return h


dx = np.zeros_like(img)
dy = np.zeros_like(img)
angles = np.zeros_like(img)

dy[1:-1,:] = img[2:,:] - img[0:-2,:]
dx[:,1:-1] = img[:,2:] - img[:,0:-2]
angles = floor(arctan2(dy,dx)/pi*180 + 90) % 180
angles = angles.astype(int)

dy = dy[1:-1, 1:-1]
dx = dx[1:-1, 1:-1]
mag = sqrt(dx**2 + dy**2)
mag = mag/mag.max()/2
angles = angles[1:-1, 1:-1]


    

In [146]:

    

fig = figure(figsize=(11,3))

imshow_num(img, 131)
title('(a) Fiber from PMT image')

imshow_num(dx, 132)
title('(b) dx')

imshow_num(dy, 133)
title('(c) dy')

tight_layout()


    

In [148]:

    

fig = figure(figsize=(12,6))
imshow_arrow(angles, mag, 121)
title('(a) Gradient illustrated as arrows')

#h = exposure.histogram(angles, nbins=angles.max())
#bar(h[1], h[0])
h = hog(angles, mag)
ax = subplot(122)
ax.plot(h)
title('(b) Histogram of weighted angles')

pass;


    

In [ ]: