PCA (eig and svd approach)

pca.py

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


def eig_decomposition(img, p=0.9):
    A = np.array(img, dtype=float)
    mean_A = np.mean(A)
    A -= mean_A
    # calculate covarience matrix
    cov = np.cov(A.T)
    # calculate eig values and eig vectors
    eig_val, eig_vec = np.linalg.eig(cov)
    # get indexs of top-k eig values
    topk = np.where(np.cumsum(sorted(eig_val, reverse=True) / np.sum(eig_val)) >= p)[0][0]
    topk = 100
    topk_idx = np.argsort(eig_val)[::-1][:topk]
    # construct projection matrix
    projection_matrix = eig_vec[:, topk_idx]
    low_dim_A = np.dot(A, projection_matrix)
    reconstructed_A = np.dot(low_dim_A, projection_matrix.T) + mean_A

    new_img = Image.fromarray((reconstructed_A).astype('uint8'))
    return new_img


def svd(img, p=0.8):
    A = np.array(img, dtype=float)
    U, sigma, VT = np.linalg.svd(A)
    topk = np.where(np.cumsum(sigma / np.sum(sigma)) >= p)[0][0]
    low_dim_A = U[:, :topk].dot(np.diag(sigma[:topk]))
    # V is actually eig_vec. Thus VT[:topk, :] is the projection_matrix
    reconstruted_A = np.dot(low_dim_A, VT[:topk, :])
    new_img = Image.fromarray(reconstruted_A.astype('uint8'))
    return new_img


if __name__ == '__main__':
    img = Image.open('./666.png')
    print 'Image loaded. Size: ', img.size
    # img.show()
    # svd(img).show()
    eig_decomposition(img).show()