No BSD License  

Highlights from
Simple Particle Tracking

from Simple Particle Tracking by Martin Pastor
Detection of bright spots into a sequence of images and track its trajectory

enlazar_original(s_raiz,L,epsilon,salto)
function enlazar_original(s_raiz,L,epsilon,salto)

% Este programa sigue la posicin de una bola en un archivo de datos que
% est ordenado tal como _fxya (mirar el programa track_ca_01): elije un
% punto y lo va encontrando el punto ms cercano del siguiente frame
% construyendo as la trayectoria. Se permite que haya frames donde se
% pierda la partcula

% L: longitud mnima de la trayectoria
% epsilon: separacin en pixel de la posicion de un spot de un frame a otro
% salto: nmero de perdidas permitido (recomiendo cero)

% Atencin: el formato con el que se guardan los datos es crudo, es decir,
% posicin "x" e "y" de los centroides en pixeles, es decir, el spot de
% coordenadas (1,1) es el vertice superior izquierdo y el (n_rows,n_columns)
% el vertice inferior derecho.

dir_0='C:\programas';           % origen del programa
dir_1='C:\origen_datos';        %origen de la datos
dir_2='C:\guardar_datos';       % carpeta donde se guardan los datos
dir_3='C:\origen_imagenes';     % origen de la carpeta de imagenes

close all;

n_number=6;
formato='.tif';
coletilla='_fxyA.dat';
pausar=1;

% carga algunos parmetros

param=load(strcat(dir_1,s_raiz,'_param.dat'));
n_rows=param(8);
n_columns=param(9);

% borra los archivos anteriores

delete(strcat(dir_2,s_raiz,'_t*.dat'));
delete(strcat(dir_2,s_raiz,'_cut.dat'));
delete(strcat(dir_2,s_raiz,'_Les.dat'));
delete(strcat(dir_2,s_raiz,'_ind.dat'));

set(0,'defaulttextinterpreter','none');

% carga de los datos de el reconocimiento de espots

dat=load(strcat(dir_1,s_raiz,coletilla));       

tam=length(dat);
frame=dat(:,1);
x=dat(:,2);
y=dat(:,3);
area=dat(:,4);
hecho=zeros(tam+1,1);     % columna que dice que punto ha sido seguido ya


% busca que indices corresponde a cada brame en el archivo de espots, el
% tamao de cut debe ser igual al nmero de frames.

n=1;    % contador absoluto
m=1;    % contador de frames
     
while (n<tam),    
    cut(m,:)=[frame(n) n 0];                % variable que almacena donde empieza y donde termina los datos de cada frame:
                                            % frame     indice_empieza      indice termina
    while((n<tam)&(frame(n)==frame(n+1))),
        n=n+1;
    end
    cut(m,3)=n;
    m=m+1;
    n=n+1;
end

save(strcat(dir_2,s_raiz,'_cut.dat'),'cut','-ASCII');        % guarda el archivo en el que empizan y terminan los frames

% construye las trayectorias: elige un punto y busca dentro del archivo
% cual es el ms prximo en el siguiente frame guardando el indice absoluto
% en la variable indices "ind".

n=1;        % contador absoluto
t=1;        % contador de trayectorias
save_ind=[]; % guardar el vector indices

while (n<tam)
    
        if(hecho(n)==0)

        i=1;        % contador de la variable de indices de trayectoria (ind)
        ind(i)=n;   % variable que contiene los indices
        hecho(n)=1; % cuando un punto se ha hecho se marca con uno
        dat_save=[];    % variable que guardar los datos encontrados
        dat_falso=[];   % variable que guardar los datos que estn interpolados
        p=[x(ind(i)) y(ind(i))];    %punto a seguir
        
        ind_f_ini=find(cut(:,3)>n); % busca el primer frame que todava no ha sido analizado
        f_ini=cut(ind_f_ini(1),1)+1;

        terminar=0;             % cuando se pierde un punto durante mas de "saltos" teminar es 1 y cierra la trayectoria 

        for f=f_ini:length(cut),    % busca en los datos correspondientes a cada frame

            if(terminar<=salto)

                i=i+1;
                exito=0;                % exito es 1 si se encuentra en este frame un punto a una distacia menor que epsilon del punto que estabamos siguiendo
                d_min=2*epsilon;        % inicializa la menor distancia
                ind_ok_prev=tam+1;      % inicializa el incide bueno anterior

                n_follow_ini=cut(f,2);  % se busca desde donde empieza hasta que termina cada frame
                n_follow_fin=cut(f,3);

                for m=n_follow_ini:n_follow_fin,    % busca en los datos en cada frame

                    if(hecho(m)==0)

                        pbis=[x(m) y(m)];
                        d=sqrt((p(1)-pbis(1)).^2+(p(2)-pbis(2))^2);

                        if((d_min>d)&(d<=(epsilon*(terminar+1))))      % si la distancia es menor que la minima y menor que epsilon

                            hecho(ind_ok_prev)=0;       % vuelve a hbilitar el punto anterior
                            ind(i)=m;                   % guarda este punto como el siguiente en la trayectoria
                            ind_ok_prev=m;              % este ser el indice bueno anterior para el siguente frame
                            hecho(m)=1;
                            d_min=d;                    % nueva menor distancia
                            exito=1;

                        end

                    end

                end

                if(exito==1)

                    p=[x(ind(i)) y(ind(i))];        % si no se ha perdido el punto el nuevo punto a seguir es el de menor distancia
                    
                else        % si se pierde el punto se pone 0 en la variable de indices de la trayectoria

                    hecho(ind_ok_prev)=0;
                    terminar=terminar+1;
                    ind(i)=0;               % los puntos que se ha saltado se marcan con 0 en la variable ind
                end

                hasta=i; % numero de frames de la trayectoria

                while (ind(end)==0)    % se eliminan los puntos donde se ha perdido el espot
                    ind(end)=[];
                    hasta=hasta-1;
                    i=i-1;
                end


                % si salto es mayor que cero se hace una interpolacin lineal
                % con los puntos que no se han encontrado

                if((i>salto)&(salto~=0)&(ind([i-salto:i-1])==zeros(1,salto))&(ind(i)~=0))

                    ind_fit=[(i-salto-1) (i)];
                    x_fit=x(ind_fit);
                    y_fit=y(ind_fit);

                    param_x=polyfit(ind_fit,x_fit',1);
                    param_y=polyfit(ind_fit,y_fit',1);

                    x_falso=polyval(param_x,[i-salto:i-1]);
                    y_falso=polyval(param_y,[i-salto:i-1]);
                    matriz_falso=[[(frame(ind(i-salto-1))+1):(frame(ind(i))-1)]' x_falso' y_falso' zeros(salto,1)];
                    dat_falso=[dat_falso;matriz_falso];
                    % los puntos interpolados tienen rea cero
                end


            end

        end
        
        % gurdar las trayectorias si el nmero de puntos es mayor que "L"

        if(length(ind)<L)

            ind=[];

        else

            for o=1:hasta,

                % si el valor del indice encontrado es cero se sustituye por el
                % dato interpolado
                if(ind(o)==0),
                    dat_save(o,:)=dat_falso(1,:);
                    dat_falso(1,:)=[];
                else
                    dat_save(o,:)=dat(ind(o),:);
                end

            end

            
            save(strcat(dir_2,s_raiz,'_t',num2str(t),'.dat'),'dat_save','-ASCII');
            disp(strcat(s_raiz,'_t',num2str(t)));
            n_image=ind(1);
            save_ind=[save_ind ind];
            ind=[];
            t=t+1;

            % representa trayectorias
            t_p=dat_save(:,1);
            x_p=dat_save(:,2);
            y_p=dat_save(:,3);
            A_p=dat_save(:,4);

            rojos=[];
            rojos=find(A_p==0);

            subplot(2,2,1)
            plot(t_p,n_columns-y_p,'.');
            if ~isempty(rojos)
                hold on
                plot(t_p(rojos),n_columns-y_p(rojos),'o');
                hold off
            end
            xlabel('frame');
            ylabel('pixel');
            title(strcat('trajectory ',num2str(t-1),'. y coordinate'));
            subplot(2,2,3)
            plot(t_p,n_rows-x_p,'.');
            if ~isempty(rojos)
                hold on
                plot(t_p(rojos),n_rows-x_p(rojos),'o');
                hold off
            end
            xlabel('frame');
            ylabel('pixel');title(strcat('trajectory ',num2str(t-1),'. x coordinate'));
            
                     % represento el punto sobre la imagen

            n_zeros=n_number-length(num2str(frame(n_image)));    
            cero='0';
            while (length(cero)<n_zeros),
                cero=strcat(cero,'0');
            end
            s=strcat(s_raiz,cero,num2str(frame(n_image)),formato);

            A=imread(strcat(dir_3,s));
            
            subplot(1,2,2)
            imshow(A);
            hold on;
            plot(x(n_image),y(n_image),'or','MarkerSize',5)
            title(s);
            hold off;
            drawnow
            
            pause(pausar);

        end

    end

    n=n+1;
    
end


% guarda el archivo con los parmetros utilizados.

param=[L,epsilon,salto];
save_ind=sort(save_ind');
save(strcat(dir_2,s_raiz,'_Les.dat'),'param','-ASCII');
save(strcat(dir_2,s_raiz,'_ind.dat'),'save_ind','-ASCII');


% representar trayectorias todas juntas

files=dir(strcat(dir_2,s_raiz,'_t*.dat'));
nfiles=length(files);
colores=hsv(nfiles);
clear files;

close
figure
hold on

for n=1:nfiles,
    dat=[];
    name=strcat(dir_2,s_raiz,'_t',num2str(n),'.dat');
    dat=load(name);
    frame=dat(:,1);
    x=dat(:,2);
    y=n_rows-dat(:,3);
    plot(x,y,'Color',colores(n,:));
end

xlabel('pixel');
ylabel('pixel');
title(strcat('Trajectories'));

cd(dir_0)

Contact us at files@mathworks.com