diff --git a/BIOP_Operetta_Import.groovy b/BIOP_Operetta_Import.groovy
index 7310530..fc6dbb3 100644
--- a/BIOP_Operetta_Import.groovy
+++ b/BIOP_Operetta_Import.groovy
@@ -1,522 +1,562 @@
//@File(label="Select your directory with your exported images", style="directory") theDir
//@Integer(label="Resize Factor", value=1) resize
-//@Boolean(label="Only Process Selected Wells", value=false) is_select_wells
-
+//@Boolean(label="Only Process Selected Wells", value=false) is_select_wells //
+//@String(label="X Y W H of box to extract", value="") str_xywh
// ---------------- DESCRIPTION ----------------- //
/*
* PERKINELMER OPERETTA STITCHER
* v3.1, October 2017
* This tool allows for the tiling and saving of tiffs exported
* with the Operetta Symphony software to be viewed and processed
* with Fiji
*
* So far this tool will tile all fields in each well
* The output is one tiled hyperstack per well (CZT + fields )
* We benefit from Gpars for parallel processing, so there are a
* few dependent libraries.
* See https://c4science.ch/w/bioimaging_and_optics_platform_biop/image-processing/imagej_tools/perkinelmer-stitching/
* For dependencies and instructions
*
* Authors: Olivier Burri, Romain Guiet
* BioImaging and Optics Platform (BIOP)
* Ecole Polytechnique Fédérale de Lausanne
*
* Copyright 2017 Olivier Burri, Romain Guiet
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
// ------------------ IMPORTS ------------------ //
import System.*
import groovy.util.XmlSlurper
import ij.*
import ij.gui.*
import ij.plugin.*
import ij.process.*
import ij.measure.Calibration
// Play with parallel stuff
import groovyx.gpars.GParsPool
import groovyx.gpars.GParsExecutorsPool
// Number converter
import java.text.DecimalFormat
// Log
import ij.text.TextPanel
import ij.text.TextWindow
import groovy.transform.Synchronized
import java.awt.Font;
+// GUI goodness
import groovy.swing.SwingBuilder
import javax.swing.*
import java.awt.*
// ------------------ SCRIPT ------------------ //
def pe = new PerkinElmerStitcher()
pe.isTestMode(false)
+if (str_xywh.size() > 0) {
+ def coordinates = str_xywh.tokenize(' ').collect{ it.toInteger() }
+
+ pe.setRoi(new Roi(coordinates[0], coordinates[1], coordinates[2], coordinates[3]))
+}
+
+
+
// Parse the file
def parser = new Timer("Parser")
parser.tic()
pe.parseXML(theDir)
parser.toc()
if(is_select_wells) {
def selproc = new Timer("Processing selected")
selproc.tic()
pe.selectWellsGUI(resize)
selproc.toc()
} else {
// Process all images in parallel, see IJ log
def processor = new Timer("Processing")
processor.tic()
pe.process(resize)
processor.toc()
}
// ----------------- CLASSES ----------------- //
/*
* The big boy that does everything
*/
class PerkinElmerStitcher {
def xml_name = "Index.idx.xml"
def path
def savePath
def ims
def voxelSize
def imageSize
def is_test = false
+
+ def roi = null
+
String experimentName
HashSet rc_indexes
HashSet rc_selected
def cal
void isTestMode(boolean is_test) {
this.is_test = is_test
}
void parseXML(File dir) {
path = dir.getAbsolutePath()+"\\"
// Prepare save directory
savePath = path+"output\\"
new File(savePath).mkdir()
// Parsing XML file
String xml_file = path.concat(xml_name)
def experiment = new XmlSlurper().parse(xml_file)
ims = new ArrayList()
// Get the experiment name
experimentName = experiment.Plates.Plate.PlateID
// from the experiement file, create list of informations for each Image
experiment.Images.Image.each{
def im = new Image()
im.rowcol.add(it.Row.toInteger())
im.rowcol.add(it.Col.toInteger())
im.field = it.FieldID.toInteger()
im.channel = it.ChannelID.toInteger()
im.slice = it.PlaneID.toInteger()
im.timepoint= it.TimepointID.toInteger()
im.posx = it.PositionX.toFloat()
im.posy = it.PositionY.toFloat()
im.posz = it.PositionZ.toFloat()
im.toffset = it.MeasurementTimeOffset.toFloat()
im.image = it.URL.text()
ims.add(im)
}
// Get the size of the image
this.imageSize = [x: experiment.Images.Image[0].ImageSizeX.toInteger(), y: experiment.Images.Image[0].ImageSizeY.toInteger() ]
//Get the voxel size
this.voxelSize = experiment.Images.Image[0].ImageResolutionX.toFloat()
//Get the calibration
this.cal = computeCalibration()
// It's not necessarly consecutive wells that are acquired
// so we have to check the indexes of each existing one !
// fortunatly a function "HasSet" that does exactly this kind of job for "Set"
// you give list/collection and return the list of unique !
this.rc_indexes = new HashSet(ims.rowcol)
this.rc_selected = rc_indexes
}
+ void setRoi(Roi roi) {
+ this.roi = roi;
+ }
+
void process(int resize) {
// Get the extents of this Dataset, all other dimensions are assumed to be identical per well
def c_xtents = [ start: ims.min{ it.channel }.channel,
end: ims.max{ it.channel }.channel
]
def t_xtents = [ start: ims.min{ it.timepoint }.timepoint,
end: ims.max{ it.timepoint }.timepoint
]
def z_xtents = [ start: ims.min{ it.slice }.slice,
end: ims.max{ it.slice }.slice
]
def f_xtents = [ start: ims.min{ it.field }.field,
end: ims.max{ it.field }.field
]
// Get minimum extent of position in xy
def minxy = [ x: ims.min { it.posx }.posx, y: ims.min{ it.posy }.posy ]
def maxxy = [ x: ims.max { it.posx }.posx, y: ims.max{ it.posy }.posy ]
// Get the full size of the TILED image, made of individual TILEs
def tiledXYsize = [x: Math.round((maxxy['x'] - minxy['x']) / voxelSize), y: Math.round((maxxy['y'] - minxy['y']) / voxelSize) ]
tiledXYsize['x'] += imageSize['x']
tiledXYsize['y'] += imageSize['y']
tiledXYsize['x'] /= resize
tiledXYsize['y'] /= resize
//Test mode only tried to process 2 wells
if (is_test) {
rc_selected = rc_selected.take(2)
}
// Build list with all unique combinations of channels, times, and z and work for each field, that way we are parallel for anything in a given well
// Thanks, Groovy!
def allCZT = new ArrayList()
(c_xtents['start']..c_xtents['end']).each{ c -> (z_xtents['start']..z_xtents['end']).each{ z -> (t_xtents['start']..t_xtents['end']).each{ t -> allCZT.add(new CZT(c,z,t)) } } }
// Get max memory, but be a bit conservative
def maxMemory = IJ.maxMemory() / 1e9 * 0.90;
// Calculate the size on RAM needed for the data...
// A TILED Well image is (in CZTF)
def czt_tiled_stack = 16 * allCZT.size() * tiledXYsize['x'] * tiledXYsize['y'] / 8 / 1e9
IJ.log("We will be processing "+rc_selected.size()+" images, each "+(new DecimalFormat("##.#").format(czt_tiled_stack))+" GB...")
// Get Number of Cores
def cores = Runtime.getRuntime().availableProcessors();
// Define number of images we can work on in parallel
def nImages = (int) Math.floor(maxMemory/czt_tiled_stack) >= (cores-3) ? (cores-3) : (int) Math.floor(maxMemory/czt_tiled_stack)
// Inform the user
IJ.log("One Tiled CZT image stack is "+ (new DecimalFormat("##.##").format(czt_tiled_stack)) + " GB.\n"+
"You have "+ (new DecimalFormat("##.##").format(maxMemory) ) + " GB of RAM\n" +
"We will try to work on "+nImages+" images in parallel to process your data, totalling "+(10*nImages)+" Threads")
print("\n\n SELECTED:"+rc_selected)
// Log window to show concurrent processing
def log = new TxtLog(rc_selected, allCZT.size())
// ExecutorsPool is less optimized than GParsPool but this way we can nest calls :)
GParsExecutorsPool.withPool(nImages) {
rc_selected.eachWithIndexParallel{ rc_idx, i ->
// find all images matching this row and column, which will be a stack
// So channels, timepoints, slices and fields
print("\nProcessing field "+rc_idx+", (Active Threads: "+Thread.activeCount()+")")
// Prepare the hyperstack that will contain the full row-column data
- def largeField = ImageStack.create((int) tiledXYsize['x'], (int) tiledXYsize['y'], allCZT.size(), 16 )
-
+ def largeField
+ if(roi != null) {
+ def bounds = roi.getBounds()
+ largeField = ImageStack.create((int)bounds.width, (int)bounds.height, allCZT.size(), 16 )
+ } else {
+ largeField = ImageStack.create((int) tiledXYsize['x'], (int) tiledXYsize['y'], allCZT.size(), 16 )
+ }
// Prepare final ImagePlus here so we can access the getStackIndex function
final def rcImage = new ImagePlus("FinalImage", largeField);
rcImage.setDimensions(c_xtents['end'], z_xtents['end'], t_xtents['end']+1)
// Image name
// def czt_tiled_stack_name = experimentName+" - R"+IJ.pad(rc_idx[0],2)+"-C"+IJ.pad(rc_idx[1],2)
def czt_tiled_stack_name = experimentName+" - R"+IJ.pad(rc_idx[0],2)+"-C"+IJ.pad(rc_idx[1],2)
if(resize != 1) {
czt_tiled_stack_name+="-Resized "+resize
}
// Process slices in parallel, as what is really slow is read times for each tiff
// The Oli Rule of thumb implementation: zipped tiffs take 10 times longer to open than uncompressed tiffs
// So in terms of disk bandwidth, we have room for 10 parallel reads. Woo!
GParsPool.withPool(10) {
allCZT.eachParallel { czt ->
// Make slice that will contain all we will be stitching
ImageProcessor slice_ip = new ShortProcessor((int) tiledXYsize['x'], (int) tiledXYsize['y'])
// Work on each field
for(int f = f_xtents['start']; f<= f_xtents['end']; f++) {
// Get the image matching this CZT
def current_image = ims.find { it.rowcol == rc_idx && it.channel == czt.getC() && it.timepoint == czt.getT() && it.field == f && it.slice == czt.getZ()}
if(current_image != null) {
// Get the filename
def imageName = current_image.image.toString()
// Open the image
def current_imp = IJ.openImage(path + imageName)
// Compute Field Position
float current_x = current_image.posx.toFloat() - minxy['x']
float current_y = maxxy['y'] - (current_image.posy.toFloat())
// in Pixels
int current_x_pixel = current_x / voxelSize
int current_y_pixel = current_y / voxelSize
// Had issues with some being null once in a while...? concurrency issue of IJ.openImage()?
if(current_imp != null) {
// resize the image as requested and add it to the large slice
def current_ip = current_imp.getProcessor().resize((int) (imageSize['x'] / resize))
+
slice_ip.copyBits(current_ip, (int) (current_x_pixel / resize), (int) (current_y_pixel / resize), Blitter.COPY)
current_imp.close()
}
+
} else {
IJ.log("No Image at c:"+czt.getC()+" z:"+czt.getZ()+" t:"+czt.getT()+" for field:" +f)
}
}
// Now add this image to the hyperstack
def stack_position = rcImage.getStackIndex(czt.getC(), czt.getZ(), czt.getT()+1)
log.update(rc_idx, stack_position)
+ // if a ROI was defined, crop it before adding it
+ if (roi != null) {
+ slice_ip.setRoi(roi)
+ slice_ip = slice_ip.crop()
+ }
+
// We have the position, we can now place the data
largeField.setProcessor(slice_ip, stack_position)
}
}
// Machen-Sie-pretty
def finalstackImp = HyperStackConverter.toHyperStack(rcImage, c_xtents['end'], z_xtents['end'], t_xtents['end']+1, "xyczt", "Composite")
-
- // Set Calibration
- finalstackImp.setCalibration(this.cal)
+
+ if(resize != 1) {
+ def newCal = this.cal
+ newCal.pixelWidth *= resize
+ newCal.pixelHeight *= resize
+ finalstackImp.setCalibration(newCal)
+
+ } else {
+ // Set Calibration
+ finalstackImp.setCalibration(this.cal)
// Save
IJ.saveAs(finalstackImp, "Tiff", savePath+czt_tiled_stack_name+".tif");
rcImage.close()
}
}
}
+ }
/*
* xy size is straightformward but time and Z are not stored as intervals but absolute values
* So we need to compute their values from two subsequent frames or slices
*/
Calibration computeCalibration() {
def z_xtents = (ims.min{ it.slice }.slice)..(ims.max{ it.slice }.slice)
def t_xtents = (ims.min{ it.timepoint }.timepoint)..(ims.max{ it.timepoint }.timepoint)
// Need to compute voxelDepth
def voxelDepth = 0.0
if(z_xtents.size() > 1) {
def z1_image = ims.find { it.rowcol == ims[0].rowcol && it.channel == ims[0].channel && it.timepoint == ims[0].timepoint && it.slice == z_xtents[0] }
def z2_image = ims.find { it.rowcol == ims[0].rowcol && it.channel == ims[0].channel && it.timepoint == ims[0].timepoint && it.slice == z_xtents[1] }
voxelDepth = z2_image.posz - z1_image.posz
}
// Need to compute frameInterval
def timeDelta = 1.0
if(t_xtents.size() > 1) {
def t1_image = ims.find { it.rowcol == ims[0].rowcol && it.channel == ims[0].channel && it.timepoint == t_xtents[0] && it.slice == ims[0].slice }
def t2_image = ims.find { it.rowcol == ims[0].rowcol && it.channel == ims[0].channel && it.timepoint == t_xtents[1] && it.slice == ims[0].slice }
timeDelta = t2_image.toffset - t1_image.toffset
}
def cal = new Calibration()
cal.pixelWidth = voxelSize * 1e6
cal.pixelHeight = voxelSize * 1e6
cal.pixelDepth = voxelDepth * 1e6
cal.setUnit("um")
cal.frameInterval = (double) timeDelta
cal.setTimeUnit("s")
return cal
}
/*
* GUI for selecting wells in case this is requested
*/
Boolean selectWellsGUI(int resize) {
+
def peGUI = new SwingBuilder()
+
def positionList = {
peGUI.panel() {
scrollPane() {
list(id: "wells",
listData: rc_indexes,
selectionMode: ListSelectionModel.MULTIPLE_INTERVAL_SELECTION,
preferredSize: [150, 300]
)
}
+
}
}
def myframe = peGUI.frame(title : 'Choose Wells',
location : [100, 400],
size : [200, 300],
defaultCloseOperation : WindowConstants.DISPOSE_ON_CLOSE,
) {
panel() {
boxLayout(axis : BoxLayout.Y_AXIS)
label(text : 'Select multiple with Shift or Ctrlt',
horizontalAlignment : JLabel.CENTER
)
positionList()
button(text : 'Run',
horizontalAlignment : JLabel.CENTER,
actionPerformed : { act ->
this.rc_selected = new HashSet(peGUI.wells.getSelectedIndices().collect{ val -> rc_indexes[val] })
IJ.log("Processing Selection")
this.process(resize)
dispose()
} )
}
}
myframe.setVisible(true)
}
}
/*
* Time class to 'tic-toc' a few steps and check time spent.
*/
class Timer{
Long startTime
Long endTime
def name
public Timer(String name){
this.name = name
}
public void tic(){
this.startTime = System.nanoTime()
}
public void toc(){
this.endTime = System.nanoTime()
IJ.log("'"+name+"' took : "+((endTime-startTime)/1e9)+" s")
}
}
/*
* Image class containing important imformation about each image file
*/
class Image {
ArrayList rowcol = new ArrayList(2)
int field
int channel
int slice
int timepoint
float posx
float posy
float posz
float toffset
String image
}
/*
* Small class to store the CZT Indexes
*/
class CZT {
int c
int z
int t
CZT(int c, int z, int t) {
this.c = c
this.z = z
this.t = t
}
int getC() { return c }
int getZ() { return z }
int getT() { return t }
}
/*
* Text log to show progress of all the fun things happening during execution
*/
class TxtLog {
def rc_idx
TextWindow tw
TextPanel tp
def logList
TxtLog(HashSet rc_idx, int nSeries) {
tw = WindowManager.getWindow("Operetta Log");
if( tw == null) tw = new TextWindow("Operetta Log", "Field\tSeries", "", 300, 800)
tp = tw.getTextPanel()
tp.clear()
tp.setFont(new Font("monospaced", Font.PLAIN, 12), true)
this.rc_idx = rc_idx
logList = Collections.synchronizedList(new ArrayList(rc_idx.size()))
rc_idx.each{
def serTxt = ""
(1..nSeries).each{ serTxt+="-" }
logList.add(serTxt)
tp.append(it.toString()+":\t"+serTxt)
}
}
@Synchronized
void update(ArrayList rc_i, int pos) {
def idx = (int) rc_idx.findIndexOf{rc_i == it}
def ser = logList.get(idx)
def c = ser.toCharArray()
c[pos-1] ='+'
logList.set(idx, String.valueOf(c))
tp.setLine(idx, rc_i.toString()+":\t"+logList.get(idx))
tp.updateDisplay()
}
-}
+}
\ No newline at end of file