From 506b449f4dc91ef172a18e59230ffc64457ab47b Mon Sep 17 00:00:00 2001 From: Wim Pomp Date: Thu, 29 Jun 2023 14:23:03 +0200 Subject: [PATCH] - read metadata into ome structure - pytest - use pathlib - series as part of the path: path/PosN - summary only shows some available metadata - allow dict in Imread[dict(c=c, z=z, t=t)] - bfread in different process so the user can start another jvm - deal with multiple images (series/positions) in czi files - use jpype instead of javabridge/bioformats - poetry for install --- ndbioimage/.gitignore => .gitignore | 0 README.md | 45 +- ndbioimage/__init__.py | 898 ++++++++++++---------------- ndbioimage/_version.py | 2 - ndbioimage/jvm.py | 38 +- ndbioimage/readers/__init__.py | 5 +- ndbioimage/readers/bfread.py | 246 +++++--- ndbioimage/readers/cziread.py | 479 ++++++++++++--- ndbioimage/readers/ndread.py | 54 +- ndbioimage/readers/seqread.py | 161 ++--- ndbioimage/readers/tifread.py | 77 ++- ndbioimage/transforms.py | 88 +-- pyproject.toml | 37 ++ setup.py | 48 -- tests/test_open.py | 10 + 15 files changed, 1268 insertions(+), 920 deletions(-) rename ndbioimage/.gitignore => .gitignore (100%) delete mode 100644 ndbioimage/_version.py create mode 100644 pyproject.toml delete mode 100644 setup.py create mode 100644 tests/test_open.py diff --git a/ndbioimage/.gitignore b/.gitignore similarity index 100% rename from ndbioimage/.gitignore rename to .gitignore diff --git a/README.md b/README.md index d306ea4..668d564 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,9 @@ -# ndbioimage +# ndbioimage - Work in progress -Exposes (bio) images as a numpy ndarray like object, but without loading the whole -image into memory, reading from the file only when needed. Some metadata is read -and exposed as attributes to the Imread object (TODO: structure data in OME format). -Additionally, it can automatically calculate an affine transform that corrects for -chromatic abberrations etc. and apply it on the fly to the image. +Exposes (bio) images as a numpy ndarray-like-object, but without loading the whole +image into memory, reading from the file only when needed. Some metadata is read and +and stored in an ome structure. Additionally, it can automatically calculate an affine +transform that corrects for chromatic abberrations etc. and apply it on the fly to the image. Currently supports imagej tif files, czi files, micromanager tif sequences and anything bioformats can handle. @@ -13,13 +12,8 @@ bioformats can handle. pip install ndbioimage@git+https://github.com/wimpomp/ndbioimage.git -### With bioformats (if java is properly installed) - - pip install ndbioimage[bioformats]@git+https://github.com/wimpomp/ndbioimage.git - -### With affine transforms (only for python 3.8, 3.9 and 3.10) - - pip install ndbioimage[transforms]@git+https://github.com/wimpomp/ndbioimage.git +Optionally: +https://downloads.openmicroscopy.org/bio-formats/latest/artifacts/bioformats_package.jar ## Usage @@ -27,34 +21,34 @@ bioformats can handle. import matplotlib.pyplot as plt - from ndbioimage import imread - with imread('image_file.tif', axes='ctxy', dtype=int) as im: + from ndbioimage import Imread + with Imread('image_file.tif', axes='ctxy', dtype=int) as im: plt.imshow(im[2, 1]) - Showing some image metadata - from ndbioimage import imread + from ndbioimage import Imread from pprint import pprint - with imread('image_file.tif') as im: + with Imread('image_file.tif') as im: pprint(im) - Slicing the image without loading the image into memory - from ndbioimage import imread - with imread('image_file.tif', axes='cztxy') as im: + from ndbioimage import Imread + with Imread('image_file.tif', axes='cztxy') as im: sliced_im = im[1, :, :, 100:200, 100:200] -sliced_im is an instance of imread which will load any image data from file only when needed +sliced_im is an instance of Imread which will load any image data from file only when needed - Converting (part) of the image to a numpy ndarray - from ndbioimage import imread + from ndbioimage import Imread import numpy as np - with imread('image_file.tif', axes='cztxy') as im: + with Imread('image_file.tif', axes='cztxy') as im: array = np.asarray(im[0, 0]) ## Adding more formats @@ -63,9 +57,8 @@ automatically recognize it and use it to open the appropriate file format. Image subclass Imread and are required to implement the following methods: - staticmethod _can_open(path): return True if path can be opened by this reader -- \_\_metadata__(self): reads metadata from file and adds them to self as attributes, - - the shape of the data in the file needs to be set as self.shape = (X, Y, C, Z, T) - - other attributes like pxsize, acquisitiontime and title can be set here as well +- property ome: reads metadata from file and adds them to an OME object imported +from the ome-types library - \_\_frame__(self, c, z, t): return the frame at channel=c, z-slice=z, time=t from the file Optional methods: @@ -78,5 +71,5 @@ Optional fields: for example: any file handles # TODO -- structure the metadata in OME format tree +- more image formats - re-implement transforms diff --git a/ndbioimage/__init__.py b/ndbioimage/__init__.py index d99e7bc..50a2388 100755 --- a/ndbioimage/__init__.py +++ b/ndbioimage/__init__.py @@ -1,9 +1,12 @@ -import os import re import inspect import pandas import yaml import numpy as np +import multiprocessing +import ome_types +from ome_types import ureg, model +from pint import set_application_registry from datetime import datetime from tqdm.auto import tqdm from itertools import product @@ -14,51 +17,56 @@ from parfor import parfor from tiffwrite import IJTiffFile from numbers import Number from argparse import ArgumentParser -from warnings import warn from ndbioimage.transforms import Transform, Transforms from ndbioimage.jvm import JVM -from ndbioimage._version import __version__, __git_commit_hash__ +from typing import List +from pathlib import Path +from importlib.metadata import version +from traceback import print_exc -class ImTransformsBase(Transforms): - def coords(self, array, colums=None): - if isinstance(array, pandas.DataFrame): - return pandas.concat([self(int(row['C']), int(row['T'])).coords(row, colums) - for _, row in array.iterrows()], axis=1).T - elif isinstance(array, pandas.Series): - return self(int(array['C']), int(array['T'])).coords(array, colums) - else: - raise TypeError('Not a pandas DataFrame or Series.') +try: + __version__ = version(Path(__file__).parent.name) +except (Exception,): + __version__ = 'unknown' + +try: + with open(Path(__file__).parent.parent / '.git' / 'HEAD') as g: + head = g.read().split(':')[1].strip() + with open(Path(__file__).parent.parent / '.git' / head) as h: + __git_commit_hash__ = h.read().rstrip('\n') +except (Exception,): + __git_commit_hash__ = 'unknown' + +ureg.default_format = '~P' +set_application_registry(ureg) -class ImTransforms(ImTransformsBase): - """ Transforms class with methods to calculate channel transforms from bead files etc. - """ - def __init__(self, path, cyllens, tracks=None, detectors=None, file=None, transforms=None): +class ImTransforms(Transforms): + """ Transforms class with methods to calculate channel transforms from bead files etc. """ + def __init__(self, path, cyllens, file=None, transforms=None): super().__init__() self.cyllens = cyllens - self.tracks = tracks - self.detectors = detectors if transforms is None: # TODO: check this - if re.search(r'^Pos\d+', os.path.basename(path.rstrip(os.path.sep))): - self.path = os.path.dirname(os.path.dirname(path)) + if re.search(r'^Pos\d+', path.name): + self.path = path.parent.parent else: - self.path = os.path.dirname(path) + self.path = path.parent if file is not None: if isinstance(file, str) and file.lower().endswith('.yml'): self.ymlpath = file self.beadfile = None else: - self.ymlpath = os.path.join(self.path, 'transform.yml') + self.ymlpath = self.path / 'transform.yml' self.beadfile = file else: - self.ymlpath = os.path.join(self.path, 'transform.yml') + self.ymlpath = self.path / 'transform.yml' self.beadfile = None - self.tifpath = self.ymlpath[:-3] + 'tif' + self.tifpath = self.ymlpath.with_suffix('.tif') try: self.load(self.ymlpath) - except Exception: + except (Exception,): print('No transform file found, trying to generate one.') if not self.files: raise FileNotFoundError('No bead files found to calculate the transform from.') @@ -70,8 +78,17 @@ class ImTransforms(ImTransformsBase): else: # load from dict transforms self.path = path self.beadfile = file - for key, value in transforms.items(): - self[tuple([int(i) for i in key.split(':')])] = Transform(value) + for i, (key, value) in enumerate(transforms.items()): + self[key] = Transform(value) + + def coords_pandas(self, array, cnamelist, colums=None): + if isinstance(array, pandas.DataFrame): + return pandas.concat([self[(cnamelist[int(row['C'])],)].coords(row, colums) + for _, row in array.iterrows()], axis=1).T + elif isinstance(array, pandas.Series): + return self[(cnamelist[int(array['C'])],)].coords(array, colums) + else: + raise TypeError('Not a pandas DataFrame or Series.') @cached_property def files(self): @@ -81,118 +98,110 @@ class ImTransforms(ImTransformsBase): else: files = self.beadfile if isinstance(files, str): + files = (Path(files),) + elif isinstance(files, Path): files = (files,) return files - except Exception: + except (Exception,): return () - def __reduce__(self): - return self.__class__, (self.path, self.cyllens, self.tracks, self.detectors, self.files, self.asdict()) - - def __call__(self, channel, time=None, tracks=None, detectors=None): - tracks = tracks or self.tracks - detectors = detectors or self.detectors - return super().__call__(channel, time, tracks, detectors) - def get_bead_files(self): - files = sorted([os.path.join(self.path, f) for f in os.listdir(self.path) if f.lower().startswith('beads') - and not f.lower().endswith('.pdf') and not f.lower().endswith('pkl')]) + files = sorted([f for f in self.path.parent.iterdir() if f.name.lower().startswith('beads') + and not f.suffix.lower() == '.pdf' and not f.suffix.lower() == 'pkl']) if not files: raise Exception('No bead file found!') - Files = [] + checked_files = [] for file in files: try: - if os.path.isdir(file): - file = os.path.join(file, 'Pos0') - with Imread(file) as im: # check for errors opening the file - pass - Files.append(file) - except Exception: + if file.is_dir(): + file /= 'Pos0' + with Imread(file): # check for errors opening the file + checked_files.append(file) + except (Exception,): continue - if not Files: + if not checked_files: raise Exception('No bead file found!') - return Files + return checked_files def calculate_transform(self, file): """ When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162 - in the horizontal direction - """ - with Imread(file) as im: - ims = [im.max(c) for c in range(im.shape[2])] - goodch = [c for c in range(im.shape[2]) if not im.isnoise(im.max(c))] - untransformed = [c for c in range(im.shape[2]) if self.cyllens[im.detector[c]].lower() == 'none'] + in the horizontal direction """ + with Imread(file, axes='zcxy') as im: + max_ims = im.max('z') + goodch = [c for c, max_im in enumerate(max_ims) if not im.is_noise(max_im)] + if not goodch: + goodch = list(range(len(max_ims))) + untransformed = [c for c in range(im.shape['c']) if self.cyllens[im.detector[c]].lower() == 'none'] good_and_untrans = sorted(set(goodch) & set(untransformed)) if good_and_untrans: masterch = good_and_untrans[0] else: masterch = goodch[0] - print(f'{untransformed = }, {masterch = }, {goodch = }') - C = Transform() + transform = Transform() if not np.any(good_and_untrans): - M = C.matrix - M[0, 0] = 0.86 - C.matrix = M - Tr = Transforms() + matrix = transform.matrix + matrix[0, 0] = 0.86 + transform.matrix = matrix + transforms = Transforms() for c in tqdm(goodch): if c == masterch: - Tr[im.track[c], im.detector[c]] = C + transforms[(im.cnamelist[c],)] = transform else: - Tr[im.track[c], im.detector[c]] = Transform(ims[masterch], ims[c]) * C - return Tr + transforms[(im.cnamelist[c],)] = Transform(max_ims[masterch], max_ims[c]) * transform + return transforms def calculate_transforms(self): - Tq = [self.calculate_transform(file) for file in self.files] - for key in set([key for t in Tq for key in t.keys()]): - T = [t[key] for t in Tq if key in t] - if len(T) == 1: - self[key] = T[0] + transforms = [self.calculate_transform(file) for file in self.files] + for key in set([key for transform in transforms for key in transform.keys()]): + new_transforms = [transform[key] for transform in transforms if key in transform] + if len(new_transforms) == 1: + self[key] = new_transforms[0] else: self[key] = Transform() - self[key].parameters = np.mean([t.parameters for t in T], 0) - self[key].dparameters = (np.std([t.parameters for t in T], 0) / np.sqrt(len(T))).tolist() + self[key].parameters = np.mean([t.parameters for t in new_transforms], 0) + self[key].dparameters = (np.std([t.parameters for t in new_transforms], 0) / + np.sqrt(len(new_transforms))).tolist() def save_transform_tiff(self): - C = 0 + n_channels = 0 for file in self.files: with Imread(file) as im: - C = max(C, im.shape[2]) - with IJTiffFile(self.tifpath, (C, 1, len(self.files))) as tif: + n_channels = max(n_channels, im.shape['c']) + with IJTiffFile(self.tifpath, (n_channels, 1, len(self.files))) as tif: for t, file in enumerate(self.files): with Imread(file) as im: with Imread(file, transform=True) as jm: - for c in range(im.shape[2]): - tif.save(np.hstack((im.max(c), jm.max(c))), c, 0, t) + for c in range(im.shape['c']): + tif.save(np.hstack((im(c=c, t=0).max('z'), jm(c=c, t=0).max('z'))), c, 0, t) -class ImShiftTransforms(ImTransformsBase): +class ImShiftTransforms(Transforms): """ Class to handle drift in xy. The image this is applied to must have a channeltransform already, which is then - replaced by this class. - """ + replaced by this class. """ def __init__(self, im, shifts=None): """ im: Calculate shifts from channel-transformed images im, t x 2 array Sets shifts from array, one row per frame im, dict {frame: shift} Sets shifts from dict, each key is a frame number from where a new shift is applied - im, file Loads shifts from a saved file - """ + im, file Loads shifts from a saved file """ super().__init__() with (Imread(im, transform=True, drift=False) if isinstance(im, str) else im.new(transform=True, drift=False)) as im: self.impath = im.path - self.path = os.path.splitext(self.impath)[0] + '_shifts.txt' + self.path = self.impath.parent / self.impath.stem + '_shifts.txt' self.tracks, self.detectors, self.files = im.track, im.detector, im.beadfile if shifts is not None: if isinstance(shifts, np.ndarray): self.shifts = shifts self.shifts2transforms(im) elif isinstance(shifts, dict): - self.shifts = np.zeros((im.shape[4], 2)) + self.shifts = np.zeros((im.shape['t'], 2)) for k in sorted(shifts.keys()): self.shifts[k:] = shifts[k] self.shifts2transforms(im) elif isinstance(shifts, str): self.load(im, shifts) - elif os.path.exists(self.path): + elif self.path.exists(): self.load(im, self.path) else: self.calulate_shifts(im) @@ -209,9 +218,6 @@ class ImShiftTransforms(ImTransformsBase): else: return Transform() - def __reduce__(self): - return self.__class__, (self.impath, self.shifts) - def load(self, im, file): self.shifts = np.loadtxt(file) self.shifts2transforms(im) @@ -220,32 +226,41 @@ class ImShiftTransforms(ImTransformsBase): self.path = file or self.path np.savetxt(self.path, self.shifts) + def coords(self, array, colums=None): + if isinstance(array, pandas.DataFrame): + return pandas.concat([self(int(row['C']), int(row['T'])).coords(row, colums) + for _, row in array.iterrows()], axis=1).T + elif isinstance(array, pandas.Series): + return self(int(array['C']), int(array['T'])).coords(array, colums) + else: + raise TypeError('Not a pandas DataFrame or Series.') + def calulate_shifts0(self, im): """ Calculate shifts relative to the first frame """ im0 = im[:, 0, 0].squeeze().transpose(2, 0, 1) - @parfor(range(1, im.shape[4]), (im, im0), desc='Calculating image shifts.') + @parfor(range(1, im.shape['t']), (im, im0), desc='Calculating image shifts.') def fun(t, im, im0): return Transform(im0, im[:, 0, t].squeeze().transpose(2, 0, 1), 'translation') transforms = [Transform()] + fun self.shifts = np.array([t.parameters[4:] for t in transforms]) - self.setTransforms(transforms, im.transform) + self.set_transforms(transforms, im.transform) def calulate_shifts(self, im): """ Calculate shifts relative to the previous frame """ - @parfor(range(1, im.shape[4]), (im,), desc='Calculating image shifts.') + @parfor(range(1, im.shape['t']), (im,), desc='Calculating image shifts.') def fun(t, im): return Transform(im[:, 0, t-1].squeeze().transpose(2, 0, 1), im[:, 0, t].squeeze().transpose(2, 0, 1), 'translation') transforms = [Transform()] + fun self.shifts = np.cumsum([t.parameters[4:] for t in transforms]) - self.setTransforms(transforms, im.transform) + self.set_transforms(transforms, im.transform) def shifts2transforms(self, im): - self.setTransforms([Transform(np.array(((1, 0, s[0]), (0, 1, s[1]), (0, 0, 1)))) - for s in self.shifts], im.transform) + self.set_transforms([Transform(np.array(((1, 0, s[0]), (0, 1, s[1]), (0, 0, 1)))) + for s in self.shifts], im.transform) - def setTransforms(self, shift_transforms, channel_transforms): + def set_transforms(self, shift_transforms, channel_transforms): for key, value in channel_transforms.items(): for t, T in enumerate(shift_transforms): self[key[0], key[1], t] = T * channel_transforms[key] @@ -270,10 +285,8 @@ class DequeDict(OrderedDict): self.__truncate__() -def tolist(item): - if isinstance(item, XmlData): - return [item] - elif hasattr(item, 'items'): +def tolist(item) -> List: + if hasattr(item, 'items'): return item elif isinstance(item, str): return [item] @@ -284,6 +297,51 @@ def tolist(item): return list((item,)) +def find(obj, **kwargs): + for item in obj: + if all([getattr(item, key) == value for key, value in kwargs.items()]): + return item + + +def find_rec(obj, **kwargs): + if isinstance(obj, list): + for item in obj: + ans = find_rec(item, **kwargs) + if ans: + return ans + elif isinstance(obj, ome_types._base_type.OMEType): + if all([hasattr(obj, key) for key in kwargs.keys()]) and all( + [getattr(obj, key) == value for key, value in kwargs.items()]): + return obj + for k, v in obj.__dict__.items(): + ans = find_rec(v, **kwargs) + if ans: + return ans + + +def try_default(fun, default, *args, **kwargs): + try: + return fun(*args, **kwargs) + except (Exception,): + return default + + +def ome_subprocess(path): + try: + jvm = JVM() + ome_meta = jvm.metadata_tools.createOMEXMLMetadata() + reader = jvm.image_reader() + reader.setMetadataStore(ome_meta) + reader.setId(str(path)) + ome = ome_types.from_xml(str(ome_meta.dumpXML()), parser='lxml') + except (Exception,): + print_exc() + ome = model.OME() + finally: + jvm.kill_vm() + return ome + + class Shape(tuple): def __new__(cls, shape, axes='xyczt'): instance = super().__new__(cls, shape) @@ -303,185 +361,6 @@ class Shape(tuple): return tuple(self[i] for i in 'xyczt') -class XmlData(OrderedDict): - def __init__(self, elem=None): - super(XmlData, self).__init__() - if elem: - if isinstance(elem, dict): - self.update(elem) - else: - self.update(XmlData._todict(elem)[1]) - - def re_search(self, reg, default=None, *args, **kwargs): - return tolist(XmlData._output(XmlData._search(self, reg, True, default, *args, **kwargs)[1])) - - def search(self, key, default=None): - return tolist(XmlData._output(XmlData._search(self, key, False, default)[1])) - - def re_search_all(self, reg, *args, **kwargs): - K, V = XmlData._search_all(self, reg, True, *args, **kwargs) - return {k: XmlData._output(v) for k, v in zip(K, V)} - - def search_all(self, key): - K, V = XmlData._search_all(self, key, False) - return {k: XmlData._output(v) for k, v in zip(K, V)} - - @staticmethod - def _search(d, key, regex=False, default=None, *args, **kwargs): - if isinstance(key, (list, tuple)): - if len(key) == 1: - key = key[0] - else: - for v in XmlData._search_all(d, key[0], regex, *args, **kwargs)[1]: - found, value = XmlData._search(v, key[1:], regex, default, *args, **kwargs) - if found: - return True, value - return False, default - - if hasattr(d, 'items'): - for k, v in d.items(): - if isinstance(k, str): - if (not regex and k == key) or (regex and re.findall(key, k, *args, **kwargs)): - return True, v - elif isinstance(v, dict): - found, value = XmlData._search(v, key, regex, default, *args, **kwargs) - if found: - return True, value - elif isinstance(v, (list, tuple)): - for w in v: - found, value = XmlData._search(w, key, regex, default, *args, **kwargs) - if found: - return True, value - else: - found, value = XmlData._search(v, key, regex, default, *args, **kwargs) - if found: - return True, value - return False, default - - @staticmethod - def _search_all(d, key, regex=False, *args, **kwargs): - K = [] - V = [] - if hasattr(d, 'items'): - for k, v in d.items(): - if isinstance(k, str): - if (not regex and k == key) or (regex and re.findall(key, k, *args, **kwargs)): - K.append(k) - V.append(v) - elif isinstance(v, dict): - q, w = XmlData._search_all(v, key, regex, *args, **kwargs) - K.extend([str(k) + '|' + i for i in q]) - V.extend(w) - elif isinstance(v, (list, tuple)): - for j, val in enumerate(v): - q, w = XmlData._search_all(val, key, regex, *args, **kwargs) - K.extend([str(k) + '|' + str(j) + '|' + i for i in q]) - V.extend(w) - else: - q, w = XmlData._search_all(v, key, regex, *args, **kwargs) - K.extend([str(k) + '|' + i for i in q]) - V.extend(w) - return K, V - - @staticmethod - def _enumdict(d): - d2 = {} - for k, v in d.items(): - idx = [int(i) for i in re.findall(r'(?<=:)\d+$', k)] - if idx: - key = re.findall(r'^.*(?=:\d+$)', k)[0] - if key not in d2: - d2[key] = {} - d2[key][idx[0]] = d['{}:{}'.format(key, idx[0])] - else: - d2[k] = v - rec = False - for k, v in d2.items(): - if [int(i) for i in re.findall(r'(?<=:)\d+$', k)]: - rec = True - break - if rec: - return XmlData._enumdict(d2) - else: - return d2 - - @staticmethod - def _unique_children(l): - if l: - keys, values = zip(*l) - d = {} - for k in set(keys): - value = [v for m, v in zip(keys, values) if k == m] - if len(value) == 1: - d[k] = value[0] - else: - d[k] = value - return d - else: - return {} - - @staticmethod - def _todict(elem): - d = {} - if hasattr(elem, 'Key') and hasattr(elem, 'Value'): - name = elem.Key.cdata - d = elem.Value.cdata - return name, d - - if hasattr(elem, '_attributes') and elem._attributes is not None and 'ID' in elem._attributes: - name = elem._attributes['ID'] - elem._attributes.pop('ID') - elif hasattr(elem, '_name'): - name = elem._name - else: - name = 'none' - - if name == 'Value': - if hasattr(elem, 'children') and len(elem.children): - return XmlData._todict(elem.children[0]) - - if hasattr(elem, 'children'): - children = [XmlData._todict(child) for child in elem.children] - children = XmlData._unique_children(children) - if children: - d = OrderedDict(d, **children) - if hasattr(elem, '_attributes'): - children = elem._attributes - if children: - d = OrderedDict(d, **children) - if not len(d.keys()) and hasattr(elem, 'cdata'): - return name, elem.cdata - - return name, XmlData._enumdict(d) - - @staticmethod - def _output(s): - if isinstance(s, dict): - return XmlData(s) - elif isinstance(s, (tuple, list)): - return [XmlData._output(i) for i in s] - elif not isinstance(s, str): - return s - elif len(s) > 1 and s[0] == '[' and s[-1] == ']': - return [XmlData._output(i) for i in s[1:-1].split(', ')] - elif re.search(r'^[-+]?\d+$', s): - return int(s) - elif re.search(r'^[-+]?\d?\d*\.?\d+([eE][-+]?\d+)?$', s): - return float(s) - elif s.lower() == 'true': - return True - elif s.lower() == 'false': - return False - elif s.lower() == 'none': - return None - else: - return s - - def __getitem__(self, item): - value = super().__getitem__(item) - return XmlData(value) if isinstance(value, dict) else value - - class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): """ class to read image files, while taking good care of important metadata, currently optimized for .czi files, but can open anything that bioformats can handle @@ -535,8 +414,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): optional close(self): close the file in a proper way optional field priority: subclasses with lower priority will be tried first, default = 99 Any other method can be overridden as needed - wp@tl2019-2021 - """ + wp@tl2019-2023 """ # TODO: more numpy.ndarray methods as needed # TODO: imread.std @@ -545,20 +423,21 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): priority = 99 do_not_pickle = 'base', 'copies', 'cache' do_not_copy = 'extrametadata' + ureg = ureg @staticmethod @abstractmethod def _can_open(path): # Override this method, and return true when the subclass can open the file return False - @abstractmethod - def __metadata__(self): - return - @abstractmethod def __frame__(self, c, z, t): # Override this, return the frame at c, z, t return np.random.randint(0, 255, self.shape['xy']) + @cached_property + def ome(self): + return self.get_ome(self.path) + def open(self): # Optionally override this, open file handles etc. """ filehandles cannot be pickled and should be marked such by setting do_not_pickle = 'file_handle_name' """ return @@ -566,6 +445,13 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): def close(self): # Optionally override this, close file handles etc. return + @staticmethod + def get_ome(path): + """ Use java BioFormats to make an ome metadata structure. """ + with multiprocessing.get_context('spawn').Pool(1) as pool: + ome = pool.map(ome_subprocess, (path,))[0] + return ome + def __new__(cls, path=None, *args, **kwargs): if cls is not Imread: return super().__new__(cls) @@ -573,7 +459,8 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): raise Exception('Restart python kernel please!') if isinstance(path, Imread): return path - for subclass in sorted(cls.__subclasses__(), key=lambda subclass: subclass.priority): + path, _ = Imread.split_path_series(path) + for subclass in sorted(cls.__subclasses__(), key=lambda subclass_: subclass_.priority): if subclass._can_open(path): do_not_pickle = (cls.do_not_pickle,) if isinstance(cls.do_not_pickle, str) else cls.do_not_pickle subclass_do_not_pickle = (subclass.do_not_pickle,) if isinstance(subclass.do_not_pickle, str) \ @@ -587,8 +474,15 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): return super().__new__(subclass) - def __init__(self, path, series=0, transform=False, drift=False, beadfile=None, sigma=None, dtype=None, - axes='cztxy'): + @staticmethod + def split_path_series(path): + if isinstance(path, str): + path = Path(path) + if isinstance(path, Path) and path.name.startswith('Pos'): + return path.parent, int(path.name.lstrip('Pos')) + return path, 0 + + def __init__(self, path, transform=False, drift=False, beadfile=None, dtype=None, axes=None): if isinstance(path, Imread): return self.isclosed = False @@ -596,48 +490,88 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self.base = None self.copies = [] if isinstance(path, str): - self.path = os.path.abspath(path) - self.title = os.path.splitext(os.path.basename(self.path))[0] - self.acquisitiondate = datetime.fromtimestamp(os.path.getmtime(self.path)).strftime('%y-%m-%dT%H:%M:%S') - else: - self.path = path # ndarray + path = Path(path) + self.path, self.series = self.split_path_series(path) + if isinstance(path, Path): + self.title = self.path.name + self.acquisitiondate = datetime.fromtimestamp(self.path.stat().st_mtime).strftime('%y-%m-%dT%H:%M:%S') + else: # ndarray + self.title = 'ndarray' + self.acquisitiondate = 'now' self.transform = transform self.drift = drift self.beadfile = beadfile self.dtype = dtype - self.series = series - self.pxsize = 1e-1 - self.settimeinterval = 0 - self.pxsizecam = 0 - self.magnification = 0 - self.exposuretime = (0,) - self.deltaz = 1 - self.pcf = (1, 1) - self.laserwavelengths = [[]] - self.laserpowers = [[]] - self.powermode = 'normal' - self.optovar = (1,) - self.binning = 1 - self.collimator = (1,) - self.tirfangle = (0,) - self.gain = (100, 100) - self.objective = 'unknown' - self.filter = 'unknown' - self.NA = 1 + self.reader = None + + self.pcf = None + self.powermode = None + self.collimator = None + self.tirfangle = None self.cyllens = ['None', 'None'] self.duolink = 'None' self.detector = [0, 1] self.track = [0] - self.metadata = {} self.cache = DequeDict(16) self._frame_decorator = None self.frameoffset = 0, 0 # how far apart the centers of frame and sensor are self.open() - self.__metadata__() + + # extract some metadata from ome + instrument = self.ome.instruments[0] + image = self.ome.images[0] + pixels = image.pixels + self.shape = pixels.size_x, pixels.size_y, pixels.size_c, pixels.size_z, pixels.size_t + self.pxsize = pixels.physical_size_x_quantity + self.exposuretime = tuple(find(image.pixels.planes, the_c=c).exposure_time_quantity + for c in range(self.shape['c'])) + if self.zstack: + self.deltaz = image.pixels.physical_size_z_quantity + self.deltaz_um = None if self.deltaz is None else self.deltaz.to(self.ureg.um).m + else: + self.deltaz = self.deltaz_um = None + if self.ome.images[0].objective_settings: + self.objective = find(instrument.objectives, id=self.ome.images[0].objective_settings.id) + else: + self.objective = None + self.timeval = [find(image.pixels.planes, the_c=0, the_t=t, the_z=0).delta_t for t in range(self.shape['t'])] + self.timeinterval = np.diff(self.timeval).mean() if len(self.timeval) > 1 else 0 + try: + self.binning = [int(i) for i in image.pixels.channels[0].detector_settings.binning.value.split('x')] + except (Exception,): + self.binning = None + self.cnamelist = [channel.name for channel in image.pixels.channels] + optovars = [objective for objective in instrument.objectives if 'tubelens' in objective.id.lower()] + if len(optovars) == 0: + self.tubelens = None + else: + self.tubelens = optovars[0] + if self.objective: + if self.tubelens: + self.magnification = self.objective.nominal_magnification * self.tubelens.nominal_magnification + else: + self.magnification = self.objective.nominal_magnification + self.NA = self.objective.lens_na + else: + self.magnification = None + self.NA = None + + self.gain = [find(instrument.detectors, id=channel.detector_settings.id).amplification_gain + for channel in image.pixels.channels + if channel.detector_settings + and find(instrument.detectors, id=channel.detector_settings.id).amplification_gain] + self.laserwavelengths = [(channel.emission_wavelength_quantity.to(self.ureg.nm).m,) + for channel in pixels.channels if channel.emission_wavelength_quantity] + self.laserpowers = try_default(lambda channel: [(1 - channel.light_source_settings.attenuation,) + for channel in pixels.channels], []) + self.filter = try_default(lambda: [find(instrument.filter_sets, id=channel.filter_set_ref.id).model + for channel in image.pixels.channels], None) + self.pxsize_um = None if self.pxsize is None else self.pxsize.to(self.ureg.um).m + self.exposuretime_s = [None if i is None else i.to(self.ureg.s).m for i in self.exposuretime] self.file_shape = self.shape.xyczt - if axes.lower() == 'squeeze': + if axes is None: self.axes = ''.join(i for i in 'cztxy' if self.shape[i] > 1) elif axes.lower() == 'full': self.axes = 'cztxy' @@ -645,9 +579,6 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self.axes = axes self.slice = [np.arange(s, dtype=int) for s in self.shape.xyczt] - if not hasattr(self, 'cnamelist'): - self.cnamelist = 'abcdefghijklmnopqrstuvwxyz'[:self.shape['c']] - m = self.extrametadata if m is not None: try: @@ -657,35 +588,34 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self.feedback = m['FeedbackChannels'] else: self.feedback = m['FeedbackChannel'] - except Exception: + except (Exception,): self.cyllens = ['None', 'None'] self.duolink = 'None' self.feedback = [] try: self.cyllenschannels = np.where([self.cyllens[self.detector[c]].lower() != 'none' for c in range(self.shape['c'])])[0].tolist() - except Exception: + except (Exception,): pass self.set_transform() try: - s = int(re.findall(r'_(\d{3})_', self.duolink)[0]) - except Exception: - s = 561 - if sigma is None: - try: - sigma = [] - for t, d in zip(self.track, self.detector): - l = np.array(self.laserwavelengths[t]) + 22 - sigma.append([l[l > s].max(initial=0), l[l < s].max(initial=0)][d]) - sigma = np.array(sigma) - sigma[sigma == 0] = 600 - sigma /= 2 * self.NA * self.pxsize * 1000 - self.sigma = sigma.tolist() - except Exception: - self.sigma = [2] * self.shape['c'] - else: - self.sigma = sigma - if 1.5 < self.NA: + s = int(re.findall(r'_(\d{3})_', self.duolink)[0]) * ureg.nm + except (Exception,): + s = 561 * ureg.nm + try: + sigma = [] + for c, d in enumerate(self.detector): + emission = np.hstack(self.laserwavelengths[c]) + 22 * ureg.nm + sigma.append([emission[emission > s].max(initial=0), emission[emission < s].max(initial=0)][d]) + sigma = np.hstack(sigma) + sigma[sigma == 0] = 600 * ureg.nm + sigma /= 2 * self.NA * self.pxsize + self.sigma = sigma.magnitude.tolist() + except (Exception,): + self.sigma = [2] * self.shape['c'] + if not self.NA: + self.immersionN = 1 + elif 1.5 < self.NA: self.immersionN = 1.661 elif 1.3 < self.NA < 1.5: self.immersionN = 1.518 @@ -711,9 +641,9 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self.transform = self.transform else: if isinstance(self.transform, str): - self.transform = ImTransforms(self.path, self.cyllens, self.track, self.detector, self.transform) + self.transform = ImTransforms(self.path, self.cyllens, self.transform) else: - self.transform = ImTransforms(self.path, self.cyllens, self.track, self.detector, self.beadfile) + self.transform = ImTransforms(self.path, self.cyllens, self.beadfile) if self.drift is True: self.transform = ImShiftTransforms(self) elif not (self.drift is False or self.drift is None): @@ -734,8 +664,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): @staticmethod def get_config(file): - """ Open a yml parameter file - """ + """ Open a yml config file """ loader = yaml.SafeLoader loader.add_implicit_resolver( r'tag:yaml.org,2002:float', @@ -763,73 +692,72 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self._frame_decorator = decorator self.cache = DequeDict(self.cache.maxlen) - def __repr__(self): - """ gives a helpful summary of the recorded experiment - """ - s = [100 * '#'] - s.append('path/filename: {}'.format(self.path)) - s.append('shape ({}): '.format(self.axes) + ' ' * (5 - self.ndim) + - ' x '.join(('{}',) * self.ndim).format(*self.shape)) - s.append('pixelsize: {:.2f} nm'.format(self.pxsize * 1000)) - if self.zstack: - s.append('z-interval: {:.2f} nm'.format(self.deltaz * 1000)) - s.append('Exposuretime: ' + ('{:.2f} ' * len(self.exposuretime)).format( - *(np.array(self.exposuretime) * 1000)) + 'ms') - if self.timeseries: - if self.timeval and np.diff(self.timeval).shape[0]: - s.append('t-interval: {:.3f} ± {:.3f} s'.format( - np.diff(self.timeval).mean(), np.diff(self.timeval).std())) - else: - s.append('t-interval: {:.2f} s'.format(self.settimeinterval)) - s.append('binning: {}x{}'.format(self.binning, self.binning)) - s.append('laser colors: ' + ' | '.join([' & '.join(len(l)*('{:.0f}',)).format(*l) - for l in self.laserwavelengths]) + ' nm') - s.append('laser powers: ' + ' | '.join([' & '.join(len(l)*('{}',)).format(*[100 * i for i in l]) - for l in self.laserpowers]) + ' %') - s.append('objective: {}'.format(self.objective)) - s.append('magnification: {}x'.format(self.magnification)) - s.append('optovar: ' + (' {}' * len(self.optovar)).format(*self.optovar) + 'x') - s.append('filterset: {}'.format(self.filter)) - s.append('powermode: {}'.format(self.powermode)) - s.append('collimator: ' + (' {}' * len(self.collimator)).format(*self.collimator)) - s.append('TIRF angle: ' + (' {:.2f}°' * len(self.tirfangle)).format(*self.tirfangle)) - s.append('gain: ' + (' {:.0f}' * len(self.gain)).format(*self.gain)) - s.append('pcf: ' + (' {:.2f}' * len(self.pcf)).format(*self.pcf)) + @property + def summary(self): + """ gives a helpful summary of the recorded experiment """ + s = [f"path/filename: {self.path}", + f"reader: {self.__class__.__module__.split('.')[-1]}", + f"shape ({self.axes}):".ljust(15) + f"{' x '.join(str(i) for i in self.shape)}"] + if self.pxsize_um: + s.append(f'pixel size: {1000 * self.pxsize_um:.2f} nm') + if self.zstack and self.deltaz_um: + s.append(f'z-interval: {1000 * self.deltaz_um:.2f} nm') + if self.exposuretime_s[0]: + s.append(f'exposuretime: {self.exposuretime_s[0]:.2f} s') + if self.timeseries and self.timeinterval: + s.append(f'time interval: {self.timeinterval:.3f} s') + if self.binning: + s.append('binning: {}x{}'.format(*self.binning)) + if self.laserwavelengths: + s.append('laser colors: ' + ' | '.join([' & '.join(len(w)*('{:.0f}',)).format(*w) + for w in self.laserwavelengths]) + ' nm') + if self.laserpowers: + s.append('laser powers: ' + ' | '.join([' & '.join(len(p)*('{:.3g}',)).format(*[100 * i for i in p]) + for p in self.laserpowers]) + ' %') + if self.objective: + s.append('objective: {}'.format(self.objective.model)) + if self.magnification: + s.append('magnification: {}x'.format(self.magnification)) + if self.tubelens: + s.append('tubelens: {}'.format(self.tubelens.model)) + if self.filter: + s.append('filterset: {}'.format(self.filter)) + if self.powermode: + s.append('powermode: {}'.format(self.powermode)) + if self.collimator: + s.append('collimator: ' + (' {}' * len(self.collimator)).format(*self.collimator)) + if self.tirfangle: + s.append('TIRF angle: ' + (' {:.2f}°' * len(self.tirfangle)).format(*self.tirfangle)) + if self.gain: + s.append('gain: ' + (' {:.0f}' * len(self.gain)).format(*self.gain)) + if self.pcf: + s.append('pcf: ' + (' {:.2f}' * len(self.pcf)).format(*self.pcf)) return '\n'.join(s) def __str__(self): - return self.path + return str(self.path) def __len__(self): return self.shape[0] - def __call__(self, *n): - """ returns single 2D frame - im(n): index linearly in czt order - im(c,z): return im(c,z,t=0) - im(c,z,t): return im(c,z,t) - """ - if len(n) == 1: - n = self.get_channel(n[0]) - c = int(n % self.shape['c']) - z = int((n // self.shape['c']) % self.shape['z']) - t = int((n // (self.shape['c'] * self.shape['z'])) % self.shape['t']) - return self.frame(c, z, t) - else: - return self.frame(*[int(i) for i in n]) + def __call__(self, c=None, z=None, t=None, x=None, y=None): + """ same as im[] but allowing keyword axes, but slices need to made with slice() or np.s_ """ + return self[{k: slice(v) if v is None else v for k, v in dict(c=c, z=z, t=t, x=x, y=y).items()}] def __getitem__(self, n): - # None = : - if isinstance(n, (slice, Number)): + """ slice like a numpy array but return an Imread instance """ + if isinstance(n, (slice, Number)): # None = : n = (n,) elif isinstance(n, type(Ellipsis)): n = (None,) * len(self.axes) + elif isinstance(n, dict): # allow im[dict(z=0)] etc. + n = [n.get(i, slice(None)) for i in self.axes] n = list(n) # deal with ... ell = [i for i, e in enumerate(n) if isinstance(e, type(Ellipsis))] if len(ell) > 1: - raise IndexError("an index can only have a single ellipsis (...)") + raise IndexError('an index can only have a single ellipsis (...)') if len(ell): if len(n) > self.ndim: n.remove(Ellipsis) @@ -840,8 +768,8 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): while len(n) < self.ndim: n.append(None) - T = [self.shape.axes.find(i) for i in 'xyczt'] - n = [n[j] if 0 <= j < len(n) else None for j in T] # reorder n + axes_idx = [self.shape.axes.find(i) for i in 'xyczt'] + n = [n[j] if 0 <= j < len(n) else None for j in axes_idx] # reorder n new_slice = [] for s, e in zip(self.slice, n): @@ -862,11 +790,11 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): return new def __contains__(self, item): - def unique_yield(l, i): - for k in l: + def unique_yield(a, b): + for k in a: yield k - for k in i: - if k not in l: + for k in b: + if k not in a: yield k for idx in unique_yield(list(self.cache.keys()), product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t']))): @@ -878,15 +806,15 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): def __array__(self, dtype=None): block = self.block(*self.slice) - T = [self.shape.axes.find(i) for i in 'xyczt'] - S = tuple({i for i, j in enumerate(T) if j == -1}.union( + axes_idx = [self.shape.axes.find(i) for i in 'xyczt'] + axes_squeeze = tuple({i for i, j in enumerate(axes_idx) if j == -1}.union( {i for i, j in enumerate(self.slice) if isinstance(j, Number)})) - block = block.squeeze(S) + block = block.squeeze(axes_squeeze) if dtype is not None: block = block.astype(dtype) if block.ndim == 0: return block.item() - axes = ''.join(j for i, j in enumerate('xyczt') if i not in S) + axes = ''.join(j for i, j in enumerate('xyczt') if i not in axes_squeeze) return block.transpose([axes.find(i) for i in self.shape.axes if i in axes]) def asarray(self): @@ -925,7 +853,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): self.copies = [] self.cache = DequeDict(16) - # TODO: this is causing problems when multiprocessing + # TODO: this is causing problems when multiprocessing and doesn't work anyway # def __del__(self): # if not self.copies: # if self.base is None: @@ -941,7 +869,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): new.copies = [] new.base = self for key, value in self.__dict__.items(): - if not key in self.do_not_copy and not hasattr(new, key): + if key not in self.do_not_copy and not hasattr(new, key): new.__dict__[key] = value self.copies.append(new) return new @@ -1008,46 +936,11 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): def moveaxis(self, source, destination): raise NotImplementedError('moveaxis it not implemented') - def czt(self, n): - """ returns indices c, z, t used when calling im(n) - """ - if not isinstance(n, tuple): - c = n % self.shape['c'] - z = (n // self.shape['c']) % self.shape['z'] - t = (n // (self.shape['c'] * self.shape['z'])) % self.shape['t'] - return c, z, t - n = list(n) - if len(n) == 2 or len(n) == 4: - n.append(slice(0, -1, 1)) - if len(n) == 3: - n = list(n) - for i, (ax, e) in enumerate(zip('czt', n)): - if isinstance(e, slice): - a = [e.start, e.stop, e.step] - if a[0] is None: - a[0] = 0 - if a[1] is None: - a[1] = -1 - if a[2] is None: - a[2] = 1 - for j in range(2): - if a[j] < 0: - a[j] %= self.shape[ax] - a[j] += 1 - n[i] = np.arange(*a) - n = [np.array(i) for i in n] - return tuple(n) - if len(n) == 5: - return tuple(n[2:5]) - - def czt2n(self, c, z, t): - return c + z * self.shape['c'] + t * self.shape['c'] * self.shape['z'] - def transform_frame(self, frame, c, t=0): if self.transform is None: return frame else: - return self.transform(c, t, self.track, self.detector).frame(frame) + return self.transform[(self.cnamelist[c],)].frame(frame) def get_czt(self, c, z, t): czt = [] @@ -1068,35 +961,6 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): czt.append([k % self.shape[i] for k in n]) return [self.get_channel(c) for c in czt[0]], *czt[1:] - def _stats(self, fun, c=None, z=None, t=None, ffun=None): - """ fun = np.min, np.max, np.sum or their nan varieties """ - warn('Warning: _stats is deprecated.') - c, z, t = self.get_czt(c, z, t) - if fun in (np.min, np.nanmin): - val = np.inf - elif fun in (np.max, np.nanmax): - val = -np.inf - else: - val = 0 - if ffun is None: - ffun = lambda im: im - T = np.full(self.shape['xy'], val, self.dtype) - for ic in c: - m = np.full(self.shape['xy'], val, self.dtype) - if isinstance(self.transform, ImShiftTransforms): - for it in t: - n = np.full(self.shape['xy'], val, self.dtype) - for iz in z: - n = fun((n, ffun(self.__frame__(ic, iz, it))), 0) - m = self.transform_frame(n, ic, it) - else: - for it, iz in product(t, z): - m = fun((m, ffun(self.__frame__(ic, iz, it))), 0) - if isinstance(self.transform, ImTransforms): - m = self.transform_frame(m, ic, 0) - T = fun((T, m), 0) - return T - @staticmethod def __array_arg_fun__(self, fun, axis=None, out=None): """ frame-wise application of np.argmin and np.argmax """ @@ -1167,9 +1031,10 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): return Imread.__array_arg_fun__(self, np.argmax, *args, **kwargs) def __array_fun__(self, fun, axis=None, dtype=None, out=None, keepdims=False, initial=None, where=True, ffun=None): - """ frame-wise application of np.min, np.max, np.sum, np.mean and theis nan equivalents """ + """ frame-wise application of np.min, np.max, np.sum, np.mean and their nan equivalents """ if ffun is None: - ffun = lambda im: im + def ffun(im): + return im if dtype is None: dtype = self.dtype if out is None else out.dtype @@ -1275,12 +1140,13 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): def reshape(self, *args, **kwargs): return np.asarray(self).reshape(*args, **kwargs) - @property - def is_noise(self): + def is_noise(self, volume=None): """ True if volume only has noise """ - F = np.fft.fftn(self) - S = np.fft.fftshift(np.fft.ifftn(F * F.conj()).real / np.sum(self ** 2)) - return -np.log(1 - S[tuple([0] * S.ndim)]) > 5 + if volume is None: + volume = self + fft = np.fft.fftn(volume) + corr = np.fft.fftshift(np.fft.ifftn(fft * fft.conj()).real / np.sum(volume ** 2)) + return -np.log(1 - corr[tuple([0] * corr.ndim)]) > 5 @property def dtype(self): @@ -1300,8 +1166,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): return c[0] def frame(self, c=0, z=0, t=0): - """ returns single 2D frame - """ + """ returns single 2D frame """ c = self.get_channel(c) c %= self.file_shape[2] z %= self.file_shape[3] @@ -1322,14 +1187,12 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): return f.copy() def data(self, c=0, z=0, t=0): - """ returns 3D stack of frames - """ + """ returns 3D stack of frames """ c, z, t = [np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) for i, e in zip('czt', (c, z, t))] return np.dstack([self.frame(ci, zi, ti) for ci, zi, ti in product(c, z, t)]) def block(self, x=None, y=None, c=None, z=None, t=None): - """ returns 5D block of frames - """ + """ returns 5D block of frames """ x, y, c, z, t = [np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) for i, e in zip('xyczt', (x, y, c, z, t))] d = np.full((len(x), len(y), len(c), len(z), len(t)), np.nan, self.dtype) @@ -1337,24 +1200,9 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): d[:, :, ci, zi, ti] = self.frame(cj, zj, tj)[x][:, y] return d - @cached_property - def timeval(self): - if hasattr(self, 'metadata') and isinstance(self.metadata, XmlData): - image = self.metadata.search('Image') - if (isinstance(image, dict) and self.series in image) or (isinstance(image, list) and len(image)): - image = XmlData(image[0]) - return sorted(np.unique(image.search_all('DeltaT').values()))[:self.shape['t']] - else: - return (np.arange(self.shape['t']) * self.settimeinterval).tolist() - - @cached_property - def timeinterval(self): - return float(np.diff(self.timeval).mean()) if len(self.timeval) > 1 else 1 - @cached_property def piezoval(self): - """ gives the height of the piezo and focus motor, only available when CylLensGUI was used - """ + """ gives the height of the piezo and focus motor, only available when CylLensGUI was used """ def upack(idx): time = list() val = list() @@ -1407,36 +1255,36 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): @cached_property def extrametadata(self): - if isinstance(self.path, str) and len(self.path) > 3: - if os.path.isfile(self.path[:-3] + 'pzl2'): - pname = self.path[:-3] + 'pzl2' - elif os.path.isfile(self.path[:-3] + 'pzl'): - pname = self.path[:-3] + 'pzl' + if isinstance(self.path, Path): + if self.path.with_suffix('.pzl2').exists(): + pname = self.path.with_suffix('.pzl2') + elif self.path.with_suffix('.pzl').exists(): + pname = self.path.with_suffix('.pzl') else: return try: return self.get_config(pname) - except Exception: + except (Exception,): return return def save_as_tiff(self, fname=None, c=None, z=None, t=None, split=False, bar=True, pixel_type='uint16'): """ saves the image as a tiff-file - split: split channels into different files - """ + split: split channels into different files """ if fname is None: - if isinstance(self.path, str): - fname = self.path[:-3] + 'tif' - else: - raise Exception('No filename given.') - elif not fname[-3:] == 'tif': - fname += '.tif' + fname = self.path.with_suffix('.tif') + if fname == self.path: + raise FileExistsError(f'File {fname} exists already.') + if not isinstance(fname, Path): + fname = Path(fname) if split: for i in range(self.shape['c']): if self.timeseries: - self.save_as_tiff(fname[:-3] + '_C{:01d}.tif'.format(i), i, 0, None, False, bar, pixel_type) + self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, 0, None, False, + bar, pixel_type) else: - self.save_as_tiff(fname[:-3] + '_C{:01d}.tif'.format(i), i, None, 0, False, bar, pixel_type) + self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, None, 0, False, + bar, pixel_type) else: n = [c, z, t] for i, ax in enumerate('czt'): @@ -1447,16 +1295,16 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, metaclass=ABCMeta): shape = [len(i) for i in n] at_least_one = False - with IJTiffFile(fname, shape, pixel_type, pxsize=self.pxsize, deltaz=self.deltaz) as tif: + with IJTiffFile(fname.with_suffix('.tif'), shape, pixel_type, + pxsize=self.pxsize_um, deltaz=self.deltaz_um) as tif: for i, m in tqdm(zip(product(*[range(s) for s in shape]), product(*n)), total=np.prod(shape), desc='Saving tiff', disable=not bar): if np.any(self(*m)) or not at_least_one: tif.save(self(*m), *i) at_least_one = True - @cached_property - def summary(self): - return self.__repr__() + def __repr__(self): + return self.summary def main(): @@ -1471,19 +1319,15 @@ def main(): parser.add_argument('-f', '--force', help='force overwrite', action='store_true') args = parser.parse_args() - if os.path.exists(args.file): - with Imread(args.file, transform=args.register) as im: - print(im.summary) - if args.out: - out = os.path.abspath(args.out) - if not os.path.exists(os.path.dirname(out)): - os.makedirs(os.path.dirname(out)) - if os.path.exists(out) and not args.force: - print('File {} exists already, add the -f flag if you want to overwrite it.'.format(args.out)) - else: - im.save_as_tiff(out, args.channel, args.zslice, args.time, args.split) - else: - print('File does not exist.') + with Imread(args.file, transform=args.register) as im: + print(im.summary) + if args.out: + out = Path(args.out).absolute() + out.parent.mkdir(parents=True, exist_ok=True) + if out.exists() and not args.force: + print('File {} exists already, add the -f flag if you want to overwrite it.'.format(args.out)) + else: + im.save_as_tiff(out, args.channel, args.zslice, args.time, args.split) from ndbioimage.readers import * diff --git a/ndbioimage/_version.py b/ndbioimage/_version.py deleted file mode 100644 index e1d65be..0000000 --- a/ndbioimage/_version.py +++ /dev/null @@ -1,2 +0,0 @@ -__version__ = '2022.7.0' -__git_commit_hash__ = 'unknown' diff --git a/ndbioimage/jvm.py b/ndbioimage/jvm.py index 349a5cc..0b85bb3 100644 --- a/ndbioimage/jvm.py +++ b/ndbioimage/jvm.py @@ -1,7 +1,6 @@ -try: - import javabridge - import bioformats +from pathlib import Path +try: class JVM: """ There can be only one java virtual machine per python process, so this is a singleton class to manage the jvm. @@ -9,30 +8,43 @@ try: _instance = None vm_started = False vm_killed = False + success = True def __new__(cls, *args): if cls._instance is None: cls._instance = object.__new__(cls) return cls._instance - def start_vm(self): + def __init__(self, classpath=None): if not self.vm_started and not self.vm_killed: - javabridge.start_vm(class_path=bioformats.JARS, run_headless=True) - outputstream = javabridge.make_instance('java/io/ByteArrayOutputStream', "()V") - printstream = javabridge.make_instance('java/io/PrintStream', "(Ljava/io/OutputStream;)V", outputstream) - javabridge.static_call('Ljava/lang/System;', "setOut", "(Ljava/io/PrintStream;)V", printstream) - javabridge.static_call('Ljava/lang/System;', "setErr", "(Ljava/io/PrintStream;)V", printstream) + if classpath is None: + classpath = [str(Path(__file__).parent / 'jars' / '*')] + + import jpype + jpype.startJVM(classpath=classpath) + import jpype.imports + from loci.common import DebugTools + from loci.formats import ImageReader + from loci.formats import ChannelSeparator + from loci.formats import FormatTools + from loci.formats import MetadataTools + + DebugTools.setRootLevel("ERROR") self.vm_started = True - log4j = javabridge.JClassWrapper("loci.common.Log4jTools") - log4j.enableLogging() - log4j.setRootLevel("ERROR") + self.image_reader = ImageReader + self.channel_separator = ChannelSeparator + self.format_tools = FormatTools + self.metadata_tools = MetadataTools if self.vm_killed: raise Exception('The JVM was killed before, and cannot be restarted in this Python process.') def kill_vm(self): - javabridge.kill_vm() + if self.vm_started and not self.vm_killed: + import jpype + jpype.shutdownJVM() self.vm_started = False self.vm_killed = True + except ImportError: JVM = None diff --git a/ndbioimage/readers/__init__.py b/ndbioimage/readers/__init__.py index 7d6663a..30d03de 100644 --- a/ndbioimage/readers/__init__.py +++ b/ndbioimage/readers/__init__.py @@ -1,3 +1,2 @@ -import os -__all__ = [os.path.splitext(os.path.basename(file))[0] for file in os.listdir(os.path.dirname(__file__)) - if file.endswith('.py') and not file == os.path.basename(__file__)] +from pathlib import Path +__all__ = [file.stem for file in Path(__file__).parent.iterdir() if file.suffix == ".py" and not file == Path(__file__)] diff --git a/ndbioimage/readers/bfread.py b/ndbioimage/readers/bfread.py index 1a00a48..ada5bc3 100644 --- a/ndbioimage/readers/bfread.py +++ b/ndbioimage/readers/bfread.py @@ -1,89 +1,187 @@ -from ndbioimage import Imread, XmlData, JVM -import os +import multiprocessing import numpy as np -import untangle +from abc import ABC +from ndbioimage import Imread, JVM +from multiprocessing import queues +from traceback import print_exc +from urllib import request +from pathlib import Path -if JVM is not None: - import bioformats - class Reader(Imread): - """ This class is used as a last resort, when we don't have another way to open the file. We don't like it - because it requires the java vm. +class JVMReader: + def __init__(self, path, series): + bf_jar = Path(__file__).parent.parent / 'jars' / 'bioformats_package.jar' + if not bf_jar.exists(): + print('Downloading bioformats_package.jar.') + url = 'https://downloads.openmicroscopy.org/bio-formats/latest/artifacts/bioformats_package.jar' + bf_jar.write_bytes(request.urlopen(url).read()) + + mp = multiprocessing.get_context('spawn') + self.path = path + self.series = series + self.queue_in = mp.Queue() + self.queue_out = mp.Queue() + self.queue_error = mp.Queue() + self.done = mp.Event() + self.process = mp.Process(target=self.run) + self.process.start() + self.is_alive = True + + def close(self): + if self.is_alive: + self.done.set() + while not self.queue_in.empty(): + self.queue_in.get() + self.queue_in.close() + self.queue_in.join_thread() + while not self.queue_out.empty(): + print(self.queue_out.get()) + self.queue_out.close() + self.process.join() + self.process.close() + self.is_alive = False + + def frame(self, c, z, t): + self.queue_in.put((c, z, t)) + return self.queue_out.get() + + def run(self): + """ Read planes from the image reader file. + adapted from python-bioformats/bioformats/formatreader.py """ - priority = 99 # panic and open with BioFormats - do_not_pickle = 'reader', 'key', 'jvm' + jvm = JVM() + reader = jvm.image_reader() + ome_meta = jvm.metadata_tools.createOMEXMLMetadata() + reader.setMetadataStore(ome_meta) + reader.setId(str(self.path)) + reader.setSeries(self.series) - @staticmethod - def _can_open(path): - return True + open_bytes_func = reader.openBytes + width, height = int(reader.getSizeX()), int(reader.getSizeY()) - def open(self): - self.jvm = JVM() - self.jvm.start_vm() - self.key = np.random.randint(1e9) - self.reader = bioformats.get_image_reader(self.key, self.path) + pixel_type = reader.getPixelType() + little_endian = reader.isLittleEndian() - def __metadata__(self): - s = self.reader.rdr.getSeriesCount() - if self.series >= s: - print('Series {} does not exist.'.format(self.series)) - self.reader.rdr.setSeries(self.series) + if pixel_type == jvm.format_tools.INT8: + dtype = np.int8 + elif pixel_type == jvm.format_tools.UINT8: + dtype = np.uint8 + elif pixel_type == jvm.format_tools.UINT16: + dtype = 'u2' + elif pixel_type == jvm.format_tools.INT16: + dtype = 'i2' + elif pixel_type == jvm.format_tools.UINT32: + dtype = 'u4' + elif pixel_type == jvm.format_tools.INT32: + dtype = 'i4' + elif pixel_type == jvm.format_tools.FLOAT: + dtype = 'f4' + elif pixel_type == jvm.format_tools.DOUBLE: + dtype = 'f8' + else: + dtype = None - X = self.reader.rdr.getSizeX() - Y = self.reader.rdr.getSizeY() - C = self.reader.rdr.getSizeC() - Z = self.reader.rdr.getSizeZ() - T = self.reader.rdr.getSizeT() - self.shape = (X, Y, C, Z, T) + try: + while not self.done.is_set(): + try: + c, z, t = self.queue_in.get(True, 0.02) + if reader.isRGB() and reader.isInterleaved(): + index = reader.getIndex(z, 0, t) + image = np.frombuffer(open_bytes_func(index), dtype) + image.shape = (height, width, reader.getSizeC()) + if image.shape[2] > 3: + image = image[:, :, :3] + elif c is not None and reader.getRGBChannelCount() == 1: + index = reader.getIndex(z, c, t) + image = np.frombuffer(open_bytes_func(index), dtype) + image.shape = (height, width) + elif reader.getRGBChannelCount() > 1: + n_planes = reader.getRGBChannelCount() + rdr = jvm.channel_separator(reader) + planes = [np.frombuffer(rdr.openBytes(rdr.getIndex(z, i, t)), dtype) for i in range(n_planes)] + if len(planes) > 3: + planes = planes[:3] + elif len(planes) < 3: + # > 1 and < 3 means must be 2 + # see issue #775 + planes.append(np.zeros(planes[0].shape, planes[0].dtype)) + image = np.dstack(planes) + image.shape = (height, width, 3) + del rdr + elif reader.getSizeC() > 1: + images = [np.frombuffer(open_bytes_func(reader.getIndex(z, i, t)), dtype) + for i in range(reader.getSizeC())] + image = np.dstack(images) + image.shape = (height, width, reader.getSizeC()) + # if not channel_names is None: + # metadata = MetadataRetrieve(self.metadata) + # for i in range(self.reader.getSizeC()): + # index = self.reader.getIndex(z, 0, t) + # channel_name = metadata.getChannelName(index, i) + # if channel_name is None: + # channel_name = metadata.getChannelID(index, i) + # channel_names.append(channel_name) + elif reader.isIndexed(): + # + # The image data is indexes into a color lookup-table + # But sometimes the table is the identity table and just generates + # a monochrome RGB image + # + index = reader.getIndex(z, 0, t) + image = np.frombuffer(open_bytes_func(index), dtype) + if pixel_type in (jvm.format_tools.INT16, jvm.format_tools.UINT16): + lut = reader.get16BitLookupTable() + if lut is not None: + lut = np.array(lut) + # lut = np.array( + # [env.get_short_array_elements(d) + # for d in env.get_object_array_elements(lut)]) \ + # .transpose() + else: + lut = reader.get8BitLookupTable() + if lut is not None: + lut = np.array(lut) + # lut = np.array( + # [env.get_byte_array_elements(d) + # for d in env.get_object_array_elements(lut)]) \ + # .transpose() + image.shape = (height, width) + if (lut is not None) and not np.all(lut == np.arange(lut.shape[0])[:, np.newaxis]): + image = lut[image, :] + else: + index = reader.getIndex(z, 0, t) + image = np.frombuffer(open_bytes_func(index), dtype) + image.shape = (height, width) - omexml = bioformats.get_omexml_metadata(self.path) - self.metadata = XmlData(untangle.parse(omexml)) + if image.ndim == 3: + self.queue_out.put(image[..., c]) + else: + self.queue_out.put(image) + except queues.Empty: + continue + except (Exception,): + print_exc() + self.queue_out.put(np.zeros((32, 32))) + finally: + jvm.kill_vm() - image = list(self.metadata.search_all('Image').values()) - if len(image) and self.series in image[0]: - image = XmlData(image[0][self.series]) - else: - image = self.metadata - unit = lambda u: 10 ** {'nm': 9, 'µm': 6, 'um': 6, 'mm': 3, 'm': 0}[u] +class Reader(Imread, ABC): + """ This class is used as a last resort, when we don't have another way to open the file. We don't like it + because it requires the java vm. + """ + priority = 99 # panic and open with BioFormats + do_not_pickle = 'reader', 'key', 'jvm' - pxsizeunit = image.search('PhysicalSizeXUnit')[0] - pxsize = image.search('PhysicalSizeX')[0] - if pxsize is not None: - self.pxsize = pxsize / unit(pxsizeunit) * 1e6 + @staticmethod + def _can_open(path): + return not path.is_dir() - if self.zstack: - deltazunit = image.search('PhysicalSizeZUnit')[0] - deltaz = image.search('PhysicalSizeZ')[0] - if deltaz is not None: - self.deltaz = deltaz / unit(deltazunit) * 1e6 + def open(self): + self.reader = JVMReader(self.path, self.series) - if self.path.endswith('.lif'): - self.title = os.path.splitext(os.path.basename(self.path))[0] - self.exposuretime = self.metadata.re_search(r'WideFieldChannelInfo\|ExposureTime', self.exposuretime) - if self.timeseries: - self.settimeinterval = \ - self.metadata.re_search(r'ATLCameraSettingDefinition\|CycleTime', self.settimeinterval * 1e3)[ - 0] / 1000 - if not self.settimeinterval: - self.settimeinterval = self.exposuretime[0] - self.pxsizecam = self.metadata.re_search(r'ATLCameraSettingDefinition\|TheoCamSensorPixelSizeX', - self.pxsizecam) - self.objective = self.metadata.re_search(r'ATLCameraSettingDefinition\|ObjectiveName', 'none')[0] - self.magnification = \ - self.metadata.re_search(r'ATLCameraSettingDefinition\|Magnification', self.magnification)[0] - elif self.path.endswith('.ims'): - self.magnification = self.metadata.search('LensPower', 100)[0] - self.NA = self.metadata.search('NumericalAperture', 1.47)[0] - self.title = self.metadata.search('Name', self.title) - self.binning = self.metadata.search('BinningX', 1)[0] + def __frame__(self, c, z, t): + return self.reader.frame(c, z, t) - def __frame__(self, *args): - frame = self.reader.read(*args, rescale=False).astype('float') - if frame.ndim == 3: - return frame[..., args[0]] - else: - return frame - - def close(self): - bioformats.release_image_reader(self.key) + def close(self): + self.reader.close() diff --git a/ndbioimage/readers/cziread.py b/ndbioimage/readers/cziread.py index 86e5a8b..fe1c074 100644 --- a/ndbioimage/readers/cziread.py +++ b/ndbioimage/readers/cziread.py @@ -1,115 +1,432 @@ -from ndbioimage import Imread, XmlData, tolist import czifile -import untangle import numpy as np import re +from lxml import etree +from ome_types import model +from abc import ABC +from ndbioimage import Imread from functools import cached_property +from itertools import product +from pathlib import Path -class Reader(Imread): +class Reader(Imread, ABC): priority = 0 - do_not_pickle = 'reader', 'filedict', 'extrametadata' + do_not_pickle = 'reader', 'filedict' @staticmethod def _can_open(path): - return isinstance(path, str) and path.endswith('.czi') + return isinstance(path, Path) and path.suffix == '.czi' def open(self): self.reader = czifile.CziFile(self.path) filedict = {} for directory_entry in self.reader.filtered_subblock_directory: idx = self.get_index(directory_entry, self.reader.start) - for c in range(*idx[self.reader.axes.index('C')]): - for z in range(*idx[self.reader.axes.index('Z')]): - for t in range(*idx[self.reader.axes.index('T')]): - if (c, z, t) in filedict: - filedict[(c, z, t)].append(directory_entry) - else: - filedict[(c, z, t)] = [directory_entry] + if 'S' not in self.reader.axes or self.series in range(*idx[self.reader.axes.index('S')]): + for c in range(*idx[self.reader.axes.index('C')]): + for z in range(*idx[self.reader.axes.index('Z')]): + for t in range(*idx[self.reader.axes.index('T')]): + if (c, z, t) in filedict: + filedict[c, z, t].append(directory_entry) + else: + filedict[c, z, t] = [directory_entry] self.filedict = filedict def close(self): self.reader.close() - def __metadata__(self): - # TODO: make sure frame function still works when a subblock has data from more than one frame - self.shape = tuple([self.reader.shape[self.reader.axes.index(directory_entry)] for directory_entry in 'XYCZT']) - self.metadata = XmlData(untangle.parse(self.reader.metadata())) - - image = [i for i in self.metadata.search_all('Image').values() if i] - if len(image) and self.series in image[0]: - image = XmlData(image[0][self.series]) + @cached_property + def ome(self): + xml = self.reader.metadata() + attachments = {i.attachment_entry.name: i.attachment_entry.data_segment() + for i in self.reader.attachments()} + tree = etree.fromstring(xml) + metadata = tree.find("Metadata") + version = metadata.find("Version") + if version is not None: + version = version.text else: - image = self.metadata + version = metadata.find("Experiment").attrib["Version"] - pxsize = image.search('ScalingX')[0] - if pxsize is not None: - self.pxsize = pxsize * 1e6 - if self.zstack: - deltaz = image.search('ScalingZ')[0] - if deltaz is not None: - self.deltaz = deltaz * 1e6 + if version == '1.0': + return self.ome_10(tree, attachments) + elif version == '1.2': + return self.ome_12(tree, attachments) - self.title = self.metadata.re_search(('Information', 'Document', 'Name'), self.title)[0] - self.acquisitiondate = self.metadata.re_search(('Information', 'Document', 'CreationDate'), - self.acquisitiondate)[0] - self.exposuretime = self.metadata.re_search(('TrackSetup', 'CameraIntegrationTime'), self.exposuretime) - if self.timeseries: - self.settimeinterval = self.metadata.re_search(('Interval', 'TimeSpan', 'Value'), - self.settimeinterval * 1e3)[0] / 1000 - if not self.settimeinterval: - self.settimeinterval = self.exposuretime[0] - self.pxsizecam = self.metadata.re_search(('AcquisitionModeSetup', 'PixelPeriod'), self.pxsizecam) - self.magnification = self.metadata.re_search('NominalMagnification', self.magnification)[0] - attenuators = self.metadata.search_all('Attenuator') - self.laserwavelengths = [[1e9 * float(i['Wavelength']) for i in tolist(attenuator)] - for attenuator in attenuators.values()] - self.laserpowers = [[float(i['Transmission']) for i in tolist(attenuator)] - for attenuator in attenuators.values()] - self.collimator = self.metadata.re_search(('Collimator', 'Position')) - detector = self.metadata.search(('Instrument', 'Detector')) - self.gain = [int(i.get('AmplificationGain', 1)) for i in detector] - self.powermode = self.metadata.re_search(('TrackSetup', 'FWFOVPosition'))[0] - optovar = self.metadata.re_search(('TrackSetup', 'TubeLensPosition'), '1x') - self.optovar = [] - for o in optovar: - a = re.search(r'\d?\d*[,.]?\d+(?=x$)', o) - if hasattr(a, 'group'): - self.optovar.append(float(a.group(0).replace(',', '.'))) - self.pcf = [2 ** self.metadata.re_search(('Image', 'ComponentBitCount'), 14)[0] / float(i) - for i in self.metadata.re_search(('Channel', 'PhotonConversionFactor'), 1)] - self.binning = self.metadata.re_search(('AcquisitionModeSetup', 'CameraBinning'), 1)[0] - self.objective = self.metadata.re_search(('AcquisitionModeSetup', 'Objective'))[0] - self.NA = self.metadata.re_search(('Instrument', 'Objective', 'LensNA'))[0] - self.filter = self.metadata.re_search(('TrackSetup', 'BeamSplitter', 'Filter'))[0] - self.tirfangle = [50 * i for i in self.metadata.re_search(('TrackSetup', 'TirfAngle'), 0)] - self.frameoffset = [self.metadata.re_search(('AcquisitionModeSetup', 'CameraFrameOffsetX'))[0], - self.metadata.re_search(('AcquisitionModeSetup', 'CameraFrameOffsetY'))[0]] - self.cnamelist = [c['DetectorSettings']['Detector']['Id'] for c in - self.metadata['ImageDocument']['Metadata']['Information']['Image'].search('Channel')] - try: - self.track, self.detector = zip(*[[int(i) for i in re.findall(r'\d', c)] for c in self.cnamelist]) - except ValueError: - self.track = tuple(range(len(self.cnamelist))) - self.detector = (0,) * len(self.cnamelist) + def ome_12(self, tree, attachments): + def text(item, default=""): + return default if item is None else item.text + + ome = model.OME() + + metadata = tree.find("Metadata") + + information = metadata.find("Information") + display_setting = metadata.find("DisplaySetting") + ome.experimenters = [model.Experimenter(id="Experimenter:0", + user_name=information.find("Document").find("UserName").text)] + + instrument = information.find("Instrument") + for _ in instrument.find("Microscopes"): + ome.instruments.append(model.Instrument()) + + for detector in instrument.find("Detectors"): + try: + detector_type = model.detector.Type(text(detector.find("Type")).upper() or "") + except ValueError: + detector_type = model.detector.Type.OTHER + + ome.instruments[0].detectors.append( + model.Detector( + id=detector.attrib["Id"].replace(' ', ''), model=text(detector.find("Manufacturer").find("Model")), + type=detector_type + )) + + for objective in instrument.find("Objectives"): + ome.instruments[0].objectives.append( + model.Objective( + id=objective.attrib["Id"], + model=text(objective.find("Manufacturer").find("Model")), + immersion=text(objective.find("Immersion")), + lens_na=float(text(objective.find("LensNA"))), + nominal_magnification=float(text(objective.find("NominalMagnification"))))) + + for tubelens in instrument.find("TubeLenses"): + ome.instruments[0].objectives.append( + model.Objective( + id=f'Objective:{tubelens.attrib["Id"]}', + model=tubelens.attrib["Name"], + nominal_magnification=1.0)) # TODO: nominal_magnification + + for light_source in instrument.find("LightSources"): + if light_source.find("LightSourceType").find("Laser") is not None: + ome.instruments[0].light_source_group.append( + model.Laser( + id=f'LightSource:{light_source.attrib["Id"]}', + power=float(text(light_source.find("Power"))), + wavelength=float(light_source.attrib["Id"][-3:]))) + + x_min = min([f.start[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_min = min([f.start[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + x_max = max([f.start[f.axes.index('X')] + f.shape[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_max = max([f.start[f.axes.index('Y')] + f.shape[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + size_x = x_max - x_min + size_y = y_max - y_min + size_c, size_z, size_t = [self.reader.shape[self.reader.axes.index(directory_entry)] + for directory_entry in 'CZT'] + + image = information.find("Image") + pixel_type = text(image.find("PixelType"), "Gray16") + if pixel_type.startswith("Gray"): + pixel_type = "uint" + pixel_type[4:] + objective_settings = image.find("ObjectiveSettings") + scenes = image.find("Dimensions").find("S").find("Scenes") + center_position = [float(pos) for pos in text(scenes[0].find("CenterPosition")).split(',')] + um = model.simple_types.UnitsLength.MICROMETER + nm = model.simple_types.UnitsLength.NANOMETER + + ome.images.append( + model.Image( + id="Image:0", + name=f'{text(information.find("Document").find("Name"))} #1', + pixels=model.Pixels( + id="Pixels:0", size_x=size_x, size_y=size_y, + size_c=size_c, size_z=size_z, size_t=size_t, + dimension_order="XYCZT", type=pixel_type, + significant_bits=int(text(image.find("ComponentBitCount"))), + big_endian=False, interleaved=False, metadata_only=True), + experimenter_ref=model.ExperimenterRef(id='Experimenter:0'), + instrument_ref=model.InstrumentRef(id='Instrument:0'), + objective_settings=model.ObjectiveSettings( + id=objective_settings.find("ObjectiveRef").attrib["Id"], + medium=text(objective_settings.find("Medium")), + refractive_index=float(text(objective_settings.find("RefractiveIndex")))), + stage_label=model.StageLabel( + name=f"Scene position #0", + x=center_position[0], x_unit=um, + y=center_position[1], y_unit=um))) + + for distance in metadata.find("Scaling").find("Items"): + if distance.attrib["Id"] == "X": + ome.images[0].pixels.physical_size_x = float(text(distance.find("Value"))) * 1e6 + elif distance.attrib["Id"] == "Y": + ome.images[0].pixels.physical_size_y = float(text(distance.find("Value"))) * 1e6 + elif size_z > 1 and distance.attrib["Id"] == "Z": + ome.images[0].pixels.physical_size_z = float(text(distance.find("Value"))) * 1e6 + + channels_im = {channel.attrib["Id"]: channel for channel in image.find("Dimensions").find("Channels")} + channels_ds = {channel.attrib["Id"]: channel for channel in display_setting.find("Channels")} + + for idx, (key, channel) in enumerate(channels_im.items()): + detector_settings = channel.find("DetectorSettings") + laser_scan_info = channel.find("LaserScanInfo") + detector = detector_settings.find("Detector") + try: + binning = model.simple_types.Binning(text(detector_settings.find("Binning"))) + except ValueError: + binning = model.simple_types.Binning.OTHER + + light_sources_settings = channel.find("LightSourcesSettings") + # no space in ome for multiple lightsources simultaneously + light_source_settings = light_sources_settings[0] + light_source_settings = model.LightSourceSettings( + id="LightSource:" + "_".join([light_source_settings.find("LightSource").attrib["Id"] + for light_source_settings in light_sources_settings]), + attenuation=float(text(light_source_settings.find("Attenuation"))), + wavelength=float(text(light_source_settings.find("Wavelength"))), + wavelength_unit=nm) + + ome.images[0].pixels.channels.append( + model.Channel( + id=f"Channel:0:{idx}", + name=channel.attrib["Name"], + acquisition_mode=text(channel.find("AcquisitionMode")), + color=model.simple_types.Color(text(channels_ds[channel.attrib["Id"]].find("Color"))), + detector_settings=model.DetectorSettings( + id=detector.attrib["Id"].replace(" ", ""), + binning=binning), + emission_wavelength=text(channel.find("EmissionWavelength")), + excitation_wavelength=text(channel.find("ExcitationWavelength")), + # filter_set_ref=model.FilterSetRef(id=ome.instruments[0].filter_sets[filterset_idx].id), + illumination_type=text(channel.find("IlluminationType")), + light_source_settings=light_source_settings, + samples_per_pixel=int(text(laser_scan_info.find("Averaging"))))) + + exposure_times = [float(text(channel.find("LaserScanInfo").find("FrameTime"))) for channel in + channels_im.values()] + delta_ts = attachments['TimeStamps'].data() + for t, z, c in product(range(size_t), range(size_z), range(size_c)): + ome.images[0].pixels.planes.append( + model.Plane(the_c=c, the_z=z, the_t=t, delta_t=delta_ts[t], exposure_time=exposure_times[c])) + + idx = 0 + for layer in metadata.find("Layers"): + rectangle = layer.find("Elements").find("Rectangle") + if rectangle is not None: + geometry = rectangle.find("Geometry") + roi = model.ROI(id=f"ROI:{idx}", description=text(layer.find("Usage"))) + roi.union.append( + model.Rectangle( + id='Shape:0:0', + height=float(text(geometry.find("Height"))), + width=float(text(geometry.find("Width"))), + x=float(text(geometry.find("Left"))), + y=float(text(geometry.find("Top"))))) + ome.rois.append(roi) + ome.images[0].roi_ref.append(model.ROIRef(id=f"ROI:{idx}")) + idx += 1 + return ome + + def ome_10(self, tree, attachments): + def text(item, default=""): + return default if item is None else item.text + + ome = model.OME() + + metadata = tree.find("Metadata") + + information = metadata.find("Information") + display_setting = metadata.find("DisplaySetting") + experiment = metadata.find("Experiment") + acquisition_block = experiment.find("ExperimentBlocks").find("AcquisitionBlock") + + ome.experimenters = [model.Experimenter(id="Experimenter:0", + user_name=information.find("User").find("DisplayName").text)] + + instrument = information.find("Instrument") + ome.instruments.append(model.Instrument(id=instrument.attrib["Id"])) + + for detector in instrument.find("Detectors"): + try: + detector_type = model.detector.Type(text(detector.find("Type")).upper() or "") + except ValueError: + detector_type = model.detector.Type.OTHER + + ome.instruments[0].detectors.append( + model.Detector( + id=detector.attrib["Id"], model=text(detector.find("Manufacturer").find("Model")), + amplification_gain=float(text(detector.find("AmplificationGain"))), + gain=float(text(detector.find("Gain"))), zoom=float(text(detector.find("Zoom"))), + type=detector_type + )) + + for objective in instrument.find("Objectives"): + ome.instruments[0].objectives.append( + model.Objective( + id=objective.attrib["Id"], + model=text(objective.find("Manufacturer").find("Model")), + immersion=text(objective.find("Immersion")), + lens_na=float(text(objective.find("LensNA"))), + nominal_magnification=float(text(objective.find("NominalMagnification"))))) + + for light_source in instrument.find("LightSources"): + if light_source.find("LightSourceType").find("Laser") is not None: + ome.instruments[0].light_source_group.append( + model.Laser( + id=light_source.attrib["Id"], + model=text(light_source.find("Manufacturer").find("Model")), + power=float(text(light_source.find("Power"))), + wavelength=float( + text(light_source.find("LightSourceType").find("Laser").find("Wavelength"))))) + + multi_track_setup = acquisition_block.find("MultiTrackSetup") + for idx, tube_lens in enumerate(set(text(track_setup.find("TubeLensPosition")) + for track_setup in multi_track_setup)): + ome.instruments[0].objectives.append( + model.Objective(id=f"Objective:Tubelens:{idx}", model=tube_lens, + nominal_magnification=float( + re.findall(r'\d+[,.]\d*', tube_lens)[0].replace(',', '.')) + )) + + for idx, filter_ in enumerate(set(text(beam_splitter.find("Filter")) + for track_setup in multi_track_setup + for beam_splitter in track_setup.find("BeamSplitters"))): + ome.instruments[0].filter_sets.append( + model.FilterSet(id=f"FilterSet:{idx}", model=filter_) + ) + + for idx, collimator in enumerate(set(text(track_setup.find("FWFOVPosition")) + for track_setup in multi_track_setup)): + ome.instruments[0].filters.append(model.Filter(id=f"Filter:Collimator:{idx}", model=collimator)) + + x_min = min([f.start[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_min = min([f.start[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + x_max = max([f.start[f.axes.index('X')] + f.shape[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_max = max([f.start[f.axes.index('Y')] + f.shape[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + size_x = x_max - x_min + size_y = y_max - y_min + size_c, size_z, size_t = [self.reader.shape[self.reader.axes.index(directory_entry)] + for directory_entry in 'CZT'] + + image = information.find("Image") + pixel_type = text(image.find("PixelType"), "Gray16") + if pixel_type.startswith("Gray"): + pixel_type = "uint" + pixel_type[4:] + objective_settings = image.find("ObjectiveSettings") + scenes = image.find("Dimensions").find("S").find("Scenes") + positions = scenes[0].find("Positions")[0] + um = model.simple_types.UnitsLength.MICROMETER + nm = model.simple_types.UnitsLength.NANOMETER + + ome.images.append( + model.Image( + id="Image:0", + name=f'{text(information.find("Document").find("Name"))} #1', + pixels=model.Pixels( + id="Pixels:0", size_x=size_x, size_y=size_y, + size_c=size_c, size_z=size_z, size_t=size_t, + dimension_order="XYCZT", type=pixel_type, + significant_bits=int(text(image.find("ComponentBitCount"))), + big_endian=False, interleaved=False, metadata_only=True), + experimenter_ref=model.ExperimenterRef(id='Experimenter:0'), + instrument_ref=model.InstrumentRef(id='Instrument:0'), + objective_settings=model.ObjectiveSettings( + id=objective_settings.find("ObjectiveRef").attrib["Id"], + medium=text(objective_settings.find("Medium")), + refractive_index=float(text(objective_settings.find("RefractiveIndex")))), + stage_label=model.StageLabel( + name=f"Scene position #0", + x=float(positions.attrib["X"]), x_unit=um, + y=float(positions.attrib["Y"]), y_unit=um, + z=float(positions.attrib["Z"]), z_unit=um))) + + for distance in metadata.find("Scaling").find("Items"): + if distance.attrib["Id"] == "X": + ome.images[0].pixels.physical_size_x = float(text(distance.find("Value"))) * 1e6 + elif distance.attrib["Id"] == "Y": + ome.images[0].pixels.physical_size_y = float(text(distance.find("Value"))) * 1e6 + elif size_z > 1 and distance.attrib["Id"] == "Z": + ome.images[0].pixels.physical_size_z = float(text(distance.find("Value"))) * 1e6 + + channels_im = {channel.attrib["Id"]: channel for channel in image.find("Dimensions").find("Channels")} + channels_ds = {channel.attrib["Id"]: channel for channel in display_setting.find("Channels")} + channels_ts = {detector.attrib["Id"]: track_setup + for track_setup in + experiment.find("ExperimentBlocks").find("AcquisitionBlock").find("MultiTrackSetup") + for detector in track_setup.find("Detectors")} + + for idx, (key, channel) in enumerate(channels_im.items()): + detector_settings = channel.find("DetectorSettings") + laser_scan_info = channel.find("LaserScanInfo") + detector = detector_settings.find("Detector") + try: + binning = model.simple_types.Binning(text(detector_settings.find("Binning"))) + except ValueError: + binning = model.simple_types.Binning.OTHER + + filterset = text(channels_ts[key].find("BeamSplitters")[0].find("Filter")) + filterset_idx = [filterset.model for filterset in ome.instruments[0].filter_sets].index(filterset) + + light_sources_settings = channel.find("LightSourcesSettings") + # no space in ome for multiple lightsources simultaneously + light_source_settings = light_sources_settings[0] + light_source_settings = model.LightSourceSettings( + id="_".join([light_source_settings.find("LightSource").attrib["Id"] + for light_source_settings in light_sources_settings]), + attenuation=float(text(light_source_settings.find("Attenuation"))), + wavelength=float(text(light_source_settings.find("Wavelength"))), + wavelength_unit=nm) + + ome.images[0].pixels.channels.append( + model.Channel( + id=f"Channel:0:{idx}", + name=channel.attrib["Name"], + acquisition_mode=text(channel.find("AcquisitionMode")), + color=model.simple_types.Color(text(channels_ds[channel.attrib["Id"]].find("Color"))), + detector_settings=model.DetectorSettings(id=detector.attrib["Id"], binning=binning), + emission_wavelength=text(channel.find("EmissionWavelength")), + excitation_wavelength=text(channel.find("ExcitationWavelength")), + filter_set_ref=model.FilterSetRef(id=ome.instruments[0].filter_sets[filterset_idx].id), + illumination_type=text(channel.find("IlluminationType")), + light_source_settings=light_source_settings, + samples_per_pixel=int(text(laser_scan_info.find("Averaging"))))) + + exposure_times = [float(text(channel.find("LaserScanInfo").find("FrameTime"))) for channel in + channels_im.values()] + delta_ts = attachments['TimeStamps'].data() + for t, z, c in product(range(size_t), range(size_z), range(size_c)): + ome.images[0].pixels.planes.append( + model.Plane(the_c=c, the_z=z, the_t=t, delta_t=delta_ts[t], + exposure_time=exposure_times[c], + position_x=float(positions.attrib["X"]), position_x_unit=um, + position_y=float(positions.attrib["Y"]), position_y_unit=um, + position_z=float(positions.attrib["Z"]), position_z_unit=um)) + + idx = 0 + for layer in metadata.find("Layers"): + rectangle = layer.find("Elements").find("Rectangle") + if rectangle is not None: + geometry = rectangle.find("Geometry") + roi = model.ROI(id=f"ROI:{idx}", description=text(layer.find("Usage"))) + roi.union.append( + model.Rectangle( + id='Shape:0:0', + height=float(text(geometry.find("Height"))), + width=float(text(geometry.find("Width"))), + x=float(text(geometry.find("Left"))), + y=float(text(geometry.find("Top"))))) + ome.rois.append(roi) + ome.images[0].roi_ref.append(model.ROIRef(id=f"ROI:{idx}")) + idx += 1 + return ome def __frame__(self, c=0, z=0, t=0): - f = np.zeros(self.file_shape[:2], self.dtype) - for directory_entry in self.filedict[(c, z, t)]: + f = np.zeros(self.shape['xy'], self.dtype) + directory_entries = self.filedict[c, z, t] + x_min = min([f.start[f.axes.index('X')] for f in directory_entries]) + y_min = min([f.start[f.axes.index('Y')] for f in directory_entries]) + xy_min = {'X': x_min, 'Y': y_min} + for directory_entry in directory_entries: subblock = directory_entry.data_segment() tile = subblock.data(resize=True, order=0) - index = [slice(i - j, i - j + k) for i, j, k in - zip(directory_entry.start, self.reader.start, tile.shape)] - index = tuple([index[self.reader.axes.index(i)] for i in 'XY']) + axes_min = [xy_min.get(ax, 0) for ax in directory_entry.axes] + index = [slice(i - j - m, i - j + k) + for i, j, k, m in zip(directory_entry.start, self.reader.start, tile.shape, axes_min)] + index = tuple(index[self.reader.axes.index(i)] for i in 'XY') f[index] = tile.squeeze() return f @staticmethod def get_index(directory_entry, start): return [(i - j, i - j + k) for i, j, k in zip(directory_entry.start, start, directory_entry.shape)] - - @cached_property - def timeval(self): - tval = np.unique(list(filter(lambda x: x.attachment_entry.filename.startswith('TimeStamp'), - self.reader.attachments()))[0].data()) - return sorted(tval[tval > 0])[:self.shape['t']] \ No newline at end of file diff --git a/ndbioimage/readers/ndread.py b/ndbioimage/readers/ndread.py index ab4d46b..f73db50 100644 --- a/ndbioimage/readers/ndread.py +++ b/ndbioimage/readers/ndread.py @@ -1,29 +1,55 @@ -from ndbioimage import Imread import numpy as np +from ome_types import model +from functools import cached_property +from abc import ABC +from ndbioimage import Imread +from itertools import product -class Reader(Imread): +class Reader(Imread, ABC): priority = 20 + do_not_pickle = 'ome' + do_not_copy = 'ome' @staticmethod def _can_open(path): return isinstance(path, np.ndarray) and 1 <= path.ndim <= 5 - def __metadata__(self): - self.base = np.array(self.path, ndmin=5) - self.title = self.path = 'numpy array' - self.axes = self.axes[:self.base.ndim] - self.shape = self.base.shape - self.acquisitiondate = 'now' + @cached_property + def ome(self): + def shape(size_x=1, size_y=1, size_c=1, size_z=1, size_t=1): + return size_x, size_y, size_c, size_z, size_t + size_x, size_y, size_c, size_z, size_t = shape(*self.array.shape) + try: + pixel_type = model.simple_types.PixelType(self.array.dtype.name) + except ValueError: + if self.array.dtype.name.startswith('int'): + pixel_type = model.simple_types.PixelType('int32') + else: + pixel_type = model.simple_types.PixelType('float') + + ome = model.OME() + ome.instruments.append(model.Instrument()) + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, size_x=size_x, size_y=size_y, + dimension_order="XYCZT", type=pixel_type, + objective_settings=model.ObjectiveSettings(id="Objective:0")))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append(model.Plane(the_c=c, the_z=z, the_t=t, delta_t=0)) + return ome + + def open(self): + self.array = np.array(self.path) + self.path = 'numpy array' def __frame__(self, c, z, t): - xyczt = (slice(None), slice(None), c, z, t) - in_idx = tuple(xyczt['xyczt'.find(i)] for i in self.axes) - frame = self.base[in_idx] + # xyczt = (slice(None), slice(None), c, z, t) + # in_idx = tuple(xyczt['xyczt'.find(i)] for i in self.axes) + # print(f'{in_idx = }') + frame = self.array[:, :, c, z, t] if self.axes.find('y') < self.axes.find('x'): return frame.T else: return frame - - def __str__(self): - return self.path diff --git a/ndbioimage/readers/seqread.py b/ndbioimage/readers/seqread.py index 63296fc..a40692e 100644 --- a/ndbioimage/readers/seqread.py +++ b/ndbioimage/readers/seqread.py @@ -1,85 +1,110 @@ -from ndbioimage import Imread, XmlData -import os +from abc import ABC + import tifffile import yaml -import json import re +from ndbioimage import Imread +from pathlib import Path +from functools import cached_property +from ome_types import model +from itertools import product +from datetime import datetime -class Reader(Imread): +class Reader(Imread, ABC): priority = 10 @staticmethod def _can_open(path): - return isinstance(path, str) and os.path.splitext(path)[1] == '' + return isinstance(path, Path) and path.suffix == "" - def __metadata__(self): - filelist = sorted([file for file in os.listdir(self.path) if re.search(r'^img_\d{3,}.*\d{3,}.*\.tif$', file)]) + def get_metadata(self, c, z, t): + with tifffile.TiffFile(self.filedict[c, z, t]) as tif: + return {key: yaml.safe_load(value) for key, value in tif.pages[0].tags[50839].value.items()} - try: - with tifffile.TiffFile(os.path.join(self.path, filelist[0])) as tif: - self.metadata = XmlData({key: yaml.safe_load(value) - for key, value in tif.pages[0].tags[50839].value.items()}) - except Exception: # fallback - with open(os.path.join(self.path, 'metadata.txt'), 'r') as metadatafile: - self.metadata = XmlData(json.loads(metadatafile.read())) + @cached_property + def ome(self): + ome = model.OME() + metadata = self.get_metadata(0, 0, 0) + ome.experimenters.append( + model.Experimenter(id="Experimenter:0", user_name=metadata["Info"]["Summary"]["UserName"])) + objective_str = metadata["Info"]["ZeissObjectiveTurret-Label"] + ome.instruments.append(model.Instrument()) + ome.instruments[0].objectives.append( + model.Objective( + id="Objective:0", manufacturer="Zeiss", model=objective_str, + nominal_magnification=float(re.findall(r"(\d+)x", objective_str)[0]), + lens_na=float(re.findall(r"/(\d\.\d+)", objective_str)[0]), + immersion=model.objective.Immersion.OIL if 'oil' in objective_str.lower() else None)) + tubelens_str = metadata["Info"]["ZeissOptovar-Label"] + ome.instruments[0].objectives.append( + model.Objective( + id="Objective:Tubelens:0", manufacturer="Zeiss", model=tubelens_str, + nominal_magnification=float(re.findall(r"\d?\d*[,.]?\d+(?=x$)", tubelens_str)[0].replace(",", ".")))) + ome.instruments[0].detectors.append( + model.Detector( + id="Detector:0", amplification_gain=100)) + ome.instruments[0].filter_sets.append( + model.FilterSet(id='FilterSet:0', model=metadata["Info"]["ZeissReflectorTurret-Label"])) + + pxsize = metadata["Info"]["PixelSizeUm"] + pxsize_cam = 6.5 if 'Hamamatsu' in metadata["Info"]["Core-Camera"] else None + if pxsize == 0: + pxsize = pxsize_cam / ome.instruments[0].objectives[0].nominal_magnification + pixel_type = metadata["Info"]["PixelType"].lower() + if pixel_type.startswith("gray"): + pixel_type = "uint" + pixel_type[4:] + else: + pixel_type = "uint16" # assume + + size_c, size_z, size_t = [max(i) + 1 for i in zip(*self.filedict.keys())] + t0 = datetime.strptime(metadata["Info"]["Time"], "%Y-%m-%d %H:%M:%S %z") + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, + size_x=metadata['Info']['Width'], size_y=metadata['Info']['Height'], + dimension_order="XYCZT", type=pixel_type, physical_size_x=pxsize, physical_size_y=pxsize, + physical_size_z=metadata["Info"]["Summary"]["z-step_um"]), + objective_settings=model.ObjectiveSettings(id="Objective:0"))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append( + model.Plane( + the_c=c, the_z=z, the_t=t, exposure_time=metadata["Info"]["Exposure-ms"] / 1000, + delta_t=(datetime.strptime(self.get_metadata(c, z, t)["Info"]["Time"], + "%Y-%m-%d %H:%M:%S %z") - t0).seconds)) # compare channel names from metadata with filenames - cnamelist = self.metadata.search('ChNames') - cnamelist = [c for c in cnamelist if any([c in f for f in filelist])] + pattern_c = re.compile(r"img_\d{3,}_(.*)_\d{3,}$") + for c in range(size_c): + ome.images[0].pixels.channels.append( + model.Channel( + id=f"Channel:{c}", name=pattern_c.findall(self.filedict[c, 0, 0].stem)[0], + detector_settings=model.DetectorSettings( + id="Detector:0", binning=metadata["Info"]["Hamamatsu_sCMOS-Binning"]), + filter_set_ref=model.FilterSetRef(id='FilterSet:0'))) + return ome - self.filedict = {} - maxc = 0 - maxz = 0 - maxt = 0 - for file in filelist: - T = re.search(r'(?<=img_)\d{3,}', file) - Z = re.search(r'\d{3,}(?=\.tif$)', file) - C = file[T.end() + 1:Z.start() - 1] - t = int(T.group(0)) - z = int(Z.group(0)) - if C in cnamelist: - c = cnamelist.index(C) - else: - c = len(cnamelist) - cnamelist.append(C) - - self.filedict[(c, z, t)] = file - if c > maxc: - maxc = c - if z > maxz: - maxz = z - if t > maxt: - maxt = t - self.cnamelist = [str(cname) for cname in cnamelist] - - X = self.metadata.search('Width')[0] - Y = self.metadata.search('Height')[0] - self.shape = (int(X), int(Y), maxc + 1, maxz + 1, maxt + 1) - - self.pxsize = self.metadata.re_search(r'(?i)pixelsize_?um', 0)[0] - if self.zstack: - self.deltaz = self.metadata.re_search(r'(?i)z-step_?um', 0)[0] - if self.timeseries: - self.settimeinterval = self.metadata.re_search(r'(?i)interval_?ms', 0)[0] / 1000 - if 'Hamamatsu' in self.metadata.search('Core-Camera', '')[0]: - self.pxsizecam = 6.5 - self.title = self.metadata.search('Prefix')[0] - self.acquisitiondate = self.metadata.search('Time')[0] - self.exposuretime = [i / 1000 for i in self.metadata.search('Exposure-ms')] - self.objective = self.metadata.search('ZeissObjectiveTurret-Label')[0] - self.optovar = [] - for o in self.metadata.search('ZeissOptovar-Label'): - a = re.search(r'\d?\d*[,.]?\d+(?=x$)', o) - if hasattr(a, 'group'): - self.optovar.append(float(a.group(0).replace(',', '.'))) - if self.pxsize == 0: - self.magnification = int(re.findall(r'(\d+)x', self.objective)[0]) * self.optovar[0] - self.pxsize = self.pxsizecam / self.magnification + def open(self): + if not self.path.name.startswith("Pos"): + path = self.path / f"Pos{self.series}" else: - self.magnification = self.pxsizecam / self.pxsize - self.pcf = self.shape[2] * self.metadata.re_search(r'(?i)conversion\sfactor\scoeff', 1) - self.filter = self.metadata.search('ZeissReflectorTurret-Label', self.filter)[0] + path = self.path + + filelist = sorted([file for file in path.iterdir() if re.search(r'^img_\d{3,}.*\d{3,}.*\.tif$', file.name)]) + with tifffile.TiffFile(self.path / filelist[0]) as tif: + metadata = {key: yaml.safe_load(value) for key, value in tif.pages[0].tags[50839].value.items()} + + # compare channel names from metadata with filenames + cnamelist = metadata["Info"]["Summary"]["ChNames"] + cnamelist = [c for c in cnamelist if any([c in f.name for f in filelist])] + + pattern_t = re.compile(r"img_(\d{3,})") + pattern_c = re.compile(r"img_\d{3,}_(.*)_\d{3,}$") + pattern_z = re.compile(r"(\d{3,})$") + self.filedict = {(int(pattern_t.findall(file.stem)[0]), + int(pattern_z.findall(file.stem)[0]), + cnamelist.index(pattern_c.findall(file.stem)[0])): file for file in filelist} def __frame__(self, c=0, z=0, t=0): - return tifffile.imread(os.path.join(self.path, self.filedict[(c, z, t)])) + return tifffile.imread(self.path / self.filedict[(c, z, t)]) diff --git a/ndbioimage/readers/tifread.py b/ndbioimage/readers/tifread.py index cfdff88..f03a3c4 100644 --- a/ndbioimage/readers/tifread.py +++ b/ndbioimage/readers/tifread.py @@ -1,50 +1,73 @@ -from ndbioimage import Imread, XmlData +from abc import ABC + +from ndbioimage import Imread import numpy as np import tifffile import yaml +from functools import cached_property +from ome_types import model +from pathlib import Path +from itertools import product -class Reader(Imread): +class Reader(Imread, ABC): priority = 0 do_not_pickle = 'reader' @staticmethod def _can_open(path): - if isinstance(path, str) and (path.endswith('.tif') or path.endswith('.tiff')): + if isinstance(path, Path) and path.suffix in ('.tif', '.tiff'): with tifffile.TiffFile(path) as tif: return tif.is_imagej else: return False + @cached_property + def ome(self): + metadata = {key: yaml.safe_load(value) if isinstance(value, str) else value + for key, value in self.reader.imagej_metadata.items()} + + page = self.reader.pages[0] + self.p_ndim = page.ndim + size_x = page.imagelength + size_y = page.imagewidth + if self.p_ndim == 3: + size_c = page.samplesperpixel + self.p_transpose = [i for i in [page.axes.find(j) for j in 'SYX'] if i >= 0] + size_t = metadata.get('frames', 1) # // C + else: + size_c = metadata.get('channels', 1) + size_t = metadata.get('frames', 1) + size_z = metadata.get('slices', 1) + if 282 in page.tags and 296 in page.tags and page.tags[296].value == 1: + f = page.tags[282].value + pxsize = f[1] / f[0] + else: + pxsize = None + + ome = model.OME() + ome.instruments.append(model.Instrument(id='Instrument:0')) + ome.instruments[0].objectives.append(model.Objective(id='Objective:0')) + ome.images.append( + model.Image( + id='Image:0', + pixels=model.Pixels( + id='Pixels:0', + size_c=size_c, size_z=size_z, size_t=size_t, size_x=size_x, size_y=size_y, + dimension_order="XYCZT", type=page.dtype.name, physical_size_x=pxsize, physical_size_y=pxsize), + objective_settings=model.ObjectiveSettings(id="Objective:0"))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append(model.Plane(the_c=c, the_z=z, the_t=t, delta_t=0)) + return ome + def open(self): self.reader = tifffile.TiffFile(self.path) def close(self): self.reader.close() - def __metadata__(self): - self.metadata = XmlData({key: yaml.safe_load(value) if isinstance(value, str) else value - for key, value in self.reader.imagej_metadata.items()}) - P = self.reader.pages[0] - self.pndim = P.ndim - X = P.imagelength - Y = P.imagewidth - if self.pndim == 3: - C = P.samplesperpixel - self.transpose = [i for i in [P.axes.find(j) for j in 'SYX'] if i >= 0] - T = self.metadata.get('frames', 1) # // C - else: - C = self.metadata.get('channels', 1) - T = self.metadata.get('frames', 1) - Z = self.metadata.get('slices', 1) - self.shape = (X, Y, C, Z, T) - if 282 in P.tags and 296 in P.tags and P.tags[296].value == 1: - f = P.tags[282].value - self.pxsize = f[1] / f[0] - # TODO: more metadata - def __frame__(self, c, z, t): - if self.pndim == 3: - return np.transpose(self.reader.asarray(z + t * self.shape[3]), self.transpose)[c] + if self.p_ndim == 3: + return np.transpose(self.reader.asarray(z + t * self.file_shape[3]), self.p_transpose)[c] else: - return self.reader.asarray(c + z * self.shape[2] + t * self.shape[2] * self.shape[3]) + return self.reader.asarray(c + z * self.file_shape[2] + t * self.file_shape[2] * self.file_shape[3]) diff --git a/ndbioimage/transforms.py b/ndbioimage/transforms.py index 2b10c1a..dc31b3b 100644 --- a/ndbioimage/transforms.py +++ b/ndbioimage/transforms.py @@ -1,21 +1,19 @@ import yaml -import os +import re import numpy as np from copy import deepcopy -from collections import OrderedDict +from pathlib import Path try: # best if SimpleElastix is installed: https://simpleelastix.readthedocs.io/GettingStarted.html import SimpleITK as sitk - installed = True except ImportError: - installed = False + sitk = None try: - pp = True from pandas import DataFrame, Series except ImportError: - pp = False + DataFrame, Series = None, None if hasattr(yaml, 'full_load'): @@ -24,42 +22,52 @@ else: yamlload = yaml.load -class Transforms(OrderedDict): +class Transforms(dict): def __init__(self, *args, **kwargs): super().__init__(*args[1:], **kwargs) + self.default = Transform() if len(args): self.load(args[0]) + def __mul__(self, other): + new = Transforms() + if isinstance(other, Transforms): + for key0, value0 in self.items(): + for key1, value1 in other.items(): + new[key0 + key1] = value0 * value1 + return new + elif other is None: + return self + else: + for key in self.keys(): + new[key] = self[key] * other + return new + def asdict(self): - return {f'{key[0]:.0f}:{key[1]:.0f}': value.asdict() for key, value in self.items()} + return {':'.join(str(i).replace('\\', '\\\\').replace(':', r'\:') for i in key): value.asdict() + for key, value in self.items()} def load(self, file): if isinstance(file, dict): d = file else: - if not file[-3:] == 'yml': - file += '.yml' - with open(file, 'r') as f: + with open(file.with_suffix(".yml"), 'r') as f: d = yamlload(f) + pattern = re.compile(r'[^\\]:') for key, value in d.items(): - self[tuple([int(k) for k in key.split(':')])] = Transform(value) + self[tuple(i.replace(r'\:', ':').replace('\\\\', '\\') for i in pattern.split(key))] = Transform(value) - def __call__(self, channel, time, tracks, detectors): - track, detector = tracks[channel], detectors[channel] - if (track, detector) in self: - return self[track, detector] - elif (0, detector) in self: - return self[0, detector] - else: - return Transform() + def __missing__(self, key): + return self.default - def __reduce__(self): - return self.__class__, (self.asdict(),) + def __getstate__(self): + return self.__dict__ + + def __setstate__(self, state): + self.__dict__.update(state) def save(self, file): - if not file[-3:] == 'yml': - file += '.yml' - with open(file, 'w') as f: + with open(file.with_suffix(".yml"), 'w') as f: yaml.safe_dump(self.asdict(), f, default_flow_style=None) def copy(self): @@ -68,6 +76,7 @@ class Transforms(OrderedDict): def adapt(self, origin, shape): for value in self.values(): value.adapt(origin, shape) + self.default.adapt(origin, shape) @property def inverse(self): @@ -76,15 +85,20 @@ class Transforms(OrderedDict): inverse[key] = value.inverse return inverse + @property + def ndim(self): + return len(list(self.keys())[0]) + class Transform: def __init__(self, *args): - if not installed: - raise ImportError('SimpleElastix is not installed: https://simpleelastix.readthedocs.io/GettingStarted.html') - self.transform = sitk.ReadTransform(os.path.join(os.path.dirname(__file__), 'transform.txt')) - self.dparameters = (0, 0, 0, 0, 0, 0) - self.shape = (512, 512) - self.origin = (255.5, 255.5) + if sitk is None: + raise ImportError('SimpleElastix is not installed: ' + 'https://simpleelastix.readthedocs.io/GettingStarted.html') + self.transform = sitk.ReadTransform(str(Path(__file__).parent / 'transform.txt')) + self.dparameters = 0, 0, 0, 0, 0, 0 + self.shape = 512, 512 + self.origin = 255.5, 255.5 if len(args) == 1: # load from file or dict if isinstance(args[0], np.ndarray): self.matrix = args[0] @@ -92,7 +106,7 @@ class Transform: self.load(*args) elif len(args) > 1: # make new transform using fixed and moving image self.register(*args) - self._last = None + self._last, self._inverse = None, None def __mul__(self, other): # TODO: take care of dmatrix result = self.copy() @@ -114,13 +128,13 @@ class Transform: return deepcopy(self) @staticmethod - def castImage(im): + def cast_image(im): if not isinstance(im, sitk.Image): im = sitk.GetImageFromArray(im) return im @staticmethod - def castArray(im): + def cast_array(im): if isinstance(im, sitk.Image): im = sitk.GetArrayFromImage(im) return im @@ -190,7 +204,7 @@ class Transform: dtype = im.dtype im = im.astype('float') intp = sitk.sitkBSpline if np.issubdtype(dtype, np.floating) else sitk.sitkNearestNeighbor - return self.castArray(sitk.Resample(self.castImage(im), self.transform, intp, default)).astype(dtype) + return self.cast_array(sitk.Resample(self.cast_image(im), self.transform, intp, default)).astype(dtype) def coords(self, array, columns=None): """ Transform coordinates in 2 column numpy array, @@ -198,7 +212,7 @@ class Transform: """ if self.is_unity(): return array.copy() - elif pp and isinstance(array, (DataFrame, Series)): + elif DataFrame is not None and isinstance(array, (DataFrame, Series)): columns = columns or ['x', 'y'] array = array.copy() if isinstance(array, DataFrame): @@ -239,7 +253,7 @@ class Transform: """ kind = kind or 'affine' self.shape = fix.shape - fix, mov = self.castImage(fix), self.castImage(mov) + fix, mov = self.cast_image(fix), self.cast_image(mov) # TODO: implement RigidTransform tfilter = sitk.ElastixImageFilter() tfilter.LogToConsoleOff() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..a0ba74b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,37 @@ +[tool.poetry] +name = "ndbioimage" +version = "2023.6.0" +description = "Bio image reading, metadata and some affine registration." +authors = ["W. Pomp "] +license = "GPLv3" +readme = "README.md" +keywords = ["bioformats", "imread", "numpy", "metadata"] +include = ["transform.txt"] +repository = "https://github.com/wimpomp/ndbioimage" + +[tool.poetry.dependencies] +python = "^3.8" +numpy = "*" +pandas = "*" +tifffile = "*" +czifile = "*" +tiffwrite = "*" +ome-types = "*" +pint = "*" +tqdm = "*" +lxml = "*" +pyyaml = "*" +parfor = "*" +JPype1 = "*" +SimpleITK-SimpleElastix = "*" +pytest = { version = "*", optional = true } + +[tool.poetry.extras] +test = ["pytest-xdist"] + +[tool.poetry.scripts] +ndbioimage = "ndbioimage:main" + +[build-system] +requires = ["poetry-core"] +build-backend = "poetry.core.masonry.api" diff --git a/setup.py b/setup.py deleted file mode 100644 index abf1674..0000000 --- a/setup.py +++ /dev/null @@ -1,48 +0,0 @@ -import setuptools -import platform -import os - -version = '2022.7.1' - -if platform.system().lower() == 'linux': - import pkg_resources - pkg_resources.require(['pip >= 20.3']) - -with open('README.md', 'r') as fh: - long_description = fh.read() - - -with open(os.path.join(os.path.dirname(__file__), 'ndbioimage', '_version.py'), 'w') as f: - f.write(f"__version__ = '{version}'\n") - try: - with open(os.path.join(os.path.dirname(__file__), '.git', 'refs', 'heads', 'master')) as h: - f.write("__git_commit_hash__ = '{}'\n".format(h.read().rstrip('\n'))) - except Exception: - f.write(f"__git_commit_hash__ = 'unknown'\n") - - -setuptools.setup( - name='ndbioimage', - version=version, - author='Wim Pomp', - author_email='w.pomp@nki.nl', - description='Bio image reading, metadata and some affine registration.', - long_description=long_description, - long_description_content_type='text/markdown', - url='https://github.com/wimpomp/ndbioimage', - packages=['ndbioimage', 'ndbioimage.readers'], - classifiers=[ - 'Programming Language :: Python :: 3', - 'License :: OSI Approved :: GNU General Public License v3 (GPLv3)', - 'Operating System :: OS Independent', - ], - python_requires='>=3.7', - install_requires=['untangle', 'pandas', 'psutil', 'numpy', 'tqdm', 'tifffile', 'czifile', 'pyyaml', 'dill', - 'tiffwrite'], - extras_require={'transforms': 'SimpleITK-SimpleElastix', - 'bioformats': ['python-javabridge', 'python-bioformats']}, - tests_require=['pytest-xdist'], - entry_points={'console_scripts': ['ndbioimage=ndbioimage.ndbioimage:main']}, - package_data={'': ['transform.txt']}, - include_package_data=True, -) diff --git a/tests/test_open.py b/tests/test_open.py new file mode 100644 index 0000000..cd3883b --- /dev/null +++ b/tests/test_open.py @@ -0,0 +1,10 @@ +import pytest +from pathlib import Path +from ndbioimage import Imread + + +@pytest.mark.parametrize("file", + [file for file in (Path(__file__).parent / 'files').iterdir() if file.suffix != '.pzl']) +def test_open(file): + with Imread(file) as im: + print(im[dict(c=0, z=0, t=0)].mean())