Homec4science
Diffusion FlatField And Stitching Multiseries files (master) Browse

ijs-flatfield-and-stitch/BIOP_Flatfield_And_Stitch.groovy
89694b6aa336master

Last Change
Raw File

BIOP_Flatfield_And_Stitch.groovy
View Options

#@File(label="Multiposition Files Directory",style="directory") image_directory
#@String(value="Select a file for flat field correction or write none", visibility="MESSAGE") textFF
#@File(label="Flatfield Stack File") ff_file
#@int(label="Downsample Factor", style = "slider", min = 1, max = 16, stepSize = 1) downsample
#@Boolean(label="Perform Stitching") is_do_stitch
#@Boolean(label="Compute Overlap") is_compute
#@Boolean(label="Keep Original Bit Depth") is_keep_bit_depth
#@String(label="Macro Mode",choices={"Fuse and display","Write to disk"}) output_type
#@Boolean(label="Save final stitched as *.ids (if write to disk)") save_Stiched_asIDS
#@Boolean(label="Delete temporary files (if write to disk)") delete_temp_files
#@LogService log
/*
* 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: January 2018
*
* 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
*
* March 2020: Harmonized code with other version and fixes new issues with BioFormats
*
*/
// Starting the Script, we need to test whether autostitch is active in bioformats
def was_auto_stitch = Prefs.get( LociPrefs.PREF_CZI_AUTOSTITCH, false )
// In case of CZI autostitch, set the preference to false
Prefs.set( LociPrefs.PREF_CZI_AUTOSTITCH, false )
// Open the flatfield image as needed
def ff_image = ff_file.getName().matches("[Nn]one") ? null : IJ.openImage( ff_file.getPath() )
// List all files. We are focusing only on lif, czi and lsm
def all_files = image_directory.listFiles().findAll{ it.getAbsolutePath().endsWith( ".lif" ) || it.getAbsolutePath().endsWith( ".czi" ) || it.getAbsolutePath().endsWith( ".lsm" ) }
// Process each file
all_files.each{ f ->
exportAndStitch( f, ff_image, downsample )
}
//Reset CZI Autostitch Preference at the end of the script
Prefs.set( LociPrefs.PREF_CZI_AUTOSTITCH, was_auto_stitch );
/**
* Helper method to divide an image with a flatfield
* Channels are split in case there are Z stacks
*/
ImagePlus divideImages( ImagePlus image, ImagePlus flatfield ) {
// This handles the logig of whether we stitch or not
if( flatfield == null ) return image
def image_channels = ChannelSplitter.split( image )
def flatfield_channels = ChannelSplitter.split( flatfield )
def ic = new ImageCalculator();
if( image_channels.length != flatfield_channels.length ) return image
def corrected_channels = [image_channels, flatfield_channels].transpose().collect{ i, f -> def corr_imp = ic.run( "Divide create 32-bit stack", i, f ) }
if(is_keep_bit_depth) {
corrected_channels.each {
def bd = it.getBitDepth()
def i_con = new ImageConverter( it )
if(bd == 16) i_con.convertToGray16()
else if(bd==8) i_con.convertToGray8()
}
}
return RGBStackMerge.mergeChannels( corrected_channels as ImagePlus[], false )
}
void exportAndStitch(File image_file, ImagePlus ff_image, int downsample) {
// Get the base directory for the multiposition file
def base_dir = image_file.getParent()
// Use it to create a save directory
// The regular expression gets rid of the extension. Read: "a dot followed by (not a dot) at least one time (+) then the string should end ($)
def exportLif_saveDir = new File( base_dir, image_file.getName().replaceFirst("[.][^.]+\$", "") )
exportLif_saveDir.mkdir()
IJ.log(exportLif_saveDir.getAbsolutePath())
// file to export image positions
File positions = new File( exportLif_saveDir, "positions.txt" )
def stitched_saveDir = new File( exportLif_saveDir , "stitched" )
if(is_do_stitch) stitched_saveDir.mkdir()
// Get the full path of the file.
def theFile = image_file.getPath()
// lif files have their coordiante metadata stored in another unit...
// or at least bioformats has issues with them...
def corr_factor = 1.0
if(theFile.endsWith(".lif")) corr_factor= 1e6
// Prepare BioFormats Importer options
def opts = new ImporterOptions()
opts.setId( theFile )
opts.setWindowless( true )
opts.setQuiet( true )
opts.setUngroupFiles( true )
opts.clearSeries()
//set up import process
def process = new ImportProcess( opts )
process.execute()
// Get how many files we will be working on
n_series = process.getSeriesCount()
//reader belonging to the import process
def i_reader = process.getReader()
def imp_reader = new ImagePlusReader( process )
// Prepare all the metadata
def factory = new ServiceFactory()
def service = factory.getInstance( OMEXMLService.class )
MetadataRetrieve retrieve = service.asRetrieve( i_reader.getMetadataStore() )
// Pixel size is assumed to be the same for all series in the file!
// This is because getPixelsPhysicalSizeX(int) sometimes returns null
def vx = retrieve.getPixelsPhysicalSizeX(0).value()
// to re-order virtual stack and export as a single file
def ch_count = retrieve.getChannelCount(0)
def plane_count = retrieve.getPlaneCount(0)
//def timepoint_count = TO DO !
//println (ch_count)
//println ( plane_count)
//println ( timepoint_count)
def posX = []
def posY = []
def names = []
// Go through each serie within the file
(1..n_series).each{ idx ->
// Store positions and names based on metadata
posX.add( retrieve.getPlanePositionX( idx - 1, 0 ).value() )
posY.add( retrieve.getPlanePositionY( idx - 1, 0 ).value() )
names.add( retrieve.getImageName( idx - 1 )+"_"+idx+".tif" )
// If the file was not saved already, do the processing
def file_name = new File( exportLif_saveDir, names.last() )
log.info("File: "+file_name)
if ( !file_name.exists() ) {
// Open the actual image
opts.setSeriesOn( idx - 1, true )
i_reader.setSeries( idx - 1 )
def imps = imp_reader.openImagePlus()
// Because there is no grouping, there is always only one element
def imp = imps[0]
dims = imp.getDimensions()
// Do flatfield: if ff_image is null, it will return imp
imp = divideImages( imp, ff_image )
// Downsample if desired
if ( downsample > 1 ) {
def size_x = Math.round( imp.getWidth() / downsample ) as int
def size_y = Math.round( imp.getHeight() / downsample ) as int
def stack = imp.getStack()
def slices = stack.size()
def stack_down = ImageStack.create( size_x, size_x, slices, imp.getBitDepth() )
(1..slices).each{ stack_down.setProcessor( stack.getProcessor( it ).resize( size_x, size_y, true ), it ) }
imp.setStack( stack_down )
}
// Save image as TIFF
IJ.saveAs( imp, "Tiff", file_name.getAbsolutePath() )
dims = imp.getDimensions()
// Some cleanup
imp.changes=false
imp.close()
opts.setSeriesOn(idx-1, false)
}
}
i_reader.close()
// Get min in X and Y
def minX = posX.min()
def minY = posY.min()
log.warn(minX)
// open the first image again to check the dimensions
def image = IJ.openImage( new File( exportLif_saveDir, names.first() ).getAbsolutePath() )
def dims = image.getDimensions()
def dim=2
def z = ""
if(dims[3] > 1) {
dim = 3
z = ", 0.0"
}
positions.write( "#Define the number of dimensions we are working on:\n" )
positions << "dim = "+dim+"\n"
positions << "# Define the image coordinates\n"
[posX, posY, names].transpose().each{ px, py, name ->
// Ignore the merged images in LIF files
if ( !name.contains( "Merging" ) ) {
def fposX = Math.round( ( px - minX ) * corr_factor / vx / downsample )
def fposY = Math.round( ( py - minY ) * corr_factor / vx / downsample )
positions << name+"; ; ("+fposX+", "+fposY + z+")\n"
}
}
// Prepare to run the stitching
if(is_do_stitch) {
// compute overlap or not ?
compute = is_compute ? compute = "compute overlap " : ""
// Define an output directory if needed
output_dir = (output_type != "Fuse and display") ? output_dir = "output_directory=["+stitched_saveDir.getAbsolutePath()+File.separator+"]" : ""
// Do the stitching here with defined option
IJ.run("Grid/Collection stitching", "type=[Positions from file] "+
"order=[Defined by TileConfiguration] "+
"directory=["+exportLif_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=["+output_type+"] "+
output_dir);
// Save image in case output type was set to display
if ( output_type.equals( "Fuse and display" ) ) {
def final_image = IJ.getImage();
//IJ.saveAs(final_image, "Tiff", new File( stitched_saveDir, image_file.getName() + suffix + "-fused.tif" ).getAbsolutePath() )
IJ.saveAs(final_image, "Tiff", new File( stitched_saveDir, image_file.getName() + "-fused.tif" ).getAbsolutePath() )
final_image.close();
log.info( "Fused image Saved" )
} else if (save_Stiched_asIDS) {
// re-open exported indiviual stitched plane as a virtual stack
def impV = FolderOpener.open(stitched_saveDir.getAbsolutePath(), "virtual");
// re-order using ch_count and total plane number ! yeah doesn't work (yet) with time-series
def reordered_impV = HyperStackConverter.toHyperStack(impV, ch_count, (plane_count/ch_count as int), 1, "Color");
//save as ids file
def idsFile_path = new File(image_directory, image_file.getName()+"_stiched.ids").getAbsolutePath()
IJ.run(reordered_impV, "Bio-Formats Exporter", "save=["+idsFile_path+"]");
// delete individual stitched images if asked to do so
if (delete_temp_files){
stitched_saveDir.listFiles().each{ it.delete() }
stitched_saveDir.delete()
}
}
}
}
/**
* Imports below
*/
import loci.formats.ImageReader
import loci.formats.MetadataTools
import loci.formats.meta.IMetadata
import loci.formats.meta.MetadataRetrieve
import loci.common.services.ServiceFactory
import loci.formats.services.OMEXMLService
import loci.plugins.in.ImagePlusReader
import loci.plugins.in.ImporterOptions
import loci.plugins.in.ImportProcess
import loci.plugins.util.LociPrefs
import ij.IJ
import ij.Prefs
import ij.ImagePlus
import ij.ImageStack
import ij.plugin.ChannelSplitter
import ij.plugin.RGBStackMerge
import ij.plugin.ImageCalculator
import ij.process.ImageConverter
import ij.plugin.FolderOpener
import ij.plugin.HyperStackConverter
c4science · Help