Added Resolution info to exported TIFF.

CHANGELOG.md

@@ -21,6 +21,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 - `tiff_auto3dseg`
     + Now using `ggc.check_threads()`.
 - Fixed `--compressed` option label (now compatible with ImageJ).
+- `tools.plot.save_tif`
+    + Added support to retain voxel resolution in TIFF metadata (ImageJ compatible).
+- `nd2_to_tiff`
+    + Now saves Resolution (XYZ) metadata when exporting.
 
 
 

bin/czi_to_tiff

@@ -46,6 +46,7 @@ import czifile
 import os
 import sys
 import tifffile
+from tqdm import tqdm
 import warnings
 
 from pygpseq.tools import image as imt
@@ -142,6 +143,12 @@ channel_path = "Metadata/DisplaySetting/Channels/Channel/DyeName"
 channel_names = [x.text for x in images.metadata.findall(channel_path)]
 channel_names = [x.replace(" ", "").lower() for x in channel_names]
 
+res_path = "Metadata/Scaling/Items/Distance"
+resolution = [x for x in images.metadata.findall(res_path)
+    if x.attrib['Id'] in ["X", "Y", "Z"]]
+resolution = dict([(x.attrib['Id'], float(x[0].text)) for x in resolution])
+resolutionXY = (1e-6/resolution["X"], 1e-6/resolution["Y"])
+
 # Remove previous axes
 if 0 != Si:
     for i in range(Si):
@@ -156,7 +163,7 @@ if axes.endswith("0"):
     pixels = pixels.reshape(pixels.shape[:-1])
     axes = axes[:-1]
 
-for si in range(images.shape[Si]):
+for si in tqdm(range(images.shape[Si])):
     for ci in range(images.shape[Ci]):
         stack = pixels[si, ci]
 
@@ -182,6 +189,8 @@ for si in range(images.shape[Si]):
 
         plot.save_tif(os.path.join(args.outdir, outpath),
             stack, imt.get_dtype(stack.max()), args.doCompress,
-            bundled_axes = stack_axes)
+            bundled_axes = stack_axes, resolution = resolutionXY,
+            inMicrons = True, ResolutionZ = resolution["Z"]*1e6,
+            forImageJ = True)
 
 ################################################################################

bin/nd2_to_tiff

@@ -43,9 +43,12 @@
 
 import argparse
 from nd2reader import ND2Reader
+from nd2reader.parser import Parser as ND2Parser
+import numpy as np
 import os
 import sys
 import tifffile
+from tqdm import tqdm
 
 from pygpseq.tools import image as imt
 from pygpseq.tools import plot
@@ -78,6 +81,9 @@ parser.add_argument('-m', '--mode', type = str,
     choices = ("DOTTER", "GPSeq"), metavar = 'mode',
     help = """Output filename notation. Default: GPSeq.""",
     default = "GPSeq")
+parser.add_argument('-Z', '--deltaZ', type = float, metavar = 'dZ',
+    help = """If provided (in um), the script does not check delta Z consistency
+    and instead uses the provided one.""", default = None)
 
 # Add flags
 parser.add_argument('--compressed',
@@ -86,7 +92,7 @@ parser.add_argument('--compressed',
     help = 'Force compressed TIFF as output.')
 
 # Version flag
-version = "2.0.0"
+version = "2.1.0"
 parser.add_argument('--version', action = 'version',
     version = '%s %s' % (sys.argv[0], version,))
 
@@ -109,6 +115,29 @@ if not os.path.isdir(args.outdir): os.mkdir(args.outdir)
 
 # FUNCTIONS ====================================================================
 
+def getOutpath(args, metadata, cid, fid):
+    outpath = None
+    # Identify ouytput file name notation
+    if "GPSeq" == args.mode:
+        outpath = "%s.channel%03d.series%03d.tif" % (
+            metadata['channels'][cid].lower(), cid + 1, fid + 1)
+    elif "DOTTER" == args.mode:
+        outpath = "%s_%03d.tif" % (
+            metadata['channels'][cid].lower(), fid + 1)
+    return(outpath)
+
+def export_channel(args, ch, outpath, metadata,
+    bundled_axes, resolutionZ):
+    resolutionXY = (1/metadata['pixel_microns'],
+        1/metadata['pixel_microns'])
+
+    plot.save_tif(os.path.join(args.outdir, outpath),
+        ch, imt.get_dtype(ch.max()), args.doCompress,
+        bundled_axes = bundled_axes.upper(),
+        resolution = resolutionXY,
+        inMicrons = True, forImageJ = True,
+        ResolutionZ = resolutionZ)
+
 def export_fov(fov, metadata, fid, bundled_axes):
     '''Export a field of view after bundling the axes to separate channel TIFF.
 
@@ -118,68 +147,63 @@ def export_fov(fov, metadata, fid, bundled_axes):
         fid (int): field of view 0-based index.
     '''
 
-    if "c" in bundled_axes:
-        print(" Field #%d: found %d channels %s." % (
-            fid + 1, fov.shape[3], str(metadata['channels'])))
 
+    # Get Z resolution
+    if not type(None) == type(args.deltaZ):
+        resolutionZ = args.deltaZ
+    else:
+        with open(args.input, "rb") as fh:
+            p = ND2Parser(fh)
+
+            Zdata = np.array(p._raw_metadata.z_data)
+            Zlevels = np.array(p.metadata['z_levels']).astype('int')
+            Zlevels = Zlevels + len(Zlevels) * fid
+            Zdata = Zdata[Zlevels]
+
+            resolutionZ = set(np.round(np.diff(Zdata), 3))
+
+            assert_msg = "Z resolution is not constant: %s" % (str(resolutionZ))
+            assert 1 == len(resolutionZ), assert_msg
+
+            resolutionZ = list(resolutionZ)[0]
+
+    if "c" in bundled_axes:
         # Iterate over channels
         for cid in range(fov.shape[3]):
+            outpath = getOutpath(args, metadata, cid, fid)
+            if type(None) == type(outpath): return()
             ch = fov[:, :, :, cid]
-
-            # Identify ouytput file name notation
-            if "GPSeq" == args.mode:
-                outpath = "%s.channel%03d.series%03d.tif" % (
-                    metadata['channels'][cid].lower(), cid + 1, fid + 1)
-            elif "DOTTER" == args.mode:
-                outpath = "%s_%03d.tif" % (
-                    metadata['channels'][cid].lower(), fid + 1)
-            else: return()
-
-            print("  Exporting %s..." % outpath)
-            plot.save_tif(os.path.join(args.outdir, outpath),
-                ch, imt.get_dtype(ch.max()), args.doCompress,
-                bundled_axes = bundled_axes.upper()[:-1])
+            export_channel(args, ch, outpath, metadata,
+                bundled_axes, resolutionZ)
     else:
-        print(" Field #%d: found 1 channel (%s)." % (
-            fid + 1, metadata['channels'][0]))
-
-        # Identify ouytput file name notation
-        if "GPSeq" == args.mode:
-            outpath = "%s.channel%03d.series%03d.tif" % (
-                metadata['channels'][0].lower(), 1, fid + 1)
-        elif "DOTTER" == args.mode:
-            outpath = "%s_%03d.tif" % (
-                metadata['channels'][0].lower(), fid + 1)
-        else: return()
-
-        print("  Exporting %s..." % outpath)
-        plot.save_tif(os.path.join(args.outdir, outpath),
-            fov, imt.get_dtype(fov.max()), args.doCompress,
-            bundled_axes = bundled_axes.upper())
+        outpath = getOutpath(args, metadata, 1, fid)
+        if type(None) == type(outpath): return()
+        export_channel(args, fov, outpath, metadata,
+            bundled_axes, resolutionZ)
 
 # RUN ==========================================================================
 
+if not type(None) == type(args.deltaZ):
+    print("Enforcing a deltaZ of %.3f um." % args.deltaZ)
+
 # Create buffer pointer to nd2 image
 images = ND2Reader(args.input)
 
 # Check if multiple fields of view are present
 axes_for_bundling = "zyxc" if "c" in images.axes else "zyx"
 if 'v' not in images.axes:
-    print("Single field of view found.")
-
+    print("Found 1 field of view, with %d channels." % (images.sizes['c']))
     # Export single field of view
     images.bundle_axes = axes_for_bundling
     export_fov(images[0], images.metadata, 0, axes_for_bundling)
 else:
-    print("Multiple fields of view found (%d)." % images.sizes['v'])
-
+    print("Found %d fields of view, with %d channels each." % (
+        images.sizes['v'], images.sizes['c']))
     # Export multiple fields of view
     images.iter_axes = 'v'
     images.bundle_axes = axes_for_bundling
-    for fid in range(images.sizes['v']):
+    for fid in tqdm(range(images.sizes['v'])):
         fov = images[fid]
         export_fov(fov, images.metadata, fid, axes_for_bundling)
 
-print("DONE")
-
 ################################################################################

pygpseq/tools/plot.py

@@ -869,7 +869,9 @@ def profile_density_scatterplot(pp, field_label, nbins, xyrange,
     if new_figure:
         return(fig)
 
-def save_tif(path, img, dtype, compressed, bundled_axes = "CZYX"):
+def save_tif(path, img, dtype, compressed, bundled_axes = "CZYX",
+    inMicrons = False, ResolutionZ = None, forImageJ = False,
+    **kwargs):
     # Add channel axis for ImageJ compatibility
     if not "C" in bundled_axes:
         bundled_axes = "C" + bundled_axes
@@ -883,14 +885,20 @@ def save_tif(path, img, dtype, compressed, bundled_axes = "CZYX"):
     assert_msg = "shape mismatch between bundled axes and image."
     assert len(bundled_axes) == len(img.shape), assert_msg
 
+    metadata = {'axes' : bundled_axes}
+    if inMicrons: metadata['unit'] = "um"
+    if not type(None) == ResolutionZ: metadata['spacing'] = ResolutionZ
+
     if compressed:
         tifffile.imsave(path, img.astype(dtype),
             shape = img.shape, compress = 9,
-            dtype = dtype, imagej = False, metadata = {'axes' : bundled_axes})
+            dtype = dtype, imagej = forImageJ,
+            metadata = metadata, **kwargs)
     else:
         tifffile.imsave(path, img.astype(dtype),
             shape = img.shape, compress = 0,
-            dtype = dtype, imagej = False, metadata = {'axes' : bundled_axes})
+            dtype = dtype, imagej = forImageJ,
+            metadata = metadata, **kwargs)
 
 def set_font_size(size):
     """Set plot font size."""

setup.py

@@ -42,12 +42,12 @@
 		'jinja2==2.10',
 		'joblib==0.11',
 		'matplotlib==2.2.2',
-		'nd2reader==3.0.8',
+		'nd2reader>=3.1.0',
 		'numpy>=1.14.2',
 		'pandas>=0.22.0',
 		'scipy>=1.0.0',
 		'scikit-image==0.14.0',
-		'tifffile==0.14.0',
+		'tifffile>=0.15.1',
 		'tqdm>=4.23.4',
 		'weasyprint>=0.42.2'],
 	scripts=[os.path.join(bindir, fp) for fp in os.listdir(bindir)],