track.m

%;
% ; see http://glinda.lrsm.upenn.edu/~weeks/idl
% ;   for more information
% ;
% ;+
% ; NAME:
% ;	track
% ; PURPOSE:
% ;	Constructs n-dimensional trajectories from a scrambled list of
% ;	particle coordinates determined at discrete times (e.g. in
% ;	consecutive video frames).
% ; CATEGORY:
% ;	Image Processing
% ; CALLING SEQUENCE:
% ;	result = track( positionlist, maxdisp [, inipos=inipos, 
% ;	   	memory=memory, dim=dim, verbose=verbose, 
% ;		goodenough=goodenough ] )
% ; INPUTS:
% ;	positionlist: an array listing the scrambled coordinates and data 
% ;		of the different particles at different times, such that:
% ;	 positionlist(0:d-1,*):	contains the d coordinates and
% ;		data for all the particles, at the different times.
% ;	 positionlist(d,*):	contains the time t that the position 
% ;		was determined (e.g. frame number.  These values must 
% ;               be monotonically increasing and uniformly gridded in time.
% ;	maxdisp: an estimate of the maximum distance that a particle 
% ;		would move in a single time interval.(see Restrictions)
% ;	inipos: if the user wants to track only a subset of the 
% ;		particles, this argument is set to an array (d,n)
% ;		which contains the d dimensional initial positions of
% ;		the n particles to be tracked. Other 'new' particles
% ;		will then NOT be added.
% ;	memory: this is the number of time steps that a particle can be
% ;		'lost' and then recovered again.  If the particle reappears
% ;		after this number of frames has elapsed, it will be
% ;		tracked as a new particle. The default setting is zero.
% ;		this is useful if particles occasionally 'drop out' of
% ;		the data.
% ;	dim: if the user would like to unscramble non-coordinate data
% ;		for the particles (e.g. apparent radius of gyration for
% ;		the particle images), then positionlist should
% ;		contain the position data in positionlist(0:dim-1,*)
% ;		and the extra data in positionlist(dim:d-1,*). It is then
% ;		necessary to set dim equal to the dimensionality of the
% ;		coordinate data to so that the track knows to ignore the
% ;		non-coordinate data in the construction of the 
% ;		trajectories. The default value is two.
% ;	verbose: set this keyword for more informational messages.
% ;	goodenough: set this keyword to eliminate all trajectories with
% ;		fewer than goodenough valid positions.  This is useful
% ;		for eliminating very short, mostly 'lost' trajectories
% ;		due to blinking 'noise' particles in the data stream.
% ; OUTPUTS:
% ;	result:  a list containing the original data rows sorted 
% ;		into a series of trajectories.  To the original input 
% ;		data structure there is appended an additional column 
% ;		containing a unique 'id number' for each identified 
% ;		particle trajectory.  The result array is sorted so 
% ;		rows with corresponding id numbers are in contiguous 
% ;		blocks, with the time variable a monotonically
% ;		increasing function inside each block.  For example:
% ;		
% ;		For the input data structure (positionlist):
% ;			(x)	     (y)	  (t)
% ;		pos = 3.60000      5.00000      0.00000
% ;		      15.1000      22.6000      0.00000
% ;		      4.10000      5.50000      1.00000	
% ;		      15.9000      20.7000      2.00000
% ;		      6.20000      4.30000      2.00000
% ;
% ;		IDL> res = track(pos,5,mem=2)
% ;
% ;		track will return the result 'res'
% ;			(x)	     (y)	  (t) 	       (id)
% ;		res = 3.60000      5.00000      0.00000      0.00000
% ;		      4.10000      5.50000      1.00000      0.00000
% ;		      6.20000      4.30000      2.00000      0.00000
% ;		      15.1000      22.6000      0.00000      1.00000
% ;		      15.9000      20.7000      2.00000      1.00000
% ;
% ;		NB: for t=1 in the example above, one particle temporarily
% ;		vanished.  As a result, the trajectory id=1 has one time
% ;		missing, i.e. particle loss can cause time gaps to occur 
% ;		in the corresponding trajectory list. In contrast:
% ;
% ;		IDL> res = track(pos,5)
% ;
% ;		track will return the result 'res'
% ;			(x)	     (y)	  (t) 	       (id)
% ;		res = 15.1000      22.6000      0.00000      0.00000
% ;      		      3.60000      5.00000      0.00000      1.00000
% ;  		      4.10000      5.50000      1.00000      1.00000
% ; 		      6.20000      4.30000      2.00000      1.00000
% ; 		      15.9000      20.7000      2.00000      2.00000
% ;	
% ;		where the reappeared 'particle' will be labelled as new
% ;		rather than as a continuation of an old particle since
% ;		mem=0.  It is up to the user to decide what setting of 
% ;		'mem' will yeild the highest fidelity tracking.
% ;	
% ; SIDE EFFECTS:
% ;	Produces informational messages.  Can be memory intensive for
% ;	extremely large data sets.
% ; RESTRICTIONS:
% ;	maxdisp should be set to a value somewhat less than the mean 
% ;	spacing between the particles. As maxdisp approaches the mean
% ;	spacing the runtime will increase significantly. The function 
% ;	will produce an error message: "Excessive Combinatorics!" if
% ;	the run time would be too long, and the user should respond 
% ;	by re-executing the function with a smaller value of maxdisp.
% ;	Obviously, if the particles being tracked are frequently moving
% ;	as much as their mean separation in a single time step, this
% ;	function will not return acceptable trajectories.
% ; PROCEDURE:
% ;	Given the positions for n particles at time t(i), and m possible
% ;	new positions at time t(i+1), this function considers all possible 
% ;	identifications of the n old positions with the m new positions,
% ;	and chooses that identification which results in the minimal total
% ;	squared displacement. Those identifications which don't associate
% ;	a new position within maxdisp of an old position ( particle loss )
% ;	penalize the total squared displacement by maxdisp^2. For non-
% ;	interacting Brownian particles with the same diffusivity, this
% ;	algorithm will produce the most probable set of identifications 
% ;	( provided maxdisp >> RMS displacement between frames ).
% ;	In practice it works reasonably well for systems with oscillatory,
% ;	ballistic, correlated and random hopping motion, so long as single 
% ;	time step displacements are reasonably small.  NB: multidimensional
% ;	functionality is intended to facilitate tracking when additional
% ;	information regarding target identity is available (e.g. size or 
% ;	color).  At present, this information should be rescaled by the
% ;	user to have a comparable or smaller (measurement) variance than 
% ;	the spatial displacements.
% ;
% ; MODIFICATION HISTORY:
% ;	 2/93 Written by John C. Crocker, University of Chicago (JFI).
% ;	 7/93 JCC fixed bug causing particle loss and improved performance
% ;		for large numbers of (>100) particles.
% ;	11/93 JCC improved speed and memory performance for large
% ;		numbers of (>1000) particles (added subnetwork code).
% ;	 3/94 JCC optimized run time for trivial bonds and d<7. (Added
% ;		d-dimensional raster metric code.)
% ;	 8/94 JCC added functionality to unscramble non-position data
% ;		along with position data.
% ;	 9/94 JCC rewrote subnetwork code and wrote new, more efficient 
% ;		permutation code.
% ;	 5/95 JCC debugged subnetwork and excessive combinatorics code.
% ;	12/95 JCC added memory keyword, and enabled the tracking of
% ;		newly appeared particles.
% ;	 3/96 JCC made inipos a keyword, and disabled the adding of 'new'
% ;		particles when inipos was set.
% ;	 3/97 JCC added 'add' keyword, since Chicago users didn't like 
% ;		having particle addition be the default. 
% ;	 9/97 JCC added 'goodenough' keyword to improve memory efficiency
% ;		when using the 'add' keyword and to filter out bad tracks.
% ;       10/97 JCC streamlined data structure to speed runtime for >200 
% ;               timesteps.  Changed 'quiet' keyword to 'verbose'. Made
% ;               time labelling more flexible (uniform and sorted is ok).
% ;	 9/98 JCC switched trajectory data structure to a 'list' form,
% ;		resolving memory issue for large, noisy datasets.
% ;	 2/99 JCC added Eric Weeks's 'uberize' code to post-facto 
% ;		rationalize the particle id numbers, removed 'add' keyword.
% ;
% ;	This code 'track.pro' is copyright 1999, by John C. Crocker. 
% ;	It should be considered 'freeware'- and may be distributed freely 
% ;	(outside of the military-industrial complex) in its original form 
% ;	when properly attributed.
% ;
% ;-


function info = track(xyzs,maxdisp)
%xyzs = read_gdf('/Users/blair/data/test_data.gdf');
%maxdisp = 4.5
dd = length(xyzs(1,:));
dim = dd - 1;
% checking the input time vector
t = xyzs(:,dd);
st = shiftr(t,-1);
st = t(2:end) - st(2:end);
if  sum(st(find(st < 0))) ~= 0
    disp('The time vectors is not in order')
    return
end
info = 1;
memory_b = 1300;
w = find(st > 0);
z = length(w);
z = z +1;
if isempty(w)
    disp('All positions are at the same time... go back!')
    return
end

% partitioning the data with unique times

res = unq(t);
res = [1,res',length(t)];
ngood = res(2) - res(1) + 1;
eyes = 1:ngood;
pos = xyzs(eyes,1:dim);
istart = 2;
n = ngood;

zspan = 50;
if n > 200 
    zspan = 20;
end
if n > 500 
    zspan = 10;
end
resx = zeros(zspan,n) - 1;

bigresx = zeros(z,n) - 1;
mem = zeros(n,1);
%  whos resx
%  whos bigresx
uniqid = 1:n;
maxid = n;
olist = [0.,0.];
%  whos eyes;
resx(1,:) = eyes;
% setting up constants
maxdisq = maxdisp^2;

% John calls this the setup for "fancy code" ???
notnsqrd = (sqrt(n*ngood) > 200) & (dim < 7);
notnsqrd = notnsqrd(1);

if notnsqrd
    %;   construct the vertices of a 3x3x3... d-dimensional hypercube
    
    cube = zeros(3^dim,dim);
    
    
    for d=0:dim-1,
        numb = 0;
        for j=0:(3^d):(3^dim)-1,
            cube(j+1:j+(3^(d)),d+1) = numb;
            numb = mod(numb+1,3);
        end
    end    
    
    %   calculate a blocksize which may be greater than maxdisp, but which
    %   keeps nblocks reasonably small.  
    
    volume = 1;
    for d = 0:dim-1
        minn = min(xyzs(w,d+1));
        maxx = max(xyzs(w,d+1));
        volume = volume * (maxx-minn);
    end
    volume;
    blocksize = max( [maxdisp,((volume)/(20*ngood))^(1.0/dim)] );
end
%   Start the main loop over the frames.
for i=istart:z
    ispan = mod(i-1,zspan)+1;
    %disp(ispan)
    % get new particle positions
    m = res(i+1) - res(i);
    % res(i);
    eyes = 1:m;
    eyes = eyes + res(i);
    
    if m > 0
        
        xyi = xyzs(eyes,1:dim);
        found = zeros(m,1);
        
        % THE TRIVIAL BOND CODE BEGINS   
        
        if notnsqrd
            %Use the raster metric code to do trivial bonds
            
            % construct "s", a one dimensional parameterization of the space 
            % which consists of the d-dimensional raster scan of the volume.)
            
            abi = fix(xyi./blocksize);
            abpos = fix(pos./blocksize);
            si = zeros(m,1);
            spos = zeros(n,1);
            dimm = zeros(dim,1);
            coff = 1.;
            
            for j=1:dim
                minn = min([abi(:,j);abpos(:,j)]);
                maxx = max([abi(:,j);abpos(:,j)]);
                abi(:,j) = abi(:,j) - minn;
                abpos(:,j) = abpos(:,j) - minn;
                dimm(j) = maxx-minn + 1;
                si = si + abi(:,j).*coff;
                spos = spos + abpos(:,j).*coff;
                coff = dimm(j).*coff;
            end
            nblocks = coff;
            % trim down (intersect) the hypercube if its too big to fit in the
            % particle volume. (i.e. if dimm(j) lt 3)
            
            cub = cube;
            deg = find( dimm < 3);
            if ~isempty(deg)
                for j = 0:length(deg)-1
                    cub = cub(find(cub(:,deg(j+1)) < dimm(deg(j+1))));
                end
            end	
            
            % calculate the "s" coordinates of hypercube (with a corner @ the origin)
            scube = zeros(length(cub(:,1)),1);
            coff = 1;
            for j=1:dim
                scube = scube + cub(:,j).*coff;
                coff = coff*dimm(j);	  
            end
            
            % shift the hypercube "s" coordinates to be centered around the origin
            
            coff = 1;
            for j=1:dim
                if dimm(j) > 3
                    scube = scube - coff;
                end
                coff = dimm(j).* coff;
            end
            scube = mod((scube + nblocks),nblocks);
            % get the sorting for the particles by their "s" positions.
            [ed,isort] = sort(si);
            
            % make a hash table which will allow us to know which new particles
            % are at a given si.
            strt = zeros(nblocks,1) -1;
            fnsh = zeros(nblocks,1);
            h = find(si == 0);
            lh = length(h);
            if lh > 0
                si(h) = 1;  
            end
            
            for j=1:m
                if strt(si(isort(j))) == -1
                    strt(si(isort(j))) = j;
                    fnsh(si(isort(j))) = j;
                else
                    fnsh(si(isort(j))) = j;
                end
            end
            if lh > 0
            si(h) = 0;   
            end
            coltot = zeros(m,1);
            rowtot = zeros(n,1);
            which1 = 1:n;%zeros(n,1);
            for j=1:n
                map = fix(-1);
                
                scub_spos = scube + spos(j);
                s = mod(scub_spos,nblocks);
                whzero = find(s == 0 );
                if ~isempty(whzero)
                    s(whzero) = 1;
                end
               
                w = find(strt(s) ~= -1);
               
                ngood = length(w);
                if ngood ~= 0 
                    s = s(w);
                    for k=1:ngood
                        map = [map;isort( strt(s(k)):fnsh(s(k)))];
                    end
                    map = map(2:end);
                    if length(map) == 2
                        if (map(1) - map(2)) == 0
                            map = unique(map);
                         end
                    end
                    %   map = map(umap);
                    %end
                    % find those trival bonds
                    distq = zeros(length(map),1);
                    for d=1:dim
                        
                        distq = distq + (xyi(map,d) - pos(j,d)).^2;
                        
                    end
                    ltmax = distq < maxdisq;
                  
                    rowtot(j) = sum(ltmax);
                    
                    if rowtot(j) >= 1 
                        w = find(ltmax);
                        coltot( map(w) ) = coltot( map(w)) +1;
                        which1(j) = map( w(1) );
                    end
                end
            end
          
            ntrk = fix(n - sum(rowtot == 0));
        
            w = find( rowtot == 1);
            ngood = length(w);
            %disp('bonds')
           
            if ngood ~= 0 
                ww = find(coltot( which1(w) ) == 1);
                ngood = length(ww);
                if ngood ~= 0 
                     %disp(size(w(ww)))
                    resx(ispan,w(ww)) = eyes( which1(w(ww)));
                    found( which1( w(ww))) = 1;
                    rowtot( w(ww)) = 0;
                    coltot( which1(w(ww))) = 0;
                end
            end
            
            labely = find( rowtot > 0);
            ngood = length(labely);
            if ngood ~= 0 
                labelx = find( coltot > 0);
                
                nontrivial = 1;
            else
                nontrivial = 0;
            end
          %disp(size(resx))
        else 
            % 'this is in the trival code'           
            %   or: Use simple N^2 time routine to calculate trivial bonds      
            
            % let's try a nice, loopless way!
            % don't bother tracking perm. lost guys.
            wh = find( pos(:,1) >= 0);
            ntrack = length(wh);
            if ntrack == 0 
                'There are no valid particles to track idiot!'
                break
            end
            %%%%damon changes here 10/4/04 -- to end of next comment
            xmat=ones(ntrack,1)*[1:m];
            ymat=[1:ntrack]'*ones(1,m);
            %xmat = zeros(ntrack,m);
            %count = 0;
            %for k=1:ntrack
            %    for j=1:m
            %        xmat(k,j) = count;
            %        count = count+1;
            %    end
            %end
            
            %xmat = mod(xmat,m) + 1;
            %ymat = xmat';
            for d=1:dim
                x = xyi(:,d);
                y = pos(wh,d);
                
                if d == 1
                    %keyboard;
                    dq = (x(xmat)-y(ymat)).^2;
                else
                    dq = dq + (x(fix(xmat))-y(fix(ymat))).^2;
                end
            end
            
            ltmax = dq < maxdisq;
            
            % figure out which trivial bonds go with which
            
            rowtot = zeros(n,1);
            %keyboard;
            rowtot(wh) = sum(ltmax,2); % damon changed '1' to '2' in this line, 10/4/04
            
            
            if ntrack > 1 
                coltot = sum(ltmax,1);  % damon changed '2' to '1' in this line, 10/4/04
            else
                coltot = ltmax;
            end
            which1 = zeros(n,1);
            for j=1:ntrack 
                [mx, w] = max(ltmax(j,:));
                which1(wh(j)) = w;
            end
            
            ntrk = fix( n - sum(rowtot(find(rowtot == 0))));
            w= find( rowtot == 1) ;
            ngood = length(w);
            
            if ngood ~= 0
                %keyboard;
                ww = find(coltot(which1(w)) == 1);
                ngood = length(ww);
                if ngood ~= 0 
                    resx( ispan, w(ww) ) = eyes( which1( w(ww)));
                    found(which1( w(ww))) = 1;
                    rowtot(w(ww)) = 0;
                    coltot(which1(w(ww))) = 0;
                end
            end
            
            labely = find( rowtot > 0);
            ngood = length(labely);
            
            if ngood ~= 0
                labelx = find( coltot > 0);
                nontrivial = 1;
            else
                nontrivial = 0;
            end
            
        end
        
        %THE TRIVIAL BOND CODE ENDS
        
        if nontrivial
            
            xdim = length(labelx);
            ydim = length(labely);
            
            %  make a list of the non-trivial bonds            
            
            bonds = zeros(1,2);
            bondlen = 0;
            
            for j=0:ydim-1
                distq = zeros(xdim,1);
                
                for d=1:dim
                    %distq
                    distq = distq + (xyi(labelx,d) - pos(labely(j+1),d)).^2;
                    %distq    
                end
                
                w= find(distq <  maxdisq)' - 1;
                ngood = length(w);
                newb = [w;(zeros(1,ngood)+j)];
                
                
                bonds = [bonds;newb'];
               
                bondlen = [ bondlen;distq( w + 1) ];
                
            end
            bonds = bonds(2:end,:);
            
            bondlen = bondlen(2:end);
            numbonds = length(bonds(:,1));
            mbonds = bonds;
            max([xdim,ydim]);
                
                
            if max([xdim,ydim]) < 4
                nclust = 1;
                maxsz = 0;
                mxsz = xdim;
                mysz = ydim;
                bmap = zeros(length(bonds(:,1)+1),1) - 1;
               
            else
           

                %   THE SUBNETWORK CODE BEGINS            
                lista = zeros(numbonds,1);
                listb = zeros(numbonds,1);
                nclust = 0;
                maxsz = 0;
                thru = xdim;
                
                while thru ~= 0
                    %  the following code extracts connected 
                    %   sub-networks of the non-trivial 
                    %   bonds.  NB: lista/b can have redundant entries due to 
                    %   multiple-connected subnetworks	    
                    
                    bonds;
                    w = find(bonds(:,2) >= 0);
                    size(w);
                    
                    lista(1) = bonds(w(1),2);
                    listb(1) = bonds(w(1),1);
                    bonds(w(1),:) = -(nclust+1);
                    bonds;
                    adda = 1;
                    addb = 1;
                    donea = 0;
                    doneb = 0;
                    if (donea ~= adda) | (doneb ~= addb)
                        true = 0;
                    else
                    true = 1;   
                    end
                                        
                    while ~true
                        
                        if (donea ~= adda)
                            w = find(bonds(:,2) == lista(donea+1));
                            ngood = length(w);
                            if ngood ~= 0 
                                listb(addb+1:addb+ngood,1) = bonds(w,1);
                                bonds(w,:) = -(nclust+1);
                                addb = addb+ngood;
                            end
                            donea = donea+1;
                        end
                        if (doneb ~= addb) 
                            w = find(bonds(:,1) == listb(doneb+1));
                            ngood = length(w);
                            if ngood ~= 0
                                lista(adda+1:adda+ngood,1) = bonds(w,2);
                                bonds(w,:) = -(nclust+1);
                                adda = adda+ngood;
                            end
                            doneb = doneb+1;
                        end
                      if (donea ~= adda) | (doneb ~= addb) 
                          true = 0;
                      else  
                          true = 1;
                      end
                    end
                    
                    [pp,pqx] = sort(listb(1:doneb));
                    unx =  unq(listb(1:doneb),pqx);
                    xsz = length(unx);
                    [pp,pqy] = sort(lista(1:donea));
                    uny =  unq(lista(1:donea),pqy);
                    ysz = length(uny);
                    if xsz*ysz > maxsz
                        maxsz = xsz*ysz;
                        mxsz = xsz;
                        mysz = ysz; 
                    end
                    
                    
                    thru = thru -xsz;
                    nclust = nclust + 1;
                end
                bmap = bonds(:,2);                    
            end
            
            % THE SUBNETWORK CODE ENDS
            % put verbose in for Jaci
            
            %   THE PERMUTATION CODE BEGINS
            
            for nc =1:nclust
                w = find( bmap == -1*(nc));
                
                nbonds = length(w);
                bonds = mbonds(w,:);
                lensq = bondlen(w);
                [pq,st] = sort( bonds(:,1));
                un = unq(bonds(:,1),st);
                
                uold = bonds(un,1);
                
                nold = length(uold);
                
                un = unq(bonds(:,2));
                
                unew = bonds(un,2);
                nnew = length(unew);
                
                if nnew > 5
                    rnsteps = 1;
                    for ii =1:nnew
                        rnsteps = rnsteps * length( find(bonds(:,2) == ...
                            unew(ii)));
                        if rnsteps > 5.e+4
                            disp('Warning: difficult combinatorics encountered.')
                        end
                        if rnsteps > 2.e+5
                            disp(['Excessive Combinitorics you FOOL LOOK WHAT YOU HAVE' ...
                                    ' DONE TO ME!!!'])
                            return
                        end
                    end
                end
                st = zeros(nnew,1);
                fi = zeros(nnew,1);
                h = zeros(nbonds,1);
                ok = ones(nold,1);
                nlost = (nnew - nold) > 0;
                
                
                for ii=1:nold 
                    h(find(bonds(:,1) == uold(ii))) = ii;
                end
                st(1) = 1 ;
                fi(nnew) = nbonds; % check this later
                if nnew > 1 
                    sb = bonds(:,2);
                    sbr = shiftr(sb,-1);
                    sbl = shiftr(sb,1);
                    st(2:end) = find( sb(2:end) ~= sbr(2:end)) + 1;
                    fi(1:nnew-1) = find( sb(1:nbonds-1) ~= sbl(1:nbonds-1));
                end
%                if i-1 == 13
%                    hi
%                end
                checkflag = 0;
                while checkflag ~= 2
                    
                    pt = st -1;
                    lost = zeros(nnew,1);
                    who = 0;
                    losttot = 0;
                    mndisq = nnew*maxdisq;
                    
                    
                    while who ~= -1
                    
                        if pt(who+1) ~= fi(who+1)
                            
                            
                            w = find( ok( h( pt( who+1 )+1:fi( who+1 ) ) ) ); % check this -1
                            ngood = length(w);
                            if ngood > 0
                                if pt(who+1) ~= st(who+1)-1
                                    ok(h(pt(who+1))) = 1;
                                end
                                pt(who+1) = pt(who+1) + w(1);
                                ok(h(pt(who+1))) = 0;
                                if who == nnew -1
                                    ww = find( lost == 0);
                                    dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                    
                                    if dsq < mndisq 
                                        minbonds = pt(ww);
                                        mndisq = dsq;
                                    end
                                else
                                    who = who+1;
                                end
                            else
                                if ~lost(who+1) & (losttot ~= nlost)
                                    lost(who+1) = 1;
                                    losttot = losttot + 1;
                                    if pt(who+1) ~= st(who+1) -1;
                                        ok(h(pt(who+1))) = 1;
                                    end
                                    if who == nnew-1
                                        ww = find( lost == 0);
                                        dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                        if dsq < mndisq
                                            minbonds = pt(ww);
                                            mndisq = dsq;
                                        end
                                    else    
                                       who = who + 1;
                                    end
                                   
                                else
                                    if pt(who+1) ~= (st(who+1) -1) 
                                        ok(h(pt(who+1))) = 1;
                                    end
                                    pt(who+1) = st(who+1) -1;
                                    if lost(who+1) 
                                        lost(who+1) = 0;
                                        losttot = losttot -1;
                                    end
                                    who = who -1;
                                end
                            end
                        else  
                            if ~lost(who+1) & (losttot ~= nlost)
                                lost(who+1) = 1;
                                losttot = losttot + 1;
                                if pt(who+1) ~= st(who+1)-1
                                    ok(h(pt(who+1))) = 1;
                                end
                                if who == nnew -1
                                    ww = find( lost == 0);
                                    dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                    
                                    if dsq < mndisq
                                        minbonds = pt(ww);
                                        mndisq = dsq;
                                    end
                                else
                                    who = who + 1;
                                end
                            else
                                if pt(who+1) ~= st(who+1) -1
                                    ok(h(pt(who+1))) = 1;
                                end
                                pt(who+1) = st(who+1) -1;
                                if lost(who+1) 
                                    lost(who+1) = 0;
                                    losttot = losttot -1;
                                end
                                who = who -1;
                            end
                        end
                    end
                    
                    checkflag = checkflag + 1;
                    if checkflag == 1
                        plost = min([fix(mndisq/maxdisq) , (nnew -1)]);
                        if plost > nlost 
                            nlost = plost; 
                        else
                            checkflag = 2;
                        end
                    end
                    
                end  
                %   update resx using the minimum bond configuration               
                
                resx(ispan,labely(bonds(minbonds,2)+1)) = eyes(labelx(bonds(minbonds,1)+1));
                found(labelx(bonds(minbonds,1)+1)) = 1;
                %disp(size(resx))
                %disp('first')
                %disp(ispan)
            end

            %   THE PERMUTATION CODE ENDS
        end
        
        w = find(resx(ispan,:) >= 0);
        nww = length(w);
        
        if nww > 0 
            pos(w,:) = xyzs( resx(ispan,w) , 1:dim);
        end  %go back and add goodenough keyword thing   
        newguys = find(found == 0);
        
       nnew = length(newguys);
      
        if (nnew > 0) % & another keyword to workout inipos
%            disp(nnew);
%            disp(size(resx));
            newarr = zeros(zspan,nnew) -1;
            resx = [resx,newarr];
            %disp(size(resx)) 
            resx(ispan,n+1:end) = eyes(newguys);
            pos = [[pos];[xyzs(eyes(newguys),1:dim)]];
            nmem = zeros(nnew,1);
            mem = [mem;nmem];
            nun = 1:nnew;
            uniqid = [uniqid,((nun) + maxid)];
            maxid = maxid + nnew;
            % put in goodenough 
            n = n + nnew;
                
        end
 
    else
        ' Warning- No positions found for t='
    end
    w = find( resx(ispan,:) ~= -1);
    nok = length(w);
  %  nok
    if nok ~= 0
        mem(w) =0;
    end
    
    %r_eq_no = (resx(ispan,:) == -1);
    nadd = length(mem)-length(resx(ispan,:));
    if nadd ~= 0
        mem = [mem,zeros(nnew,1)];
    end    

    num = find(resx(ispan,:) == -1);
    if ~isempty(num)
        mem(num(:)) = 1;   
    end
 
    wlost = find(mem == memory_b+1);
    nlost = length(wlost);

    if nlost > 0 
        pos(wlost,:) = -maxdisp;
        % put in goodenough keyword stuff if 
    end
    if (ispan == zspan) | (i == z)
        nold = length(bigresx(1,:));
        nnew = n-nold;
        if nnew > 0
            newarr = zeros(z,nnew) -1;
            bigresx = [bigresx,newarr];
        end
        
        % again goodenough keyword

        %disp(strcat(num2str(i), ' of ' ,num2str(z), ' done.  Tracking  ',num2str(ntrk),' particles ', num2str(n),' tracks total'));
  
        bigresx(i-(ispan)+1:i,:) = resx(1:ispan,:);
        resx = zeros(zspan,n) - 1;
        %disp(size(resx))  
 
        wpull = find(pos(:,1) == -maxdisp);
        npull = length(wpull);
        
        if npull > 0
            lillist = zeros(1,2);
            for ipull=1:npull
                wpull2 = find(bigresx(:,wpull(ipull)) ~= -1);
                npull2 = length(wpull2);
                thing = [bigresx(wpull2,wpull(ipull)),zeros(npull2,1)+uniqid(wpull(ipull))];
                lillist = [lillist;thing];
                
            end
            olist = [[olist];[lillist(2:end,:)]];
             
        end
        %disp( olist);
        
        %if nkeep == 0 
        %    'Were going to crash now, no particles....'
        %end
        wkeep = find(pos(:,1) >= 0);
        nkeep = length(wkeep);
        
        resx = resx(:,wkeep);
        bigresx = bigresx(:,wkeep);
        pos = pos(wkeep,:);
        mem = mem(wkeep);
        uniqid = uniqid(wkeep);
        n = nkeep;
        dumphash = zeros(nkeep);
    end
   
end

wpull = find( pos(:,1) ~= -2*maxdisp);
npull = length(wpull);
if npull > 0
    lillist = zeros(1,2);
    for ipull=1:npull
        wpull2 = find(bigresx(:,wpull(ipull)) ~= -1);
        npull2 = length(wpull2);   
        thing = [bigresx(wpull2,wpull(ipull)),zeros(npull2,1)+uniqid(wpull(ipull))];
        lillist = [lillist;thing];
    end
    olist = [olist;lillist(2:end,:)];
end

olist = olist(2:end,:);
%bigresx = 0;
%resx = 0;

nolist = length(olist(:,1));
res = zeros(nolist,dd+1);
for j=1:dd
    res(:,j) = xyzs(olist(:,1),j);
end
res(:,dd+1) = olist(:,2);

info = uberize(res);