/*This code performs the normalized median test on ImageJ PIV output*/


function median(a){
	Array.sort(a);
	med = 0.0;
	if(a.length%2 == 0){
		med = (0 + a[(a.length/2 - 1) ] + a[(a.length/2)])/2;

	}else{
		med = 0 + a[a.length/2 -1];
	}
	
	return med;
}
		

function setElement(myArray, height, row, col, value){
//hashes 2-D array values in a 1-D array in row major order 

	myArray[col * height + row] = (0 + value);

}

function getElement(myArray, height, row, col){
//returns values from the hashed 1-D array
	return (0 + myArray[col * height + row]);


}


			
Dialog.create("Normalized Median Test");
Dialog.addNumber("eta (Noise Parameter)", 0.2);
Dialog.addNumber("threshold", 2);
Dialog.addNumber("magnitude threshold (lower bound)", 0);
Dialog.addNumber("magnitude threshold (upper bound)", 99999);
Dialog.addNumber("Video Frames Per Second", 2000);
Dialog.addNumber("Pixels per mm", 475);
Dialog.addNumber("frames", 2000);
Dialog.show;

eta = Dialog.getNumber();
rthresh = Dialog.getNumber();
magThreshLower = Dialog.getNumber();
magThreshUpper = Dialog.getNumber();
fps = Dialog.getNumber();
pixelspermm = Dialog.getNumber();
slices = Dialog.getNumber();
replacedVectors = 0;
totalVectors = 0;

path = getDirectory("select a folder with Multiple PIV Results");

for(z=1; z <=slices; z++){
	call("java.lang.System.gc");
	dataPath = File.openAsString(path + "seq_" + z + "_PIV2_disp.txt");
	outputStr = "";
	rows = split(dataPath, "\n");
	cols = split(rows[0], " ");
	temp = split(rows[rows.length-1], " ");

	xlength = 0 + temp[0];
	ylength = 0 + temp[1];

	//initialize 2-D arrays 
	u = newArray((xlength+1) * (ylength+1));
	Array.fill(u, 0.0);
	v = newArray((xlength+1) * (ylength+1));
	Array.fill(v, 0.0);
	height = ylength;
	width = xlength;


	xcoords = newArray(rows.length);
	Array.fill(xcoords, 0);		
	ycoords = newArray(rows.length);
	Array.fill(ycoords, 0);


	for(i=0; i< rows.length; i++){

		temp = split(rows[i], " ");
		xcoords[i] = temp[0];
		ycoords[i] = temp[1];

		setElement(u, height, temp[0], temp[1], temp[2]);
		setElement(v, height, temp[0], temp[1], temp[3]);

	}





	step = 0 + xcoords[1] - xcoords[0];
	

	uNeighbs = newArray(8);
	vNeighbs =newArray(8);
	uResiduals = newArray(8);
	vResiduals = newArray(8);

	//
	for(x = xcoords[0] ; x<=xlength; x+= step){
		for(y = ycoords[0]; y<=ylength; y+=step){

			utmp = getElement(u,height,x,y);
			//if vectors are on edges, set the appropriate neighbors to zero
			if(x-step <= 0){
				uNeighbs[0] = 0;
				vNeighbs[0] = 0;
				uNeighbs[3] = 0;
				vNeighbs[3] = 0;
				uNeighbs[5] = 0;
				vNeighbs[5] = 0;
			}

			if(y - step <= 0){
				uNeighbs[0] = 0;
				vNeighbs[0] = 0;
				uNeighbs[1] = 0;
				vNeighbs[1] = 0;
				uNeighbs[2] = 0;
				vNeighbs[2] = 0;
			}
		
			if(x+step >= xlength){
				uNeighbs[2] = 0;
				vNeighbs[2] = 0;
				uNeighbs[4] = 0;
				vNeighbs[4] = 0;
				uNeighbs[7] = 0;
				vNeighbs[7] = 0;
			}

			if(y + step >= ylength){
				uNeighbs[5] = 0;
				vNeighbs[5] = 0;
				uNeighbs[6] = 0;
				vNeighbs[6] = 0;
				uNeighbs[7] = 0;
				vNeighbs[7] = 0;
			}

			//if vector is in the middle, populate an array of all its neighbors
			if(x-step > 0 && y-step > 0 && x+step < xlength && y+step < ylength){
				uNeighbs[0] = 0 + getElement(u,height, x-step, y-step);
				uNeighbs[1] = 0 + getElement(u,height,x, y-step);
				uNeighbs[2] = 0 + getElement(u,height,x+step, y-step);
				uNeighbs[3] = 0 + getElement(u,height,x - step, y);
				//print("getElement (" + (x + step) + ", " + (y+step) + "): " + getElement(u,height, x + step, y));
				uNeighbs[4] = 0 + getElement(u,height, x + step, y);
				uNeighbs[5] = 0 + getElement(u,height,x-step, y+step);
				uNeighbs[6] = 0 + getElement(u,height,x, y+step);
				uNeighbs[7] = 0 + getElement(u,height,x+step, y+step);
				uMedian = 0 + median(uNeighbs);

			}


			vtmp = getElement(v,height,x,y);

				vNeighbs[0] = 0 + getElement(v,height, x-step, y-step);
				vNeighbs[1] = 0 + getElement(v,height,x, y-step);
				vNeighbs[2] = 0 + getElement(v,height,x+step, y-step);
				vNeighbs[3] = 0 + getElement(v,height,x - step, y);
				vNeighbs[4] = 0 + getElement(v,height,x+step, y);
				vNeighbs[5] = 0 + getElement(v,height,x-step, y+step);
				vNeighbs[6] = 0 + getElement(v,height,x, y+step);
				vNeighbs[7] = 0 + getElement(v,height,x+step, y+step);
				vMedian = 0 + median(vNeighbs);

			for(p=0; p<8; p++){
				//calculate array of residuals
				uResiduals[p] = abs(uNeighbs[p] - uMedian);
				vResiduals[p] = abs(vNeighbs[p] - vMedian);
			}

			if( abs(utmp - uMedian)/(median(uResiduals) + eta) > rthresh || abs(vtmp - vMedian)/(median(vResiduals) + eta) > rthresh){
				//perform normalized median test

				setElement(u,height, x,y, 0);
				setElement(v,height, x,y, 0);

				//Uncomment the following four lines to replace the vector with a bilinear interpolation
				//uinterp = 0.25* (uNeighbs[0] + uNeighbs[2] + uNeighbs[5] + uNeighbs[7]);	uncomment this line to 
				//setElement(u,height, x,y, uinterp);
				//vinterp = 0.25* (vNeighbs[0] + vNeighbs[2] + vNeighbs[5] + vNeighbs[7]);
				//setElement(v,height, x,y, vinterp);

				replacedVectors = replacedVectors + 1;

			}
		}
	}

	for(j=0; j< rows.length; j++){
		//selectWindow("U");
		uout = getElement(u,height, xcoords[j], ycoords[j]);

		//selectWindow("V");
		vout = getElement(v,height, xcoords[j], ycoords[j]);

		if(uout != 0 || vout != 0){
			totalVectors = totalVectors + 1;
		}

		outputStr += xcoords[j] + " " + ycoords[j] + " " + 0+uout*fps/pixelspermm + " " + 0+vout*fps/pixelspermm + " " + "\n";

	}

	File.saveString(outputStr, path + "\\NMT" + z + ".txt");
	print("Iteration " + z + " complete");
	call("java.lang.System.gc"); //garbage collection
}

outputStr2 = '';

//open each NMT'ed frame and perform a temporal average
for(q=1;q<slices; q++){
	call("java.lang.System.gc"); //garbage collection
	NMTdata = File.openAsString(path + "NMT" + q + ".txt");
	NMTrows = split(NMTdata, "\n");

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

	}

	for(h=0; h<NMTrows.length; h++){
			cols = split(NMTrows[h], " ");
			

			
			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[h]++;
				yZeros[h]++;
			}
			
			if(mag != 0){			
				xdisp[h] = xdisp[h] + cols[2];
				ydisp[h] = ydisp[h] + cols[3];
			}


		if(q == slices-1){


			if(i-xZeros[h] !=0 && q != xZeros[h]){
				//if there are zeros for some, but not all vectors at this location throughout the stack then only average non-zero vectors
				xdisp[h] = xdisp[h]/abs(q-xZeros[h]);
			} else{
				xdisp[h] = 0;
			}

			if(i-yZeros[h]!= 0 && q != yZeros[h]){
				//if there are zeros for some, but not all vectors at this location throughout the stack then only average non-zero vectors
				ydisp[h] = ydisp[h]/abs(q-yZeros[h]);
			} else{
				ydisp[h] = 0;
			}
			NMTxcoords[h] = cols[0];
			NMTycoords[h] = cols[1];

			mags[h] = sqrt( xdisp[h]*xdisp[h] + ydisp[h]*ydisp[h] );

			outputStr2 += NMTxcoords[h] + " " + NMTycoords[h] + " " + xdisp[h] + " " + ydisp[h] + " " + mags[h] + " " + "\n";


		}
	}



}


File.saveString(outputStr2, path + "\\NMTaveraged(No Interpolation).txt");
//uncomment to report vector calidation statistics
//File.saveString("Number of Replaced Vectors:" + replacedVectors + "\n Number of Total Vectors: " + totalVectors + " Percent Replaced: " + replacedVectors/totalVectors, path + "\\replacedVectors.txt");

print("Done!!!!");