function [dbratio,uniqidx,bmask,nn,totalarea,iso]=mcxnuvoxel(vol, varargin)
%
% Format:
%    newvol=mcxnuvoxel(vol)
%       or
%    [dbratio,uniqidx,nn,totalarea,iso]=mcxnuvoxel(vol)
%    [dbratio,uniqidx,nn,totalarea,iso]=mcxnuvoxel(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:
%    dbratio: a vector of length as uniqidx, the approx. partial volume of 
%         the higher-valued label in a mixed label voxel (between 0-1)
%    uniqidx: the 1-D index of all voxels that are made of mixed labels
%    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
%    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, dbratio, bmask, nn, and totalarea are
%    merged into a single 4-D volume, where the indices in the 1st
%    dimension represents
%      1: dbratio, 2-3: bmask, 4-6: nn, 7: totalarea
%
% 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',zeros(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(vol==labels(i),'g',ksize,kstd);
    else
        volsmooth=(vol>=labels(i));
    end
    fv0=isosurface(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(bmask(uniqidx)==0);
    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))
    dbratio=[];
    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);

totalarea=zeros(size(vol));
maxvid=max(voxid);
totalarea(1:maxvid)=accumarray(voxid,areas); % total areas of cross-sections in each boundary voxel
totalarea=totalarea(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,voxid)=nn(i,voxid)./nnlen(voxid)';
end
nn=nn(:,uniqidx)';

%% calculating distances between 8 nodes of the enclosing voxel and the
% centroid of the surface patch
p0=[floor(c0(vox2patch,1)) floor(c0(vox2patch,2)) floor(c0(vox2patch,3))];
p0=repmat(p0,1,8)+repmat([0 0 0  1 0 0  0 1 0  0 0 1  1 1 0  1 0 1  0 1 1  1 1 1],size(p0,1),1);

pdiff=reshape(p0-repmat(c0(vox2patch,:),1,8), length(vox2patch), 3, 8);
pdiff=squeeze(sum(pdiff.*repmat(nn,[1,1,8]),2)); % pdiff=dr .* n: unique voxels x3

dbvox=zeros(size(pdiff,1),2);
[s1,s2]=find(pdiff>0);
dbvox(:,1)=accumarray(s1, pdiff(pdiff>0));    % summing positive distances
[s1,s2]=find(pdiff<0);
dbvox(:,2)=-accumarray(s1, pdiff(pdiff<0));   % dbvox: unique voxels x2, summing negative distances

dbratio=dbvox(:,2)./(dbvox(:,1)+dbvox(:,2));  % dbratio: unique voxels

%% remove x/y/z oriented 
if(curveonly)
    [ix,iy]=find(abs(nn)==1);
    boxmask=zeros(size(vol));
    boxmask(uniqidx(ix))=1;
    boxmask=smooth3(boxmask,'b',3);
    boxmask=(boxmask>0);
    ix=find(boxmask(uniqidx)==1);

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

%% assemble the final volume

if(nargout==1)
    newvol=zeros([7,size(vol)]);
    newvol(1,uniqidx)=dbratio;
    newvol(2:3,uniqidx)=bmask';
    newvol(4:6,uniqidx)=nn';
    newvol(7,uniqidx)=totalarea';
    dbratio=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=find(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
    disp(dbratio(cc2(patid(1))))     % mix ratio of the selected voxel
end