一、簡介

基於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
完整程式碼或代寫加1564658423