《数字图像处理原理与实践(MATLAB版)》一书之代码Part5

时间:2022-11-04 18:46:47

本文系《数字图像处理原理与实践(MATLAB版)》一书之代码系列的Part5,辑录该书第281至第405页之代码,供有须要读者下载研究使用。代码运行结果请參见原书配图。另外,书中实验用图片文件及其他资源可以从图像处理书籍读者QQ群中获取(入群通道请见文末)。


P338

i=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-1

X1 = imread('cathe1.bmp');X2 = imread('cathe2.bmp');XFUS = wfusimg(X1,X2,'sym4',5,'mean','max');imshow(XFUS,[]);

P341-2

X1 = 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,[]);

P344

I = 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');

P361

function 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 );end

P363

num_iter=50; delta_t=0.125;k=4; option=2;i = imread('noise_lena.bmp');diff = anisodiff(i, num_iter, delta_t, k, option);

P370

function 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);end

P371

function 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');end

P372

Img = 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

P385-1

function y = tv(X)    [M,N] = size(X);    Dh = diff(X,[],1);    Dh = [Dh;zeros(1,N)];    Dv = diff(X,[],2);    Dv = [Dv zeros(M,1)];    y = sum(sum(sqrt(Dh.^2+Dv.^2)));end


P385-2

function y = atv(X)    [M,N] = size(X);    Dh = -diff(X,[],1);    Dh = [Dh;zeros(1,N)];    Dv = -diff(X,[],2);    Dv = [Dv zeros(M,1)];    y = sum(sum(abs(Dh)+abs(Dv)));end


P387

I = double(rgb2gray(imread('Lena.bmp')));I0 = I;ep=1; dt=0.25; lam=0;ep2=ep^2; [ny,nx]=size(I);iter = 80;for i=1:iter,    % 中心差法计算梯度和微分    % WN  N  EN    % W   O  E    % WS  S  ES    I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2; % Ix = (E-W)/2    I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2; % Iy = (S-N)/2    I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I; % Ixx = E+W-2*O    I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I; % Iyy = S+N-2*O    Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);    Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);    I_xy = (Dp-Dm)/4; % Ixy = Iyx = ((ES+WN)-(EN+WS))/4        Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);    Den = (ep2+I_x.^2+I_y.^2).^(3/2);    I_t = Num./Den + lam.*(I0-I);    I=I+dt*I_t;  %% evolve image by dtendimshow(I,[]);

P400

RGB = imread('moon.jpg');I = rgb2gray(RGB);J = imnoise(I,'gaussian',0,0.025);imshow(J)K = wiener2(J,[5 5]);figure, imshow(K)

P401

I = im2double(imread('cameraman.tif'));imshow(I), title('Original Image');%模拟运动模糊,生成一个点扩散函数PSF,相应的线性运动长度为21个像素,角度为11LEN = 21;THETA = 11;PSF = fspecial('motion', LEN, THETA);blurred = imfilter(I, PSF, 'conv', 'circular');figure, imshow(blurred), title('Blurred Image');%对运动模糊图像进行复原result_w1= deconvwnr(blurred, PSF);figure, imshow(result_w1), title('Restoration of Blurred Image')%模拟加性噪声noise_mean = 0;noise_var = 0.0001;blurred_noisy = imnoise(blurred, 'gaussian', noise_mean, noise_var);figure, imshow(blurred_noisy), title('Simulate Blur and Noise')%假设没有噪声的情况下复原图像estimated_nsr = 0;wnr2 = deconvwnr(blurred_noisy, PSF, estimated_nsr);figure, imshow(wnr2)title('Restoration of Blurred, Noisy Image Using NSR = 0')%使用一个更好的信噪功率比评估值来复原图像estimated_nsr = noise_var / var(I(:));wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);figure, imshow(wnr3)title('Restoration of Blurred, Noisy Image Using NSR');

P404

%载入原始图像I = imread('board.tif');I = I(50+[1:256],2+[1:256],:);figure, imshow(I)title('Original Image')%模糊处理PSF = fspecial('gaussian',5,5);Blurred = imfilter(I,PSF,'symmetric','conv');figure, imshow(Blurred);title('Blurred Image')%添加噪声V = 0.002;BlurredNoisy = imnoise(Blurred,'gaussian',0,V);figure, imshow(BlurredNoisy)title('Blurred and Noisy Image')

P405

%利用露茜-理查德森算法复原图像(5次迭代)luc1 = deconvlucy(BlurredNoisy, PSF, 5);figure, imshow(luc1)%利用露茜-理查德森算法复原图像luc1_cell = deconvlucy({BlurredNoisy}, PSF, 5);luc2_cell = deconvlucy(luc1_cell, PSF);luc2 = im2uint8(luc2_cell{2});figure, imshow(luc2);%控制阈值的复原效果DAMPAR = im2uint8(3*sqrt(V));luc3 = deconvlucy(BlurredNoisy, PSF, 15, DAMPAR);figure, imshow(luc3);



(代码发布未完,请待后续...)


《数字图像处理原理与实践(MATLAB版)》一书之代码Part5