[Flatfield correction followed by stitching] Performs a flatfield correction and stitching from Raw LIF, CZI and LSM files using plugins in Fiji #Fiji #Stitching #Groovy #Flatfield #LIF #LSM #CZI #BIOP

FlatField_And_Stitch.groovy

/*
 * Flatfield correction followed by stitching
 * 
 * It is often the case that tiled microscopy images contain artifacts due to uneven illumination.
 * This is due to the shape of the illumination, the quality of the optics and the size of the field of view, but
 * losses of up to 50% can be seen between the periphery and the center of a field of view. This loss usually follows
 * a parabolic function, but is not necessarily centered, causing artifacts when stitching tiled acquisitions
 * 
 * While many methods exist to compensate for this, the simplest consists in acquiring a 'flat field' image by
 * taking a homogeneous sample (Chroma slides or even better, dye solutions [https://www.ncbi.nlm.nih.gov/pubmed/18173639]) and acquiring it for each channel of the dataset
 * 
 * This image can then be used to divide the original one and thus help 'flatten' the illumination. 
 * 
 * Furthermore, excellent tools already exist to perform image stitching such as the Stitch Grid/Collection Plugin
 * 
 * By leveraging the scriptability of that plugin, this script 
 * 1. Flattens the intensities based on a user-provided Flat Field image stack (1 slice per channel)
 * 2. Saves a tiff (and possiby downsampled) version of the images
 * 3. Reads the XY Position of each tile to create a TileConfiguration file readable by the Stitch Grid/Collection Plugin
 * 4. Eventually calls said plugin
 * 
 * This script has been tested and is known to work with .lif, .czi and .lsm multiseries files.
 * These are due to the fact they are the most common formats to come out of our microscopes. 
 * 
 * Authors: Olivier Burri, Romain Guiet, Joao Firmino
 * BioImaging & Optics Platform (BIOP)
 * 
 * Last modification: November 2020
 * 
 * Update List
 * January 2018:  Added the possibility to work with tiled CZI files. Make sure that you disable 'Automatically Stitch Tiled Images'
 * 				  is checked under Plugins > Bio-Formats > Bio-Formats Plugin Configuration > Formats > Zeiss CZI
 * 				  User no longer needs to check the box, this is done automarically
 * November 2020: Updated how to get position units using Bioformats.
 * 				  Added hack to flip and swap coordinates because of some lif files 
 * 				  removed corr_factor
 * 				  Refactoring in more Groovy-esque way
 * 
 */

#@ File image_directory ( label = "Multiposition Files Directory", style = "directory" )
#@ File ff_file ( label = "Flatfield Stack File" )
#@ Integer downsample ( label = "Downsample Factor", style = "slider", min = 1, max = 16, stepSize = 1, value = 1 ) 
#@ Boolean is_do_stitch ( label = "Perform Stitching" ) 
#@ Boolean is_compute ( label = "Compute Overlap" ) 
#@ Boolean is_not_flatfield( label = "Don't Perform Flatfield" ) 
#@ Boolean is_keep_bit_depth( label = "Keep Original Bit Depth" )
#@ Boolean is_flip_swap ( label = "Needed for some LIF Files: Flip and Swap Axes?", value = false ) 
#@ LogService log

// In case of CZI autostitch, set the preference to false
def was_auto_stitch = Prefs.get( LociPrefs.PREF_CZI_AUTOSTITCH, false )

Prefs.set( LociPrefs.PREF_CZI_AUTOSTITCH, false )

def ff_image = null

if (!is_not_flatfield) {
	ff_image = IJ.openImage(ff_file.getPath())
}

def all_files = image_directory.listFiles().findAll{ it.getAbsolutePath().endsWith(".lif") || it.getAbsolutePath().endsWith(".czi") || it.getAbsolutePath().endsWith(".lsm")}

all_files.each{ f ->
	exportAndStitch(f, ff_image, downsample)
}

//Reset CZI Autostitch Preference
Prefs.set(LociPrefs.PREF_CZI_AUTOSTITCH, was_auto_stitch);


/* 
 *  This divides the image array (each index == 1 channel) with a flatfield image
*/
ImagePlus divideImages( ImagePlus image, ImagePlus flatfield ) {

	// split them 
	def imageC = ChannelSplitter.split( image )
	def flatfieldC = ChannelSplitter.split( flatfield )
	
	if( imageC.length != flatfieldC.length ) return null
		
	def ic = new ImageCalculator();
		
	def flattenedC = [imageC, flatfieldC].transpose().collect{ img, ff ->		
   			
   		def temp_image = ic.run( "Divide create 32-bit stack", img, ff );
   			
   		if(is_keep_bit_depth) {
   			def i_c = new ImageConverter( temp_image )
   			def bd = img.getBitDepth()
   			if( bd == 16 ) i_c.convertToGray16()
   			else if( bd==8 ) i_c.convertToGray8()
		}

		return temp_image
	}
		
	return RGBStackMerge.mergeChannels(flattenedC, false);
}

void exportAndStitch( File image_file, ImagePlus ff_image, Integer downsample ) {
	
	// Get the base directory for the multiposition file
	def base_dir = image_file.getParent()
	
	// Use it to create  a save directory
	def saveDir = new File( base_dir, image_file.getName().replaceFirst( "[.][^.]+\$", "" ) + File.separator )
	IJ.log( "Save Directory: " + saveDir.getAbsolutePath() )
	saveDir.mkdir()
	
	// Get the full path of the file. 
	def theFile = image_file.getPath()
	
	// Setup import process
	def opts = new ImporterOptions()
	opts.setId( theFile )
	opts.setWindowless( true )
	opts.setQuiet( true )
	opts.setUngroupFiles( true )
	opts.clearSeries()
	
	def process = new ImportProcess( opts )
	process.execute()
	nseries = process.getSeriesCount()
	
	//reader belonging to the import process
	def i_reader = process.getReader()
	def impReader = new ImagePlusReader( process )
	
	// Get all the metadata
	def factory = new ServiceFactory()
	def service = factory.getInstance( OMEXMLService.class )
	
	// And the retreive function will help us get all the metadata back
	def retrieve = service.asRetrieve( i_reader.getMetadataStore() )
	
	def posX  = []
	def posY  = []
	def names = []
	
	def vx = retrieve.getPixelsPhysicalSizeX( 0 ).value( UNITS.MICROMETER ).doubleValue()
	
	all_dims = null
	for (int i=0; i<nseries; i++) {
		// for (int i=0; i<2; i++) {
		def px = retrieve.getPlanePositionX( i, 0 ).value( UNITS.MICROMETER ).doubleValue()
		def py = retrieve.getPlanePositionY( i, 0 ).value( UNITS.MICROMETER ).doubleValue()

		if ( is_flip_swap ) {
			// Flip
			px *= -1
			py *= -1
			
			// Swap
			def temp = px
			px = py
			py = temp
		}
		
		posX.add( px )
		posY.add( py )
		
		def i_name = retrieve.getImageName( i ) + "_" + (i+1) + ".tif"
		
		names.add( i_name )
		
		opts.setSeriesOn( i, true )
		i_reader.setSeries( i )
		
		//read and process all images in series
		def imp = impReader.openImagePlus()[0]
		 
		def bd= imp.getBitDepth()


		// Perform the flatfield, if there is a file
		
		def ff_imp = imp	
		if( ff_image != null ) {
			ff_imp = divideImages( imp, ff_image )
		}

		// Downsample if requested
		small_ff_imp = ff_imp;
	
		if ( downsample > 1 ) {
		  	IJ.run( ff_imp, "Scale...", "x=" + ( 1.0 / downsample ) + " y=" + ( 1.0 / downsample ) + " z=" + ( 1.0 / downsample ) + " create interpolation=Bilinear average" )
		  	small_ff_imp = IJ.getImage()
		}

		// Save image as TIFF
		IJ.saveAs( small_ff_imp, "Tiff", new File( saveDir, i_name ).getAbsolutePath() )
		
		
		// Some cleanup	  
		imp.changes=false
		imp.close()
		ff_imp.close()
		all_dims = small_ff_imp.getDimensions()
		small_ff_imp.close()
		opts.setSeriesOn( i, false )
	}
	
	i_reader.close()
	
	
	// Get min in X and Y
	def minX = Collections.min( posX )
	def minY = Collections.min( posY )
	
	def dim=2
	def z = ""
	if(all_dims[3] > 1) {
		dim = 3
		z   = ", 0.0"
	}

	// Prepare a positions file for Grid/Collection Stitching
	def out = new File( saveDir, "positions.txt" )
	out.write( "#Define the number of dimensions we are working on:\n" )
	out << "dim = " << dim << "\n"
	out << "# Define the image coordinates" << "\n"

	[ names, posX, posY ].transpose().each{ name, x, y ->
		if ( !name.contains( "Merging" ) ) {
			def fposX = Math.round( ( ( x - minX ) / vx / downsample ) )
			def fposY = Math.round( ( ( y - minY ) / vx / downsample ) )		
			out << name + ";	;	(" + fposX + ", " + fposY + z+ ")" << "\n"
		}
	}

	// Call Grid/Collection Stitching
	compute = ""
	if( is_compute ) { compute = "compute overlap " }
	if( is_do_stitch ) {
		// run the stitcher
		IJ.run( "Grid/Collection stitching", "type=[Positions from file] "+
		"order=[Defined by TileConfiguration] "+
		"directory=[" + saveDir.getAbsolutePath() + "] "+
		"layout_file=positions.txt "+
		"fusion_method=[Linear Blending] "+
		"regression_threshold=0.30 "+
		"max/avg_displacement_threshold=2.50 "+
		"absolute_displacement_threshold=3.50 "+
		compute +
		"computation_parameters=[Save memory (but be slower)] "+
		"image_output=[Fuse and display]")
		
		final_image = IJ.getImage()
		IJ.saveAs( final_image, "Tiff", saveDir.getAbsolutePath()+File.separator+image_file.getName()+"-fused.tif" )
		final_image.close()
	}
}

// Imports
import ij.Prefs
import ij.IJ
import ij.ImagePlus
import ij.plugin.ChannelSplitter
import ij.plugin.ImageCalculator
import ij.plugin.RGBStackMerge
import ij.process.ImageConverter

import ome.units.UNITS

import loci.formats.ImageReader
import loci.formats.meta.IMetadata
import loci.formats.meta.MetadataRetrieve
import loci.formats.MetadataTools
import loci.formats.services.OMEXMLService
import loci.common.services.ServiceFactory
import loci.plugins.in.*
import loci.plugins.util.LociPrefs