function [delta_x, delta_y, delta_xy, cont_error, seg_lengths] = get_tree_statistics(tree, max_deltas, lossless_cont)
% -------------------------------------------------------------------------
% MSc Project - MSc in Communications and Signal Processing
% Department of Electrical and Electronic Engineering
% Imperial College London, 2011
%
% Supervisor  : Dr Pier Luigi Dragotti
% Contributor : Andriy Gelman
% Author      : Jon Onativia Bravo
%
% File        : get_tree_statistics.m
% -------------------------------------------------------------------------
% 
%
% [delta_x, delta_y, delta_xy, error, seg_lengths] = get_tree_statistics(tree, max_deltas, lossless_cont)
%
% INPUT:
%  - 
%
% OUTPUT:
%  - 
%

% Obtain tree parameters
num_layers      = length(tree);
max_depth       = length(tree{1}{1});

% Allocate output variables
delta_x     = cell(max_depth, 1);
delta_y     = cell(max_depth, 1);
delta_xy    = cell(max_depth, 1);
cont_error  = cell(max_depth, 1);
seg_lengths = cell(max_depth, 1);

% Reconstruct the original lossless contour
cont_inds      = cell(num_layers, 1);
cont_lengths   = cell(num_layers, 1);
for ith_layer = 1 : num_layers
    cont_inds{ith_layer}    = find(lossless_cont{ith_layer}(1,:) == 0);
    cont_lengths{ith_layer} = lossless_cont{ith_layer}(2,cont_inds{ith_layer});
end

for ith_level = 1 : max_depth
    cur_delta    = max_deltas(ith_level);
    
    cur_delta_x  = zeros(2*cur_delta+1, 2);
    cur_delta_y  = zeros(2*cur_delta+1, 2);
    cur_delta_xy = zeros(2*cur_delta+1, 2*cur_delta+1);
    
    cur_delta_x(:,1) = (-cur_delta:cur_delta)';
    cur_delta_y(:,1) = (-cur_delta:cur_delta)';
    
    % Obtain the maximum segment length and the total number of vertices
    max_length  = 0;
    total_verts = 0;
    for ith_layer = 1 : num_layers
        layer_tree   = tree{ith_layer};
        num_contours = length(layer_tree);
        
        for ith_cont = 1 : num_contours
            
            if ~iscell(layer_tree{ith_cont}{ith_level})
                total_verts = total_verts + size(layer_tree{ith_cont}{ith_level}, 2) - 1;
                
                verts      = layer_tree{ith_cont}{ith_level};
                max_length = max([max_length size(verts,2)]);
            else
                total_verts = total_verts + size(cell2mat(layer_tree{ith_cont}{ith_level}), 2) - 1;
                
                num_segments = length(layer_tree{ith_cont}{ith_level});
                for ith_seg = 1 : num_segments
                    verts      = layer_tree{ith_cont}{ith_level}{ith_seg};
                    max_length = max([max_length size(verts,2)]);
                end
            end
            
        end % for ith_cont
    end % for ith_layer

    % Allocate the segment length vector
    cur_cont_error = zeros(total_verts, 1);
    cur_seg_len    = zeros(max_length, 1);
    
    % Compute statistics
    idx = 1;
    for ith_layer = 1 : num_layers
        layer_tree   = tree{ith_layer};
        num_contours = length(layer_tree);
        
        for ith_cont = 1 : num_contours
            
            % Retrieve the vertices of the current level
            if iscell(layer_tree{ith_cont}{ith_level})
                cur_verts = cell2mat(layer_tree{ith_cont}{ith_level});
            else
                cur_verts = layer_tree{ith_cont}{ith_level};
            end
            
            % Retrieve the original contour
            cont_start = cont_inds{ith_layer}(ith_cont) + 1;
            cont_end   = cont_inds{ith_layer}(ith_cont) + cont_lengths{ith_layer}(ith_cont);
            verts      = lossless_cont{ith_layer}(:,cont_start:cont_end);
            verts_t    = verts';
            
            len = size(cur_verts, 2);
            for ith_vert = 2 : len
                first_vert = cur_verts(:,ith_vert-1);
                last_vert  = cur_verts(:,ith_vert);
                
                d_x = last_vert(1) - first_vert(1);
                d_y = last_vert(2) - first_vert(2);

                % Retrieve the indices
                i_x = find(cur_delta_x(:,1) == d_x);
                i_y = find(cur_delta_y(:,1) == d_y);

                % Increase the corresponding values
                cur_delta_x(i_x,2)    = cur_delta_x(i_x,2) + 1;
                cur_delta_y(i_y,2)    = cur_delta_y(i_y,2) + 1;
                cur_delta_xy(i_x,i_y) = cur_delta_xy(i_x,i_y) + 1;
                
                % Measure the error of the approximated contour
                first_vert_ind = find(ismember(verts_t, first_vert', 'rows'));
                last_vert_ind  = find(ismember(verts_t, last_vert', 'rows'));
                
                % Extract the sub_segment between the two vertices
                if first_vert_ind < last_vert_ind
                    sub_seg = verts(:,first_vert_ind:last_vert_ind);
                else
                    sub_seg = [verts(:,first_vert_ind:end) verts(:,1)];
                end
                cur_cont_error(idx) = get_max_distance(first_vert, last_vert, sub_seg);
                idx = idx + 1;
            end % for ith_vert
            
            % Measure the segments lengths
            if ~iscell(layer_tree{ith_cont}{ith_level})
                len = size(layer_tree{ith_cont}{ith_level}, 2);
                cur_seg_len(len) = cur_seg_len(len) + 1;
            else
                num_segments = length(layer_tree{ith_cont}{ith_level});
                for ith_seg = 1 : num_segments
                    verts      = layer_tree{ith_cont}{ith_level}{ith_seg};
                    len        = size(verts, 2);
                    cur_seg_len(len) = cur_seg_len(len) + 1;
                end % for ith_seg
            end
            
        end % for ith_cont
    end % for ith_layer
    
    delta_x{ith_level}     = cur_delta_x;
    delta_y{ith_level}     = cur_delta_y;
    delta_xy{ith_level}    = cur_delta_xy;
    cont_error{ith_level}  = cur_cont_error;
    seg_lengths{ith_level} = cur_seg_len;
end % for ith_level
