Bilateral Filters(双边滤波算法)原理及实现

双边滤波算法原理:

双边滤波是一种非线性滤波器,它可以达到保持边缘、降噪平滑的效果。和其他滤波原理一样,双边滤波也是采用加权平均的方法,用周边像素亮度值的加权平均代表某个像素的强度,所用的加权平均基于高斯分布[1]。最重要的是,双边滤波的权重不仅考虑了像素的欧氏距离(如普通的高斯低通滤波,只考虑了位置对中心像素的影响),还考虑了像素范围域中的辐射差异(例如卷积核中像素与中心像素之间相似程度、颜色强度,深度距离等),在计算中心像素的时候同时考虑这两个权重。 公式1a,1b给出了双边滤过的操作,Iq为输入图像,Ipbf为滤波后图像:


mark下双边滤波里的两个权重域的概念:空间域(spatial domain S)和像素范围域(range domain R),这个是它跟高斯滤波等方法的最大不同点。下面是我找到的对比说明,更好地理解双边滤波,首先是高斯滤波的情况:


然后对比再看一下双边滤波的过程:


双边滤波的核函数是空间域核与像素范围域核的综合结果:在图像的平坦区域,像素值变化很小,对应的像素范围域权重接近于1,此时空间域权重起主要作用,相当于进行高斯模糊;在图像的边缘区域,像素值变化很大,像素范围域权重变大,从而保持了边缘的信息。


为了更加形象的说明两个权重的影响,作者还给出了二维图像的直观说明:


双边滤波算法实现:

在原理部分,从双边滤波的公式就可以得到该算法的实现途径。由于直接的编码实现上述过程,其时间复杂度为O(σs2) ,比较耗时,所以后来出现了一些改进算法,比较经典的有:论文《Fast O(1) bilateral filtering using trigonometric range kernels》,提出了用Raised cosines函数来逼近高斯值域函数,并利用一些特性把值域函数分解为一些列函数的叠加,从而实现函数的加速[5,8]。
这里只对原始方法进行实现,从而有助于更加清楚的了解算法的原理。

matlab实现方法,这里也附一下核心代码:

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

if ~exist( 'edge', 'var' ),
    edge = data;
end

inputHeight = size( data, 1 );
inputWidth = size( data, 2 );

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

edgeMin = min( edge( : ) );
edgeMax = max( edge( : ) );
edgeDelta = edgeMax - edgeMin;

if ~exist( 'sigmaRange', 'var' ),
    sigmaRange = 0.1 * edgeDelta;
end

if ~exist( 'samplingSpatial', 'var' ),
    samplingSpatial = sigmaSpatial;
end

if ~exist( 'samplingRange', 'var' ),
    samplingRange = sigmaRange;
end

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

% parameters
derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
derivedSigmaRange = sigmaRange / samplingRange;

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

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

gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );

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

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

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

        gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
        gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;
        
    end
end

% make gaussian kernel
kernelWidth = 2 * derivedSigmaSpatial + 1;
kernelHeight = kernelWidth;
kernelDepth = 2 * derivedSigmaRange + 1;

halfKernelWidth = floor( kernelWidth / 2 );
halfKernelHeight = floor( kernelHeight / 2 );
halfKernelDepth = floor( kernelDepth / 2 );

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

% convolve
blurredGridData = convn( gridData, kernel, 'same' );
blurredGridWeights = convn( gridWeights, kernel, 'same' );

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

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

% interpn takes rows, then cols, etc
% i.e. size(v,1), then size(v,2), ...
output = interpn( normalizedBlurredGrid, di, dj, dz );

双边滤波算法实例:

%% 测试函数,其中参数设置请参见函数注释
clc,clear all,close all;
ori=imread('D:\proMatlab\vessel_edge_extration\image\3.jpg'); 
ori=double(rgb2gray(ori))/255.0;
[width, height]=size(ori);
sigmaSpatial  = min( width, height ) / 30;
samplingSpatial=sigmaSpatial;
sigmaRange = ( max( ori( : ) ) - min( ori( : ) ) ) / 30;
samplingRange= sigmaRange;
output = bilateralFilter( ori, ori, sigmaSpatial, sigmaRange, ...
    samplingSpatial, samplingRange );
figure,
subplot(1,2,1),imshow(ori,[]);title('input image');
subplot(1,2,2),imshow(output,[]);title('output image');

参考:

  1. https://en.wikipedia.org/wiki/Bilateral_filtering
  2. http://people.csail.mit.edu/sparis/bf_course/
  3. http://people.csail.mit.edu/sparis/bf/
  4. http://blog.csdn.net/fightingforcv/article/details/52723376
  5. http://blog.csdn.net/mumusan2016/article/details/54578038
  6. http://blog.csdn.net/majinlei121/article/details/50463514
  7. http://kaiminghe.com/eccv10/index.html  (Guided Image Filtering
  8. http://blog.csdn.net/dangchangying/article/details/14451963
  • 86
    点赞
  • 452
    收藏
    觉得还不错? 一键收藏
  • 17
    评论
自适应双边滤波算法是一种常用的图像滤波算法,其可以在保持图像边缘信息的同时去除图像噪声。下面是C++实现自适应双边滤波算法的代码,具体注释见代码中: ```c++ #include <opencv2/opencv.hpp> #include <iostream> using namespace cv; using namespace std; // 定义自适应双边滤波函数 Mat adaptive_bilateral_filter(Mat input, int size, double sigma_d, double sigma_r) { Mat output = Mat::zeros(input.size(), input.type()); // 初始化输出图像 int radius = size / 2; // 计算滤波半径 int rows = input.rows; int cols = input.cols; int channels = input.channels(); double c1 = 1.0 / (2.0 * sigma_d * sigma_d); // 计算常数C1 double c2 = 1.0 / (2.0 * sigma_r * sigma_r); // 计算常数C2 // 对每个像素进行滤波处理 for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { double sum1 = 0.0; double sum2 = 0.0; for (int k = -radius; k <= radius; k++) { for (int l = -radius; l <= radius; l++) { int row = i + k; int col = j + l; if (row >= 0 && row < rows && col >= 0 && col < cols) { double d = sqrt((double)k * k + (double)l * l); // 计算欧式距离 double w1 = exp(-c1 * d * d); // 计算空域权重 double diff = input.at<uchar>(i, j) - input.at<uchar>(row, col); // 计算像素值差异 double w2 = exp(-c2 * diff * diff); // 计算灰度值权重 sum1 += w1 * w2 * input.at<uchar>(row, col); // 计算加权像素值之和 sum2 += w1 * w2; // 计算加权权重之和 } } } output.at<uchar>(i, j) = (uchar)(sum1 / sum2); // 更新输出图像像素值 } } return output; } int main() { Mat input = imread("lena.jpg", IMREAD_GRAYSCALE); // 读取灰度图像 if (input.empty()) { cout << "读取图像失败!" << endl; return -1; } Mat output = adaptive_bilateral_filter(input, 11, 10.0, 10.0); // 自适应双边滤波 imshow("input", input); imshow("output", output); waitKey(0); return 0; } ``` 其中,`adaptive_bilateral_filter()`函数是自适应双边滤波的具体实现函数,输入参数包括原始图像`input`、窗口大小`size`、空域标准差`sigma_d`和灰度值标准差`sigma_r`,输出参数为滤波后的图像`output`。 在`main()`函数中,首先读取灰度图像`input`,然后调用`adaptive_bilateral_filter()`函数进行自适应双边滤波,并将结果保存在`output`中,最后显示原始图像和滤波后的图像。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 17
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值