SHARE
TWEET

[Matlab] Bilateral filter for RGB images

LightningIsMyName Sep 23rd, 2011 501 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % Bilateral Filter for RGB images for Matlab/Octave
  2.  
  3. % Author: Barak Itkin
  4. % Author Website: http://barak-itkin.blogspot.com/
  5. % Created: 2011-09-08
  6.  
  7. % This work is licensed under a Creative Commons Attribution-ShareAlike
  8. % 3.0 Unported License.
  9. % See http://creativecommons.org/licenses/by-sa/3.0/
  10.  
  11. % The bilateral filter is a 'Smart Blur' filter that avoids smoothing
  12. % edges. It does so by considering the intensity differences between the
  13. % pixels - if the difference between the 'center pixel' and another
  14. % pixel is too high, then that pixel get's less weight in the weight
  15. % matrix of the blur. In addition to considering the intensity
  16. % difference, we also consider the actual distance (as in a regular
  17. % Gaussian blur).
  18. % The bilateral filter is most commonly used as a noise reduction filter
  19. % due to it's ability to smooth large areas without destroying edges.
  20. % For more information, see http://en.wikipedia.org/wiki/Bilateral_filter
  21. %
  22. % The equation of the bilateral filter is
  23. %
  24. %         (       dx ^ 2       )       (     dI ^ 2      )
  25. % F = exp (- ----------------- ) * exp (---------------- )
  26. %         (  sigma_spatial ^ 2 )       ( sigma_color ^ 2 )
  27. %     ~~~~~~~~~~~~~~~~~~~~~~~~~~
  28. %     This is a guassian filter!
  29. %
  30. % dx - The 'geometric' distance between the 'center pixel' and the pixel
  31. %      to sample
  32. % dI - The difference between the color of the 'center pixel' and
  33. %      the pixel to sample. To calculate it, we calculate the root of the
  34. %      avarege of the squared differences in each of the RGB channels
  35. %      
  36. % sigma_spatial and sigma_color are constants. Higher values mean
  37. % that we 'tolerate more' higher value of the distances dx and dI.
  38.  
  39. % Usage:
  40. %   im         - A RGB image
  41. %   filterSize - The size of the filter (the lookup area would be
  42. %                filterSize X filterSize). This is a real positive
  43. %                integer. This MUST BE AN ODD (NON EVEN) REAL NUMBER!
  44. %   sSpatial   - The sigma_spatial of the bilateral filter equation.
  45. %                Typically, this is half of the filterSize.
  46. %   sColor     - The sigma_color of the bilateral filter equation.
  47. %
  48. % Return:
  49. %   ret        - The resulting image.
  50.  
  51. function [ ret ] = bilateralRGB (im, filterSize, sSpatial, sColor)
  52.  
  53. R = 1;
  54. G = 2;
  55. B = 3;
  56. RGB = [R G B];
  57.  
  58. if (not (nargin == 4))
  59.         error ('Expected exactly 4 arguments!')
  60. end
  61.  
  62. if (not (ismatrix(im)) || not (ndims (im) == 3))
  63.         error ('Argument 1 (im) must be a RGB image')
  64. end
  65.  
  66. if (not (isscalar (filterSize)) || filterSize <= 0 || not (isreal (filterSize)) || not (rem (filterSize, 2) == 1))
  67.         error ('Argument 2 (filterSize) must be a positive real odd integer!')
  68. end
  69.  
  70. if (not (isscalar (sSpatial)) || not (isreal (sSpatial)))
  71.         error ('Argument 3 (sSpatial) must be a real number!')
  72. end
  73.  
  74. if (not (isscalar (sColor)) || not (isreal (sColor)))
  75.         error ('Argument 4 (sColor must be a real number!')
  76. end
  77.  
  78. % Compute the Gaussian filter part of the Bilateral filter
  79. gauss = fspecial ('gaussian', filterSize, sSpatial);
  80.  
  81. % The radius of the filter
  82. radius = floor (filterSize / 2);
  83.  
  84. % The original image dimensions
  85. [height,width,depth] = size (im);
  86.  
  87. % The resulting image
  88. ret = zeros (size (im));
  89.  
  90. % We will now compute the result by iterating over the result image,
  91. % computing one pixel at a time. For the pixel at (i,j) we want to look
  92. % at:
  93. %
  94. %      im((i - radius):(i + radius), (j - radius):(j + radius))
  95. %
  96. % But, we we must also remember to restrict ourselves to remain inside
  97. % the image boundries
  98. for i=1:height
  99.  
  100.         % The top part of the lookup area, clamped to the image
  101.         ymin  = max (i - radius, 1);
  102.         % How many rows were outside the image, on the top?
  103.         dymin = ymin - (i - radius);
  104.        
  105.         % The bottom part of the lookup area, clamped to the image
  106.         ymax  = min (i + radius, height);
  107.         % How many rows were outside the image, on the bottom?
  108.         dymax = (i + radius) - ymax;
  109.        
  110.         for j=1:width
  111.        
  112.                 % The left part of the lookup area, clamped to the image
  113.                 xmin = max (j - radius, 1);
  114.                 % How many columns were outside the image, on the left?
  115.                 dxmin = xmin - (j - radius);
  116.                
  117.                 % The right part of the lookup area, clamped to the image
  118.                 xmax = min (j + radius, width);
  119.                 % How many columns were outside the image, on the right?
  120.                 dxmax = (j + radius) - xmax;
  121.                
  122.                 % The actual area of the image we will look at
  123.                 area = im (ymin:ymax, xmin:xmax, RGB);
  124.                
  125.                 % The center pixel
  126.                 center = im (i, j, RGB);
  127.  
  128.                 % The left expression in the bilateral filter equation
  129.                 % We take only the relevant parts of the matrix of the
  130.                 % Gaussian weights - we use dxmin, dxmax, dymin, dymax to
  131.                 % ignore the parts that are outside the image
  132.                 expS = gauss((1+dymin):(filterSize-dymax),(1+dxmin):(filterSize-dxmax));
  133.  
  134.                 dR = area(:,:,R) - center (:,:,R);
  135.                 dG = area(:,:,G) - center (:,:,G);
  136.                 dB = area(:,:,B) - center (:,:,B);
  137.  
  138.                 dIsquare = ((dR .* dR) + (dG .* dG) + (dB .* dB)) / 3;
  139.  
  140.                 % The right expression in the bilateral filter equation
  141.                 expI = exp (- dIsquare / (sColor * sColor));
  142.                
  143.                 % The bilater filter (weights matrix)
  144.                 F = expI .* expS;
  145.  
  146.                 % Normalized bilateral filter
  147.                 Fnormalized = F / sum(F(:));
  148.                
  149.                 % Multiply the area by the filter
  150.                 tempR = area(:,:,R) .* Fnormalized;
  151.                 tempG = area(:,:,G) .* Fnormalized;
  152.                 tempB = area(:,:,B) .* Fnormalized;
  153.                
  154.                 % The resulting pixel is the sum of all the pixels in
  155.                 % the area, according to the weights of the filter
  156.                 ret(i,j,R) = sum (tempR(:));
  157.                 ret(i,j,G) = sum (tempG(:));
  158.                 ret(i,j,B) = sum (tempB(:));
  159.         end
  160. end
RAW Paste Data
Pastebin PRO Summer Special!
Get 40% OFF on Pastebin PRO accounts!
Top