demo_surf2vol_ex1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   demo script to convert a closed surface to a binary image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0002 %   demo script to convert a closed surface to a binary image
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 
0005 %% preparation
0006 % user must add the path of iso2mesh to matlab path list
0007 % addpath('../');
0008 
0009 %% load the sample data
0010 load rat_head.mat
0011 
0012 % first, generate a surface from the original image
0013 % similar to demo_shortcuts_ex1.m
0014 
0015 [node,face,regions,holes]=v2s(volimage,0.5,3);
0016 
0017 node=sms(node,face(:,1:3),3,0.5); % apply 3 mesh smoothing
0018 
0019 mdim=ceil(max(node)+1);
0020 dstep=0.25;
0021 zslice=15;
0022 xrange=0:dstep:mdim(1);
0023 yrange=0:dstep:mdim(2);
0024 zrange=0:dstep:mdim(3);
0025 img=surf2vol(node,face(:,1:3),xrange,yrange,zrange);
0026 
0027 imagesc(squeeze(img(:,:,zslice))); % z=10
0028 
0029 hold on
0030 
0031 z0=zslice*dstep;
0032 plane=[min(node(:,1)) min(node(:,2)) z0
0033        min(node(:,1)) max(node(:,2)) z0
0034        max(node(:,1)) min(node(:,2)) z0];
0035 
0036 % run qmeshcut to get the cross-section information at z=mean(node(:,1))
0037 % use the x-coordinates as the nodal values
0038 
0039 [bcutpos,bcutvalue,bcutedges]=qmeshcut(face(:,1:3),node,node(:,1),plane);
0040 [bcutpos,bcutedges]=removedupnodes(bcutpos,bcutedges);
0041 bcutloop=extractloops(bcutedges);
0042 bcutloop(isnan(bcutloop))=[]; % there can be multiple loops, remove the separators
0043 plot(bcutpos(bcutloop,2)*(1/dstep),bcutpos(bcutloop,1)*(1/dstep),'w');
0044 
0045 if(isoctavemesh)
0046   if(~exist('bwfill'))
0047     error('you need to install octave-image toolbox first');
0048   end
0049   img2=zeros(size(img),'uint8');
0050   for i=1:size(img,3)
0051     img2(:,:,i)=bwfill(img(:,:,i),'holes');
0052   end
0053   img2=img2+img;
0054 else
0055   img2=imfill(img,'holes')+img;
0056 end
0057 figure;
0058 imagesc(squeeze(img2(:,:,zslice))); % z=10
0059 hold on;
0060 plot(bcutpos(bcutloop,2)*(1/dstep),bcutpos(bcutloop,1)*(1/dstep),'y--');