% CIS 601 Computer Graphics and Image Processing
% Homework 4                    Spring 2007
% Riu Baring
% ==============================================
% Main function
% ==============================================

% 
clear all

% Load a Matlab data file
load SurvTestFrames2.mat;

BLOCK3D = 3 * 8 * 8;
NUMPCA = 3;

% Array pre-allocation
blockmean = zeros(1, 1200);
blockstd = zeros(1, 1200);
data = zeros(240, 320);
matrix3d = zeros(3, 192);
projected = zeros(NUMPCA, NUMPCA);
projPCA = cell(1200, 240, 320);

% Get Projection Matrix
fprintf('PCA Initialization - Start\n');
frameWindow = ag(:,:,1:5);
for zy = 1:30
    for zx = 1:40
        zz = ((zy - 1) * 40) + zx;
        yoffset = (zy - 1) * 8;
        xoffset = (zx - 1) * 8;
        
        i3Drow = 1;
        i3Dpix = 1;
        
        for i = 1:3
            i3Dpix = 1;
            for j = 1:3
                data = frameWindow(:,:,i + (j - 1));
                for py = 1:8
                    for px = 1:8
                        matrix3d(i3Drow, i3Dpix) = data(yoffset + py, xoffset + px);
                        i3Dpix = i3Dpix + 1;
                    end
                end
            end
            i3Drow = i3Drow + 1;
        end
        [projPCA{zz}, D] = eig(cov(matrix3d));
    end
end
fprintf('PCA Initialization - Done\n');

fprintf('PCA Projection Matrix Construction - Start\n');
frameWindow = ag(:,:,1:99);
for frameCnt = 1:94
    frameWindow = ag(:,:,frameCnt:frameCnt+4);
    for zy = 1:30
        for zx = 1:40
            zz = ((zy -1) * 40) + zx;
            yoffset = (zy - 1) * 8;
            xoffset = (zx - 1) * 8;

            i3Drow = 1;
            i3Dpix = 1;
            for i = 1:3
                for j = 1:3
                    data = frameWindow(:,:,i + (j - 1));
                    for py = 1:8
                        for px = 1:8
                            matrix3d(i3Drow, i3Dpix) = data(yoffset + py, xoffset + px);
                            i3Dpix = i3Dpix + 1;
                        end
                    end
                end
                i3Drow = i3Drow + 1;
            end

            for k = 1:3
                for l = 1:3
                    projected(k, l) = zeros(1, 1);
                    coloffset = BLOCK3D - l + 1;
                    for m = 1:BLOCK3D
                        projected(k, l) = projected(k, l) + (matrix3d(k, m) * projPCA{zz}(m, coloffset));
                    end
                end
            end

            [V D] = eig(cov(projected));
            blockmean(zz) = blockmean(zz) + max(D(:));
        end
    end
end
blockmean = blockmean / 99;
stdblock = 0.3 * blockmean;
fprintf('PCA Projection Matrix Construction - Done\n');