vol2surf

 [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)

 converting a 3D volumetric image to surfaces

 author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)

 input:
     img: a volumetric binary image; if img is empty, vol2surf will
          return user defined surfaces via opt.surf if it exists
     ix,iy,iz: subvolume selection indices in x,y,z directions
     opt: function parameters
       if method is 'cgalsurf' or 'cgalpoly':
         opt=a float number>1: max radius of the Delaunay sphere(element size) 
         opt.radbound: same as above, max radius of the Delaunay sphere
         opt.distbound: maximum deviation from the specified isosurfaces
         opt(1,2,...).radbound: same as above, for each levelset
       if method is 'simplify':
         opt=a float number<1: compression rate for surf. simplification
         opt.keepratio=a float less than 1: same as above, same for all surf.
         opt(1,2,..).keepratio: setting compression rate for each levelset
       opt(1,2,..).maxsurf: 1 - only use the largest disjointed surface
                0 - use all surfaces for that levelset
          opt(1,2,..).side: - 'upper': threshold at upper interface
                              'lower': threshold at lower interface
       opt(1,2,..).maxnode: - the maximum number of surface node per levelset
       opt(1,2,..).holes: user specified holes interior pt list
       opt(1,2,..).regions: user specified regions interior pt list
       opt(1,2,..).surf.{node,elem}: add additional surfaces
       opt(1,2,..).{A,B}: linear transformation for each surface
       opt.autoregion: if set to 1, vol2surf will try to determine 
              the interior points for each closed surface automatically
     dofix: 1: perform mesh validation&repair, 0: skip repairing
     method: - if method is 'simplify', iso2mesh will first call
           binsurface to generate a voxel-based surface mesh and then
           use meshresample/meshcheckrepair to create a coarser mesh;
         - if method is 'cgalsurf', iso2mesh will call the surface
           extraction program from CGAL to make surface mesh
         - if method is not specified, 'cgalsurf' is assumed by default
     isovalues: a list of isovalues where the levelset is defined

 output: 
     no:  list of nodes on the resulting suface mesh, 3 columns for x,y,z
     el:  list of trianglular elements on the surface, [n1,n2,n3,region_id]
     regions: list of interior points for all sub-region, [x,y,z]
     holes:   list of interior points for all holes, [x,y,z]

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

0001 function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
0002 %
0003 % [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
0004 %
0005 % converting a 3D volumetric image to surfaces
0006 %
0007 % author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
0008 %
0009 % input:
0010 %     img: a volumetric binary image; if img is empty, vol2surf will
0011 %          return user defined surfaces via opt.surf if it exists
0012 %     ix,iy,iz: subvolume selection indices in x,y,z directions
0013 %     opt: function parameters
0014 %       if method is 'cgalsurf' or 'cgalpoly':
0015 %         opt=a float number>1: max radius of the Delaunay sphere(element size)
0016 %         opt.radbound: same as above, max radius of the Delaunay sphere
0017 %         opt.distbound: maximum deviation from the specified isosurfaces
0018 %         opt(1,2,...).radbound: same as above, for each levelset
0019 %       if method is 'simplify':
0020 %         opt=a float number<1: compression rate for surf. simplification
0021 %         opt.keepratio=a float less than 1: same as above, same for all surf.
0022 %         opt(1,2,..).keepratio: setting compression rate for each levelset
0023 %       opt(1,2,..).maxsurf: 1 - only use the largest disjointed surface
0024 %                0 - use all surfaces for that levelset
0025 %          opt(1,2,..).side: - 'upper': threshold at upper interface
0026 %                              'lower': threshold at lower interface
0027 %       opt(1,2,..).maxnode: - the maximum number of surface node per levelset
0028 %       opt(1,2,..).holes: user specified holes interior pt list
0029 %       opt(1,2,..).regions: user specified regions interior pt list
0030 %       opt(1,2,..).surf.{node,elem}: add additional surfaces
0031 %       opt(1,2,..).{A,B}: linear transformation for each surface
0032 %       opt.autoregion: if set to 1, vol2surf will try to determine
0033 %              the interior points for each closed surface automatically
0034 %     dofix: 1: perform mesh validation&repair, 0: skip repairing
0035 %     method: - if method is 'simplify', iso2mesh will first call
0036 %           binsurface to generate a voxel-based surface mesh and then
0037 %           use meshresample/meshcheckrepair to create a coarser mesh;
0038 %         - if method is 'cgalsurf', iso2mesh will call the surface
0039 %           extraction program from CGAL to make surface mesh
0040 %         - if method is not specified, 'cgalsurf' is assumed by default
0041 %     isovalues: a list of isovalues where the levelset is defined
0042 %
0043 % output:
0044 %     no:  list of nodes on the resulting suface mesh, 3 columns for x,y,z
0045 %     el:  list of trianglular elements on the surface, [n1,n2,n3,region_id]
0046 %     regions: list of interior points for all sub-region, [x,y,z]
0047 %     holes:   list of interior points for all holes, [x,y,z]
0048 %
0049 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0050 %
0051 
0052 fprintf(1,'extracting surfaces from a volume ...\n');
0053 
0054 el=[];
0055 no=[];
0056 
0057 if(isstruct(opt) & isfield(opt,'holes')) 
0058     holes=opt.holes;
0059 else
0060     holes=[];
0061 end
0062 if(isstruct(opt) & isfield(opt,'regions')) 
0063     regions=opt.regions;
0064 else
0065     regions=[];
0066 end
0067 maxlevel=0;
0068 
0069 if(~isempty(img))
0070 
0071     img=img(ix,iy,iz);
0072     dim=size(img);
0073     newdim=dim+[2 2 2];
0074     newimg=zeros(newdim);
0075     newimg(2:end-1,2:end-1,2:end-1)=img;
0076 
0077     if(nargin<8)
0078         maxlevel=max(newimg(:));
0079         isovalues=1:maxlevel;
0080     else
0081         isovalues=unique(sort(isovalues));
0082         maxlevel=length(isovalues);
0083     end
0084 
0085     for i=1:maxlevel
0086       if(i<maxlevel)
0087           levelmask=int8(newimg>=isovalues(i) & newimg<isovalues(i+1));
0088       else
0089           levelmask=int8(newimg>=isovalues(i));
0090       end
0091       [levelno,levelel]=binsurface(levelmask);
0092       if(~isempty(levelel))
0093           if(isstruct(opt) & isfield(opt,'autoregion'))
0094               if(opt.autoregion)
0095                   seeds=surfseeds(levelno,levelel);
0096               else
0097                   seeds=surfinterior(levelno,levelel);
0098               end
0099           else
0100               seeds=surfinterior(levelno,levelel);
0101           end
0102           if(~isempty(seeds))
0103               disp([sprintf('region %d centroid :',i) sprintf('\t%f %f %f\n', seeds')]);
0104               if(~isempty(regions))
0105                   regions(end+1:end+size(seeds,1),:)=seeds;
0106               else
0107                   regions=seeds;
0108               end
0109           end
0110       end
0111     end
0112 
0113     for i=1:maxlevel
0114         fprintf(1,'processing threshold level %d...\n',i);
0115 
0116         if(nargin>=7 & strcmp(method,'simplify'))
0117 
0118           [v0,f0]=binsurface(newimg>=isovalues(i)); % not sure if binsurface works for multi-value arrays
0119           % with binsurface, I think the following line is not needed anymore
0120           %  v0(:,[1 2])=v0(:,[2 1]); % isosurface(V,th) assumes x/y transposed
0121           if(dofix)  [v0,f0]=meshcheckrepair(v0,f0);  end  
0122 
0123           if(isstruct(opt) & length(opt)==maxlevel) keepratio=opt(i).keepratio;
0124           elseif (isstruct(opt) & length(opt)==1) keepratio=opt.keepratio;
0125           else keepratio=opt;  end;
0126 
0127           % first, resample the surface mesh with cgal
0128           fprintf(1,'resampling surface mesh for level %d...\n',i);
0129           [v0,f0]=meshresample(v0,f0,keepratio);
0130 
0131           % iso2mesh is not stable for meshing small islands,remove them (max 3x3x3 voxels)
0132           f0=removeisolatedsurf(v0,f0,3);
0133 
0134           if(dofix) [v0,f0]=meshcheckrepair(v0,f0); end
0135 
0136         elseif(nargin<7 | strcmp(method,'cgalsurf') | strcmp(method,'cgalpoly'))
0137           if(isstruct(opt) & length(opt)==maxlevel) radbound=opt(i).radbound;
0138           elseif (isstruct(opt) & length(opt)==1) radbound=opt.radbound;
0139           else radbound=opt;  end;
0140 
0141           maxsurfnode=40000;  % maximum node numbers for each level
0142           if(isstruct(opt) & length(opt)==maxlevel) 
0143               if(isfield(opt(i),'maxnode')) maxsurfnode=opt(i).maxnode; end
0144           elseif (isstruct(opt) & length(opt)==1 )
0145               if(isfield(opt(1),'maxnode')) 
0146                  maxsurfnode=opt.maxnode; 
0147               end
0148           end
0149 
0150       distbound=radbound;
0151           if(isstruct(opt) & length(opt)==maxlevel)
0152               if(isfield(opt(i),'distbound')) distbound=opt(i).distbound; end
0153           elseif (isstruct(opt) & length(opt)==1 )
0154               if(isfield(opt(1),'distbound')) distbound=opt.distbound; end
0155           end
0156       surfside='';
0157           if(isstruct(opt) & length(opt)==maxlevel)
0158               if(isfield(opt(i),'side')) surfside=opt(i).side; end
0159           elseif (isstruct(opt) & length(opt)==1)
0160               if(isfield(opt(1),'side')) surfside=opt(1).side; end
0161           end
0162       if(~isempty(surfside))
0163          newimg0=newimg;
0164       end
0165           if(strcmp(surfside,'upper'))
0166           newimg(find(newimg<=isovalues(i)-1e-9))=isovalues(i)-1e-9;
0167       elseif(strcmp(surfside,'lower'))
0168           newimg(find(newimg>=isovalues(i)+1e-9))=isovalues(i)+1e-9;
0169       end
0170           perturb=1e-4*abs(max(isovalues));
0171           if(all(newimg>isovalues(i)-perturb)) perturb=-perturb;  end
0172           [v0,f0]=vol2restrictedtri(newimg,isovalues(i)-perturb,regions(i,:),...
0173                      sum(newdim.*newdim)*2,30,radbound,distbound,maxsurfnode);
0174 
0175       if(~isempty(surfside))
0176             newimg=newimg0;
0177         clear newimg0;
0178       end
0179         else
0180             error('method can only be one of "cgalsurf", "cgalpoly" or "simplify".');
0181         end
0182 
0183         % if use defines maxsurf=1, take only the largest closed surface
0184         if(isstruct(opt))
0185             if( (isfield(opt,'maxsurf') && length(opt)==1 && opt.maxsurf==1) | ...
0186                     (length(opt)==maxlevel && isfield(opt(i),'maxsurf') && opt(i).maxsurf==1))
0187                     f0=maxsurf(finddisconnsurf(f0));
0188             end
0189         end
0190 
0191         % if a transformation matrix/offset vector supplied, apply them
0192 
0193         if(isstruct(opt) & length(opt)==maxlevel) 
0194           if(isfield(opt(i),'A') & isfield(opt(i),'B'))
0195             v0=(opt(i).A*v0'+repmat(opt(i).B(:),1,size(v0,1)))';
0196           end
0197         elseif (isstruct(opt) & length(opt)==1) 
0198           if(isfield(opt,'A') & isfield(opt,'B'))
0199             v0=(opt.A*v0'+repmat(opt.B(:),1,size(v0,1)))';
0200           end
0201         end
0202 
0203         % if user specified holelist and regionlist, append them
0204         if(isstruct(opt)  & length(opt)==maxlevel)
0205         if(isfield(opt(i),'hole'))
0206                 holes=[holes;opt(i).hole]
0207         end
0208         if(isfield(opt(i),'region'))
0209                 regions=[regions;opt(i).region]
0210         end
0211         end
0212 
0213         if(i==0)
0214         el=[f0 (i+1)*ones(size(f0,1),1)];
0215         no=v0;
0216         else
0217         el=[el;f0+length(no) (i+1)*ones(size(f0,1),1)];
0218         no=[no;v0];
0219         end
0220     end
0221 
0222     %some final fix and scaling
0223     no(:,1:3)=no(:,1:3)-1; % because we padded the image with a 1 voxel thick null layer in newimg
0224 
0225     no(:,1)=no(:,1)*(max(ix)-min(ix)+1)/dim(1)+(min(ix)-1);
0226     no(:,2)=no(:,2)*(max(iy)-min(iy)+1)/dim(2)+(min(iy)-1);
0227     no(:,3)=no(:,3)*(max(iz)-min(iz)+1)/dim(3)+(min(iz)-1);
0228 
0229 end  % if not isempty(img)
0230 
0231 if(isstruct(opt) & isfield(opt,'surf'))
0232    for i=1:length(opt.surf)
0233     opt.surf(i).elem(:,4)=maxlevel+i;
0234         el=[el;opt.surf(i).elem+length(no)];
0235         no=[no;opt.surf(i).node];
0236    end
0237 end
0238 
0239 fprintf(1,'surface mesh generation is complete\n');
0240