dworkLib/bilateralFilter.m at master · ndwork/dworkLib · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

function out = bilateralFilter( img, varargin )
  % out = bilateralFilter( img, [ 'S', S, 'sigmaD', sigmaD, ...
  %   'sigmaR', sigmaR, 'verbose', verbose ] );
  %
  % Inputs
  %   img: the input image to be denoised
  %     either a MxN 2D grayscale image or a MxNx3 color image
  %
  % Optional Inputs:
  %   S: the length of each side of the kernel (must be odd)
  %   sigmaD: euclidean distance spread parameter
  %   sigmaR: photometric distance spread parameter
  %   verbose: displays status
  %
  % Outputs
  %   out: the denoised image

  if nargin < 1
    disp( 'Usage:  out = bilateralFilter( img, [ ''S'', S, ''sigmaD'', sigmaD, ... ' );
    disp( '  ''sigmaR'', sigmaR, ''verbose'', verbose ] );' );
    return
  end

  isPositiveAndOdd = @(x) (x>0) && (mod(x,2) == 1);
  isNumericAndPositive = @(x) isnumeric(x) && (x>0);

  defaultSigmaD = 3;
  defaultSigmaR = 0.3;
  p = inputParser;
  p.addRequired( 'img', @isnumeric );
  p.addParameter( 'S', [], isPositiveAndOdd );
  p.addParameter( 'sigmaD', defaultSigmaD, isNumericAndPositive );
  p.addParameter( 'sigmaR', defaultSigmaR, isNumericAndPositive );
  p.addParameter( 'verbose', false, @(x) islogical(x) || isnumeric(x) );
  p.parse( img, varargin{:} );
  s = p.Results.S;
  sigmaD = p.Results.sigmaD;
  sigmaR = p.Results.sigmaR;
  verbose = p.Results.verbose;

  if numel( s ) == 0, s = 3 * sigmaD; end

  if ismatrix(img)
    out = bilateralFilter_2D( img, s, sigmaD, sigmaR, verbose );
  elseif ndims(img)==3
    out = bilateralFilter_3D( img, s, sigmaD, sigmaR, verbose );
  end
end


function out = bilateralFilter_2D( img, s, sigmaD, sigmaR, verbose )

  dKernel = fspecial( 'gaussian', s, sigmaD );
  halfS = floor(s/2);
  varR = sigmaR * sigmaR;

  sImg = size( img );

  firstJ = ceil(s/2);   lastJ = sImg(1)-floor(s/2);
  firstI = ceil(s/2);   lastI = sImg(2)-floor(s/2);
  jCells = cell( lastJ - firstJ + 1, 1 );

  p = parforProgress( lastJ-firstJ+1 );
  parfor j=1:lastJ-firstJ+1
    if verbose == true, p.progress( j, 50 ); end

    tmp = zeros( 1, lastI-firstI+1 );
    imgLine = img( j+firstJ-1-halfS:j+firstJ-1+halfS, : );   %#ok<PFBNS>

    for i=1:lastI-firstI+1
      %subImg = img( j-halfS:j+halfS, i-halfS:i+halfS );
      subImg = imgLine( :, i+firstI-1-halfS : i+firstI-1+halfS );

      pKernel = exp( -( img(j+firstJ-1,i+firstI-1) - subImg ).^2 / ( 2 * varR ) );

      weights = pKernel .* dKernel;
      weights = weights / sum( weights(:) );

      tmp(i) = sum( subImg(:) .* weights(:) );
    end
    jCells{j} = tmp;
  end
  p.clean;
  out = img;
  out(firstJ:lastJ,firstI:lastI) = cell2mat( jCells );
end


function out = bilateralFilter_3D( img, s, sigmaD, sigmaR, verbose )

  dKernel = fspecial3d( 'gaussian', s, sigmaD );
  halfS = floor(s/2);
  varR = sigmaR * sigmaR;

  sImg = size( img );
  firstJ = ceil(s/2);   lastJ = sImg(1)-floor(s/2);
  firstI = ceil(s/2);   lastI = sImg(2)-floor(s/2);
  firstK = ceil(s/2);   lastK = sImg(3)-floor(s/2);
  sliceCells = cell(1,1,lastK-firstK+1);

  p = parforProgress( lastK-firstK+1 );
  parfor k=1:lastK-firstK+1
    if verbose == true, p.progress( k ); end   %#ok<PFBNS>

    tmp = zeros( sImg(1:2) );   %#ok<PFBNS>
    kImg = img(:,:,k+firstK-1-halfS:k+firstK-1+halfS);   %#ok<PFBNS>

    for j=1:lastJ-firstJ+1
      for i=1:lastI-firstI+1

        subImg = kImg( j+firstJ-1-halfS:j+firstJ-1+halfS, ...
                       i+firstI-1-halfS:i+firstI-1+halfS, : );

        pKernel = exp( -( img(j,i,k) - subImg ).^2 / ( 2 * varR ) );

        weights = pKernel .* dKernel;
        weights = weights / sum( weights(:) );

        tmp(j,i) = sum( subImg(:) .* weights(:) );
      end
    end
    sliceCells{1,1,k} = tmp;
  end
  p.clean;
  out = img;
  out(firstJ:lastJ,firstI:lastI,firstK:lastK) = cell2mat( sliceCells );
end