Graph-based image processing

We will analyze this dataset of images

In [78]:
from glob import glob
from PIL import Image
from numpy import *
import matplotlib.pyplot as pl
%pylab inline

image_files = glob('images/*.pgm')# set of all pictures names
print(image_files,'\n')

image = 255-array(Image.open(image_files[3]))
h,w = array(image).shape
print('size = ',h,w)

Populating the interactive namespace from numpy and matplotlib
['images/airplane.pgm', 'images/fgr.pgm', 'images/pentagon.pgm', 'images/clown.pgm', 'images/peppers.pgm', 'images/feathers.pgm', 'images/couple.pgm', 'images/goldhill.pgm', 'images/mandrill.pgm', 'images/boats.pgm', 'images/cameraman.pgm', 'images/lena.pgm', 'images/tibet.pgm', 'images/house.pgm'] 

size =  512 512

Build a patch embedding

In [79]:
# Define patch properties

s = int(9)
d = s**2
patch_filter = ones((s,s))
spacer = 9
min_pix = int((s-1)/2)
max_pix = int(h - (s-1)/2)

x_pos = range(min_pix,max_pix,spacer)
y_pos = range(min_pix,max_pix,spacer)
nn = len(x_pos)
N = nn**2
d = s**2

X = zeros((N,d))
imshow(image,cmap='Greys')
shape(image)
Out[79]:
(512, 512)
In [80]:
i = 10
j=100
i_ids = linspace(i-(s-1)/2,i+(s-1)/2,(s),dtype=int)
j_ids = linspace(j-(s-1)/2,j+(s-1)/2,(s),dtype=int)
patch = 255-(array(image)[i_ids,])[:,j_ids]
imshow(patch,cmap='Greys')
Out[80]:
<matplotlib.image.AxesImage at 0x145f93e10>
In [81]:
counter = 0
for i in x_pos:
    for j in y_pos:
        i_ids = linspace(i-(s-1)/2,i+(s-1)/2,(s),dtype=int)
        j_ids = linspace(j-(s-1)/2,j+(s-1)/2,(s),dtype=int)
        patch = (image[i_ids,])[:,j_ids]
        X[counter,:] = reshape(patch,(1,d))
        counter +=1
In [82]:
shape(X)
Out[82]:
(3136, 81)

Lets view a scatter plot of the patches

In [83]:
# figure(figsize=(190,190))
# for i in range(0,9,1):
#     for j in range(0,9,1):
#         # plot the i,j coordinate plane projections
#         subplot(d,d,i*d+j+1)
#         if i!=j: 
#             plot(X[:,j],X[:,i],'b.',alpha=0.5)
#         #xticks([])
#         #yticks([])

Not that helpful... Recall that this is a 25-dimensional space, and we are only viewing 2D projections onto pairs of the dimensions.

Instead of these 2D projections, lets consider optimal 2D projections using principle complnent analysis.

In [84]:
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
pca = PCA(n_components=3)
pca.fit(X)
pca_positions = pca.transform(X)

#plot3(pca_positions[:,0],pca_positions[:,1],pca_positions[:,2],'b.',alpha=0.5)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(pca_positions[:,0],pca_positions[:,1],pca_positions[:,2],zdir='z',s=10,c='r',marker='.')
Out[84]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x11943a780>

Make a histogram of entries

In [85]:
figure(figsize=(15,2))
for k in range(shape(pca_positions)[1]):
    subplot(1,shape(pca_positions)[1],1+k)
    a = hist(pca_positions[:,k],50)
    title('mode = '+str(k))

Let's visualize the 3 dominant modes

In [86]:
figure(figsize=(15,15))
mode = []
for k in range(shape(pca_positions)[1]):
    mode.append(reshape(pca_positions[:,k],(len(x_pos),len(x_pos))))
    ax = subplot(3,shape(pca_positions)[1],1+k)
    imshow(mode[k])
    title('mode = '+str(k)+', raw values')

    ax2 = subplot(3,shape(pca_positions)[1],1+k+shape(pca_positions)[1])
    imshow(mode[k]>80)
    title('mode = '+str(k)+', thresholded positive')
    
    ax3 = subplot(3,shape(pca_positions)[1],1+k+2*shape(pca_positions)[1])
    imshow(mode[k]<-80)
    title('mode = '+str(k)+', thresholded negative')

The 3 dominant modes pick up (i) dark vs light patches and (ii) diagonal edge detection

Let's view the extremal patches (i.e, most negative and most positive)

In [87]:
figure(figsize=(15,15))
for k in range(shape(pca_positions)[1]):
    id_max = argmax(pca_positions[:,k])
    ax = subplot(2,shape(pca_positions)[1],1+k)
    imshow(reshape(X[id_max,:],(s,s)))
    title('mode = '+str(k)+', most positive patch')

    id_min = argmin(pca_positions[:,k])
    ax = subplot(2,shape(pca_positions)[1],1+k + shape(pca_positions)[1])
    imshow(reshape(X[id_min,:],(s,s)))
    title('mode = '+str(k)+', most negative patch')

Denoising images with patch-graphs and nonlinear dimension reduction

In [88]:
import networkx as nx
from sklearn.metrics.pairwise import euclidean_distances
distances = euclidean_distances(X)
print(shape(distances))

(3136, 3136)
In [99]:
A = exp(-.005*distances) - eye(len(distances))

A[A<0.01] = 0
hist(A.max(axis=0))
#imshow(A)

d = A.sum(axis=0)
L = eye(len(d)) - dot(diag(d**-.5),dot(A,diag(d**-.5)))
min(d)
Out[99]:
0.32824367155781908
In [90]:
#G = nx.from_numpy_matrix(A)
In [100]:
import scipy.sparse as sparse
vals, vecs = sparse.linalg.eigs(sparse.csr_matrix(L), k=4,which='SR')
vecs = real(vecs[:,1:4])
vals = vals[1:4]
shape(vecs)
Out[100]:
(3136, 3)
In [101]:
figure(figsize=(15,2))
for k in range(shape(vecs)[1]):
    subplot(1,shape(vecs)[1],1+k)
    a = hist(vecs[:,k],50)
    title('mode = '+str(k))
In [93]:
figure(figsize=(15,15))
mode = []
for k in range(shape(vecs)[1]):
    mode.append(reshape(vecs[:,k],(len(x_pos),len(x_pos))))
    ax = subplot(3,shape(vecs)[1],1+k)
    imshow(mode[k])
    title('mode = '+str(k)+', raw values')

    ax2 = subplot(3,shape(vecs)[1],1+k+shape(vecs)[1])
    imshow(mode[k]>.01)
    title('mode = '+str(k)+', thresholded positive')
    
    ax3 = subplot(3,shape(vecs)[1],1+k+2*shape(vecs)[1])
    imshow(mode[k]<-.02)
    title('mode = '+str(k)+', thresholded negative')

Now, lets add noise to the image

In [109]:
noise_level = 200
noisy_image = image + noise_level*random.rand(shape(image)[0],shape(image)[1])
imshow(noisy_image,cmap='Greys')
counter = 0
N = nn**2
d = s**2
X_noisy=zeros((N,d))
for i in x_pos:
    for j in y_pos:
        i_ids = linspace(i-(s-1)/2,i+(s-1)/2,(s),dtype=int)
        j_ids = linspace(j-(s-1)/2,j+(s-1)/2,(s),dtype=int)
        patch = (noisy_image[i_ids,])[:,j_ids]
        X_noisy[counter,:] = reshape(patch,(1,d))
        counter +=1
In [110]:
distances = euclidean_distances(X_noisy)
A = exp(-.005*distances) #- eye(len(distances))
A[A<0.01] = 0
d = A.sum(axis=1)
L = eye(len(d)) - dot(diag(d**-.5),dot(A,diag(d**-.5)))

dd = 10
vals, vecs = sparse.linalg.eigs(sparse.csr_matrix(L), k=dd+1,which='SR')
vecs = real(vecs[:,1:dd+1])
vals = vals[1:dd+1]
In [111]:
projected_X = zeros((N,s**2))
for pix in range(shape(X_noisy)[1]):
    inner_prods = (dot(vecs.T,X_noisy[:,pix]))
    projected_X[:,pix] = dot(vecs,array(inner_prods))

shape(array(inner_prods))
#scatter(X[:,1],projected_X[:,1])
Out[111]:
(10,)
In [114]:
denoised_image = zeros((spacer*len(x_pos)+int((s-1)),spacer*len(y_pos)+int((s-1)) ))
print(shape(denoised_image))
counter = 0
for i in x_pos:
    for j in y_pos:
        i_ids = linspace(i-(s-1)/2,i+(s-1)/2,(s),dtype=int)
        j_ids = linspace(j-(s-1)/2,j+(s-1)/2,(s),dtype=int)
        #print(j_ids)
        temp = reshape(projected_X[counter,:],(s,s))
        shape(temp)
        denoised_image[np.ix_(i_ids, j_ids)] += temp
        counter +=1
        
imshow(denoised_image,cmap='Greys')        

(512, 512)
Out[114]:
<matplotlib.image.AxesImage at 0x11cd27da0>