Homec4science
Diffusion FlatField And Stitching Multiseries files (master) Browse

ijs-flatfield-and-stitch/BIOP_Flatfield_And_Stitch.groovy
c6c2850d1890master

Last Change
Raw File

BIOP_Flatfield_And_Stitch.groovy
View Options

//@File(label="Multiposition Files Directory",style="directory") image_directory
//@File(label="Flatfield Stack File") ff_file
//@int(label="Downsample Factor", value=1, style = "slider", min = "1", max = "16", stepSize = "1") downsample
//@Boolean(label="Perform Stitching") is_do_stitch
//@Boolean(label="Compute Overlap") is_compute
//@Boolean(label="Don't Perform Flatfield") is_not_flatfield
//@Boolean(label="Keep Original Bit Depth") is_keep_bit_depth
//@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
*
*
*/
import ij.IJ
import loci.formats.ImageReader
import loci.formats.meta.IMetadata
import loci.formats.meta.MetadataRetrieve
import loci.formats.MetadataTools
import loci.common.services.ServiceFactory
import loci.formats.services.OMEXMLService
import java.util.Arrays
import java.util.Collections
import org.apache.commons.lang.ArrayUtils
import ij.plugin.ChannelSplitter
import ij.plugin.RGBStackMerge
import ij.ImagePlus
import ij.plugin.ImageCalculator
import loci.plugins.in.ImagePlusReader
import loci.plugins.in.ImporterOptions
import loci.plugins.in.ImportProcess
import java.io.FilenameFilter
import ij.process.ImageConverter
import ij.gui.WaitForUserDialog
import loci.plugins.util.LociPrefs
import ij.Prefs
// 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)
/*
* flatfield file opening and management
* We need to open the Flatfield image and split the channels.
*/
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) {
if(image.length == flatfield.length) {
def ic = new ImageCalculator();
image.eachWithIndex{ img, idx ->
image[idx] = ic.run("Divide create 32-bit stack", img, flatfield[idx]);
if(is_keep_bit_depth) {
def i_c = new ImageConverter(image[idx])
def bd = img.getBitDepth()
if(bd == 16)
i_c.convertToGray16()
else if(bd==8)
i_c.convertToGray8()
}
}
return RGBStackMerge.mergeChannels(image, 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
def saveDir = new File(base_dir+File.separator+image_file.getName().replaceFirst("[.][^.]+\$", "")+File.separator)
IJ.log(saveDir.getAbsolutePath())
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
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()
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
MetadataRetrieve retrieve = service.asRetrieve(i_reader.getMetadataStore())
double[] posX = new double[i_reader.getSeriesCount()]
double[] posY = new double[i_reader.getSeriesCount()]
String[] names = new String[i_reader.getSeriesCount()]
vx = retrieve.getPixelsPhysicalSizeX(0).value()
all_dims = null
for (int i=0; i<nseries; i++) {
// for (int i=0; i<2; i++) {
posX[i] = (retrieve.getPlanePositionX( i, 0 ).value())
posY[i] = (retrieve.getPlanePositionY( i, 0 ).value())
// print("X:"+posX[i])
// print("Y:"+posY[i])
names[i] = retrieve.getImageName(i)
opts.setSeriesOn(i,true)
i_reader.setSeries(i)
//read and process all images in series
imps = impReader.openImagePlus()
//IJ.log("Length "+imps.length);
imp=imps[0]
bd= imp.getBitDepth()
ff_imp = imp
// Perform the flatfield, if there is a file
if( ff_image != null ) {
ff_imp = divideImages(ChannelSplitter.split(imp),ChannelSplitter.split(ff_image))
}
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", saveDir.getAbsolutePath()+File.separator+names[i]+"_"+(i+1)+".tif")
// 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
List lX = Arrays.asList(ArrayUtils.toObject(posX));
List lY = Arrays.asList(ArrayUtils.toObject(posY));
minX = Collections.min(lX)
minY = Collections.min(lY)
dim=2
z = ""
if(all_dims[3] > 1) {
dim = 3
z = ", 0.0"
}
PrintWriter out = new PrintWriter(saveDir.getAbsolutePath()+File.separator+"positions.txt")
out.println("#Define the number of dimensions we are working on:")
out.println("dim = "+dim)
out.println("# Define the image coordinates")
for (i=0; i<posX.length; i++) {
if (names[i].contains("Merging")) continue
fposX = Math.round(((posX[i]-minX)*corr_factor/vx/downsample))
fposY = Math.round(((posY[i]-minY)*corr_factor/vx/downsample))
out.println(names[i]+"_"+(i+1)+".tif; ; ("+fposX+", "+fposY + z+")")
}
out.close()
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()
}
}
c4science · Help