Rubio_etal_2021/Code/Scripts/MAIN_Estim_Long_Axon_por_Esferas_EJEMPLOS.m

clear all;

% %%
clc


% Leemos el archivo de Matlab conteniendo las neuronas
Path_Data_1 = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Neuronas/Selección de NeuroMorpho 23-2/Específicas - tipo 1/';
Path_Data_File_1 = '/Neuronas Especificas - Tipo 1.txt';

% Path_Data_1 = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Neuronas/Selección de NeuroMorpho 23-2/Multiespecíficas - tipo 2/';
% Path_Data_File_1 = '/Neuronas Multiespecificas - Tipo 2.txt';

% Path_Data_1 = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Neuronas/Selección de NeuroMorpho 23-2/Inespecíficas - tipo 3/';
% Path_Data_File_1 = '/Neuronas Inespecificas - Tipo 3.txt';

% Path_Data_1 = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Neuronas/Selección de NeuroMorpho 23-2/Locales- tipo 4/';
% Path_Data_File_1 = '/Neuronas Locales - Tipo 4.txt';

% Read the names of files with raw data:

fileID = fopen([Path_Data_1, Path_Data_File_1]);    % Spanish Data by Carlos
%
names_1 = textscan(fileID,'%s', 'delimiter', '\n', 'whitespace', '');
fclose(fileID);
names_1 = names_1{1};
%
iii = 24;
names_1 = {names_1{iii}};

h_Cortes = 1;   % Analizamos un corte cada 'h_Cortes' cortes de la muestra
Grosor_Corte = 100; % En micras

% El 'frame' es la sampling box en la que se introduce la esfera; esa caja luego se introduce en la celda de la rejilla en que se ha dividido el corte.
sampling_box_height = Grosor_Corte;   % en micras
lado_x_Frame = Grosor_Corte;   % en micras
lado_y_Frame = Grosor_Corte;   % en micras

% Diams_Sonda = [10:5:50];    % Vector con los diámetros de la sonda a considerar. En micras
% Step_Lengths = [70:10:150];   % Vector con los step lengths a considerar. En micras
Diams_Sonda = 100;    % Vector con los diámetros de la sonda a considerar. En micras
Step_Lengths = 100;   % Vector con los step lengths a considerar. En micras


MATRIZ_Q = zeros( length( Step_Lengths ), length( Diams_Sonda ), length(names_1) );   % Filas para los step_lengths y columnas para los diámetros de la sonda
MATRIZ_ERROR_LENGTH = zeros( length( Step_Lengths ), length( Diams_Sonda ), length(names_1) );   % Filas para los step_lengths y columnas para los diámetros de la sonda
AXON_REAL_LENGTH = zeros(1, length(names_1));
ESTIMATED_AXON_LENGTH = zeros( length( Step_Lengths ), length( Diams_Sonda ), length(names_1) );

for rr = 1:length(names_1)
% for rr = [101, 102, 106, 110]
    tic
    % Leemos los archivos
    fname=fullfile(Path_Data_1, [names_1{rr}]);
    fname_new = [fname(1:end-4), '.mat'];
    load( fname_new );
    
    [Matriz_Q, Matriz_Error_Length, Axon_Real_Length, Estimated_Axon_Length] = Sgript_Estim_Long_Axon_Estereo ( AXON_Cell, AXON, Diams_Sonda, Step_Lengths, Grosor_Corte, h_Cortes, lado_x_Frame, lado_y_Frame, sampling_box_height );
    
    MATRIZ_Q( :, :, rr ) = Matriz_Q;
    MATRIZ_ERROR_LENGTH( :, :, rr ) = Matriz_Error_Length;
    AXON_REAL_LENGTH(1, rr) = Axon_Real_Length;
    ESTIMATED_AXON_LENGTH(:, :, rr) = Estimated_Axon_Length;
    clc
    toc
    disp(rr)
    
end

% %% Salvamos los resultados
% 
% clc
% 
% Path_Data_Save = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Resultados/';
% fname_Save=fullfile(Path_Data_Save, ['Result_Estim_Long_Axon_por_Estereo_h_Cortes=', num2str(h_Cortes),'_', Path_Data_File_1(2:end-4)]);
% save( fname_Save );

%%



%% Exploramos la distribución de los errores e intersecciones en diferentes corrdenadas [e, d]

ind_e = 1;
ind_d = 2;

num_bins = 15;

Distrib_Error_Length = MATRIZ_ERROR_LENGTH(ind_e, ind_d, :);
Distrib_Q = MATRIZ_Q(ind_e, ind_d, :);

figure('color', 'w', 'position', [50, 200, 1100, 400]); 

subplot(1, 2, 1);
histogram( Distrib_Error_Length, num_bins, 'facecolor', 'k', 'normalization', 'probability' );
xlabel('Error Length (%)', 'fontsize', 16);
% ylim([0, 0.2]);
% xlim([-25, 25]);
title(['Steph Length = ', num2str(Step_Lengths(ind_e)), '.   Probe Diameter = ', num2str(Diams_Sonda(ind_d))]);

subplot(1, 2, 2);
histogram( Distrib_Q, num_bins, 'facecolor', 'b', 'normalization', 'probability' );
xlabel('Intersections', 'fontsize', 16);
% ylim([0, 0.2]);
% xlim([-25, 25]);
title(['Steph Length = ', num2str(Step_Lengths(ind_e)), '.   Probe Diameter = ', num2str(Diams_Sonda(ind_d))]);

%% Ploteamos las medias y std de los errores

clc

Saturacion = 100;   % A partir de este error (%) saturamos la matriz
Matriz_Error_Length_MEAN_SATURADA = mean(abs(MATRIZ_ERROR_LENGTH), 3);
Matriz_Error_Length_MEAN_SATURADA( Matriz_Error_Length_MEAN_SATURADA >= Saturacion ) = Saturacion;

Matriz_Error_Length_STD_SATURADA = std(abs(MATRIZ_ERROR_LENGTH), 0, 3);
Matriz_Error_Length_STD_SATURADA( Matriz_Error_Length_STD_SATURADA >= Saturacion ) = Saturacion;

figure('color', 'w', 'position', [50, 200, 1100, 400]); 

suptitle(Path_Data_File_1(2:end-4));

subplot(1, 2, 1);
imagesc(Matriz_Error_Length_MEAN_SATURADA);
xticklabels({ '10', '15', '20', '25', '30', '35', '40', '45', '50' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Probe Diameter', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Mean Error Length (%)', 'fontsize', 15);
colorbar

subplot(1, 2, 2);
imagesc(Matriz_Error_Length_STD_SATURADA);
xticklabels({ '10', '15', '20', '25', '30', '35', '40', '45', '50' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Probe Diameter', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('STD Error Length (%)', 'fontsize', 15);
colorbar


%%

%% Ploteamos

rr = 46;
fname=fullfile(Path_Data_1, [names_1{rr}]);
fname_new = [fname(1:end-4), '.mat'];
load( fname_new );
%%
clc
corte = 6;
colores = {'r', 'b', 'm', 'k', 'g'};

figure;
suptitle(names_1{rr});

hold on;
for i = 1:corte
Puntos_Axon = AXON_Cell{1, i};
    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-r', 'linewidth', 1.5 );
end

for i = corte+1:length(AXON_Cell)

    Puntos_Axon = AXON_Cell{1, i};
    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-k', 'linewidth', 1.5 );
    xlabel('X');
    ylabel('Y');    
    zlabel('Z');
    % zlim([150 200]);
    % view([-96, -90]);
        axis equal
end

%%
clc

h_Cortes = 1;   % Analizamos un corte cada 'h_Cortes' cortes de la muestra
Grosor_Corte = 50; % En micras

% El 'frame' es la sampling box en la que se introduce la esfera; esa caja luego se introduce en la celda de la rejilla en que se ha dividido el corte.
sampling_box_height = Grosor_Corte;   % en micras
lado_x_Frame = Grosor_Corte;   % en micras
lado_y_Frame = Grosor_Corte;   % en micras

Diams_Sonda = 50;    % Vector con los diámetros de la sonda a considerar. En micras
Step_Lengths = 150;   % Vector con los step lengths a considerar. En micras

AXON_Cell_0 = AXON_Cell(1, 7:end);
AXON_0 = AXON(88:end, :);
% for rr = [101, 102, 106, 110]
    tic
    % Leemos los archivos
    fname=fullfile(Path_Data_1, [names_1{rr}]);
    fname_new = [fname(1:end-4), '.mat'];
    load( fname_new );
    
    [Matriz_Q, Matriz_Error_Length, Axon_Real_Length, Estimated_Axon_Length, Filas_Centros, Columnas_Centros, Z_Centros] = Sgript_Estim_Long_Axon_Estereo ( AXON_Cell_0, AXON_0, Diams_Sonda, Step_Lengths, Grosor_Corte, h_Cortes, lado_x_Frame, lado_y_Frame, sampling_box_height );
    
    MATRIZ_Q( :, :, rr ) = Matriz_Q;
    MATRIZ_ERROR_LENGTH( :, :, rr ) = Matriz_Error_Length;
    AXON_REAL_LENGTH(1, rr) = Axon_Real_Length;
    clc
    toc







%% 
rr = 46;

Saturacion = 100;   % A partir de este error (%) saturamos la matriz
Matriz_Error_Length_SATURADA = abs(MATRIZ_ERROR_LENGTH(:, :, rr));
Matriz_Error_Length_SATURADA( Matriz_Error_Length_SATURADA >= Saturacion ) = Saturacion;

Matriz_Q_Dummy = abs(MATRIZ_Q(:, :, rr));

figure('color', 'w', 'position', [50, 200, 1100, 400]); 
suptitle(names_1{rr});

subplot(1, 2, 1);
imagesc(Matriz_Error_Length_SATURADA);
xticklabels({ '10', '15', '20', '25', '30', '35', '40', '45', '50' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Probe Diameter', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Error Length (%)', 'fontsize', 15);
colorbar

subplot(1, 2, 2);
imagesc(Matriz_Q_Dummy);
xticklabels({ '10', '15', '20', '25', '30', '35', '40', '45', '50' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Probe Diameter', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Number of Intersections', 'fontsize', 15);
colorbar


%% Ploteamos el axón con las esferas

clc

[x,y,z] = sphere;
rad_sonda = Diams_Sonda/2;

figure;
hold on;
for i = 1:length(AXON_Cell)
    
    Puntos_Axon = AXON_Cell{1, i};
    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-k', 'markersize', 10, 'linewidth', 1.5 );
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
    % zlim([150 200]);
    % view([-96, -90]);
        axis equal
end
% Esferas
for i_x = Columnas_Centros   % Columna de la rejilla
    for i_y = Filas_Centros   % Fila de la rejilla
        for i_z = 1:h_Cortes:length(Z_Centros)   % Corte
            surf( rad_sonda*x + i_x, rad_sonda*y + i_y, rad_sonda*z + Z_Centros(i_z),'FaceAlpha',0.9,'EdgeColor','none' );
        end
    end
end























%% Para chequear
clc

for j=ind_CORTE
% j=437
        [row_a, col_a] = find(~cellfun(@isempty, dummy_CORTE{1, j}));   % para cada rama estos son los índices de las celdas de ese corte en el que aparecen dichas ramas
% disp([row_a, col_a]);
%
k = 2;

if length(row_a)>1
dummy_a = dummy_CORTE{1, j}{row_a(k), col_a(k)};    % Coordenadas de los puntos de la rama cuya intersección estamos estudiando
            % Obtenemos las distancias al centro de la sonda esférica:
            dist_centro_sonda = sqrt( ( dummy_a(:, 1) - Columnas_Centros(col_a(k)) ).^2 + ( dummy_a(:, 2) - Filas_Centros(row_a(k)) ).^2 +...
                                  ( dummy_a(:, 3) - Z_Centros(i) ).^2 );
                              
            num_inters = sum( abs( diff( dist_centro_sonda <= rad_sonda ) ) );
            if num_inters > 3
            disp(j)
            end
end
end
%% Comprobamos que está bien, ploteando rama a rama:

[x,y,z] = sphere;

% clc
i_x = col_a(k);   % Columna de la rejilla
i_y = row_a(k);   % Fila de la rejilla
i_z = i;   % Corte

colores = {'r', 'b', 'm', 'k', 'g'};

figure('color', 'w', 'position', [250, 200, 700, 600]);

hold on
% for iii = 1:length(Cortes_XY_Axon{i_z, 1})
for iii = ind_CORTE
    
    if ~isempty( Cortes_XY_Axon{i_z, 1}{1, iii} )
        
        dummy = Cortes_XY_Axon{i_z, 1}{1, iii};
        for i_fila = 1:size(dummy, 1)
            
            for i_columna = 1:size(dummy, 2)
                if ~isempty(dummy{i_fila, i_columna})
                    Puntos_Axon = dummy{i_fila, i_columna};
                    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-o', 'color', colores{mod(i, 5)+1}, 'markersize', 7, 'linewidth', 1.5 );
                end
            end
            
        end
        
    end
    axis equal
    xlim([Columnas_Rejilla(i_x), Columnas_Rejilla(i_x) + Rejilla_X]);
    ylim([Filas_Rejilla(i_y), Filas_Rejilla(i_y) + Rejilla_Y]);
    zlim([Cortes(i_z), Cortes(i_z) + Grosor_Corte]);
    view([32, 20]);
end
surf( rad_sonda*x + Columnas_Centros(col_a(k)), rad_sonda*y + Filas_Centros(row_a(k)), rad_sonda*z + Z_Centros(i) );
alpha 0.7
% surf( rad_sonda*x + Columnas_Centros(i_x), rad_sonda*y + Filas_Centros(i_y), rad_sonda*z + Z_Centros(i_z) );

%% Comprobamos que está bien, ploteando ramas enteras:

[x,y,z] = sphere;

% clc
i_x = 2;   % Columna de la rejilla
i_y = 1;   % Fila de la rejilla
i_z = 2;   % Corte


    ind_CORTE = find(~cellfun(@isempty, Cortes_XY_Axon{i_z, 1}));
    
colores = {'r', 'b', 'm', 'k', 'g'};

figure('color', 'w', 'position', [250, 200, 700, 600]);

hold on
% for iii = 1:length(Cortes_XY_Axon{i_z, 1})
for iii = ind_CORTE
    
    if ~isempty( Cortes_XY_Axon{i_z, 1}{1, iii} )
        
        dummy = Cortes_XY_Axon{i_z, 1}{1, iii};
        for i_fila = 1:size(dummy, 1)
            
            for i_columna = 1:size(dummy, 2)
                if ~isempty(dummy{i_fila, i_columna})
                    Puntos_Axon = dummy{i_fila, i_columna};
                    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-', 'color', colores{mod(i, 5)+1}, 'markersize', 7, 'linewidth', 1.5 );
                end
            end
            
        end
        
    end
    axis equal
    view([32, 20]);
end
surf( rad_sonda*x + Columnas_Centros(i_x), rad_sonda*y + Filas_Centros(i_y), rad_sonda*z + Z_Centros(i_z) );
alpha 0.7
xlabel('X');
ylabel('Y');
zlabel('Z');
%     xlim([Columnas_Rejilla(i_x), Columnas_Rejilla(i_x) + Rejilla_X]);
%     ylim([Filas_Rejilla(i_y), Filas_Rejilla(i_y) + Rejilla_Y]);
%     zlim([Cortes(i_z), Cortes(i_z) + Grosor_Corte]);

% % Para la neurona 5
% xlim([-100, 100]);
% ylim([-200, 0]);
% zlim([-80, 60]);

% % Para la neurona 4
% xlim([-600, 400]);
% ylim([-600, 800]);
% zlim([-200, 700]);

% Para la neurona 3
xlim([-250, 100]);
ylim([-300, 200]);
zlim([-100, 80]);






%%
figure
hold on
surf(peaks(30))
alpha 0.5
plot3(10,10,10,'r*')
hold off













%%