updates for compatibility with input zstacks

Blaze microscope now outputs data with one tif file per zstack.
This adds compatibility for it.

TODO: test full end-to-end

config/config.yml

@@ -1,6 +1,5 @@
 datasets: 'config/datasets.tsv'
 
-
 root: 'bids' # can use a s3:// or gcs:// prefix to write output to cloud storage
 work: 'work' 
 
@@ -13,6 +12,7 @@ cores_per_rule: 32
 #import wildcards:  tilex, tiley, channel, zslice (and prefix - unused)
 import_blaze:
   raw_tif_pattern: "{prefix}_Blaze[{tilex} x {tiley}]_C{channel}_xyz-Table Z{zslice}.ome.tif"
+  raw_tif_pattern_zstack:  "{prefix}_Blaze[{tilex} x {tiley}]_C{channel}.ome.tif"
   intensity_rescaling: 0.5 #raw images seem to be at the upper end of uint16 (over-saturated) -- causes wrapping issues when adjusting with flatfield correction etc. this rescales the raw data as it imports it..
 
 import_prestitched:

workflow/lib/dask_image.py

@@ -0,0 +1,45 @@
+from __future__ import annotations
+
+import os
+from glob import glob
+
+try:
+    import tifffile
+except (AttributeError, ImportError):
+    pass
+
+from dask.array.core import Array
+from dask.base import tokenize
+
+
+def add_leading_dimension(x):
+    return x[None, ...]
+
+def imread(fn, page):
+    return tifffile.imread(fn,key=page)
+
+
+def imread_pages(filename, preprocess=None):
+
+
+    tif = tifffile.TiffFile(filename)
+    pages = [i for i in range(len(tif.pages))]
+    name = "imread-%s" % tokenize(filename, os.path.getmtime(filename))
+
+    sample = tif.pages[0].asarray()
+
+    if preprocess:
+        sample = preprocess(sample)
+
+    keys = [(name, i) + (0,) * len(sample.shape) for i in pages]
+    if preprocess:
+        values = [
+            (add_leading_dimension, (preprocess, (imread, filename,i))) for i in pages
+        ]
+    else:
+        values = [(add_leading_dimension, (imread, filename,i)) for i in pages]
+    dsk = dict(zip(keys, values))
+
+    chunks = ((1,) * len(pages),) + tuple((d,) for d in sample.shape)
+
+    return Array(dsk, name, chunks, sample.dtype)

workflow/rules/import.smk

@@ -79,11 +79,6 @@ rule blaze_to_metadata_gcs:
 rule blaze_to_metadata:
     input:
         ome_dir=get_input_dataset,
-    params:
-        in_tif_pattern=lambda wildcards, input: os.path.join(
-            input.ome_dir,
-            config["import_blaze"]["raw_tif_pattern"],
-        ),
     output:
         metadata_json=temp(
             bids(
@@ -197,10 +192,6 @@ rule tif_to_zarr:
         ome_dir=get_input_dataset,
         metadata_json=rules.copy_blaze_metadata.output.metadata_json,
     params:
-        in_tif_pattern=lambda wildcards, input: os.path.join(
-            input.ome_dir,
-            config["import_blaze"]["raw_tif_pattern"],
-        ),
         intensity_rescaling=config["import_blaze"]["intensity_rescaling"],
     output:
         zarr=temp(

workflow/scripts/blaze_to_metadata.py

@@ -4,25 +4,48 @@
 import re
 from itertools import product
 from snakemake.io import glob_wildcards
+from glob import glob
+from pathlib import Path
+import os
 
-in_tif_pattern = snakemake.params.in_tif_pattern
+tif_files = glob(f"{snakemake.input.ome_dir}/*.tif")
+
+#check first tif file to see if it is zstack or not:
+if 'xyz-Table Z' in Path(tif_files[0]).name:
+    is_zstack=False
+else:
+    is_zstack=True
+
+if is_zstack:
+    in_tif_pattern = os.path.join(snakemake.input.ome_dir,snakemake.config["import_blaze"]["raw_tif_pattern_zstack"])
+else:
+    in_tif_pattern = os.path.join(snakemake.input.ome_dir,snakemake.config["import_blaze"]["raw_tif_pattern"])
 
 #add a wildcard constraint to ensure no
 #subfolders get parsed (ie don't match anything with / in it):
 prefix_constraint=r'[^/]+' 
 in_tif_pattern_constrained = in_tif_pattern.replace('{prefix}',f'{{prefix,{prefix_constraint}}}')
 
+
 #parse the filenames to get number of channels, tiles etc..
-prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern_constrained)
+if is_zstack:
+    prefix, tilex, tiley, channel = glob_wildcards(in_tif_pattern_constrained,tif_files)
+else:
+    prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern_constrained,tif_files)
 
 tiles_x = sorted(list(set(tilex)))
 tiles_y = sorted(list(set(tiley)))
 channels = sorted(list(set(channel)))
-zslices = sorted(list(set(zslice)))
 prefixes = sorted(list(set(prefix)))
 
-#read in series metadata from first file
-in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+if not is_zstack:
+    zslices = sorted(list(set(zslice)))
+
+    #read in series metadata from first file
+    in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+
+else:
+    in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0])
 
 raw_tif = tifffile.TiffFile(in_tif,mode='r')
 
@@ -37,12 +60,18 @@
 custom_metadata = ome_dict['OME']['Image']['ca:CustomAttributes']
 
 
-
 #read tile configuration from the microscope metadata
 if axes == 'CZYX':
-    tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+), (?P<z>\S+)\)"
+    if is_zstack:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+)\)"
+    else:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+), (?P<z>\S+)\)"
 elif axes == 'ZYX': 
-    tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+), (?P<z>\S+)\)"
+    if is_zstack:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+)\)"
+    else:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+), (?P<z>\S+)\)"
+
 
 tile_pattern = re.compile(tile_config_pattern)
 
@@ -58,23 +87,29 @@
     chunks.append(chunk)
 
 for line in  custom_metadata['TileConfiguration']['@TileConfiguration'].split('  ')[1:]:
-
     d = re.search(tile_pattern,line).groupdict()
     chunk = map_tiles_to_chunk[d['tilex']+d['tiley']] # want the key to have chunk instad of tilex,tiley, so map to that first
-
-    #key is:  tile-{chunk}_chan-{channel}_z-{zslice} 
-    key = f"tile-{chunk}_chan-{d['channel']}_z-{d['zslice']}" 
+
+    if is_zstack:
+        key = f"tile-{chunk}_chan-{d['channel']}_z-0" 
+    else:
+        #key is:  tile-{chunk}_chan-{channel}_z-{zslice} 
+        key = f"tile-{chunk}_chan-{d['channel']}_z-{d['zslice']}" 
 
     map_x[key] = float(d['x'])
     map_y[key] = float(d['y'])
-    map_z[key] = float(d['z'])
+    if is_zstack:
+        map_z[key] = float(0)
+    else:
+        map_z[key] = float(d['z'])
 
 
 metadata={}
 metadata['tiles_x'] = tiles_x
 metadata['tiles_y'] = tiles_y
 metadata['channels'] = channels
-metadata['zslices'] = zslices
+if not is_zstack:
+    metadata['zslices'] = zslices
 metadata['prefixes'] = prefixes
 metadata['chunks'] = chunks
 metadata['axes'] = axes
@@ -88,6 +123,7 @@
 metadata['ome_full_metadata'] = ome_dict
 metadata['PixelSize'] = [ metadata['physical_size_z']/1000.0, metadata['physical_size_y']/1000.0, metadata['physical_size_x']/1000.0 ] #zyx since OME-Zarr is ZYX
 metadata['PixelSizeUnits'] = 'mm' 
+metadata['is_zstack'] = is_zstack
 
 #write metadata to json
 with open(snakemake.output.metadata_json, 'w') as fp:

workflow/scripts/blaze_to_metadata_gcs.py

@@ -5,31 +5,50 @@
 import os
 from itertools import product
 from snakemake.io import glob_wildcards
+from pathlib import Path
 import gcsfs
 from lib.cloud_io import get_fsspec
 
 dataset_uri = snakemake.params.dataset_path
-in_tif_pattern = snakemake.params.in_tif_pattern
-
 
 gcsfs_opts={'project': snakemake.params.storage_provider_settings['gcs'].get_settings().project,
                         'token': snakemake.input.creds}
 fs = gcsfs.GCSFileSystem(**gcsfs_opts)
 
 tifs = fs.glob(f"{dataset_uri}/*.tif")
 
+#check first tif file to see if it is zstack or not:
+if 'xyz-Table Z' in Path(tifs[0]).name:
+    is_zstack=False
+else:
+    is_zstack=True
+
+if is_zstack:
+    in_tif_pattern = snakemake.config["import_blaze"]["raw_tif_pattern_zstack"]
+else:
+    in_tif_pattern = snakemake.config["import_blaze"]["raw_tif_pattern"]
+
+
+
 #parse the filenames to get number of channels, tiles etc..
-prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern,files=tifs)
+if is_zstack:
+    prefix, tilex, tiley, channel = glob_wildcards(in_tif_pattern,files=tifs)
+else:
+    prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern,files=tifs)
+    zslices = sorted(list(set(zslice)))
 
 tiles_x = sorted(list(set(tilex)))
 tiles_y = sorted(list(set(tiley)))
 channels = sorted(list(set(channel)))
-zslices = sorted(list(set(zslice)))
 prefixes = sorted(list(set(prefix)))
 
+if not is_zstack:
+    zslices = sorted(list(set(zslice)))
 
-#read in series metadata from first file
-in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+    #read in series metadata from first file
+    in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+else:
+    in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0])
 
 with fs.open(f"gcs://{in_tif}", 'rb') as tif_file:
     raw_tif = tifffile.TiffFile(tif_file,mode='r')
@@ -49,9 +68,15 @@
 
 #read tile configuration from the microscope metadata
 if axes == 'CZYX':
-    tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+), (?P<z>\S+)\)"
+    if is_zstack:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+)\)"
+    else:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+),(?P<chan>\S+), (?P<z>\S+)\)"
 elif axes == 'ZYX': 
-    tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+), (?P<z>\S+)\)"
+    if is_zstack:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+)\)"
+    else:
+        tile_config_pattern=r"Blaze\[(?P<tilex>[0-9]+) x (?P<tiley>[0-9]+)\]_C(?P<channel>[0-9]+)_xyz-Table Z(?P<zslice>[0-9]+).ome.tif;;\((?P<x>\S+), (?P<y>\S+), (?P<z>\S+)\)"
 
 tile_pattern = re.compile(tile_config_pattern)
 
@@ -71,19 +96,28 @@
     d = re.search(tile_pattern,line).groupdict()
     chunk = map_tiles_to_chunk[d['tilex']+d['tiley']] # want the key to have chunk instad of tilex,tiley, so map to that first
 
-    #key is:  tile-{chunk}_chan-{channel}_z-{zslice} 
-    key = f"tile-{chunk}_chan-{d['channel']}_z-{d['zslice']}" 
+    if is_zstack:
+        key = f"tile-{chunk}_chan-{d['channel']}_z-0" 
+    else:
+        #key is:  tile-{chunk}_chan-{channel}_z-{zslice} 
+        key = f"tile-{chunk}_chan-{d['channel']}_z-{d['zslice']}" 
 
     map_x[key] = float(d['x'])
     map_y[key] = float(d['y'])
-    map_z[key] = float(d['z'])
+    if is_zstack:
+        map_z[key] = float(0)
+    else:
+        map_z[key] = float(d['z'])
 
 
+
+
 metadata={}
 metadata['tiles_x'] = tiles_x
 metadata['tiles_y'] = tiles_y
 metadata['channels'] = channels
-metadata['zslices'] = zslices
+if not is_zstack:
+    metadata['zslices'] = zslices
 metadata['prefixes'] = prefixes
 metadata['chunks'] = chunks
 metadata['axes'] = axes
@@ -97,6 +131,7 @@
 metadata['ome_full_metadata'] = ome_dict
 metadata['PixelSize'] = [ metadata['physical_size_z']/1000.0, metadata['physical_size_y']/1000.0, metadata['physical_size_x']/1000.0 ] #zyx since OME-Zarr is ZYX
 metadata['PixelSizeUnits'] = 'mm' 
+metadata['is_zstack'] = is_zstack
 
 #write metadata to json
 with open(snakemake.output.metadata_json, 'w') as fp: