function [valid selected_slice euclidean pearson all_features tot_time] = getHippoCampus(image_name, show_output, things_needed_for_ranking, my_hippocampus_output_folder, my_source_folder, my_segmented_output_folder, my_slice)

% % % parameters:
% image_name: name of the image the toolbox will extract the
% hippocampus. It must be placed in my_source_folder (if set),
% otherwise in ../dataset. The name is considered without extension, the
% file must be a .jpg
% show_output: -1 no output, 0 no output but save temporary images and files, 
% 1 output text, 2 output text and most relevant images, 3 all images
% things_needed_for_ranking the name of the mat file containing the
% information used for ranking. If it is missing, no ranking will be
% performed
% my_hippocampus_output_folder: save located hippocampus in this folder
% my_source_folder: read input images from this folder
% my_segmented_output_folder: write output images in this folder
% my_slice: look only around this slice, not in all of them


% % % return values:
% result of the validation (1=valid, 0=an error occurred)
% best slice
% euclidean distance (with reference to the prototype contained in
% things_needed_for_ranking)
% pearson distance
% all features (-1 if not valid)
% total time employed

close all

%% GLOBAL VARIABLES, CONSTANTS ...
%FOLDERS OF THE IMAGES
if exist('my_source_folder', 'var') && ~strcmp(my_source_folder,'')
    SOURCE_FOLDER = my_source_folder;
else
    SOURCE_FOLDER = '../dataset';
end
OUTPUT_FOLDER = 'images/hippocampi';
REF_FOLDER = 'images/reference';
REF_CONT_FOLDER = 'images/reference_cont';
ELAB_FOLDER = 'images/elab';
ELAB_CONT_FOLDER = 'images/elab_cont';
if exist('my_segmented_output_folder', 'var') && ~strcmp(my_segmented_output_folder,'')
    SEGMENTED_HOME_FOLDER = my_segmented_output_folder;
else
    SEGMENTED_HOME_FOLDER = 'images/segmented';
end
IMAGE_SIZE = [255 517]; %it is the size of the biggest reference image
if exist('my_slice', 'var')
    REFERENCES = [my_slice - 2 ; my_slice ; my_slice + 2];
else
    REFERENCES = [ 30; 32; 34; 36; 38; 40; 42 ];
end
%these rois has been selected manually (xmin ymin xmax ymax)
% HIPPO_REGION.('reference24') = [166 58 245 112]; % 24
% HIPPO_REGION.('reference26') = [155 49 240 105]; % 26
HIPPO_REGION.('reference28') = [195 50 273 110]; % 28
HIPPO_REGION.('reference30') = [200 54 272 105]; % 30
HIPPO_REGION.('reference32') = [206 53 272 105]; % 32
HIPPO_REGION.('reference34') = [220 45 299 104]; % 34
HIPPO_REGION.('reference36') = [215 53 278 109]; % 36
HIPPO_REGION.('reference38') = [220 61 289 110]; % 38
HIPPO_REGION.('reference40') = [212 66 271 110]; % 40
HIPPO_REGION.('reference42') = [215 76 278 118]; % 42
HIPPO_REGION.('reference44') = [224 68 289 114]; % 44
% HIPPO_REGION.('reference46') = [180 76 231 107]; % 46
% HIPPO_REGION.('reference50') = [200 75 251 106]; % 50

MY_OFFSET = [-5 -5 5 5]; %xmin, ymin, xmax, ymax: increase ROI area

% SET: if the best slice selection phase must take into account also local
% registration (slower but more accurate)
DO_LOCAL_REGISTRATION = false;

% ADD ALL FUNCTION PATHS
addpaths

%% REGISTRATION OPTIMIZATION PARAMETERS

% boundaries of registration (offset x, offset y, rotation, resize x, resize y)
lbg=[-50  -50  -0.314  0.75  0.75];
ubg=[50  50  0.314  1.5  1.5];

%parameters of Levenberg-Marquardt for global registration
if show_output <= 0
    display_lm = 'off';
else
    display_lm = 'final';
end
itermax_lm = 2000;

%% PRELIMINAR PHASES OF BEST REFERENCE SLICE SELECTION
tic
source_image_name = sprintf('%s/%s.jpg',SOURCE_FOLDER,image_name);
if ~exist(source_image_name,'file')
    error(strcat('Source file (',image_name,'.jpg) is not in the correct folder (', SOURCE_FOLDER ,'). Elaboration stopped.'));
end
if show_output > 0
    fprintf(1,'Reading Input Image...\n');
end

%read image
Iorig = cvReadImage(source_image_name,0);

if show_output > 0
    fprintf(1,'Start Contour Extraction...\n');
end
if show_output > -1
    image_filename = sprintf('%s/%s_cont',ELAB_CONT_FOLDER,image_name); %without extension
    entire_image_filename = sprintf('%s/%s',ELAB_FOLDER,image_name); %without extension
    command = ['Extract_External_Contour.exe ', source_image_name, ' ', image_filename, ' ', entire_image_filename];
    image_filename = strcat(image_filename, '.png');
    entire_image_filename = strcat(entire_image_filename, '.png');
else
    image_filename = 'out_c.png';
    entire_image_filename = 'out.png';
    command = ['Extract_External_Contour.exe ', source_image_name];%, ' ', image_filename, ' ', entire_image_filename];
end
system(command); %execute external command

res_img = im2double(imread(entire_image_filename)); %it is the original image resized and in graylevels

%output handling
if show_output > -1
    timestamp = datestr(now, 'yymmdd_HHMMSS');
    out_dir_name = sprintf('%s/%s_%s',OUTPUT_FOLDER, image_name,timestamp);
    mkdir(out_dir_name);
    out_file_name = sprintf('%s/%s_%s.txt',out_dir_name,image_name,timestamp);
    f = fopen(out_file_name,'w');
end

% at this point in image_filename there is the name of the image to register

I1=im2double(imread(image_filename));
% convert to gray scale if needed
if (size(I1,3)>1)
    I1 = rgb2gray(I1);
end;

%if the target image is larger than the largest reference, change it.
orig_width_target = size(I1, 2);
if size(I1, 2) > IMAGE_SIZE(2)
    IMAGE_SIZE(2) = size(I1, 2);
else %fix width
    I1 = fix_image_width(I1,IMAGE_SIZE(2),0);
end

%% GLOBAL REGISTRATION
siz_ref = size(REFERENCES);
fitnesses = zeros(2,siz_ref(1)); %this will keep the fitnesses functions (global and local)
sum_variance = 0.0; %avg variance of selected rois
orig_width_ref = zeros(siz_ref(1),1);
P = zeros(siz_ref(1),5); %registration values
scale = [1 1 1 1 1];
if show_output > 0
    fprintf(1,'Start Global Registration...\n');
end
ORIG_HIPPO_REGION = zeros(siz_ref(1),4,3);
CROP_HIPPO_REGION = zeros(siz_ref(1),5,3);
T_MATRIX = zeros(siz_ref(1),3,3);

for ref = 1:siz_ref(1)  %for each reference atlas
    actual_reference = sprintf('reference%i',REFERENCES(ref));
    if show_output > 0
        fprintf(1,'%s\n',actual_reference);
    end

    reference_filename = sprintf('%s/%s_cont.png',REF_CONT_FOLDER,actual_reference);
    I2=im2double(imread(reference_filename));

    if (size(I2,3)>1)
        I2 = rgb2gray(I2);
    end;

    orig_width_ref(ref) = size(I2, 2);
    if size(I2, 2) < IMAGE_SIZE(2)
        I2 = fix_image_width(I2, IMAGE_SIZE(2),0);
    end
    
    
    doAffine(0,I1, I2, lbg, ubg, 'sd'); %initialize by passing the two images, the boundaries and the distance method

    %optimization with levenberg-marquardt
    [xmin xxx tfxmin]=lsqnonlin(@(x)affine_mask(x,scale,I1,I2,'sd',lbg,ubg),[0.5 0.5 0.5 0.5 0.5],[],[],optimset('Display',display_lm,'TolX',0,'MaxIter',30000,'TolFun',0,'MaxFunEvals', itermax_lm, 'Algorithm','levenberg-marquardt'));
    xmin = xmin.*(ubg-lbg) + lbg;
    
    fitnesses(1,ref) = tfxmin;
    P(ref,:)=xmin;
    M=make_transformation_matrix(xmin(1:2),xmin(3),xmin(4:5));
    T_MATRIX(ref,:,:) = M;
    Minv=make_transformation_matrix([xmin(2) xmin(1)],-xmin(3),[xmin(5) xmin(4)]);

    Icor=affine_transform(I1,M,3); % 3 stands for cubic interpolation
    
    ref_hippo = HIPPO_REGION.(actual_reference); %load ROI coordinates
    
    %add offset to enlarge the roi
    ref_hippo = ref_hippo + MY_OFFSET;
    ref_hippo = [ref_hippo(1) ref_hippo(2) 1; ref_hippo(1) ref_hippo(4) 1; ref_hippo(3) ref_hippo(2) 1; ref_hippo(3) ref_hippo(4) 1]; %from coordinates to points
    ref_hippo(:,1) = ref_hippo(:,1) + (IMAGE_SIZE(2) - orig_width_ref(ref)) / 2.0;  %move if you added bands on side
    ref_to_draw = [ref_hippo(1,:);ref_hippo(2,:);ref_hippo(4,:);ref_hippo(3,:);ref_hippo(1,:)]; %to draw a rectangle on images
    change_coords = repmat([IMAGE_SIZE(2)./2 , IMAGE_SIZE(1)./2 , 0],4,1);
    ref_hippo = ref_hippo - change_coords;
    orig_hippo = (Minv*ref_hippo')'; %back transformation
    orig_hippo = orig_hippo + change_coords;
    orig_to_draw = [orig_hippo(1,:);orig_hippo(2,:);orig_hippo(4,:);orig_hippo(3,:);orig_hippo(1,:)];
    ORIG_HIPPO_REGION(ref,:,:) = orig_hippo;  %remember original and trasformed roi for the next cycle
    CROP_HIPPO_REGION(ref,:,:) = orig_to_draw;
    if show_output > 2 %show some images if desired
        entire_ref_filename = sprintf('%s/%s.png',REF_FOLDER,actual_reference);
        I1show = res_img; %load grayscale images and adjust the width
        I2show = im2double(rgb2gray(imread(entire_ref_filename)));
        I1show = fix_image_width(I1show, IMAGE_SIZE(2),1);
        I2show = fix_image_width(I2show, IMAGE_SIZE(2),1);
        Icorshow=affine_transform(I1show,M,3);
        figure_title = sprintf('Contour Not Registered, %d',REFERENCES(ref));
        figure('Name',figure_title,'NumberTitle','off'); imshow(I2)
        hold on
        h = imshow(I1);
        set(h, 'AlphaData', 0.5)

        figure_title = sprintf('Contour Registered, %d',REFERENCES(ref));
        figure('Name',figure_title,'NumberTitle','off'); imshow(I2)
        hold on
        h = imshow(Icor);
        set(h, 'AlphaData', 0.5)

        figure_title = sprintf('Image Not Registered, %d',REFERENCES(ref));
        figure('Name',figure_title,'NumberTitle','off'); imshow(I2show)
        hold on
        line(orig_to_draw(:,1),orig_to_draw(:,2));
        h = imshow(I1show);
        set(h, 'AlphaData', 0.5)

        figure_title = sprintf('Image Registered, %d',REFERENCES(ref));
        figure('Name',figure_title,'NumberTitle','off'); imshow(I2show)
        hold on
        line(ref_to_draw(:,1),ref_to_draw(:,2));
        h = imshow(Icorshow);
        set(h, 'AlphaData', 0.5)
    end         
    clear I2
end %for global registration


%% LOCAL REGISTRATION
if DO_LOCAL_REGISTRATION
    PERC_MIN_SAT = 0.7; %saturate background
    PERC_MAX_SAT = 0.95; %saturate possible parts darker than hippocampus
    WORST_GLOBAL_FITNESS_FACTOR = 1.3; %to discard slices with very bad global registration
    %don't do the local registration on images with a very low global fitness
    calculate_local = fitnesses(1,:)<(min(fitnesses(1,:))*WORST_GLOBAL_FITNESS_FACTOR);

    % pso parameters for local registration
    pso_options = get_local_pso_params(show_output);
    % boundaries for local registration
    lbl = [-14  -14  -0.314/2  0.9  0.9];
    ubl = [14  14  0.314/2  1.1  1.1];
    
    if show_output > 0
        fprintf(1,'Start Local Registration...\n');
    end
    for ref = 1:siz_ref(1)  %for all reference atlases
        if calculate_local(ref) %only for the ones with good global fitness
            actual_reference = sprintf('reference%i',REFERENCES(ref));
            if show_output > 0
                fprintf(1,'%s\n',actual_reference);
            end
            entire_ref_filename = sprintf('%s/%s.png',REF_FOLDER,actual_reference);
            ref_to_draw = reshape(CROP_HIPPO_REGION(ref,:,:),5,3);
            minx = min(ref_to_draw(:,1));
            miny = min(ref_to_draw(:,2));
            width = max(ref_to_draw(:,1))-minx;
            height = max(ref_to_draw(:,2))-miny;

            %load images and crop hippocampi
            I1show = res_img;%im2double(rgb2gray(imread(entire_image_filename)));
            I2show = im2double(rgb2gray(imread(entire_ref_filename)));
            I1show = fix_image_width(I1show, IMAGE_SIZE(2),1);
            I2show = fix_image_width(I2show, IMAGE_SIZE(2),1);
            M = reshape(T_MATRIX(ref,:,:),3,3);
            Icorshow=affine_transform(I1show,M,1);
            hip1 = imcrop(Icorshow,[minx miny width height]);
            hip2 = imcrop(I2show,[minx miny width height]);
            sum_variance = sum_variance + std(hip1(:));
            %here we can set a normalization-equalization of images
            hip1 = preprocess_small_hippo(hip1, PERC_MIN_SAT, PERC_MAX_SAT);
            hip2 = preprocess_small_hippo(hip2, PERC_MIN_SAT, PERC_MAX_SAT);

            doAffine(0,hip1, hip2, lbl, ubl, 'sd'); %initialize by passing the two images and the distance method
            [tfxminh, xminh] = pso_proposed(pso_options); %optimization
            xminh = xminh.*(ubl-lbl) + lbl;
            LOCAL_TRANSFORM.(actual_reference) = xminh;
            fitnesses(2,ref) = tfxminh;
        else
            fitnesses(2,ref) = 100;
            xminh = [0 0 0 0 0];
        end
        if show_output > -1
            fprintf(f,'%d\t%f\t%f\t[%f %f,%f,%f %f]\t[%f %f,%f,%f %f]\n',REFERENCES(ref), fitnesses(1,ref), fitnesses(2,ref),P(ref,1),P(ref,2),P(ref,3),P(ref,4),P(ref,5),xminh(1),xminh(2),xminh(3),xminh(4),xminh(5));
        end
    end %local registration
end
clear I1
%% BEST REFERENCE SLICE SELECTION

if DO_LOCAL_REGISTRATION %weighted mean of the two registrations
    MIN_TRUST = 0.3;
    MAX_TRUST = 0.7;
    MIN_STD = 0.5;
    MAX_STD = 1.5;
    NORM_FACTOR = 10;
    %calculate reliability measure and fitnesses
    %it's a linear relation
    sum_variance = (sum_variance/sum(calculate_local))*NORM_FACTOR;
    trust_roi = (sum_variance-MIN_STD)/((MAX_STD-MIN_STD)/(MAX_TRUST-MIN_TRUST)) + MIN_TRUST;
    if trust_roi < MIN_TRUST
        trust_roi = MIN_TRUST;
    elseif trust_roi > MAX_TRUST
        trust_roi = MAX_TRUST;
    end
    fitn = (1-trust_roi)*fitnesses(1,:)+trust_roi*fitnesses(2,:);
    if show_output > -1
        fprintf(f,'TRUST: %f\n',trust_roi);
        for ref = 1:siz_ref(1)
            fprintf(f,'%d: %f\n',REFERENCES(ref),fitn(ref));
        end
    end
else
    fitn = fitnesses(1,:);
end
[best_fitn best_ref] = min(fitn); %find best reference atlas
selected_slice=REFERENCES(best_ref);
best_orig_hippo = reshape(ORIG_HIPPO_REGION(best_ref,:,:),4,3);


% do operations on the best: 
% best_orig_hippo contains the coordinates of the best hippocampus. From
% here you can obtain the image of the hippocampus in the original image.

if show_output > 0
    fprintf(1,'Saving the hippocampus...\n');
end

%reproject to original image 
ratio = repmat([size(Iorig,2), size(Iorig,1), 1]./[orig_width_target,IMAGE_SIZE(1),1],5,1);
orig_to_draw = [best_orig_hippo(1,:);best_orig_hippo(2,:);best_orig_hippo(4,:);best_orig_hippo(3,:);best_orig_hippo(1,:)];
if show_output > 2
    figure('Name','Hippocampus','NumberTitle','Off');imshow(I1show);
    hold on
    line(orig_to_draw(:,1),orig_to_draw(:,2));     
end
best_orig_hippo(:,1) = best_orig_hippo(:,1) - (IMAGE_SIZE(2) - orig_width_target) / 2.0; %to remove the enlargement of image
big_to_draw = [best_orig_hippo(1,:);best_orig_hippo(2,:);best_orig_hippo(4,:);best_orig_hippo(3,:);best_orig_hippo(1,:)].*ratio;

minx = min(big_to_draw(:,1)); %because they can be rotated...
miny = min(big_to_draw(:,2));
width = max(big_to_draw(:,1))-minx;
height = max(big_to_draw(:,2))-miny;
hip_img = imcrop(Iorig,[minx miny width height]);
if show_output > -1
    hippo_filename = sprintf('%s/%s_hippo_ref%d_%s.png',out_dir_name,image_name, REFERENCES(best_ref),timestamp); %save extracted ROI
    imwrite(hip_img, hippo_filename,'png');
end

if show_output > 2
    figure('Name','Original Hippocampus','NumberTitle','Off'); imshow(Iorig);
    hold on
    line(big_to_draw(:,1),big_to_draw(:,2));
end

clear Iorig

if show_output > 0
    fprintf(1,'Best Reference Slice Selection complete.\n');
end

%% ROI ELABORATION
%in hip_img there is the extracted ROI
image_side = 800; %roi_size
hip_img = im2double(hip_img);
%preprocess_roi is a method based on histogram equalization and connected components
roi_img = preprocess_roi(hip_img);
if show_output > 2
    roi_filename = sprintf('%s/%s_roi_ref%d_%s.png',out_dir_name,image_name, REFERENCES(best_ref), timestamp);
    imwrite(roi_img, roi_filename,'png'); %save reduced and binarized ROI
end

%% PRELIMINAR PHASES OF HIPPOCAMPUS LOCALIZATION 
output_hippo_bound = 30; % when you find a model, how much space leave around in the saved image
ERROR_DG = 16.0;%threshold of acceptable error on DG model
ERROR_CA = 32.0;%threshold of acceptable error on CA model
N_MAX_TEST_MODEL = min([3, numel(REFERENCES)]); %don't even try with very bad models!
n_m = 1; %index of fitn_sort
if show_output > 0
    fprintf(1,'Hippocampus Localization starting...\n');
end
target = roi_img;
target([1 size(target,1)],:) = 0;%because i don't want the model to go to the bounds
target(:,[1 size(target,2)]) = 0;
[xxx fitn_sort] = sort(fitn); %fitn_sort will contain the order of the models to be used
model_used = REFERENCES(fitn_sort(n_m));

[movements_min] = get_models_std(model_used,1); %called just to set n_points
n_points_dg = size(movements_min,2)/2;
PENALTIES_DG = [27 0.055 1]; %penalty factors for the internal forces of DG NOTE: now the last element is not used

%% DIFFERENTIAL EVOLUTION PARAMETERS
S_struct = get_localization_de_params(show_output, n_points_dg);

model = zeros(n_points_dg,2); %it will contain the cartesian coordinates of the points
movements_cart = zeros(size(model)); %it will contain the cartesian coordinates of the shifts

%% DG LOCALIZATION
if show_output > 0
    fprintf(1,'Dentate Gyrus Localization.\n');
end

best_dg_error = 5*n_points_dg + 20*(n_points_dg-1) + 100; %maximum possible error
dg_error = best_dg_error;
while dg_error > ERROR_DG && n_m <= N_MAX_TEST_MODEL
    model_used = REFERENCES(fitn_sort(n_m));
    if show_output > 0
        fprintf(1,'Using model %d\n', model_used);
    end
    [movements_min movements_max movements_mean outer_shift] = get_models_std(model_used,1);
    outer_shift_cart = zeros(n_points_dg,2);
    [outer_shift_cart(:,1) outer_shift_cart(:,2)] = pol2cart(outer_shift(:,1),outer_shift(:,2));
    
    outer_shift_cart_op = ceil(reshape(outer_shift_cart',n_points_dg*2,1));
    de_find(0,target,outer_shift_cart_op,movements_min, movements_max, movements_mean, PENALTIES_DG);
    
    %optimization
    [xmin,tfxmin]=deopt('deopt_find_mask',S_struct);
    
    xmin = xmin.*(movements_max-movements_min)+movements_min;
    xmin = reshape(xmin,n_points_dg,2);
    [movements_cart(:,1) movements_cart(:,2)] = pol2cart(xmin(:,1),xmin(:,2));
    model(1,:) = ceil(movements_cart(1,:));
    for i=2:size(model,1) %obtain points from the shifts
        model(i,:) = model(i-1,:) + ceil(movements_cart(i,:));
    end
    %move points outside the image to the boundaries
    model(model<1)=1; model(model>image_side)=image_side;
    model_out = model + ceil(outer_shift_cart);
    model_out(model_out<1)=1; model_out(model_out>image_side)=image_side;

    dg_error = 5*n_points_dg + 20*(n_points_dg-1) + tfxmin;
    if show_output > 0
        fprintf(1,'DG Displacement: %f\n',dg_error);
    end
    if n_m == 1 %in this way i am sure that the model exists
        modeldg = model;
    end
    if dg_error < best_dg_error
        modeldg = model;
        modeldg_out = model_out;
        best_dg_error = dg_error;
    end
    n_m = n_m + 1;
    
end %while

if show_output > 2
    figure('Name','Hippocampus on ROI','NumberTitle','off');
    imshow(target);
    hold on
    line(modeldg(:,1),modeldg(:,2),'Color','r','LineWidth',4);
    line(modeldg_out(:,1),modeldg_out(:,2),'Color','g','LineWidth',2);
end

%% CA LOCALIZATION
if show_output > 0
    fprintf(1,'Ammon''s Horn Localization.\n');
end
model_used = REFERENCES(fitn_sort(1));
[movements_min] = get_models_std(model_used,2);
n_points = size(movements_min,2)/2;
PENALTIES_CA = [20 0.025 0.065];

 

best_ca_error = 5*n_points + 20*(n_points-1) + 100;
ca_error = best_ca_error;
n_m = 1; %index of model to try
while ca_error > ERROR_CA && n_m <= N_MAX_TEST_MODEL
    model_used = REFERENCES(fitn_sort(n_m));
    if show_output > 0
        fprintf(1,'Using model %d\n', model_used);
    end
    [movements_min movements_max movements_mean outer_shift] = get_models_std(model_used,2);
    n_points = size(movements_min,2)/2;
    %change differential evolution params for CA
    S_struct = get_localization_de_params(show_output, n_points);
    outer_shift_cart = zeros(n_points,2);
    [outer_shift_cart(:,1) outer_shift_cart(:,2)] = pol2cart(outer_shift(:,1),outer_shift(:,2));

    %change the possible starting position of CA model based on DG model results
    [movements_max(1) movements_max(n_points+1)] = cart2pol(modeldg(4,1),modeldg(4,2));%this part is above the other
     movements_max(n_points+1) =  movements_max(n_points+1)*1.1; %start position vector module
     movements_max(1) =  movements_max(1)*0.9; %it can not start from the very same point
     movements_min(1) =  max([movements_max(1)-0.314 movements_min(1)]);
     
    outer_shift_cart_op = ceil(reshape(outer_shift_cart',n_points*2,1));
    de_find(0,target,outer_shift_cart_op,movements_min, movements_max, movements_mean, PENALTIES_CA);%[27 0.055 1]);
    [xmin,tfxmin]=deopt('deopt_find_mask',S_struct);
    model = zeros(n_points,2);
    movements_cart = zeros(size(model));
    xmin = xmin.*(movements_max-movements_min)+movements_min;
    xmin = reshape(xmin, n_points,2);
    [movements_cart(:,1) movements_cart(:,2)] = pol2cart(xmin(:,1),xmin(:,2));
    model(1,:) = ceil(movements_cart(1,:));
    for i=2:size(model,1)
        model(i,:) = model(i-1,:) + ceil(movements_cart(i,:));
    end
    model(model<1)=1; model(model>image_side)=image_side;
    model_out = model + ceil(outer_shift_cart);
    model_out(model_out<1)=1; model_out(model_out>image_side)=image_side;

    ca_error = 5*n_points + 20*(n_points-1) + tfxmin;%pso_hippo_fitness([model; model_out], target, reshape(movements_mean,n_points-1,2),[0 0 0]); % here i want to see only how many points are outside their positions, i don't care about distortion
    if show_output > 0
        fprintf(1,'CA Displacement: %f\n',ca_error);
    end
    if n_m == 1 %in this way i am sure that the model exists
        modelca = model;
    end
    if ca_error < best_ca_error
        modelca = model;
        modelca_out = model_out;
        best_ca_error = ca_error;
    end
    n_m = n_m + 1;
    
end %while

if show_output > 2
    line(modelca(:,1),modelca(:,2),'Color','r','LineWidth',4);
    line(modelca_out(:,1),modelca_out(:,2),'Color','g','LineWidth',2);
end

%% OUTPUT LOCALIZATION
n_points = length(modelca);
if show_output > 0
    fprintf(1,'Saving results of localization...\n');
end
ratio = size(hip_img)./size(target);
ratio = [ratio(2) ratio(1)];
roi_points = ceil(repmat(ratio,n_points,1).*modelca);
roi_pointsdg = ceil(repmat(ratio,n_points_dg,1).*modeldg);
roi_points_out = ceil(repmat(ratio,n_points,1).*modelca_out);
roi_pointsdg_out = ceil(repmat(ratio,n_points_dg,1).*modeldg_out);

xmin = min([roi_pointsdg_out(:,1); roi_points_out(:,1); roi_pointsdg(:,1); roi_points(:,1)])-output_hippo_bound;
ymin = min([roi_pointsdg_out(:,2); roi_points_out(:,2); roi_pointsdg(:,2); roi_points(:,2)])-output_hippo_bound;
w = max([roi_pointsdg_out(:,1); roi_points_out(:,1); roi_pointsdg(:,1); roi_points(:,1)]) - xmin+output_hippo_bound*2;
h = max([roi_pointsdg_out(:,2); roi_points_out(:,2); roi_pointsdg(:,2); roi_points(:,2)]) - ymin+output_hippo_bound*2;
loc_hip = imcrop(hip_img, [xmin ymin w h]);
if show_output > -1
    loc_filename = sprintf('%s/%s_loc_ref%d_%s.png',out_dir_name,image_name, REFERENCES(best_ref), timestamp);
    imwrite(loc_hip,loc_filename,'png');
end

if show_output > 1
    figure('Name','Hippocampus Localization','NumberTitle','off');
    imshow(hip_img)
    hold on
    line(roi_points(:,1),roi_points(:,2),'Color','y','LineWidth',4);
    line(roi_pointsdg(:,1),roi_pointsdg(:,2),'Color','r','LineWidth',4);
    tmp_filename = sprintf('%s/%s_onmodel_ref%d_%s.png',out_dir_name,image_name, REFERENCES(best_ref), timestamp);
    saveas(gcf,tmp_filename,'png');
end

if exist('my_hippocampus_output_folder', 'var') && ~strcmp(my_hippocampus_output_folder,'') %the user wants to save the hippocampus
    roi_filename = sprintf('%s/%s.png',my_hippocampus_output_folder,image_name);
    imwrite(hip_img, roi_filename,'png'); %save reduced and binarized ROI    
end


%to write the points in a file
roi_points(:,1) = roi_points(:,1) - xmin;
roi_points(:,2) = roi_points(:,2) - ymin;
roi_pointsdg(:,1) = roi_pointsdg(:,1) - xmin;
roi_pointsdg(:,2) = roi_pointsdg(:,2) - ymin;
roi_points_out(:,1) = roi_points_out(:,1) - xmin;
roi_points_out(:,2) = roi_points_out(:,2) - ymin;
roi_pointsdg_out(:,1) = roi_pointsdg_out(:,1) - xmin;
roi_pointsdg_out(:,2) = roi_pointsdg_out(:,2) - ymin;

if show_output > -1
    fprintf(f,'DG model\n');
    for i = 1:n_points_dg
        fprintf(f,'%d %d\n',roi_pointsdg(i,1),roi_pointsdg(i,2));
    end
    fprintf(f,'CA model\n');
    for i = 1:n_points
        fprintf(f,'%d %d\n',roi_points(i,1),roi_points(i,2));
    end
    fprintf(f,'External DG model\n');
    for i = 1:n_points_dg
        fprintf(f,'%d %d\n',roi_pointsdg_out(i,1),roi_pointsdg_out(i,2));
    end
    fprintf(f,'External CA model\n');
    for i = 1:n_points
        fprintf(f,'%d %d\n',roi_points_out(i,1),roi_points_out(i,2));
    end
    fclose(f);
end

clear hip_img;

%% EXTRACTION OF THE PARTS
%loc_hip contains the image
%roi_points, roi_pointsdg (+out) contain the coordinates of the model on loc_hip
if show_output > 0
    fprintf(1,'Extracting Parts of the Hippocampus (CA1, CA3, DG)...\n');
end
[ca1im ca3im dgim] = simpleSegmentation(loc_hip, roi_pointsdg, roi_points, roi_pointsdg_out, roi_points_out);

if show_output > -1
    tmp_filename = sprintf('%s/%s_ca1_%s.png',out_dir_name,image_name, timestamp);
    imwrite(ca1im, tmp_filename, 'png');
    tmp_filename = sprintf('%s/%s_ca3_%s.png',out_dir_name,image_name, timestamp);
    imwrite(ca3im, tmp_filename, 'png');
    tmp_filename = sprintf('%s/%s_dg_%s.png',out_dir_name,image_name, timestamp);
    imwrite(dgim, tmp_filename, 'png');
end
if show_output > 1
    figure('Name','Segmented Parts','NumberTitle','off');
    subplot(1,3,1), imshow(ca1im);
    subplot(1,3,2), imshow(ca3im);
    subplot(1,3,3), imshow(dgim);
end

[dg_in dg_out_int dg_out_ext dg_out ca1_in ca1_out ca3_in ca3_out ...
    ca1_border_int ca1_border_ext ca3_border_int ca3_border_ext dg_border_int dg_border_ext] = ...
    get_segmented_images(ca1im, ca3im, dgim);
clear ca1im
clear ca3im
clear dgim

if any([size(ca1_in) size(ca3_in) size(dg_in)]==1) %surely there is an error
    valid = 0;
else
    valid = validateSegmentationResults(dg_in, dg_out_int, dg_out_ext, dg_out, ca1_in, ca1_out, ca3_in, ca3_out);
end
if valid
    if show_output > 0
        fprintf(1,'Saving Parts of the Hippocampus (CA1, CA3, DG)...\n');
    end
    store_segmented_images(dg_in, dg_out_int ,dg_out_ext ,ca1_in ,ca1_out ,ca3_in ,ca3_out ,...
    ca1_border_int ,ca1_border_ext ,ca3_border_int ,ca3_border_ext ,dg_border_int ,dg_border_ext, image_name, SEGMENTED_HOME_FOLDER);
    % VERY IMPORTANT:
    % things_needed_for_ranking must contains:
    % features: features to use (LOGICAL 1x220)
    % prototype: reference vector (DOUBLE 1xnumber_of_selected_features)
    % mymu, mysigma: normalization factors (DOUBLE 1x220)
    % to generate this file use GeneticFeatureSelectionToolbox
    if exist('things_needed_for_ranking', 'var') && ~strcmp(things_needed_for_ranking,'')
        if exist(things_needed_for_ranking, 'file')
            load(things_needed_for_ranking)
            [euclidean pearson all_features] = rankGeneFromImages(ca1_border_int ,ca1_border_ext, ca1_in ,ca1_out , ...
                                                     ca3_border_int ,ca3_border_ext ,ca3_in ,ca3_out , ...
                                                     dg_border_int ,dg_border_ext, dg_in, dg_out_int ,dg_out_ext , ...
                                                     selected_features, prototype, mymu, mysigma);
        else
            if show_output > -1
                fprintf('WARNING!!!! %s not found!!!', things_needed_for_ranking);
            end
            euclidean = -1;
            pearson = -1;
            all_features = -1;    
        end
    else
        euclidean = -1;
        pearson = -1;
        all_features = -1;
    end
    clear features
    clear prototype
    clear mymu
    clear mysigma

else
    euclidean = -1;
    pearson = -1;
    all_features = -1;
    if show_output > 0
        fprintf(1,'Results validation not passed. Try again!\n');
    end
end
clear ca1_border_int
clear ca1_border_ext
clear ca1_in
clear ca1_out
clear ca3_border_int
clear ca3_border_ext
clear ca3_in
clear ca3_out
clear dg_border_int
clear dg_border_ext
clear dg_in
clear dg_out_int
clear dg_out_ext
tot_time = toc;

end