Viveca Stitching Tool

/*
VIVECA STITCHING
Started in January 2019
Macro build and maintained by Vasco Fachada (vasco.fachada@gmail.com)

At the moment this macro is using the stitching plugin "Grid/Collection stitching" used in Preibisch et al. Bioinformatics (2009).
This macro basically goes throught the folder of the user's choice and opens .czi files.
The metadata from the .czi is read and the individual tiles of the file are saved as orderly tiff files in a temp folder. 
The read metada is used to feed the mentioned plugin the parameters for stitching. These include number of tiles in each axis.
The laser excitation data is also read and the image should display: 405->Blue; 488->Green; 555->Orange; 639->Red.
The tile overlap data is also read from the metdata and fed into the stitching plugin
Z-stack stiching is supported
Both modes of tile scanning in ZEN black are supported, uni and bidirectional
*/

close("*");
Dialog.create("Ok, grab your needle. This is just like knitting");
Dialog.addMessage("Welcome to our Viveca stitching tool. \n");
Dialog.addMessage("\nThis macro is mainly aimed for stitching tiled-images produced by our Zeiss LSM700. \nAt the moment, stitching is only working with .czi files which cannot contain spaces in their file name.");
Dialog.addCheckbox("  Display your images once they are stitched?", true); 
Dialog.addCheckbox("\n  My tiles were scanned in 'zigzag / snake / bidirection' mode", false);
Dialog.addMessage("Constantly under development. Please consider acknowledging the author in publications. \nAny issues should be reported in order to improve the tool. \n \nauthor: Vasco Fachada (vasco.fachada(at)gmail.com)");
Dialog.show();
visualize_checkBox = Dialog.getCheckbox();
tilingMode_checkBox = Dialog.getCheckbox();

//Finding the .czi files to stitch
dir = getDirectory("Select a directory containing your CZI files you want to stitch");
start = getTime();
files = getFileList(dir);
outputFold = dir+"Stitchings";
File.makeDirectory(outputFold);

setBatchMode(true);

//Deleting any previously produced folders in the case the macro broke half-way  
for(f=0; f<files.length; f++) {
	if(matches(files[f], ".*stitch_temp.*")){
		old_temp_folder = dir+File.separator+files[f];
		old_content = getFileList(old_temp_folder);
		for (li=0; li<old_content.length; li++){
			old_images = getFileList(old_temp_folder+File.separator+old_content[li]);
			for (lii=0; lii<old_images.length; lii++){
				File.delete(old_temp_folder+File.separator+old_content[li]+File.separator+old_images[lii]);
			} 
			File.delete(old_temp_folder+File.separator+old_content[li]);
		}
  		File.delete(old_temp_folder);	
		print("successfully deleted previously produced temp files");
	}
}


k=0;
n=0;
Tot_nmrOfTilesOpened = 0;
nmrOfFilesStitched = 0;
resultArray = newArray(0);
pixelTotal = 0;
sizeLoaded = 0;

	run("Bio-Formats Macro Extensions");
	for(f=0; f<files.length; f++) {
		if(endsWith(files[f], ".czi")) {
			nmrOfTilesOpened = 0; 
			k++;
			id = dir+files[f];
			//print(files[f]);
			//waitForUser;
			Ext.setId(id);
			Ext.getSeriesCount(seriesCount);
			n+=seriesCount;


			//READING METADATA
			run("Bio-Formats Importer", "open=["+id+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack series_"+(1));
			fullMeta=getMetadata("Info");
	
			//Determining if 2D or 3D
			zSteps = parseInt(substring(fullMeta, indexOf(fullMeta, "SizeZ =")+8, indexOf(fullMeta, "SizeZ =")+10)); 
			print("z-steps: "+zSteps); 
				
			//Determining if timeseries
			tPoints = parseInt(substring(fullMeta, indexOf(fullMeta, "SizeT =")+8, indexOf(fullMeta, "SizeT =")+9));
			print("timepoints: "+tPoints);
			
			//Determining number of positions. Ine the case of the 12-well plate container functions in LSM700 Zen black, each position corresponds to each well.
			if (indexOf(fullMeta, "Information|Image|SizeS #1 =") > -1){ 
				positionsNum = substring(fullMeta, lastIndexOf(fullMeta, "Information|Image|SizeS #1 =")+lengthOf("Information|Image|SizeS #1 =")+0, lastIndexOf(fullMeta, "Information|Image|SizeS #1 =")+lengthOf("Information|Image|SizeS #1 =")+3);
				positionsNum = parseInt(positionsNum);				
			}else
				positionsNum = 1;
			print("number of positions: "+positionsNum);
								
			//Getting the number of tiles in each axis for later stitching steps
			posStartInd = 3+lengthOf(toString(positionsNum));
			posLastInd = 4+lengthOf(toString(positionsNum));	
			gridMetaX = substring(fullMeta, lastIndexOf(fullMeta, "PositionGroup|TilesX #")+lengthOf("PositionGroup|TilesX #")+posStartInd, lastIndexOf(fullMeta, "PositionGroup|TilesX #")+lengthOf("PositionGroup|TilesX #")+posLastInd);
			gridMetaY = substring(fullMeta, lastIndexOf(fullMeta, "PositionGroup|TilesY #")+lengthOf("PositionGroup|TilesY #")+posStartInd, lastIndexOf(fullMeta, "PositionGroup|TilesY #")+lengthOf("PositionGroup|TilesY #")+posLastInd);
			gridX = parseInt(gridMetaX);
			print("X grid: "+gridX);
			gridY = parseInt(gridMetaY);
			print("Y grid: "+gridY);
			TilesPerPos = gridX * gridY;
			print("Tiles per position: "+TilesPerPos);
			TotalNumTiles = TilesPerPos * positionsNum * tPoints;
			print("total number of tiles: "+TotalNumTiles);
			if (TotalNumTiles < 10)
				mosaicSize = "{i}";
			else if (TotalNumTiles > 9 && TotalNumTiles < 100)
				mosaicSize = "{ii}";
			else if (TotalNumTiles > 99 && TotalNumTiles < 1000)
				mosaicSize = "{iii}";
			else
				mosaicSize = "{iiii}";

			print("mosaic size: "+mosaicSize);
			//Reading tile scanning direction
			if (tilingMode_checkBox==1)
				tilingDirection = "[Grid: snake by rows] order=[Right & Down                ]";
			else
				tilingDirection = "[Grid: row-by-row] order=[Right & Down                ]";
			print("Tiling direction: "+tilingDirection);	
			
			//Reading tile overlap
			tlOv = substring(fullMeta, indexOf(fullMeta, "TileAcquisitionOverlap #")+28, indexOf(fullMeta, "TileAcquisitionOverlap #")+30);
			if (tlOv == '0.')
				tlOverlap = 100*parseFloat(substring(fullMeta, indexOf(fullMeta, "TileAcquisitionOverlap #")+28, indexOf(fullMeta, "TileAcquisitionOverlap #")+39));
			else
				tlOverlap = 0;
			print("tile overlap: "+tlOverlap);

			//Looking and determining number of channels in data
			channels_nmr = 	parseInt(substring(fullMeta, indexOf(fullMeta, "SizeC =")+8, indexOf(fullMeta, "SizeC =")+9)); 
			channels = newArray(channels_nmr);
			colors = newArray(channels_nmr);
			print(channels_nmr+" channels:");
				
			for (ii=0; ii<channels_nmr; ii++){
				laserUsed = substring(fullMeta, indexOf(fullMeta, "Attenuator|Laser #"+ii+1+" =")+27, indexOf(fullMeta, "Attenuator|Laser #"+ii+1+" =")+31);
				channels[ii] = parseInt(laserUsed);
				//print(channels[ii]);	
				if (channels[ii] == 405) colors[ii] = "Blue";
				else if (channels[ii] == 488) colors[ii] = "Green";
				else if (channels[ii] == 555) colors[ii] = "Yellow";
				else if (channels[ii] == 639){
					colors[ii] = "Red";
					print(colors[ii]+" with "+laserUsed+" laser");
				//else print(colors[ii]);
				}else exit("While reading the metadata we have been unable to find which laser was used to scan one or more of the channels. \nPlease contact vasco.fachada@gmail.com for help.");	
			}
			
		//waitForUser;
		
			getVoxelSize(width, height, depth, unit); //detecting the voxel size of the current raw file, so the final image also has the same scale.
			close();
			
			//Creating temporary tiff files for analysis and organization of stiched file
			for (l=1; l<positionsNum+1; l++){
							
				savingDir_temp = outputFold+File.separator+"stitch_temp";
				File.makeDirectory(savingDir_temp);
				positionDir = savingDir_temp+File.separator+"position"+l;
				File.makeDirectory(positionDir);

				if (seriesCount > 1){
					for (i=1; i<TilesPerPos+1; i++) {
						run("Bio-Formats Importer", "open=["+id+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack series_"+(i+nmrOfTilesOpened));

						fullName	= getTitle();
						//dirName 	= substring(fullName, 0,lastIndexOf(fullName, ".czi"));
						fileName 	= substring(fullName, lastIndexOf(fullName, " - ")+3, lengthOf(fullName));
						//5print(fileName);
						fileSeries  = substring(fileName, indexOf(fileName, "#"), lengthOf(fileName)); //collects the number of series starting from the hashtag. #1 or #0001
				
						saveAs("tiff", positionDir+File.separator+fileName+".tif");									
						close();
					}
					
				
					//STICHING AND VISUALIZATION
		//			setBatchMode(false);
					//"Grid/Collection stitching" used in Preibisch et al. Bioinformatics (2009).
					run("Grid/Collection stitching", "type="+tilingDirection+" grid_size_x="+gridX+" grid_size_y="+gridY+" tile_overlap="+tlOverlap+" first_file_index_i="+1+nmrOfTilesOpened+" directory="+positionDir+" file_names=["+substring(fileName, 0, lengthOf(fileName)-lengthOf(fileSeries))+"#"+mosaicSize+".tif] output_textfile_name=TileConfiguration.txt fusion_method=[Linear Blending] regression_threshold=0.30 max/avg_displacement_threshold=2.50 absolute_displacement_threshold=3.50 ignore_z_stage computation_parameters=[Save memory (but be slower)] image_output=[Fuse and display]");
					//run("Grid/Collection stitching", "type= grid_size_x="+gridX+" grid_size_y="+gridY+" tile_overlap="+tlOverlap+" first_file_index_i="+1+nmrOfTilesOpened+" directory="+positionDir+" file_names=["+substring(fileName, 0, lengthOf(fileName)-lengthOf(fileSeries))+"#"+mosaicSize+".tif] output_textfile_name=TileConfiguration.txt fusion_method=[Linear Blending] regression_threshold=0.30 max/avg_displacement_threshold=2.50 absolute_displacement_threshold=3.50 ignore_z_stage computation_parameters=[Save memory (but be slower)] image_output=[Fuse and display]");
					nmrOfTilesOpened = nmrOfTilesOpened + TilesPerPos;

					//setBatchMode(false);
					//for (xx=0;xx<colors.length;xx++){
						//print(colors[xx]);
					//}
					if (isOpen('Fused') == 1)
						print('Stiching successful!');
					else{
						waitForUser('There was a problem. \nCheck for instance if you have spaces in your file path, there should be none. \nRestart the stitching.');
					}
					
					if (colors.length<2) 
						run(colors[0]);
					else{
						for (ii=0; ii<colors.length; ii++){ //Assigning colors to final image corresponding to excitation laser wavelengths
							Stack.setDisplayMode("color");
							Stack.setChannel(ii+1);
							run(colors[ii]);
						}
						run("Make Composite", "stack");	
					}
			
					//Setting the native scale
					run("Properties...", "channels="+channels_nmr+" slices="+zSteps+" frames="+tPoints+" unit="+unit+" pixel_width="+width+" pixel_height="+height+" voxel_depth="+depth+"");
					//Stack.setDisplayMode("composite");

					savingDir = outputFold+File.separator+"STITCHED_"+substring(files[f], 0, lengthOf(files[f])-4);
					File.makeDirectory(savingDir);
					saveAs("tiff", savingDir+File.separator+"STITCHED_"+substring(files[f], 0, lengthOf(files[f])-4)+"_pos"+l+".tif"); //saving stiched image
					saveType = "CHED_";
									
					pixelTotal = pixelTotal+ (getHeight()*getWidth());
					imgInfo=getImageInfo();
					sizeLoaded = sizeLoaded + parseInt(substring(imgInfo, lastIndexOf(imgInfo, "Size:")+7, lastIndexOf(imgInfo, "MB")));
					nmrOfFilesStitched = nmrOfFilesStitched+1;
					//close("Fused");list[st]
				
					list = getFileList(positionDir);
  					for (li=0; li<list.length; li++){
  						File.delete(positionDir+File.separator+list[li]); //deleting individual files created. The original data is still intact in the .czi files
  						File.delete(positionDir);	//deletes folder
  						File.delete(savingDir_temp);	//deletes folder										
					}
				}
			}
		Tot_nmrOfTilesOpened = Tot_nmrOfTilesOpened + TotalNumTiles;
		resultArray = append(resultArray, savingDir);
		//print("result array: "+resultArray[k-1]);
		//print("result array length: "+resultArray.length);
		}
	}
//print(f);
						
if (k==0) 
	exit("Sorry, there seems to be no .czi files in \n"+dir+".");
else if (seriesCount == 1) 
	exit("Your .czi file(s) in \n"+dir+" \ndon't seem to have more than 1 tile.");

close("*");

	for (rA=0; rA<resultArray.length; rA++){
		//print(resultArray[rA]);
		//print("array length is: "+resultArray.length);
		//print("current array is: "+resultArray[rA]);
		list = getFileList(resultArray[rA]);
		for (st=0; st<list.length; st++){
  			if (matches(list[st], ".*STIT"+saveType+".*")){
  				//print(resultArray[rA]+File.separator+list[st]);
  				curTil = resultArray[rA]+File.separator+list[st];
  				open(curTil);
				//run("Color Balance...");
				run("RGB Color", "slices");
				title=getTitle();
				//print(title);
				title=substring(title,0,lengthOf(title)-13);
				//print(title);
				//print(curTil);
  				//print(list[st]);
  				close(title);
  			}
		}
		if (positionsNum>1){
			//setBatchMode(false);
			run("Images to Stack", "name=Stack title="+title+" keep");
			colNum = list.length/2;
			rowNum = list.length/colNum;
			run("Make Montage...", "columns="+colNum+" rows="+rowNum+" scale="+1/positionsNum+" border=4 font=30 label");
			saveAs("tiff", resultArray[rA]+File.separator+"Montage.tif"); //saving stiched image
			close("stack");
			if(visualize_checkBox==1){	
				setBatchMode(false);
				//open(resultArray[rA]+File.separator+"Montage.tif");
				//run("Despeckle", "stack");
     			//run("Enhance Contrast...", "saturated=0.1 process_all use");
			}	
		}
		else{
			if(visualize_checkBox==1){
				//print(resultArray[rA]);
				//print(list[0]);
				setBatchMode(false);
				//print(curTil);
				open(curTil);
			}
		}
	}


end = getTime();
Dialog.create("This is not good for Zeiss business...");
Dialog.addMessage("You just stitched "+Tot_nmrOfTilesOpened+" Tile(s) into "+nmrOfFilesStitched+" Mosaic position(s) from "+k+" files.\n \nThat's a total of "+pixelTotal+" pixels and "+sizeLoaded+" MB... \n...within "+(end-start)/1000+" seconds ("+(end-start)/1000/60+" minutes).\n \nYou know, that's Super Computer kind of stuff...clap clap!");
Dialog.show();

//function to append data to expanding arrays
function append(arr, value) {
     arr2 = newArray(arr.length+1);
     for (i=0; i<arr.length; i++)
        arr2[i] = arr[i];
     arr2[arr.length] = value;
     return arr2;
  }