// Get_Chroma_Coefs macro by Christophe Leterrier
// 27-01-2014
// last version 15-04-2021
//
// Processes a folder full of multichannel tif, register each channel relative to the first one using the Image Stabilzer plugin;
// (requires Image Stabilizer plugin, with the 2010 bug patch for Log Applier batch run)
// http://www.kangli.org/code/Image_Stabilizer.html
// then generates TransformJ matrices from the shifts stored in the Image Stabilizer log
//
// Correction of stacks is then made using the Correct_STORM or Correct_Apo macros, applying the TrasnformJ matrices to each channel

macro "Get_Chroma_Coefs" {

//	IN_PATH = "/Volumes/Bushwick/Unison/Data Raw/140203 Etalonnage/Tetraspeck (chroma correction)/040X_1.00X_4c_10f/040X_1.00X_4c_10f.tif - T=0.tif";
//	IN_NAME = File.getName(IN_PATH);

	CORR_MODE = "Translation";
//	CORR_MODE = "Affine";
	// channel number as reference (1 = channel 0, etc.)
	REF_CHANNEL = 1;

	INPUT_DIR=getDirectory("Select the input stacks directory");
	ALL_NAMES=getFileList(INPUT_DIR);
	Array.sort(ALL_NAMES);

	print("\n\n\n*** Get_Chroma_Coefs Log ***");
	print("");
	print("INPUT_DIR :"+INPUT_DIR);

	// Stack counter for average
	f = 0;

	// Loop on tif images
	for (n=0; n<ALL_NAMES.length; n++) {
		NAME_PARTS = split(ALL_NAMES[n], ".");
		NAME_EXT = NAME_PARTS[NAME_PARTS.length - 1];

		if (NAME_EXT == "tif") {

			// increment stack counter
			f++;

			// Open stack
			open(INPUT_DIR + ALL_NAMES[n]);

			// Create flags & sum arrays at first stack
			if (f == 1) {
				getPixelSize(unit, pw, ph, pd);
				N_CHANNELS = nSlices;

				F1 = newArray(N_CHANNELS);
				F2 = newArray(N_CHANNELS);
				C3 = newArray(N_CHANNELS);
				C4 = newArray(N_CHANNELS);

				if (CORR_MODE == "Affine") {
					C5 = newArray(N_CHANNELS);
					C6 = newArray(N_CHANNELS);
					C7 = newArray(N_CHANNELS);
					C8 = newArray(N_CHANNELS);

				}
			}

			print("\nStack name: " + ALL_NAMES[n]);
			print("Number of channels: " + N_CHANNELS);
			print("Pixelsize: " + pw);

			//Put stack at reference channel before launching image stabilizer
			setSlice(REF_CHANNEL);

			// Run image stabilizer plugin
			run("Image Stabilizer", "transformation=" + CORR_MODE + " maximum_pyramid_levels=1 template_update_coefficient=0.90 maximum_iterations=200 error_tolerance=0.0000001 log_transformation_coefficients");
			close();

			// Get image stabilizer log
			IN_LOG = getLogName(ALL_NAMES[n]);
			selectWindow(IN_LOG + ".log");
			LOG_STRING = getInfo("window.contents");
			run("Close");

			// Split lines of log
			LOG_SPLIT = split(LOG_STRING, "\n");

			// Loop on channel lines from log
			N_SLICES = LOG_SPLIT.length - 2;
			for (j = 0; j < N_SLICES; j++) {

				print("Coefs for channel #" + (j+1) + ": " + LOG_SPLIT[j+2]);
				// Split channel lines to get coefs
				LINE_SPLIT = split(LOG_SPLIT[j+2], ",");

				// Initialize flag and sum arrays values at first stack
				if (f==1) {
					F1[j] = parseInt(LINE_SPLIT[0]);
					F2[j] = parseInt(LINE_SPLIT[1]);
					C3[j] = parseFloat(LINE_SPLIT[2]);
					C4[j] = parseFloat(LINE_SPLIT[3]);
					if (CORR_MODE == "Affine") {
						C5[j] = parseFloat(LINE_SPLIT[4]);
						C6[j] = parseFloat(LINE_SPLIT[5]);
						C7[j] = parseFloat(LINE_SPLIT[6]);
						C8[j] = parseFloat(LINE_SPLIT[7]);
					}
				}

				// Add value to sum arrays for subsequent stacks
				else {
					C3[j] = C3[j] + parseFloat(LINE_SPLIT[2]);
					C4[j] = C4[j] + parseFloat(LINE_SPLIT[3]);
					if (CORR_MODE == "Affine") {
						C5[j] = C5[j] + parseFloat(LINE_SPLIT[4]);
						C6[j] = C6[j] + parseFloat(LINE_SPLIT[5]);
						C7[j] = C7[j] + parseFloat(LINE_SPLIT[6]);
						C8[j] = C8[j] + parseFloat(LINE_SPLIT[7]);
					}
				}
			} // end of for loop on channels

		} // end of if loop on tif images
	} // end of for loop on all files


	// Create average arrays for coefs
	A3 = newArray(N_CHANNELS);
	A4 = newArray(N_CHANNELS);

	if (CORR_MODE == "Affine") {
		A5 = newArray(N_CHANNELS);
		A6 = newArray(N_CHANNELS);
		A7 = newArray(N_CHANNELS);
		A8 = newArray(N_CHANNELS);
	}

	// Create string array for coefs string
	A_ARRAY = newArray(N_CHANNELS);
	M_ARRAY = newArray(N_CHANNELS);

	// Calculate average coefs for each channel

	print ("\n");
	for (j = 0; j < N_CHANNELS; j++) {

		A3[j] = C3[j] / f;
		A4[j] = C4[j] / f;

		L_STRING =	"Average coefficients for channel #" + (j+1) + ": ";
		A_STRING = "" + F1[j] + "," + F2[j] + "," + A3[j] + "," + A4[j];


//		M_STRING = "1\t0\t0\t" + (A3[j]*pw) + "\n0\t1\t0\t"+ (A4[j]*pw) +"\n0\t0\t1\t0\n0\t0\t0\t1\n";
		M_STRING = "1\t0\t0\t" + (-A3[j]*pw) + "\n0\t1\t0\t"+ (-A4[j]*pw) +"\n0\t0\t1\t0\n0\t0\t0\t1\n";

		if (CORR_MODE == "Affine") {
			A5[j] = C5[j] / f;
			A6[j] = C6[j] / f;
			A7[j] = C7[j] / f;
			A8[j] = C8[j] / f;

			A_STRING = A_STRING + "," + A5[j] + "," + A6[j] + "," + A7[j] + "," + A8[j];
//			M_STRING = "" + (1-A3[j]) + "\t" + (-A4[j]) + "\t0\t" + (-A5[j]*pw) + "\n" + (-A6[j]) + "\t" + (1-A7[j]) + "\t0\t" + (-A8[j]*pw) +"\n0\t0\t1\t0\n0\t0\t0\t1\n";
			M_STRING = "" + (1-A3[j]) + "\t" + (-A4[j]) + "\t0\t" + (-A5[j]) + "\n" + (-A6[j]) + "\t" + (1-A7[j]) + "\t0\t" + (-A8[j]) +"\n0\t0\t1\t0\n0\t0\t0\t1\n";
		}

	print(L_STRING + A_STRING);
	print("Matrix:\n" + M_STRING);

	A_ARRAY[j] = A_STRING;
	M_ARRAY[j] = M_STRING;

	File.saveString(M_STRING, INPUT_DIR + "matrix_CH" + j + ".txt");

	}


	// Build the log file for averaged coefficients
	OUT_STRING = "Image Stabilizer Log File for \"averaged coefs\"\n";
	if (CORR_MODE == "Translation") OUT_STRING += "0\n";
	else OUT_STRING += "1\n";
	for (j = 0; j < N_CHANNELS; j++) {
		OUT_STRING += A_ARRAY[j] + "\n";
	}

	print("\nAverage Image Stabilizer Log:");
	print(OUT_STRING);


	File.saveString(OUT_STRING, INPUT_DIR + "average_coefs.log");
	print("\nAveraged log saved to: "+ INPUT_DIR + "average_coefs.log");

	print("\n*** Get_Chroma_Coefs end ***");
	showStatus("Get_Chroma_Coefs finished");
}




//*************** Functions ***************

function getLogName(Name) {
	tifindex = indexOf(Name, ".tif");
	zviindex = indexOf(Name, ".zvi");
	if (zviindex > 0) index = minOf(tifindex, zviindex);
	else index = tifindex;
	short = substring(Name, 0, index);
//	print(short);
	return short;
}