//
// Displacement vector analysis
//
// ver 2019.12.09 Ken Takahashi
//

// Get folder list
folderList = File.openAsString("C:/vector_analysis/movie_temp/joblist.txt");
folderNames = split(folderList,"\n");
folderNo = folderNames.length;
print(folderNo)

// Call analysis function and repeat for all folders
for (i=0; i<folderNo; i++) {
	print("***********************************");
	print("Analyzing "+folderNames[i]);
	print("***********************************");
	print("");
	vector_analysis(folderNames[i]);
}

print("Vector analyses of all folders completed!");



//
// Analysis function
//
function vector_analysis(folderName) {
	// Load AVI video
	open(folderName+"/analyze.avi");
	frameNo = nSlices();
	saveAs("PNG", folderName+"/Phase_contrast.png");
	run("Close");
	print("Number of frames in analyze.avi: "+frameNo);

	// Run PIV analysis
	for (i=2; i<frameNo+1; i++) {
		run("AVI...", "open="+folderName+"/analyze.avi first=1 last=1 use convert");
		run("AVI...", "open="+folderName+"/analyze.avi first="+i+" last="+i+" use convert");
		run("Images to Stack", "name=analyze.avi title=[] use");
		print("PIV analysis between frames 1 and "+i);
		run("iterative PIV(Advanced)...", "  piv1=128 sw1=128 vs1=64 piv2=64 sw2=64 vs2=32 piv3=32 sw3=32 vs3=16 correlation=0.60 debug_x=-1 debug_y=-1 path=/ what=[Accept this PIV and output] noise=0.20 threshold=5 c1=3 c2=1 save="+folderName+"/vec_"+i-1+".txt");
	
		selectWindow("analyze.avi_PIV1");
		close();
		selectWindow("analyze.avi_PIV2");
		close();
		selectWindow("analyze.avi_PIV3");
		close();
		selectWindow("analyze.avi");
		close();
		if ( isOpen("Scale Graph") ) { 
			selectWindow("Scale Graph"); 
			run("Close"); 
		} 
	}

	// find maximum vector size at each pixel
	for (i=1; i<frameNo; i++) {
		run("Text Image... ", "open="+folderName+"/vec_"+i+".txt");
	}

	run("Images to Stack", "name=Stack title=[] use");
	run("Reslice [/]...", "output=1.000 start=Left avoid");
	setSlice(5);
	selectWindow("Stack");
	close();
	run("Green Fire Blue");
	run("Enhance Contrast", "saturated=0.35");

	vectorNo = getWidth();
	traj1 = newArray(vectorNo);

	for (k=0; k<vectorNo; k++){	
		setSlice(5);
		for (i=0; i<frameNo; i++) {
			traj1[i] = getPixel(k, i);
		}
		Array.getStatistics(traj1, min, max, mean, dummy);

	    for (j=0; j<vectorNo; j++) { 
    	    if (traj1[j]==max) { 
        	    index=j; 
	        } 
    	}

		// set max vector size at column 5
		setPixel(k,0,max);

		// set x vector value at column 3
		setSlice(3);
		max_x = getPixel(k,index);
		setPixel(k,0,max_x);

		// set y vector value at column 4
		setSlice(4);
		max_y = getPixel(k,index);
		setPixel(k,0,max_y);
	}

	// Calculate mean vector size and write
	setSlice(5);
	for (k=0; k<vectorNo; k++){	
		traj1[k] = getPixel(k, 0);
	}
	Array.getStatistics(traj1, min, max, mean, dummy);
	setSlice(7);
	setPixel(0,0,mean);

	// Save max vector image
	run("Reslice [/]...", "output=1.000 start=Top rotate avoid");
	run("Slice Remover", "first=2 last="+frameNo+" increment=1");
	saveAs("Text Image", folderName+"/Max_vector.txt");

	// Display max vector map and save as PNG file
	run("plot...", "select="+folderName+"/Max_vector.txt vector_scale=2 max=3 plot_width=0 plot_height=0 show lut=S_Pet");
	selectWindow("Vector plot_Max_vector.txt");
	saveAs("PNG", folderName+"/Max_vector.png");

	// Create overlay image and save
	open(folderName+"/Phase_contrast.png");
	run("RGB Color");
	run("Calculator Plus", "i1=Max_vector.png i2=Phase_contrast.png operation=[Add: i2 = (i1+i2) x k1 + k2] k1=1 k2=0");
	selectWindow("Phase_contrast.png");
	saveAs("PNG", folderName+"/Overlaid.png");

	print("Max vector map generated successfully!");
	print("");
	print("");
	print("");

	close();
	selectWindow("Max_vector.png");
	close();
	selectWindow("Scale Graph");
	close();
	selectWindow("Reslice of Reslice");
	close();
	selectWindow("Reslice of Stack");
	close();
}