An Implementation of the robust subspace projection technique described in: A. Leonardis and H. Bischof. Robust recognition using eigenimages. Computer Vision and Image Understanding, 78:99--118, 2000. based off of source provided by the original authors

robust_pca_project.py

#!/usr/bin/env python

# robust_pca_project.py
#
# Author: Kevin Hughes
# Date: May 2013 

"""
An Implementation of the robust subspace projection technique described in:

A. Leonardis and H. Bischof. Robust recognition using eigenimages.
Computer Vision and Image Understanding, 78:99--118, 2000.

Implementation based off of source provided by the original authors
"""

import numpy as np

def robust_project(X, Xm, U, k=300, r=0.85, s=100, n=3):
	""" robust_project """
	# X	... image vector
	# Xm ... mean vector
	# U	... eigenvectors
	# k	... k*K is the initial number of selected pixels
	# r	... reduction factor in alpha trimming
	# s	... s*K is stopping criterion
	# n	... number of iterations 
	#return ... resulting coefficient vector
	
	""" authors params: k=48, r=0.85, s=24, n=3 """
	
	if X.shape != Xm.shape:
		raise ValueError("robust_project X and Xm must be the same shape")
	
	if X.size != U.shape[1]:
		raise ValueError("robust_project X must have the same number of cols as U")		
	
	K,M = U.shape
	
	if k*K > X.size:
		raise ValueError("robust_project k is too large, k*K is > than len(X)")
	
	if s > k:
		raise ValueError("robust_project s must be less than k")
	
	if r <= 0 or r >= 1:
		raise ValueError("robust_project r must be between 0 and 1")
	
	# mean center X
	Xd = X - Xm
	
	# run n iterations to get n candidate projection coeffs
	a_cands = []
	errors = []
	for i in range(n):
	
		# get k*K random points
		idx = np.arange(M)
		np.random.shuffle(idx)
		idx = idx[:(k*K)]
		nbp = len(idx)
	
		B = Xd[:,idx]
		A = U[:,idx]
	
		# solve lstq squares for the random subset
		x, residuals, rank, s_vals = np.linalg.lstsq(A.transpose(),B)	
		err = np.absolute(np.dot(A.transpose(),x) - B)

		#alpha trimming
		while(nbp > s*K): 
			errh = np.sort(err)
			err_marg = errh[np.ceil(r*nbp)]			
		
			# trim
			idx = np.where(err <= err_marg)[0]
			nbp = len(idx)
		
			B = B[:,idx]
			A = A[:,idx]
		
			# solve lstq squares for the trimmed random subset
			x, residuals, rank, s_vals = np.linalg.lstsq(A.transpose(),B)	
			err = np.absolute(np.dot(A.transpose(),x) - B)
	
		a_cands.append(x.copy())
		errors.append(np.sum(err) / nbp)
	
	# Select best the projection coefficients
	idx = np.argsort(errors)
	return a_cands[idx[0]]
	
def face_example():
	"""face_example
	demonstrates the drawbacks of naive lst sqrs projection
	when corrupted data is present. Shows how the robust projection
	function handles this problem.
        This example will download and extract the orl faces dataset. Then it will
        perform PCA on a set of faces. The mean is reshaped, normalized and then 
        displayed. The same is done for the first 6 eigenvectors. The first face is
        projected and reconstructed, the result is then visualized.  

        Next the data is corrupted to demonstate the robust projection method. The
        corrupted face is projected and reconstructed using both the standard projection
        and the robust projection. The results are displayed and also saved to file along
        with the corrupted_face.
	"""	
	from os.path import join, exists
	import urllib
	import zipfile

	import cv2
	from sklearn.decomposition import PCA

	ORL_FACES_URL = "http://www.cl.cam.ac.uk/Research/DTG/attarchive/pub/data/att_faces.zip"
	
	# get the data
	if not exists('orl_faces'):
		print 'downloading orl faces from %s' % ORL_FACES_URL
		urllib.urlretrieve(ORL_FACES_URL, 'orl_faces.zip')
		print 'extracing orl_faces.zip'
		zfile = zipfile.ZipFile("orl_faces.zip",)
		zfile.extractall('orl_faces')

	pathToORL = 'orl_faces/'

	faces = []
	for i in range(1,20+1):
		for j in range(1,1+1):
			path =	pathToORL + 's' + str(i) + '/' + str(j) + '.pgm'
			face = cv2.imread(path,0)
			faces.append(face)
			
			#cv2.imshow("face",face)
			#cv2.waitKey()
		
	imgW = faces[0].shape[1]
	imgH = faces[0].shape[0]

	faces_rowvecs = []
	for (i,face) in enumerate (faces):
		faces_rowvecs.append(face.flatten("C").copy())
						
	faceData = np.vstack(faces_rowvecs)

	k = 50
	pca = PCA(k)
	pca.fit(faceData)
	
	# mean
	mean = pca.mean_
	mean = mean.reshape(imgH,imgW)
	mean = cv2.normalize(mean,norm_type=cv2.NORM_MINMAX)

	cv2.imshow("pca mean",mean)
	cv2.moveWindow("pca mean",(1)*(imgW+10), 100)

	# eigenvectors
	#for i in range(k):
	for i in range(6):
	
		eigenvector = pca.components_[i,:]
		eigenvector = eigenvector.reshape(imgH,imgW)
		eigenvector = cv2.normalize(eigenvector,alpha=0,beta=255,norm_type=cv2.NORM_MINMAX)
		eigenvector = eigenvector.astype(np.uint8)		

		cv2.imshow("pca"+str(i),eigenvector)
		cv2.moveWindow("pca"+str(i),(i+3)*(imgW+10), 100)
		#cv2.imwrite("eigenface_"+str(i)+".png", eigenvector)
	
	# reconstructions
	face = faces[0]
	
	cv2.imshow("face", face)
	cv2.moveWindow("face",(0+3)*(imgW+10), 300)
	
	face = face.flatten("C")
	proj = pca.transform(face)
	recon = pca.inverse_transform(proj)
	recon = recon.reshape(imgH,imgW)
	recon = cv2.normalize(recon,norm_type=cv2.NORM_MINMAX)
	
	cv2.imshow("recon", recon)
	cv2.moveWindow("recon",(1+3)*(imgW+10), 300)
	
	# corrupted data
	face = faces[0]
	#cv2.imwrite("face.png", face) 
	
	# top quarter white
	#face[0:imgH/2, 0:imgW/2] = 255
	
	# salt and pepper slash just salt
	idx = np.arange(face.size)
	np.random.shuffle(idx)
	idx = idx[:(int(face.size * 0.25))]
	face = face.flatten("C")
	face[idx] = 255
	face = face.reshape(imgH,imgW)
	
	cv2.imshow("face D", face)
	cv2.moveWindow("face D",(2+3)*(imgW+10), 300)
	cv2.imwrite("corrupted_face.png", face) 	

	# non robust projection
	face = face.flatten("C")
	proj = pca.transform(face)
	recon = pca.inverse_transform(proj)
	recon = recon.reshape(imgH,imgW)
	recon = cv2.normalize(recon,alpha=0,beta=255,norm_type=cv2.NORM_MINMAX)
	recon = recon.astype(np.uint8)	

	cv2.imshow("recon D", recon)
	cv2.moveWindow("recon D",(3+3)*(imgW+10), 300)
	cv2.imwrite("pca_reconstruction.png", recon)	

	# robust projection
	proj = robust_project(face, pca.mean_, pca.components_, 48, 0.85, 24, 100)
	recon = pca.inverse_transform(proj)
	recon = recon.reshape(imgH,imgW)
	recon = cv2.normalize(recon,alpha=0,beta=255,norm_type=cv2.NORM_MINMAX)
	recon = recon.astype(np.uint8)	
	
	cv2.imshow("recon R", recon)
	cv2.moveWindow("recon R",(4+3)*(imgW+10), 300)
	cv2.imwrite("robust_pca_reconstruction.png", recon)
	
	cv2.waitKey()		
	
if __name__ == "__main__":
	face_example()