本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part6，辑录该书第281至第374页之代码，供有需要读者下载研究使用。代码执行结果请参见原书配图，建议下载代码前阅读下文：关于《数字图像处理原理与实践（MATLAB版）》一书代码发布的说明 P338i=double(imread(‘vase.tif’));[C,S]=wavedec2(i,2,’db1′);a2=appcoef2(C,S,’db1′,2);dh1=detcoef2(‘h’,C,S,1);dv1=detcoef2(‘v’,C,S,1);dd1=detcoef2(‘d’,C,S,1);dh2=detcoef2(‘h’,C,S,2);dv2=detcoef2(‘v’,C,S,2);dd2=detcoef2(‘d’,C,S,2);[x,y]=size(i);img = zeros(x,y);img(1:x/4,1:y/4) =im2uint8(mat2gray(a2));img(((x/4)+1):y/2,1:y/4) = im2uint8(mat2gray(dv2));img(((x/4)+1):x/2,1:y/4) = im2uint8(mat2gray(dv2));img(1:x/4,((y/4)+1):y/2) = im2uint8(mat2gray(dh2));img(((x/4)+1):x/2,((y/4)+1):y/2) = im2uint8(mat2gray(dd2));img(((x/2)+1):x,1:y/2) = im2uint8(mat2gray(dv1));img(1:x/2,((y/2)+1):y) = im2uint8(mat2gray(dh1));img(((x/2)+1):x,((y/2)+1):y) = im2uint8(mat2gray(dd1));imshow(img,[]);P341-1X1 = imread(‘cathe1.bmp’);X2 = imread(‘cathe2.bmp’);XFUS = wfusimg(X1,X2,’sym4′,5,’mean’,’max’);imshow(XFUS,[]);P341-2X1 = imread(‘cathe1.bmp’);X2 = imread(‘cathe2.bmp’);M1 = double(X1) / 256;M2 = double(X2) / 256;N = 4;wtype = ‘sym4’;[c0,s0] = wavedec2(M1, N, wtype);[c1,s1] = wavedec2(M2, N, wtype);length = size(c1);Coef_Fusion = zeros(1,length(2));%低频系数的处理，取平均值Coef_Fusion(1:s1(1,1)) = (c0(1:s1(1,1))+c1(1:s1(1,1)))/2;%处理高频系数，，取绝对值大者，这里用到了矩阵乘法MM1 = c0(s1(1,1)+1:length(2));MM2 = c1(s1(1,1)+1:length(2));mm = (abs(MM1)) > (abs(MM2));Y = (mm.*MM1) + ((~mm).*MM2);Coef_Fusion(s1(1,1)+1:length(2)) = Y;%重构Y = waverec2(Coef_Fusion,s0,wtype);imshow(Y,[]);P344I = imread(‘noise_lena.bmp’);[thr,sorh,keepapp] = ddencmp(‘den’,’wv’,I);de_I = wdencmp(‘gbl’,I,’sym4′,2,thr,sorh,keepapp);imwrite(im2uint8(mat2gray(de_I)), ‘denoise_lena.bmp’);P361function diff_im = anisodiff(im, num_iter, delta_t, k, option)im = double(im);% 赋初值diff_im = im;% 用以计算方向梯度的卷积模板hN = [0 1 0; 0 -1 0; 0 0 0];hS = [0 0 0; 0 -1 0; 0 1 0];hE = [0 0 0; 0 -1 1; 0 0 0];hW = [0 0 0; 1 -1 0; 0 0 0];hNE = [0 0 1; 0 -1 0; 0 0 0];hSE = [0 0 0; 0 -1 0; 0 0 1];hSW = [0 0 0; 0 -1 0; 1 0 0];hNW = [1 0 0; 0 -1 0; 0 0 0];% 各向异性扩散滤波for t = 1:num_iter% 计算梯度 nablaN = conv2(diff_im,hN,’same’); nablaS = conv2(diff_im,hS,’same’); nablaW = conv2(diff_im,hW,’same’); nablaE = conv2(diff_im,hE,’same’); nablaNE = conv2(diff_im,hNE,’same’); nablaSE = conv2(diff_im,hSE,’same’); nablaSW = conv2(diff_im,hSW,’same’); nablaNW = conv2(diff_im,hNW,’same’); % 计算扩散系数 % OPTION 1： c(x,y,t) = exp(-(nablaI/kappa).^2) if option == 1 cN = exp(-(nablaN/k).^2); cS = exp(-(nablaS/k).^2); cW = exp(-(nablaW/k).^2); cE = exp(-(nablaE/k).^2); cNE = exp(-(nablaNE/k).^2); cSE = exp(-(nablaSE/k).^2); cSW = exp(-(nablaSW/k).^2); cNW = exp(-(nablaNW/k).^2); % OPTION 2： c(x,y,t) = 1./(1 + (nablaI/kappa).^2) elseif option == 2 cN = 1./(1 + (nablaN/k).^2); cS = 1./(1 + (nablaS/k).^2); cW = 1./(1 + (nablaW/k).^2); cE = 1./(1 + (nablaE/k).^2); cNE = 1./(1 + (nablaNE/k).^2); cSE = 1./(1 + (nablaSE/k).^2); cSW = 1./(1 + (nablaSW/k).^2); cNW = 1./(1 + (nablaNW/k).^2); end % 计算一次迭代结果 diff_im = diff_im + delta_t*(… cN.*nablaN + cS.*nablaS + cW.*nablaW + cE.*nablaE + … cNE.*nablaNE + cSE.*nablaSE + cSW.*nablaSW + cNW.*nablaNW );endP363num_iter=50; delta_t=0.125;k=4; option=2;i = imread(‘noise_lena.bmp’);diff = anisodiff(i, num_iter, delta_t, k, option);P370function x=Thomas(N, alpha, beta, gama, d)x=d;m=zeros(1,N); l=zeros(1,N);m(1)=alpha(1);for i=2:N l(i)=gama(i)/m(i-1); m(i)=alpha(i)-l(i)*beta(i-1);endy=zeros(1,N);y(1)=d(1);for i=2:N y(i)=d(i)-l(i)*y(i-1);endx=zeros(1,N);x(N)=y(N)/m(N);for i=N-1:-1:1 x(i)=(y(i)-beta(i)*x(i+1))/m(i);endP371function Ig=gauss(I,ks,sigma2)[Ny,Nx]=size(I);hks=(ks-1)/2; % 高斯核的一半%%- 一维卷积if (Ny<ks) x=(-hks:hks); flt=exp(-(x.^2)/(2*sigma2)); % 一维高斯函数 flt=flt/sum(sum(flt)); % 归一化 x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx)); eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)]; Ig=conv(eI,flt); Ig=Ig(:,ks+hks+1:Nx+ks+hks);else %%- 二维卷积 x=ones(ks,1)*(-hks:hks); y=x’; flt=exp(-(x.^2+y.^2)/(2*sigma2)); % 二维高斯函数 flt=flt/sum(sum(flt)); % 归一化 if (hks>1) xL=mean(I(:,1:hks)’)’; xR=mean(I(:,Nx-hks+1:Nx)’)’; else xL=I(:,1); xR=I(:,Nx); end eI=[xL*ones(1,hks) I xR*ones(1,hks)]; if (hks>1) xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:)); else xU=eI(1,:); xD=eI(Ny,:); end eI=[ones(hks,1)*xU; eI; ones(hks,1)*xD]; Ig=conv2(eI,flt,’valid’);endP372Img = imread(‘Lena.bmp’);Img = double(Img);[nrow, ncol] = size(Img);N=max(nrow, ncol);%储存三对角矩阵alpha=zeros(1,N); beta=zeros(1,N); gama=zeros(1,N);%储存中间结果u1=zeros([nrow, ncol]);u2=zeros([nrow, ncol]);timestep=5;%用以控制迭代次数%iterations = 2;%for times = 1:iterations I_temp=gauss(Img,3,1); Ix = 0.5*(I_temp(:,[2:ncol,ncol])-I_temp(:,[1,1:ncol-1])); Iy = 0.5*(I_temp([2:nrow,nrow],:)-I_temp([1,1:nrow-1],:)); K = 10 grad=Ix.^2+Iy.^2; g=1./(1+grad/K*K); %边缘压迫因子 % 使用Thomas算法逐行求解u1 for i=1:nrow beta(1)=-0.5*timestep*(g(i,2)+g(i,1)); alpha(1)=1-beta(1); for j=2:ncol-1 beta(j)=-0.5*timestep*(g(i,j+1)+g(i,j)); gama(j)=-0.5*timestep*(g(i,j-1)+g(i,j)); alpha(j)=1-beta(j)-gama(j); end gama(ncol)=-0.5*timestep*(g(i,ncol)+g(i,ncol-1)); alpha(ncol)=1- gama(ncol); u1(i,:)=Thomas(ncol,alpha,beta,gama,Img(i,:)); end % 使用Thomas算法逐列求解u2 for j=1:ncol beta(1)=-0.5*timestep*(g(2,j)+g(1,j)); alpha(1)=1-beta(1); for i=2:nrow-1 beta(j)=-0.5*timestep*(g(i+1,j)+g(i,j)); gama(j)=-0.5*timestep*(g(i-1,j)+g(i,j)); alpha(j)=1-beta(j)-gama(j); end gama(nrow)=-0.5*timestep*(g(nrow,j)+g(nrow-1,j)); alpha(nrow)=1- gama(nrow); u2(:,j)=Thomas(nrow,alpha,beta,gama,Img(:,j)); end Img=0.5*(u1+u2); % 显示处理结果 imshow(uint8(Img));%end（代码发布未完，请待后续…）
人，总是很难改正自己的缺点，
