Rubio_etal_2021/Code/Scripts/MAIN_Estim_Long_Axon_por_Planos_Virtuales.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};
%%
clc

num_Repet = 1;   % Cada vez que se muestrea sale diferente porque los planos cambian. Se repite 5 veces apra obtener cierta estadística

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
Sampling_Box_X = 50;   % en micras
Sampling_Box_Y = 50;   % en micras

% Distancias_Entre_Planos = [3:3:30];    % Vector con las distancia entre planos virtuales en micras
Distancias_Entre_Planos = [1:1:10];    % Vector con las distancia entre planos virtuales en micras
Step_Lengths = [70:10:150];   % Vector con los step lengths a considerar. En micras
% Distancias_Entre_Planos = 5;    % Vector con los diámetros de la sonda a considerar. En micras
% Step_Lengths = 75;   % Vector con los step lengths a considerar. En micras

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

for kk = 1:num_Repet
    for rr = 1:length(names_1)
%     for rr = 1:230
%     for rr = 13
        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, Matriz_Estimated_Length] = Sgript_Estim_Long_Axon_Estereo_PLANOS ( AXON_Cell, AXON, Distancias_Entre_Planos, Step_Lengths, Grosor_Corte, h_Cortes, Sampling_Box_X, Sampling_Box_Y, sampling_box_height );
        
        MATRIZ_Q( :, :, rr, kk ) = Matriz_Q;
        MATRIZ_ERROR_LENGTH( :, :, rr, kk ) = Matriz_Error_Length;
        AXON_REAL_LENGTH(1, rr) = Axon_Real_Length;
        MATRIZ_ESTIMATED_AXON_LENGTH(:, :, rr, kk) = Matriz_Estimated_Length;
        clc
        toc
        
        % Salvamos los resultados
        Path_Data_Save = '/Users/pepo/Google Drive (1)/Cesar Mario Axones/Resultados/';
        fname_Save=fullfile(Path_Data_Save, ['Result_Estim_Long_Axon_por_Planos_h_Cortes=', num2str(h_Cortes),'_', Path_Data_File_1(2:end-4)]);
        save( fname_Save, 'h_Cortes', 'Grosor_Corte', 'sampling_box_height', 'Sampling_Box_X', 'Sampling_Box_Y', 'Distancias_Entre_Planos', 'Step_Lengths', 'MATRIZ_Q', 'MATRIZ_ERROR_LENGTH', 'AXON_REAL_LENGTH', 'MATRIZ_ESTIMATED_AXON_LENGTH' );
        
        disp(rr)
        
    end
    
    disp(kk)
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_Planos_h_Cortes=', num2str(h_Cortes),'_', Path_Data_File_1(2:end-4)]);
% save( fname_Save );


%%

%% 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(squeeze(MATRIZ_ERROR_LENGTH)), 3);
Matriz_Error_Length_MEAN_SATURADA( Matriz_Error_Length_MEAN_SATURADA >= Saturacion ) = Saturacion;

Matriz_Error_Length_STD_SATURADA = std(abs(squeeze(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]);
% figure('color', 'w', 'position', [50, 200, 560, 500]);

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

subplot(1, 2, 1);
% imagesc(Matriz_Error_Length_MEAN_SATURADA);
% yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
contourf(flipud(Matriz_Error_Length_MEAN_SATURADA), 10);
yticklabels(fliplr({ '70', '80', '90', '100', '110', '120', '130', '140', '150' }));

% xticklabels({ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10' });
% xticklabels({ '2', '4', '6', '8', '10', '12', '14', '16', '18', '20' });
xticklabels({ '3', '6', '9', '12', '15', '18', '21', '24', '27', '30' });
xlabel('Dist between Planes', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Virtual Plane Mean Error Length (%)', 'fontsize', 15);
% title('Length Estim Error (Virtual Planes)  (%)', 'fontsize', 15);
colorbar

subplot(1, 2, 2);
% imagesc(Matriz_Error_Length_STD_SATURADA);
% yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
contourf(flipud(Matriz_Error_Length_STD_SATURADA), 10);
yticklabels(fliplr({ '70', '80', '90', '100', '110', '120', '130', '140', '150' }));

% xticklabels({ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10' });
% xticklabels({ '2', '4', '6', '8', '10', '12', '14', '16', '18', '20' });
xticklabels({ '3', '6', '9', '12', '15', '18', '21', '24', '27', '30' });
xlabel('Dist between Planes', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Virtual Plane STD Error Length (%)', 'fontsize', 15);
colorbar

%% Ploteamos las medias y std de las intersecciones

clc

Matriz_Inters_MEAN = mean((MATRIZ_Q), 3);
Matriz_Inters_STD = std((MATRIZ_Q), 0, 3);

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

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

subplot(1, 2, 1);
imagesc(Matriz_Inters_MEAN);
xticklabels({ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10' });
% xticklabels({ '2', '4', '6', '8', '10', '12', '14', '16', '18', '20' });
% xticklabels({ '3', '6', '9', '12', '15', '18', '21', '24', '27', '30' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Dist between Planes', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Mean Intersections (%)', 'fontsize', 15);
colorbar

subplot(1, 2, 2);
imagesc(Matriz_Inters_STD);
xticklabels({ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10' });
% xticklabels({ '2', '4', '6', '8', '10', '12', '14', '16', '18', '20' });
% xticklabels({ '3', '6', '9', '12', '15', '18', '21', '24', '27', '30' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Dist between Planes', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('STD Intersections (%)', 'fontsize', 15);
colorbar


%%
figure;
imagesc(abs(Matriz_Error_Length));
xticklabels({ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10' });
yticklabels({ '70', '80', '90', '100', '110', '120', '130', '140', '150' });
xlabel('Dist between Planes', 'fontsize', 18);
ylabel('Step Length', 'fontsize', 18);
title('Virtual Plane Mean Error Length (%)', 'fontsize', 15);
% title('Length Estim Error (Virtual Planes)  (%)', 'fontsize', 15);
colorbar


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

ind_e = 9;
ind_d = 10;

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)), '.   Dist between Planes = ', num2str(Distancias_Entre_Planos(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)), '.   Dist between Planes = ', num2str(Distancias_Entre_Planos(ind_d))]);


%% 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
Sampling_Box_X = Grosor_Corte;   % en micras
Sampling_Box_Y = Grosor_Corte;   % en micras

Distancias_Entre_Planos = 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, Distancias_Entre_Planos, Step_Lengths, Grosor_Corte, h_Cortes, Sampling_Box_X, Sampling_Box_Y, 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 = Distancias_Entre_Planos/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













%%