Created
June 16, 2025 14:56
-
-
Save romainGuiet/7a32a36cd1e602fcf99e7e479c45f1b1 to your computer and use it in GitHub Desktop.
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
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #@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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment