Bilateral Filter - MATLAB

bilateral_filter.m

clc, close all, clear all;

img = imread('puppy.jpg');
A = rgb2gray(img);
[m,n] = size(A);

%Filter Parameters setup

windowSize = input('Enter window size (3/5/7):\n ');
stDevSpace = input('Enter space standard deviation (1-20):\n ');
stDevColor = input('Enter color standard deviation (1-20):\n ');

switch windowSize
    case 3
        disp('Window size 3x3 is set')
        wNum = 1;
        centerPixelX = 1 + wNum;
        centerPixelY = 1 + wNum;
        xStartBoundary = 1 + wNum;
        xEndBoundary = m - wNum;
        yStartBoundary = 1 + wNum;
        yEndBoundary = n - wNum;
    case 5
        disp('Window size 5x5 is set')
        wNum = 2;
        centerPixelX = 1 + wNum;
        centerPixelY = 1 + wNum;
        xStartBoundary = 1 + wNum;
        xEndBoundary = m - wNum;
        yStartBoundary = 1 + wNum;
        yEndBoundary = n - wNum;
    case 7
        disp('Window size 7x7 is set')
        wNum = 3;
        centerPixelX = 1 + wNum;
        centerPixelY = 1 + wNum;
        xStartBoundary = 1 + wNum;
        xEndBoundary = m - wNum;
        yStartBoundary = 1 + wNum;
        yEndBoundary = n - wNum;
    otherwise
        disp('Please choose something')
end


for i = xStartBoundary:xEndBoundary
    for j = yStartBoundary:yEndBoundary
        windowM = A([i-wNum:i+wNum],[j-wNum:j+wNum]);

        sumTogether = 0.0;
        sumOutput = 0;

        %% Before Normalization
        for runX = 1:windowSize
            for runY = 1:windowSize
                meanDifX = runX-centerPixelX;
                meanDifY = runY-centerPixelY;
                
                difSpace = double(sqrt(meanDifX^2 + meanDifY^2));
                difColor = double(abs(windowM(runX,runY) - A(i,j)));
                
                beforeNormSpace(runX,runY) = (0.399/stDevSpace) * exp(-difSpace/(2*stDevSpace^2)); %0.399 = 1/sqrt(2pi)
                beforeNormColor(runX,runY) = (0.399/stDevColor) * exp(-difColor/(2*stDevColor^2)); %0.399 = 1/sqrt(2pi)  
                
                sumTogether = sumTogether + beforeNormSpace(runX,runY)*beforeNormColor(runX,runY);
            end
        end
        %% After Normalization
        for runXX = 1:windowSize
            for runYY = 1:windowSize
                afterNormMatrix(runXX,runYY) = (beforeNormSpace(runXX,runYY)*beforeNormColor(runXX,runYY))/sumTogether;
            end
        end 
        %% Convolution with input
        for runXXX = 1:windowSize
            for runYYY = 1:windowSize
                sumOutput = sumOutput + afterNormMatrix(runXXX,runYYY) * windowM(runXXX,runYYY);
            end
        end
        %%
        outputImage(i,j) = sumOutput;
    end
end

%Ignoring last rows and columns
for i = 1:m
    for j = 1:n
        if (i == 1 || j == 1 || i == m || j == n)
            outputImage(i,j) = A(i,j);
        end
    end
end

%Result
figure(1)
subplot(1,2,1), imshow(A);
title('Input:');
subplot(1,2,2), imshow(outputImage);
str = sprintf('wndw = %d, stdspc = %d, stdclr = %d', windowSize, stDevSpace, stDevColor);
title(str)