function [tracks, fitP] = track(xyzs,maxdisp,param, tempfitP)

% ; xyzs=PList(:,1:3);maxdisp = param.stepLength;
% ; 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, param )
% ;  set all keywords in the space below
% ; 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. must be positve
% ;  positionlist(d,*): contains the time t that the position 
% ;     was determined, must be integers (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)
%  OPTIONAL INPUT:
%   param:  a structure containing a few tracking parameters that are
%       needed for many applications.  If param is not included in the
%       function call, then default values are used.  If you set one value
%       make sure you set them all:
% ;         param.mem: 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.
% ;         param.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:param.dim-1,*)
% ;             and the extra data in positionlist(param.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.
% ;         param.good: set this keyword to eliminate all trajectories with
% ;             fewer than param.good valid positions.  This is useful
% ;             for eliminating very short, mostly 'lost' trajectories
% ;             due to blinking 'noise' particles in the data stream.
%;          param.quiet: set this keyword to 1 if you don't want any text
% ; 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 .
% ; 
% ; 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.
% ;  1/05 Transmuted to MATLAB by D. Blair
% ;  5/05  ERD Added the param structure to simplify calling.
%    6/05  ERD Added quiet to param structure
%    7/05  DLB Fixed slight bug in trivial bond code
%    3/07  DLB Fixed bug with max disp pointed out by Helene Delanoe-Ayari
%
% ; 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.
% ;
% ;-

dd = length(xyzs(1,:));

%use default parameters if none given
if nargin==2
    %default values
    memory_b=0; % if mem is not needed set to zero
    goodenough = 0;  % if goodenough is not wanted set to zero
    dim = dd - 1;
    quiet=0;
else
    memory_b    =   param.mem;
    goodenough  =   param.good;
    dim         =   param.dim;
    quiet       =   param.quiet;
end


% checking the input time vector
t = xyzs(:,dd);
st = circshift(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;

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);
% implanting unq directly
    indices = find(t ~= circshift(t,-1));
        count = length(indices);
        if count > 0
            res = indices;
        else  
            res = length(t) -1;
        end
 %%%%%%%%%%%%%%%%%%%%%%%       
        
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.];

if goodenough > 0 
    dumphash = zeros(n,1);
    nvalid = ones(n,1);
end

%  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);



%   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   
        
            %   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); %number of old positions ew
            if ntrack == 0 
                'There are no valid particles to track idiot!'
                break
            end
            xmat = zeros(ntrack,m);
            count = 0;
            for kk=1:ntrack %old ew
                for ll=1:m %new ew
                    xmat(kk,ll) = count;
                    count = count+1;
                end
            end
            % array: xmat: durchzaehlen: 1 2 3; 4 5 6; 7 8 9 ... (n zeilen, m spalten)
            count = 0;
            for kk=1:m
                for ll=1:ntrack
                    ymat(kk,ll) = count;
                    count = count+1;
                end
            end
            % array: xmat: durchzaehlen: 1 2 3 4; 5 6 7 8; 9 ... (m zeilen, n spalten)

            xmat = (mod(xmat,m) + 1); % 1 2 3; 1 2 3; 1 2 3;...
            ymat = (mod(ymat,ntrack) +1)'; % 1 1 1 1; 2 2 2 2; 3 3 3 3;...
            [lenxn,lenxm] = size(xmat);
%            whos ymat
%            whos xmat
%            disp(m)

            for d=1:dim
                x = xyi(:,d);  % new ew
                y = pos(wh,d); % old ew
                xm = x(xmat);
                ym = y(ymat(1:lenxn,1:lenxm));
                if size(xm) ~= size(ym)
                    xm = xm';
                end
                
                if d == 1
                    dq = (xm -ym).^2;
                    %dq = (x(xmat)-y(ymat(1:lenxn,1:lenxm))).^2;
                else
                    dq = dq + (xm-ym).^2;
                    %dq = dq + (x(xmat)-y(ymat(1:lenxn,1:lenxm)) ).^2;
                end
            end
            % dq: array with all distances ew
            % logisches array mit 1 ueberall wo distances kleiner als max
            % dist --> moegliche bonds
            ltmax = dq < maxdisq;
            
            % figure out which trivial bonds go with which
            
            rowtot = zeros(n,1); % Spaltenvektor in Laenge der Anzahl der alten Pos
            rowtot(wh) = sum(ltmax,2); % ich verstehe: rowtot speichert die Anzahl der moeglichen bonds per alter Position
            
            
            if ntrack > 1 % ntrack: Anzahl der alten Positionen 
                coltot = sum(ltmax,1);
            else
                coltot = ltmax;
            end
            which1 = zeros(n,1);
            for j=1:ntrack %loop ueber alle alten Positionen
                [mx, w] = max(ltmax(j,:)); 
                which1(wh(j)) = w;
            end
            
            ntrk = fix( n - sum(rowtot == 0));
            w= find( rowtot == 1) ;
            ngood = length(w);
            if ngood ~= 0
                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
        %THE TRIVIAL BOND CODE ENDS
        
        w = find(resx(ispan,:) >= 0);
        nww = length(w);
        
        if nww > 0 
            pos(w,:) = xyzs( resx(ispan,w) , 1:dim);
            if goodenough > 0 
                nvalid(w) = nvalid(w) + 1;
            end
        end  %go back and add goodenough keyword thing   
        newguys = find(found == 0);
        
       nnew = length(newguys);
      
        if (nnew > 0) % & another keyword to workout inipos
            newarr = zeros(zspan,nnew) -1;
            resx = [resx,newarr];

            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;
            if goodenough > 0 
                dumphash = [dumphash;zeros(1,nnew)'];
                nvalid = [nvalid;zeros(1,nnew)'+1];
            end
            % put in goodenough 
            n = n + nnew;
                
        end
    end
 %------------------------------------------------------------------------------

    w = find( resx(ispan,:) ~= -1);
    nok = length(w);
    if nok ~= 0
        mem(w) =0;
    end
    
    mem = mem + (resx(ispan,:)' == -1);
    wlost = find(mem == memory_b+1);
    nlost =length(wlost);

    if nlost > 0 
        pos(wlost,:) = -maxdisp;
        if goodenough > 0
            wdump = find(nvalid(wlost) < goodenough);
            ndump = length(wdump);
            if ndump > 0
                dumphash(wlost(wdump)) = 1;
            end
        end
        % 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
        if goodenough > 0  
            if (sum(dumphash)) > 0
                wkeep = find(dumphash == 0);
                nkeep = length(wkeep);
                resx = resx(:,wkeep);
                bigresx = bigresx(:,wkeep);
                pos = pos(wkeep,:);
                mem = mem(wkeep);
                uniqid = uniqid(wkeep);
                nvalid = nvalid(wkeep);
                n = nkeep;
                dumphash = zeros(nkeep,1);
            end
        end
        
        % again goodenough keyword
        if quiet~=1
            disp(strcat(num2str(i), ' of ' ,num2str(z), ' done.  Tracking  ',num2str(ntrk),' particles  ', num2str(n),' tracks total'));
        end
        bigresx(i-(ispan)+1:i,:) = resx(1:ispan,:);
        resx = zeros(zspan,n) - 1;

 
        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

        
 
        wkeep = find(pos(:,1) >= 0);
        nkeep = length(wkeep);
        if nkeep == 0 
                'Were going to crash now, no particles....'
        end
        resx = resx(:,wkeep);
        bigresx = bigresx(:,wkeep);
        pos = pos(wkeep,:);
        mem = mem(wkeep);
        uniqid = uniqid(wkeep);
        n = nkeep;
        dumphash = zeros(nkeep,1);
        if goodenough > 0
            nvalid = nvalid(wkeep);
        end
    end
   
end

if goodenough > 0 
    nvalid = sum(bigresx >= 0 ,1);
    wkeep = find(nvalid >= goodenough);
    nkeep = length(wkeep);
    if nkeep == 0
        for i=1:10
        disp('You are not going any further, check your params and data')
        end
        disp('the code broke at line 1045')
        %return
    end
    if nkeep < n
        bigresx = bigresx(:,wkeep);
        n = nkeep;
        uniqid = uniqid(wkeep);
        pos = pos(wkeep,:);
    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);

fitP = tempfitP(olist(:,1),:);

% this is uberize included for simplicity of a single monolithic code

ndat=length(res(1,:));
newtracks=res;

%u=unq(newtracks(:,ndat));

% inserting unq
indices = find(newtracks(:,ndat) ~= circshift(newtracks(:,ndat),-1));
        count = length(indices);
        if count > 0
            u = indices;
        else  
            u = length(newtracks(:,ndat)) -1;
        end


ntracks=length(u);
u=[0;u];
for i=2:ntracks+1
    newtracks(u(i-1)+1:u(i),ndat) = i-1;
end

% end of uberize code

tracks = newtracks;



 
