eigenface

eigenface.py

import argparse
import numpy as np
from scipy import linalg

def ReadTillWhite(fd):
    data = b''
    while True:
        ch = fd.read(1)
        if ch in b"\r\t\n ":
            break
        data += ch
    return data

def LoadPGM(pgm):
    f = open(pgm, "rb")
    magic = ReadTillWhite(f)
    w = int(ReadTillWhite(f))
    h = int(ReadTillWhite(f))
    g = int(ReadTillWhite(f))
    pixels = []
    while True:
        ch = f.read(1)
        if ch == b"": break
        pixels.append(ord(ch))
    f.close()
    return w, h, g, pixels

def main():
    parser = argparse.ArgumentParser()

    parser.add_argument(dest="imgs",
            nargs="+",
            metavar="1.pgm 2.pgm ...", help="ORL Image Data")
    args = parser.parse_args()

    ntrain = 5
    ntest = 10 - ntrain
    trainset = []
    testset = []
    for img in args.imgs:
        num = img.split("/")[-1].split(".")[0]
        if int(num) <= ntrain:
            trainset.append(img)
        else:
            testset.append(img)
    print("total images: %d" % len(args.imgs))
    print("train images: %d" % len(trainset))
    print("test images:  %d" % len(testset))

    w, h = 0, 0
    vecs = []
    mean = np.array([])
    for img in trainset:
        w, h, g, pixels = LoadPGM(img)
        pixels = np.array(pixels)
        vecs.append(pixels)
        if len(mean) != 0:
            mean = mean + pixels
        else:
            mean = pixels
    mean = np.array([float(pixel/len(trainset)) for pixel in mean])
    vecs = [(vec - mean) for vec in vecs]

    # play a little mathematical trick
    A_T = np.asmatrix(vecs)
    A = A_T.T
    C_prime = A_T * A
    la, V_prime = linalg.eig(C_prime)
    la = np.real(la)
    idxlst = np.argsort(-1 * la)
    V_prime = V_prime[idxlst]
    la = la[idxlst]
    # real eigenvectors and eigenmatrix
    V = [np.asarray(np.dot(A, v))[0] for v in V_prime]
    U_T = np.asmatrix(V)[10:190,:]

    # project original training vectors into a new space
    x_stars = []
    for x in vecs:
        x_star = np.asarray(np.dot(U_T, x)[0])
        x_stars.append(x_star)

    # test rest images
    ncorrect = 0
    for im in range(len(testset)):
        w, h, g, pixels = LoadPGM(testset[im])
        pixels = np.array(pixels)
        t_star = np.asarray(np.dot(U_T, (pixels-mean))[0])
        
        midx = -1
        mdst = float("inf")
        for ist in range(len(x_stars)):
            dst = linalg.norm(t_star - x_stars[ist])
            if mdst > dst:
                mdst = dst
                midx = ist
        if int(midx/ntrain) == int(im/ntest):
            ncorrect += 1

    print("correct:  %d" % ncorrect)
    print("accuracy: %.2f%%" % (100.0*ncorrect/len(testset)))

if __name__ == "__main__":
    main()