Partial Solutions
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?
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()
#### 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()
# 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()
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()
ind1 = np.squeeze(np.where(img_erosion==1))
np.shape(ind1)
(2, 16900)
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()
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()
# 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')
Text(0.5, 1.0, 'Reconstructed Image')
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()
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))
plt.semilogy(criterion_mu+criterion_sigma+criterion_phi)
plt.title('Convergence of the EM Algorithm')
plt.show()
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
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()
print(np.max(gamma[:,2]))
4.565117792090513e-93
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?
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()
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()
Given your result from the previous question, compress the image by computing its principal component factorization for increasing numbers of components. Compare the reconstructions.
subspace = np.real(v_sorted[:,:8])
reconstructed_image = np.matmul(subspace, np.matmul(subspace.T, img))
plt.imshow(reconstructed_image)
plt.show()
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.
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()
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
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()