Rubio_etal_2021/Code/Scripts/PRUEBA_I_Script_Analisis_PLANOS_Neuronas_1.m

function [Q, Intersecciones_Detalladas, Error_Length, Estimated_Axon_Length, Filas_Centros, Columnas_Centros, Z_Centros] = PRUEBA_I_Script_Analisis_PLANOS_Neuronas_1 ( Cortes_XY_Axon, h_Cortes, Columnas_Rejilla, Rejilla_X, Filas_Rejilla, Rejilla_Y, Cortes, Grosor_Corte,       Sampling_Box_X, Sampling_Box_Y, Sampling_box_height,       Axon_Real_Length,  dist_Planos,     AZIM_0, ELEV_0, ONSET_0 )

% Definimos los parámetros estereológicos:

% h_Cortes = 2;   % Analizamos un corte cada 'h_Cortes' cortes de la muestra

% Creamos unos vectores con los centros de las celdas de la rejilla en que hemos dividido el corte; en esos centros irá la sonda esférica:
Filas_Centros = Filas_Rejilla + Rejilla_Y/2;
Columnas_Centros = Columnas_Rejilla + Rejilla_X/2;
Z_Centros = Cortes + Grosor_Corte/2;

Intersecciones_Detalladas = cell(size(Cortes_XY_Axon));
Q = 0;   % Intersecciones totales
% Recorremos los cortes para analizar:
for i = 1:h_Cortes:length(Cortes_XY_Axon)
    
    % Para cada corte generamos tres matrices, las dos primeras con los ángulos de azimutal y elevación de los planos en cada sampling box (o en cada celda de la rejilla):
    AZIM = AZIM_0(:, :, i);
    ELEV = ELEV_0(:, :, i);
%     AZIM = pi - 2*pi*rand( length(Filas_Centros), length(Columnas_Centros) );
%     ELEV = pi/2 - pi*rand( length(Filas_Centros), length(Columnas_Centros) );
    % La tercera con el onset, i.e. 
%     ONSET = dist_Planos*rand( length(Filas_Centros), length(Columnas_Centros) );
    ONSET = dist_Planos*ONSET_0(:, :, i);
    
    % 'dummy_CORTE' contiene la distribución de ramas axónicas en cada celda de la rejilla de ese corte:
    dummy_CORTE = Cortes_XY_Axon{i, 1};
    % 'dummy_INTERSECCIONES' contiene las intersecciones en cada celda de la rejilla en que se dividió el corte. En cada celda, cada 
    % intersección está codificada como un número que es la rama que en esa celda intersecta con la esfera:
    dummy_INTERSECCIONES = cell( length( Filas_Rejilla ), length( Columnas_Rejilla ) );
    ind_CORTE = find(~cellfun(@isempty, dummy_CORTE));     % índices de las ramas que aparecen en este corte
    for j = ind_CORTE   % Recorremos las ramas de ese corte
        
        [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
%         row_a_Y = length( Filas_Rejilla ) - row_a + 1;   % Le damos la vuela porque los índices empiezan por arriba (índices matriciales), pero hemos generado la rejilla empezando por abajo en la Y, i.e. por los negativos.
        for k = 1:length(row_a)   % Recorremos los índices de las celdas de ese corte en el que aparece la presente rama
            
            dummy_a = dummy_CORTE{1, j}{row_a(k), col_a(k)};    % Coordenadas de los puntos de la rama cuya intersección estamos estudiando
            % Nos quedamos sólo con los puntos de la rama que están dentro de la sampling box:
            ind_sampling_box = (abs( dummy_a(:, 1) - Columnas_Centros(col_a(k)) ) <= Sampling_Box_X/2) &...
                                              (abs( dummy_a(:, 2) - Filas_Centros(row_a(k)) ) <= Sampling_Box_Y/2) &...
                                              (abs( dummy_a(:, 3) - Z_Centros(i) )) <= Sampling_box_height/2;
            dummy_b = dummy_a( ind_sampling_box, : );
            % Obtenemos los cortes entre ese trozo de rama y la sampling box, incluidos los puntos de dentro, para coger toda la longitud de la rama:
%             dummy_b = Corte_Rama_Samplng_Box( dummy_a, Columnas_Centros(col_a(k)), Filas_Centros(row_a(k)), Z_Centros(i), Sampling_Box_X, Sampling_Box_Y, Sampling_box_height );
            
            if ~isempty( dummy_b )    % Si hay algún punto de la rama en la sampling box pasamos a los cortes con los planos
                % Num intersecciones con los planos
                num_inters = Cortes_con_Planos_Virtuales( AZIM( row_a(k), col_a(k) ), ELEV( row_a(k), col_a(k) ), ONSET( row_a(k), col_a(k) ), ( dummy_b(:, 1:3) )', dist_Planos );
                
                % Metemos los índices de las ramas que intersectan los planos, en la celda adecuada de 'dummy_INTERSECCIONES':
                if num_inters > 0
                    dummy_INTERSECCIONES{row_a(k), col_a(k)} = [dummy_INTERSECCIONES{row_a(k), col_a(k)}; repmat(j, num_inters, 1) ];
                    % Añadimos estas intersecciones al número total:
                    Q = Q + num_inters;
                end
            end
            
        end
        
    end
    Intersecciones_Detalladas{i, 1} = dummy_INTERSECCIONES;
    
end

% %% Estimación de la longuitud del axón

% clc

% Section Sampling Fraction ( ssf ) = the fraction of the total number of sections containing the reference space 
% that were actually used in the analysis, es decir, cada cuantos cortes tomamos uno para analizar con las esferas.
% Por ejemplo, si tomamos 1 corte de cada 8 --> ssf = 1/8;
ssf = 1/h_Cortes;

% Area Sampling Fraction ( asf ) = the fraction of the area of the profiles of the reference volume seen in the sections 
% that were sampled by the sampling boxes that contained the spherical probes (profiles son los que aparece del axón 
% cuando éste se corta con un plano de enfoque paralelo a XY); ver Figs. 6, 7 y 8 de "Stereological length estimation using spherical probes";
% El 'frame' es la caja 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.
% Es el area_frame / area x-y step
                            % lado_x_Frame = 50;   % en micras
                            % lado_y_Frame = 50;   % en micras
asf = Sampling_Box_X*Sampling_Box_Y / ( Rejilla_X*Rejilla_Y );

% Thickness Sampling Fraction ( tsf ) = the fraction of the final, i.e. after-processing, section thickness sampled by 
% the sampling boxes = 1//(sampling box height/section thickness).
                            % sampling_box_height = 50;   % en micras
hsf = Sampling_box_height / Grosor_Corte;

% Longitud axón estimada:
Estimated_Axon_Length = 2*Q*dist_Planos*( 1/ssf )*( 1/asf )*( 1/hsf );

% %% Error de la esitimación

Error_Length = ( Axon_Real_Length - Estimated_Axon_Length )*100/Axon_Real_Length;