Rubio_etal_2021/Code/Scripts/Export_Axons_Neuromorph_3.m

%

function Export_Axons_Neuromorph_3( filename )
% Leemos un archivo .swc. Obtenemos una tabla

LL = importdata(filename, ' ');
if ~isfield(LL, 'data')
    
    warning('¡Hay que importar a mano!');
    disp(filename(end-13:end));
    
else
    
    A = LL.data;
    if size(A, 2) ~= 7 && size(A, 2) ~= 3
        warning('¡Hay que importar a mano!');
        disp(filename(end-13:end));
    else
        if size(A, 2) == 3
            A = Preparar_SWC_Raro( LL );
        end
        A = A(:, 2:end);
        
        AXON = A(A(:, 1) == 2, :);
        
        % Organizamos el axon por ramas como una cell.
        
        % Buscamos los ínidices de la columna 6 que estén repetidos; esos indicarán branching y por lo tanto señalarán ramas distintas:
        
        [~, idx ] = unique(AXON(:, 6), 'stable');%don't sort
        fila = setxor(1:numel(AXON(:, 6)), idx);
        % ind2 = AXON(idx, 6);
        
        AXON_Cell = {};
        % Recorremos de mos a más los índices de los puntos en el axón:
        rama = 1;
        dummy = [];
        Ind_Repet = [];
        for i = 1:size(AXON, 1)
            
            if ~ismember( AXON(i, 6), AXON(fila, 6) )   % El punto del AXON en el que estamos NO está repetido; seguimos almacenando los puntos de la rama
                dummy = [dummy; [AXON(i, 2), AXON(i, 3), AXON(i, 4), AXON(i, 6)]];
                if i == size(AXON, 1)    % Se acaban los índices y guardamos la rama
                    % Pegamos la rama al árbol:
                    ind0 = find( AXON(:, 2) == dummy(1, 1) & AXON(:, 3) == dummy(1, 2) & AXON(:, 4) == dummy(1, 3) );   % fila del primer elemento de la rama
                    if ind0 > 1
                        if AXON(ind0-1, 6) + 1 == dummy(1, 4)
                            ind1 = ind0 - 1;    % Tomamos el índice del punto anterior al primero de la rama
                        else  % se ramifica.
                            ind1 = find( AXON(:, 6) == dummy(1, 4) );   % debe tener dos elementos;
                            ind1 = ind1(1);  % tomamos el primero
                        end
                        dummy = [[AXON(ind1, 2), AXON(ind1, 3), AXON(ind1, 4), AXON(ind1, 6)]; dummy];
                    end
                    
                    AXON_Cell{1, rama} = dummy;  % Guardamos la rama que acaba
                end
            else  % El punto de AXON en el que estamos SÍ está repetido; creamos una rama nueva
                %                 AXON_Cell{1, rama} = [dummy; [AXON(i, 2), AXON(i, 3), AXON(i, 4), AXON(i, 6)]];  % Guardamos la rama que acaba
                % Empalmamos esta rama al árbol. Tomamos el índice del primer punto de la rama. Si ese índice es mayor que 1 y el índice del anterior
                % punto + 1 es igual a ese, entonces ponemos el anterior punto primero en la rama. Si no lo es, significa que el anterior punto de la rama
                % es una bifurcación; la buscamos y la ponemos en primer lugar en la rama:
%                 if AXON(i, 6) == 1
                if i == 1
                    dummy = [AXON(i, 2), AXON(i, 3), AXON(i, 4), AXON(i, 6)];
                end
                ind0 = find( AXON(:, 2) == dummy(1, 1) & AXON(:, 3) == dummy(1, 2) & AXON(:, 4) == dummy(1, 3) );   % fila del primer elemento de la rama
                if ind0 > 1
                    if AXON(ind0-1, 6) + 1 == dummy(1, 4)
                        ind1 = ind0 - 1;    % Tomamos el índice del punto anterior al primero de la rama
                    else  % se ramifica.
                        ind1 = find( AXON(:, 6) == dummy(1, 4) );   % debe tener dos elementos;
                        ind1 = ind1(1);  % tomamos el primero
                    end
                    dummy = [[AXON(ind1, 2), AXON(ind1, 3), AXON(ind1, 4), AXON(ind1, 6)]; dummy];
                end
                
                AXON_Cell{1, rama} = dummy;  % Guardamos la rama que acaba
                
                % Comenzamos una nueva rama
                rama = rama + 1;   % índice de la rama en la que nos encontramos
                dummy = [AXON(i, 2), AXON(i, 3), AXON(i, 4), AXON(i, 6)];   % Guardamos en la rama como primer punto, el punto repetido
                 if i == size(AXON, 1)    % Se acaban los índices y guardamos la rama
                    % Pegamos la rama al árbol:
                    ind0 = find( AXON(:, 2) == dummy(1, 1) & AXON(:, 3) == dummy(1, 2) & AXON(:, 4) == dummy(1, 3) );   % fila del primer elemento de la rama
                    if ind0 > 1
                        if AXON(ind0-1, 6) + 1 == dummy(1, 4)    % En este caso el punto anterior continúa a este punto (no hay ramificación)
                            ind1 = ind0 - 1;    % Tomamos el índice del punto anterior al primero de la rama
                        else                            
                            if AXON(ind0-1, 6) == dummy(1, 4)  % Este punto y el anterior tienen el mísmo parent index
                                ind1 = ind0-2;
                            else% se ramifica.
                                ind1 = find( AXON(:, 6) == dummy(1, 4) );   % debe tener dos elementos;
                                ind1 = ind1(1);  % tomamos el primero
                            end
                        end
                        dummy = [[AXON(ind1, 2), AXON(ind1, 3), AXON(ind1, 4), AXON(ind1, 6)]; dummy];
                    end
                    
                    AXON_Cell{1, rama} = dummy;  % Guardamos la rama que acaba                     
                end
                Ind_Repet = [Ind_Repet; i];
                
            end
            
        end
    
        % Salvamos los datos
        filename_new = filename(1:end-4); % Quitamos la extensión .swc
        save([filename_new, '.mat']);
        
    end
end