《数字图像处理原理与实践（MATLAB版）》一书之代码Part6

时间：2014-12-17 00:23:01 阅读：180 评论：0 收藏：0 [点我收藏+]

本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part6，辑录该书第281至第374页之代码，供有需要读者下载研究使用。代码执行结果请参见原书配图，建议下载代码前阅读下文：

关于《数字图像处理原理与实践（MATLAB版）》一书代码发布的说明

http://blog.csdn.net/baimafujinji/article/details/40987807

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);
end
y=zeros(1,N);
y(1)=d(1);
for i=2:N
   y(i)=d(i)-l(i)*y(i-1);
end

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

（代码发布未完，请待后续...）

标签：图像处理 matlab aos perona-malik方程

原文地址：http://blog.csdn.net/baimafujinji/article/details/41971057

踩

(0)

评论一句话评论（0）

分享档案

更多>

2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)

周排行