// @string(label = "XY pixels", value = 256,persist=false) cubeDim
// @string(label = "Pixel Size (um)", value = 0.1,persist=false) pixSize
// @string(label = "Number of Features", min=1,value = 5,persist=false) numFeats
// @string(label = "Feature Size (pixels)", value = 5,persist=false) featSize
// @string(label = "Number of frames", value = 100,persist=false) timeFrames
// @string(label = "Interval (s)", value = 0.1,persist=false) dT
// @string(label = "Diffusion Coefficient (um^2/s)", value = 20,persist=false) Diffusion
// @Boolean(label = "Avoid Edges?", description = "If checked, will not create spots within 10% of the edge", value = true,persist=false) avoidEdges
//
//
//									- Dave Mason [dnmason@liv.ac.uk]
//									- Centre for Cell Imaging [http://cci.liv.ac.uk]
//--------------------------------------------------------------------------------
// CREATE BROWNIAN DATA FOR TEACHING TRACKING 
//
setBatchMode(true);
run("Options...", "iterations=1 count=1 black do=Nothing");

if (1==2){
//-- Set Cube dimension in pixels
cubeDim=512;
//-- What is the pixel size?
pixSize=0.1; //--um/pixel
//-- Number of features per slice:
numFeats=15;
//-- How big do you want the features in pixels?
featSize=5;
//-- How many time frames?
timeFrames=500;

//-- Movement Parameters
//D=0.003; // um^2/s
D=0.1; // pixels^2/s
D=D/(pixSize*pixSize); // um^2/s
dT=0.05; // s
//dT=1; // s
}
cubeDim=parseFloat(cubeDim);
pixSize=parseFloat(pixSize);
numFeats=parseFloat(numFeats);
featSize=parseFloat(featSize);
timeFrames=parseFloat(timeFrames);
dT=parseFloat(dT);
D=parseFloat(D);

k=sqrt(2*D*dT);


//-- Make a new stack with identical x, y dimensions
newImage("Tracks", "8-bit grayscale-mode", cubeDim, cubeDim, 1, 1, timeFrames);
//-- Set foreground colour so points are drawn in white
setForegroundColor(255, 255, 255);


//-- Make an array to hold the values
sX=newArray(timeFrames+1);
sY=newArray(timeFrames+1);

//-- For each particle make starting coordinates
for (i=1;i<numFeats+1;i++){
setSlice(1);	


//-- Set starting coordinates and optionally bias middle of frame
if (avoidEdges==true){
sX[1]=floor(random*(cubeDim*0.8))+(cubeDim*0.1);
sY[1]=floor(random*(cubeDim*0.8))+(cubeDim*0.1);
}else{
sX[1]=floor(random*cubeDim);
sY[1]=floor(random*cubeDim);
}

//-- Draw in starting positions
makeOval(sX[1], sY[1], featSize, featSize);
//-- Draw the point into the image (using foreground colour)
run("Draw", "slice");
//-- Deselect everything
run("Select None");

//-- Make the track
for (j=2;j<timeFrames+1;j++){
	
	sX[j]=sX[j-1]+(k*randn());
	sY[j]=sY[j-1]+(k*randn());

//-- Draw in features
setSlice(j);
//makeOval(round(sX[j]), round(sY[j]), featSize, featSize);
makeOval(sX[j], sY[j], featSize, featSize);
		//-- Draw the point into the image (using foreground colour)
		run("Draw", "slice");
		//-- Deselect everything
		run("Select None");
		
} // j loop (timepoints)
//-- Fill the holes
run("Fill Holes", "stack");	


}
 // i loop (features)


setSlice(1);

//-- Do some post processing

run("Gaussian Blur...", "sigma="+(featSize/2)+" stack");
run("Add Specified Noise...", "stack standard=10");

//-- Set the calibration to match the Motion Parameters
run("Properties...", "unit=µm pixel_width="+pixSize+" pixel_height="+pixSize+" voxel_depth=1.0000000 frame=["+dT+" sec]");

//-- Print a summary to the log
print("----------------------------------");
print("Image Properties:");
print("----------------------------------");
print("Dimensions: "+cubeDim);
print("Pixel size: "+pixSize);
print("Timeframes: "+timeFrames);
print("Interval: "+dT+" s");
print("----------------------------------");
print(" Trajectory Properties:");
print("----------------------------------");
print("Creating "+numFeats+" features");
print("Feature size: "+featSize+" um");
print("Diffusion Coefficient: "+D+" um^2/s");
print("----------------------------------");

setBatchMode(false);


//-- Functions

function randn() {
	//-- Use the Box-Muller method to produce normally distributed random numbers
	U=random();
	V=random();
	return sqrt(-2*log(U))*cos(2*3.14*V);
}