port the script for ACAMAS

ImageJ/Load_MM_files_Virtual_.py

+#@ File (label="Root directory of image", style="directory") rootdir
+#@ String (label="file pattern", description="Give the file pattern mydata_{z:d}_{x:003d}_{filter}") filepattern
+#@ ConvertService convertService
+
+"""
+A way to load any series of .tif in several folders as a virtualstack.
+"""
+
+import os
+from ij import ImagePlus, VirtualStack, IJ  
+from net.imglib2.img.display.imagej import ImageJFunctions as IL
+import sys
+
+if __name__ == '__main__':
+    # Add basedir to path
+    basedir = os.path.dirname(__file__)
+    sys.path.append(basedir)
+    # Load the image sequence loader
+    from utils.ImagesLoaders import load_image_sequence
+
+    #filepath = '/home/hugo/Documents/Boulo/ligne_DISCO/20200160/Pn_long_100x_1530tiles_1/'
+    filepath = rootdir.absolutePath
+    #path_pattern = 'roi6_tile{tile:d}/img_000000000_{filter}_{z:003d}.tif'
+    imp, fn = load_image_sequence(filepath, filepattern, return_filenames=True)
+
+    imp.show()

ImageJ/Test_imgJ2_ops.py

+#@ OpService ops
+#@ Dataset data
+
+# to run this tutorial run 'file->Open Samples->Confocal Series' and make sure that
+# confocal-series.tif is the active image
+
+
+from net.imglib2.util import Intervals
+from net.imagej.axis import Axes
+from net.imglib2.img.display.imagej import ImageJFunctions as IL
+from net.imglib2.algorithm.math.ImgMath import computeIntoFloat, sub, div, GT, LT, IF, THEN, ELSE, minimum
+from net.imglib2.view import Views 
+from net.imglib2.converter import Converters
+from net.imglib2.img.array import ArrayImgs  
+from net.imglib2.view import Views 
+
+from ij import IJ, ImagePlus
+
+# first take a look at the size and type of each dimension
+for d in range(data.numDimensions()):
+	print "axis d: type: "+str(data.axis(d).type())+" length: "+str(data.dimension(d))
+
+img=data.getImgPlus()
+print(img.numDimensions())
+print(img.dimensionIndex(Axes.CHANNEL), img.dimensionIndex(Axes.Z), img.dimensionIndex(Axes.TIME))
+
+xLen = data.dimension(data.dimensionIndex(Axes.X))
+yLen = data.dimension(data.dimensionIndex(Axes.Y))
+zLen = data.dimension(data.dimensionIndex(Axes.Z))
+tLen = data.dimension(data.dimensionIndex(Axes.TIME))
+cLen = data.dimension(data.dimensionIndex(Axes.CHANNEL))
+
+# crop a channel
+c0=ops.transform().crop(img, Intervals.createMinMax(0, 0, 3, 5, 0, xLen-1, yLen-1, 3, 5, tLen-1))
+converted = ops.convert().float32(c0)
+print(converted.numDimensions())
+
+minvalue = 6  
+op1 = sub(converted, 750)
+op = IF(LT(op1, 0), THEN(minvalue), ELSE(op1))
+
+res = computeIntoFloat(op)
+res = ops.convert().uint16(res)
+
+#print(ops.op("stat.min",res))
+#IL.show(converted)
+aaa=IL.wrap(res, "ImgMath view")
+aaa.show() 
+print(type(aaa))
+one = Views.interval(res, Intervals.createMinMax(0, 0, 3, xLen-1, yLen-1, 3))
+
+resa = float(ops.run("stats.max", one).toString())
+print(resa)
+IL.show(one)
\ No newline at end of file

ImageJ/script_mosaic_telemos.py

@@ -41,12 +41,85 @@ from utils.ImagesLoaders import load_image_sequence
 # Import to do operation with ImageJ2 and ImgMath
 from net.imglib2.util import Intervals
 from net.imglib2.img.display.imagej import ImageJFunctions as IL
-from net.imglib2.algorithm.math.ImgMath import computeIntoFloat, sub, div, GT, LT, IF, THEN, ELSE, minimum
+from net.imglib2.algorithm.math.ImgMath import computeIntoFloats, sub, div, GT, LT, IF, THEN, ELSE, minimum
 from net.imglib2.type.numeric.real import FloatType
 from net.imglib2.view import Views 
 from net.imagej import Dataset
 from net.imagej.axis import Axes
 
+def getFirstMetadata(metadata):
+    """
+    Return the first keys starting with Metadata
+    """
+    ml = list(metadata)
+    cpt = 0
+    for l in ml:
+        cpt += 1
+        if l.startswith('Metadata'):
+            break
+    
+    return cpt
+    
+def getMMversion(metadata):
+    """
+    Get MM version from metadata 
+    """
+
+    if int(metadata['Summary']["MicroManagerVersion"].split('.')[0]) >= 2:
+        version = 2
+    else:
+        version = 1
+
+    return version
+
+
+def getRoi(metadata, mm_version=2):
+    """
+    Return the ROI rectangle from metadata file
+    """
+
+    if mm_version == 2:
+        # Use the first image to get the ROI
+        i = getFirstMetadata(metadata)
+        roi = metadata[metadata.keys()[i]]['ROI']
+    else:
+        roi = metadata['Summary']['ROI']
+
+    roi = [int(r) for r in roi.split('-')]
+    return roi
+
+
+def getScale(metadata):
+    """
+    Return the scale from metadata of MM
+    """
+
+    try:
+        # MM2 version
+        i = getFirstMetadata(metadata)
+        scale=metadata[metadata.keys()[i]]['PixelSizeUm']
+    except:
+        # MM1 version
+        scale = metadata['Summary']['PixelSize_um']
+        
+    return scale
+
+def getWidthHeight(metadata):
+    """
+    Get the width and height from metadata
+    """
+    try:
+        # MM2 version
+        i = getFirstMetadata(metadata)
+        height = metadata[metadata.keys()[i]]['Height']
+        width = metadata[metadata.keys()[i]]['Width']
+    except:
+        # MM1 version
+        height = metadata['Summary']['Height']
+        width = metadata['Summary']['Width']
+        
+    return width, height
+    
 def load_image_tiles_stack(images_dir):
     """
     Fonction pour charge le stack d'images depuis
@@ -117,12 +190,20 @@ def load_image_tiles_stack(images_dir):
         file_pattern = selected_roi+'-Pos_{rows:003d}_{cols:003d}'
         mmetapath = os.path.join(base_dir,file_pattern.format(rows=0, cols=0),'metadata.txt')
 
+    # Chargement des métadonnées MM (juste le resumée de la premiere tuile suifit)
+    metadata = json.load(open(mmetapath))
+    mmversion = getMMversion(metadata)
+    
+    # Test de version sur MicroManager
+    if mmversion >= 2:
+        # Add MM2.0 file structure
+        file_pattern = os.path.join(file_pattern, 'img_channel{filter}_position{tile-1:003d}_time000000000_z{z:003d}.tif')
+        hstack, good_names = load_image_sequence(base_dir, file_pattern, return_filenames=True)
+    else:
     	# Add MM1.4 file structure
     	file_pattern = os.path.join(file_pattern, 'img_000000000_{filter}_{z:003d}.tif')
     	hstack, good_names = load_image_sequence(base_dir, file_pattern, return_filenames=True)
 
-    # Chargement des métadonnées MM (juste le resumée de la premiere tuile suifit)
-    metadata = json.load(open(mmetapath))['Summary']
     
     hstack.setTitle("Data")
 
@@ -130,13 +211,18 @@ def load_image_tiles_stack(images_dir):
 
 def load_images(images_dir, roi=None, virtual=False):
     d = os.path.abspath(images_dir)
-    # MM1.4 pattern for white
+    # MM1.4 pattern for white has a pos0
+    if os.path.isdir(os.path.join(d,'Pos0')):
         white_pattern = os.path.join('Pos0','img_000000000_{filter}_{z:003d}.tif')
-    imp = load_image_sequence(d, white_pattern)
+    else:
+        # Use MM2 pattern img_channel000_position000_time000000000_z000.tif
+        white_pattern = os.path.join('Default','img_channel{filter}_position000_time000000000_z{z:003d}.tif')
         
+    imp = load_image_sequence(d, white_pattern)
     if roi is not None and len(roi) == 4 and imp.width > roi[2] and imp.height > roi[3]:
+        IJ.log('Roi %s on image size (%i, %i)' % (str(roi[2])+','+str(roi[3]), imp.width, imp.height))
         imp.setRoi(*roi)
-       #IJ.run(img, "Crop", "")
+        #IJ.run(imp, "Crop", "")
 
     return imp
 
@@ -208,7 +294,7 @@ def process_white(implus_white, implus_darkofwhite, zofwhite, roi=None):
 
     # Substract white and dark
     opsus = sub(img_whiteC, img_darkC)
-    whitemindark = computeIntoFloat(opsus)
+    whitemindark = computeIntoFloats(opsus)
     maxwhite = float(ops.run("stats.max", whitemindark).toString())
 
     # Compute the minimum value of the ratio
@@ -220,7 +306,7 @@ def process_white(implus_white, implus_darkofwhite, zofwhite, roi=None):
             THEN(div(opsus,maxwhite)), 
             ELSE(ratiomin))
 
-    result = computeIntoFloat(op)
+    result = computeIntoFloats(op)
     # Convert back to ImagePlus
     img_white_corr = IL.wrap(result, 
                              "white corrected from dark (z: %i)" % zofwhite)
@@ -253,7 +339,7 @@ def process_tiles(implus_tiles_stack, channel, z, implus_white_corr, implus_dark
         The stack containing the dark images for each channels
 
     - metadata: dict,
-        The metadata of MM 1.4 (only summary)
+        The metadata of MM 1.4 or MM2
 
     return
     ------
@@ -266,11 +352,31 @@ def process_tiles(implus_tiles_stack, channel, z, implus_white_corr, implus_dark
     # Extract the dark for the color position of the dataset:
     # TODO: Fix selecting channel when dark and white does not
     # have the same amount of filter
-    dark = IL.wrap(implus_dark_stack)
+    dark = convertService.convert(implus_dark_stack, Dataset)
+    if dark.numDimensions() > 2:
         dark = Views.hyperSlice(dark, 2, channel-1)
+    else:
+        # Need to apply crop on dark
+        roi = getRoi(metadata)
+        if roi is not None and len(roi) == 4 and implus_dark_stack.width > roi[2] and implus_dark_stack.height > roi[3]:
+            xmin, ymin, xmax, ymax = roi
+            xmax = xmin + xmax 
+            ymax = ymin + ymax 
+        else:
+            xmin = ymin = 0
+            xmax = implus_dark_stack.width
+            ymax = implus_dark_stack.height
+            
+        
+        dark = ops.transform().crop(dark, 
+                              Intervals.createMinMax(xmin, ymin,
+                                                           xmax-1, ymax-1))
+                                                           
+        IJ.log('dark is not an hyperslice, only one dark, keep this one')
+
     
     # Prepare the white
-    white = IL.wrap(implus_white_corr)
+    white = convertService.convert(implus_white_corr, Dataset)
 
     # Extract the zposition and channel of the data-to-procced
     data = convertService.convert(implus_tiles_stack, Dataset)
@@ -297,22 +403,24 @@ def process_tiles(implus_tiles_stack, channel, z, implus_white_corr, implus_dark
     if inter is not None:
         data = ops.transform().crop(data, inter)
 
+
     # Operation
     op = [None] * Ntiles
     cpt = 0
 
     for tile in xrange(Ntiles):
         IJ.showProgress(tile, Ntiles)
+        IJ.log('Compute %i / %i' % (tile, Ntiles))
         data_views = Views.hyperSlice(data, 2, tile)
         opsus = sub(data_views, dark)
-        op[cpt] = IF(GT(opsus, 0), THEN(div(opsus, white)), ELSE(0))
+        op[cpt] = IF(GT(opsus, 0.0), THEN(div(opsus, white)), ELSE(0.0))
         cpt += 1
 
     result = datasetService.create(Views.stack([o.view(FloatType()) for o in op]))
     result = ops.convert().uint16(result)
 
     # Create the ImagePlus from the result converted back to uint16
-    name = metadata['ChNames'][channel-1]
+    name = metadata['Summary']['ChNames'][channel-1]
     imp_title = "Processed Tiles (channel: %i-%s, z: %i)" % (channel, name, z)
     img_tiles_corr = IL.wrap(result, imp_title)
     # Convert axis to Time and se the resolution
@@ -336,7 +444,7 @@ def run_processing(channel, zimg, zwhite, run_basic=True, optimizeXY=True, savep
     if not os.path.exists(savepath):
        os.makedirs(savepath)
 
-    chname = metadata['ChNames'][channel-1]
+    chname = metadata['Summary']['ChNames'][channel-1]
     # A folder for the processed tiles
     procfolder = savepath+'/DW_processed_tiles_%s_z%i_whitez%i' % (chname, zimg, zwhite)
     if not os.path.exists(procfolder):
@@ -402,8 +510,9 @@ def run_processing(channel, zimg, zwhite, run_basic=True, optimizeXY=True, savep
     proctiles.close()
 
 def run_gui(default_values=None):
-    channels = metadata['ChNames']
-    dataZmax = metadata['Slices']
+
+    channels = metadata['Summary']['ChNames']
+    dataZmax = metadata['Summary']['Slices']
     whiteZmax = white.getNSlices()
 
     if default_values is None:
@@ -476,30 +585,32 @@ def run_gui(default_values=None):
        white.close()
 
 
-def extract_tile_positions(metadata, selected_roi, tilesondisk, tilename='tile_'):
+def getMMtiles(metadata, selected_roi, tilesondisk):
     """
-    Extraction des positions des tuiles depuis les métadonnées
-    pour les rendres compatible avec le plugin de stitching.
+    Extract tiles position from micromanager metadata
+
+    tilesondisk, list of tiles name found on disk
     """
     
-    IJ.log("Extract tiles coordinates from metadata")
-    # Extract pixel size
-    pix2um = metadata['PixelSize_um']
-    # The roi width
-    roiwidth = metadata['Width']
+    mmversion = getMMversion(metadata)
 
+    mmsum = metadata['Summary']
+    xt = None
+    yt = None
+    
+    if mmversion == 1:
         # Positions
-    xt = [None] * len(metadata['InitialPositionList'])
-    yt = [None] * len(metadata['InitialPositionList'])
+        xt = [None] * len(mmsum['InitialPositionList'])
+        yt = [None] * len(mmsum['InitialPositionList'])
         test_nameSE = 'roi'+selected_roi
         test_nameGC = '%s-Pos_' % selected_roi
 
-    uniquecols = set(a['GridColumnIndex'] for a in metadata['InitialPositionList'])
-    uniquerows = set(a['GridRowIndex'] for a in metadata['InitialPositionList'])
+        uniquecols = set(a['GridColumnIndex'] for a in mmsum['InitialPositionList'])
+        uniquerows = set(a['GridRowIndex'] for a in mmsum['InitialPositionList'])
         Nrow = max(uniquerows) + 1
         Ncol = max(uniquecols) + 1
 
-    for i, pl in enumerate(metadata['InitialPositionList']):
+        for i, pl in enumerate(mmsum['InitialPositionList']):
             if pl['Label'] in tilesondisk:
                 # print(pl['Label'])
                 if test_nameSE in pl['Label']:
@@ -519,11 +630,77 @@ def extract_tile_positions(metadata, selected_roi, tilesondisk, tilename='tile_'
             else:
                 IJ.log("[WARNING]: File %s listed on metadata does not exist on disk" % pl['Label'])
 
+    else:
+        # for MM2
+        xt = [None] * len(mmsum['StagePositions'])
+        yt = [None] * len(mmsum['StagePositions'])
+        test_nameSE = 'roi'+selected_roi
+        test_nameGC = '%s-Pos_' % selected_roi
+
+        uniquecols = set(a['GridCol'] for a in mmsum['StagePositions'])
+        uniquerows = set(a['GridRow'] for a in mmsum['StagePositions'])
+        Nrow = max(uniquerows) + 1
+        Ncol = max(uniquecols) + 1
+
+        for i, pl in enumerate(mmsum['StagePositions']):
+            if pl['Label'] in tilesondisk:
+                # print(pl['Label'])
+                if test_nameSE in pl['Label']:
+                    ind = i
+
+                if test_nameGC in pl['Label']:
+                    # Compute the globel indice ind = col*(Nrow)+row
+                    # as the order is not the same as the one loaded by bio-format reader
+                    c = pl['GridCol']
+                    r = pl['GridRow']
+                    ind = c * Nrow + r
+                    # print(ind, pl['Label'])
+
+                xt[ind] = pl['DevicePositions'][0]['Position_um'][0]
+                yt[ind] = pl['DevicePositions'][0]['Position_um'][1]
+
+            else:
+                IJ.log("[WARNING]: File %s listed on metadata does not exist on disk" % pl['Label'])
+
+    return xt, yt
+    
+def extract_tile_positions(metadata, selected_roi, tilesondisk, tilename='tile_'):
+    """
+    Extraction des positions des tuiles depuis les métadonnées
+    pour les rendres compatible avec le plugin de stitching.
+    """
+    
+    IJ.log("Extract tiles coordinates from metadata")
+    # Extract pixel size
+    pix2um = getScale(metadata)
+    # The roi width
+    roiwidth, _ = getWidthHeight(metadata)
+
+    # Is the microscope use Andor camera ?
+    i = getFirstMetadata(metadata)
+    if metadata[metadata.keys()[i]]['Camera'] == 'Andor':
+        ACAMAS = True
+    else:
+        ACAMAS = False
+        
+    xt, yt = getMMtiles(metadata, selected_roi, tilesondisk)
+    
+    if ACAMAS:
+        xtmp = yt
+        ytmp = xt
+        xt = xtmp
+        yt = ytmp
+        #xt = [-x for x in xtmp]
+        #yt = [-y for y in ytmp]        
+        
     # Remove None values in xt and yt (when the tile is not on disk)
     igood = [i for i,v in enumerate(xt) if v is not None]
     xt = [xt[i] for i in igood]
     yt = [yt[i] for i in igood]
 
+    if ACAMAS:
+        xt = [-x for x in xt]
+        
     xt = map(lambda x: x/pix2um, xt)
     yt = map(lambda x: x/pix2um, yt)
 
@@ -533,6 +710,7 @@ def extract_tile_positions(metadata, selected_roi, tilesondisk, tilename='tile_'
     xt = [x-xt[0] for x in xt]
     yt = [y-yt[0] for y in yt]
 
+        
     outstr = []
     for i in range(len(yt)):
         outstr += ['%s%04i.tif; ; (%0.2f, %0.2f)'%(tilename, i+1, yt[i], -xt[i])]
@@ -564,7 +742,11 @@ if __name__ in ['__builtin__','__main__']:
     hstack, metadata, rootdir, selectedroi, file_names = load_image_tiles_stack(fdata)
     hstack.show()
 
-    roi = metadata['ROI']
+	# GET the version of MM
+    mmversion = getMMversion(metadata)
+
+    # LOAD ROI info
+    roi = getRoi(metadata, mmversion)
     IJ.log('ROI: %s' % str(roi))
     
     # Make a list of tiles names

ImageJ/utils/ImagesLoaders.py

-# -*- coding: utf-8 -*-
-"""
-Contain function or class to make import of image sequences easier in ImageJ
-"""
-import os
-from path_tools import find_groups_values
-from string import Formatter
-from itertools import product
-
-from loci.formats import ChannelSeparator
-from ij import ImagePlus, VirtualStack, IJ
-
-
-def load_image_sequence(rootdir, path_pattern, order=None, channel_group=None,
-                        z_group=None, return_filenames=False):
-    """
-    Load an image sequence (several tiff files in several sub-folders) as an ImageJ VirtualStack.
-    It uses the python string format to define the pattern of file names
-    If order is specified as a list of pattern keys, image will be loaded in that order,
-    otherwise the order are the one in the pattern file names.
-
-    Parameters:
-    -----------
-
-    rootdir: string,
-        The fixed part of the path to image folder and subfolders
-    path_pattern: string,
-        The moving part of the folder/file_names.ext where you can use
-        python string format convention.
-    order: list of string or None,
-        The order of pattern keys in which images should be loaded
-    channel_group: string or None,
-        The key which contains the color-band dimension of your stack
-    z_group: string or None,
-        The key of the z dimension of your stack
-    return_filenames: boolean,
-        If true return the list of file names
-
-    Examples:
-    ---------
-
-    TODO
-    """
-
-    # Check that path_pattern doest not start with a /
-    if path_pattern.startswith(os.path.sep):
-        path_pattern = path_pattern[len(os.path.sep):]
-
-    # Find group values
-    gp_values = find_groups_values(os.path.join(rootdir, path_pattern))
-
-    if order is None:
-        keys = list(i[1] for i in Formatter().parse(path_pattern))
-        keys_order = []
-        for k in keys:
-            if k is not None and k not in keys_order:
-                keys_order += [k]
-
-    else:
-        keys_order = order
-
-    # Set some smath things (if key is c or channel or filter set as channel_group, or if key is z set as z_group)
-    if z_group is None:
-        for t in ['z', 'Z']:
-            if t in keys_order:
-                z_group = t
-                break
-
-    if channel_group is None:
-        for t in ['c', 'channel', 'channels', 'filter', 'filters']:
-            if t in keys_order:
-                channel_group = t
-                break
-
-            if t.upper() in keys_order:
-                channel_group = t.upper()
-                break
-
-            if t.capitalize() in keys_order:
-                channel_group = t.capitalize()
-                break
-
-    # Z should be before channel group but at the end (if given for hyperstack dimension)
-    if z_group is not None:
-        keys_order.remove(z_group)
-        keys_order += [z_group]
-
-    # If a channel group is given (need to be at the end of keys_order)
-    if channel_group is not None:
-        keys_order.remove(channel_group)
-        keys_order += [channel_group]
-
-    # Create the combination of all possible groups values
-    allpattern_values = list(product(*(gp_values[k] for k in keys_order)))
-
-    # Generate the list of file names to load
-    allnames = tuple(path_pattern.format(**dict(zip(keys_order, p))) for p in allpattern_values)
-
-    # We need to load one image to get it's size
-    fr = ChannelSeparator()
-    fr.setGroupFiles(False)
-    fr.setId(os.path.join(rootdir, allnames[0]))
-    metadata = fr.getGlobalMetadata()
-    width, height = fr.getSizeX(), fr.getSizeY()
-    # Create a virtualStack
-    vstack = VirtualStack(width, height, None, rootdir)
-    for imgname in allnames:
-        vstack.addSlice(imgname)
-
-    # Create the ImagePlus to display the stack
-    imp = ImagePlus("%s" % rootdir, vstack)
-
-    # Set the dimension of hyperstack if keys for color channels
-    # and z channels are given
-    if channel_group is not None:
-        nchannel = len(gp_values[channel_group])
-    else:
-        nchannel = 1
-
-    if z_group is not None:
-        nslice = len(gp_values[z_group])
-    else:
-        nslice = 1
-
-    # Compute the time from the product of not used channel
-    not_used_keys = []
-    for k in keys_order:
-        if k != channel_group and k != z_group:
-            not_used_keys += [k]
-
-    ntimes = 1
-    for k in not_used_keys:
-        ntimes *= len(gp_values[k])
-
-    imp.setDimensions(nchannel, nslice, ntimes)
-    imp.setOpenAsHyperStack(True)
-
-    # Try to set the physical size from Metadata
-    try:
-        IJ.run(imp,
-               "Properties...",
-               "pixel_width={resox} pixel_height={resoy} pixel_unit='um'".format(resox=1/metadata['XResolution'],
-                                                                                 resoy=1/metadata['YResolution']));
-    except:
-        pass
-
-    # Set auto contrast
-    IJ.run(imp, "Enhance Contrast", "saturated=0.35");
-
-    if return_filenames:
-        return imp, allnames
-    else:
-        return imp