function IMG_DDT = DDTriangulate( IMG )
% function IMG_DDT = DDTriangulate( IMG )
% IMG: A three color image matrix sampled by following CFA:
%   R G R G...
%   G B G B...
%   R G R G...
%   ..........
% IMG_DDT: The reconstructed Image

%Ref: Dan Su and Philip Wills, "Image Interpolation by Pixel-Leve
%Data-Dependent Traiangulation", Volume 23(2004, number 2 pp. 189-201,
%Computer Graphic Forum.

% %RGGB Test Image Load
% IMG_ORI = double( imread('Lenna-1.bmp') );
% IMG = double( imread('Lenna-1_RGGB.bmp') );
IMG = double(IMG);
IMG_SIZE = size(IMG(:,:,1));

%Initialization & Parameter declaration
IMG_DDT = zeros( [IMG_SIZE, 3]);
IsYsizeOdd = mod( IMG_SIZE(1), 2);
IsXsizeOdd = mod( IMG_SIZE(2), 2);

%We need 2 types of interpolation:
% R  1  R  According to the paper, interpolation on pixel "1" only depends
% 1  2  1  on the adjacent R pixels; only pixel "2" involves with the
% R  1  R  direction of diagonal.

% R channel and B channel Interpolation on type "1"
conv_win = [0 1 0; 1 2 1; 0 1 0] /2;
IMG_DDT(:,:,1) = conv2(IMG(:,:,1), conv_win, 'same');
IMG_DDT(:,:,2) = IMG(:,:,2);
IMG_DDT(:,:,3) = conv2(IMG(:,:,3), conv_win, 'same');

% R channel and B channel diagonal map generation, 0 means "\" diagonal is chosen, 1
% means "/" diagonal is chosen.
Ycount_R = floor( (IMG_SIZE(1)-1)/2);
Xcount_R = floor( (IMG_SIZE(2)-1)/2);
Rmap = zeros( Ycount_R, Xcount_R );
for y = 1:Ycount_R
    for x = 1:Xcount_R
        Diag_0 = abs( IMG(2*y-1, 2*x-1,1) - IMG(2*y+1, 2*x+1,1) );
        Diag_1 = abs( IMG(2*y-1, 2*x+1,1) - IMG(2*y+1, 2*x-1,1) );
        if Diag_0 < Diag_1
            Rmap(y,x) = 1;
        end
    end
end

Ycount_B = floor( IMG_SIZE(1)/2 ) -1;
Xcount_B = floor( IMG_SIZE(2)/2 ) -1;
Bmap = zeros( Ycount_B, Xcount_B);
for y = 1:Ycount_B
    for x = 1:Xcount_B
        Diag_0 = abs( IMG(2*y, 2*x, 3) - IMG(2*y+2, 2*x+2, 3) );
        Diag_1 = abs( IMG(2*y, 2*x+2, 3) - IMG(2*y+2, 2*x, 3) );
        if Diag_0 < Diag_1
            Bmap(y,x) = 1;
        end
    end
end

%R chanel and B channel Interpolation on type "2"
for y = 1:Ycount_R
    for x = 1:Xcount_R
        if Rmap(y,x) %Choose "/" diagonal
            IMG_DDT(2*y, 2*x,1) = mean( [IMG(2*y-1, 2*x+1,1) IMG(2*y+1, 2*x-1,1)] );
        else %Choose "\" diagonal
            IMG_DDT(2*y, 2*x,1) = mean( [IMG(2*y-1, 2*x-1,1) IMG(2*y+1, 2*x+1,1)] );
        end
    end
end

for y = 1:Ycount_B
    for x = 1:Xcount_B
        if Bmap(y,x) %Choose "/" diagonal
            IMG_DDT(2*y+1, 2*x+1,3) = mean( [IMG(2*y, 2*x+2,3) IMG(2*y+2, 2*x,3)] );
        else %Choose "\" diagonal
            IMG_DDT(2*y+1, 2*x+1,3) = mean( [IMG(2*y, 2*x,3) IMG(2*y+2, 2*x+2,3)] );
        end
    end
end

%G channel Interpolation
for y = 2:IMG_SIZE(1)-1
    IsYOdd = mod(y,2);
    for x = 2+IsYOdd:2:IMG_SIZE(2)-1
        if( abs( diff(IMG(y, x-1:2:x+1, 2)) ) < abs( diff(IMG(y-1:2:y+1, x, 2))) )
            IMG_DDT(y,x,2) = mean( IMG(y, x-1:2:x+1,2) );
        else
            IMG_DDT(y,x,2) = mean( IMG(y-1:2:y+1, x,2) );
        end
    end
end

%Boundary Refinement, R
if ~IsXsizeOdd
    IMG_DDT(:,IMG_SIZE(2),1) = IMG_DDT(:,IMG_SIZE(2)-1,1);
end
if ~IsYsizeOdd
    IMG_DDT(IMG_SIZE(1),:,1) = IMG_DDT(IMG_SIZE(1)-1,:,1);
end

%Boundary Refinement, G
conv_win = [1 2 1]/2;
IMG_DDT(1,:,2) = conv2( IMG(1,:,2), conv_win, 'same');
IMG_DDT(IMG_SIZE(1),:,2) = conv2( IMG(IMG_SIZE(1),:,2), conv_win, 'same');
%IMG_DDT(1,1,2) = ( IMG(1,2,2) + IMG(2,1,2) )/2;
IMG_DDT(:,1,2) = conv2( IMG(:,1,2), conv_win', 'same');
IMG_DDT(:,IMG_SIZE(2),2) = conv2( IMG(:,IMG_SIZE(2),2), conv_win', 'same');

%Boundary Refinement, B
IMG_DDT(:,1,3) = IMG_DDT(:,2,3);
IMG_DDT(1,:,3) = IMG_DDT(2,:,3);
if IsXsizeOdd
    IMG_DDT(:,IMG_SIZE(2),3) = IMG_DDT(:,IMG_SIZE(2)-1,3);
end
if IsYsizeOdd
    IMG_DDT(IMG_SIZE(1),:,3) = IMG_DDT(IMG_SIZE(1)-1,:,3);
end

IMG_DDT = uint8(IMG_DDT);