Merge pull request #4 from constantinpape/main

Add initial data in ome.zarr format

.gitignore

@@ -0,0 +1,3 @@
+__pycache__/
+*.n5
+*.ome.zarr

create_data.py

@@ -0,0 +1,120 @@
+import json
+import os
+
+import z5py
+from data_conversion import convert_bdv_n5
+from pybdv.metadata import (get_size, get_resolution,
+                            write_size_and_resolution,
+                            write_affine)
+
+from mobie.xml_utils import copy_xml_as_n5_s3
+from mobie.metadata.image_dict import default_layer_setting
+
+IMAGE_DICT = './data/images.json'
+
+
+def write_metadata(in_xml, out_xml, out_path):
+    bucket_name = 'i2k-2020'
+    path_in_bucket = os.path.split(out_path)[1]
+    copy_xml_as_n5_s3(in_xml, out_xml,
+                      service_endpoint='https://s3.embl.de',
+                      bucket_name=bucket_name,
+                      path_in_bucket=path_in_bucket,
+                      authentication='Anonymous',
+                      bdv_type='bdv.zarr.s3')
+
+    with z5py.File(out_path, 'r') as f:
+        shape = f['setup0/timepoint0/s0'].shape[2:]
+
+    # check if we need to update the shape and resolution
+    exp_shape = get_size(out_xml, setup_id=0)
+    if shape != exp_shape:
+        resolution = get_resolution(out_xml, setup_id=0)
+        scale_factor = [float(esh) / sh for sh, esh in zip(shape, exp_shape)]
+        resolution = [round(res * sf, 2) for res, sf in zip(resolution, scale_factor)]
+        print("Updating shape and resolution to:")
+        print(shape)
+        print(resolution)
+
+        write_size_and_resolution(out_xml, setup_id=0,
+                                  size=shape, resolution=resolution)
+
+        # make transformation the hacky way ...
+        dz, dy, dx = resolution
+        oz, oy, ox = 0., 0., 0.
+        trafo = '{} 0.0 0.0 {} 0.0 {} 0.0 {} 0.0 0.0 {} {}'.format(dx, ox,
+                                                                   dy, oy,
+                                                                   dz, oz)
+        trafo = list(map(float, trafo.split(' ')))
+        write_affine(out_xml, setup_id=0, affine=trafo, overwrite=True)
+
+
+def add_to_image_dict(name, layer_type, xml_path):
+    settings = default_layer_setting(layer_type)
+    storage = {"remote": os.path.split(xml_path)[1]}
+    settings.update({"storage": storage})
+
+    if os.path.exists(IMAGE_DICT):
+        with open(IMAGE_DICT) as f:
+            image_dict = json.load(f)
+    else:
+        image_dict = {}
+
+    image_dict[name] = settings
+    with open(IMAGE_DICT, 'w') as f:
+        json.dump(image_dict, f, indent=2, sort_keys=True)
+
+
+def add_volume(in_path, vol_name, layer_type, start_scale=0):
+    out_path = os.path.join('data', f'{vol_name}.ome.zarr')
+
+    # convert to ome zarr
+    convert_bdv_n5(in_path=in_path,
+                   out_path=out_path,
+                   out_key='setup0/timepoint0',
+                   vol_name=vol_name,
+                   use_nested_store=False,
+                   n_threads=8,
+                   start_scale=start_scale)
+
+    # create the bdv.xml
+    in_xml = in_path.replace('.n5', '.xml')
+    out_xml = os.path.join('data', f'{vol_name}.xml')
+    write_metadata(in_xml, out_xml, out_path)
+
+    add_to_image_dict(vol_name, layer_type, out_xml)
+
+
+# def convert_bdv_n5(in_path, out_path, use_nested_store, n_threads):
+# add the myosin prospr data
+def add_myosin():
+    print("Add myosin")
+    in_path = os.path.join('/g/arendt/EM_6dpf_segmentation/platy-browser-data/data/0.6.3',
+                           'images/local/prospr-6dpf-1-whole-mhcl4.n5')
+    add_volume(in_path, vol_name='prospr-myosin', layer_type='image')
+
+
+# add the em raw data
+def add_raw():
+    print("Add raw")
+    in_path = '/g/arendt/EM_6dpf_segmentation/platy-browser-data/data/rawdata/sbem-6dpf-1-whole-raw.n5'
+    add_volume(in_path, vol_name='em-raw', layer_type='image', start_scale=3)
+
+
+# add the em cell segmentation
+def add_seg():
+    print("Add cells")
+    in_path = os.path.join('/g/arendt/EM_6dpf_segmentation/platy-browser-data/data/1.0.1',
+                           'images/local/sbem-6dpf-1-whole-segmented-cells.n5')
+    add_volume(in_path, vol_name='em-cells', layer_type='segmentation', start_scale=2)
+
+
+def add_all_volumes():
+    os.makedirs('./data', exist_ok=True)
+    add_myosin()
+    add_raw()
+    add_seg()
+
+
+if __name__ == '__main__':
+    add_all_volumes()

data/em-cells.xml

@@ -0,0 +1,43 @@
+<SpimData version="0.2">
+  <BasePath type="relative">.</BasePath>
+  <SequenceDescription>
+    <ViewSetups>
+      <Attributes name="channel">
+        <Channel>
+          <id>0</id>
+          <name>0</name>
+        </Channel>
+      </Attributes>
+      <ViewSetup>
+        <id>0</id>
+        <name>Setup0</name>
+        <size>3438 3240 2854</size>
+        <voxelSize>
+          <unit>micrometer</unit>
+          <size>0.08 0.08 0.1</size>
+        </voxelSize>
+        <attributes>
+          <channel>0</channel>
+        </attributes>
+      </ViewSetup>
+    </ViewSetups>
+    <Timepoints type="range">
+      <first>0</first>
+      <last>0</last>
+    </Timepoints>
+    <ImageLoader format="bdv.zarr.s3">
+      <Key>em-cells.ome.zarr</Key>
+      <SigningRegion>us-west-2</SigningRegion>
+      <ServiceEndpoint>https://s3.embl.de</ServiceEndpoint>
+      <BucketName>i2k-2020</BucketName>
+      <Authentication>Anonymous</Authentication>
+    </ImageLoader>
+  </SequenceDescription>
+  <ViewRegistrations>
+    <ViewRegistration setup="0" timepoint="0">
+      <ViewTransform type="affine">
+        <affine>0.08 0.0 0.0 0.0 0.0 0.08 0.0 0.0 0.0 0.0 0.1 0.0</affine>
+      </ViewTransform>
+    </ViewRegistration>
+  </ViewRegistrations>
+</SpimData>

data/em-raw.xml

@@ -0,0 +1,43 @@
+<SpimData version="0.2">
+  <BasePath type="relative">.</BasePath>
+  <SequenceDescription>
+    <ViewSetups>
+      <Attributes name="channel">
+        <Channel>
+          <id>1</id>
+          <name>1</name>
+        </Channel>
+      </Attributes>
+      <ViewSetup>
+        <id>0</id>
+        <name>channel 1</name>
+        <size>3438 3240 2854</size>
+        <voxelSize>
+          <unit>micrometer</unit>
+          <size>0.08 0.08 0.1</size>
+        </voxelSize>
+        <attributes>
+          <channel>1</channel>
+        </attributes>
+      </ViewSetup>
+    </ViewSetups>
+    <Timepoints type="range">
+      <first>0</first>
+      <last>0</last>
+    </Timepoints>
+    <ImageLoader format="bdv.zarr.s3">
+      <Key>em-raw.ome.zarr</Key>
+      <SigningRegion>us-west-2</SigningRegion>
+      <ServiceEndpoint>https://s3.embl.de</ServiceEndpoint>
+      <BucketName>i2k-2020</BucketName>
+      <Authentication>Anonymous</Authentication>
+    </ImageLoader>
+  </SequenceDescription>
+  <ViewRegistrations>
+    <ViewRegistration setup="0" timepoint="0">
+      <ViewTransform type="affine">
+        <affine>0.08 0.0 0.0 0.0 0.0 0.08 0.0 0.0 0.0 0.0 0.1 0.0</affine>
+      </ViewTransform>
+    </ViewRegistration>
+  </ViewRegistrations>
+</SpimData>

data/images.json

@@ -0,0 +1,35 @@
+{
+  "em-cells": {
+    "color": "randomFromGlasbey",
+    "contrastLimits": [
+      0.0,
+      1000.0
+    ],
+    "storage": {
+      "remote": "em-cells.xml"
+    },
+    "type": "segmentation"
+  },
+  "em-raw": {
+    "color": "white",
+    "contrastLimits": [
+      0.0,
+      255.0
+    ],
+    "storage": {
+      "remote": "em-raw.xml"
+    },
+    "type": "image"
+  },
+  "prospr-myosin": {
+    "color": "white",
+    "contrastLimits": [
+      0.0,
+      255.0
+    ],
+    "storage": {
+      "remote": "prospr-myosin.xml"
+    },
+    "type": "image"
+  }
+}

data/prospr-myosin.xml

@@ -0,0 +1,43 @@
+<SpimData version="0.2">
+  <BasePath type="relative">.</BasePath>
+  <SequenceDescription>
+    <ViewSetups>
+      <Attributes name="channel">
+        <Channel>
+          <id>0</id>
+          <name>0</name>
+        </Channel>
+      </Attributes>
+      <ViewSetup>
+        <id>0</id>
+        <name>Setup0</name>
+        <size>500 471 519</size>
+        <voxelSize>
+          <unit>micrometer</unit>
+          <size>0.55 0.55 0.55</size>
+        </voxelSize>
+        <attributes>
+          <channel>0</channel>
+        </attributes>
+      </ViewSetup>
+    </ViewSetups>
+    <Timepoints type="range">
+      <first>0</first>
+      <last>0</last>
+    </Timepoints>
+    <ImageLoader format="bdv.zarr.s3">
+      <Key>prospr-myosin.ome.zarr</Key>
+      <SigningRegion>us-west-2</SigningRegion>
+      <ServiceEndpoint>https://s3.embl.de</ServiceEndpoint>
+      <BucketName>i2k-2020</BucketName>
+      <Authentication>Anonymous</Authentication>
+    </ImageLoader>
+  </SequenceDescription>
+  <ViewRegistrations>
+    <ViewRegistration setup="0" timepoint="0">
+      <ViewTransform type="affine">
+        <affine>0.55 0.0 0.0 0.0 0.0 0.55 0.0 0.0 0.0 0.0 0.55 0.0</affine>
+      </ViewTransform>
+    </ViewRegistration>
+  </ViewRegistrations>
+</SpimData>

data_conversion/__init__.py

@@ -0,0 +1 @@
+from .to_ome_zarr import convert_bdv_n5

data_conversion/check_result.py

@@ -0,0 +1,22 @@
+import sys
+import zarr
+
+
+def check_result(path, check_data):
+    with zarr.open(path, mode='r') as f:
+        for name, ds in f.items():
+            shape = ds.shape
+            chunks = ds.chunks
+            assert len(shape) == len(chunks) == 5
+            print(name, shape, chunks)
+
+            if check_data:
+                data = ds[:]
+                # print(data[0, 0, :10, :10, :10])
+                assert data.shape == shape
+
+    print("All tests passed")
+
+
+path = sys.argv[1]
+check_result(path, True)

data_conversion/joshs_script.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+# This assumes that n5-copy has already been used
+
+import argparse
+import zarr
+
+parser = argparse.ArgumentParser()
+parser.add_argument("input")
+parser.add_argument("output")
+ns = parser.parse_args()
+
+zin = zarr.open(ns.input)
+
+sizes = []
+
+def groups(z):
+    rv = sorted(list(z.groups()))
+    assert rv
+    assert not list(z.arrays())
+    return rv
+
+def arrays(z):
+    rv = sorted(list(z.arrays()))
+    assert rv
+    assert not list(z.groups())
+    return rv
+
+setups = groups(zin)
+assert len(setups) == 1  # TODO: multiple channels?
+for sname, setup in setups:
+    timepoints = groups(setup)
+    for tname, timepoint in timepoints:
+        resolutions = arrays(timepoint)
+        for idx, rtuple in enumerate(resolutions):
+            rname, resolution = rtuple
+            try:
+                expected = sizes[idx]
+                assert expected[0] == rname
+                assert expected[1] == resolution.shape
+                assert expected[2] == resolution.chunks
+                assert expected[3] == resolution.dtype
+            except:
+                sizes.append((rname,
+                              resolution.shape,
+                              resolution.chunks,
+                              resolution.dtype))
+
+
+datasets = []
+out = zarr.open(ns.output, mode="w")
+
+for idx, size in enumerate(sizes):
+    name, shape, chunks, dtype = size
+    shape = tuple([len(timepoints), len(setups)] + list(shape))
+    chunks = tuple([1, 1] + list(chunks))
+    a = out.create_dataset(name, shape=shape, chunks=chunks, dtype=dtype)
+    datasets.append({"path": name})
+    for sidx, stuple in enumerate(groups(zin)):
+        for tidx, ttuple in enumerate(groups(stuple[1])):
+            resolutions = arrays(ttuple[1])
+            a[tidx, sidx, :, :, :] = resolutions[idx][1]
+out.attrs["multiscales"] = [
+    {
+        "version": "0.1",
+        "datasets": datasets,
+    }
+]
+