// Calib_EMGain macro by Christophe Leterrier
// 02/02/2012
// automate the EMCCD EM gain measurement explained at http://www.photometrics.com/resources/learningzone/calculating-electron-multiplication-gain.php
// keep the illumination constant when acquiring all images !
// for each gain/exposure time couple, acquire two images:
// - a "bright" images named gXXXX_YYYY.tif with XXXX = gain value (padded to 4 digits with zeros) and YYYY = exposition value (in ms padded to 4 digits with zeros)
// - a "dark" image with no illumination light names gXXXX_dark.tif (same exposure time).
// The first couple will be "no gain" i.e. gain = 1 g0001_YYYY.tif and g0000_dark.tif (a point at EM gain = zero is necessary)
// just point the macro to the folder containing the images (and no other files!) and it will output a table with the theoretical gain (slider value) vs experimental gain (measured).
// Then you can save the table as a txt or xls file.

macro "Calib_EMGain" {
	// Choose a folder
	InputDir = getDirectory("Please Choose Acquired Stacks Directory");
	Files = getFileList(InputDir);

	// initialize results arrays
	nCond = Files.length/2;
	GainsTheo = newArray(nCond);
	Expo = newArray(nCond);
	BrightMean = newArray(nCond);
	DarkMean = newArray(nCond);
	GainsExp = newArray(nCond);

setBatchMode(true);
	
	// loops on all files in the folder (only tif files from the experiment!)
	for (i=0; i<Files.length; i++) {
		// if the file corresponds to a bright image (i.e. named gXXXX_YYYY.tif).
		if (substring(Files[i], lengthOf(Files[i])-8, lengthOf(Files[i])) != "dark.tif") {
			// relates current i to index of results
			index = i / 2;
			// get the theoretical gain = XXXX
			GainsTheo[index] = parseInt(substring(Files[i], 1, 5));
			// get the exposure time = YYYY
			Expo[index] = parseInt(substring(Files[i], lengthOf(Files[i])-8, lengthOf(Files[i])-4));

			// open the bright image
			open(InputDir + Files[i]);
			// restrict the measurement to the center of the image
			run("Specify...", "width=256 height=256 x=128 y=128");
			// get the mean intensity from the center area as "Bright Mean"
			getStatistics(area, BrightMean[index], min, max, sd, histo);
			close();	

//			print("i=" + i + " / GainsTheo[" + index + "]=" + GainsTheo[index] + " / Expo[" + index + "]=" + Expo[index] + " / BrightMean[" + index + "]=" + BrightMean[index]);
		}	
		// else the file is a dark image (i.e. named gXXXX_dark.tif)
		else {
			// relates current i to index of results
			index = (i-1) / 2;
			// open the dark image
			open(InputDir + Files[i]);
			// restrict the measurement to the center of the image
			run("Specify...", "width=256 height=256 x=128 y=128");
			// get the mean intensity from the center area as "Dark Mean"
			getStatistics(area, DarkMean[index], min, max, sd, histo);
			close();

			// Calculate the experimental gain value based on the Photometrics formula
			// Beware ! Formula modified to subtract the count of the dark image without gain (DarkMean[0]) to the count of the bright image without gain (BrightMean[0])
			// so the experimental gain = theoretical gain = 1 for the "no gain" situation
			GainsExp[index] = ((BrightMean[index] - DarkMean[index]) / (BrightMean[0] - DarkMean[0])) * (Expo[0]/Expo[index]);

//			print("i=" + i + " / DarkMean[" + index + "]=" + DarkMean[index] + " / GainsExp[" + index + "]=" + GainsExp[index]);
		}
		
	}


	// Generate the results table
	title1= "EM Gains Calibration Results";
	title2 = "["+title1+"]";
	f = title2;
	if (isOpen(title1))
		print(f, "\\Clear");
	else
		run("New... ", "name="+title2+" type=Table");
	Headings="\\Headings:n\tTheoric Gain\tExpo Time\tBright Mean\tDark Mean\tExp Gain";
	print(f,Headings);
	for (i=0; i<nCond; i++) {
		ResultsLine=d2s(i+1,0)+"\t"+GainsTheo[i]+"\t"+Expo[i]+"\t"+BrightMean[i]+"\t"+DarkMean[i]+"\t"+GainsExp[i];
		print(f,ResultsLine);
	}

setBatchMode("exit and display");	
}