clear
% PCA for classification with cross validation.

% Constant
MnumComponents = 20;
sx = 85; sy = 127; % original image size (this would need much memory !) 
% sx = 20; sy = 30; % resized image size

addpath(['../']);

for iRun = 1:25
    disp(['run ' num2str(iRun)]);
    load(sprintf('../RandomSplits/split%d',iRun),'posTrainImages');
    % load training positives
    numComponents = min([length(posTrainImages)-1, MnumComponents]);
    dataMatrix = zeros(length(posTrainImages),sx*sy);
    PC = zeros(numComponents,sx*sy);
    for i = 1:length(posTrainImages)
        im = imread(['../posAll/' posTrainImages{i}]);
        if size(im,3) == 3
            im = rgb2gray(im);
        end
        im = double(imresize(im,[sx sy]));
        im = im - mean(im(:));
        im = im / sqrt(var(im(:)));
        dataMatrix(i,:) = reshape(im,1,sx*sy);
    end
    % substracting the mean image
    meanIm = mean(dataMatrix, 1);
    dataMatrix = dataMatrix - repmat(meanIm, size(dataMatrix,1), 1);
    % perform PCA
    [COEFF, SCORE, latent] = princomp(dataMatrix, 'econ');
    PC(1:numComponents,:) = COEFF(:,1:numComponents)';
    
    % testing
    load(sprintf('../RandomSplits/split%d',iRun),'posTestImages','negTestImages');
    
    scorePos = zeros(length(posTestImages),numComponents);
    for i = 1:length(posTestImages)
        im = imread(['../posAll/' posTestImages{i}]);
        if size(im,3) == 3
            im = rgb2gray(im);
        end
        im = double(imresize(im,[sx sy]));
        im = im - mean(im(:));
        im = im / sqrt(var(im(:)));
        im = reshape(im,1,sx*sy);
        % substract mean
        im = im - meanIm;
        % compute the L2 distance
        for j = 1:numComponents
            im = im - sum(im.*PC(j,:)) * PC(j,:);
            score = exp(-sum(im.^2)/sx/sy);
            scorePos(i,j) = score;
        end
        
    end
    
    scoreNeg = zeros(length(negTestImages),numComponents);
    for i = 1:length(negTestImages)
        im = imread(['../negAll/' negTestImages{i}]);
        if size(im,3) == 3
            im = rgb2gray(im);
        end
        im = double(imresize(im,[sx sy]));
        im = im - mean(im(:));
        im = im / sqrt(var(im(:)));
        im = reshape(im,1,sx*sy);
        % substract mean
        im = im - meanIm;
        % compute the L2 distance
        for j = 1:numComponents
            im = im - sum(im.*PC(j,:)) * PC(j,:);
            score = exp(-sum(im.^2)/sx/sy);
            scoreNeg(i,j) = score;
        end

    end

    [ap testAUC tp fp] = evaluate(scorePos,scoreNeg);
    save(sprintf('pca_result_run%d.mat',iRun),...
        'scoreNeg','scorePos','PC','tp','fp','testAUC','ap','meanIm');
end
