
function [R, parent] = tree(f)    
    parent = zpar = nan(size(f));

    [_, s] = sort(f(:), "descend");
    [_, R] = intersect(s, [1:numel(f)]);
    
    for p = s'
        parent(p) = p;
        zpar(p) = p;
        
        for n = N(p, f)
            if isnan(zpar(n))
                continue;
            end
            
            [r, zpar] = root(n, zpar);
            
            if r != p
                parent(r) = p;
                zpar(r) = p;
            end
        end
    end
    
    R = reshape(R, size(f));
end
