基于小波变换的图像修复算法

clear close all clc %************************************************************************** %************************************************************************** x=double(imread('F:\课题有关\张闻霞\张文霞设计资料\张文霞程序\基于小波变换的图像修复算法\LENA.bmp')); figure,imshow(uint8(x)); X=x; %N=size(x); %************************************************************************** %************************************************************************* c=0; L =1;% 小波分解次数 R11=300;%分解前受损区域 R21=310; C11=230; C21=340; R12=450;%分解前受损区域 R22=460; C12=260; C22=360; R13=80;%分解前受损区域 R23=100; C13=80; C23=110; R14=330;%分解前受损区域 R24=400; C14=340; C24=360; R15=360;%分解前受损区域 R25=380; C15=410; C25=480; R16=16;%分解前受损区域 R26=32; C16=300; C26=380; F=zeros(512,512); %************************************************************************** %对图像x进行人工破损 %************************************************************************** for i=R11:R21 for j=C11:C21 x(i,j)=0; F(i,j)=1; c=c+1; end end for i=R12:R22 for j=C12:C22 x(i,j)=0; F(i,j)=1; c=c+1; end end for i=R13:R23 for j=C13:C23 x(i,j)=255; F(i,j)=1; c=c+1; end end for i=R14:R24 for j=C14:C24 x(i,j)=255; F(i,j)=1; c=c+1; end end for i=R15:R25 for j=C15:C25 x(i,j)=0; F(i,j)=1; c=c+1; end end for i=R16:R26 for j=C16:C26 x(i,j)=0; F(i,j)=1; c=c+1; end end %************************************************************************** %************************************************************************** c x3=x; tic figure,imshow(uint8(x)); % forward transform %[af, sf] = farras; [lo,hi]=wfilters('db1'); af=[lo',hi']; lo=lo';hi=hi'; sf=[lo(end:-1:1,:),hi(end:-1:1,:)]; W = dwt2D(x,L,af); % figure,imshow(W{L}{1},[]); % figure,imshow(W{L}{2},[]); % figure,imshow(W{L}{3},[]); k5=6; MSE=0; MSE1=0; x0=W{L+1}; N=size(x0); dist=zeros(N(1),N(2)); dirc_real=zeros(N(1),N(2)); dirc_im=zeros(N(1),N(2)); f=zeros(N(1),N(2)); T=zeros(N(1),N(2)); R1=R11;%分解前受损区域 R2=R21; C1=C11; C2=C21; while(k5~=0) k5=k5-1; min1=10000;% 寻找T值最小的像素 k=1;%判断破损区域是否修复完毕 m=0; n=0;%锁定T值最小像素的坐标 p=1;%窗口大小调整参数 sum1=0; sum2=0; sum_Ex0n=0; sum_Ex0d=0; sum_Ey0n=0; sum_Ey0d=0; c2=0;%判断循环次数 %************************************************************************** % 去除高频子带中的受损边界 %************************************************************************** W{L}{1}((fix(R1/(2^L))-4:fix(R2/(2^L)))+4,(fix(C1/(2^L))-4:fix(C2/(2^L))+4))=0; W{L}{2}((fix(R1/(2^L))-4:fix(R2/(2^L)))+4,(fix(C1/(2^L))-4:fix(C2/(2^L))+4))=0; W{L}{3}((fix(R1/(2^L))-4:fix(R2/(2^L)))+4,(fix(C1/(2^L))-4:fix(C2/(2^L))+4))=0; %figure,imshow(W{L}{1},[]); %figure,imshow(W{L}{2},[]); %figure,imshow(W{L}{3},[]); %************************************************************************** %************************************************************************** %************************************************************************** %修复低频子带 %************************************************************************** f(fix(R1/(2^L)):fix(R2/(2^L)),fix(C1/(2^L)):fix(C2/(2^L)))=255; T(fix(R1/(2^L)):fix(R2/(2^L)),fix(C1/(2^L)):fix(C2/(2^L)))=1e6; f(fix(R1/(2^L))-1,fix(C1/(2^L))-1:fix(C2/(2^L))+1)=1; f(fix(R2/(2^L))+1,fix(C1/(2^L))-1:fix(C2/(2^L))+1)=1; f(fix(R1/(2^L))-1:fix(R2/(2^L))+1,fix(C2/(2^L))+1)=1; f(fix(R1/(2^L))-1:fix(R2/(2^L))+1,fix(C1/(2^L))-1)=1; while(k~=0) k=0; c2=c2+1; min1=1000; %************************************************************************** %look for the least T %************************************************************************** for i=fix(R1/(2^L))-1:fix(R2/(2^L))+1 for j=fix(C1/(2^L))-1:fix(C2/(2^L))+1 if((f(i,j)==1)&(T(i,j)<min1)) min1=T(i,j); k=k+1; m=i; n=j; end end end f(m,n)=0; %************************************************************************** % advance boundary %************************************************************************** if(k~=0) i=m-1; j=n; if (f(i,j)~=0) if(f(i,j)==255) f(i,j)=1; sum1=0; sum2=0; sum_Ex0=0; sum_Ey0=0; sum_Ex0n=0; sum_Ex0d=0; sum_Ey0n=0; sum_Ey0d=0; r=0; for k=i-p:i+p for l=j-p:j+p if f(k,l)==0 %Gx=x0(k,l)-x0(k-1,l-1); %Gy=x0(k-1,l)-x0(k,l-1); Gx=x0(k-1,l+1)-2*x0(k,l-1)-x0(k+1,l-1)-x0(k-1,l-1)+2*x0(k,l+1)+x0(k+1,l+1); Gy=x0(k+1,l-1)-x0(k-1,l-1)+2*x0(k+1,l)-2*x0(k-1,l)+x0(k+1,l+1)-x0(k-1,l+1); Ex=-Gy; Ey=Gx; d0=1/(1+sqrt((k-i)^2+(l-j)^2)); sum_Ex0n=sum_Ex0n+d0*Ex; sum_Ex0d=sum_Ex0d+d0; sum_Ey0n=sum_Ey0n+d0*Ey; sum_Ey0d=sum_Ey0d+d0; dist(k,l)=d0; Ex=Ex/(1+sqrt(Ex^2+Ey^2)); Ey=Ey/(1+sqrt(Ex^2+Ey^2)); dirc_real(k,l)=Ex; dirc_im(k,l)=Ey; end end end Ex0=sum_Ex0n/sum_Ex0d;% (i,j)的等照度线方向。 Ey0= sum_Ey0n/sum_Ey0d; %************************************************************************** %归一化 %************************************************************************** Ex0=Ex0/(1+sqrt(Ex0^2+Ey0^2)); Ey0=Ey0/(1+sqrt(Ex0^2+Ey0^2)); %************************************************************************** %求出权重W(k,l) %************************************************************************** for k=i-p:i+p for l=j-p:j+p if f(k,l)==0 Q=exp(Ex0*dirc_real(k,l)+Ey0* dirc_im(k,l));%cosQ w(k,l)=Q*dist(k,l)*(1/(1+T(k,l))); end end end %************************************************************************** %估计(i,j)灰度值 %************************************************************************** for k=i-p:i+p for l=j-p:j+p if f(k,l)==0 sum1=sum1+x0(k,l)*w(k,l); sum2=sum2+w(k,l); end end end x0(i,j)=sum1/sum2; end %************************************************************************** %解扩散方程 %************************************************************************ %************************************************************************** %sol=solve(i-1,j,i,j-1,f,T);Sol(1)=sol; %************************************************************************* Sol(1)=1e+6; s=0; r=0; d=0; if(f(i-1,j)==0) if(f(i,j-1)==0) d=T(i-1,j)-T(i,j-1); r=sqrt(abs(2-d*d)); s=(T(i-1,j)+T(i,j-1)-r)/2; if((s>=T(i-1,j))&(s>=T(i,j-1))) Sol(1)=s; else s=s+r; if((s>=T(i-1,j))&(s>=T(i,j-1))) Sol(1)=s; end end end if(f(i,j-1)~=0) Sol(1)=1+T(i-1,j); end end if(f(i-1,j)~=0) if(f(i,j-1)==0) Sol(1)=1+T(i,j-1); end end %************************************************************************ %sol=solve(i+1,j,i,j-1,f,T);Sol(2)=sol; %************************************************************************ Sol(2)=1e+6; s=0; r=0; d=0; if(f(i+1,j)==0) if(f(i,j-1)==0) d=T(i+1,j)-T(i,j-1); r=sqrt(abs(2-d*d)); s=(T(i+1,j)+T(i,j-1)-r)/2; if((s>=T(i+1,j))&(s>=T(i,j-1)))