clear all;

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

Num_File = 6;
load(['Neurona_Matlab_', num2str(Num_File), '.mat']);

%% Ploteamos

clc

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

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), '-', 'color', colores{mod(i, 5)+1}, 'markersize', 10, 'linewidth', 1.5 );
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
    % zlim([150 200]);
    % view([-96, -90]);
    %     axis equal
end

%% Generamos cortes de en el plano xy

clc

% Preparamos la muestra en Z para que en ella quepa exactamente la rejilla de ancho "Grosor_C"
Grosor_Corte = 50;   % Espesor del corte = 50 micras
sobra_Z = Grosor_Corte - mod( max( AXON(:, 4) ) - min( AXON(:, 4) ), Grosor_Corte );
Max_Z = max( AXON(:, 4) ) + sobra_Z/2;    % Coordenada máxima de z
min_Z = min( AXON(:, 4) ) - sobra_Z/2;    % Coordenada mínima de z

Cortes = min_Z:Grosor_Corte:Max_Z - Grosor_Corte;

% Preparamos la muestra en Y para que en ella quepa exactamente la rejilla de ancho "Rejilla_Y"
Rejilla_Y = 100;
sobra_Y = Rejilla_Y - mod( max( AXON(:, 3) ) - min( AXON(:, 3) ), Rejilla_Y );
Max_Y = max( AXON(:, 3) ) + sobra_Y/2;    % Coordenada máxima de z
min_Y = min( AXON(:, 3) ) - sobra_Y/2;    % Coordenada mínima de z

Filas_Rejilla = min_Y:Rejilla_Y:Max_Y - Rejilla_Y;

% Preparamos la muestra en X para que en ella quepa exactamente la rejilla de ancho "Rejilla_X":
Rejilla_X = 100;
sobra_X = Rejilla_X - mod( max( AXON(:, 2) ) - min( AXON(:, 2) ), Rejilla_X );
Max_X = max( AXON(:, 2) ) + sobra_X/2;    % Coordenada máxima de z
min_X = min( AXON(:, 2) ) - sobra_X/2;    % Coordenada mínima de z

Columnas_Rejilla = min_X:Rejilla_X:Max_X - Rejilla_X;

tic
% Separamos la muestra en cortes en XY y en la rejilla:
Cortes_XY_Axon = cell( length( Cortes ), 1 );
for i = 1:length( Cortes )  % Recorremos los cortes
    
    dummy_CORTE = {};
    dummy_RAMA_en_REJILLA = cell( 1, length(AXON_Cell) );
    for j = 1:length(AXON_Cell)
        dummy_2 = AXON_Cell{1, j};   % Puntos en esa rama del axón y en ese corte
        ind_Corte = dummy_2(:, 3) > Cortes(i) & dummy_2(:, 3) < Cortes(i) + Grosor_Corte;
        
        if sum(ind_Corte) > 0    % En este corte hay puntos de esa rama
            
            dummy_REJILLA = cell( length(Filas_Rejilla), length(Columnas_Rejilla) );  % Esta celda tendrá NxM, con N y M el númeor de filas y columnas respect. de la rejilla en que se divide cada corte
            % Vemos en qué celda de la rejilla están:
            for filas = 1:length(Filas_Rejilla)
                ind_Rejilla_Y = dummy_2(:, 2) > Filas_Rejilla( filas ) & dummy_2(:, 2) < Filas_Rejilla( filas ) + Rejilla_Y;
                
                for columnas = 1:length( Columnas_Rejilla )
                    ind_Rejilla_X = dummy_2(:, 1) > Columnas_Rejilla( columnas ) & dummy_2(:, 1) < Columnas_Rejilla( columnas ) + Rejilla_X;
                    
                    Rama_en_Rejilla = dummy_2( ind_Corte & ind_Rejilla_X & ind_Rejilla_Y, : );    % Puntos de esa rama que están en esa celda de la rejilla
                    dummy_REJILLA{filas, columnas} = Rama_en_Rejilla;
                end
                
            end
            dummy_RAMA_en_REJILLA{1, j} = dummy_REJILLA;   % Cada celda contiene NxM celdas, codificando la presencia de cada rama correspondiente en la rejilla. Las que no tienen NxM es que están vacías porque sus ramas no están en ese corte
            
        end
        
    end
    Cortes_XY_Axon{i, 1} = dummy_RAMA_en_REJILLA;
    disp(i)
end
toc

%% Salvamos los datos:

save(['Cortes_Muestra_Axon_', num2str(Num_File), '.mat']);

%% Comprobamos que está bien ploteando:

clc
i_x = 16;   % Columna de la rejilla
i_y = 14;   % Fila de la rejilla
i_z = 14;   % Corte

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

figure('color', 'w', 'position', [50, 200, 1300, 500]);

subplot(1, 2, 1);
hold on;
for i = 1:length(AXON_Cell)
    
    Puntos_Axon = AXON_Cell{1, i};
    plot3( Puntos_Axon(:, 1), Puntos_Axon(:, 2), Puntos_Axon(:, 3), '-o', 'color', colores{mod(i, 5)+1}, 'markersize', 7, 'linewidth', 1.5 );
    
    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]);
    %     axis equal
end

subplot(1, 2, 2);
hold on
for i = 1:length(Cortes_XY_Axon{i_z, 1})
    
    if ~isempty( Cortes_XY_Axon{i_z, 1}{1, i} )
        
        dummy = Cortes_XY_Axon{i_z, 1}{1, i};
        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
    
    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]);
    %     axis equal
end