% measureRestore.m

clc;
printTitle('Measuring restoration quality ...');

% Loop over all input images
if ~exist('queryFilenames', 'var')
    load(SURFQueryMAT);
end
if ~exist('databaseFilenames', 'var')
    load(SURFDatabaseMAT);
end
if exist('queryImagesGray', 'var')
    clear queryImagesGray;
end
numImages = length(queryFilenames);
allImages = 1:numImages;
MSERGBOrig = NaN * ones(size(allImages));
MSERGB1stRGB = NaN * ones(size(allImages));
MSERGB1stLAB = NaN * ones(size(allImages));
MSELABOrig = NaN * ones(size(allImages));
MSELAB1stRGB = NaN * ones(size(allImages));
MSELAB1stLAB = NaN * ones(size(allImages));
for nImage = allImages
    
    disp(sprintf('\nMeasuring image %d ... \n', nImage));
        
    % Load database image
    nImageDB = queryToDatabaseNumbers(nImage);
    imgColorDB = double( imread(databaseFilenames{nImageDB}) );
    
    % Load projected images
    if SURFEnhanced
        enhancedTag = 'enhanced';
    else
        enhancedTag = 'basic';
    end
    projName = sprintf('%s%.2d-%s-orig-proj.jpg', restoreDir, nImage, enhancedTag);
    projName1stRGB = sprintf('%s%.2d-%s-1stRGB-proj.jpg', restoreDir, nImage, enhancedTag);
    projName1stLAB = sprintf('%s%.2d-%s-1stLAB-proj.jpg', restoreDir, nImage, enhancedTag);
    if ~exist(projName, 'file') || ~exist(projName1stRGB, 'file') || ~exist(projName1stLAB, 'file')
        continue;
    end
    imgColorProj = double( imread(projName) );
    imgColor1stRGBProj = double( imread(projName1stRGB) );
    imgColor1stLABProj = double( imread(projName1stLAB) );
    
    % Crop to center region
    crop = 1;
    if crop
        imgColorDB = imgColorDB(50:450, 50:450, 1:3);
        imgColorProj = imgColorProj(50:450, 50:450, 1:3);
        imgColor1stRGBProj = imgColor1stRGBProj(50:450, 50:450, 1:3);
        imgColor1stLABProj = imgColor1stLABProj(50:450, 50:450, 1:3);
    end
    [height, width, channels] = size(imgColorDB);
    area = height * width;
    
    % Measure MSE in RGB space
    MSERGBOrig(nImage) = sum( (imgColorProj(:) - imgColorDB(:)).^2 ) / area;
    MSERGB1stRGB(nImage) = sum( (imgColor1stRGBProj(:) - imgColorDB(:)).^2 ) / area;
    MSERGB1stLAB(nImage) = sum( (imgColor1stLABProj(:) - imgColorDB(:)).^2 ) / area;
    
    % Measure MSE in LAB space
    imgColorDBFlat = reshape(imgColorDB, area, channels);
    imgColorProjFlat = reshape(imgColorProj, area, channels);
    imgColor1stRGBProjFlat = reshape(imgColor1stRGBProj, area, channels);
    imgColor1stLABProjFlat = reshape(imgColor1stLABProj, area, channels);
    whiteRGB = [255 255 255];
    imgColorDBFlat = rgb2xyz(imgColorDBFlat, whiteRGB);
    imgColorProjFlat = rgb2xyz(imgColorProjFlat, whiteRGB);
    imgColor1stRGBProjFlat = rgb2xyz(imgColor1stRGBProjFlat, whiteRGB);
    [imgColor1stLABProjFlat, whiteXYZ] = rgb2xyz(imgColor1stLABProjFlat, whiteRGB);
    imgColorDBFlat = vcXYZ2lab(imgColorDBFlat, whiteXYZ);
    imgColorProjFlat = vcXYZ2lab(imgColorProjFlat, whiteXYZ);
    imgColor1stRGBProjFlat = vcXYZ2lab(imgColor1stRGBProjFlat, whiteXYZ);
    imgColor1stLABProjFlat = vcXYZ2lab(imgColor1stLABProjFlat, whiteXYZ);
    MSELABOrig(nImage) = sum( (imgColorProjFlat(:) - imgColorDBFlat(:)).^2 ) / area;
    MSELAB1stRGB(nImage) = sum( (imgColor1stRGBProjFlat(:) - imgColorDBFlat(:)).^2 ) / area;
    MSELAB1stLAB(nImage) = sum( (imgColor1stLABProjFlat(:) - imgColorDBFlat(:)).^2 ) / area;
    
    % Clean up memory
    clear imgColor*;
    
end % nImage

figure(1); clf;
subplot(1,2,1);
h = plot(allImages, MSERGBOrig ./ MSERGBOrig, 'b-o', ...
    allImages, MSERGB1stRGB ./ MSERGBOrig, 'r-x');
xlabel('image number'); ylabel('normalized MSE');
title('RGB MSE Metric');
legend('original', '1st order RGB');
axis([1 numImages 0 2]);
set(h(1), 'MarkerSize', 5, 'MarkerFaceColor', 'blue', 'LineWidth', 2);
set(h(2), 'MarkerSize', 7, 'MarkerFaceColor', 'red', 'LineWidth', 2);
subplot(1,2,2);
h = plot(allImages, MSERGBOrig ./ MSERGBOrig, 'b-o', ...
    allImages, MSERGB1stLAB ./ MSERGBOrig, 'r-x');
xlabel('image number'); ylabel('normalized MSE');
title('RGB MSE Metric');
legend('original', '1st order LAB');
axis([1 numImages 0 2]);
set(h(1), 'MarkerSize', 5, 'MarkerFaceColor', 'blue', 'LineWidth', 2);
set(h(2), 'MarkerSize', 7, 'MarkerFaceColor', 'red', 'LineWidth', 2);
disp('Average RGB MSE:');
disp(round( [mean(MSERGBOrig) mean(MSERGB1stRGB) mean(MSERGB1stLAB)] ));

figure(2); clf;
subplot(1,2,1);
h = plot(allImages, MSELABOrig ./ MSELABOrig, 'b-o', ...
    allImages, MSELAB1stRGB ./ MSELABOrig, 'r-x');
xlabel('image number'); ylabel('normalized MSE');
title('LAB MSE Metric');
legend('original', '1st order RGB');
axis([1 numImages 0 2]);
set(h(1), 'MarkerSize', 5, 'MarkerFaceColor', 'blue', 'LineWidth', 2);
set(h(2), 'MarkerSize', 7, 'MarkerFaceColor', 'red', 'LineWidth', 2);
subplot(1,2,2);
h = plot(allImages, MSELABOrig ./ MSELABOrig, 'b-o', ...
    allImages, MSELAB1stLAB ./ MSELABOrig, 'r-x');
xlabel('image number'); ylabel('normalized MSE');
title('LAB MSE Metric');
legend('original', '1st order LAB');
axis([1 numImages 0 2]);
set(h(1), 'MarkerSize', 5, 'MarkerFaceColor', 'blue', 'LineWidth', 2);
set(h(2), 'MarkerSize', 7, 'MarkerFaceColor', 'red', 'LineWidth', 2);
disp('Average RGB MSE:');
disp(round( [mean(MSELABOrig) mean(MSELAB1stRGB) mean(MSELAB1stLAB)] ));

% Save results
save(SURFRestoreMAT, 'MSE*');