// Mean_vs_SD macro by Christophe Leterrier
// 02/02/2012
// Perform Photon Transfer Analysis as explained by Nico Stuurman here: http://valelab.ucsf.edu/~MM/MMwiki/index.php/Measuring_camera_specifications
// For each condition, take any number of images (2 according to the uManager website, but it can be more) at different exposure levels, and save it as a stack named XXXX.tif where XXXX is the exposure time (in ms padded to 4 digits with zeros)
// Point this macro to the folder containing all the acquired XXXX.tif images for this condition
// You will have to enter the offset in digital number units (as calculated from a dark stack of images), you can put zero to bypass offset correction
// It will output a table containing the image name, exposure time, mean, standard deviation and variance for each stack that you can save as a .txt or .xls.


macro "Mean_vs_SD" {

	InputDir = getDirectory("Please Choose Acquired Stacks Directory");
	Files = getFileList(InputDir);
	Dialog.create("Please Enter Offset Value in Digital Units");
	Dialog.addNumber("Offset", 1000, 0, 4, "ADU");
	Dialog.show();
	Offset = Dialog.getNumber();
	
	run("Set Measurements...", "area mean standard display redirect=None decimal=3");

	Expo = newArray(Files.length);
	Mean = newArray(Files.length);
	SD = newArray(Files.length);
	Var = newArray(Files.length);

setBatchMode(true);
	
	for (i=0; i<Files.length; i++) {

		// Gets the exposure time (from the name of the file)
		Expo[i] = parseInt(substring(Files[i], 0, lengthOf(Files[i])-4));
		
		// open the acquired stack
		open(InputDir + Files[i]);
		StackID = getImageID();
		// Sometimes the first slice is different from the rest
		run("Delete Slice");
		// subtract the Offset value
		run("Subtract...", "value=" + Offset + " stack");


		// Measure the mean intensity on a centered 256x256 region
		run("Z Project...", " projection=[Average Intensity]");
		run("Specify...", "width=256 height=256 x=128 y=128");
		// This is the "Mean"
		getStatistics(area, Mean[i], min, max, sd, histo);
		close();

		run("Duplicate...", "duplicate");
		DupID = getImageID();
		run("Delete Slice");
		selectImage(StackID);
		setSlice(nSlices);
		run("Delete Slice");
		// create a (nSlice-1) difference stack by substracting SLice (i) with Slice (i-1) 
		imageCalculator("Subtract create 32-bit stack", StackID, DupID);
		DiffID = getImageID();
		// get the SD of the whole difference stack
		run("Z Project...", "start=1 stop=98 projection=[Standard Deviation]");
		// Measure the Mean of this image as the SD of the difference stack
		// this is the "SD"
		getStatistics(area, SD[i], min, max, sd, histo);
		close();
		
		selectImage(DiffID);
		close();
		selectImage(DupID);
		close();
		selectImage(StackID);
		close();

		// Calculates the variance as square of SD
		Var[i] = SD[i] * SD[i];
	}
	
	// Generate the results table
	title1= "Mean vs SD Results";
	title2 = "["+title1+"]";
	f = title2;
	if (isOpen(title1))
		print(f, "\\Clear");
	else
		run("New... ", "name="+title2+" type=Table");
	Headings="\\Headings:n\tImage Name\tExpo Time (ms)\tMean Intensity\tStandard Deviation\tVariance";
	print(f,Headings);
	for (i=0; i<Files.length; i++) {
		ResultsLine=d2s(i+1,0)+"\t"+Files[i]+"\t"+Expo[i]+"\t"+Mean[i]+"\t"+SD[i]+"\t"+Var[i];
		print(f,ResultsLine);
	}

setBatchMode("exit and display");
}