function [cut, map, short] = vl_aibcut(parents, n)
% VL_AIBCUT  Cut VL_AIB tree
%  CUT = VL_AIBCUT(PARENTS, N) cuts the binary merge tree PARENTS and
%  returns a cut CUT of N nodes. The format of PARENTS is the same
%  used by the VL_AIB() function.
%
%  A cut is a set of N nodes such that no node is a descendant of any
%  other node in the CUT and such that any leaf descend from a node in
%  the cut. CUT lists the nodes of the binary merge tree PARENT that
%  form the cut.
%
%  Nodes with null parent (as defined by PARENTS) are comprised in the
%  cut if the other nodes are not enough to fill a cut of N elements.
%
%  [CUT, MAP] = VL_AIBCUT(...) returns a vector MAP with the same size as
%  PARENTS. MAP assign each node below or in the cut to the
%  corresponding element in the CUT vector (each element above the cut
%  or with null parent is mapped to 0). To get the index of the
%  corresponding cut nodes use CUT(MAP). MAP by itself is useful to
%  quantize the leaves in a sequences of N contiguous indexes,
%  starting from one (see also VL_AIBCUTPUSH()).
%
%  [CUT, MAP, SHORT] = VL_AIBCUT(...) returns also a vector SHORT that
%  represents a merge tree compressed to the cut. This is obtained by
%  mapping all nodes below to the cut element above them. Nodes in or
%  above the cut are mapped to themselves. Null parents are left
%  unchanged, except if the corresponding node is in the cut (in which
%  case the map-to-itself rule has the precedence).
%
%  See also: VL_HELP(), VL_AIB().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% --------------------------------------------------------------------
%                                           Determine nodes in the cut
% --------------------------------------------------------------------

if n > 1
  root    = max(parents) ;

  % count number of null nodes
  z = sum(parents(1:root) == 0) ;

  % determine number of leves
  nleaves = (root - z + 1) / 2 ;

  % find first node of the cut
  mu   = root - min(n, nleaves) + 1 ;

  % correction for presence of null nodes
  nz   = find(parents(1:mu) > 0) ;
  mu   = nz(end) ;

  % find node belnoging to the cut
  cut  = find(parents(1:mu) > mu) ;

  % In the presence of null nodes, the cut size might exceed
  % nleaves, which is the maximum cut size we can obtain with the
  % specified tree. The additional nodes have to be picked up from
  % the null nodes.

  if length(cut) < n
    sel_z = find(parents == 0) ;
    cut = [sel_z(1:n-length(cut)) cut] ;
  end

  % aesthetic reasons only
  cut = sort(cut) ;

else
  mu   = max(parents) ;
  cut  = mu ;
end

% --------------------------------------------------------------------
%                                       Short-circuit nodes to the cut
% --------------------------------------------------------------------

stop = [cut find(parents == 0)] ;
short = 1:length(parents) ;

while 1
  [drop,sel] = setdiff(short(1:mu), stop)  ;
  sel = setdiff(sel, stop) ;
  if isempty(sel), break ; end
  short(sel) = parents(short(sel))  ;
end

short(setdiff(find(parents == 0), cut)) = 0 ;

% --------------------------------------------------------------------
%                                                  Build quantizer map
% --------------------------------------------------------------------

map             = 1:numel(parents) ;
map(cut)        = 1:n ;
map(short >  0) = map(short(short > 0)) ;
map(short == 0) = 0 ;
map(mu+1:end)   = 0 ;
