diff --git a/BIOP_Flatfield_And_Stitch.groovy b/BIOP_Flatfield_And_Stitch.groovy index 7fa2895..3abb802 100644 --- a/BIOP_Flatfield_And_Stitch.groovy +++ b/BIOP_Flatfield_And_Stitch.groovy @@ -1,343 +1,339 @@ #@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="Delete temporary files") delete_temp_files +#@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("[.][^.]+\$", "") ) - def stitched_saveDir = new File( exportLif_saveDir , "stitched" ) - + exportLif_saveDir.mkdir() IJ.log(exportLif_saveDir.getAbsolutePath()) + + // file to export image positions + File positions = new File( exportLif_saveDir, "positions.txt" ) - exportLif_saveDir.mkdir() + 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 = retrieve.getTimestampAnnotationCount() - + //def timepoint_count = TO DO ! //println (ch_count) //println ( plane_count) //println ( timepoint_count) def posX = [] def posY = [] def names = [] - - File positions = new File( exportLif_saveDir, "positions.txt" ) - - - // Go through each file + // 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 - compute = "" - if(is_compute) { compute = "compute overlap " } - if(is_do_stitch) { - // Define an Output directory if needed + // 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() + 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() + } } - } - - if (output_type != "Fuse and display"){ - // re-open exported indiviual stitched plane as a virtual stack - def impV = FolderOpener.open(stitched_saveDir.getAbsolutePath(), "virtual"); - //impV.show() - // 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"); - //reordered_impV.show() - //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_fileList = stitched_saveDir.listFiles(); - stitched_fileList.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 \ No newline at end of file