/**
 * Make montage stacks from MetaMorph 7 High Throughput Screening.
 * 
 * MetaMorph 7's Screen Acquisition (screenacq) does not natively
 * support time and Z acquisitions from it's UI or journal functions.
 * Therefore it's easier to assemble the experiment with FIJI
 * including stitch TIF files together with overlap.  It is
 * assumed the grid of all sites in a well are imaged, i.e. supports
 * Metamorph's "Well Selection" but not "Site Selection".
 * 
 * Developer Notes:
 * - loci_tools does not index Metamorph's HTD TIF files correctly
 *   (it reads all image leaves as index 1).  The stitching plugin
 *   uses loci_tools. Thus the files cannot be used directly and the
 *   HTD file needs to be need to be renamed so that loci does not
 *   recognize it.
 * - Written against Version 1.0 of HTSInfoFile
 * 
 * Written by Pariksheet Nanda <omsai@member.fsf.org> Oct, 2012
 * License: Public Domain
 * Version: 0.7
 */


// Global variables (for use by functions)
var datadir;  // Path to HTS files containing HTD and TIF files.
var htd_string;  // HTD file fully read in as a string.


macro 'Unused Tool-2 - ' {}  // empty slot


macro 'Stitch Raw Data Action Tool - C000 R00ffL505fLa0afL05f5L0afa' {
    // Boilerplate
    setBatchMode(true);
    saveSettings();

    datadir = getDirectory('Choose your raw "Screen Acquisition" data folder');
    if (unhide_htd_files(datadir))
	exit('Exiting... Restored .HTD file naming from previous failed run.');
    if (is_filenaming_invalid(datadir)) {
	ok_to_rename =
	    getBoolean('Extra "." characters found in filenames which will ' +
		       'confuse the Stitching / Bioformats plugins.\n' +
		       'Ok to replace them with "_" characters?');
	if (ok_to_rename == 0)
	    exit('Exiting... filenames contain invalid characters, ' +
		 'which will crash the macro');
	else {
	    IJ.log('Renaming files...');
	    IJ.log('... ' + fix_filenaming(datadir) + ' files renamed.');
	}
    }
    Dialog.create('Hyperstack');
    Dialog.addNumber('Z slices:', 1);
    Dialog.show();
    slices = Dialog.getNumber();
    outputdir = datadir + 'montages/';
    finaldir = datadir + 'hyperstacks/';
    if (! File.exists(outputdir))
	File.makeDirectory(outputdir);
    if (! File.exists(finaldir))
	File.makeDirectory(finaldir);
    htd_files = get_files(datadir, 'HTD');
    // make sure number of z slices is sane, otherwise the hyperstack
    // convertor will fail later.
    extra_frames = htd_files.length % slices;
    if (extra_frames != 0) {
	Dialog.create('Warning');
	Dialog.addMessage('You have entered ' + slices + ' Z slices, but ' +
			  'there are ' + htd_files.length + ' timepoints. ' +
			  'This does not divide evenly.');
	Dialog.addMessage('Click OK to ignore the ' + extra_frames + 
			  ' extra time points ' +
			  '(if, say, you ended the acquisition prematurely) ' +
			  'or cancel, to restart the macro and enter the ' +
			  'correct number of z slices.');
	Dialog.show();
    }
    // Hide the HTD file from Bioformats since it does not read in the
    // well metadata and repeatedly stitches the first frame only
    if (hide_htd_files(datadir) == 0)
	exit('Exiting... No HTD files found.');
    
    // Iterate over all HTD acquisition primitives to create stack
    // of experiment wells
    htd_string = File.openAsString(datadir + htd_files[0] + '.bak');
    wells = get_well_coordinates();
    x_sites = parseInt(get_htd_value('XSites'));
    y_sites = parseInt(get_htd_value('YSites'));
    total_sites = x_sites * y_sites;
    channels = get_wavelengths();

    for (well_index = 0; well_index < wells.length; well_index++) {
    	well = wells[well_index];
	file_prefix = '';

	// save individual files
	for (wavelength = 0; wavelength < channels.length; wavelength++) {
	    for (acq = 0; acq < htd_files.length; acq++) {
		show_progress(1 +
			      acq +
			      wavelength * htd_files.length +
			      well_index * channels.length * htd_files.length,
			      wells.length * channels.length * htd_files.length);
		base = prefix(htd_files[acq]);
		file_prefix = base + '_' + well;
		if (channels.length == 1)
		    channel_prefix = '';
		else
		    channel_prefix = '_w' + wavelength+1;
		tif_file_name = file_prefix + '_s{i}' + channel_prefix + '.TIF';
		// Plugins > Stitching > Grid/Collection stitching
		run('Grid/Collection stitching',
		    'type=[Grid: row-by-row] ' +
		    'order=[Right & Down                ] ' +
		    'grid_size_x=&x_sites grid_size_y=&y_sites ' +
		    'tile_overlap=20 ' +
		    'first_file_index_i=1 ' +
		    'directory=&datadir ' +
		    'file_names=[&tif_file_name] ' +
		    'fusion_method=[Linear Blending] ' +
		    'regression_threshold=0.30 ' +
		    'max/avg_displacement_threshold=2.50 ' +
		    'absolute_displacement_threshold=3.50 ' +
		    'subpixel_accuracy ' +
		    'computation_parameters=[Save computation time ' +
		    '(but use more RAM)] ' +
		    'image_output=[Fuse and display]');
		save_file_name = file_prefix + '_w' + wavelength+1;
		saveAs('tif', outputdir + save_file_name);
		IJ.log('Saved montage ' + save_file_name);
		close();
	    }
	}
	
	if (htd_files.length > 1 || channels.length > 1) {
	    // assemble hyperstack
	    tif_file_name = outputdir + file_prefix + '_w' + 1 + '.tif';
	    // File > Import > Image Sequence
	    run('Image Sequence...',
		'open=[&tif_file_name] ' + 
		'starting=1 ' + 
		'increment=1 ' +
		'scale=100 ' + 
		'file=[&well] ' + 
		'or=[] ' +
		'sort');
	    // strip extra frames if experiment ended prematurely
	    frames = floor(htd_files.length / slices);
	    total_channels = channels.length;
	    extra_frames = htd_files.length % slices;
	    if (extra_frames != 0) {
		IJ.log('Deleting last ' + extra_frames * total_channels +
		       ' frames from stack of ' + nSlices + '...');
		Stack.setSlice(htd_files.length * total_channels);
		for (i = 0; i < extra_frames * total_channels; i++)
		    run('Delete Slice');
		IJ.log('...new stack size is ' + nSlices);
	    }
	
	    if (nSlices > 1) {
		// convert to hyperstack

		// Running 'Stack to Hyperstack...' with batch mode on in
		// ImageJ 1.47c goves the following exception when turning off
		// batch mode, so as a workaround we turn off batch mode.
		//
		// java.lang.NullPointerException
		// at ij.CompositeImage.setupLuts(CompositeImage.java:139)
		// at ij.CompositeImage.updateImage(CompositeImage.java:201)
		// at ij.CompositeImage.getImage(CompositeImage.java:97)
		// at ij.ImagePlus.show(ImagePlus.java:364)
		// at ij.ImagePlus.show(ImagePlus.java:347)
		// at ij.macro.Functions.displayBatchModeImage(Functions.java:2605)
		// at ij.macro.Functions.setBatchMode(Functions.java:2588)
		// at ij.macro.Functions.doFunction(Functions.java:131)
		// at ij.macro.Interpreter.doStatement(Interpreter.java:216)
		// at ij.macro.Interpreter.doBlock(Interpreter.java:539)
		// at ij.macro.Interpreter.runMacro(Interpreter.java:131)
		// at ij.macro.MacroRunner.run(MacroRunner.java:143)
		// at java.lang.Thread.run(Thread.java:662)
		setBatchMode(false);
		run('Stack to HyperStack...',
		    'order=xyczt(default) ' +
		    'channels=&total_channels ' +
		    'slices=&slices ' + 
		    'frames=&frames ' + 
		    'display=Color');
		// using slices=1 above since one can always reslice the stack
		// to set Z planes, i.e. converting to a stack and then back
		// to hyperstack.
		// Image > Hyperstacks > Stack to Hyperstack
		setBatchMode(true);
		base = prefix(htd_files[0]);
		file_prefix = base + '_';
		save_file_name = file_prefix + well;
		saveAs('tif', finaldir + save_file_name);
		close();
		IJ.log('Saved hyperstack ' + well);
	    }
	}
    }
    IJ.log('Restoring .HTD files...');
    unhide_htd_files(datadir);
    IJ.log('...restored');
    
    // Boilerplate
    setBatchMode(false);
    restoreSettings();
    IJ.log('Finished HTS Grid Stitching');
}


//macro 'Review Stitched Wells Action Tool - C111 O0966 O0066 O9966 O9066' {}


// Easy to do this from ImageJ, but saves a few mouse clicks
macro 'Hyperstack Z Projection Action Tool - C000 P62f2a71762 Pc5f5aa1a47 Pc8f8ad1d4a' {
    setBatchMode(true);
    Stack.getDimensions(width, height, channels, slices, frames);
    if (slices == 1)
	exit('Image only has 1 Z slice');
    // Reduce depth of hyperstack to a single MIP for each frame
    run('Z Project...',
	'start=1 ' +
	'stop=&slices ' +
	'projection=[Max Intensity] ' +
	'all');
    setBatchMode(false);
}


/**
 * Extract basename without extension.
 * 
 * @param  directory path where files are located
 * @return           string of file list
 */
function prefix(directory) {
    file_name = File.getName(directory);
    extension_index = lastIndexOf(file_name, '.');
    base_name = substring(file_name, 0, extension_index);
    return base_name;
}

/**
 * Find files with a specific extension.
 * 
 * @param  directory path where files are located
 * @param  extension desired files
 * @return           array file list
 */
function get_files(directory, extension) {
    all_file_list = getFileList(directory);
    match_file_list = newArray();
    for (i = 0; i < all_file_list.length; i++) {
	file_split = split(File.getName(all_file_list[i]), '.');
	file_extension = Array.slice(file_split, file_split.length-1,
				     file_split.length);
	if (matches(file_extension[0], extension)) {
	    match_file = Array.slice(all_file_list, i, i+1);
	    match_file_list = Array.concat(match_file_list, match_file);
	}
    }
    return match_file_list;
}

/**
 * Check for filename characters that upset the Stitching plugin.
 * 
 * @param  directory path where files are located
 * @return           1 if any filename invalid, otherwise 0
 */
function is_filenaming_invalid(directory) {
    htd_files = get_files(directory, 'HTD');
    for (i = 0; i < htd_files.length; i++)
	if (indexOf(prefix(htd_files[i]), '.') != -1)
	    return 1;
    return 0;
}	    

/**
 * Replace invalid TIF filename characters with underscores in directory.
 * 
 * @param  directory path where files are located
 * @return           number of files renamed
 */
function fix_filenaming(directory) {
    types_to_rename = newArray('HTD', 'TIF');
    files_renamed = 0;
    for (t = 0; t < types_to_rename.length; t++) {
	extension = types_to_rename[t];
	files = get_files(directory, extension);
	for (i = 0; i < files.length; i++) {
	    file = files[i];
	    base = prefix(file);
	    if (indexOf(base, '.') != -1) {
		new_name = replace(base, '.', '_') + '.' + extension;
		File.rename(directory + File.separator + file,
			    directory + File.separator + new_name);
		files_renamed++;
	    }
	}
    }
    return files_renamed;
}

/**
 * Read HTD metadata using key names.
 * 
 * Lines in the HTD file format consist of 'key_name', value
 * 
 * @param  key_name key contained in first column of HTD file
 * @return          corresponding value read from HTD file
 */
function get_htd_value(key_name) {
    key_index = indexOf(htd_string, key_name);
    if (key_index == -1)
	return 'Error: key "' + key_name + '" not found in HTD file';
    line_end_index = indexOf(htd_string, '\n', key_index);
    KEY_VALUE_SEPARATOR = '", ';
    value_index = key_index + lengthOf(key_name) + lengthOf(KEY_VALUE_SEPARATOR);
    value = substring(htd_string, value_index, line_end_index);
    return value;
}

/**
 * Well coordinates of images from reading HTD metadata.
 * 
 * @return array of coordinates like A01,A02,B03,...
 */
function get_well_coordinates() {
    x_wells = parseInt(get_htd_value('XWells'));
    y_wells = parseInt(get_htd_value('YWells'));
    wells = newArray();
    for (i = 1; i <= y_wells; i++) {
	row = get_htd_value('WellsSelection' + i);
	row = replace(row, ' ', '');  // strip spaces
	row = split(row, ',');
	for (j = 1; j <= x_wells; j++) {
	    // Populate array with strings like A01, A02, etc
	    if (row[j-1] == 'TRUE') {
		coordinates = fromCharCode(64 + i) + IJ.pad(j,2);
		wells = Array.concat(wells, coordinates);
	    }
	}
    }
    return wells;
}

/**
 * Wavelengths of images from reading HTD metadata.
 * 
 * @return array of wavelengths like Confocal 488,Cam Trans...
 */
function get_wavelengths() {
    wavelengths = newArray();
    total_wavelengths = parseInt(get_htd_value('NWavelengths'));
    for (i = 1; i <= total_wavelengths; i++) {
	wavelength = get_htd_value('WaveName' + i);
	wavelengths = Array.concat(wavelengths, wavelength);
    }
    return wavelengths;
}

/**
 * Strips .bak from all directory .HTD files.
 *
 * @return 1 if any file renamed, otherwise 0
 */
function unhide_htd_files(directory) {
    // this can be fail if there is a non-HTD bak file because it will
    // only rename HTD bak files
    bak_files = get_files(directory, 'bak');
    for (bak = 0; bak < bak_files.length; bak++) {
 	ret = File.rename(directory + bak_files[bak],
			  directory + replace(bak_files[bak],
					      '.HTD.bak', '.HTD'));
    }
    if (bak > 0)
	return 1;
    return 0;
}

/**
 * Adds .bak to all directory .HTD files.
 *
 * @return 1 if any file renamed, otherwise 0
 */
function hide_htd_files(directory) {
    htd_files = get_files(directory, 'HTD');
    for (acq = 0; acq < htd_files.length; acq++) {
	ret = File.rename(directory + htd_files[acq],
			  directory + htd_files[acq] + '.bak');
    }
    if (acq > 0)
	return 1;
    return 0;
}

/**
 * Show progress string.
 * 
 * @param  current_value index data point
 * @param  total	 total data points
 * @return string indicating progress, like ****---------------
 */
function _get_bar(current_value, total) {
    n = 20;
    empty_bar = '--------------------';
    complete_bar = '********************';
    index = round(n * (current_value/total));
    if (index < 1)
	index = 1;
    if (index > n-1)
	index = n-1;
    return substring(complete_bar, 0, index) + substring(empty_bar, index+1, n);
}

/**
 * Show progress in window.
 *
 * Adapted from http://imagej.nih.gov/ij/macros/ProgressBar.txt
 *
 * @param  current_value index data point
 * @param  total	 total data points
 */
function show_progress(current_value, total) {
    title = '[Progress]';
    // create progress window for first run
    windows = getList('window.titles');
    found_progress_window = 0;
    clean_title = replace(title, '\\[', '');
    clean_title = replace(clean_title, '\\]', '');
    for (i = 0; i < windows.length; i++) {
	if (matches(windows[i], clean_title))
	    found_progress_window = 1;
    }
    if (found_progress_window == 0)
	run('Text Window...',
	    'name='+ title +' width=27 height=3 monospaced');
    // show progress
    if (current_value >= total) {
    	print(title, '\\Close');
    }
    else {
	print(title,
	      '\\Update:' + current_value + '/' + total +
	      ' (' + (current_value / total) * 100 + '%)\n' +
	      _get_bar(current_value, total));
    }
}