function [bmask,uniqidx,nc,nn,totalarea,iso]=mcxsvmc(vol, varargin)
%
% Format:
%    newvol=mcxsvmc(vol)
%       or
%    [bmask,uniqidx,nn,nc,totalarea,iso]=mcxsvmc(vol)
%    [bmask,uniqidx,nn,nc,totalarea,iso]=mcxsvmc(vol,'option1',v1,'option2',v2,...)
%
% Preprocessing voxelated labels to find and incorporate curved boundary
% information and improve accuracy
%
% Author: Qianqian Fang <q.fang at neu.edu>
%
% Input:
%    vol: a 3D volume with integer labels, a label 0 is assumed to be
%         background and tissue labels are non-zero positive numbers
%    options (optional): one can add additional 'name', value pairs to
%         specify user options, these options include
%         'debug': [0], if set to 1, print example surfaces for debugging
%         'smoothing': [1], if set to 1, do a Gaussian smoothing for the
%                 labels before extracting isosurface using marching cube
%         'kernelsize': [5], the 3D gaussian kernel size
%         'kernelstd': [2], the 3D gaussian kernel standard-deviation(sigma)
%         'threshold': [0.5], default threshold when extracting isosurfaces
%         'debugpatch': [1], ID of the patch to plot
%         'bmask': [], a pre-computed mask of the same size as vol; a
%                 voxel of non-zero value prevents adding new boundary info
%         'curveonly': [1] if set to 1, remove voxels where patch normals
%                 are align with x/y/z/ axes.
%
% Output:
%    bmask: an Nx2 array, where N is the length of uniqidx, with the 1st
%         column records the lower-valued label ID, and the 2nd column
%         records the higher-valued label ID in a mixed label voxel
%    uniqidx: the 1-D index of all voxels that are made of mixed labels
%    nc: nc is a Nx3 array, where N is the length of uniqidx, representing
%         the reference point of the intravoxel interface in a mixed label 
%         voxel
%    nn: nn is a Nx3 array, where N is the length of uniqidx, representing
%         the normalized normal vector in a mixed label voxel, pointing
%         from low to high label numbers
%    totalarea: an Nx1 array, where N is the length of uniqidx, representing
%         the total cross-sectional area of the boundary in each mixed
%         label voxel
%    iso: a struct, iso.vertices denotes the nodes and iso.faces denotes
%         the triangle patches of the combined isosurfaces from all labels.
%
%    if a single output is given, bmask, nn, nc, are merged into a single
%    4-D volume, where the indices in the 1st dimension represents
%    1-2: bmask, 3-5: nn, 6-8: nc
%
% Example:
%    [xi,yi,zi]=ndgrid(0.5:59.5,0.5:59.5,0.5:59.5);
%    dist=(xi-30).*(xi-30)+(yi-30).*(yi-30)+(zi-30).*(zi-30);
%    vol=zeros(size(dist));
%    vol(dist<10*10)=1;
%    vol(20:40,20:40,1:30)=2;
%    nuvox=mcxnuvoxel(vol);
%
% Dependency:
%    this function depends on the Iso2Mesh toolbox (http://iso2mesh.sf.net)
%
% This function is part of Monte Carlo eXtreme (MCX) URL: http://mcx.space
%
% License: GNU General Public License version 3, please read LICENSE.txt for details
%

%% parse user options
opt=varargin2struct(varargin{:});
dodebug=jsonopt('debug',0,opt);
ksize=jsonopt('kernelsize',3,opt);
kstd=jsonopt('kernelstd',1,opt);
level=jsonopt('threshold',0.5,opt);
debugpatch=jsonopt('debugpatch',1,opt);
bmask=jsonopt('bmask',NaN(size(vol)),opt);
dosmooth=jsonopt('smoothing',1,opt);
curveonly=jsonopt('curveonly',1,opt);

bmask2=zeros(size(vol));

%% read unique labels

labels=sort(unique(vol(:)));
% labels(labels==0)=[];

if(length(labels)>255)
    error('MCX currently supports up to 255 labels for this function');
end

%% loop over unique labels in ascending order

iso=struct('vertices',[],'faces',[]);

for i=1:length(labels)
    % convert each label into a binary mask, smooth it, then extract the
    % isosurface using marching cube algorithm (matlab builtin)
    if(dosmooth)
        volsmooth=smooth3(double(vol==labels(i)),'g',ksize,kstd);
    else
        volsmooth=(vol==labels(i));
    end
    [xi,yi,zi]=ndgrid(1:size(volsmooth,1),1:size(volsmooth,2),1:size(volsmooth,3));
    fv0=isosurface(xi,yi,zi,volsmooth,level);
    if(isempty(fv0.vertices))
        continue;
    end

    % get the containing voxel linear id
    c0=meshcentroid(fv0.vertices,fv0.faces);
    voxid=sub2ind(size(vol),floor(c0(:,1))+1,floor(c0(:,2))+1,floor(c0(:,3))+1);
    % identify unique voxels
    uniqidx=unique(voxid);
    bmask2(uniqidx)=labels(i);
    
    % find new boundary voxels that are not covered by previous levelsets
    uniqidx=uniqidx(isnan(bmask(uniqidx)));
    goodpatchidx=ismember(voxid,uniqidx);
    
    % merge surface patches located inside new boundary voxels
    iso.faces=[iso.faces; size(iso.vertices,1)+fv0.faces(goodpatchidx==1,:)];
    iso.vertices=[iso.vertices; fv0.vertices];
    
    % label those voxels as covered
    bmask(uniqidx)=labels(i);
end

%% handle uniform domains
if(isempty(iso.vertices))
    uniqidx=[];
    nn=[];
    totalarea=[];
    return;
end

%% get the voxel mapping for the final combined isosurface
[iso.vertices,iso.faces]=removeisolatednode(iso.vertices,iso.faces);
c0=meshcentroid(iso.vertices,iso.faces);
voxid=sub2ind(size(vol),floor(c0(:,1))+1,floor(c0(:,2))+1,floor(c0(:,3))+1);
[uniqidx, vox2patch]=unique(voxid);
    
[cc1,cc2]=histc(voxid,uniqidx);
cc=cc1(cc2);

if(dodebug)
    plotmesh(iso.vertices,iso.faces,'facealpha',0.4,'facecolor','b','edgealpha',0.2)
    disp(max(iso.vertices)-min(iso.vertices))
end

%% obtain the low and high labels in all mix-label voxels
bmask=[bmask(uniqidx),bmask2(uniqidx)];
bmask(bmask(:,1)==bmask(:,2),2)=0;

%% computing total area and normal vector for each boundary voxel
areas=elemvolume(iso.vertices,iso.faces);
normals=surfacenorm(iso.vertices,iso.faces);
centroids=meshcentroid(iso.vertices,iso.faces);

totalarea=zeros(size(vol));
maxvid=max(voxid);
totalarea(1:maxvid)=accumarray(voxid,areas); % total areas of cross-sections in each boundary voxel
% totalarea_unique=totalarea(uniqidx);

%% compute weighted average surface centroid per boundary voxel;
nc=zeros([3 size(vol)]);
nc(1,1:maxvid)=accumarray(voxid,centroids(:,1).*areas); % total normal_x of cross-sections in each boundary voxel
nc(2,1:maxvid)=accumarray(voxid,centroids(:,2).*areas); % total normal_y of cross-sections in each boundary voxel
nc(3,1:maxvid)=accumarray(voxid,centroids(:,3).*areas); % total normal_z of cross-sections in each boundary voxel
for i=1:3
    nc(i,uniqidx)=nc(i,uniqidx)./totalarea(uniqidx)';
end
nc=nc(:,uniqidx)';

%% creating weighted average surface patch normals per boundary voxel
nn=zeros([3 size(vol)]);
nn(1,1:maxvid)=accumarray(voxid,normals(:,1).*areas); % total normal_x of cross-sections in each boundary voxel
nn(2,1:maxvid)=accumarray(voxid,normals(:,2).*areas); % total normal_y of cross-sections in each boundary voxel
nn(3,1:maxvid)=accumarray(voxid,normals(:,3).*areas); % total normal_z of cross-sections in each boundary voxel
nnlen=sqrt(sum(nn.*nn,1));
for i=1:3
    nn(i,uniqidx)=nn(i,uniqidx)./nnlen(uniqidx)';
end
nn=nn(:,uniqidx)';

%% remove x/y/z oriented 
if(curveonly)
    [ixc,iyc]=find((nc-floor(nc))<1e-6);
    [ixn,iyn]=find(abs(nn)==1);
    ix=intersect(ixc,ixn);
    boxmask=zeros(size(vol));
    boxmask(uniqidx(ix))=1;
    boxmask=smooth3(boxmask,'b',3);
    boxmask=(boxmask>0);
    ix=find(boxmask(uniqidx)==1);

    nn(ix,:)=[];
    nc(ix,:)=[];
    uniqidx(ix)=[];
    bmask(ix,:)=[];
    totalarea(ix)=[];
end

%% discretize nn and nc vector components to 0-255 gray-scale numbers

nc=nc-floor(nc);
nc=floor(nc*255);
nn=min(floor((nn+1)*255/2),254);

%% assemble the final volume
if(nargout==1)
    newvol=zeros([8,size(vol)]);
    newvol(1,:,:,:)=vol;
    newvol(1:2,uniqidx)=bmask';
    newvol(3:5,uniqidx)=nc';
    newvol(6:8,uniqidx)=nn';
    bmask=newvol;
end

%% plotting for verification

if(dodebug)
    figure;

    pcidx=find(cc>1);         % find 4-patch that blong to the same voxel (in patch idx)
    pidx0=pcidx(debugpatch);          % pick one such patch group to debug
    idx1=voxid(pidx0);        % find the corresponding voxel id idx1 for the patch pidx0
    patid=voxid==idx1;  % these patches are within voxel linear id idx1

    testmask=zeros(size(volsmooth));
    testmask(idx1)=1;
    [no3,fc3]=binsurface(testmask,4);

    figure;
    plotmesh(no3,fc3,'facealpha',0.3,'facecolor','none')
    hold on;
    plotmesh(iso.vertices,iso.faces(patid,:),'facealpha',0.3,'facecolor','b','edgealpha',0.2)
    disp(nn(cc2(pidx0),:))           % normal in voxel id idx1
    plotmesh([c0(pidx0,:); c0(pidx0,:)+nn(cc2(pidx0),:)],'ro-');   % plot centroid of pidx0-th patch to the normal in voxel id idx1
end