function [contours] = estimate_contours_tree(bin_image, max_errors, ...
                                             max_vert_deltas)
% -------------------------------------------------------------------------
% 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        : estimate_contours_tree.m
% -------------------------------------------------------------------------
% Estimate contour vertices for a given binary shape allowing a maximum
% distance of max_errors(i) between the real contour and the approximating
% polygon and a maximum Dx and Dy of max_vert_deltas(i) between 
% consecutive vertices for each decomposition level.
%
% [contours] = estimate_contours_tree(bin_image, max_errors, ...
%                                     max_vert_deltas)
%
% INPUT:
%  - bin_image       : Binary shape image
%  - max_errors      : Maximum allowed distance between real contour and
%                      approximated polygon. Each element of the vector
%                      corresponds to the allowed distance for each level.
%  - max_vert_deltas : Maximum offset in pixels for X and Y coordinates
%                      between two consecutive vertices. The vector has to 
%                      be of the same size as max_errors argument.
%
% OUTPUT:
%  - contours        : Tree of contours containing the trees of vertices.
%


%% Read input parameters
num_of_levels = length(max_errors);


%% Binarize bin_image image
bin_image(bin_image ~= 0) = 1;


%% Find the contour
figure
[c, ~] = contour(bin_image, [0 0]);
close gcf;
c = round(c);


%% Determine the number of segments.
% The c vector is 2xN matrix, with first column [0; len_of_segment] and
% then len_of_segment columns  with first row contains x coordinate 
% and second row contains y coordinate of the points of the contour. If
% there are more contours another [0; len_of_segment] element will appear.
segment_inds = find(c(1,:) == 0);
seg_lengths  = c(2,segment_inds);
num_of_segs  = length(segment_inds);


%% Remove duplicated values
c1 = [];
for i_seg = 1 : num_of_segs
    % Retrieve the segment
    seg_len   = seg_lengths(i_seg);
    seg_start = segment_inds(i_seg) + 1;
    seg_end   = seg_start + seg_len - 1;
    seg       = c(:,seg_start:seg_end);
    
    % Remove duplicated values
%     [~, l]    = unique(seg', 'rows');
%     seg       = seg(:,sort(l));
    new_seg = seg(:,1);
    for i = 2 : seg_len-1
        if seg(1,i) ~= new_seg(1,end) || seg(2,i) ~= new_seg(2,end)
            new_seg = [new_seg, seg(:,i)];
        end
    end
    % Insert the last vertex only if it's different from previous and first
    if (seg(1,end) ~= new_seg(1,end) || seg(2,end) ~= new_seg(2,end)) ...
            && (seg(1,end) ~= new_seg(1,1)  || seg(2,end) ~= new_seg(2,1))
    	new_seg = [new_seg, seg(:,end)];
    end
    % Remove last vertex if it is equal to the first one
    if size(new_seg,2) > 1 ...
            && new_seg(1,end) == new_seg(1,1) ...
            && new_seg(2,end) == new_seg(2,1)
        new_seg = new_seg(:,1:end-1);
    end
    seg = new_seg;
    
    % Insert segment in vertices array
    c1  = [c1, [0; length(seg(1,:))]];
    c1  = [c1, seg];
end


%% Get the different new segments
segment_inds = find(c1(1,:) == 0);
seg_lengths  = c1(2,segment_inds);
num_of_segs  = length(segment_inds);


%% Estimate vertices

% The number of contours will determine the number of cells of the output
% parameter 'contours'
contours = cell(1, num_of_segs);

for i_seg = 1 : num_of_segs
    
    % For each contour, a tree of vertices of depth num_of_levels will be
    % created
    vertices = cell(num_of_levels, 1);

    % Retrieve the segment of the current contour
    seg_len   = seg_lengths(i_seg);
    seg_start = segment_inds(i_seg) + 1;
    seg_end   = seg_start + seg_len - 1;
    seg       = c1(:,seg_start:seg_end);
    seg_t     = seg';

    % Check if it is an open or a closed contour
    contour_is_closed = true;
    if abs(seg(1,end) - seg(1,1)) > 1 || abs(seg(2,end) - seg(2,1)) > 1
        contour_is_closed = false;
    end
    
    % Compute the first level of decomposition
    max_offset = max_vert_deltas(1);
    max_distance = max_errors(1);
    
    verts = find_vertices(seg, max_distance, max_offset);
    
    if ~isempty(verts)
        vertices{1} = verts;
    end

    % Continue in the next levels of decomposition
    for l = 2 : num_of_levels
        
        max_distance = max_errors(l);
        max_offset   = max_vert_deltas(l);
        
        prev_verts        = vertices{l-1};
        if iscell(prev_verts)
            prev_verts    = cell2mat(prev_verts);
        end
        num_of_prev_verts = size(prev_verts, 2);
        vertices{l}       = cell(1, num_of_prev_verts);
        
        for i = 1 : num_of_prev_verts
            
            first_vert = prev_verts(:,i);
            if i < num_of_prev_verts
                last_vert = prev_verts(:,i+1);
            else
                if contour_is_closed == true
                    last_vert = prev_verts(:,1);
                else
                    % If contour is open, insert the last vertex and stop
                    vertices{l}{i} = first_vert;
                    break;
                end
            end
            
            % Find the indices of the extreme vertices in the original
            % segment
            first_vert_ind = find(ismember(seg_t, first_vert', 'rows'));
            last_vert_ind  = find(ismember(seg_t, last_vert', 'rows'));
            
            % Extract the sub_segment between the two vertices
            if first_vert_ind < last_vert_ind
                sub_seg = seg(:,first_vert_ind:last_vert_ind);
            else
                sub_seg = [seg(:,first_vert_ind:end) seg(:,1)];
            end
            verts = find_vertices(sub_seg, max_distance, max_offset);
            
            % Remove the last vertex, only if it's not the last vertex
            if size(verts, 2) > 1
                if verts(1,end) == last_vert(1,1) ...
                        && verts(2,end) == last_vert(2,1)
                    verts = verts(:, 1:end-1);
                end
            end
            
            % Insert the vertices in the tree
            vertices{l}{i} = verts;
        end
        
    end

    contours{i_seg} = vertices;
end
    

end



%%
function vertices = find_vertices(seg, max_distance, max_offset)

cur_vert   = 1;
next_vert  = 2;

seg_len    = length(seg(1,:));

if seg_len == 1
    new_seg = seg;
else
    new_seg = [];
    while cur_vert < seg_len

        new_seg = [ new_seg, [seg(1,cur_vert); seg(2,cur_vert)] ];

        if next_vert < seg_len
            d = get_max_distance([seg(1,cur_vert); seg(2,cur_vert)], ...
                           [seg(1,next_vert+1); seg(2,next_vert+1)], ...
                           seg(:,cur_vert+1:next_vert));
        else
            break
        end

        while d < max_distance ...
              && next_vert+1 < seg_len ...
              && abs(seg(1,next_vert+1)-seg(1,cur_vert)) <= max_offset ...
              && abs(seg(2,next_vert+1)-seg(2,cur_vert)) <= max_offset
            next_vert = next_vert + 1;
            d = get_max_distance([seg(1,cur_vert); seg(2,cur_vert)], ...
                           [seg(1,next_vert+1); seg(2,next_vert+1)], ...
                           seg(:,cur_vert+1:next_vert));
        end

        cur_vert  = next_vert;
        next_vert = next_vert + 1;

    end % end of while
end

% If the original contour is not closed (the contour ends in the
% border of the image), force the insertion of the last vertex so
% the contour contains the vertices of the image border.
if abs(seg(1,end) - seg(1,1)) > 1 || abs(seg(2,end) - seg(2,1)) > 1
    % Check if the last vertex of the original contour has already
    % been inserted
    if seg(1,end) ~= new_seg(1,end) || seg(2,end) ~= new_seg(2,end)
        new_seg = [ new_seg, [seg(1,end); seg(2,end)] ];
    end
end

vertices = new_seg;

% Check that the max delta between vertices satisfies the restriction
max_delta = get_max_offset([ [0; size(vertices, 2)] , vertices ]);
if max_delta > max_offset
    disp('The maximum delta is bigger than the maximum allowed')
end

end
%%