//MUSICAL time-lapse from multi-colour stack 

img = getTitle();
img_without_extension = File.nameWithoutExtension;
dir = getDirectory("image");

getDimensions(width, height, n_colours, n_slices, n_frames);
if(n_frames < n_slices){
	tmp = n_frames;
	n_frames = n_slices;
	n_slices = tmp;
}

// make dialog box
Dialog.create("MUSICAL");
Dialog.addMessage("Channel specifications:");
default_colours = newArray(520, 520, 520, 520, 520, 520, 520);
default_thresholds = newArray(-0.5, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7);
for (c = 1; c <= n_colours; c++) {
	Dialog.addNumber("C" + c + " emission [nm]:", default_colours[c-1]);
	Dialog.addNumber("C" + c + " threshold:", default_thresholds[c-1]);
}

Dialog.addMessage("Time-lapse specifications:");
Dialog.addNumber("Batch size:", 150);
Dialog.addNumber("Slide by:", 50);

Dialog.addMessage("Microscope specifications:");
Dialog.addNumber("Pixel size [nm]:", 65);
Dialog.addNumber("Subpixels per pixel:", 4);
Dialog.addNumber("Magnification:", 1);
Dialog.addNumber("Alpha:", 4.0);
Dialog.addNumber("Numerical aperture:", 1.42);

Dialog.addCheckbox("Multithreading", true);
Dialog.addNumber("Threads", 4);
Dialog.show();

// read in user in-put
wavelengths = newArray(n_colours);
thresholds = newArray(n_colours);
for (c = 1; c <= n_colours; c++) {
	wavelengths[c-1] = Dialog.getNumber();
	thresholds[c-1] = Dialog.getNumber();
}
batch_size = Dialog.getNumber();
slide_by = Dialog.getNumber();
pixel_size = Dialog.getNumber();
sub_pixels = Dialog.getNumber();
magnification = Dialog.getNumber();
alpha = Dialog.getNumber();
numerical_aperture = Dialog.getNumber();
multi_threading = Dialog.getCheckbox();
threads = Dialog.getNumber();

n_MUSICAL_frames = floor((n_frames-batch_size)/slide_by) + 1;

if (n_colours>1) {
	run("Split Channels");
}
else {
	run("Duplicate...", "duplicate title=" + "C1-" + img);
}

for (c = 1; c <= n_colours; c++) {
	current_wavelength = wavelengths[c-1];
	current_threshold = thresholds[c-1];
	current_img = "C" + c + "-" + img;
	
	MUSICAL_frame_ID = 0; // initializing
	for (i = 1; i <= n_MUSICAL_frames; i++) {
		
		// define batch (substack numbers)
		batch_name = "batch_C" + c + "_" + i + ".tif";
		batch_path = dir + batch_name;
		current_batch_start = (i-1)*slide_by + 1;
		current_batch_end = current_batch_start + batch_size - 1;
		print("Current batch C" + c + " = [" + current_batch_start + ", " + current_batch_end + "]");
	
		// make current batch
		selectImage(current_img);
		run("Duplicate...", "duplicate title=" + batch_name + " range=" + current_batch_start + "-" + current_batch_end);
		saveAs("tif", batch_path); // temporary image stack (MusiJ v0.9 needs the raw file to be saved in a defined location)
		selectImage(current_img);
		current_frame = File.nameWithoutExtension + "_C" + c + "_" + i + "_[" + current_batch_start + "-" + current_batch_end + "].tif";

		// make MUSICAL image
		selectImage(batch_name);
		rename(current_frame);
		run("MUSICAL Image", "em=" + current_wavelength + 
							" na=" + numerical_aperture + 
							" mag=" + magnification +
							" pixsize=" + pixel_size +
							" subpixels=" + sub_pixels + 
							" threshold=" + current_threshold +
							" alpha=" + alpha +
							" multithreading=" + multi_threading +
							" threads=" + threads +
							" save=" + true);
							
		// add MUSICAL image to output stack
		run("Copy");
		w = getWidth; h = getHeight; 
	  	close();
	  	if (MUSICAL_frame_ID == 0) {
	  		selectImage(current_img);
	  		MUSICAL_stack_name = File.nameWithoutExtension  + "_C" + c + "_MuImS";
		    newImage(MUSICAL_stack_name, "Grey", w, h, 1); 
		    MUSICAL_frame_ID = getImageID; 
		} 
		else { 
		   selectImage(MUSICAL_frame_ID); 
		   run("Add Slice"); 
		}
		run("Paste");
	
	  	// tidy up
	  	selectImage(current_frame);
		close();
		dummy_variable = File.delete(batch_path);
		
	}
}

// save combined colour stack
merge_string = "";
for (c = 1; c <= n_colours; c++) {
	merge_string = merge_string + "c" + c + "=" + img_without_extension + "_C" + c + "_MuImS ";
}
run("Merge Channels...", merge_string + " create ignore");
run("Re-order Hyperstack ...", "channels=[Channels (c)] slices=[Frames (t)] frames=[Slices (z)]");
final_MUSICAL_stack_path = dir + img_without_extension + "_MuImS.tif";
saveAs("tif", final_MUSICAL_stack_path);

// tidy up
for (c = 1; c <= n_colours; c++) {
	selectImage("C" + c + "-" + img);
	close();
}
if (n_colours==1) {
	selectImage(img);
	close();
}

// save parameters to file
settings_name = img_without_extension + "_settings.txt";
settings_file = File.open(dir + settings_name);
for (c = 1; c <= n_colours; c++) {
	print(settings_file, "Emission wavelength " + c +" [nm]: " + wavelengths[c-1]);
	print(settings_file, "Threshold " + c +": " + thresholds[c-1]);
}
print(settings_file, "Batch size: " + batch_size);
print(settings_file, "Slide by: " + slide_by);
print(settings_file, "Pixel size [nm]: " + pixel_size);
print(settings_file, "Subpixels [nm]: " + sub_pixels);
print(settings_file, "Magnification: " + magnification);
print(settings_file, "Alpha: " + alpha);
print(settings_file, "Numerical aperture: " + numerical_aperture);
if (multi_threading) {
	print(settings_file, "Number of threads in multi-threading: " + threads);
}
File.close(settings_file);