一、简介
基于matlab Wiener+Non-Local Means+Lucy_Richardson+Lee+kuwahara+Bilateral滤波器图像去噪
二、源代码
function [out, psn]=bif_filter(im,sigd,sigr)
% bilateral filter双边滤波器
% 函数输入:
% im 输入的图像
% sigd 空间内核的时域参数
% sigr 内核参数强度变化范围
% 函数输出:
% out 滤波图像 = output imagespatial kernel
w=(2*sigd)+1;
% sigr=(n*100)^2/(.003*(sigd^2)); % 自适应R值,n为高斯噪声强度,n=0.001
% 高斯滤波器
gw=zeros(w,w); % 高斯权值矩阵初始化
c=ceil(w/2); % 向前取整
c=[c c]; % 中心元素位置
for i=1:w
for j=1:w
q=[i,j]; % 记录相连像素位置标识位
gw(i,j)=norm(c-q); % 欧氏距离
end
end
Gwd=(exp(-(gw.^2)/(2*(sigd^2)))); % 高斯函数
% Padding 扩展图像的边界,防止滑动窗口边界值溢出
proci=padarray(im,[sigd sigd],'replicate');
% A = [1 2; 3 4];
% B = padarray(A,[3 2],'replicate','post')
% B =
% 1 2 2 2
% 3 4 4 4
% 3 4 4 4
% 3 4 4 4
% 3 4 4 4
[row clm]=size(proci); % Size of image
if ~isa(proci,'double')
proci = double(proci)/255; % 转换为double类型
end
K=sigd;
L=[-K:K];
c=K+1; % 中心元素位置
iter=length(L); % 迭代次数
ind=1;
for r=(1+K):(row-K) % 行
for s=(1+K):(clm-K) % 列
for i=1:iter % 窗口大小 行
for j=1:iter % 窗口大小 列
win(i,j)=proci((r+L(i)),(s+L(j))); % 获取窗口
end
end
I=win; % 灰度矩阵
win=win(c,c)-win; % 相对中心点处的强度差异,中心点为参考灰度值
win=sqrt(win.^2); % 保证win中的每一个元素为正
Gwi=exp(-(win.^2)/(2*(sigr^2))); % 高斯函数
weights=(Gwi.*Gwd)/sum(sum(Gwi.*Gwd)); % 高斯权值
Ii=sum(sum(weights.*I)); % 得到当前双边滤波值
proci(r,s)=Ii; % 替换当前灰度值
win=[];
end
end
% 移除边界扩展值
proci=rpadd(proci,K); % 移除边界扩展值
out=im2uint8(proci); % 类型转换
%% 滤波重建后,图像峰值信噪比计算
if ~isa(out,'double')
dimg = double(out)/255; % 转换为double类型
end
psn = PSN(im,dimg); % PSNR,峰值信噪比
end
% PSF: 退化函数的空域模板
% NP: 表示噪声的功率
%函数输出:
% J: 约束最小平方滤波图像
% LAGRA: 可以为一个数值,表示指定约束最小平方的最佳复原参数y,
% 也可以为[min,max]形式,表示参数y的搜索范围
% 若此参数省略,则表示搜索范围为[1e-9,1e9] 。
% 约束最小平方滤波
if ~isa(I,'double')
I = double(I)/255;
end
LR = [1e-9 1e9]; % 复原参数搜索范围
% 求解输入图像维数
sizeI = size(I);
% PSF 矩阵
sizePSF = size(PSF);
% 输入图像的维数求解
numNSdim = find(sizePSF~=1);
NSD = length(numNSdim);
% 转化PSF函数到期望的维数 光传递函数OTF
otf = psf2otf(PSF,sizeI);
% regop:通过计算拉普拉斯算子计算图像的平滑性
% 具体见表达式(10.23)
if NSD == 1,
regop = [1 -2 1];
else % 二维矩阵
% 3x3 Laplacian matrixes
regop = repmat(zeros(3),[1 1 3*ones(1,NSD-2)]);
% center matrix of Laplacian
idx = [{':'}, {':'}, repmat({2},[1 NSD-2])];
regop(idx{:}) = [0 1 0;1 -NSD*2 1;0 1 0]; % 模板算子
end
% regop =
% 0 1 0
% 1 -4 1
% 0 1 0
% 改变REGOP折返回原始维数
idx1 = repmat({1},[1 length(sizePSF)]);
idx1(numNSdim) = repmat({':'},[1 NSD]);
REGOP(idx1{:}) = regop;
% 转化PSF函数到期望的维数 光传递函数OTF
REGOP = psf2otf(REGOP,sizeI);
fftnI = fftn(I);
R2 = abs(REGOP).^2;
H2 = abs(otf).^2;
% 计算LAGRA值
if (numel(LR)==1) || isequal(diff(LR),0),% LAGRA is given
LAGRA = LR(1);
else % 采用fminbnd在[1e-9,1e9]优化,加速计算
R4G2 = (R2.*abs(fftnI)).^2;
H4 = H2.^2;
R4 = R2.^2;
H2R22 = 2*H2.*R2;
ScaledNP = NP*prod(sizeI);
LAGRA = fminbnd(@ResOffset,LR(1),LR(2),[],R4G2,H4,R4,H2R22,ScaledNP);
end;
三、运行结果
四、备注
版本:2014a