关于数字稳像的matlab程序

tic N=16; r=32; m=7; mov=aviread('foreman'); D=aviinfo('foreman'); for i=1:D.NumFrames A(:,:,:,i)=mov(i).cdata; end B=mov.colormap; for j=1:D.NumFrames F(j)=im2frame(A(:,:,:,j)); end %movie(F); BR2(:,:,:)=A(1:8*N,1:8*N,:,1); Sf=size(A); e2=3; f2=3; % Block motion estimation g=zeros(Sf(4),2); g1=zeros(Sf(4),2); g2=zeros(Sf(4),2); for i=1:D.NumFrames-1 S=size(BR2); h=S(1)/N; w=S(2)/N; b=0; for q=1:h for k=1:w refblk=BR2((q-1)*N+1:q*N,(k-1)*N+1:k*N,:); rt=(2*(f2+k-6)-1)*N/2+r-1; if rt>Sf(2)-N/2 rt=Sf(2)-N/2; end lt=(2*(f2+k+2)-1)*N/2-r; if lt<N/2 lt=N/2; end tp=(2*(e2+q-6)-1)*N/2-r; if tp<N/2 tp=N/2; end bt=(2*(e2+q+2)-1)*N/2+r-1; if bt>Sf(1)-N/2 bt=Sf(1)-N/2; end e=tp; f=lt; for l=tp:bt for n=lt:rt if i+1<=Sf(4) slidblk=A(l-N/2+1:l+N/2,n-N/2+1:n+N/2,:,i+1); else slidblk=A(l-N/2+1:l+N/2,n-N/2+1:n+N/2,:,i); end % calculating SAD sad(l,n)=0; for o=1:N for p=1:N a=slidblk(o,p)-refblk(o,p); sad(l,n)=sad(l,n)+abs(a); end end if sad(l,n)<sad(e,f) e=l; f=n; end end end b=b+1; mv(b,1)=e-(2*(e2+q-3)-1)*N/2; mv(b,2)=f-(2*(f2+k-3)-1)*N/2; end end % GLOBAL MOTION DETECTION edges=cell(2,1); edges={-r:r-1,-r:r-1}; H=hist3(mv,'Edges',edges); px(i)=1; py(i)=1; for k=1:2*r for l=1:2*r H1(k,l)=0; for n=-2:2 for o=-2:2 if (1<=(k+n))&&((k+n)<=2*r)&&(1<=(l+o))&&((l+o)<=2*r) H1(k,l)=H1(k,l)+H(k+n,l+o); else H1(k,l)=H1(k,l); end end end if H1(k,l)>H1(px(i),py(i)) px(i)=k; py(i)=l; end end end for n=-1:1 for o=-1:1 if (1<=(px(i)+n))&&((px(i)+n)<=2*r)&&(1<=(py(i)+o))&&((py(i)+o)<=2*r) g(i,1)=g(i,1)+H(px(i)+n,py(i)+o)*(px(i)-r-1-n); g(i,2)=g(i,2)+H(px(i)+n,py(i)+o)*(py(i)-r-1-o); else g(i,1)=g(i,1); g(i,2)=g(i,2); end end end g(i,1)= g(i,1)/H1(px(i),py(i)); g(i,2)= g(i,2)/H1(px(i),py(i)); % GLOBAL MOTION SMOOTHING z = gausswin(m); z=z/norm(z); for n=1:m if (i-n+1)>=1 g1(i,1)=g1(i,1)+z(n)*g(i-n+1,1); g1(i,2)=g1(i,2)+z(n)*g(i-n+1,2); else g1(i,1)=g1(i,1); g1(i,2)=g1(i,2); end end % global motion compensation d(i,1)=convergent(g1(i,1)-g(i,1)); d(i,2)=convergent(g1(i,2)-g(i,2)); for j=1:S(1) for k=1:S(2) if (1<=(j+d(i,1)))&&((j+d(i,1))<=S(1))&&(1<=(k+d(i,2)))&&((k+d(i,2))<=S(2)) BR2(j,k,:)=BR2(j+d(i,1),k+d(i,2),:); else BR2(j,k,:)=0; end end end % ENCODER % BLOCK MOTION ESTIMATION FOR ENTIRE FRAME h1=Sf(1)/N; w1=Sf(2)/N; b1=0; for q1=1:h1 for k1=1:w1 refblk1=A((q1-1)*N+1:q1*N,(k1-1)*N+1:k1*N,:,i); rt1=(2*k1-1)*N/2+r-1; if rt1>Sf(2)-N/2 rt1=Sf(2)-N/2; end lt1=(2*k1-1)*N/2-r; if lt1<N/2 lt1=N/2; end tp1=(2*q1-1)*N/2-r; if tp1<N/2 tp1=N/2; end bt1=(2*q1-1)*N/2+r-1; if bt1>Sf(1)-N/2 bt1=Sf(1)-N/2; end e1=tp1; f1=lt1; for l1=tp1:bt1 for n1=lt1:rt1 if i+1<=Sf(4) slidblk1=A(l1-N/2+1:l1+N/2,n1-N/2+1:n1+N/2,:,i+1); else slidblk1=A(l1-N/2+1:l1+N/2,n1-N/2+1:n1+N/2,:,i); end % calculating SAD sad1(l1,n1)=0; for o1=1:N for p1=1:N a1=slidblk1(o1,p1)-refblk1(o1,p1); sad1(l1,n1)=sad1(l1,n1)+abs(a1); end end if sad1(l1,n1)<sad1(e1,f1) e1=l1; f1=n1; end end end b1=b1+1; mv1(b1,1)=e1-(2*q1-1)*N/2; mv1(b1,2)=f1-(2*k1-1)*N/2; end end % GLOBAL MOTION DETECTION edges=cell(2,1); edges={-r:r-1,-r:r-1}; H2=hist3(mv1,'Edges',edges); px1(i)=1; py1(i)=1; for k1=1:2*r for l1=1:2*r H3(k1,l1)=0; for n1=-2:2 for o1=-2:2 if (1<=(k1+n1))&&((k1+n1)<=2*r)&&(1<=(l1+o1))&&((l1+o1)<=2*r) H3(k1,l1)=H3(k1,l1)+H2(k1+n1,l1+o1); else H3(k1,l1)=H3(k1,l1); end end end if H3(k1,l1)>H3(px1(i),py1(i)) px1(i)=k1; py1(i)=l1; end end end for n1=-2:2 for o1=-2:2 if (1<=(px1(i)+n1))&&((px1(i)+n1)<=2*r)&&(1<=(py1(i)+o1))&&((py1(i)+o1)<=2*r) g2(i,1)=g2(i,1)+H2(px1(i)+n1,py1(i)+o1)*(px1(i)-r/4-1-n1); g2(i,2)=g2(i,2)+H2(px1(i)+n1,py1(i)+o1)*(py1(i)-r/4-1-o1); else g2(i,1)=g2(i,1); g2(i,2)=g2(i,2); end end end g2(i,1)= g2(i,1)/H3(px1(i),py1(i)); g2(i,2)= g2(i,2)/H3(px1(i),py1(i)); % BR Estimation n2 = wgn(h1,w1,2); e2=3; f2=3; for u=3:h1-2 for v=3:w1-2 br{u,v}=0; h3=mv1(u*h1+v,:)+n2(u,v); for a2=-2:2 for b2=-2:2 br{u,v}=br{u,v}+abs(mv1((u+a2)*h1+v+b2,:)-h3); end end if br{u,v}<br{e2,f2} e2=u; f2=v; end end end if i+1<Sf(4) BR2(:,:,:)=A((e2-5)*N+1:(e2+3)*N,(f2-5)*N+1:(f2+3)*N,:,i+1); else BR2(:,:,:)=A((e2-5)*N+1:(e2+3)*N,(f2-5)*N+1:(f2+3)*N,:,i); end end % block motion compensation for j=1:h1 for k=1:w1 d1=mv1(j*h1+k,1); d2=mv1(j*h1+k,2); for u1=1:N for v1=1:N if (1<=(j*N+u1+d1))&&((j*N+u1+d1)<=S(1))&&(1<=(k*N+v1+d2))&&((k*N+v1+d2)<=S(2)) A(j*N+u1,k*N+v1,:,i)=A(j*N+u1+d1,k*N+v1+d2,:,i); else A(j*N+u1,k*N+v1,:,i)=0; end end end end end plot(1:Sf(4),g3(:,2)); toc