//Example macro that does a diffusion analysis on a results table of tracking 
//data. The required format is as generated by the Manual Tracking plugin 
//http://rsbweb.nih.gov/ij/plugins/track/track.html
//assuming Results table column headings: Track > Slice > X > Y
//It requires many tracks and many time points!
//The analysis uses the time ensemble average method described in Charsooghi, MA et al 2011
//http://www.sciencedirect.com/science/article/pii/S0010465510003620

//Get number of tracks (nTracks)
nTracks = 0;
for (a=0; a<nResults(); a++) {
    if (getResult("Track",a)>nTracks)
    {
     nTracks = getResult("Track",a);
    	}
    	else{};
}

//Find the track length for each track - write track length to results table
//some variables
Track=1;
L_Track=0;

//work though tracks and determine length
for (i=0; i<nResults(); i++){
	if (getResult("Track", i)==Track) {L_Track++;} else {

//Do this if its the first track
	if(Track==1){		for (j=(i-L_Track); j<i; j++){
				setResult("T_Length", j, L_Track);			
		}
				L_Track=0;
				Track++;
				
//Do this if its in the middle
	L_Track=1;	
	} else {for (j=(i-L_Track); j<i; j++){
				setResult("T_Length", j, L_Track);
		}
				L_Track=1;
				Track++;
			}
		}
	}
//Do this to get the last track
	L_Track=0;
	for (j=0; j<nResults; j++) {
		if (getResult("Track", j)==nTracks) {L_Track++;}
	}

	for (k=0; k<nResults; k++) {
		if (getResult("Track", k)==nTracks) {setResult("T_Length", k, L_Track);
	}
}
updateResults();

//get last slice
maxslice = 0;
for (b=0; b<nResults(); b++) {
    if (getResult("Slice",b)>maxslice)
    {
     maxslice = getResult("Slice",b);
    	}
    	else{};
}

//get first slice
minslice = maxslice;
for (c=0; c<nResults(); c++) {
    if (getResult("Slice",c)<minslice)
    {
     minslice = getResult("Slice",c);
    	}
    	else{};
}

//The window sizes for analysis range from 1-step to - maxslice-1
//Calculate squared dispalcement from tracking data for all possible window sizes 

MSD = newArray();
time = newArray();
divide=0;
r_total=0;
dis2 = 0;

for (u=1; u<(maxslice); u++) {

for (i=0; i<nResults(); i++){
	if (getResult("Slice", i)<=u) {}
	
	else{ if (getResult("Track", i)>getResult("Track", i-u)) {}
	
	else { if (getResult("T_Length", i)>=u && getResult("Track", i-u)==getResult("Track", i)) {
	B9 = getResult("X", i);
	B8 = getResult("X", i-u);
	C9 = getResult("Y", i);
	C8 = getResult("Y", i-u);
	disx2 = (B9-B8)*(B9-B8);
	disy2 = (C9-C8)*(C9-C8);
	dis2 = (disx2 + disy2);
	r_total = r_total+dis2;	
	divide++;
			}	
		}
	}
}
time = Array.concat(time, u);	
MSD = Array.concat(MSD, (r_total)/divide);
r_total=0;
divide=0;
}

Fit.doFit("Straight Line", time, MSD);
intercept = d2s(Fit.p(0),6);
slope = d2s(Fit.p(1),6);
r2 = d2s(Fit.rSquared,3);
print("slope = "+slope);
print("intercept = "+intercept);
print("R^2 = "+r2);
Fit.plot();

//Richard Mort 26/01/2013