function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
%
% [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
%
% converting a 3D volumetric image to surfaces
%

fprintf(1,'extracting surfaces from a volume ...\n');

el=[];
no=[];

if(isstruct(opt) & isfield(opt,'holes')) 
    holes=opt.holes;
else
    holes=[];
end
if(isstruct(opt) & isfield(opt,'regions')) 
    regions=opt.regions;
else
    regions=[];
end
maxlevel=0;

if(~isempty(img))

    img=img(ix,iy,iz);
    dim=size(img);
    newdim=dim+[2 2 2];
    newimg=zeros(newdim);
    newimg(2:end-1,2:end-1,2:end-1)=img;

    if(nargin<8)
        maxlevel=max(newimg(:));
        isovalues=1:maxlevel;
    else
        isovalues=unique(sort(isovalues));
        maxlevel=length(isovalues);
    end

    for i=1:maxlevel
      if(i<maxlevel)
          levelmask=int8(newimg>=isovalues(i) & newimg<isovalues(i+1));
      else
          levelmask=int8(newimg>=isovalues(i));
      end
      [levelno,levelel]=binsurface(levelmask);
      if(~isempty(levelel))
          if(isstruct(opt) & isfield(opt,'autoregion'))
              if(opt.autoregion)
                  seeds=surfseeds(levelno,levelel);
              else
                  seeds=surfinterior(levelno,levelel);
              end
          else
              seeds=surfinterior(levelno,levelel);
          end
          if(~isempty(seeds))
              disp([sprintf('region %d centroid :',i) sprintf('\t%f %f %f\n', seeds')]);
              if(~isempty(regions))
                  regions(end+1:end+size(seeds,1),:)=seeds;
              else
                  regions=seeds;
              end
          end
      end
    end

    for i=1:maxlevel
        fprintf(1,'processing threshold level %d...\n',i);

        if(nargin>=7 & strcmp(method,'simplify'))

          [v0,f0]=binsurface(newimg>=isovalues(i)); % not sure if binsurface works for multi-value arrays
          % with binsurface, I think the following line is not needed anymore
          %  v0(:,[1 2])=v0(:,[2 1]); % isosurface(V,th) assumes x/y transposed
          if(dofix)  [v0,f0]=meshcheckrepair(v0,f0);  end  

          if(isstruct(opt) & length(opt)==maxlevel) keepratio=opt(i).keepratio;
          elseif (isstruct(opt) & length(opt)==1) keepratio=opt.keepratio;
          else keepratio=opt;  end;

          % first, resample the surface mesh with cgal
          fprintf(1,'resampling surface mesh for level %d...\n',i);
          [v0,f0]=meshresample(v0,f0,keepratio);

          % iso2mesh is not stable for meshing small islands,remove them (max 3x3x3 voxels)
          f0=removeisolatedsurf(v0,f0,3);

          if(dofix) [v0,f0]=meshcheckrepair(v0,f0); end

        elseif(nargin<7 | strcmp(method,'cgalsurf') | strcmp(method,'cgalpoly'))
          if(isstruct(opt) & length(opt)==maxlevel) radbound=opt(i).radbound;
          elseif (isstruct(opt) & length(opt)==1) radbound=opt.radbound;
          else radbound=opt;  end;

          maxsurfnode=1000%40000;  % maximum node numbers for each level
          if(isstruct(opt) & length(opt)==maxlevel) 
              if(isfield(opt(i),'maxnode')) maxsurfnode=opt(i).maxnode; end
          elseif (isstruct(opt) & length(opt)==1 )
              if(isfield(opt(1),'maxnode')) 
                 maxsurfnode=opt.maxnode; 
              end
          end

	  distbound=radbound;
          if(isstruct(opt) & length(opt)==maxlevel)
              if(isfield(opt(i),'distbound')) distbound=opt(i).distbound; end
          elseif (isstruct(opt) & length(opt)==1 )
              if(isfield(opt(1),'distbound')) distbound=opt.distbound; end
          end
	  surfside='';
          if(isstruct(opt) & length(opt)==maxlevel)
              if(isfield(opt(i),'side')) surfside=opt(i).side; end
          elseif (isstruct(opt) & length(opt)==1)
              if(isfield(opt(1),'side')) surfside=opt(1).side; end
          end
	  if(~isempty(surfside))
	     newimg0=newimg;
	  end
          if(strcmp(surfside,'upper'))
	      newimg(find(newimg<=isovalues(i)-1e-9))=isovalues(i)-1e-9;
	  elseif(strcmp(surfside,'lower'))
	      newimg(find(newimg>=isovalues(i)+1e-9))=isovalues(i)+1e-9;
	  end
          perturb=1e-4*abs(max(isovalues));
          if(all(newimg>isovalues(i)-perturb)) perturb=-perturb;  end
          [v0,f0]=vol2restrictedtri(newimg,isovalues(i)-perturb,regions(i,:),...
                     sum(newdim.*newdim)*2,30,radbound,distbound,maxsurfnode);

	  if(~isempty(surfside))
            newimg=newimg0;
	    clear newimg0;
	  end
        else
            error('method can only be one of "cgalsurf", "cgalpoly" or "simplify".');
        end

        % if use defines maxsurf=1, take only the largest closed surface
        if(isstruct(opt))
            if( (isfield(opt,'maxsurf') && length(opt)==1 && opt.maxsurf==1) | ...
                    (length(opt)==maxlevel && isfield(opt(i),'maxsurf') && opt(i).maxsurf==1))
                    f0=maxsurf(finddisconnsurf(f0));
            end
        end

        % if a transformation matrix/offset vector supplied, apply them

        if(isstruct(opt) & length(opt)==maxlevel) 
          if(isfield(opt(i),'A') & isfield(opt(i),'B'))
            v0=(opt(i).A*v0'+repmat(opt(i).B(:),1,size(v0,1)))';
          end
        elseif (isstruct(opt) & length(opt)==1) 
          if(isfield(opt,'A') & isfield(opt,'B'))
            v0=(opt.A*v0'+repmat(opt.B(:),1,size(v0,1)))';
          end
        end

        % if user specified holelist and regionlist, append them
        if(isstruct(opt)  & length(opt)==maxlevel)
        if(isfield(opt(i),'hole'))
                holes=[holes;opt(i).hole]
        end
        if(isfield(opt(i),'region'))
                regions=[regions;opt(i).region]
        end
        end

        if(i==0)
        el=[f0 (i+1)*ones(size(f0,1),1)];
        no=v0;
        else
        el=[el;f0+length(no) (i+1)*ones(size(f0,1),1)];
        no=[no;v0];
        end
    end

    %some final fix and scaling
    no(:,1:3)=no(:,1:3)-1; % because we padded the image with a 1 voxel thick null layer in newimg

    no(:,1)=no(:,1)*(max(ix)-min(ix)+1)/dim(1)+(min(ix)-1);
    no(:,2)=no(:,2)*(max(iy)-min(iy)+1)/dim(2)+(min(iy)-1);
    no(:,3)=no(:,3)*(max(iz)-min(iz)+1)/dim(3)+(min(iz)-1);

end  % if not isempty(img)

if(isstruct(opt) & isfield(opt,'surf'))
   for i=1:length(opt.surf)
	opt.surf(i).elem(:,4)=maxlevel+i;
        el=[el;opt.surf(i).elem+length(no)];
        no=[no;opt.surf(i).node];
   end
end

fprintf(1,'surface mesh generation is complete\n');

