• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# [Matlab] Bilateral Filter

LightningIsMyName Sep 8th, 2011 2,588 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. % Bilateral Filter for Matlab/Octave
2.
3. % Author: Barak Itkin
4. % Author Website: http://barak-itkin.blogspot.com/
5. % Created: 2011-09-08
6.
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.
21. %
22. % The equation of the bilateral filter is
23. %
24. %         (       dx ^ 2       )       (         dI ^2        )
25. % F = exp (- ----------------- ) * exp (- ------------------- )
26. %         (  sigma_spatial ^ 2 )       (  sigma_Intensity ^ 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 intensity of the 'center pixel' and
33. %      the pixel to sample
34. % sigma_spatial and sigma_Intesity are constants. Higher values mean
35. % that we 'tolerate more' higher value of the distances dx and dI.
36.
37. % Usage:
38. %   im         - A GRAY image
39. %   filterSize - The size of the filter (the lookup area would be
40. %                filterSize X filterSize). This is a real positive
41. %                integer. This MUST BE AN ODD (NON EVEN) REAL NUMBER!
42. %   sSpatial   - The sigma_spatial of the bilateral filter equation.
43. %                Typically, this is half of the filterSize.
44. %   SIntensity - The sigma_Intensity of the bilateral filter equation.
45. %
46. % Return:
47. %   ret        - The resulting image.
48.
49. function [ ret ] = bilateral (im, filterSize, sSpatial, sIntensity)
50.
51. if (not (nargin == 4))
52.         error ('Expected exactly 4 arguments!')
53. end
54.
55. if (not (ismatrix(im)) || not (ndims (im) == 2))
56.         error ('Argument 1 (im) must be a GRAY image')
57. end
58.
59. if (not (isscalar (filterSize)) || filterSize <= 0 || not (isreal (filterSize)) || not (rem (filterSize, 2) == 1))
60.         error ('Argument 2 (filterSize) must be a positive real odd integer!')
61. end
62.
63. if (not (isscalar (sSpatial)) || not (isreal (sSpatial)))
64.         error ('Argument 3 (sSpatial) must be a real number!')
65. end
66.
67. if (not (isscalar (sIntensity)) || not (isreal (sIntensity)))
68.         error ('Argument 4 (sIntensity must be a real number!')
69. end
70.
71. % Compute the Gaussian filter part of the Bilateral filter
72. gauss = fspecial ('gaussian', filterSize, sSpatial);
73.
74. % The radius of the filter
75. radius = floor (filterSize / 2);
76.
77. % The original image dimensions
78. [height,width] = size (im);
79.
80. % The resulting image
81. ret = zeros (size (im));
82.
83. % We will now compute the result by iterating over the result image,
84. % computing one pixel at a time. For the pixel at (i,j) we want to look
85. % at:
86. %
88. %
89. % But, we we must also remember to restrict ourselves to remain inside
90. % the image boundries
91. for i=1:height
92.
93.         % The top part of the lookup area, clamped to the image
94.         ymin  = max (i - radius, 1);
95.         % How many rows were outside the image, on the top?
96.         dymin = ymin - (i - radius);
97.
98.         % The bottom part of the lookup area, clamped to the image
99.         ymax  = min (i + radius, height);
100.         % How many rows were outside the image, on the bottom?
101.         dymax = (i + radius) - ymax;
102.
103.         for j=1:width
104.
105.                 % The left part of the lookup area, clamped to the image
106.                 xmin = max (j - radius, 1);
107.                 % How many columns were outside the image, on the left?
108.                 dxmin = xmin - (j - radius);
109.
110.                 % The right part of the lookup area, clamped to the image
111.                 xmax = min (j + radius, width);
112.                 % How many columns were outside the image, on the right?
113.                 dxmax = (j + radius) - xmax;
114.
115.                 % The actual area of the image we will look at
116.                 area = im (ymin:ymax, xmin:xmax);
117.
118.                 % The center pixel
119.                 center = im (i,j);
120.
121.                 % The left expression in the bilateral filter equation
122.                 % We take only the relevant parts of the matrix of the
123.                 % Gaussian weights - we use dxmin, dxmax, dymin, dymax to
124.                 % ignore the parts that are outside the image
125.                 expS = gauss((1+dymin):(filterSize-dymax),(1+dxmin):(filterSize-dxmax));
126.
127.                 % The intensity difference from the center pixel
128.                 dI = area - center;
129.                 % The right expression in the bilateral filter equation
130.                 expI = exp((-dI .* dI) / (sIntensity * sIntensity));
131.
132.                 % The bilater filter (weights matrix)
133.                 F = expI .* expS;
134.                 % Normalized bilateral filter
135.                 Fnormalized = F / sum(F(:));
136.
137.                 % Multiply the area by the filter
138.                 temp = area .* Fnormalized;
139.
140.                 % The resulting pixel is the sum of all the pixels in
141.                 % the area, according to the weights of the filter
142.                 ret(i,j) = sum (temp(:));
143.         end
144. end
RAW Paste Data
Top