Week 10: Unsupervised Learning (Part II)

CSCI-UA 9473 - Introduction to Machine Learning

Partial Solutions

Exercise 1. White matter, Gray matter and cerebrospinal fluid¶

For this exercise, you will need to install nibabel. You can do this by typing the line 'pip install nibabel' in your terminal.

We consider the MR image shown below. We want to segment this image into the white matter, gray matter and cerebrospinal fluid.

Start by plotting the historgram of all pixel values (you can do this with 'numpy.histogram'). How does the distribution look like ?
Use your histogram to initialize the E-M algorithm (means and variance). The plot the resulting distributions on top ofthe histogram. Does your algorithm correctly capture the three clusters?

In [214]:

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt

img = nib.load('T1.nii')
plt.imshow(img.dataobj[:,80,:], cmap = 'gray')

np_array = np.asarray(img.dataobj[:,80,:])

plt.show()

In [40]:

#### Exercise 2. Spectral clustering

pixel_list = np.reshape(np_array,np.shape(np_array)[0]*np.shape(np_array)[1])

pixel_list2 = np.delete(pixel_list, np.argwhere(pixel_list==0))
plt.hist(pixel_list2, bins=np.linspace(np.min(pixel_list2), np.max(pixel_list2),100), density=True)
plt.title('Original pixels histogram')
plt.show()

In [47]:

# finding the mask

mask = 0*np.ones(np.shape(np_array))

for i in np.arange(np.shape(np_array)[0]):
    for j in np.arange(np.shape(np_array)[1]):
    
        if np_array[i,j] >0:
            
            mask[i,j]=1
            
plt.imshow(mask,cmap = 'gray')
plt.show()

In [98]:

import cv2

kernel = np.ones((5,5), np.uint8)

img_erosion = cv2.erode(mask, kernel, iterations=2)
#img_dilation = cv2.dilate(mask, kernel, iterations=1)

# for display
kernel2 = np.ones((4,4), np.uint8)
img_erosion2 = cv2.erode(img_erosion, kernel2, iterations=1)
border = img_erosion - img_erosion2 

#image_contour = np.zeros(np.shape)
color_tmp = np.zeros((np.shape(np_array)[0], \
                      np.shape(np_array)[1], \
                     3))

for i in np.arange(3):
    
    color_tmp[:,:,i] = np_array[:,:]/np.max(np_array)

border_c = np.squeeze(np.where(border==1))

for k in np.arange(len(border_c[0,:])):
    
    color_tmp[border_c[0,k],border_c[1,k],0:1] = 0 
    color_tmp[border_c[0,k],border_c[1,k],2] = 1 
    
plt.imshow(color_tmp)
plt.show()

In [99]:

ind1 = np.squeeze(np.where(img_erosion==1))
np.shape(ind1)

Out[99]:

(2, 16900)

In [100]:

cropped_image = np.zeros(np.shape(img_erosion))

for k in np.arange(len(ind1[0,:])):

    if img_erosion[ind1[0,k], ind1[1,k]]==1:
        cropped_image[ind1[0,k], ind1[1,k]] = np_array[ind1[0,k], ind1[1,k]]

plt.imshow(cropped_image, cmap = 'gray')
plt.show()

In [223]:

pixel_list = np.reshape(cropped_image,np.shape(cropped_image)[0]*np.shape(cropped_image)[1])

pixel_list2 = np.delete(pixel_list, np.argwhere(pixel_list==0))
plt.hist(pixel_list2, bins=np.linspace(np.min(pixel_list2), np.max(pixel_list2),100), density=True)
plt.title('Original pixels histogram')
plt.show()

In [195]:

# optional median filter 
from scipy import ndimage, misc 
filtered_image = ndimage.median_filter(cropped_image,[2,2])

plt.figure()
plt.imshow(filtered_image, plt.cm.gray)
plt.title('Reconstructed Image')

Out[195]:

Text(0.5, 1.0, 'Reconstructed Image')

In [200]:

pixel_list = np.reshape(filtered_image,np.shape(cropped_image)[0]*np.shape(cropped_image)[1])

pixel_list2 = np.delete(pixel_list, np.argwhere(pixel_list<100))
plt.hist(pixel_list2, bins=np.linspace(np.min(pixel_list2), np.max(pixel_list2),50), density=True)
plt.title('Original pixels histogram')
plt.show()

In [254]:

import copy 

# Initialization

mu1 = 100
mu2 = 200
mu3 = 300

mu = np.asarray([mu1, mu2, mu3])

sigma1 = 20
sigma2 = 20
sigma3 = 10

sigma = np.asarray([sigma1, sigma2, sigma3])

phi1 = .8
phi2 = .1
phi3 = .1

phi = np.asarray([phi1, phi2, phi3])

num_clusters = 3

criterion = 10

maxIter = 500

criterion_mu = np.zeros((maxIter,))
criterion_sigma = np.zeros((maxIter,))
criterion_phi = np.zeros((maxIter,))
    
for kk in np.arange(maxIter):

    # E-step 
    
    gamma = np.zeros((len(pixel_list2), num_clusters))
    numerator = np.zeros((len(pixel_list2), num_clusters))
    denominator = np.zeros((len(pixel_list2),))
    
    for j in np.arange(num_clusters):

        numerator[:,j] = np.squeeze(np.exp(-np.power(pixel_list2- mu[j], 2.) / (2 * np.power(sigma[j], 2.)))*phi[j])
        denominator+= numerator[:,j]

    for j in np.arange(num_clusters):
            
        gamma[:,j] = np.divide(numerator[:,j],denominator)

    # M-step 
    
    mu_old = copy.deepcopy(mu)
    sigma_old = copy.deepcopy(sigma)
    phi_old = copy.deepcopy(phi)
    
    for j in np.arange(num_clusters):

        tmp = np.dot(gamma[:,j],pixel_list2)
        mu[j] = tmp/np.sum(gamma[:,j])

        tmp = np.dot(gamma[:,j],(pixel_list2 - mu_old[j])**2)/np.sum(gamma[:,j])
        sigma[j] = np.sqrt(tmp)

        phi[j] = np.sum(gamma[:,j])/len(pixel_list2)
        
    criterion_mu[kk] = np.sqrt(np.sum((mu.reshape(-1,1)-mu_old.reshape(-1,1))**2))
    criterion_sigma[kk] = np.sqrt(np.sum((sigma.reshape(-1,1)-sigma_old.reshape(-1,1))**2))
    criterion_phi[kk] = np.sqrt(np.sum((phi.reshape(-1,1)-phi_old.reshape(-1,1))**2))

In [255]:

plt.semilogy(criterion_mu+criterion_sigma+criterion_phi)
plt.title('Convergence of the EM Algorithm')
plt.show()

In [256]:

import scipy.stats as stats
import math

numClasses = 3
plt.hist(pixel_list2, bins=np.linspace(np.min(pixel_list2), np.max(pixel_list2),100), density=True, alpha=0.5)

for i in np.arange(numClasses):
    
    
    x = np.linspace(mu[i] - 3*sigma[i], mu[i] + 3*sigma[i], 100)
    plt.plot(x, stats.norm.pdf(x, mu[i], sigma[i]))
    
plt.title('Recovered Gaussian Mixture')
plt.show()

To recover the independent regions, we take for each pixel the Gaussian that returns the highest probability

In [259]:

import numpy.matlib

max_row = np.arange(np.shape(np_array)[0]).reshape(-1,1)
max_col = np.arange(np.shape(np_array)[1]).reshape(-1,1).T

row_ind = np.tile(max_row, [1,np.shape(np_array)[1]])
col_ind = np.tile(max_col, [np.shape(np_array)[0],1])

list_rind = np.reshape(row_ind,np.shape(np_array)[0]*np.shape(np_array)[1])
list_cind = np.reshape(col_ind,np.shape(np_array)[0]*np.shape(np_array)[1])


list_rind = np.delete(list_rind,np.argwhere(pixel_list<5))
list_cind = np.delete(list_cind,np.argwhere(pixel_list<5))

segmented_image = 3*np.ones(np.shape(np_array))

for i in np.arange(len(list_rind)):
    
    segmented_image[list_rind[i], list_cind[i]] = np.argmax(gamma[i,:])
    
plt.imshow(segmented_image, cmap='gray')
plt.show()

In [130]:

print(np.max(gamma[:,2]))

4.565117792090513e-93

Exercise 2. Image compression¶

2.1. plotting the eigenvalues¶

Consider the simple image given below. We would like to compress this image by means of principal component factorization. For this purpose, we will see the columns as our feature vectors. Build the covariance matrix associated to the columns of the image (don't forget to center them first) and plot the eigenvalues associated to this covariance matrix. What do you notice?

In [2]:

import numpy as np
from PIL import Image
from matplotlib import pyplot as plt

img = Image.open('shapes.png')
data1 = np.array(img)

plt.imshow(data1)
plt.show()

In [3]:

mean_col = np.mean(img, axis=1)

centered_image = img - np.matmul(mean_col.reshape(-1,1), np.ones((1, np.shape(img)[1])))

covariance = (1/np.shape(img)[1]) * np.matmul(centered_image, centered_image.T)

w, v = np.linalg.eig(covariance)

# sort the eigenvalues in decreasing order

wsorted_indices = np.argsort(w)[::-1]

w_sorted = w[wsorted_indices]
v_sorted = v[:,wsorted_indices]

plt.semilogy(np.real(w_sorted))
plt.show()

2.2. projection onto the principal vectors¶

Given your result from the previous question, compress the image by computing its principal component factorization for increasing numbers of components. Compare the reconstructions.

In [4]:

subspace = np.real(v_sorted[:,:8])

reconstructed_image = np.matmul(subspace, np.matmul(subspace.T, img))

plt.imshow(reconstructed_image)
plt.show()

Exercise 3: Face recognition¶

In this exercise, we consider the 'Labeled Faces in the Wild' dataset from UM Amherst. We would like to train a Maximum margin classifier on this dataset (scikit learn 'svm') with a Gaussian kernel. In its original form, the dimensionality of the data is too high (if we used the raw images as our feature vectors, each vector would encode 75$\times$56 features). Although it can sometimes be a blessing, high dimensionality can also be a curse. Working in a high dimensional space means that the feature vectors will be very distant from each other which could result in a meaningless classifier. To avoid this, we will start by computing a principal component decomposition of our set of images. The SVM classifier will then be applied on the PCA representation of the training images.

In [5]:

import numpy as np

from sklearn.datasets import fetch_lfw_people
from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt

lfw_people = fetch_lfw_people(min_faces_per_person=70, resize=.6)

n_samples, h, w = lfw_people.images.shape

images = lfw_people.data
targets = lfw_people.target

training_images, test_images, train_targets, test_targets = \
train_test_split(images,targets, test_size = .10, random_state=10)

plt.imshow(training_images[0,:].reshape(h,w), cmap='gray')
plt.axis('off')
plt.show()

Question 3.1. PCA decomposition¶

Extract the first 50, 100, 200 and 300 principal components from the image dataset and display the corresponding reconstruction for one of the faces.
Display the plot of the reconstruction error (deviation between the true image and its PCA reconstruction for K=1 to h*w)
Finally, plot each of the principal images for K=1,...10

In [6]:

import numpy as np
from sklearn.decomposition import PCA

# put your code here 

import numpy as np
from sklearn.decomposition import PCA


pca = PCA(n_components=300, whiten=True)
latent_rep = pca.fit_transform(training_images)

reconstructed_image0 = pca.inverse_transform(latent_rep[2,:])

plt.subplot(1, 2, 1)
plt.imshow(reconstructed_image0.reshape(h, w))

plt.subplot(1, 2, 2)
plt.imshow(training_images[2,:].reshape(h, w))

plt.show()