clear all; close all; clc;
for img = 5:-1:1

% In this routine, we assume the CFA image is of the format RGGB, where the
% 1st pixel in the 1st line is R.
Bilinear = 1;
FFT_TimeDomainRECOVER = 0;
FFT_FreqDomainRECOVER = 0;
TwoLine = 1;
DDT = 1;
ED = 1;
HQL = 1;
ERROR_METHOD = 'PSNR';

%RGGB Test Image Load
%filename = '.\Test Images\Mackay-GoodOptics_ori';
switch img
    case 1
filename = '.\Test Images\FREQ-ORIENT_ori';
filefmt = 'tif';
    case 2
filename = '.\Test Images\pointArrayPoorOptics_ori';
filefmt = 'tif';
    case 3
filename = '.\Test Images\Mackay-PoorOptics_ori';
filefmt = 'tif';        
    case 4
filename = '.\Test Images\peppers';
filefmt = 'tif';
    case 5
filename = '.\Test Images\baboon';
filefmt = 'tif';
end
filename_fmt = strcat( filename, '.', filefmt );
IMG_ORI = double( imread( filename_fmt ) );
[tmp, IMG] = RGB2CFA(IMG_ORI);
clear tmp;
IMG = double(IMG);
IMG_SIZE = size(IMG(:,:,1));
filefmt = 'tif';

%Bilinear Interpolation
if Bilinear
    IMG_Bilinear = BilinearInterpolate( IMG );
    figure; imshow(IMG_Bilinear);
    imwrite(IMG_Bilinear, strcat(filename, '_bilinear.', filefmt), filefmt);
    Error_Bilinear = ImgComp( IMG_Bilinear, IMG_ORI, ERROR_METHOD )
end

%2D FFT Interpolation
if ( FFT_TimeDomainRECOVER ==1 | FFT_FreqDomainRECOVER == 1)
	barX_halfsize = floor(IMG_SIZE(2)/2);
	barY_halfsize = floor(IMG_SIZE(1)/2);
	if (FFT_TimeDomainRECOVER ==1)
% 		sinc_barX = sinc( 1/2*[-barX_halfsize:barX_halfsize ]);
% 		sinc_barY = sinc( 1/2*[-barY_halfsize:barY_halfsize ]);
		sinc_barX = sinc( 1/2*[-124:124 ]);
		sinc_barY = sinc( 1/2*[-124:124 ]);
		sinc_win = sinc_barX' * sinc_barY;
		IMG_upsample = zeros( [IMG_SIZE, 3] );
		IMG_upsample(:,:,1) = conv2(IMG(:,:,1), sinc_win, 'same');
		IMG_upsample(:,:,2) = conv2(IMG(:,:,2), sinc_win, 'same');
		IMG_upsample(:,:,3) = conv2(IMG(:,:,3), sinc_win, 'same');
    else
		IMG_f = fft2(IMG);
		IMG_upsample_f = zeros( [IMG_SIZE, 3] );
		IMG_upsample_f(1:barY_halfsize, 1:barX_halfsize, :) = 4*IMG_f(1:barY_halfsize, 1:barX_halfsize, :);
		IMG_upsample_f(:,:,2) = IMG_upsample_f(:,:,2)/2;
		IMG_upsample = ifft2(IMG_upsample_f);
    end
    
	IMG_2DFFTRECOVER = uint8( abs(IMG_upsample( 1:IMG_SIZE(1),  1:IMG_SIZE(2), :)) );
	figure; imshow(IMG_2DFFTRECOVER);
    ERROR_2DFFT = ImgComp(IMG_2DFFTRECOVER, IMG_ORI, ERROR_METHOD);
    
	if ~DEBUG
        clear IMG_upsample_f IMG_upsample IMG_f;
	end
    
end

%Tow Line Interpolation [Ref: Color demosaicing with constrained buffering]
%Functions: Interpolate2L, Interpolate2Lfactor
if TwoLine
    IMG_TwoLine = TwoLineInterpolate( IMG );
    figure;    imshow(IMG_TwoLine);
    imwrite(IMG_TwoLine, strcat(filename, '_2line.', filefmt), filefmt);
    ERROR_TwoLine = ImgComp(IMG_TwoLine, IMG_ORI, ERROR_METHOD)
end

% Data-Dependent Triangulation Method
if DDT
	IMG_DDT = DDTriangulate( IMG );
	figure; imshow(IMG_DDT);
    imwrite(IMG_DDT, strcat(filename, '_DDT.', filefmt), filefmt);
	Error_DDT = ImgComp( IMG_DDT, IMG_ORI, ERROR_METHOD )
end

%Edge Directed Interpolation.
if ED
    [IMG_ED EdgDet_map] = EdgDirInterpolate( IMG );
    figure; imshow(IMG_ED);
    imwrite(IMG_ED, strcat(filename, '_ED.', filefmt), filefmt);
    Error_ED = ImgComp( IMG_ED, IMG_ORI, ERROR_METHOD )
end

if HQL
    IMG_HQL = HiQualLinearInterpolate( IMG );
    figure; imshow(IMG_HQL);
    imwrite(IMG_HQL, strcat(filename, '_HQL.', filefmt), filefmt);
    ERROR_HQL = ImgComp( IMG_HQL, IMG_ORI, 'psnr')
end

end