MATLAB---非线性复扩散滤波

%

非线性复扩散滤波

clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic
[filename ,pathname]=…
uigetfile({‘.bmp’;'.tif’;‘*.jpg’;},‘选择图片’); %选择图片路径
str=[pathname filename]; % 合成路径+文件名
im = imread(str); % 读图
im = imnoise(im,‘gaussian’,0,1e-3); % 原图像 + 白噪声

%非线性复扩散滤波参数设置
TMAX = 0.80; % 扩散时间
[im_e, nIter, dTT] = NCD_filter(im, TMAX); % 非线性复扩散滤波

figure,
subplot(121),imshow(im);title(‘原始图像’)
colormap(jet) % 颜色
shading interp % 消隐
subplot(122),imshow(im_e,[]);title(‘非线性复扩散滤波图像’)
colormap(jet) % 颜色
shading interp % 消隐
toc
function [imgOutput, nIter, dtt] = NCD_filter(imgInput, Tmax)
%　非线性复扩散滤波
% 函数输入：
% imgInput - 含噪声的图像
% Tmax - 扩散时间
% 函数输出:
% imgOutput - f非线性复扩散滤波图像
% nIter - 滤波迭代处理次数
% dtt - 每一次迭代的时间步长
% 设定默认的迭代扩散时间
if nargin < 2 ; % 输入个数小于2
Tmax = 3.0 ; % Tmax默认赋值
end

% 初始化操作
theta = pi/30; % 初始化
j = sqrt(-1); % 初始化
e_jxtheta = exp(j*theta); % 初始化
kappamin = 2.0; % 初始化
kappamax = 28.0; % 初始化

% 高斯滤波器
wsize = 3; % 窗口大小 3 x 3
wsigma = 10; % 方差
filter_gaussian = fspecial(‘gaussian’, wsize, wsigma); % 滤波掩模
% fspecial(‘gaussian’, 3, 0.001)
% ans =
% 0 0 0
% 0 1 0
% 0 0 0
% fspecial(‘gaussian’, 5, 0.001)
% ans =
% 0 0 0 0 0
% 0 0 0 0 0
% 0 0 1 0 0
% 0 0 0 0 0
% 0 0 0 0 0
% 扩散滤波器系数
wsigma2 = 0.5; % 方差
filter_gaussian2 = fspecial(‘gaussian’, wsize, wsigma2); % 滤波掩模

lapKernel = [0,1,0; 1,-4,1; 0,1,0]; % laplacian kernel
gradKernel = [-1/2,0,1/2]; % 梯度kernel
[xSize, ySize] = size(imgInput); % 图像维数

Border = 2; % 图像边界2个像素点之间不进行梯度计算
indexX = 1+Border:xSize+Border; % (1+Border):(xSize+Border)
indexY = 1+Border:ySize+Border; % (1+Border):(xSize+Border)

if ~isa(imgInput,‘double’)
yNCDF = double(imgInput); % 图像数据类型转换
end

t_iter = 0; % 迭代时间累加和
nIter = 0; % 迭代次数

while (Tmax - t_iter) > eps % 扩散时间
nIter = nIter + 1;
ryNCDF = real(yNCDF); % 实部
iyNCDF = imag(yNCDF); % 虚部

% 滤波,见表达式（10.32）
localAvg = filter2(filter_gaussian, ryNCDF, 'same');
minlocalAvg = min(localAvg(:));   % 最小值

% 自适应系数 k，见表达式（10.31）
k_mod = kappamax + (kappamin-kappamax)* ...
    (localAvg - minlocalAvg) / (max(localAvg(:)) - minlocalAvg);
% 非线性复扩散函数
coefDif = e_jxtheta ./ (1 + ( (iyNCDF/theta) ./ k_mod ).^2);
coefDif = filter2(filter_gaussian, coefDif, 'same');

% 计算 laplacian and gradient
% lap(yNCDF)
yAux  = zeros(xSize+2*Border, ySize+2*Border);
yAux(indexX, indexY) = yNCDF;
del2Y = conv2(yAux, lapKernel, 'same');
del2Y = del2Y(indexX, indexY);

% grad(yNCDF)
dAux  = conv2(yAux, gradKernel, 'same');  % 卷积
dYx   = dAux(indexX, indexY);
dAux  = conv2(yAux, gradKernel', 'same'); % 卷积
dYy   = dAux(indexX, indexY);

% grad(coefDif)
dDx   = conv2(coefDif, gradKernel, 'same');  % 卷积
dDy   = conv2(coefDif, gradKernel', 'same'); % 卷积

dIdt  = coefDif.*del2Y + dDx.*dYx + dDy.*dYy;

% 计算自适应时间步长
dt = (1/4)*( 0.25 + 0.75*exp(- max(max( abs(real(dIdt)) ./ ryNCDF ))) );

dtt(nIter) = dt;

% 约束每一步的最大步长
if t_iter + dt > Tmax ;
    dt = Tmax - t_iter ; 
end
% 更新已经处理已经扩散的时间消耗
t_iter = t_iter + dt;

% 更新 yNCDF
yNCDF = yNCDF + dt.*dIdt;

end % 结束，对应while (Tmax - t_iter) > eps

imgOutput = real(yNCDF); % 实部，图像输出

end