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;
没有合适的资源?快使用搜索试试~ 我知道了~
matlab-医学切片图片三维重建得到三维模型matlab仿真-源码
共771个文件
jpg:760个
m:6个
db:5个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
5星 · 超过95%的资源 1 下载量 169 浏览量
2021-09-15
00:00:14
上传
评论 3
收藏 12.3MB RAR 举报
温馨提示
matlab_医学切片图片三维重建得到三维模型matlab仿真_源码
资源推荐
资源详情
资源评论
收起资源包目录
matlab-医学切片图片三维重建得到三维模型matlab仿真-源码 (771个子文件)
Thumbs.db 315KB
Thumbs.db 315KB
Thumbs.db 313KB
Thumbs.db 313KB
Thumbs.db 313KB
CT0304.jpg 21KB
CT0684.jpg 21KB
CT0076.jpg 21KB
CT0152.jpg 20KB
CT0608.jpg 20KB
CT0760.jpg 20KB
CT0228.jpg 20KB
CT0380.jpg 20KB
CT0151.jpg 20KB
CT0075.jpg 20KB
CT0607.jpg 20KB
CT0683.jpg 20KB
CT0227.jpg 20KB
CT0759.jpg 20KB
CT0531.jpg 20KB
CT0074.jpg 20KB
CT0454.jpg 20KB
CT0682.jpg 20KB
CT0758.jpg 20KB
CT0303.jpg 20KB
CT0606.jpg 20KB
CT0456.jpg 19KB
CT0455.jpg 19KB
CT0530.jpg 19KB
CT0532.jpg 19KB
CT0378.jpg 19KB
CT0377.jpg 19KB
CT0681.jpg 19KB
CT0453.jpg 19KB
CT0605.jpg 19KB
CT0301.jpg 19KB
CT0150.jpg 19KB
CT0757.jpg 19KB
CT0529.jpg 19KB
CT0379.jpg 19KB
CT0226.jpg 19KB
CT0225.jpg 19KB
CT0302.jpg 19KB
CT0224.jpg 19KB
CT0376.jpg 19KB
CT0073.jpg 19KB
CT0604.jpg 19KB
CT0300.jpg 18KB
CT0452.jpg 18KB
CT0528.jpg 18KB
CT0680.jpg 18KB
CT0148.jpg 18KB
CT0149.jpg 18KB
CT0375.jpg 18KB
CT0223.jpg 18KB
CT0147.jpg 18KB
CT0299.jpg 18KB
CT0451.jpg 18KB
CT0071.jpg 18KB
CT0527.jpg 18KB
CT0756.jpg 18KB
CT0072.jpg 18KB
CT0374.jpg 18KB
CT0754.jpg 18KB
CT0146.jpg 18KB
CT0755.jpg 18KB
CT0070.jpg 18KB
CT0298.jpg 18KB
CT0222.jpg 17KB
CT0651.jpg 17KB
CT0603.jpg 17KB
CT0575.jpg 17KB
CT0043.jpg 17KB
CT0450.jpg 17KB
CT0727.jpg 17KB
CT0679.jpg 17KB
CT0119.jpg 17KB
CT0678.jpg 17KB
CT0195.jpg 17KB
CT0271.jpg 17KB
CT0344.jpg 17KB
CT0274.jpg 17KB
CT0420.jpg 17KB
CT0499.jpg 17KB
CT0044.jpg 17KB
CT0574.jpg 17KB
CT0728.jpg 17KB
CT0350.jpg 17KB
CT0347.jpg 17KB
CT0045.jpg 17KB
CT0650.jpg 17KB
CT0652.jpg 17KB
CT0426.jpg 17KB
CT0578.jpg 17KB
CT0602.jpg 17KB
CT0275.jpg 17KB
CT0502.jpg 17KB
CT0046.jpg 17KB
CT0198.jpg 17KB
CT0345.jpg 17KB
共 771 条
- 1
- 2
- 3
- 4
- 5
- 6
- 8
资源评论
- 秃头小嘿嘿2023-04-25资源有一定的参考价值,与资源描述一致,很实用,能够借鉴的部分挺多的,值得下载。
mYlEaVeiSmVp
- 粉丝: 1908
- 资源: 19万+
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功