Rubio_etal_2021/Code/Scripts/PRELIMINAR_Analisis_Esferas_Neuronas_1.m

clear all;

% %%
clc
% Leemos el archivo de Matlab conteniendo las neuronas

% Num_File = 3;
% load(['Cortes_Muestra_Axon_', num2str(Num_File), '.mat']);

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
rr = 46;

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

Diams_Sonda = 45;    % Vector con los diámetros de la sonda a considerar. En micras
Step_Lengths = 140;   % Vector con los step lengths a considerar. En micras
lado_x_Frame = 50;   % en micras
lado_y_Frame = 50;   % en micras
sampling_box_height = 50;   % en micras

% Leemos los archivos
fname=fullfile(Path_Data_1, [names_1{rr}]);
fname_new = [fname(1:end-4), '.mat'];
load( fname_new );

[Q, Error_Length, Axon_Real_Length, Estimated_Axon_Length, Filas_Centros, Columnas_Centros, Z_Centros] = 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 );




% %% Longitud total del axón
% 
% branch_length = zeros( 1, length(AXON_Cell) );
% aa= 0;
% for i = 1:length(AXON_Cell)
%     
%     dummy = AXON_Cell{1, i};
%     diff_dummy = diff( dummy, 1, 1 );
%     branch_length( i ) = sum( sqrt( sum( diff_dummy(:, 1:3).^2, 2 ) ) );
%     aa = aa + length(dummy);
% end
% Axon_Length = sum( branch_length );
% 
% %% Obtenemos las intersecciones de la esfera para estimar la longitud total:
% 
% clc
% 
% h_Cortes = 1;   % Analizamos un corte cada 'h_Cortes' cortes de la muestra
% 
% % 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.
% lado_x_Frame = 50;   % en micras
% lado_y_Frame = 50;   % en micras
% sampling_box_height = 50;   % en micras
% 
% % Radio sonda esférica:
% rad_sonda = sampling_box_height/2;
% 
% [Q, Intersecciones_Detalladas, Error_Length, Filas_Centros, Columnas_Centros, Z_Centros] = Script_Analisis_Esferas_Neuronas_1 ( Cortes_XY_Axon, h_Cortes, rad_sonda, Columnas_Rejilla, Rejilla_X, Filas_Rejilla, Rejilla_Y, Cortes, Grosor_Corte, lado_x_Frame, lado_y_Frame, sampling_box_height, Axon_Length );
% 


%% Ploteamos el axón con las esferas

clc

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

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













%%