/*This code runs PIV on an image sequence and performs a temporal average of each velocity vector*/

macro "Run PIV on Stack [f1]"{

id0 = getImageID();
slices = nSlices;
outputStr = "";
xcomp = 0;
ycomp = 0;
mag = 0;


//creates dialog boxes to enter important parameters, these may vary by experimental setup
Dialog.create("Multi-slice PIV");
Dialog.addNumber("Interrogation window 1", 256);
Dialog.addNumber("search window 1", 512);
Dialog.addNumber("vector spacing 1", 128);
Dialog.addNumber("Interrogation window 2", 128);
Dialog.addNumber("search window 2", 256);
Dialog.addNumber("vector spacing 2", 64);
Dialog.addNumber("Interrogation window 3", 0);
Dialog.addNumber("search window 3", 0);
Dialog.addNumber("vector spacing 3", 0);
Dialog.addNumber("correlation threshold", 0.8);
Dialog.addNumber("Video Frames Per Second", 2000);
Dialog.addNumber("Pixels per mm", 475);
Dialog.show;

piv1 = Dialog.getNumber();
sw1 = Dialog.getNumber();
vs1 = Dialog.getNumber();
piv2 = Dialog.getNumber();
sw2 = Dialog.getNumber();
vs2 = Dialog.getNumber();
piv3 = Dialog.getNumber();
sw3 = Dialog.getNumber();
vs3 = Dialog.getNumber();
corr = Dialog.getNumber();
fps = Dialog.getNumber();
pixelspermm = Dialog.getNumber();

path = getDirectory("select a folder to save PIV results");

for(s=1;s<slices;s++){
	selectImage(id0);
	call("java.lang.System.gc");
	run("Duplicate...", "title=[seq_"+s+"] duplicate range="+s+"-"+s+1+"");
	run("iterative PIV(Advanced)...", "  piv1="+piv1+" sw1="+sw1+" vs1="+vs1+" piv2="+piv2+" sw2="+sw2+" vs2="+vs2+" piv3="+piv3+" sw3="+sw3+" vs3="+vs3+" correlation="+corr+" batch path=["+path+"]");
	for(p=1; p<=2; p++){
		selectWindow("seq_"+ s + "_PIV" + p);
		close();
	}
		selectWindow("seq_"+ s);
		close();
}


for(i=1;i<slices; i++){
	
	dataPath = File.openAsString(path + "seq_"+i+"_PIV2_disp.txt");
	rows = split(dataPath, "\n");

	//On first pass initialize displacement arrays to be same size as PIV data arrays
	if(i ==1){
		print ("outputting average...");
		xcoords = newArray(rows.length);
		Array.fill(xcoords, 0);		
		ycoords = newArray(rows.length);
		Array.fill(ycoords, 0);
		mags = newArray(rows.length);
		Array.fill(mags,0);
		xdisp = newArray(rows.length);
		Array.fill(xdisp, 0);
		ydisp = newArray(rows.length);
		Array.fill(ydisp, 0);
		xZeros = newArray(rows.length);
		Array.fill(xZeros, 0);		
		yZeros = newArray(rows.length);
		Array.fill(yZeros,0);

	}
	
	for(j=0; j<rows.length; j++){
			cols = split(rows[j], " ");

			xcomp = 0 + cols[2];
			ycomp = 0 + cols[3];
			mag = sqrt(xcomp*xcomp + ycomp*ycomp);

			//if no velocity vector at point, exclude it from the noise averaging
 			if(mag == 0){
				xZeros[j]++;
				yZeros[j]++;
			}else{	
			//otherwise, sum in time to be averaged later		
				xdisp[j] = xdisp[j] + cols[2];
				ydisp[j] = ydisp[j] + cols[3];
			}


		if(i == slices-1){
			//on last pass divide summed velocities by number of frames to average
			if(i-xZeros[j] !=0){
				xdisp[j] = xdisp[j]/abs(i-xZeros[j]) * fps/pixelspermm;
			} else{
				xdisp[j] = 0;
			}

			if(i-yZeros[j]!= 0){
				ydisp[j] = ydisp[j]/abs(i-yZeros[j]) * fps/pixelspermm;
			} else{
				ydisp[j] = 0;
			}
			xcoords[j] = cols[0];
			ycoords[j] = cols[1];

			mags[j] = sqrt( xdisp[j]*xdisp[j] + ydisp[j]*ydisp[j] );
			outputStr += xcoords[j] + " " + ycoords[j] + " " + xdisp[j] + " " + ydisp[j] + " " + mags[j] + " " + "\n";
			
		}
	}



}


File.saveString(outputStr, path + "\\averaged.txt");
print("Done!!!!");


function closeAllWindows(){

	for(i=1; i< nImages; i++){
		selectImage(i);
		if(getImageID()!=id0){
			close();
		}
	}
}