image thumbnail

Fast Bilateral Filter

by

 

11 May 2012 (Updated )

Fast O(1) implementation of the bilateral and the joint (or cross) bilateral filter.

shiftable_jointBF(inImg, rangeImg, sigmaS, sigmaR, w, tol)
function [ outImg , param ]  =  shiftable_jointBF(inImg, rangeImg, sigmaS, sigmaR, w, tol)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Filtering operation
%
%   inImg      : MxNxB image to filter (single/double)
%   rangeImg   : MxN image used to for pixel similarity (single/double)
%   outImg     : filtered image of size MxNxB
%   sigmaS     : width of spatial Gaussian
%   sigmaR     : width of range Gaussian
%   [-w, w]^2  : domain of spatial Gaussian
%   tol        : truncation error
%
%   Author: Kunal N. Chaudhury
%   Date:   March 1, 2012
%
%   Modified by: Derek Hoiem
%   Date: May 18, 2012
%   Modifications: Allow inImg to be multiband; allow image for computing
%   pixel similarity to be different than filtered image.
%
%   Reference:  
%
%   [1] K.N. Chaudhury, D. Sage, and M. Unser, "Fast O(1) bilateral  filtering using 
%   trigonometric range kernels," IEEE Transactions on Image Processing, vol. 20, 
%   no. 11, 2011.
%
% [2] K.N. Chaudhury, "Acceleration of the shiftable O(1) algorithm for
% bilateral filtering and non-local means,"  IEEE Transactions on Image Proc., 
% vol. 22, no. 4, 2013. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% create spatial filter
filt     = fspecial('gaussian', [w w], sigmaS);

% set range interval and the order of raised cosine
T  =  maxFilter(rangeImg, w);
N  =  ceil( 0.405 * (T / sigmaR)^2 );

gamma    =  1 / (sqrt(N) * sigmaR);
twoN     =  2^N;

% compute truncation
if tol == 0
    M = 0;
else
    if sigmaR > 40
        M = 0;
    elseif sigmaR > 10
        sumCoeffs = 0;
        for k = 0 : round(N/2)
            sumCoeffs = sumCoeffs + nchoosek(N,k)/twoN;
            if sumCoeffs > tol/2
                M = k;
                break;
            end
        end
    else
        M = ceil( 0.5 * ( N - sqrt(4*N*log(2/tol)) ) );
    end
end


% main filter
[m, n, b]  =  size(inImg);
outImg1    =  zeros(m, n, b);
outImg2    =  zeros(m, n, b);
outImg     =  zeros(m, n, b);

warning off; %#ok<WNOFF>

for k = M : N - M
    
    coeff = nchoosek(N,k) / twoN;
    
    temp1  = cos( (2*k - N) * gamma * rangeImg);
    temp2  = sin( (2*k - N) * gamma * rangeImg);
    
    if size(inImg, 3) > 1
      temp1 = repmat(temp1, [1 1 size(inImg, 3)]);
      temp2 = repmat(temp2, [1 1 size(inImg, 3)]);
    end
    
    phi1 =  imfilter(inImg .* temp1, filt);
    phi2 =  imfilter(inImg .* temp2, filt);
    phi3 =  imfilter(temp1, filt);
    phi4 =  imfilter(temp2, filt);
    
    outImg1 = outImg1 + coeff * ( temp1 .* phi1 +  temp2 .* phi2 );
    outImg2 = outImg2 + coeff * ( temp1 .* phi3 +  temp2 .* phi4 );
    
end

idx1 = find( outImg2 < 0.0001);
idx2 = find( outImg2 > 0.0001);

outImg( idx1 ) = inImg( idx1 );
outImg( idx2 ) = outImg1( idx2 ) ./ outImg2 (idx2 );

% save parameters
param.T  = T;
param.N  = N;
param.M  = M;

Contact us