非线性复扩散滤波

Halcom

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

1. function [imgOutput, nIter, dtt] = NCD_filter(imgInput, Tmax)
   
2. %　非线性复扩散滤波
   
3. % 函数输入：
   
4. %        imgInput - 含噪声的图像
   
5. %        Tmax  - 扩散时间
   
6. % 函数输出:
   
7. %        imgOutput - f非线性复扩散滤波图像
   
8. %        nIter  - 滤波迭代处理次数
   
9. %        dtt    - 每一次迭代的时间步长
   
10. % 设定默认的迭代扩散时间
    
11. if nargin < 2 ;   % 输入个数小于2
    
12.     Tmax  = 3.0 ; % Tmax默认赋值
    
13. end
    
14. 
    
15. % 初始化操作
    
16. theta     = pi/30;           % 初始化
    
17. j         = sqrt(-1);        % 初始化
    
18. e_jxtheta = exp(j*theta);    % 初始化
    
19. kappamin  =   2.0;           % 初始化
    
20. kappamax  =  28.0;           % 初始化
    
21. 
    
22. % 高斯滤波器
    
23. wsize     = 3;  % 窗口大小 3 x 3
    
24. wsigma    = 10; % 方差
    
25. filter_gaussian    = fspecial('gaussian', wsize, wsigma);  % 滤波掩模
    
26. % fspecial('gaussian', 3, 0.001)
    
27. % ans =
    
28. %      0     0     0
    
29. %      0     1     0
    
30. %      0     0     0
    
31. % fspecial('gaussian', 5, 0.001)
    
32. % ans =
    
33. %      0     0     0     0     0
    
34. %      0     0     0     0     0
    
35. %      0     0     1     0     0
    
36. %      0     0     0     0     0
    
37. %      0     0     0     0     0
    
38. % 扩散滤波器系数
    
39. wsigma2    = 0.5; % 方差
    
40. filter_gaussian2  = fspecial('gaussian', wsize, wsigma2);  % 滤波掩模
    
41. 
    
42. lapKernel      = [0,1,0; 1,-4,1; 0,1,0]; % laplacian kernel
    
43. gradKernel     = [-1/2,0,1/2];          % 梯度kernel
    
44. [xSize, ySize] = size(imgInput);         % 图像维数
    
45. 
    
46. Border = 2; % 图像边界2个像素点之间不进行梯度计算
    
47. indexX = 1+Border:xSize+Border; % (1+Border):(xSize+Border)
    
48. indexY = 1+Border:ySize+Border; % (1+Border):(xSize+Border)
    
49. 
    
50. if ~isa(imgInput,'double')
    
51.     yNCDF = double(imgInput); % 图像数据类型转换
    
52. end
    
53. 
    
54. t_iter = 0;  % 迭代时间累加和
    
55. nIter = 0;  % 迭代次数
    
56. 
    
57. while (Tmax - t_iter) > eps % 扩散时间
    
58.     nIter = nIter + 1;
    
59.     ryNCDF = real(yNCDF);  % 实部
    
60.     iyNCDF = imag(yNCDF);  % 虚部
    
61.    
    
62.     % 滤波,见表达式（10.32）
    
63.     localAvg = filter2(filter_gaussian, ryNCDF, 'same');
    
64.     minlocalAvg = min(localAvg(:));   % 最小值
    
65.    
    
66.     % 自适应系数 k，见表达式（10.31）
    
67.     k_mod = kappamax + (kappamin-kappamax)* ...
    
68.         (localAvg - minlocalAvg) / (max(localAvg(:)) - minlocalAvg);
    
69.     % 非线性复扩散函数
    
70.     coefDif = e_jxtheta ./ (1 + ( (iyNCDF/theta) ./ k_mod ).^2);
    
71.     coefDif = filter2(filter_gaussian, coefDif, 'same');
    
72. 
    
73.     % 计算 laplacian and gradient
    
74.     % lap(yNCDF)
    
75.     yAux  = zeros(xSize+2*Border, ySize+2*Border);
    
76.     yAux(indexX, indexY) = yNCDF;
    
77.     del2Y = conv2(yAux, lapKernel, 'same');
    
78.     del2Y = del2Y(indexX, indexY);
    
79.    
    
80.     % grad(yNCDF)
    
81.     dAux  = conv2(yAux, gradKernel, 'same');  % 卷积
    
82.     dYx   = dAux(indexX, indexY);
    
83.     dAux  = conv2(yAux, gradKernel', 'same'); % 卷积
    
84.     dYy   = dAux(indexX, indexY);
    
85.    
    
86.     % grad(coefDif)
    
87.     dDx   = conv2(coefDif, gradKernel, 'same');  % 卷积
    
88.     dDy   = conv2(coefDif, gradKernel', 'same'); % 卷积
    
89.    
    
90.     dIdt  = coefDif.*del2Y + dDx.*dYx + dDy.*dYy;
    
91.    
    
92.     % 计算自适应时间步长
    
93.     dt = (1/4)*( 0.25 + 0.75*exp(- max(max( abs(real(dIdt)) ./ ryNCDF ))) );
    
94.    
    
95.     dtt(nIter) = dt;
    
96.    
    
97.     % 约束每一步的最大步长
    
98.     if t_iter + dt > Tmax ;
    
99.         dt = Tmax - t_iter ;
    
100.     end
     
101.     % 更新已经处理已经扩散的时间消耗
     
102.     t_iter = t_iter + dt;
     
103.    
     
104.     % 更新 yNCDF
     
105.     yNCDF = yNCDF + dt.*dIdt;
     
106.    
     
107. end % 结束，对应while (Tmax - t_iter) > eps
     
108. 
     
109. imgOutput = real(yNCDF);  % 实部，图像输出
     
110. 
     
111. end

复制代码