Draw a image of famous Lenna and filtered image.

01_draw_lenna_and_filter.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
import mahotas as mh
import os.path

# get Lenna image
if os.path.isfile('./lenna.jpg'):
    im = mh.demos.load('lena')
    mh.imsave('lenna.jpg', im)

plt.style.use('fivethirtyeight')

def draw_image(data, size):
    Z = data.reshape(size,size)   # convert from vector to 28x28 matrix
    Z = np.array(Z[::-1,:])       # flip vertical
    plt.pcolor(Z)
    plt.gray()
    plt.tick_params(labelbottom="off")
    plt.tick_params(labelleft="off")

def gen_cos_filter(size, theta=0, freq=1, phase=0):
    if (size %2) != 0:
        tuning = 0.5
    else:
        tuning = 0

    _min = size - math.floor(size/2) - tuning
    _max = size + math.floor(size/2) + tuning
    
    xx = np.linspace(_min, _max, size)
    yy = np.linspace(_min, _max, size)
    X, Y = np.meshgrid(xx, yy)

    # Rotation 
    x_theta =  (X*freq+phase)*np.cos(theta) + (Y*freq+phase)*np.sin(theta)
    ret = np.cos(x_theta)
    return np.array(ret)*2

def add_padding_0(data, pad_length):
    if type(data) != np.ndarray:
        raise Exception('data must be numpy.ndarray')
    return np.lib.pad(data, ((0, pad_length), (0, pad_length)), 'constant', constant_values=0)

im = mh.imread('./lenna.jpg', as_grey=True)
im_size = im.shape[0]
im_lenna = im.astype(np.float32)
im_lenna = (im_lenna-np.min(im_lenna)) / float(np.max(im_lenna)-np.min(im_lenna))
im_lenna = add_padding_0(im_lenna, 16)

plt.figure(figsize=(7,7))

filter_size = 17
plt.figure(figsize=(3,3))

filter_type = 1

if filter_type == 0:
    #
    sin_filter = gen_cos_filter(filter_size, theta=np.pi/4, freq=.346, phase=-1.5)
else:
    #
    sin_filter = gen_cos_filter(filter_size, theta=-np.pi/4, freq=.646, phase=-1.5)
    

draw_image(sin_filter, filter_size)
sin_filter = sin_filter.reshape(filter_size, filter_size)

# forで回しすぎて怒られそうなコードだ・・・
converted = np.zeros(im_size*im_size).reshape(im_size, im_size)
for i in range(im_size):
    for j in range(im_size):
        converted[i,j] = np.sum([im_lenna[i+p, j+q] * sin_filter[p, q] \
                                 for p in range(filter_size) for q in range(filter_size)])
        
plt.figure(figsize=(7,7))
draw_image(converted.flatten(), im_size)
draw_image(im_lenna.flatten(), im_lenna.shape[0])

02_Lp_pooling.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('ggplot')
 
np.random.seed(0)
x_max = 1000
num = x_max * 100
P = np.linspace(1,x_max,num).astype(np.float128)
h_range = [1, 5, 10, 15]

Lp_list = []
for h in h_range:
    z = np.random.uniform(low=0,high=h,size=500)
    z = z.astype(np.float128)
    H = float(len(z))
    Lp_list.append([((1/H)*np.sum(z**p))**(1/p) for p in P])

plt.figure(figsize=(15,8))
for Lp in Lp_list:
    plt.plot(P, Lp)
plt.legend(h_range)
plt.ylim(0,15)
plt.xlabel("P")
plt.ylabel("Lp")
plt.title("Lp pooling plot with changing P")
plt.show()

03_anim_Lp_pooling_unif.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

plt.style.use('ggplot')
 
np.random.seed(0)
x_max = 1000
num = x_max * 100
# caluculation takes a certain time due to using float128
P = np.linspace(1,x_max,num).astype(np.float128)
h_range = xrange(1, 16)
z = np.array([np.random.uniform(low=0,high=h,size=500) for h in h_range])

Lp_list = []
for h in h_range:
    zz = z[h-1].astype(np.float128)
    H = float(len(zz))
    Lp_list.append([((1/H)*np.sum(zz**p))**(1/p) for p in P])

def animate(nframe):
    plt.clf()
    zmax = np.max(z[nframe])
    col = "bgrcmyk"*3
    plt.style.use('ggplot')
    plt.subplot(2,1,1)
    Lp = Lp_list[nframe]
    plt.ylim(0,15)
    plt.xlim(0,300)
    plt.xlabel("P")
    plt.ylabel("Lp")
    plt.title("Lp pooling plot with changing P, max={0:.3f}".format(zmax))
    plt.plot([0,500],[zmax,zmax],"k--")
    plt.plot(P[:30000], Lp[:30000], c=col[nframe])
    
    plt.subplot(2,1,2)
    plt.ylim(0,15.5)
    plt.title("Uniform random variable, min=0, max={0:.3f}, ave={1:.3f}".format(zmax,np.mean(z[nframe])))
    plt.plot([0,500],[zmax,zmax],"k--")
    plt.plot(z[nframe], c=col[nframe])
    
fig = plt.figure(figsize=(10,10))

anim = ani.FuncAnimation(fig, animate, frames=len(Lp_list), blit=True)
anim.save('anim_unif_Lp.gif', writer='imagemagick', fps=2, dpi=64)

04_anim_Lp_pooling_beta.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as st

plt.style.use('ggplot')
 
np.random.seed(0)
x_max = 1000
num = x_max * 100
# caluculation takes a certain time due to using float128
P = np.linspace(1,x_max,num).astype(np.float128)
h_range = xrange(1, 16)

# Beta Random Variable
a=1.5
b=20
z = np.array([st.beta.rvs(a, b, loc=0, scale=1, size=500)*h for h in h_range])
    
Lp_list = []
for h in h_range:
    zz = z[h-1].astype(np.float128)
    H = float(len(zz))
    Lp_list.append([((1/H)*np.sum(zz**p))**(1/p) for p in P])
    
def animate(nframe):
    sys.stdout.write("{}, ".format(nframe))
    plt.clf()
    zmax = np.max(z[nframe])
    col = "bgrcmyk"*3
    plt.style.use('ggplot')
    plt.subplot(2,1,1)
    Lp = Lp_list[nframe]
    plt.ylim(0,6)
    plt.xlim(0,500)
    plt.xlabel("P")
    plt.ylabel("Lp")
    plt.title("Lp pooling plot with changing P, max={0:.3f}".format(zmax))
    plt.plot([0,500],[zmax,zmax],"k--")
    plt.plot(P[:50000], Lp[:50000], c=col[nframe])
    
    plt.subplot(2,1,2)
    plt.ylim(0,6)
    plt.title("Beta({2},{3}) random variable, min=0, max={0:.3f}, ave={1:.3f}".format(zmax,np.mean(z[nframe]),a,b))
    plt.plot([0,500],[zmax,zmax],"k--")
    plt.plot(z[nframe], c=col[nframe])
    
fig = plt.figure(figsize=(10,10))

anim = ani.FuncAnimation(fig, animate, frames=len(Lp_list), blit=True)
anim.save('anim_beta_Lp.gif', writer='imagemagick', fps=2, dpi=64)

05_local_contrast_norm.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math, sys
import mahotas as mh
import scipy.stats as st

plt.style.use('fivethirtyeight')

def draw_image(data, size):
    Z = data.reshape(size,size)   # convert from vector to 28x28 matrix
    Z = np.array(Z[::-1,:])       # flip vertical
    plt.pcolor(Z)
    plt.gray()
    plt.tick_params(labelbottom="off")
    plt.tick_params(labelleft="off")

def norm_filter(size, sigma=1):
    r= math.floor(size/2)
    if (size %2) != 0:
        pad = 0
    else:
        pad = 1
    x = np.linspace(-r+pad,r, size)
    y = np.linspace(-r+pad,r, size)
    X, Y = np.meshgrid(x, y)
    
    Z = st.multivariate_normal.pdf(np.dstack((X,Y)), mean=[0,0], cov=[[sigma,0],[0,sigma]])
    plus = (1. - np.sum(Z)) / float(size*size)
    return Z + plus

def add_padding_center(data, filter_size):
    if type(data) != np.ndarray:
        raise Exception('data must be numpy.ndarray')  
        
    if (filter_size %2) != 0:
        # odd
        pad_size = (filter_size-1)/2
    else:
        #even
        pad_size = filter_size/2

    return np.lib.pad(data, ((pad_size,pad_size),(pad_size,pad_size)), 'edge')

# prepare filter
filter_size = 17
plt.figure(figsize=(7,7))
draw_image(padded_lenna.flatten(), im_lenna.shape[0] + filter_size-1)
norm_f = norm_filter(filter_size, sigma=15)
plt.figure(figsize=(5,5))
draw_image(norm_f.flatten(), filter_size)

# prepare Lenna image
im = mh.imread('./lenna-ring.jpg', as_grey=True)
im_size = im.shape[0]
im_lenna = im.astype(np.float32)
im_lenna = (im_lenna-np.min(im_lenna)) / float(np.max(im_lenna)-np.min(im_lenna))
padded_lenna = add_padding_center(im_lenna, filter_size)

# forで回しすぎて怒られそうなコードだ・・・
for i in range(im_size):
    sys.stdout.write("i:{}".format(i))
    for j in range(im_size):
        im_lenna[i,j] -= np.sum([padded_lenna[i+p, j+q] * norm_f[p, q] \
                                 for p in range(filter_size) for q in range(filter_size)])
        
plt.figure(figsize=(7,7))
draw_image(im_lenna.flatten(), im_size)