matlab-(含教程)基于matlab的医学图片序列的三维光流场提取仿真

function [fieldx,fieldy,fieldz]=func_cal_demons(img1,img2,initfieldx, initfieldy,initfieldz,iter) global clims; img1_ori=img1; img2_ori=img2; imgsize=size(img1); %gaussian filter gvariance=2; width=ceil(3*gvariance); if floor(width/2)== width/2 width=width+1; end whalf=floor(width/2); gsum=0; h=zeros(width,width,width); for i=-whalf:whalf for j=-whalf:whalf for k=-whalf:whalf h(i+whalf+1,j+whalf+1,k+whalf+1)=exp(-(i*i+j*j+k*k)/2/gvariance/gvariance);%/gvariance); gsum=gsum+h(i+whalf+1,j+whalf+1,k+whalf+1); end end end if gsum~=0 h=h/gsum; end %smooth image with gaussian filter size(img1) img1=convn(img1,h); img1=img1(3:imgsize(1)+2,3:imgsize(2)+2,3:imgsize(3)+2); size(img1) img2=convn(img2,h); img2=img2(3:imgsize(1)+2,3:imgsize(2)+2,3:imgsize(3)+2); hv=h; %if initial field is empty, set to zeros tag=0; if isequal(initfieldx,[]) u=zeros(imgsize); else u=zeros(size(initfieldx)); u=initfieldx; tag=tag+1; end if isequal(initfieldy,[]) v=zeros(imgsize); else v=zeros(size(initfieldy)); v=initfieldy; tag=tag+1; end if isequal(initfieldz,[]) w=zeros(imgsize); else w=zeros(size(initfieldz)); w=initfieldz; tag=tag+1; end img=img2; %calculate gradient for reference image for i=1:imgsize(3) gx1(:,:,i)=circshift(img1(:,:,i),[0,1])-circshift(img1(:,:,i),[0,-1]); gx1(:,:,i)=gx1(:,:,i)/2; gy1(:,:,i)=circshift(img1(:,:,i),[1,0])-circshift(img1(:,:,i),[-1,0]); gy1(:,:,i)=gy1(:,:,i)/2; if i>1 & i<imgsize(3) gz1(:,:,i)=-(img1(:,:,i+1)-img1(:,:,i-1))/2; else gz1(:,:,i)=zeros(imgsize(1:2)); end end grad12=gx1.*gx1+gy1.*gy1+gz1.*gz1; % if initial fields not empty, use them to find initial deformed image 2 if ~(tag==0) for k=1:imgsize(3) for i=1:imgsize(1) for j=1:imgsize(2) ni=i+v(i,j,k); nj=j+u(i,j,k); nk=k+w(i,j,k); if ni<1 | ni>imgsize(1) | nj<1 | nj>imgsize(2) | nk<1 | nk>imgsize(3) img(i,j,k)=0; else img(i,j,k)=img2(round(ni),round(nj),round(nk)); end end end end end ndisp=8; set(gcf,'doublebuffer','on'); subplot(2,3,3); imshow(img(:,:,ndisp),clims); subplot(2,3,1); imshow(img1(:,:,ndisp),clims); subplot(2,3,2); imshow(img2(:,:,ndisp),clims); subplot(2,3,4); imshow(grad12(:,:,ndisp),[]); subplot(2,3,5); % imshow(grad22(:,:,ndisp),[]); subplot(2,3,6); imshow((img2(:,:,ndisp)-img1(:,:,ndisp)),[]); defimg=img2_ori; currdir=pwd; currdir=[currdir '\results_deform\']; for n=1:iter n stepfactor=0.4; di=img-img1; di2=di.*di; %new gradient for image 2 for i=1:imgsize(3) gx2(:,:,i)=circshift(img(:,:,i),[0,1])-circshift(img(:,:,i),[0,-1]); gx2(:,:,i)=gx2(:,:,i)/2; gy2(:,:,i)=circshift(img(:,:,i),[1,0])-circshift(img(:,:,i),[-1,0]); gy2(:,:,i)=gy2(:,:,i)/2; if i>1 & i<imgsize(3) gz2(:,:,i)=-(img(:,:,i+1)-img(:,:,i-1))/2; else gz2(:,:,i)=zeros(imgsize(1:2)); end end grad22=gx2.*gx2+gy2.*gy2+gz2.*gz2; for k=1:imgsize(3) %for each voxel for i=1:imgsize(1) for j=1:imgsize(2) du(i,j,k)=0; dv(i,j,k)=0; dw(i,j,k)=0; if grad12(i,j,k)>10 %passive force tmp=grad12(i,j,k)+stepfactor*di2(i,j,k); if ~isequal(tmp,0) du(i,j,k)=di(i,j,k)*gx1(i,j,k)/tmp; dv(i,j,k)=di(i,j,k)*gy1(i,j,k)/tmp; dw(i,j,k)=di(i,j,k)*gz1(i,j,k)/tmp; end end if grad22(i,j,k)>10 %active force tmp2=grad22(i,j,k)+stepfactor*di2(i,j,k); if ~isequal(tmp2,0) du(i,j,k)=du(i,j,k)+di(i,j,k)*gx2(i,j,k)/tmp2; dv(i,j,k)=dv(i,j,k)+di(i,j,k)*gy2(i,j,k)/tmp2; dw(i,j,k)=dw(i,j,k)+di(i,j,k)*gz2(i,j,k)/tmp2; end end %find average force du(i,j,k)=du(i,j,k)/2; dv(i,j,k)=dv(i,j,k)/2; dw(i,j,k)=dw(i,j,k)/2; end end end u=u+du; v=v+dv; w=w+dw; %smooth field(force) u=convn(u,hv); u=u(3:imgsize(1)+2,3:imgsize(2)+2,3:imgsize(3)+2); v=convn(v,hv); v=v(3:imgsize(1)+2,3:imgsize(2)+2,3:imgsize(3)+2); w=convn(w,hv); w=w(3:imgsize(1)+2,3:imgsize(2)+2,3:imgsize(3)+2); %deform image2 using updated force for k=1:imgsize(3) for i=1:imgsize(1) for j=1:imgsize(2) ni=i+v(i,j,k); nj=j+u(i,j,k); nk=k+w(i,j,k); if ni<1 | ni>imgsize(1) | nj<1 | nj>imgsize(2) | nk<1 | nk>imgsize(3) img(i,j,k)=0; else ny1=floor(ni); nx1=floor(nj); nz1=floor(nk); ny2=round(ny1+1); nx2=round(nx1+1); nz2=round(nz1+1); if isequal(ny1,ni) ny2=ny1; end if isequal(nx1,nj) nx2=nx1; end if isequal(nz1,nk) nz2=nz1; end dy1=ni-ny1; dx1=nj-nx1; dz1=nk-nz1; dy2=1-dy1; dx2=1-dx1; dz2=1-dz1; tmp1=(img2(ny1,nx1,nz1)*dx2*dy2+img2(ny1,nx2,nz1)*dx1*dy2+img2(ny2,nx1,nz1)*dx2*dy1+img2(ny2,nx2,nz1)*dx1*dy1); tmp2=(img2(ny1,nx1,nz2)*dx2*dy2+img2(ny1,nx2,nz2)*dx1*dy2+img2(ny2,nx1,nz2)*dx2*dy1+img2(ny2,nx2,nz2)*dx1*dy1); img(i,j,k)=tmp1*dz2+tmp2*dz1; tmp1=round(img2_ori(ny1,nx1,nz1)*dx2*dy2+img2_ori(ny1,nx2,nz1)*dx1*dy2+img2_ori(ny2,nx1,nz1)*dx2*dy1+img2_ori(ny2,nx2,nz1)*dx1*dy1); tmp2=round(img2_ori(ny1,nx1,nz2)*dx2*dy2+img2_ori(ny1,nx2,nz2)*dx1*dy2+img2_ori(ny2,nx1,nz2)*dx2*dy1+img2_ori(ny2,nx2,nz2)*dx1*dy1); defimg(i,j,k)=tmp1*dz2+tmp2*dz1; end end end end subplot(2,3,3); imshow(img(:,:,ndisp),clims); subplot(2,3,5); g=squeeze(grad22(:,32,:)); imshow(g,[]); subplot(2,3,6); imshow((img(:,:,ndisp)-img1(:,:,ndisp)),[-500,500]); pause(0.5); clear fname; end %for function return purpose fieldx=u; fieldy=v; fieldz=w;