binsurface

 [node,elem]=binsurface(img,nface)

 fast isosurface extraction from 3D binary images

 author: Qianqian Fang, <q.fang at neu.edu>

 input: 
   img:  a 3D binary image
   nface: nface=3 or ignored - for triangular faces, 
          nface=4 - square faces
          nface=0 - return a boundary mask image via node

 output:
   elem: integer array with dimensions of NE x nface, each row represents
         a surface mesh face element 
   node: node coordinates, 3 columns for x, y and z respectively

 the outputs of this subroutine can be easily plotted using 
     patch('Vertices',node,'faces',elem,'FaceVertexCData',node(:,3),
           'FaceColor','interp');
 if the surface mesh has triangular faces, one can plot it with
     trisurf(elem,node(:,1),node(:,2),node(:,3))

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [node,elem]=binsurface(img,nface)
0002 %
0003 % [node,elem]=binsurface(img,nface)
0004 %
0005 % fast isosurface extraction from 3D binary images
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 %
0009 % input:
0010 %   img:  a 3D binary image
0011 %   nface: nface=3 or ignored - for triangular faces,
0012 %          nface=4 - square faces
0013 %          nface=0 - return a boundary mask image via node
0014 %
0015 % output:
0016 %   elem: integer array with dimensions of NE x nface, each row represents
0017 %         a surface mesh face element
0018 %   node: node coordinates, 3 columns for x, y and z respectively
0019 %
0020 % the outputs of this subroutine can be easily plotted using
0021 %     patch('Vertices',node,'faces',elem,'FaceVertexCData',node(:,3),
0022 %           'FaceColor','interp');
0023 % if the surface mesh has triangular faces, one can plot it with
0024 %     trisurf(elem,node(:,1),node(:,2),node(:,3))
0025 %
0026 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0027 %
0028 
0029 dim=size(img);
0030 if(length(dim)<3)
0031     dim(3)=1;
0032 end
0033 newdim=dim+[1 1 1];
0034 
0035 % find the jumps (0->1 or 1->0) for all directions
0036 d1=diff(img,1,1);
0037 d2=diff(img,1,2);
0038 d3=diff(img,1,3);
0039 [ix,iy]=find(d1==1|d1==-1);
0040 [jx,jy]=find(d2==1|d2==-1);
0041 [kx,ky]=find(d3==1|d3==-1);
0042 
0043 % compensate the dim. reduction due to diff, and
0044 % wrap the indices in a bigger array (newdim)
0045 ix=ix+1;
0046 [iy,iz]=ind2sub(dim(2:3),iy);
0047 iy=sub2ind(newdim(2:3),iy,iz);
0048 
0049 [jy,jz]=ind2sub([dim(2)-1,dim(3)],jy);
0050 jy=jy+1;
0051 jy=sub2ind(newdim(2:3),jy,jz);
0052 
0053 [ky,kz]=ind2sub([dim(2),dim(3)-1],ky);
0054 kz=kz+1;
0055 ky=sub2ind(newdim(2:3),ky,kz);
0056 
0057 id1=sub2ind(newdim,ix,iy);
0058 id2=sub2ind(newdim,jx,jy);
0059 id3=sub2ind(newdim,kx,ky);
0060 
0061 if(nargin==2 && nface==0)
0062     elem=[id1 id2 id3];
0063     node=zeros(newdim);
0064     node(elem)=1;
0065     node=node(2:end-1,2:end-1,2:end-1)-1;
0066     return
0067 end
0068 
0069 % populate all the triangles
0070 xy=newdim(1)*newdim(2);
0071 
0072 if(nargin==1 || nface==3)  % create triangular elements
0073     elem=[id1 id1+newdim(1) id1+newdim(1)+xy; id1 id1+newdim(1)+xy id1+xy];
0074     elem=[elem; id2 id2+1 id2+1+xy; id2 id2+1+xy id2+xy];
0075     elem=[elem; id3 id3+1 id3+1+newdim(1); id3 id3+1+newdim(1) id3+newdim(1)];
0076 else                       % create box elements
0077     elem=[id1 id1+newdim(1) id1+newdim(1)+xy id1+xy];
0078     elem=[elem; id2 id2+1 id2+1+xy id2+xy];
0079     elem=[elem; id3 id3+1 id3+1+newdim(1) id3+newdim(1)];
0080 end
0081 % compress the node indices
0082 nodemap=zeros(max(elem(:)),1);
0083 nodemap(elem(:))=1;
0084 id=find(nodemap);
0085 
0086 nodemap=0;
0087 nodemap(id)=1:length(id);
0088 elem=nodemap(elem);
0089 
0090 % create the coordiniates
0091 [xi,yi,zi]=ind2sub(newdim,id);
0092 
0093 % assuming the origin [0 0 0] is located at the lower-bottom corner of the image
0094 node=[xi(:),yi(:),zi(:)]-1;