% -------------------------------------------------------------------------
% 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       : compute_freeman_huffman.m
% -------------------------------------------------------------------------


clear all
close all
clc


% MATLAB data files
files = {'animal_farm.mat', 'cones.mat', 'teddy.mat', 'tsukuba_extracted4.mat'};
% files = {'tsukuba_extracted4.mat'};


% Loop over all files
for file_name_cell = files
    
    file_name         = file_name_cell{1};
    [~, data_name, ~] = fileparts(file_name);

    % Load data

    % Disparity vector
    load(file_name, 'D');
    % Images array
    load(file_name, 'in');
    % Labels array
    load(file_name, 'label');

    % Size of the data
    [r, c, num_images] = size(in);

    % Extract layers and merge

    labels     = unique(label(:,:,:));
    labels     = labels(labels > 0);
    num_layers = length(labels);
    disparity  = ceil( (0:num_images-1)' * D' );

    % The disparity measures the displacement of a layer with respect to the
    % first image (EPI images, camera only moves in the horizontal axis)
    pos_disp = max(disparity(:));
    neg_disp = -min(disparity(:));
    ext_cols = c + pos_disp + neg_disp;

    % The layers_stack will contain one label per image, adding up the layer
    % among all the images disparity compensated
    layers_stack = zeros(r, ext_cols, num_layers);

    for l = 1 : num_layers
        ith_layer                      = zeros(r, ext_cols);

        for i = 1 : num_images
            cur_disp  = disparity(i, l);
            cur_label = label(:,:,i);

            % Indices of the location corresponding to the current layer in the
            % current image
            ind_orig = (cur_label == l);

            % Indices corresponding to the location in the extended layer
            ind_dest  = false(r, ext_cols);
            first_col = neg_disp + 1 + cur_disp;
            last_col  = first_col + c - 1;
            ind_dest(:, first_col:last_col) = ind_orig;

            % Copy 
            ith_layer(ind_dest) = cur_label(ind_orig) / l;
        end

        layers_stack(:,:,l) = ith_layer;    
    end
    
    % ====> Data is available here in "layer_stack" variable
    
    % Tree parameters
    max_distances = 0;
    max_offsets   = 1;
    
    % Labels parameters
    [height, width, num_layer] = size(layers_stack);
    approx_layers              = zeros(height, width, num_layers-1);
    
    label_contours = cell(num_layers-1, 1);
    for l = 1 : num_layers-1
        label_contours{l} = estimate_contours_vertices(layers_stack(:,:,l), max_distances, max_offsets);
    end % end of l = 1 : num_layers

    bitstream = contour2bitstream(label_contours, width, height);
    
    label_contours_recov = bitstream2contour(bitstream);
    for l = 1 : num_layers-1
        approx_layers(:,:,l) = contour2shape(label_contours_recov{l}, width, height);
    end
    
    distortion = get_distortion(layers_stack(:,:,1:num_layers-1), approx_layers);
    
    disp(['File : ' file_name])
    disp(['Bits                      : ' num2str(length(bitstream))])
    disp(['Obtained distortion       : ' num2str(distortion)])
    disp('')
end

