matlab

关于Sylvain Paris ,Fr´edo Durand的快速双边滤波算法

2019-07-17 13:37发布

站内问答 / matlab
1663 2 987
本帖最后由 lxzxj_chan 于 2015-10-23 12:07 编辑

这段代总体思路是什么，怎样实现双边滤波？




  4. %
  5. % output = bilateralFilter( data, edge, sigmaSpatial, sigmaRange, ...
  6. %                          samplingSpatial, samplingRange )
  7. %
  8. % Bilateral and Cross-Bilateral Filter
  9. %
  10. % Bilaterally filters the image 'data' using the edges in the image 'edge'.
  11. % If 'data' == 'edge', then it the normal bilateral filter.
  12. % Else, then it is the "cross" or "joint" bilateral filter.
  13. %
  14. % Note that for the cross bilateral filter, data does not need to be
  15. % defined everywhere.  Undefined values can be set to 'NaN'.  However, edge
  16. % *does* need to be defined everywhere.
  17. %
  18. % data and edge should be of the same size and greyscale.
  19. % (i.e. they should be ( height x width x 1 matrices ))
  20. %
  21. % data is the only required argument
  22. %
  23. % By default:
  24. % edge = data
  25. % sigmaSpatial = samplingSpatial = min( width, height ) / 16;
  26. % sigmaRange = samplingRange = ( max( edge( : ) ) - min( edge( : ) ) ) / 10
  27. %
  28. %

  30. % Copyright (c) <2007> <Jiawen Chen, Sylvain Paris, and Fredo Durand>
  31. %
  32. % Permission is hereby granted, free of charge, to any person obtaining a copy
  33. % of this software and associated documentation files (the "Software"), to deal
  34. % in the Software without restriction, including without limitation the rights
  35. % to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  36. % copies of the Software, and to permit persons to whom the Software is
  37. % furnished to do so, subject to the following conditions:
  38. %
  39. % The above copyright notice and this permission notice shall be included in
  40. % all copies or substantial portions of the Software.
  41. %
  42. % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  43. % IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  44. % FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  45. % AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  46. % LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  47. % OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  48. % THE SOFTWARE.
  49. %

  51. function output = bilateralFilter( data, edge, sigmaSpatial, sigmaRange, ...
  52.     samplingSpatial, samplingRange )

  54. if ~exist( 'edge', 'var' ),
  55.     edge = data;
  56. end

  58. inputHeight = size( data, 1 );
  59. inputWidth = size( data, 2 );

  61. if ~exist( 'sigmaSpatial', 'var' ),
  62.     sigmaSpatial = min( inputWidth, inputHeight ) / 16;
  63. end

  65. edgeMin = min( edge( : ) );
  66. edgeMax = max( edge( : ) );
  67. edgeDelta = edgeMax - edgeMin;

  69. if ~exist( 'sigmaRange', 'var' ),
  70.     sigmaRange = 0.1 * edgeDelta;
  71. end

  73. if ~exist( 'samplingSpatial', 'var' ),
  74.     samplingSpatial = sigmaSpatial;
  75. end

  77. if ~exist( 'samplingRange', 'var' ),
  78.     samplingRange = sigmaRange;
  79. end

  81. if size( data ) ~= size( edge ),
  82.     error( 'data and edge must be of the same size' );
  83. end

  85. % parameters
  86. derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
  87. derivedSigmaRange = sigmaRange / samplingRange;

  89. paddingXY = floor( 2 * derivedSigmaSpatial ) + 1;
  90. paddingZ = floor( 2 * derivedSigmaRange ) + 1;

  92. % allocate 3D grid
  93. downsampledWidth = floor( ( inputWidth - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
  94. downsampledHeight = floor( ( inputHeight - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
  95. downsampledDepth = floor( edgeDelta / samplingRange ) + 1 + 2 * paddingZ;

  97. gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
  98. gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );

  100. % compute downsampled indices
  101. [ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 );

  103. % ii =
  104. % 0 0 0 0 0
  105. % 1 1 1 1 1
  106. % 2 2 2 2 2

  108. % jj =
  109. % 0 1 2 3 4
  110. % 0 1 2 3 4
  111. % 0 1 2 3 4

  113. % so when iterating over ii( k ), jj( k )
  114. % get: ( 0, 0 ), ( 1, 0 ), ( 2, 0 ), ... (down columns first)

  116. di = round( ii / samplingSpatial ) + paddingXY + 1;
  117. dj = round( jj / samplingSpatial ) + paddingXY + 1;
  118. dz = round( ( edge - edgeMin ) / samplingRange ) + paddingZ + 1;

  120. % perform scatter (there's probably a faster way than this)
  121. % normally would do downsampledWeights( di, dj, dk ) = 1, but we have to
  122. % perform a summation to do box downsampling
  123. for k = 1 : numel( dz ),
  124.       
  125.     dataZ = data( k ); % traverses the image column wise, same as di( k )
  126.     if ~isnan( dataZ  ),
  127.         
  128.         dik = di( k );
  129.         djk = dj( k );
  130.         dzk = dz( k );

  132.         gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
  133.         gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;
  134.         
  135.     end
  136. end

  138. % make gaussian kernel
  139. kernelWidth = 2 * derivedSigmaSpatial + 1;
  140. kernelHeight = kernelWidth;
  141. kernelDepth = 2 * derivedSigmaRange + 1;

  143. halfKernelWidth = floor( kernelWidth / 2 );
  144. halfKernelHeight = floor( kernelHeight / 2 );
  145. halfKernelDepth = floor( kernelDepth / 2 );

  147. [gridX, gridY, gridZ] = meshgrid( 0 : kernelWidth - 1, 0 : kernelHeight - 1, 0 : kernelDepth - 1 );
  148. gridX = gridX - halfKernelWidth;
  149. gridY = gridY - halfKernelHeight;
  150. gridZ = gridZ - halfKernelDepth;
  151. gridRSquared = ( gridX .* gridX + gridY .* gridY ) / ( derivedSigmaSpatial * derivedSigmaSpatial ) + ( gridZ .* gridZ ) / ( derivedSigmaRange * derivedSigmaRange );
  152. kernel = exp( -0.5 * gridRSquared );

  154. % convolve
  155. blurredGridData = convn( gridData, kernel, 'same' );
  156. blurredGridWeights = convn( gridWeights, kernel, 'same' );

  158. % divide
  159. blurredGridWeights( blurredGridWeights == 0 ) = -2; % avoid divide by 0, won't read there anyway
  160. normalizedBlurredGrid = blurredGridData ./ blurredGridWeights;
  161. normalizedBlurredGrid( blurredGridWeights < -1 ) = 0; % put 0s where it's undefined
  162. blurredGridWeights( blurredGridWeights < -1 ) = 0; % put zeros back

  164. % upsample
  165. [ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 ); % meshgrid does x, then y, so output arguments need to be reversed
  166. % no rounding
  167. di = ( ii / samplingSpatial ) + paddingXY + 1;
  168. dj = ( jj / samplingSpatial ) + paddingXY + 1;
  169. dz = ( edge - edgeMin ) / samplingRange + paddingZ + 1;

  171. % interpn takes rows, then cols, etc
  172. % i.e. size(v,1), then size(v,2), ...
  173. output = interpn( normalizedBlurredGrid, di, dj, dz );
复制代码

友情提示: 此问题已得到解决，问题已经关闭,关闭后问题禁止继续编辑，回答。
2条回答
lxzxj_chan
2019-07-17 17:46
haohaoren 发表于 2015-10-23 10:15
我觉得啊，作为一个提问题的人，应该不要那么懒，你提出一个问题，直接问人家怎么理解，你都没交代这个是什么问题，你都没告诉我你怎么理解的，我也没有义乌给你说吧，所以，我觉得现代人最大的问题啊，就是不懂得怎样虚心求教。

从我个人角度看，你这种学习态度，也很难学习很好吧！

好的，谢谢提醒

一周热门 更多>