maketree.m

function [tr] = maketree(rp,varargin)
%MAKETREE assemble an AABB search tree for a collection of 
%(hyper-)rectangles.
%   [TR] = MAKETREE(RP) returns an axis-aligned bounding-box 
%   (AABB) tree for a collection of d-rectangles embedded in 
%   R^d. The rectangles are defined as an NR-by-NDIM*2 array 
%   of min/max coordinates RP = [PMIN,PMAX], where PMIN = 
%   [A1,A2,...,ANDIM] and PMAX = [B1,B2,...,BNDIM] are the 
%   minimum/maximum coordinates associated with each rectan-
%   gle. 
%
%   The resulting tree is returned as a structure containing
%   an NT-by-NDIM*2 array of tree-node coordinates TR.XX = 
%   [PMIN,PMAX], an NT-by-2 array of tree-node indexing 
%   TR.II = [PI,CI], and an NT-by-1 cell array of node lists
%   TR.LL. PI,CI are the parent/child pointers associated 
%   with each node in the tree, and TR.LL{II} is the list of
%   rectangles associated with the II-th node.
%
%   MAKETREE forms a binary search tree, such that each 
%   "leaf" node in the tree encloses a maximum number of re-
%   ctangles. The tree is formed by recursively subdividing 
%   the bounding-box of the collection. At each division, a
%   simple heuristic is used to determine a splitting axis
%   and to position an axis-aligned splitting (hyper-)plane. 
%   The associated collection of rectangles is partitioned 
%   between two new child nodes. The dimensions of each node 
%   in the tree are selected to provide a minimal enclosure 
%   of the rectangles in its associated sub-tree. Tree nodes 
%   may overlap as a result. 
%
%   [...] = MAKETREE(RP,OP) also passes an additional user-
%   defined options structure. OP.NOBJ = {32} is the maximum 
%   allowable number of rectangles per tree-node. OP.LONG = 
%   {.75} is a relative length tolerance for "long" rectang-
%   les, such that any rectangles with RMAX(IX)-RMIN(IX) >=
%   OP.LONG * (NMAX(IX)-NMIN(IX)) remain in the parent node.
%   Here RMIN,RMAX are the coordinates of the rectangle, 
%   NMIN,NMAX are the coordinates of the enclosing node in
%   the tree, and IX is the splitting axis. Nodes that beco-
%   me "full" of "long" items may exceed their maximum rect-
%   angle capacity. OP.VTOL = {.55} is a "volume" splitting
%   criteria, designed to continue subdivision while the net
%   node volume is reducing. Specifically, a node is split 
%   if V1+V2 <= OP.VTOL*VP, where VP is the d-dim. "volume"
%   of the parent node, and V1,V2 are the volumes associated 
%   with its children.
%
%   See also DRAWTREE, QUERYSET, MAPVERT, MAPRECT

% Please see the following for additional information:
%
%   Darren Engwirda, "Locally-optimal Delaunay-refinement & 
%   optimisation-based mesh generation". Ph.D. Thesis, Scho-
%   ol of Mathematics and Statistics, Univ. of Sydney, 2014:
%   http://hdl.handle.net/2123/13148

%   Darren Engwirda : 2014 --
%   Email           : de2363@columbia.edu
%   Last updated    : 08/10/2017

    tr.xx = []; tr.ii = []; tr.ll = {}; op = [];
    
%------------------------------ quick return on empty inputs
    if (isempty(rp)), return; end
    
%---------------------------------------------- basic checks    
    if (~isnumeric(rp))
        error('maketree:incorrectInputClass', ...
            'Incorrect input class.');
    end
    
%---------------------------------------------- basic checks
    if (nargin < +1 || nargin > +2)
        error('maketree:incorrectNoOfInputs', ...
            'Incorrect number of inputs.');
    end
    if (ndims(rp) ~= +2 || ...
            mod(size(rp,2),2)~= +0)
        error('maketree:incorrectDimensions', ...
            'Incorrect input dimensions.');
    end
    
%------------------------------- extract user-defined inputs
    if (nargin >= +2), op = varargin{1}; end
    
    isoctave = exist( ...
        'OCTAVE_VERSION','builtin') > 0;
    
    if (isoctave)
    
    %-- faster for OCTAVE with large tree block size; slower
    %-- loop execution...
    
        NOBJ = +1024 ; 
        
    else
    
    %-- faster for MATLAB with small tree block size; better
    %-- loop execution... 
    
        NOBJ =  + 32 ;
        
    end
    
%--------------------------------------- user-defined inputs
    if (~isstruct(op))
        
        op.nobj = NOBJ ;
        op.long = .75;
        op.vtol = .55;
        
    else
    
    %-------------------------------- bound population count
        if (isfield(op,'nobj'))
            if (op.nobj <= +0 )
                error('Invalid options.') ;
            end
        else
            op.nobj = NOBJ ;
        end
        
    %-------------------------------- bound "long" tolerance
        if (isfield(op,'long'))
            if (op.long < +.0 || op.long > +1.)
                error('Invalid options.') ;
            end
        else
            op.long = .75; 
        end
        
    %-------------------------------- bound "long" tolerance
        if (isfield(op,'vtol'))
            if (op.vtol < +.0 || op.vtol > +1.)
                error('Invalid options.') ;
            end
        else
            op.vtol = .55; 
        end
        
    end
    
%---------------------------------- dimensions of rectangles
    nd = size(rp,2) / +2 ;
    ni = size(rp,1) ;
    
%------------------------------------------ alloc. workspace
    xl = zeros(ni*1,1*nd);
    xr = zeros(ni*1,1*nd);
    ii = zeros(ni*1,2);
    ll = cell (ni*1,1);
    ss = zeros(ni*1,1);
    
%------------------------------------ min & max coord. masks
    lv = false(size(rp,2),1);
    rv = false(size(rp,2),1);
    lv((1:nd)+nd*+0) = true ;
    rv((1:nd)+nd*+1) = true ;

%----------------------------------------- inflate rectangle
    r0 = min(rp(:,lv),[],1) ;
    r1 = max(rp(:,rv),[],1) ;
   
    rd = repmat(r1-r0,ni,1) ;
    
    rp(:,lv) = ...
    rp(:,lv) - rd * eps^.8;
    rp(:,rv) = ...
    rp(:,rv) + rd * eps^.8;
    
%----------------------------------------- rectangle centres
    rc = rp(:,lv)+rp(:,rv);
    rc = rc * .5 ;
%----------------------------------------- rectangle lengths
    rd = rp(:,rv)-rp(:,lv);

%------------------------------ root contains all rectangles
    ll{1} = (+1:ni)' ;
%------------------------------------ indexing for root node
    ii(1,1) = +0 ;
    ii(1,2) = +0 ;
%------------------------------ root contains all rectangles
    xl(1,:) = min(rp(:,lv),[],1);
    xr(1,:) = max(rp(:,rv),[],1);
    
%-- main loop : divide nodes until all constraints satisfied
    ss(+1) = +1; ns = +1; nn = +1;   
    while (ns ~= +0)
    %----------------------------------- pop node from stack
        ni = ss(ns) ; 
        ns = ns - 1 ;
    %----------------------------------- push child indexing
        n1 = nn + 1 ; 
        n2 = nn + 2 ;           
        
    %--------------------------- set of rectangles in parent
        li = ll{ni} ;    
    %--------------------------- split plane on longest axis
        dd = xr(ni,:) ...
           - xl(ni,:) ;
       [dd,ia] = sort(dd);
  
        for id = nd : -1 : +1
    %--------------------------- push rectangles to children         
            ax = ia (id) ;
            mx = dd (id) ;
    
            il = rd(li,ax) > ...
                op.long * mx ;
            lp = li( il) ;          %  "long" rectangles
            ls = li(~il) ;          % "short" rectangles
        
            if (length(lp) < ...
            0.5*length(ls)&& ...
                length(lp) < ...
            0.5 * op.nobj)
                break ;
            end 
        end
        
        if (isempty(ls) )
    %-------------------------------- partition empty, done!
            continue ;
        end
        
    % select the split position: take the mean of the set of
    % (non-"long") rectangle centres along axis AX
    %-------------------------------------------------------
        sp = sum(rc(ls,ax))/length(ls);
       
    %---------------------------- partition based on centres
        i2 = rc(ls,ax)>sp ;
        l1 = ls(~i2) ;              %  "left" rectangles
        l2 = ls( i2) ;              % "right" rectangles
    
        if (isempty(l1) || ...
            isempty(l2) )
    %-------------------------------- partition empty, done!
            continue ;
        end
        
    %-------------------------------- finalise node position
        xl(n1,:) = ...
            min(rp(l1,lv),[],1) ;
        xr(n1,:) = ...
            max(rp(l1,rv),[],1) ;
        xl(n2,:) = ...
            min(rp(l2,lv),[],1) ;
        xr(n2,:) = ...
            max(rp(l2,rv),[],1) ;
            
    %--------------------------- push child nodes onto stack        
        if (length(li) <= op.nobj )
        
            vi = prod(xr(ni,:) ...  % upper d-dim "vol."
                    - xl(ni,:) ) ;
            v1 = prod(xr(n1,:) ...  % lower d-dim "vol."
                    - xl(n1,:) ) ; 
            v2 = prod(xr(n2,:) ...
                    - xl(n2,:) ) ;
        
            if (v1+v2 < op.vtol*vi)
               
    %-------------------------------- parent--child indexing
            ii(n1,1) = ni ; 
            ii(n2,1) = ni ;
            ii(ni,2) = n1 ;
            
    %-------------------------------- finalise list indexing
            ll{ni,1} = lp ;
            ll{n1,1} = l1 ;
            ll{n2,1} = l2 ;
        
            ss(ns+1) = n1 ;
            ss(ns+2) = n2 ;
            ns = ns+2; nn = nn+2;
            
            end
            
        else
    %-------------------------------- parent--child indexing
            ii(n1,1) = ni ; 
            ii(n2,1) = ni ;
            ii(ni,2) = n1 ;
            
    %-------------------------------- finalise list indexing
            ll{ni,1} = lp ;
            ll{n1,1} = l1 ;
            ll{n2,1} = l2 ;
        
            ss(ns+1) = n1 ;
            ss(ns+2) = n2 ;
            ns = ns+2; nn = nn+2;
        
        end
        
    end
%----------------------------------------------- trim alloc.
    xl = xl(1:nn,:) ;
    xr = xr(1:nn,:) ;
    ii = ii(1:nn,:) ;
    ll = ll(1:nn,:) ;
    
%----------------------------------------------- pack struct
    tr.xx =[xl, xr] ;     
    tr.ii = ii ; 
    tr.ll = ll ;
    
end