function [approx_layers, total_segmentation_cost] = lossy_segmentation_encoding_msc(layer_stack, target_segmentation_rate)
% -------------------------------------------------------------------------
% 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        : lossy_segmentation_encoding_msc.m
% -------------------------------------------------------------------------
% 
%
% [layer_stack, total_segmentation_cost] = lossy_segmentation_encoding_msc(layer_stack, target_segmentation_rate)
%
% INPUT:
%  - 
%
% OUTPUT:
%  - 
%

% Tree parameters
max_distances = [ 3  2 1];
max_offsets   = [31 15 7];
max_depth     = length(max_distances);

% Labels parameters
[height, width, num_layers] = size(layer_stack);
approx_layers               = zeros(height, width, num_layers);

% Compute the contour trees for each layer
label_verts_tree = cell(num_layers, 1);
for l = 1 : num_layers
    if l < num_layers
        label_verts_tree{l} = estimate_contours_tree(layer_stack(:,:,l), max_distances, max_offsets);
    else
        last_layer          = padarray(layer_stack(:,:,l), [1 1]);
        label_verts_tree{l} = estimate_contours_tree(last_layer, [20 20 20], max_offsets);
    end
end % end of l = 1 : num_layers

% Obtain the bitstreams
[bitstream, layer_limits] = tree2bitstream_v3(label_verts_tree, max_offsets, width, height);

field_bits = log2(max_offsets + 1) + 1;

% Compress the bitstream
[substreams] = arithm_enc_stream(bitstream, layer_limits, field_bits);

required_bits = zeros(1, max_depth);
for ith_level = 1 : max_depth
    for remaining_levels = ith_level : max_depth
        required_bits(remaining_levels) ...
            = required_bits(remaining_levels) + length(substreams{ith_level})*8;
    end
end

% Find the level that fits best the target_segmentation_rate
req_depth               = 1;
total_segmentation_cost = required_bits(1);
for ith_level = 1 : max_depth
    if required_bits(ith_level) > target_segmentation_rate
        break;
    else
        req_depth               = ith_level;
        total_segmentation_cost = required_bits(ith_level);
    end
end

% Reconstruct the layers for this given bitrate
for l = 1 : num_layers
    if l < num_layers
        verts                = tree2contour(label_verts_tree{l}, req_depth);
        approx_layers(:,:,l) = contour2shape(verts, width, height);
    else
        verts                = tree2contour(label_verts_tree{l}, req_depth);
        last_layer           = contour2shape(verts, width+2, height+2);
        approx_layers(:,:,l) = last_layer(2:end-1,2:end-1);
    end
end

