% bfilter (Bilateral Filter) by Dan Piehl (2006)
%
% This function used the bilateral filter (Tomasi & Manduchi 1998) to filter
% a one or two dimensional signal.  A third dimension, if supplied, will
% function as color channels, with "range distances" determined by
% a Euclidean distance measurement in color space.  This function requires two
% sigma values (domain and range) which determine the strength of the filter.
% (eg. for an average color photo start with sdomain ~ 4 and srange ~ 0.2)

function b = bfilter(a,sdomain,srange)

  % Examine parameters and return error messages, if needed
  numDimsA = ndims(a);          % Get the number of input dimensions
  S1=size(a,1);                 % Get the number of rows
  S2=size(a,2);                 % Get the number of columns
  S3=size(a,3);                 % Get the number of color planes
  if (nargin < 3 )
      error('MATLAB:wavelet:NoInputs',['No input arguments specified.'])
  end
  if (numDimsA<1 || numDimsA>3) % Only 1-3 dimensional arrays are allowed
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.']) 
  end
  if not (sdomain>0 && sdomain<10 && srange >0 && srange<=256)
      error('MATLAB:wavelet:NoInputs',['Invalid filter parameters.'])
  end
  if not ( (S1==1 && S2>1) || (S1>1 && S2==1) || (S1>1 && S2>1) )
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.'])
  end
  if (S1>1 && S2>1 && sdomain>10)      % Set a 2-D limit as a precaution
      error('MATLAB:wavelet:NoInputs',['Possibly excessive domain sigma.'])
  end  

  % Determine necessary filtering variables
  R=ceil(3*sdomain);    % Compute radius (domain wt. <1% outside this radius)
  R1=R;                 % Set row and column radius limit
  R2=R;                 % (these are practical limitations to speed up the process)
  if (S1==1)            % If input is a single row
    R1=0;               % Keep filter isolated to that row
  end
  if (S2==1)            % If input is a single column
    R2=0;               % Keep filter isolated to that column 
  end
  WT=zeros(S1,S2);      % Reset weight total array
 
  % Run the bilateral filter 
  b=zeros(S1,S2,S3);    % Reset output matrix
  for C1=-R1:R1         % Loop through potential circular shift values
    for C2=-R2:R2         
      DS=C1^2+C2^2;     % Compute pixel distance squared
      if (DS < R^2)     % If within radius limit...
        t=circshift(a,[C1 C2 0]);               % Compute shifted matrix
        W = exp(-DS/(2*sdomain^2)) * exp(-sum((t-a).^2,3)/(2*srange^2));
        WT = WT + W;                            % Compute weights and weight total
        for color=1:S3                          % Loop through color planes
          b(:,:,color) = b(:,:,color) + t(:,:,color).*W;   % Add to each color total
        end
      end
    end
  end
  for color=1:S3                      % Loop through color planes
    b(:,:,color)=b(:,:,color)./WT;    % Normalize the values
  end