function [newnode,newelem,newelem0]=surfboolean(node,elem,varargin)
%
% [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)
%
% merge two or more triangular meshes and resolve intersecting elements
% 
% author: Qianqian Fang <q.fang at neu.edu>
%
% input:
%      node: node coordinates, dimension (nn,3)
%      elem: triangle surfaces (ne,3)
%      op:  a string of a boolean operator, possible op values include
%           'union' or 'or': the outter surface of the union of the enclosed space
%           'inter' or 'and': the surface of the domain contained by both meshes
%           'diff' or '-': the surface of the domain in mesh 1 excluding that of
%                   mesh 2
%           'all', 'resolve', or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
%                  column of newelem:
%                    1: mesh 1 outside of mesh 2
%                    2: mesh 2 outside of mesh 1
%                    3: mesh 1 inside of mesh 2
%                    4: mesh 2 inside of mesh 1
%                  you can use newelem(find(mod(newelem(:,4),2)==1),:) to
%                  get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:) 
%                  to get mesh 2 cut by mesh 1;
%           'first': combine 1 and 3 from the output of 'all'
%           'second': combine 2 and 4 from the output of 'all'
%           'self': test for self-intersections; only the first mesh is
%                   tested; other inputs are ignored.
%           'remesh': remesh the first input mesh (2nd input mesh is ignored)
%           'decouple': separate two shells and make sure there is no intersection;
%                   the input surfaces must be closed and ordered from outer to inner
%
%           if a dash '-' is preceding a command, including
%              -diff, -union, -isct, -resolve, -first, -second, -xor
%           a 'remesh' operation is applied to the output mesh to remove poorly shaped 
%           triangles
%
% output:
%      newnode: the node coordinates after boolean operations, dimension (nn,3)
%      newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
%      newelem0: when the operator is 'self', return the intersecting
%               element list in terms of the input node list (experimental)
%
% example:
%
%   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
%   [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
%   [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
%   plotmesh(newnode,newface);
%   [newnode,newface]=surfboolean(node1,face1,'-union',node2,face2);
%   figure; plotmesh(newnode,newface);
%   [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
%   figure;plotmesh(newnode,newface,'x>5');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=length(varargin);
newnode=node;
newelem=elem;
if(len>0 && mod(len,3)~=0)
   error('you must give operator, node and element in a triplet form');
end

try
    exename=evalin('caller','ISO2MESH_SURFBOOLEAN');
catch
    try
        exename=evalin('base','ISO2MESH_SURFBOOLEAN');
    catch
        exename='cork';
    end
end

exesuff=fallbackexeext(getexeext,exename);
randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"
if(~isempty(getvarfrom({'caller','base'},'ISO2MESH_RANDSEED')))
    randseed=getvarfrom({'caller','base'},'ISO2MESH_RANDSEED');
end

for i=1:3:len
   op=varargin{i};
   no=varargin{i+1};
   el=varargin{i+2};
   opstr=op;
   if(strcmp(op,'or'))
       opstr='union';
   end
   if(strcmp(op,'xor'))
       opstr='all';
   end
   if(strcmp(op,'and'))
       opstr='isct';
   end
   if(strcmp(op,'-'))
       opstr='diff';
   end
   if(strcmp(op,'self'))
       opstr='solid';
   end
   if(strcmp(opstr,'all'))
       opstr='resolve';
   end
   tempsuff='off';
   deletemeshfile(mwpath(['pre_surfbool*.' tempsuff]));
   deletemeshfile(mwpath('post_surfbool.off'));
   if(strcmp(opstr,'all'))
      deletemeshfile(mwpath('s1out2.off'));
      deletemeshfile(mwpath('s1in2.off'));
      deletemeshfile(mwpath('s2out1.off'));
      deletemeshfile(mwpath('s2in1.off'));
   end
   if(strcmp(op,'decouple'))
       if(exist('node1','var')==0)
          node1=node;
          elem1=elem;
          newnode(:,4)=1;
          newelem(:,4)=1;
       end
       opstr=' --decouple-inin 1 --shells 2'; %-q
       saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
       if(size(no,2)~=3)
           opstr=['-q --shells ' num2str(no)];
           cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),opstr);
       else
           saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
           cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
       end
   elseif(strcmp(op,'decoupleout'))
       if(exist('node1','var')==0)
          node1=node;
          elem1=elem;
          newnode(:,4)=1;
          newelem(:,4)=1;
       end
       opstr=[' --decouple-outout 1 --shells 2']; %-q
       saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
       if(size(no,2)~=3)
           opstr=['-q --shells ' num2str(no)];
           cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),opstr);
       else
           saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
           cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
       end
   elseif(strcmp(op,'separate'))
       if(exist('node1','var')==0)
          node1=node;
          elem1=elem;
          newnode(:,4)=1;
          newelem(:,4)=1;
       end
       opstr=' --shells 2'; %-q
       saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
       if(size(no,2)~=3)
           opstr=['-q --shells ' num2str(no)];
           cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),opstr);
       else
           saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
           cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
       end       
   elseif(strcmp(op,'remesh'))
       saveoff(newnode(:,1:3),newelem(:,1:3),mwpath(['pre_surfbool1.' tempsuff]));
       cmd=sprintf('cd "%s" && "%s%s" %s%s "%s" "%s" -%d',mwpath,mcpath(exename),exesuff,'-',...
            opstr,mwpath(['pre_surfbool1.' tempsuff]),mwpath('post_surfbool.off'),randseed);
   else
       saveoff(newnode(:,1:3),newelem(:,1:3),mwpath(['pre_surfbool1.' tempsuff]));
       saveoff(no(:,1:3),el(:,1:3),mwpath(['pre_surfbool2.' tempsuff]));
       cmd=sprintf('cd "%s" && "%s%s" %s%s "%s" "%s" "%s" -%d',mwpath,mcpath(exename),exesuff,'-',...
            opstr,mwpath(['pre_surfbool1.' tempsuff]),mwpath(['pre_surfbool2.' tempsuff]),mwpath('post_surfbool.off'),randseed);
   end
   [status, outstr]=system(cmd);
   if(status~=0 && strcmp(op,'self')==0)
       error('surface boolean command failed:\n%s\nERROR: %s\n',cmd,outstr);
   end
   if(strcmp(op,'self'))
       if(isempty(strfind(outstr,'NOT SOLID')))
           fprintf(1,'no self-intersection was found!\n');
           newnode=[];
           newelem=[];
           newelem0=[];
       else
           fprintf(1,'input mesh is self-intersecting\n');
           newnode=1;
           newelem=[];
           newelem0=1;
       end
       return
   end
   if(strcmp(opstr,'all'))
      % tag the 4 piceses of meshes, this tag do not propagate to the next boolean operation
      [nnode, nelem]=readoff(mwpath('s1out2.off'));
      newelem=[nelem ones(size(nelem,1),1)];
      newnode=[nnode ones(size(nnode,1),1)];

      [nnode, nelem]=readoff(mwpath('s1in2.off'));
      newelem=[newelem; nelem+size(newnode,1) 3*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 3*ones(size(nnode,1),1)];

      [nnode, nelem]=readoff(mwpath('s2out1.off'));
      newelem=[newelem; nelem+size(newnode,1) 2*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 2*ones(size(nnode,1),1)];

      [nnode, nelem]=readoff(mwpath('s2in1.off'));
      newelem=[newelem; nelem+size(newnode,1) 4*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 4*ones(size(nnode,1),1)];
   elseif(strcmp(op,'decouple'))
      [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off')); %[node1,elem1]
      %newelem=[newelem;elem1+size(newnode,1) (i+1)*ones(size(elem1,1),1)];
      %newnode=[newnode;node1 (i+1)*ones(size(node1,1),1)];
   elseif(strcmp(op,'separate'))
       [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off'));
   elseif(strcmp(op,'decoupleout'))
      [newnode,newelem]=readoff(mwpath('pre_decouple1_fixed.off')); %[node1,elem1]
   else
      [newnode,newelem]=readoff(mwpath('post_surfbool.off'));
      if(strcmp(op,'self'))
          fprintf(1,'a total of %d self-intersecting elements were found\n',size(newelem,1));
          if(nargout>=3)
              [found,newelem0]=ismember(newnode,node,'rows');
              if(~all(found))
                  error('self intersecting elements contain new nodes');
              end
              newelem0=newelem0(newelem);
          end
          return;
      end
   end
end
