function IMG_HQL = HiQualLinearInterpolate( IMG );

IMG = double(IMG);
[size_y size_x color] = size(IMG);

%Empirical Gain factor
alpha = 1/2;
beta = 5/8;
gamma = 3/4;

%Delta R/G/B evaluation filter banks.
filter1 = [ 0 0 -1 0 0;
            0 0  0 0 0;
           -1 0  4 0 -1;
            0 0  0 0 0;
            0 0 -1 0 0]/4;

filter2 = [ 0  0 .5  0 0;
            0 -1  0 -1 0;
            -1 0  5  0 -1;
            0 -1  0 -1 0;
            0  0 .5  0 0]/5;
filter3 = filter2';

filter4 = [ 0  0 -3/2 0  0;
            0  0   0  0  0;
          -3/2 0   6  0 -3/2;
            0  0   0  0  0;
            0  0 -3/2 0  0]/6;

IMG_tmp = IMG;
IMG_CFA = sum(IMG,3);

for y = 3:size_y-2
    IsYOdd = mod(y,2);
    for x = 3:size_x-2
        IsXOdd = mod(x,2);
        
        Neighbors1 = [ IMG_CFA(y-1,x) IMG_CFA(y,x-1) IMG_CFA(y,x+1) IMG_CFA(y+1,x) ]; % adjacent neighbors
        Neighbors2 = [ IMG_CFA(y-1,x-1) IMG_CFA(y-1,x+1)  IMG_CFA(y+1,x-1)  IMG_CFA(y+1,x+1) ]; % diagonal neighbors
        Neighbors3 = [ IMG_CFA(y-1,x)  IMG_CFA(y+1,x) ]; % vertical neighbors
        Neighbors4 = [ IMG_CFA(y,x-1)  IMG_CFA(y,x+1) ]; % horizontal neighbors
        Local_CFA = IMG_CFA(y-2:y+2, x-2:x+2);
        
        switch ( 2*IsYOdd + IsXOdd)
            case 3 % at R pixel, Interpolate G & B, deltaG uses filter1, deltaB uses filter4
                IMG_tmp(y,x,2) = mean( Neighbors1 ) + beta* sum(sum( Local_CFA  .* filter1 ));
                IMG_tmp(y,x,3) = mean( Neighbors2 ) + gamma* sum(sum( Local_CFA .* filter4 ));
            case 2 % at Gr pixel, Interpolate R & B, deltaR uses filter2, deltaB uses filter3
                IMG_tmp(y,x,1) = mean( Neighbors4 ) + alpha* sum(sum( Local_CFA .* filter2 ));
                IMG_tmp(y,x,3) = mean( Neighbors3 ) + gamma* sum(sum( Local_CFA .* filter3 ));
            case 1 % at Gb pixel, Interpolate R & B, deltaR uses filter3, deltaB uses filter2
                IMG_tmp(y,x,1) = mean( Neighbors3 ) + alpha* sum(sum( Local_CFA .* filter3 ));
                IMG_tmp(y,x,3) = mean( Neighbors4 ) + gamma* sum(sum( Local_CFA .* filter2 ));
            case 0 % at B pixel, Interpolate R & G, deltaR uses filter4, deltaG uses filter1
                IMG_tmp(y,x,1) = mean( Neighbors2 ) + alpha* sum(sum( Local_CFA .* filter4 ));
                IMG_tmp(y,x,2) = mean( Neighbors1 ) + beta* sum(sum( Local_CFA .* filter1 ));
        end
    end
end

IMG_HQL = uint8(IMG_tmp);