function IMG_TwoLine = TwoLineInterpolate( IMG )

Te = 200;
Tcr = 50;
IMG_SIZE = size(IMG(:,:,1));
tmp = zeros( [IMG_SIZE, 3] );
%    x_end_minus1 = IMG_SIZE(2)-1;
IsXsizeOdd = mod( IMG_SIZE(2), 2);
conv_win = [1 2 1]/2;

%Green Channel, line 1
for x = 2:IMG_SIZE(2)-1
    tmp(1,:,2) = conv2( IMG(1,:,2), conv_win, 'same' );
end
tmp(1,1,2) = ( IMG(2,1,2) + IMG(1,2,2) )/2; % rewrite Green at 1,1
if IsXsizeOdd % rewrite the last Green if Xsize is odd.
    tmp(1,IMG_SIZE(2), 2) = ( IMG(2,IMG_SIZE(2),2) + IMG(1,IMG_SIZE(2)-1,2) )/2;
end

%Green Channel, the rest lines.
tmp(2:IMG_SIZE(1),:,2) = IMG(2:IMG_SIZE(1),:,2);
for y = 2:IMG_SIZE(1)
    IsYodd = mod( y, 2);
    if IsYodd   %for the 1st pixel in each line.
        tmp(y, 1, 2) = Interpolate2L( tmp(y-1,1,2), tmp(y,2,2), tmp(y-1,2,2) );
    end
    for x = 2+IsYodd:2:IMG_SIZE(2)-1
        tmp(y, x, 2) = ( Interpolate2L( tmp(y-1,x,2),   tmp(y,x+1,2),   tmp(y-1,x+1,2) ) + Interpolate2L( tmp(y,x-1,2), tmp(y-1,x,2), tmp(y-1,x-1,2) ) )/2;
    end
    if ~IsXsizeOdd
        tmp(y, IMG_SIZE(2),2) = Interpolate2L( IMG(y,IMG_SIZE(2)-1,2), IMG(y-1,IMG_SIZE(2),2), IMG(y-1,IMG_SIZE(2)-1,2) );
    end
end

%Red Channel, 1st line.
for x = 2:IMG_SIZE(2)-1
    tmp(1,:,1) = conv2( IMG(1,:,1), conv_win, 'same' );
end
if ~IsXsizeOdd % rewrite the last Red if Xsize is even.
    tmp(1,IMG_SIZE(2), 1) = ( IMG(2,IMG_SIZE(2),1) + IMG(1,IMG_SIZE(2)-1,1) )/2;
end

%Red Channel, the rest lines
tmp(2:IMG_SIZE(1),:,1) = IMG(2:IMG_SIZE(1),:,1);
%There are three cases for R channel:
% R   R  R
%   3 2
% R 1 R  R

%Case 1:
for y =3:2:IMG_SIZE(1)
    for x = 2:2:IMG_SIZE(2)-1
        GradG = abs( tmp(y,x-1,2) - tmp(y,x+1,2) )/2;
        factor1 = Interpolate2Lfactor( tmp(y,x+1,2), tmp(y,x,2) );
        factor2 = Interpolate2Lfactor( tmp(y,x-1,2), tmp(y,x,2) );
        if( GradG < Te & tmp(y,x-1,2) > Tcr & tmp(y,x+1,2) > Tcr )
            tmp(y,x,1) = tmp(y,x,2) * ( factor1 * tmp(y,x+1,1) / tmp(y,x+1,2) + factor2 * tmp(y,x-1,1) / tmp(y,x-1,2) ) / (factor1 + factor2);
        else
            tmp(y,x,1) = ( factor1 * tmp(y,x+1,1) + factor2 * tmp(y,x-1,1) ) / (factor1 + factor2);
        end
    end
    % The last pixel
    if ~IsXsizeOdd
        tmp(y,IMG_SIZE(2),1) = tmp(y,IMG_SIZE(2)-1,1);
    end
end

%Case 3:
for y = 2:2:IMG_SIZE(1)
    tmp(y,1,1) = tmp(y-1,1,1);
    for x = 2:2:IMG_SIZE(2)-1
        GradG = ( (tmp(y-1,x,2)-tmp(y,x,2))^2 + (tmp(y,x-1,2)-tmp(y,x,2))^2 )^.5;    
        factor1 = Interpolate2Lfactor(tmp(y-1,x-1,2), tmp(y,x,2) );
        factor2 = Interpolate2Lfactor(tmp(y-1,x+1,2), tmp(y,x,2) );
        if( GradG < Te & tmp(y-1,x-1,2) > Tcr & tmp(y-1, x+1,2) > Tcr )
            tmp(y,x,1) = (1+tmp(y,x,2)) * ( factor1 * tmp(y-1,x-1,1) / (1+tmp(y-1,x-1,2)) + factor2 * tmp(y-1,x+1,1) / (1+tmp(y-1,x+1,2)) ) / (factor1+factor2);
        else
            tmp(y,x,1) = ( tmp(y-1,x-1,1) * factor1 + tmp(y-1,x+1,1) *factor2 ) / (factor1+factor2);
        end
    end
    if ~IsXsizeOdd
        tmp(y,IMG_SIZE(2),1) = tmp(y-1, IMG_SIZE(2),1);
    end
end

%Case 2:
for y = 2:2:IMG_SIZE(1)
    for x = 3:2:IMG_SIZE(2)
        GradG = abs( tmp(y-1,x,2) - tmp(y,x,2) )/2;
        if( GradG < Te & tmp(y-1,x,2) > Tcr)
            tmp(y,x,1) = tmp(y,x,2) * tmp(y-1,x,1) / tmp(y-1,x,2);
        else
            tmp(y,x,1) = Interpolate2L(tmp(y,x-1,1), tmp(y-1,x,1), tmp(y-1,x-1,1) ) ;
        end
    end
end

%Blue Channel, 1st line.
tmp(:,:,3) = IMG(:,:,3);
for x = 1:2:IMG_SIZE(2)-1
    tmp(1,x,3) = tmp(2,x+1,3);
    tmp(1,x+1,3) = tmp(2,x+1,3);
end
if ~IsXsizeOdd
    tmp(1, IMG_SIZE(2),3) = tmp(2, IMG_SIZE(2)-1,3);
end

%Blue Channel, the rest lines,
%Case 1:
for y =2:2:IMG_SIZE(1)
    tmp(y,1,3) = tmp(y,2,3);
    for x = 3:2:IMG_SIZE(2)-1
        GradG = abs( tmp(y,x-1,2) - tmp(y,x+1,2) )/2;
        factor1 = Interpolate2Lfactor( tmp(y,x+1,2), tmp(y,x,2) );
        factor2 = Interpolate2Lfactor( tmp(y,x-1,2), tmp(y,x,2) );
        if( GradG < Te & tmp(y,x-1,2) > Tcr & tmp(y,x+1,2) > Tcr )
            tmp(y,x,3) = tmp(y,x,2) * ( factor1 * tmp(y,x+1,3) / tmp(y,x+1,2) + factor2 * tmp(y,x-1,3) / tmp(y,x-1,2) ) / (factor1 + factor2);
        else
            tmp(y,x,3) = ( factor1 * tmp(y,x+1,3) + factor2 * tmp(y,x-1,3) ) / (factor1 + factor2);
        end
    end
    % The last pixel
    if IsXsizeOdd
        tmp(y,IMG_SIZE(2),3) = tmp(y,IMG_SIZE(2)-1,3);
    end
end

%Case 3:
for y = 3:2:IMG_SIZE(1)
    tmp(y,1,3) = tmp(y-1,1,3);
    for x = 3:2:IMG_SIZE(2)-1
        GradG = ( (tmp(y-1,x,2)-tmp(y,x,2))^2 + (tmp(y,x-1,2)-tmp(y,x,2))^2 )^.5;    
        factor1 = Interpolate2Lfactor(tmp(y-1,x-1,2), tmp(y,x,2) );
        factor2 = Interpolate2Lfactor(tmp(y-1,x+1,2), tmp(y,x,2) );
        if( GradG < Te & tmp(y-1,x-1,2) > Tcr & tmp(y-1, x+1,2) > Tcr )
            tmp(y,x,3) = (1+tmp(y,x,2)) * ( factor1 * tmp(y-1,x-1,3) / (1+tmp(y-1,x-1,2)) + factor2 * tmp(y-1,x+1,3) / (1+tmp(y-1,x+1,2)) ) / (factor1+factor2);
        else
            tmp(y,x,3) = ( tmp(y-1,x-1,3) * factor1 + tmp(y-1,x+1,3) *factor2 ) / (factor1+factor2);
        end
    end
    if IsXsizeOdd
        tmp(y,IMG_SIZE(2),3) = tmp(y-1, IMG_SIZE(2),3);
    end
end

%Case 2:
for y = 3:2:IMG_SIZE(1)
    for x = 2:2:IMG_SIZE(2)
        GradG = abs( tmp(y-1,x,2) - tmp(y,x,2) )/2;
        if( GradG < Te & tmp(y-1,x,2) > Tcr)
            tmp(y,x,3) = tmp(y,x,2) * tmp(y-1,x,3) / tmp(y-1,x,2);
        else
            tmp(y,x,3) = Interpolate2L(tmp(y,x-1,3), tmp(y-1,x,3), tmp(y-1,x-1,3) );
        end
    end
end

IMG_TwoLine = uint8(tmp);