%%%%       Bakalarska práce - Segmentace optickeho disku v obrazových datech sitnice
%%%%        Autor prace : Radek Juracek
%%%%        Vedouci prace : Ing. Jan Odstrcilik, Ph.D.
%%%%        Datum odevzdani : 27.5.2022
%%%%        Ustav : Ustav biomedicinského inženyrstvi
%%%%        Fakulta : Fakulta elektrotechniky a kominikacnich technologii
%%%%        Vysoke uceni technicke v Brne

function [ final_segmentation_mask ] = BatAlgoMethod(experiment_metadata, method_params, image_id)
    %%%% BatAlgoMethod
    %%%% Function segments image via morfological operations and bat algorithm.    
    %%%%
    %%%% :param experiment_metadata: struct experiment metadata
    %%%% :param params: vector of parametrs
    %%%% :param ind: index of segmented image
    % Check method params
    if method_params(1)<method_params(2) | method_params(3)<method_params(4)
        project_paths = experiment_metadata.project_paths;
        dataset_mask = experiment_metadata.data_metadata.dataset_mask;
        image_names = experiment_metadata.data_metadata.image_names;
        bat_config = experiment_metadata.ba_config;

        % Read preproceseed image
        image_filepath = append(project_paths.preprocessed_dir, image_names{image_id}, '_preprocessed_image.png');
        preprocessed_image = imread(image_filepath);

        % Get each optimized parameter from params vector
        hist_lower_bound = method_params(1);
        hist_upper_bound = method_params(2);
        low_morf_bound = method_params(3);
        high_morf_bound = method_params(4);
        morf_step = method_params(5);

        % Make vector of morfologi radiuses
        morph_disc_radiuses = [low_morf_bound:morf_step:high_morf_bound];

        % Initial best values for loop
        best_objective_value_found = 0;
        best_threshold = 127;
        best_morph_image = preprocessed_image;    

        % Starting image
        morphed_image = preprocessed_image;

        % Restrict image histogram low and high values
        morphed_image(dataset_mask == true & morphed_image < hist_lower_bound) = hist_lower_bound;
        morphed_image(dataset_mask == true & morphed_image > hist_upper_bound) = hist_upper_bound; 

        % Start iteration of opening and closing operation
        for radius = morph_disc_radiuses

            % Apply morfological opening on closed image
            morphed_image = imopen(morphed_image, strel('disk', radius, 0));    
            % Apply morfological closing on opened image
            morphed_image = imclose(morphed_image, strel('disk', radius, 0));

            % Compute image score and treshold via Bat Algorithm
            [ solution_threshold, solution_objective_value ] = BatAlgoOptimalization(morphed_image, dataset_mask, bat_config);
            % CHECK NEW SOLUTION
            if solution_objective_value > best_objective_value_found % If is new echo higher that actual max echo
                best_objective_value_found = solution_objective_value;            % Update max echo
                best_morph_image = morphed_image;
                best_threshold = solution_threshold;         
            end
        end

        % Segment image by thresholding, filtering and filling
        segmentated_image = imbinarize(best_morph_image, best_threshold);
        segmentated_image = bwareafilt(segmentated_image, 1, 'largest'); 
        segmentated_image = imfill(segmentated_image, 'holes');
        final_segmentation_mask = approx_ellipse(segmentated_image);
    else
        final_segmentation_mask = zeros(experiment_metadata.data_metadata.scaled_image_size,'logical');
    end
end

function [ ellipse_image ] = approx_ellipse(image)
    %%%% Function get boundaries in binary input image. 
    %%%%  Function aprox elipse via leastsquare from boundaries positions. 
    %%%%  Than draw elipse in binary image.
    %%%% :param image: binary input image
    
    if sum(sum(image)) == 0
        ellipse_image = image;
        return
    end

    % Get positions of boundaries in binary image 
    boundaries = cell2mat(bwboundaries(image));
    % Keep only unique one
    boundaries = unique(boundaries, 'rows');
    try
        x = boundaries(:, 2);                                                                                            
        y = boundaries(:, 1);                        
    catch
        elipse_image = zeros(size(image), 'logical');
        disp('No boundaries found in binary image')
        return
    end
    % Aprox elipse via least square method
    ellipse = FitEllipse(x, y);  % Gets struct of elipse parametrs
    % Draw elipse in binary image
    try
        ellipse_image = draw_ellipse(ellipse, image);  % Draw binary elipse image
    catch
        ellipse_image = zeros(size(image),'logical');
    end
end

function [ elipse_image ] = draw_ellipse(el, image)
%%%%    Function DrawElipse generates binary image of ellipse
%%%%    :param el: ellipse  info
%%%%    :param image: binary input image
    xc = el.X0_in;                                    % Row of elipse centre
    yc = el.Y0_in;                                    % Col of elipse centre
    c = sqrt((el.long_axis/2)^2-(el.short_axis/2)^2); % Calculate lenght of C axis
    % % CALCULATE FOCALS POINT POSITION
    f1 = [ yc + sin(el.phi)*c  xc+cos(el.phi)*c];     % Position of focal point after rotation
    f2 = [ yc - sin(el.phi)*c  xc-cos(el.phi)*c];     % Position of focal point after rotation
    % DRAW ELIPSE
    [xx, yy] = meshgrid(1:size(image,2),1:size(image,1));   % Make meshgrit matrix
    elipse_image = logical(zeros(size(image)));             % Alocate image array background
    elipse_image(sqrt((xx - f1(2)).^2 + (yy - f1(1)).^2 ) + sqrt((xx - f2(2)).^2 + (yy - f2(1)).^2) < el.long_axis) = 1;% Draw elipse 
end

function [ best_position, best_echo ] = BatAlgoOptimalization(morphed_image, dataset_mask, bat_config);
%%%%    Function Bat Algorithm alocated bat population and run behaivor
%%%%    to find optimum between class variance of image.    
%%%%    :param images: uint8 tresholded image
%%%%    :param mask: binary mask of ROI in image
%%%%    :param conf: struct of bat algorithm configuration

    % Get population size and number of iterations
    pop_size = bat_config.population_size;
    num_of_iter = bat_config.num_of_iter;
    
    % Initialize loudness and constant
    loudness = ones([pop_size, 1]);
    loudness_beta_const = bat_config.loudness_beta;

    % Initialize pulse rate and constant
    pulse_rate = repmat(0.05, pop_size,1);
    pulse_rate_sigma_const = bat_config.pulse_rate_sigma;

    % Initialize frequency and constants
    bat_freqs = zeros(pop_size, 1);
    freq_min_const = bat_config.freq_min;
    freq_max_const = bat_config.freq_max;

    % Set search space boundaries;
	upper_space_bound = im2double(max(morphed_image(dataset_mask == 1)));
    lower_space_bound = im2double(min(morphed_image(dataset_mask == 1)));

    % Initial bat population echoes and positions
    bat_velocities = zeros(pop_size, 1);
    bat_echoes = zeros([pop_size, 1]);
    bat_positions = (upper_space_bound - lower_space_bound) .* rand(pop_size, 1) + lower_space_bound;

    % Compute initial population echoes
    for bat_id = 1:pop_size
        [ bat_echoes(bat_id) ] = compute_echo(morphed_image, bat_positions(bat_id), dataset_mask);
    end

    % Find best initial solution (position-threshold)
    [ best_echo, best_ind ] = max(bat_echoes);
    best_position = bat_positions(best_ind, :);

    %Start bat algorithm iterations loop
    for iteration = 1:num_of_iter                                                                                   
        for bat_id = 1:pop_size
   
            % Generate frequencies
            freq(bat_id, :) = freq_min_const + ((freq_max_const - freq_min_const) * rand(1,1));

            % Change bat velocity
            bat_velocities(bat_id, :) = bat_velocities(bat_id, :) + (bat_positions(bat_id) - best_position) .* bat_freqs(bat_id, :);

            % Compute new position
            new_bat_position = bat_positions(bat_id) + bat_velocities(bat_id, :);

            % Make random walk 
            if rand > pulse_rate(bat_id)                                                                                              
                new_bat_position = new_bat_position + mean(loudness) * randn(1,1) * 0.01;
            end

            % Apply space boundaries for new position
            new_bat_position(new_bat_position < lower_space_bound) = lower_space_bound;
            new_bat_position(new_bat_position > upper_space_bound) = upper_space_bound;

            % Evaluate echo of new position
            [ new_echo ] = compute_echo(morphed_image, new_bat_position, dataset_mask);

            % Update position if new echo is better and isnt too loud
            if (new_echo >= bat_echoes(bat_id)) & (rand < loudness(bat_id))
                bat_echoes(bat_id) = new_echo;
                bat_positions(bat_id) = new_bat_position;
                
                % Reduce loudness and pulse rate
                loudness(bat_id) = loudness_beta_const * loudness(bat_id);
                pulse_rate(bat_id) = pulse_rate(bat_id) * (1 - exp(-pulse_rate_sigma_const * iteration));
            end

            %Update best solution if in new echo is better 
            if new_echo > best_echo
                best_echo = new_echo;
                best_position = new_bat_position;
            end
        end
    end
end

function [echo_value] = compute_echo(input_image, threshold, mask);
%%%%  Objective function computes between classvariation of image
    bin = imbinarize(input_image, threshold);            % Binarize image with threhold(position of bat)
    bin = bwareafilt(bin, 1, 'largest');                 % Keep largest blob only
    bin = imfill(bin, 'holes');                          % Fill holes in blob
    if sum(bin(:)) > 0                                  % If filtred bin image contain blob
        class0 = input_image(bin  == 0 & mask == 1);    % Vector of intensity disk image in disk area
        class1 = input_image(bin  == 1 & mask == 1);    % Vector of intensity disk image in background
        class0_size = size(class0, 1);             % Size of intensity disk image vector - disk
        class1_size = size(class1, 1);             % Size of intensity disk image vector - background 
        u0 = mean(class0);                         % Mean of background area in disk image
        u1 = mean(class1);                             % Mean of disk area in disk image  
        w0 = class0_size / (class0_size + class1_size);  % Calculate class probability
        w1 = class1_size / (class0_size + class1_size);  % Calculate class probability          
        echo_value = w0 * w1 * ((u0 - u1)^2);                        % Calculate echo
    else
        echo_value = 0;
    end
end

