function [R,G,B]=colorInterp4(r,g,b);
% [R,G,B]=colorInterp4(r,g,b);
%
% Written by Bao-Jun Jiang
%
% NOTE THAT if we don't worry about the one or two rows and
% columns of pixels on the borders of the images, the code
% can be shrunk by a factor of 6.  But as it is now, I have
% discussed all pixels on the borders and considered the
% different cases of even and odd rows and columns.
%
% This function interpolates the colors r, g, and b for the
% bayer pattern (as created by the command bayer), and returns
% the [R,G,B] matrices with the missing color values filled
% by using the block matching algorithm described in Tsai's paper.
% The bayer pattern is assumed to be that of the following
%
%       GRGRGRGR...
%       BGBGBGBG...
%       GRGRGRGR...
%       BGBGBGBG...
%


R=r;
G=g;
B=b;
[Row, Col]=size(R);

%===============================================
%%for pixels that contain blue information only
%===============================================

%considering the first and last columns of blue pixels
for i=1:round(Row/2-1),
   R(2*i,1)=mean([R(2*i-1,2) R(2*i+1,2)]);
   G(2*i,1)=median([G(2*i-1,1) G(2*i,2) G(2*i+1,1)]);
   if Col/2==round(Col/2) %Col is even
       R(2*i,Col-1)=median([R(2*i-1,Col-2) R(2*i-1,Col)...
		R(2*i+1,Col-2) R(2*i+1,Col)]);
       G(2*i,Col-1)=median([G(2*i-1,Col-1) G(2*i+1,Col-1)...
                              G(2*i,Col-2) G(2*i,Col)]);
       else  %Col is odd
           R(2*i,Col)=median([R(2*i-1,Col-1) R(2*i+1,Col-1)]);
           G(2*i,Col)=median([G(2*i-1,Col) G(2*i+1,Col)...
                                G(2*i,Col-1)]);
   end;
end;
if round(Row/2)==Row/2
   R(Row,1)=R(Row-1,2);
   G(Row,1)=mean([G(Row-1,1) G(Row,2)]);
   if round(Col/2)==Col/2  %Col is even
       R(Row,Col-1)=median([R(Row-1,Col-2) R(Row-1,Col)]);
       G(Row,Col-1)=median([G(Row-1,Col-1) G(Row,Col-2) ...
            G(Row, Col)]);
   else %Col is odd
       R(Row, Col)=R(Row-1, Col-1);
       G(Row, Col)=median([G(Row-1,Col) G(Row, Col-1)]);
   end;
end;

%considering the first and last rows of blue pixels
for j=1:floor(Col/2)-1,
   R(2,2*j+1)=median([R(1,2*j) R(1,2*j+2) R(3,2*j) R(3,2*j+2)]);
   G(2,2*j+1)=median([G(1,2*j+1) G(3,2*j+1) G(2,2*j) G(2,2*j+2)]);
end;

for j=1:floor(Col/2)-1,
  if floor(Row/2)==Row/2 %Row is even
     R(Row,2*j+1)=median([R(Row-1,2*j) R(Row-1,2*j+2)]);
     G(Row,2*j+1)=median([G(Row-1,2*j+1) G(Row,2*j)...
                                G(Row,2*j+2)]);
     else %Row is odd
          R(Row-1,2*j+1)=median([R(Row,2*j) R(Row,2*j+2)... 
			R(Row-2,2*j) R(Row-2,2*j+2)]);
          G(Row-1,2*j+1)=median([G(Row-2,2*j+1) G(Row,2*j+1)...
                                G(Row-1,2*j) G(Row-1,2*j+2)]);
   end;
end;  

%considering all the rest of the blue pixels
for i=2:floor(Row/2)-1,
	for j=1:round(Col/2)-2,
       G(2*i,2*j+1)=median([G(2*i-1,2*j+1) G(2*i+1,2*j+1)...
                    G(2*i,2*j) G(2*i,2*j+2)]);  %Green is usual.
        %the comparison has been simplifies mathematically below
        d1=abs(G(2*i-2,2*j)+G(2*i-1,2*j-1)-G(2*i+1, 2*j+1)-...
                     G(2*i,2*j+2));
        d2=abs(G(2*i-2,2*j+2)+G(2*i-1,2*j+3)-G(2*i+1, 2*j+1)-...
                     G(2*i,2*j));
        d3=abs(G(2*i+2,2*j+2)+G(2*i+1,2*j+3)-G(2*i-1, 2*j+1)-...
                     G(2*i,2*j)); 
        d4=abs(G(2*i+2,2*j)+G(2*i+1,2*j-1)-G(2*i-1, 2*j+1)-...
                     G(2*i,2*j+2)); 
        dmin=min([d1 d2 d3 d4]);
        if d1==dmin R(2*i,2*j+1)=R(2*i-1,2*j);
            elseif d2==dmin R(2*i,2*j+1)=R(2*i-1,2*j+2);
            elseif d3==dmin R(2*i,2*j+1)=R(2*i+1,2*j+2);
            else  R(2*i,2*j+1)=R(2*i+1,2*j);
        end;
  end;
end;


%===============================================
%%for pixels that contain red information only
%===============================================

%considering the first row
for j=1:round(Col/2-1),
   G(1,2*j)=median([G(1,2*j-1) G(1,2*j+1) G(2,2*j)]);
   B(1,2*j)=median([B(2,2*j-1) B(2,2*j+1)]);
end;
if floor(Col/2)==Col/2
   G(1,Col)=median([G(1,Col-1) G(2,Col)]);
   B(1,Col)=B(2,Col-1);
end;

% considering the last row of red pixels
if round(Row/2)==(Row+1)/2  %ie, Row is odd
   for j=1:round(Col/2-1),
      G(Row,2*j)=median([G(Row-1,2*j) G(Row,2*j-1) G(Row,2*j+1)]);
      B(Row,2*j)=median([B(Row-1,2*j-1) B(Row-1,2*j+1)]);
   end;
   if floor(Col/2)==Col/2
      G(Row,Col)=median([G(Row-1,Col) G(Row,Col-1)]);
      B(Row,Col)=B(Row-1,Col-1);
    end;

    else % Row is even
        for j=1:round(Col/2-1),
              G(Row-1,2*j)=median([G(Row-2,2*j) G(Row,2*j)...
                             G(Row-1,2*j+1) G(Row-1,2*j-1)]);
              B(Row-1,2*j)=median([B(Row-2,2*j-1) B(Row-2,2*j+1)...
                             B(Row,2*j-1) B(Row,2*j+1)]);
        end;

        if floor(Col/2)==Col/2  %ie,if Col is even
            G(Row-1,Col)=median([G(Row-2,Col) G(Row-1,Col-1) G(Row,Col)]);
            B(Row-1,Col)=median([B(Row-2,Col-1) B(Row,Col-1)]);
        end;
end;


% consideringthe first and last columns of red pixels
% skipping those in either the first or last red rows 

for i=1:round(Row/2)-2,
    G(2*i+1,2)=median([G(2*i,2) G(2*i+1,1) G(2*i+1,3) G(2*i+2,2)]);
    B(2*i+1,2)=median([B(2*i,1) B(2*i,3) B(2*i+2,1) B(2*i+2,3)]);

    if floor(Col/2)==Col/2  %ie, if Col is even
       G(2*i+1,Col)=median([G(2*i,Col) G(2*i+1,Col-1) G(2*i+2,Col)]);
       B(2*i+1,Col)=median([B(2*i,Col-1) B(2*i+2,Col-1)]);
    else %Col is odd
       G(2*i+1,Col-1)=median([G(2*i,Col-1) G(2*i+1,Col-2)...
                               G(2*i+1,Col) G(2*i+2,Col-1)]);
       B(2*i+1,Col-1)=median([B(2*i,Col-2) B(2*i,Col)...
                            B(2*i+2,Col-2) B(2*i+2,Col)]);
    end;
end;

%considering all the rest of the red pixels
for i=1:round(Row/2)-2,
	for j=2:floor(Col/2)-1,
       G(2*i+1,2*j)=median([G(2*i+1,2*j-1) G(2*i+1,2*j+1)...
                   G(2*i,2*j) G(2*i+2,2*j)]);  %Green is usual.
        %the comparison has been simplifies mathematically below
        d1=abs(G(2*i-1,2*j-1)+G(2*i,2*j-2)-G(2*i+1, 2*j+1)-...
                     G(2*i+2,2*j));
        d2=abs(G(2*i-1,2*j+1)+G(2*i,2*j+2)-G(2*i+1, 2*j-1)-...
                     G(2*i+2,2*j));
        d3=abs(G(2*i+2,2*j+2)+G(2*i+3,2*j+1)-G(2*i, 2*j)-...
                     G(2*i+1,2*j-1)); 
        d4=abs(G(2*i+3,2*j-1)+G(2*i+2,2*j-2)-G(2*i, 2*j)-...
                     G(2*i+1,2*j+1)); 
        dmin=min([d1 d2 d3 d4]);
        if d1==dmin   B(2*i+1,2*j)=B(2*i,2*j-1);
            elseif d2==dmin B(2*i+1,2*j)=B(2*i,2*j+1);
            elseif d3==dmin B(2*i+1,2*j)=B(2*i+2,2*j+1);
            else  B(2*i+1,2*j)=B(2*i+2,2*j-1);
        end;
  end;
end;



%===============================================
%%for pixels that contain green information only
%===============================================

%first row
R(1,1)=R(1,2);
B(1,1)=B(2,1);
for j=1:floor(Col/2-1),
   R(1,2*j+1)=mean([R(1,2*j) R(1,2*j+2)]);
   B(1,2*j+1)=B(2,2*j+1);
end;
if floor(Col/2)==(Col-1)/2  %ie, Col is odd
   R(1,Col)=R(1,Col-1);
   B(1,Col)=B(2,Col);
end;

%other odd rows (except possibly the last row)

for i=1:round(Row/2)-2,
   R(2*i+1,1)=R(2*i+1,2);  
   %for the first column we use bilinear method to calculate blue
   B(2*i+1,1)=mean([B(2*i,1) B(2*i+2,1)]);

   for j=1:round(Col/2-2),
      d1=abs(G(2*i+1,2*j+1)+G(2*i+1,2*j-1)-G(2*i,2*j+2)...
                    -G(2*i+2,2*j+2));
      d2=abs(G(2*i+1,2*j+1)+G(2*i+1,2*j+3)-G(2*i,2*j)...
                    -G(2*i+2,2*j));
      if d1<d2  
              R(2*i+1,2*j+1)=R(2*i+1,2*j);
         else R(2*i+1,2*j+1)=R(2*i+1,2*j+2);
      end;

      d3=abs(G(2*i+1,2*j+1)+G(2*i-1,2*j+1)-G(2*i+2,2*j)...
                    -G(2*i+2,2*j+2));
      d4=abs(G(2*i+1,2*j+1)+G(2*i+3,2*j+1)-G(2*i,2*j)...
                    -G(2*i,2*j+2));
      if d3<d4  
              B(2*i+1,2*j+1)=B(2*i,2*j+1);
         else B(2*i+1,2*j+1)=B(2*i+2,2*j+1);
      end;
   end;

   if floor(Col/2)==(Col-1)/2  %ie, Col is odd
       R(2*i+1,Col)=R(2*i+1,Col-1);
       B(2*i+1,Col)=mean([B(2*i,Col) B(2*i+2,Col)]);
      else  %Col is even
       R(2*i+1,Col-1)=mean([R(2*i+1,Col-2) R(2*i+1,Col)]);
       B(2*i+1,Col-1)=mean([B(2*i,Col-1) B(2*i+2,Col-1)]);
   end;
end;

%last odd Row
if floor(Row/2)==(Row-1)/2  %ie., Row is odd
   R(Row,1)=R(Row,2);
   B(Row,1)=B(Row-1,1);
   for j=1:floor(Col/2-1),
      R(Row,2*j+1)=mean([R(Row,2*j) R(Row,2*j+2)]);
      B(Row,2*j+1)=B(Row-1,2*j+1);
   end;
   if floor(Col/2)==(Col-1)/2  %ie, Col is odd
      R(Row,Col)=R(Row,Col-1);
      B(Row,Col)=B(Row-1,Col);
   end;
else %Row is even
    R(Row-1,1)=R(Row-1,2);
    B(Row-1,1)=mean([B(Row-2,1) B(Row,1)]);
    for j=1:floor(Col/2-1),
          R(Row-1,2*j+1)=mean([R(Row-1,2*j) R(Row-1,2*j+2)]);
          B(Row-1,2*j+1)=mean([B(Row-2,2*j+1) B(Row,2*j+1)]);
     end;
     if floor(Col/2)==(Col-1)/2  %ie, Col is odd
         R(Row-1,Col)=R(Row-1,Col-1);
         B(Row-1,Col)=mean([B(Row-2,Col) B(Row,Col)]);
      end;
end;


%%even rows
% first even column.
R(2,2)=mean([R(1,2) R(3,2)]);
B(2,2)=mean([B(2,1) B(2,3)]);

for i=2:floor(Row/2-1),
    B(2*i,2)=mean([B(2*i,1) B(2*i,3)]);
    d1=abs(G(2*i,2)+G(2*i-2,2)-G(2*i+1,1)-G(2*i+1,3));
    d2=abs(G(2*i,2)+G(2*i+2,2)-G(2*i-1,1)-G(2*i-1,3));
    if d1<d2  R(2*i,2)=R(2*i-1,2);
        else R(2*i,2)=R(2*i+1,2);
    end;
end;

if floor(Row/2)==Row/2  %Row is even
    R(Row,2)=R(Row-1,2);
    B(Row,2)=mean([B(Row,1) B(Row,3)]);
   else  %Row is odd
    R(Row-1,2)=mean([R(Row-2,2) R(Row,2)]);
    B(Row-1,2)=mean([B(Row-1,1) B(Row-1,3)]);
end;

% last even column
if floor(Col/2)==Col/2 % Col is even
    for i=1:floor(Row/2-1),
        R(2*i,Col)=mean([R(2*i-1,Col) R(2*i+1,Col)]);
        B(2*i,Col)=B(2*i,Col-1);
    end;
    if floor(Row/2)==Row/2
        R(Row, Col)=R(Row-1,Col);
        B(Row, Col)=B(Row, Col-1);
    end;   

  else  %Col is odd
     R(2,Col-1)=mean([R(1,Col-1) R(3,Col-1)]);
     B(2,Col-1)=mean([B(2,Col-2) B(2,Col)]);

    for i=2:floor(Row/2-1),
        B(2*i,Col-1)=mean([B(2*i,Col-2) B(2*i,Col)]);
        d1=abs(G(2*i,Col-1)+G(2*i-2,Col-1)...
                      -G(2*i+1,Col-2)-G(2*i+1,Col));
        d2=abs(G(2*i,Col-1)+G(2*i+2,Col-1)...
                    -G(2*i-1,Col-2)-G(2*i-1,Col));
        if d1<d2  R(2*i,Col-1)=R(2*i-1,Col-1);
            else R(2*i,Col-1)=R(2*i+1,Col-1);
        end;
    end;

   if floor(Row/2)==Row/2  %Row is even
       R(Row,Col-1)=R(Row-1,Col-1);
       B(Row,Col-1)=mean([B(Row,Col-2) B(Row,Col)]);
    else  %Row is odd
       R(Row-1,Col-1)=mean([R(Row-2,Col-1) R(Row,Col-1)]);
       B(Row-1,Col-1)=mean([B(Row-1,Col-2) B(Row-1,Col)]);
   end;
end;

% first even row
for j=2:floor(Col/2-1),
    R(2,2*j)=mean([R(1,2*j) R(3,2*j)]);
    d1=abs(G(2,2*j)+G(2,2*j-2)-G(1,2*j+1)-G(3,2*j+1));
    d2=abs(G(2,2*j)+G(2,2*j+2)-G(1,2*j-1)-G(3,2*j-1));
    if d1<d2  B(2,2*j)=B(2,2*j-1);
            else B(2,2*j)=B(2,2*j+1);
    end;
end;

%last even row
if floor(Row/2)==Row/2,  %Row is even
    for j=2:round(Col/2-1),
          R(Row, 2*j)=R(Row-1, 2*j);
          B(Row, 2*j)=mean([B(Row,2*j-1) B(Row, 2*j+1)]);
    end;
   if round(Col/2)==Col/2 %Col is even
         R(Row,Col)=R(Row-1,Col);
         B(Row,Col)=B(Row, Col-1);
   end;

else %Row is odd
   for j=2:floor(Col/2-1),
        R(Row-1,2*j)=mean([R(Row-2,2*j) R(Row,2*j)]);
        d1=abs(G(Row-1,2*j)+G(Row-1,2*j-2)...
                     -G(Row-2,2*j+1)-G(Row,2*j+1));
        d2=abs(G(Row-1,2*j)+G(Row-1,2*j+2)...
                     -G(Row-2,2*j-1)-G(Row,2*j-1));
        if d1<d2  B(Row-1,2*j)=B(Row-1,2*j-1);
            else B(Row-1,2*j)=B(Row-1,2*j+1);
        end;
   end;
end;

%the rest of the even rows
for i=2:floor(Row/2-1),
    for j=2:floor(Col/2-1),
      d1=abs(G(2*i,2*j)+G(2*i-2,2*j)-G(2*i+1,2*j-1)...
                    -G(2*i+1,2*j+1));
      d2=abs(G(2*i,2*j)+G(2*i+2,2*j)-G(2*i-1,2*j-1)...
                    -G(2*i-1,2*j+1));
      if d1<d2  
              R(2*i,2*j)=R(2*i-1,2*j);
         else R(2*i,2*j)=R(2*i+1,2*j);
      end;

      d3=abs(G(2*i,2*j)+G(2*i,2*j-2)-G(2*i-1,2*j+1)...
                    -G(2*i+1,2*j+1));
      d4=abs(G(2*i,2*j)+G(2*i,2*j+2)-G(2*i-1,2*j-1)...
                    -G(2*i+1,2*j-1));
      if d3<d4  
              B(2*i,2*j)=B(2*i,2*j-1);
         else B(2*i,2*j)=B(2*i,2*j+1);
      end;
   end;
end;













