matlab-图像重建matlab仿真,包括标准正则化,空域迭代正则化方法, 基于噪声的自适应正则化方法-源码

%**********************************% %* 空域迭代正则化方法的图像重建 *% %**********************************% clc; clear; close all; warning off; addpath(genpath(pwd)); I =imread('x.png'); J = uint8(255*ones(size(rgb2gray(I))))-rgb2gray(I); [rr,cc]=size(J); % J2 = zeros(16*rr,16*cc); % for i = 1:rr % for j = 1:cc % J2(16*i-16+1:16*i,16*j-16+1:16*j) = J(i,j); % end % end I=uint8(J); figure(1);subplot(2,3,1); imshow(I); title('原始图像'); %读入图像，预处理 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %高斯模糊降质函数 siz = (size(I)-1)/2; std = 1; [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1)); arg = -((x.*x + y.*y)/(std*std)).^1.5; H = (exp(arg)); H(H<eps*max(H(:))) = 0; sumh = sum(H(:)); if sumh ~= 0, H = H/sumh; end; %或者 %H=fspecial('gaussian',10,2); %产生高斯模糊降析函数 Gaussian = imfilter(I,H,'replicate'); figure(1);subplot(2,3,2); imshow(Gaussian);title('高斯模糊图像'); %显示模糊图像 %噪声 Gaussian_G=imnoise(Gaussian,'gaussian',0,0.001);%高斯 %Gaussian_G=imnoise(Gaussian,'salt & pepper',0.005);%椒盐 figure(1);subplot(2,3,3); imshow(Gaussian_G);title('加入噪声图像'); %显示加入噪声图像 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blurimage=Gaussian_G; [row,col,raw]=size(blurimage); %读入观测图像 %中心化算子 N1=row; N2=col; N3=raw; T=zeros(N1,N2); x0=N1/2;y0=N2/2; for re=1:N1 for im=1:N2 ta=cos( 2*pi*( (re-1)*x0+(im-1)*y0 )/N1 ); tb=-sin(2*pi*( (re-1)*x0+(im-1)*y0 )/N2 ); T(re,im)=complex(ta,tb); end end %逆滤波 p=[0 -1 0;-1 4 -1;0 -1 0]; %拉普拉斯算子，实现二阶微分 FFT_p=fft2(p,row,col); FFT_blurimage=fft2(blurimage,row,col); FFT_G=fft2(Gaussian_G,row,col); FFT_H=fft2(H,row,col); FFT_H=FFT_H.*T; CONJ_FFT_H=conj(FFT_H); det_FFT_H=abs(FFT_H); det_FFT_p=abs(FFT_p); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X1=Gaussian_G; %X1获取高分辨率图像 d=0.000001; %迭代中止条件 m=0.01; %权值 c=0.00001; %移位 Y=Gaussian_G; %Y获取低分辨率图像 Q=[0,-1,0;-1,4,-1;0,-1,0]; mark=0; %迭代过程 for k=1:20 % 求正则化参数 alfa A =H; B =imfilter(Y,Q,'replicate') N =imfilter(X1,A,'replicate'); a1=norm((double(Y)-double(N)),2).^2; a2=norm(double(Y),2).^2; %a3=norm(double(X1),2).^2; alfa=log(m*a1/(a2)+1); %参数计算方法，在不同算法中选择不同 %alfa=log(m*a1/(a2+c)+1); %alfa=(1.5*a1)/a2; xx(k)=k; %迭代曲线坐标 yy(k)=alfa; %迭代形式 FFT_X1=fft2(X1,row,col); FFT_X1_t=CONJ_FFT_H.*FFT_X1; %迭代初始值 FFT_object=FFT_X1_t+CONJ_FFT_H.*FFT_blurimage-FFT_X1_t.*(det_FFT_H.^2+alfa*det_FFT_p.^2); %FFT_object=(CONJ_FFT_H.*FFT_blurimage)./(det_FFT_H.^2+alfa*det_FFT_p.^2); object=abs(ifft2(FFT_object)); %迭代结果缓存 object_max=max(max(object)); %标定 object_min=min(min(object)); object=(object-object_min)*255/(object_max-object_min); object = uint8(object); X2=object; t=X1; X1=X2; diff=(norm(double(X2)-double(t),2).^2)/((norm(double(t),2)).^2) if (mark==0) if (diff <=d) mark=k; resultim=object; end %满足迭代终止条件，取重建结果，此处为绘制迭代曲线并不终止迭代 end end if (mark==0)%迭代次数过多 mark=50; resultim=object; end %计算信噪比 [L1,L2]=size(I) dx=norm((double(resultim)-double(I)),2).^2 PSNR=10*log((255*255*L1*L2)/dx) %显示结果 subplot(2,3,4); figure(1);imshow(resultim,[0,255]); promp=['重建图像',10,'迭代次数为' num2str(mark),10,'PSNR=' num2str(PSNR)]; title(promp); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f=resultim; %锐化相加，消除寄生波纹 f=double(f); [m,n]=size(f); %读取正则化重建结果 g=f; g1=f; g0=f; %拉普拉斯锐化 for i=2:m-1 for j=2:n-1 g00(i,j)=g(i+1,j)+g(i-1,j)+g(i,j-1)+g(i,j+1)-4*g(i,j); %+g(i+1,j+1)+g(i+1,j-1)+g(i-1,j+1)+g(i-1,j-1) end end for i=1:m-2 for j=1:n-2 g0(i,j)=f(i,j)-g00(i,j); end end g0=uint8(g0); subplot(2,3,5); figure(1);imshow(g0);title('拉普拉斯锐化后图像'); %显示锐化图像 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %结果 finalim=uint8((double(resultim)+(double(g0)-double(resultim))/2)); %计算信噪比 [L1,L2]=size(I) dx=norm((double(finalim)-double(I)),2).^2 PSNR1=10*log((255*255*L1*L2)/dx) promp1=['消除寄生波纹图像',10,'PSNR=' num2str(PSNR1)]; figure(1);subplot(2,3,6); imshow(finalim);title(promp1); %显示结果 %绘制迭代曲线 figure(2); plot(xx,yy); axis([0 mark+5 0 max(yy)+0.01]) ; promp=['正则化参数迭代曲线']; title(promp); imwrite(resultim,'result01.tif','tif'); %输出结果。