www.pudn.com > imageregistration.rar > peizhunchengxu.m, change:2015-04-26,size:3121b


clear; 
clc;  
I1=dicomread('E:\考研准备课程\IM154.IM154');  
 I1=double(I1);  
I2=dicomread('E:\考研准备课程\IM70.IM70');  
I2=double(I2); 
 
%截取部分图像 
figure(1),  
subplot(2,3,1),imshow(I1,[]);  
title('浮动图像');  
image1=I1(106:106+300-1,... 
    42:42+300-1); 
%e = imrect(gca,[0 0 300 300]); 
subplot(2,3,2),imshow(I2,[]);  
title('参考图像');  
%e = imrect(gca,[0 0 300 300]); 
image2=I2(106:106+300-1,... 
    42:42+300-1); 
subplot(2,3,3),imshow(image1,[]); 
title('局部浮动图像'); 
subplot(2,3,4),imshow(image2,[]); 
title('局部参考图像'); 
%初始化K值为0      
k(1:4) = 0;      
     
%金字塔循环     
for pyrIter=2:-1:0,      
  step=2.^pyrIter;      
     
%变换参数值二倍      
  k = [k(1) k(2) 2.*k(3) 2.*k(4)];      
     
  im1 = imresize(image1,1 /step,'nearest');      
  im2 = imresize(image2, 1/step,'nearest');      
     
%figure,subplot(2,2,1),imshow(im1,[]);  
%subplot(2,2,2),imshow(im2,[]);  
 for l=1:3      
%运动补偿帧2 l是L的小写     
    warped = warp(im2,k,pyrIter);      
     
%求图像梯度      
   [Ix Iy It] = grad(im1,warped);      
     
%找剩余的K值      
  resid_k = Parametric(Ix,Iy,It);      
     
%更新K值      
  k = k + resid_k;      
  end      
end  
Xn=warped;  
Xnp1=image2; 
x=Xn;y=Xnp1; 
[m,n]=size(x); 
%m=size(x,1);  
%n=size(x,2);  
h=m;w=n;  
x0=zeros(1,n);y0=zeros(1,m); z0=zeros(1,1);  
x=[x;x0];x=x';x=[x;y0,z0];x=x';  
y=[y;x0];y=y';y=[y;y0,z0];y=y';  
% [h,w]=size(Xn)  
% hm5=h-5; wm5=w-5;  
z=zeros(h,w); v1=z; v2=z;  
evaluate=z;  
%初始化  
dsx2=v1; dsx1=v1; dst=v1;  
alpha2=625;  
imax=10;  
  
%有限差分的梯度估计   
for i=1:m   
for j=1:n   
dsx2(i,j)=(x(i+1,j)-x(i,j)+x(i+1,j+1)-x(i,j+1)+y(i+1,j)-y(i,j)+y(i+1,j+1)-y(i,j+1))./4;   
dsx1(i,j)=(x(i,j+1)-x(i,j)+x(i+1,j+1)-x(i+1,j)+y(i,j+1)-y(i,j)+y(i+1,j+1)-y(i+1,j))./4;   
dst(i,j)=(y(i,j)-x(i,j)+y(i+1,j)-x(i+1,j)+y(i,j+1)-x(i,j+1)+y(i+1,j+1)-x(i+1,j+1))./4;   
end   
end  
  
for i=1:imax  
delta=(dsx1.*v1+dsx2.*v2+dst)./(alpha2+dsx1.^2+dsx2.^2);  
v1=v1-dsx1.*delta;  
v2=v2-dsx2.*delta;  
end;  
u=z; u(1:h,1:w)=v1(1:h,1:w);  
v=z; v(1:h,1:w)=v2(1:h,1:w);  
  
x(size(x,1),:)=[];x(:,size(x,2))=[];  
y(size(y,1),:)=[];y(:,size(y,2))=[];  
  
[s,t]=meshgrid(1:size(y,2),1:size(y,1));  
l1=s+u;  
l2=t+v;  
evaluate=interp2(s,t,y,l1,l2,'bilinear');  
evaluate(isnan(evaluate)) = 0;      
  
%R3=mean(mean(abs(evaluate-y)))  
%V3=sqrt(mean(mean((evaluate-y-R3)^2)))  
  
xskip=round(h/32);%计算光流  
[hs,ws]=size(u(1:xskip:h,1:xskip:w));  
us=zeros(hs,ws); vs=us;  
  
N=xskip^2;  
for i=1:hs-1  
for j=1:ws-1  
hk=i*xskip-xskip+1;  
hl=i*xskip;  
wk=j*xskip-xskip+1;  
wl=j*xskip;  
us(i,j)=sum(sum(u(hk:hl,wk:wl)))/N;  
vs(i,j)=sum(sum(v(hk:hl,wk:wl)))/N;  
end;  
end;  
%输出光流场  
figure(5);  
quiver(us,vs);  
colormap('default');  
axis ij;  
axis tight;  
axis equal; 
title('光流场'); 
%画恢复后的图像   
figure(1);   
subplot(2,3,5)   
imshow(evaluate,[]);   
title('光流变换后的图象');  
%剪影图像 
figure(6),subplot(1,2,1),imshow((image1-image2),[]); 
title('原第一、二帧图像之差'); 
figure(6),  
subplot(1,2,2);  
imshow((evaluate-y),[]);  
title('光流变换后的图像与第二帧差')