Source code for owl.mcam_data._export

import concurrent.futures
import functools
import math
import sys
import tempfile
from contextlib import nullcontext
from functools import partial
from itertools import chain
from math import prod
from pathlib import Path
from warnings import warn

try:
    from natsort import natsorted
except ImportError:
    print(f"""Install natsort

    {sys.prefix}/bin/conda install --prefix {sys.prefix} natsort

Falling back on sorted.
""")
    natsorted = None
import h5netcdf
import imageio.v3 as imageio
import numpy as np
import ro_json as json
import tifffile
import xarray as xr
from dask.diagnostics import ProgressBar
from packaging.version import Version

try:
    from wabisabi import kwarg_name_change
except ImportError:
    def kwarg_name_change(*args, **kwargs):
        def decorator(func):
            return func

        return decorator
    print(f"""Install wabisabi

    conda install --prefix {sys.prefix} wabisabi

""")

import owl
from owl.color.illumination_colormaps import generate_luts
from owl.util import get_localzone_name
from owl.util.io_utils import check_disk_space
from owl.util.util import visual_illumination_modes

from ..array._imageio_backend import FlatImageioArray
from ..array._tiff_backend import TiledMultiPageTiffArray
from ..array._video_array import TiledVideoArray
from ..util import make_timestamp
from ..util.index_tricks import ndrange
from ._export_constants import (
    _default_imagename_format,
    _default_imagename_format_wells,
    _default_imagename_format_wells_field_id,
    _default_tile_dims,
    _expected_compression_factor,
    _valid_chroma_for_image_export_format,
    _valid_extensions,
    _valid_extensions_for_image_export_format,
)
from ._format import DatasetFormatMapping, format_string_array
from ._mp4 import _export_multi_mp4, _export_tiledmp4
from ._netcdf import write_netcdf
from ._new import new_dataset
from ._optimized_tifffile import FileHandle
from ._properties import get_bin_mode, get_chroma, get_stack_dims, is_well_dataset
from ._update import update_dataset
from ._update_export_metadata import update_metadata
from ._util import valid_bayer_patterns

kwarg_name_change = functools.partial(
    kwarg_name_change,
    library_name='owl',
    current_library_version=owl.__version__,
)

bad_h5netcdf_version = h5netcdf.__version__ in [
    # Release on Oct 15, 2025 -- broke serilization of tuples strings in attrs
    '1.7.0',
]

if bad_h5netcdf_version:
    # End users really shouldn't be hitting this
    # So this is mostly for our internal users.
    print(f"""
You are using h5netcdf version {h5netcdf.__version__} which has a bug that breaks
saving data. Contact us at help@ramonaoptics.com to resolve this issue or through slack.
""")


def _tqdm(iterable, *args, **kwargs):
    return iterable


def _export_images_or_image_stacks(
    mcam_dataset,
    *,
    directory,
    metadata_filename,
    mode,
    tile_dims,
    save_tiff_references=None,
    tqdm=None,
    **kwargs,
):
    # Metadata may be written after export (for tiff refs), so ensure the directory exists.
    directory.mkdir(parents=True, exist_ok=True)

    image_export_format = mcam_dataset.attrs['image_export_format']

    if image_export_format == 'flat_image_stacks' and tile_dims is None:
        raise ValueError("flat_image_stacks require tile_dims to be specified at this stage.")
    elif tile_dims is None:
        tile_dims = _default_tile_dims[image_export_format]

    imagename_format = mcam_dataset.attrs['imagename_format']
    extension = Path(imagename_format).suffix.lower()

    if save_tiff_references is None:
        save_tiff_references = extension in ('.tif', '.tiff')

    if save_tiff_references and extension not in ('.tif', '.tiff'):
        raise RuntimeError(
            f'save_tiff_references is only compatible with .tif files. not {extension} files')

    chunks = mcam_dataset.images.chunks
    if chunks is None:
        # If there are no chunks, lets take the trivial case,
        # where each image is randomly accessible
        # as its own chunk
        chunks = tuple(
            (1,) * mcam_dataset.sizes[d] if d not in ("y", "x", "rgb", "rgba")
            else (mcam_dataset.sizes[d],)
            for d in mcam_dataset.images.dims
        )

    tile_chunks = tuple(
        chunks[mcam_dataset.images.dims.index(d)]
        for d in tile_dims
    )

    shape = tuple(
        mcam_dataset.sizes[d]
        for d in tile_dims
    )
    this_valid_data = np.full(shape, fill_value=True, dtype='bool')
    # Sometimes reading from the data is locked to a single threaded operation
    # And once the first job has been submitted
    # It may lock up the H5PY file for access making it impossible to submit
    # other jobs
    # In this case, it isn't about performance, but more the fact that
    # the end user doesn't see any indication of progress making them feel
    # like they are waiting too long.
    if 'valid_data' in mcam_dataset:
        valid_data = mcam_dataset.valid_data.compute()
        for i in ndrange(this_valid_data.shape):
            this_valid_data[i] = np.any(np.asarray(
                np.asarray(valid_data.isel({
                    t: s
                    for t, s in zip(tile_dims, i)
                    if t in valid_data.dims
                }))
            ))

    # This optimization was introduced to help users go from
    # tiledmp4 -> flat_image_stacks
    # it may not be completely correct for other cases
    # However tests were added in
    # tests/mcam_data/test_mcam_data_export_performance.py
    # to help ensure that if this algorithm changes the requirements
    # to help ensure that the new algorithm stays performant are there
    if all(c == (1,) * len(c) for c in tile_chunks):
        # This is the trivial case where we can use our likely lazy arrays
        # and a general concurrent future pipeline to simply write out
        # the entire data
        _trivial_export_images_or_image_stacks(
            mcam_dataset,
            directory=directory,
            imagename_format=imagename_format,
            mode=mode,
            tile_dims=tile_dims,
            this_valid_data=this_valid_data,
            tqdm=tqdm,
            **kwargs,
        )
    else:
        # Here we try to ensure that we make use of data
        # locally once the data is loaded into an array
        # This algorithm isn't perfect but we can improve on it
        # as time goes on
        # It focuses on the case where the data is in a large chunk
        # Across a dimension like the frame_number, and should only
        # be read once across an entire x, y dimension
        _export_images_or_image_stacks_chunked(
            mcam_dataset,
            directory=directory,
            imagename_format=imagename_format,
            mode=mode,
            tile_dims=tile_dims,
            tqdm=tqdm,
            this_valid_data=this_valid_data,
            **kwargs,
        )
    if metadata_filename is not None:
        mcam_dataset = _add_imagename_cache_to_metadata(mcam_dataset)
        if save_tiff_references:
            mcam_dataset = _add_tiff_refs_to_metadata(mcam_dataset, directory)

        save_metadata(
            mcam_dataset,
            directory=directory,
            include_timestamp=False,
            metadata_filename=metadata_filename,
            mode=mode,
        )

    return directory


def _export_images_or_image_stacks_chunked(
    mcam_dataset,
    *,
    directory,
    imagename_format,
    mode,
    tile_dims,
    tqdm,
    this_valid_data,
    max_workers=4,
    **kwargs,
):
    total_images = prod(
        mcam_dataset.sizes[d]
        for d in tile_dims
    )

    this_valid_data_da = xr.DataArray(
        dims=tile_dims,
        data=this_valid_data,
    )

    # We consider a few cases:
    tile_stack_dims = tuple(
        d
        for d in tile_dims
        if d not in ('image_y', 'image_x')
    )
    tile_stack_shape = tuple(
        mcam_dataset.sizes[d]
        for d in tile_stack_dims
    )

    remaining_dims = tuple(
        d
        for d in tile_dims
        if d not in tile_stack_dims
    )
    remaining_shape = tuple(
        mcam_dataset.sizes[d]
        for d in remaining_dims
    )
    if tqdm is None:
        class tqdm:
            def __init__(self, iterable=None, *args, **kwargs):
                self._iterable = iterable
                return

            def update(self, *args, **kwargs):
                pass

            def close(self, *args, **kwargs):
                pass

            def __iter__(self):
                return iter(self._iterable)

    tqdm = tqdm(total=total_images)
    with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as executor:
        for i_t in ndrange(tile_stack_shape):
            this_chunk_dataset = mcam_dataset.isel({
                d: i
                for d, i in zip(tile_stack_dims, i_t)
            }).compute()
            futures = []
            for i_r in ndrange(remaining_shape):
                this_dataset = this_chunk_dataset.isel({
                    d: i
                    for d, i in zip(remaining_dims, i_r)
                })
                if not this_valid_data_da.isel({
                    d: i
                    for d, i in chain(
                        zip(tile_stack_dims, i_t),
                        zip(remaining_dims, i_r),
                    )
                }).any():
                    continue
                futures.append(executor.submit(
                    _save_one,
                    this_dataset,
                    directory=directory,
                    imagename_format=imagename_format,
                    mode=mode,
                    **kwargs,
                ))

            for future in concurrent.futures.as_completed(futures):
                future.result()
                tqdm.update()
    tqdm.close()


def _trivial_export_images_or_image_stacks(
    mcam_dataset,
    *,
    directory,
    imagename_format,
    mode,
    tile_dims,
    tqdm,
    this_valid_data,
    max_workers=4,
    **kwargs,
):
    if tqdm is None:
        class tqdm:
            def __init__(self, iterable=None, *args, **kwargs):
                self._iterable = iterable
                return

            def update(self, *args, **kwargs):
                pass

            def close(self, *args, **kwargs):
                pass

            def __iter__(self):
                return iter(self._iterable)

    shape = tuple(
        mcam_dataset.sizes[t]
        for t in tile_dims
    )

    # Setting max workers beyond 4 will likely causing thrasing
    # since SSDs will get eaten up
    # Compression algorithms also use a lot of CPU
    # so its best to have a "small number" greater than 1
    # of workers to keep the CPU busy here.
    futures = []
    with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as executor:
        for i in ndrange(shape):
            if not this_valid_data[i]:
                continue
            this_dataset = mcam_dataset.isel({
                t: s
                for t, s in zip(tile_dims, i)
            })
            futures.append(executor.submit(
                _save_one,
                this_dataset,
                directory=directory,
                imagename_format=imagename_format,
                mode=mode,
                **kwargs,
            ))

        for future in tqdm(concurrent.futures.as_completed(futures), total=len(futures)):
            future.result()


_exporters = {
    'images': _export_images_or_image_stacks,
    # 'flat_image_stacks' is available as an image_export_format but
    # it is a shortcut for 'images' + flatten all channels
    'flat_image_stacks': _export_images_or_image_stacks,
    'tiledmp4': _export_tiledmp4,
    'multi_mp4': _export_multi_mp4,
}

_supported_image_export_formats = tuple(_exporters.keys())


def _save_one(dataset, directory, imagename_format, mode='w', **kwargs):
    extension = Path(imagename_format).suffix.lower()
    if extension in ('.tif', '.tiff'):
        func = save_one_as_tif
    elif len(get_stack_dims(dataset)) > 0:
        raise ValueError(
            f"We do not support image stacks for extension {extension}")
    elif extension == '.j2c':
        func = save_one_as_j2c
    elif extension == '.png':
        func = save_one_as_png
    elif extension == '.bmp':
        func = save_one_as_bmp
    elif extension in ('.jpeg', '.jpg'):
        func = save_one_as_jpeg
    else:
        # This function is typically called from within dask
        # therefore this error doesn't get to the end user nicely
        # so it would be better if people raised the extension before  calling
        # this function.
        raise ValueError(
            f"Unknown Extension {extension} provided in imagename_format")
    return func(
        dataset,
        directory,
        imagename_format=imagename_format,
        mode=mode,
        **kwargs
    )


def save_one_as_jpeg(
    dataset,
    directory: Path,
    *,
    imagename_format="cam{image_y}_{image_x}.jpeg",
    mode='w',
    **kwargs,
):
    """Save a single image as a jpeg.

    Parameters
    ----------
    dataset: xarray.Dataset
        The image you would like to save. The image is contained in the
        ``'images'`` variable in the dataset, and should have 2 or 3
        dimensions.

    directory:
        The directory where you would like to save the image.

    imagename_format:
        A format string that dictates the filename format how the image
        should be stored. The extension is used by `imageio` to dictate the
        compression. This should contain two formatting parameters
        `r` and `c` describing the row index and the column index of the
        image. `image_y` and `image_x` are passed to this format string.

    mode: str
        The mode to open the file in. Valid options are ``'w'`` (write),
        or ```x``` exclusive creation, failing if the file already exists.

    """
    if not imagename_format.lower().endswith(('.jpeg', '.jpg')):
        raise ValueError(f'Given ``imagename_format`` {imagename_format} does '
                         'not end with a valid jpeg extension')

    filename = imagename_format.format_map(DatasetFormatMapping(dataset))
    path = (directory / filename).resolve()
    path.parent.mkdir(parents=True, exist_ok=True)

    _save_jpeg_with_metadata(
        np.asarray(dataset.images),
        path,
        dataset,
        include_timestamp=False,
        mode=mode,
    )


def save_one_as_bmp(dataset, directory: Path, *,
                    imagename_format="cam{image_y}_{image_x}.bmp", mode='w', **kwargs):
    """Save a single image as a bmp.

    Parameters
    ----------
    dataset: xarray.Dataset
        The image you would like to save. The image is contained in the
        ``'images'`` variable in the dataset, and should have 2 or 3
        dimensions.

    directory:
        The directory where you would like to save the image.

    imagename_format:
        A format string that dictates the filename format how the image
        should be stored. The extension is used by `imageio` to dictate the
        compression. This should contain two formatting parameters
        `r` and `c` describing the row index and the column index of the
        image. `image_y` and `image_x` are passed to this format string.

    """

    filename = imagename_format.format_map(DatasetFormatMapping(dataset))
    path = (directory / filename).resolve()
    path.parent.mkdir(parents=True, exist_ok=True)

    with imageio.imopen(str(path), io_mode=mode, legacy_mode=False) as image_file:
        image_file.write(np.asarray(dataset.images))


def save_one_as_png(dataset, directory: Path, *,
                    imagename_format="cam{image_y}_{image_x}.png",
                    mode='w',
                    **kwargs):
    """Save a single image as a png.

    Parameters
    ----------
    dataset: xarray.Dataset
        The image you would like to save. The image is contained in the
        ``'images'`` variable in the dataset, and should have 2 or 3
        dimensions.

    directory:
        The directory where you would like to save the image.

    imagename_format:
        A format string that dictates the filename format how the image
        should be stored. The extension is used by `imageio` to dictate the
        compression. This should contain two formatting parameters
        `r` and `c` describing the row index and the column index of the
        image. `image_y` and `image_x` are passed to this format string.

    mode: str
        The mode to open the file in. Valid options are ``'w'`` (write),
        or ```x``` exclusive creation, failing if the file already exists.

    """
    if not imagename_format.lower().endswith('.png'):
        raise ValueError(f'Given ``imagename_format`` {imagename_format} does '
                         'not end with a valid png extension')

    filename = imagename_format.format_map(DatasetFormatMapping(dataset))
    path = (directory / filename).resolve()
    path.parent.mkdir(parents=True, exist_ok=True)

    _save_png_with_metadata(
        np.asarray(dataset.images),
        path,
        dataset,
        include_timestamp=False,
        mode=mode,
    )


def save_one_as_j2c(
    dataset,
    directory: Path,
    imagename_format="cam{image_y}_{image_x}.j2c",
    mode='w',
    *,
    num_decompositions=None,
    reversible=None,
    qstep=None,
    **kwargs,
):
    # This saves a "jpeg2000 codestream", no other metadata is stored in this
    # it isn't a full file, and is typically included in a larger container
    if not imagename_format.lower().endswith('.j2c'):
        raise ValueError(f'Given ``imagename_format`` {imagename_format} does '
                         'not end with a valid .j2c extension')

    filename = imagename_format.format_map(DatasetFormatMapping(dataset))
    path = (directory / filename).resolve()
    path.parent.mkdir(parents=True, exist_ok=True)

    if mode == 'x':
        if path.exists():
            raise FileExistsError(
                f"File {path} already exists, cannot open in exclusive mode."
            )
        mode = 'w'

    import ojph
    if Version(ojph.__version__) < Version('0.5.1'):
        raise RuntimeError(f"""You need ojph 0.5.1.

This should only happen in development mode. To fix this use

    conda install "ojph>=0.5.1" --prefix {sys.prefix}

""")
    data = np.asarray(dataset.images)
    if num_decompositions is None:
        num_decompositions = max(
            1,
            min(
                31,
                # Smallest size is 16 x 16 pixels, or about that!
                int(np.log2(max(data.shape[:2]))) - 3
            )
        )

    if data.ndim == 2:
        channel_order = 'HW'
    elif data.ndim == 3:
        channel_order = 'HWC'
    else:
        raise ValueError(f"Unsupported number of dimensions: {data.ndim}")

    ojph.imwrite(
        path,
        data,
        channel_order=channel_order,
        num_decompositions=num_decompositions,
        reversible=reversible,
        qstep=qstep,
        # These two parameters were introduced in 0.4.6.
        # The if statement should catch bad installations
        # But these parameters are critical for fast reading of
        # JPEG2000 files using OpenJPH, so we re-assert them here
        #
        # tlm_marker provides us with information in the header
        # to be able to estimate the size of the file needed
        tlm_marker=True,
        # tileparts_at_resolutions provides us with one tilepart for each
        # resolution level, allowing us to read just the necessary data
        tileparts_at_resolutions=True,
    )


def save_one_as_tif(
    dataset,
    directory: Path, *,
    imagename_format="cam{image_y}_{image_x}.tif",
    mode='w',
    **kwargs,
):
    """Save a single image as a tif.

    Parameters
    ----------
    dataset: xarray.Dataset
        The image you would like to save. The image is contained in the
        ``'images'`` variable in the dataset, and should have 2 or 3
        dimensions.

    directory:
        The directory where you would like to save the image.

    imagename_format:
        A format string that dictates the filename format how the image
        should be stored. The extension must be some form of the tif file
        extension. This should contain two formatting parameters
        `r` and `c` describing the row index and the column index of the
        image. `image_y` and `image_x` are passed to this format string.

    """
    if not (imagename_format.lower().endswith('tif') or
            imagename_format.lower().endswith('tiff')):
        raise ValueError(f'Given ``imagename_format`` {imagename_format} does '
                         'not end with a valid tif extension')

    filename = imagename_format.format_map(DatasetFormatMapping(dataset))
    path = (directory / filename).resolve()
    path.parent.mkdir(parents=True, exist_ok=True)

    data = np.asarray(dataset.images)
    _save_tiff_with_metadata(
        data,
        str(path),
        dataset,
        include_timestamp=False,
        mode=mode,
        **kwargs,
    )


def _load_json_metdata(filename):
    with open(filename, encoding='utf-8') as metadata_file:
        metadata = json.load(metadata_file, ignore_comments=True)
    return xr.Dataset.from_dict(metadata)


def _load_zarr_metadata(filename, delayed=False, chunks='auto'):
    metadata = xr.open_zarr(filename, chunks=chunks)
    if delayed is None:
        delayed = metadata.nbytes > 1E9
    if not delayed:
        metadata.load()
    return metadata


def _load_netcdf_metadata(filename, delayed=False, engine='ramona'):
    # I'm not really sure, but I feel like libnetcdf doesn't close files
    # correctly and this results in bugs in HDF5 like:
    # HDF5-DIAG: Error detected in HDF5 (1.14.0) thread 1:
    #     #000: H5A.c line 679 in H5Aopen_by_name(): unable to synchronously open attribute
    #     major: Attribute
    #     minor: Can't open object
    # By using "h5netcdf" it seems to avoid this bug....
    #
    # The metadata file is small, but its HDF5 metadata is scattered across the
    # file, so a cold open issues many small random reads. The ramona engine
    # opens files with the HDF5 'direct' driver (O_DIRECT) -- intended for
    # streaming the large image-data files -- which bypasses the OS page cache
    # and read-ahead and turns each of those small reads into a physical disk
    # seek. Metadata used to be loaded with the buffered "h5netcdf" engine;
    # !9398 (Mar 2026, "Remove dependency on libnetcdf4") flipped the default to
    # the ramona engine, so metadata started inheriting O_DIRECT as a side
    # effect. Open the small metadata file with the buffered h5netcdf engine so
    # the reads go through the page cache and the OS handles caching/read-ahead;
    # the large image data keeps O_DIRECT, where it belongs.
    metadata_engine = 'h5netcdf' if engine == 'ramona' else engine
    metadata = xr.open_dataset(filename, engine=metadata_engine)
    if delayed is None:
        delayed = metadata.nbytes > 1E9
    if not delayed:
        metadata.load()
    return metadata


def _load_metadata(directory, metadata_filename, delayed=False, engine='ramona'):
    """Load metadata handling mostly for the legacy logic.

    This function should be heavily documented as it is probably going to be
    confusing to read for future developers.
    """

    if metadata_filename is None:
        if (directory / 'metadata.json').is_file():
            metadata_filename = 'metadata.json'
        elif (directory / 'metadata.nc').is_file():
            metadata_filename = 'metadata.nc'
        elif (directory / 'analysis_metadata.nc').is_file():
            metadata_filename = 'analysis_metadata.nc'
        elif (directory / 'metadata.zarr').is_dir():
            metadata_filename = 'metadata.zarr'
        elif (directory / 'analysis_metadata.zarr').is_dir():
            metadata_filename = 'analysis_metadata.zarr'
        else:
            raise RuntimeError(
                f"Could not find metadata file in provided directory ({directory}).")

    filename = directory / metadata_filename
    suffix = filename.suffix.lower()

    try:
        if suffix == ".json":
            metadata_dataset = _load_json_metdata(filename)
        elif suffix == ".nc":
            metadata_dataset = _load_netcdf_metadata(filename, delayed=delayed, engine=engine)
        elif suffix == ".zarr":
            metadata_dataset = _load_zarr_metadata(filename, delayed=delayed)
        else:
            raise ValueError(f"Unknown suffix ({suffix}) for metadata.")
        # Ideally we would return here, but unfortunately
        # We change the way the data is saved every day!
        # 2020/11/18: metadata is simply extracted from xarray dataset using to_dict().
        metadata_dataset = update_metadata(metadata_dataset)
        return metadata_dataset
    except Exception as e:
        raise RuntimeError("We are sorry that you are unable to open this dataset. "
                           "To get help, please contact Ramona Optics at help@ramonaoptics.com "
                           "with a screenshot of your screen including this error message.") from e


[docs] def save_metadata(mcam_dataset, directory, *, include_timestamp: bool = True, metadata_filename: str = 'metadata.nc', allow_nan: bool = True, mode='w', tqdm=_tqdm, vars_to_chunk=None, virtual_dataset=None, _stacklevel_increment=0, ): """Save mcam metadata to a given directory. Parameters ---------- mcam_dataset: xarray Dataset The mcam_dataset in an xarray Dataset with all the metadata. directory: pathlike The directory where to save the metadata. metadata_filename: str, optional The name of the metadata file. include_timestamp: bool, optional If set to `True`, this will append a timestamp to the provided directory. If set to false, no timestamp will be appended to the directory possibly overwriting any files currently within the previous directory. allow_nan: bool, optional If True nan values will be allowed when exporting metadata as a json. If False attempting to export metadata as a json will result in an exception. Returns ------- save_directory: Path The directory where the metadata was saved. If ``include_timestamp`` is ``True``, then this include the added timestamp. """ for dim in mcam_dataset.dims: if dim not in mcam_dataset.coords: warn( f'{dim} dimension has not been given coordinates within the dataset. ' 'Please contact Ramona Optics at help@ramonaoptics.com to get this resolved. ' 'Please include the dataset and mention #1644 in your message.', stacklevel=2 + _stacklevel_increment, ) directory = Path(directory) if include_timestamp: directory = (directory.parent / (directory.name + '_' + make_timestamp())) directory.mkdir(parents=True, exist_ok=True) if metadata_filename is None: # There is more metadata here, and therefore we typically have # found that we save the metadata as nc metadata_filename = 'metadata.nc' metadata_filename = metadata_filename.format_map( DatasetFormatMapping(mcam_dataset, {'default_filename': 'metadata'}) ) metadata = mcam_dataset.drop_vars('images', errors='ignore') _write_metadata( metadata, directory / metadata_filename, allow_nan=allow_nan, mode=mode, tqdm=tqdm, vars_to_chunk=vars_to_chunk, virtual_dataset=virtual_dataset, ) return directory
def _write_metadata(metadata, filename, allow_nan=True, mode='w', tqdm=_tqdm, vars_to_chunk=None, virtual_dataset=None): if hasattr(filename, 'write'): # BytesIO object. No suffix here. suffix = None else: # Filepath object. Get the suffix. suffix = filename.suffix.lower() if suffix == '': filename = filename.with_suffix('.nc') suffix = '.nc' if suffix == ".zarr": metadata.to_zarr(filename, mode=mode, zarr_format=3) elif suffix == ".json": metadata_dict = metadata.to_dict() with open(filename, mode=mode, encoding='utf-8') as metadata_file: json.dump(metadata_dict, metadata_file, indent=4, allow_nan=allow_nan) elif suffix == ".nc" or suffix is None: if vars_to_chunk is None: vars_to_chunk = [] # Ensure that any analysis output masks are chunked for var_name in metadata.data_vars: var_data = metadata[var_name] if (var_data.attrs.get('__owl_settings_type__') == 'analysis_output_mask' and var_name not in vars_to_chunk): vars_to_chunk.append(var_name) write_netcdf( metadata, filename, mode=mode, tqdm=tqdm, vars_to_chunk=vars_to_chunk, virtual_dataset=virtual_dataset, ) else: raise ValueError(f"Unknown suffx ({suffix}) for saving metadata")
[docs] def load_timelapse_directory( directory, metadata_filename=None, timelapse_tqdm=None, tqdm=None, sort_reverse=None, engine='ramona', ): if natsorted is None: raise RuntimeError( "Timelapse loading requires natsort" ) directory = Path(directory) if metadata_filename is None: metadata_filename = "metadata.nc" metadata_filename = Path(metadata_filename) if timelapse_tqdm is None: timelapse_tqdm = iter timelapse_dim = "timelapse_index" parents_in_metadata_filename = len(metadata_filename.parents) - 1 timepoint_dataset_paths = natsorted( [ p.parents[parents_in_metadata_filename] for p in directory.glob(f"*/{str(metadata_filename)}") if p.parents[parents_in_metadata_filename].is_dir() ] ) if not timepoint_dataset_paths: raise ValueError(f"No timepoint datasets found in {directory}") if not timepoint_dataset_paths: return None metadata_list = [] virtual_dataset = {} data_vars_to_drop = ["images"] data_vars_to_expand = [ "software_timestamp", "_images_imagename_cache", "_images_tiff_page_offsets", "_images_tiff_page_nbytes", "_images_tiff_page_compression", "_images_tiff_page_predictor", "_images_tiff_page_dtype", "_images_tiff_page_levels_downsampling", ] analysis_output_setting_types = { "analysis_output_bounding_box", "analysis_output_mask", "analysis_output_label", "analysis_output_positioned_label", "analysis_output_animal_tracking_points", } timelapse_prefix = [] for i, timepoint_dataset_path in enumerate( timelapse_tqdm(timepoint_dataset_paths) ): timelapse_prefix.append( str(timepoint_dataset_path.relative_to(directory)) ) this_metadata = _load_metadata( timepoint_dataset_path, metadata_filename, delayed=True, engine=engine ) for key, var in this_metadata.items(): settings_type = var.attrs.get("__owl_settings_type__") if settings_type == "analysis_output_mask": if key not in virtual_dataset: data_vars_to_drop.append(key) virtual_dataset[key] = { "filenames": [], "coords": { dim: var.coords[dim].data for dim in var.dims }, "attrs": dict(var.attrs), "shape": var.data.shape, "dtype": var.data.dtype, "dims": (timelapse_dim,) + var.dims, "virtual_dims": (timelapse_dim,), } virtual_dataset[key]["filenames"].append( str((timepoint_dataset_path / metadata_filename).resolve()) ) elif settings_type in analysis_output_setting_types: data_vars_to_expand.append(key) expand_vars = list(set(data_vars_to_expand) - set(data_vars_to_drop)) this_metadata = this_metadata.drop_vars(data_vars_to_drop, errors="ignore") this_metadata.coords[timelapse_dim] = i for key in expand_vars: if key not in this_metadata: continue this_metadata[key] = this_metadata[key].expand_dims(timelapse_dim) this_metadata.load() metadata_list.append(this_metadata) if sort_reverse is not None: timelapse_prefix, metadata_list = zip(*natsorted(sorted( zip(timelapse_prefix, metadata_list), key=lambda x: x[0], reverse=sort_reverse, ), key=lambda x: x[0], reverse=sort_reverse)) # Convert to numpy array after sorting to ensure correct dtype timelapse_prefix = np.asarray(timelapse_prefix, dtype="str") timelapse_metadata = xr.concat( metadata_list, dim=timelapse_dim, coords="minimal", data_vars="minimal", compat="override", join="outer", ) timelapse_metadata["timelapse_prefix"] = ( (timelapse_dim,), timelapse_prefix, ) timelapse_metadata.attrs["images_dims"] = [ timelapse_dim ] + timelapse_metadata.attrs["images_dims"] timelapse_metadata.attrs["image_export_tile_dims"] = tuple([ timelapse_dim ] + timelapse_metadata.attrs["image_export_tile_dims"]) # We likely don't need to do the whole resolve and relative # to, but it makes the directories more readable # for future inspection timelapse_metadata.attrs["imagename_format"] = str( (directory / ( "{timelapse_prefix}/" + f"{str(Path(metadata_filename).parents[0])}/" + timelapse_metadata.attrs["imagename_format"] )).resolve().relative_to( directory.resolve() ) ) timelapse_metadata = timelapse_metadata.assign_coords( {timelapse_dim: range(len(metadata_list))} ) timelapse_metadata = _add_imagename_cache_to_metadata(timelapse_metadata) # These tiff reference should have already been aggregated # timelapse_metadata = _add_tiff_refs_to_metadata(timelapse_metadata, directory) for _, var_info in virtual_dataset.items(): var_info["coords"][timelapse_dim] = np.arange(len(metadata_list)) var_info["shape"] = (len(metadata_list),) + var_info["shape"] # HDF5 virtual datasets segfault when read from in-memory BytesIO # files, so we write to a real temp file and load it into memory. tmp = tempfile.NamedTemporaryFile(suffix='.nc', delete=False) tmp.close() tmp_path = Path(tmp.name) try: _write_metadata( timelapse_metadata, tmp_path, virtual_dataset=virtual_dataset or None, ) if engine is None: engine = "h5netcdf" metadata_ds = xr.open_dataset(tmp_path, engine=engine) finally: # TODO: cleanup these files on windows? # Mark and Clay ran out of time on March 16, 2026 try: tmp_path.unlink(missing_ok=True) except Exception as e: # https://gitlab.com/ramonaoptics/mcam/python-owl/-/work_items/1735 print(f"Failed to delete temporary file {tmp_path}: {e}") mcam_dataset = _load_directory( directory, metadata_ds, update=True, tqdm=tqdm, ) return mcam_dataset
def _ensure_2d_dataset(mcam_dataset): if ('image_x' in mcam_dataset.images.dims) and ('image_y' in mcam_dataset.images.dims): return image_data = mcam_dataset.images if 'image_x' not in image_data.dims and 'image_y' not in image_data.dims: image_data = image_data.expand_dims('image_x', 0).expand_dims('image_y', 0) elif 'image_x' not in image_data.dims: image_data = image_data.expand_dims('image_x', 1) elif 'image_y' not in image_data.dims: image_data = image_data.expand_dims('image_y', 0) mcam_dataset['images'] = image_data def save( mcam_dataset, directory: Path, *, mode='w', metadata_filename=None, include_timestamp: bool = True, imagename_format=None, image_mode=None, disk_space_tolerance=0.01E9, image_export_format=None, save_filename_separator=None, tile_dims=None, save_tiff_references=None, _stacklevel_increment=0, **kwargs ) -> Path: """Save data as bmps and metadata as json. You may pass additional data to be serialized as a dictionary to experiment_data. Returns the saved directory as a Path object. Parameters ---------- mcam_dataset: The mcam_dataset you wish to save. directory: There directory where the data should be saved. metadata_filename: The filename within the directory where the metadata stored as a `json` file should be saved. If set to `None`, no metadata file will be saved. include_timestamp: If set to `True`, this will append a timestamp to the provided directory. If set to false, no timestamp will be appended to the directory possibly overwriting any files currently within the previous directory. imagename_format: Passed to `save_one_as_tif` or 'save_one_as_bmp'. image_mode: 'gray', 'rgb', 'rgba', 'bayer', 'bggr', 'rggb', 'grbg', gbrb' or None If ``mcam_dataset`` is not an ``xarray.DataArray``, then the ``image_mode`` can be used to override the code that we use to guess the dimensions of your image. disk_space_tolerance: float The difference between the free disk space in the given directory and the dataset to be save in bytes. If the directory has fewer free bytes than this tolerance plus the dataset size an error will be raised. save_filename_separator: string If not None, this will be used to create a flat structure using the provided separator, the toplevel directory name as a prefix, and the ``imagename_format``. For example, if the separator is '_', then the output structure will be ``{directory.parent}/{directory.name}_{imagename_format}``. Returns ------- directory: The name of the folder where the data has been exported. """ # ensure mcam_dataset is actually a xr.Dataset if not isinstance(mcam_dataset, xr.Dataset): warn("mcam_datasets that are not xarray.Datasets are no longer supported." "Please contract Ramona Optics at help@ramonaoptics.com with a screenshot of " "your screen including this error message if you have any questions.", stacklevel=2) if isinstance(mcam_dataset, xr.DataArray): mcam_dataset = mcam_dataset.to_dataset() else: mcam_dataarray = _infer_dims(mcam_dataset, input_mode=image_mode) mcam_dataset = mcam_dataarray.to_dataset() for dim in mcam_dataset.dims: if dim not in mcam_dataset.coords: warn( f'{dim} dimension has not been given coordinates within the dataset. ' 'Please contact Ramona Optics at help@ramonaoptics.com to get this resolved. ' 'Please include the dataset and mention #1644 in your message.', stacklevel=2 + _stacklevel_increment, ) if metadata_filename is None: # There is more metadata here, and therefore we typically have # found that we save the metadata as nc metadata_filename = 'metadata.nc' metadata_filename = metadata_filename.format_map( DatasetFormatMapping(mcam_dataset, {'default_filename': 'metadata'}) ) stack_dims = get_stack_dims(mcam_dataset) if image_export_format is None: image_export_format = 'images' if image_export_format == 'multi_tiledmp4': raise ValueError( "The multi_tiledmp4 export format is no longer supported. " "If you need this format, contact us at help@ramonaoptics.com." ) if image_export_format == 'image_stacks': warn("The image_export_format 'image_stacks' has been deprecated " "in favor of 'images'", stacklevel=2) image_export_format = 'images' if image_export_format == 'flat_image_stacks': if tile_dims is not None: raise ValueError( 'flat_image_stacks is only compatible with tile_dims=None. ' f'Got {tile_dims}' ) tile_dims = stack_dims + ('image_y', 'image_x') if tile_dims is None: tile_dims = _default_tile_dims[image_export_format] tiff_stack_dims = tuple( dim for dim in stack_dims if dim not in tile_dims ) valid_tiff_stack_dims = ('channel', 'frame_number', 'timelapse_index') if not all(dim in valid_tiff_stack_dims for dim in tiff_stack_dims): raise ValueError( f'Got unsupported tiff stack dims {tiff_stack_dims}. ' f'Only excepted tiff stack dims are {valid_tiff_stack_dims}.' ) if image_export_format not in _supported_image_export_formats: raise ValueError( f"Got unsupported image_export_format {image_export_format}. " f"Supported formats are {_supported_image_export_formats}.") expected_compression = _expected_compression_factor[image_export_format] if 'valid_data' in mcam_dataset: fraction_valid = float(mcam_dataset.valid_data.mean()) else: fraction_valid = 1. check_disk_space( directory, 1, mcam_dataset.nbytes * expected_compression * fraction_valid, absolute_tolerance=disk_space_tolerance, ) if imagename_format is None: if is_well_dataset(mcam_dataset): if 'field_id' not in mcam_dataset or mcam_dataset.field_id.max() == 0: default_imagename_format = _default_imagename_format_wells else: default_imagename_format = _default_imagename_format_wells_field_id else: default_imagename_format = _default_imagename_format imagename_format = default_imagename_format[image_export_format] if image_export_format == 'flat_image_stacks' and stack_dims: base_path = Path(imagename_format) extension = base_path.suffix stack_parts = '_'.join(f'{{{dim}}}' for dim in stack_dims) imagename_format = f'cam{{image_y}}_{{image_x}}_{stack_parts}{extension}' if save_filename_separator is not None: prefix = directory.name directory = directory.parent imagename_format = prefix + save_filename_separator + imagename_format metadata_filename = prefix + save_filename_separator + metadata_filename extension = Path(imagename_format).suffix.lower() if 'rgba' in mcam_dataset.images.dims and extension == '.bmp': # Bug in PIL # https://stackoverflow.com/questions/10453604/load-rgba-bitmap-with-pil raise ValueError('We cannot save RGBA as bmp. Try using a different ' 'format') if extension not in _valid_extensions: raise ValueError( f"Unknown Extension {extension} provided in imagename_format") if image_export_format not in _exporters: raise ValueError( f"Unknown export format {image_export_format}." ) if extension not in _valid_extensions_for_image_export_format[image_export_format]: raise ValueError( f"The file extension ({extension}) is not compatible with {image_export_format}" ) if image_export_format == 'images' and extension not in ('.tif', '.tiff'): extra_tile_dims = tuple( dim for dim in mcam_dataset.images.dims if dim not in tile_dims + ('y', 'x', 'rgb', 'rgba') ) if extra_tile_dims: needed_tile_dims = tuple( dim for dim in mcam_dataset.images.dims if dim not in ('y', 'x', 'rgb', 'rgba') ) raise ValueError( "Cannot export images with extra dimensions " f"{extra_tile_dims} as {extension} without specifying tile_dims " f"as {needed_tile_dims}." ) chroma = get_chroma(mcam_dataset) valid_chromas = _valid_chroma_for_image_export_format[image_export_format] if chroma not in valid_chromas: raise ValueError( f"The image chroma is {chroma}, but this export format only supports chromas of " f"{valid_chromas}" ) # Shallow copy so as not to affect the user's original dataset mcam_dataset = mcam_dataset.copy() _ensure_2d_dataset(mcam_dataset) images_dims = mcam_dataset.images.dims extra_tile_dims = tuple( d for d in tile_dims if d not in images_dims ) if extra_tile_dims: raise ValueError( f"Requested image format to be flatted on dimensions {extra_tile_dims}" " which are dimensions that do not exist." ) # Re-order the tile dims so that they match they occur in the # images. # This is important since the dataset of the images will be dropped # upon saving so we need to store this information somewhere tile_dims = tuple( t for t in images_dims if t in tile_dims ) if image_export_format != 'tiledmp4' and ( 'image_x' not in tile_dims or 'image_y' not in tile_dims ): raise ValueError( "tile_dims must contain image_y and image_x" ) if extension.lower() in ['.tif', '.tiff']: # This logic must match the logic in _save_tiff_with_metadata per_tiff_dims = tuple( d for d in images_dims if d not in tile_dims ) # stack_dims_per_image = tuple( # d # for d in per_tiff_dims # if d in stack_dims # ) # if stack_dims_per_image: # mcam_dataset.attrs['image_export_imagej_axes_transposer'] = axes_transposer # else: # mcam_dataset.attrs['image_export_imagej_axes'] = None # mcam_dataset.attrs['image_export_imagej_axes_transposer'] = None per_tiff_dims = _adjust_images_dims_for_z_stack(per_tiff_dims, mcam_dataset) imagej_axes, axes_transposer = _get_imagej_axes(per_tiff_dims) mcam_dataset.attrs['image_export_imagej_axes'] = imagej_axes # This transposer was added in 0.19.200 # and was used to invert the transposition of the axes # However it is likely not entirely needed anymore # https://gitlab.com/ramonaoptics/mcam/python-owl/-/commit/d6acf7d2aacd05dc4780b3e91633256e4da3f227 # but we will keep in to allow limited datasets saved with # 0.19.3XX to be loaded with these older image loaders. # We should remove this in 0.20.0. mcam_dataset.attrs['image_export_imagej_axes_transposer'] = axes_transposer mcam_dataset.attrs['imagename_format'] = imagename_format mcam_dataset.attrs['images_dims'] = images_dims mcam_dataset.attrs['image_export_format'] = image_export_format mcam_dataset.attrs['image_export_tile_dims'] = tile_dims # 2023/07/12 # In the past, we sometimes pre-applied an exif transformation in the case of # multi-mp4 datasets. We have removed this "special case" and now always # apply 1 # Always apply the default exif transformation to our dataset, that is "1" # we may allow this to change in the future mcam_dataset.attrs['preapplied_exif_transform'] = 1 directory = Path(directory) if include_timestamp: directory = (directory.parent / (directory.name + '_' + make_timestamp())) export = _exporters[image_export_format] return export( mcam_dataset, directory=directory, metadata_filename=metadata_filename, mode=mode, tile_dims=tile_dims, save_tiff_references=save_tiff_references, **kwargs ) def _infer_dims(data, input_mode): # bayer data if len(data.shape) == 4: if input_mode in valid_bayer_patterns: coords = {'bayer_pattern': input_mode} input_mode = 'bayer' elif input_mode == 'bayer': coords = {'bayer_pattern': ['bggr']} else: # grayscale coords = {} dims = ['image_y', 'image_x', 'y', 'x'] elif len(data.shape) == 5: if input_mode is None: if data.shape[4] == 3: input_mode = 'rgb' elif data.shape[4] == 4: input_mode = 'rgba' if input_mode == 'rgb': dims = ['image_y', 'image_x', 'y', 'x', 'rgb'] coords = {'rgb': ['r', 'g', 'b']} elif input_mode == 'rgba': dims = ['image_y', 'image_x', 'y', 'x', 'rgba'] coords = {'rgba': ['r', 'g', 'b', 'a']} else: raise ValueError("Unknown image input_mode") mcam_dataarray = new_dataset( N_cameras=data.shape[:2], image_shape=data.shape[2:], dtype=data.dtype, dims=dims, coords=coords, array=data )['images'] return mcam_dataarray def load( directory, *, metadata_filename=None, delayed=True, progress=True, scheduler='threading', update=True, tqdm=None, engine='ramona', _stacklevel_increment=0, **_kwargs, ): if not update: # Mark -- 2025/10 # Our data loaders have gotten more elaborate and are requiring # the full dataset to be updated # having the split "updates" for the metadata and data itself # is proving to be unmaintainable # Our updating strategy has gotten more stable and thus we can continue # With confidence # We can revisit this situation in the future # But it has been causing us trouble when trying to load # multi-channel RGB images which are not considered by the TIFF standard # https://gitlab.com/ramonaoptics/mcam/python-owl/-/merge_requests/8223 # Particularly look at results about the inability to load data from # _update_pre_0_18_253_z_stage_to_frame_number # This data was acquired with early versions of the Kestrel0101 in 2023/04 warn( "The 'update=False' can cause unpredictable behavior and is only " "recommended for debugging purposes. It may not load your data in any " "meaningful way. If you require this functionality please contact " "Ramona Optics with information about your usecase.", stacklevel=2 + _stacklevel_increment, ) metadata = _load_metadata(directory, metadata_filename, delayed=delayed, engine=engine) if update: metadata = update_dataset(metadata) mcam_dataset = _load_directory( directory, metadata, update=update, tqdm=tqdm, _stacklevel_increment=_stacklevel_increment + 1 ) if delayed is None: # Automatically load datasets that are less than 1 GB # otherwise, just delay them delayed = mcam_dataset.nbytes > 1E9 if not delayed: if progress: # Allow the user to pass us their context manager if not hasattr(progress, '__enter__'): progress = ProgressBar() else: progress = nullcontext() with progress: _mcam_dataset = mcam_dataset mcam_dataset = mcam_dataset.compute(scheduler=scheduler) _mcam_dataset.close() return mcam_dataset def _multi_file_closer(closers): for closer in closers: if closer is not None: closer() def _load_directory(directory, metadata, update, *, tqdm=None, _stacklevel_increment=0): # Really try to avoid using xarray as long as possible, indexing in # it is really slow.... attrs = metadata.attrs sizes = metadata.sizes images_dims = attrs['images_dims'] N_cameras = (sizes['image_y'], sizes['image_x']) imagename_format = attrs['imagename_format'] image_export_format = attrs['image_export_format'] if image_export_format == 'multi_tiledmp4': raise ValueError( "Loading datasets in multi_tiledmp4 format is no longer supported. " "If you need this format, contact us at help@ramonaoptics.com." ) preapplied_exif_transform = attrs['preapplied_exif_transform'] tile_dims = tuple(attrs['image_export_tile_dims']) missing_dims = tuple(d for d in images_dims if d not in metadata.dims) if missing_dims: warn( "The dataset is missing some dimensions, this may cause unexpected " f"behavior when loading the dataset. Missing dimensions: {missing_dims}.", stacklevel=2 + _stacklevel_increment, ) images_dims = tuple(d for d in images_dims if d not in missing_dims) single_image_dims = tuple( d for d in images_dims if d in ('y', 'x', 'rgb', 'rgba') ) extension = Path(imagename_format).suffix.lower() inverse_exif_rotation = { 1: 1, 2: 2, 3: 3, 4: 4, 5: 7, 6: 8, 7: 5, 8: 6, } frame_exif_rotate = inverse_exif_rotation[preapplied_exif_transform] tile_shape = tuple( metadata.sizes[dim] for dim in tile_dims ) image_names = get_image_names( metadata, tile_dims, tile_shape, imagename_format, ) if ( image_export_format in ("images", "flat_image_stacks") and extension in [".j2c", ".jph"] ): filenames = np.empty_like(image_names) for i in ndrange(tile_shape): filenames[i] = directory / image_names[i] valid_data = np.zeros(filenames.shape, dtype='bool') from ..array._j2c_backend import J2CArray final_stack = J2CArray( filenames, valid_data=valid_data, tqdm=tqdm, frame_exif_rotate=frame_exif_rotate, ) if update or 'valid_data' in metadata: metadata['valid_data'] = (tile_dims, valid_data) elif ( image_export_format in ("images", "flat_image_stacks") and extension not in [".tif", ".tiff"] ): filenames = np.empty_like(image_names) for i in ndrange(tile_shape): filenames[i] = directory / image_names[i] valid_data = np.zeros(filenames.shape, dtype='bool') final_stack = FlatImageioArray( filenames, valid_data=valid_data, tqdm=tqdm, frame_exif_rotate=frame_exif_rotate, ) if update or 'valid_data' in metadata: metadata['valid_data'] = (tile_dims, valid_data) elif ( image_export_format in ("images", "flat_image_stacks") # and extension in [".tif", ".tiff"] ): filenames = np.empty_like(image_names) for i in ndrange(tile_shape): filenames[i] = directory / image_names[i] ref_ds = get_ref_ds(metadata) if ref_ds is not None: for var_name in list(ref_ds.data_vars): ref_dims = list(ref_ds[var_name].dims) reorder = [ ref_dims.index(d) for d in tile_dims if d in ref_dims ] + [ ref_dims.index(d) for d in ref_dims if d not in tile_dims ] arr = np.transpose(np.asarray(ref_ds[var_name].values), reorder) for axis, d in enumerate(tile_dims): if d not in ref_dims: arr = np.expand_dims(arr, axis=axis) ref_ds[var_name] = ( tile_dims + tuple(dim for dim in ref_ds[var_name].dims if dim not in tile_dims), arr, ) if 'valid_data' in metadata.data_vars: valid_data_arr = np.asarray(metadata['valid_data'].values, dtype='bool') valid_data_dims = list(metadata['valid_data'].dims) if valid_data_arr.shape != tile_shape: extra_axes = tuple( i for i, d in enumerate(valid_data_dims) if d not in tile_dims ) if extra_axes: valid_data_arr = np.any(valid_data_arr, axis=extra_axes) valid_data_dims = [d for d in valid_data_dims if d in tile_dims] reorder = [ valid_data_dims.index(d) for d in tile_dims if d in valid_data_dims ] valid_data_ordered = np.transpose(valid_data_arr, reorder) for axis, d in enumerate(tile_dims): if d not in valid_data_dims: valid_data_ordered = np.expand_dims(valid_data_ordered, axis=axis) valid_data = np.broadcast_to(valid_data_ordered, tile_shape).copy() else: valid_data = valid_data_arr.copy() else: valid_data = np.ones(filenames.shape, dtype='bool') stack_dims_per_image = tuple( d for d in images_dims if d not in tile_dims + ('y', 'x', 'rgb', 'rgba') ) image_shape = tuple( metadata.sizes[d] for d in images_dims if d in ('y', 'x', 'rgb', 'rgba') ) if not stack_dims_per_image: axes_transposer = None hyperstack_dims = () else: # vvvvvvvvvvvvvvvvvvvvv Mark 2025/11/10 vvvvvvvvvvvvvvvvvvvvvvvvvvvvv # We manually recreate the variable attrs['image_export_imagej_axes_transposer'] # In case there are missing dimensions in the dataset. per_tiff_dims = stack_dims_per_image + ('y', 'x') if 'rgba' in images_dims: per_tiff_dims += ('rgba',) if 'rgb' in images_dims: per_tiff_dims += ('rgb',) per_tiff_dims = _adjust_images_dims_for_z_stack(per_tiff_dims, metadata) _, axes_transposer = _get_imagej_axes(per_tiff_dims) # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # Compte the inverse of the transposer axes_transposer = tuple( axes_transposer.index(t) for t in range(len(axes_transposer)) ) _stack_dims_per_image = _adjust_images_dims_for_z_stack( stack_dims_per_image, metadata, ) _, hyperstack_transposer = _get_imagej_axes(_stack_dims_per_image) hyperstack_dims = tuple( stack_dims_per_image[t] for t in hyperstack_transposer ) hyperstack_shape = tuple( metadata.sizes[d] for d in hyperstack_dims ) tiff_array = TiledMultiPageTiffArray( filenames, # Tiffs have a habit of just collapsing the channel dimension # https://github.com/cgohlke/tifffile/issues/241 # So we need to expand the axes if it exists. hyperstack_shape=hyperstack_shape, tiff_transpose=axes_transposer, valid_data=valid_data, frame_exif_rotate=frame_exif_rotate, tqdm=tqdm, ref_ds=ref_ds, ref_image_shape=image_shape, ) if update or 'valid_data' in metadata: metadata['valid_data'] = (tile_dims, valid_data) # At this stage, we have our stack dims at the front of the TIFF # array and the images at the back, we must retranspose them tiff_dims = tile_dims + stack_dims_per_image tiff_dims += tuple(images_dims[len(tiff_dims):]) transpose = tuple( tiff_dims.index(d) for d in images_dims ) final_stack = tiff_array.transpose(transpose) elif image_export_format == "multi_mp4": filenames = np.empty_like(image_names) for i in ndrange(tile_shape): filenames[i] = directory / image_names[i] if single_image_dims[-1] == "rgb": output_pixel_format = 'rgb24' else: output_pixel_format = 'gray8' valid_data = np.zeros(filenames.shape, dtype='bool') video_array = TiledVideoArray( filenames, output_pixel_format=output_pixel_format, frame_exif_rotate=frame_exif_rotate, valid_data=valid_data, tqdm=tqdm, ) if update or 'valid_data' in metadata: metadata['valid_data'] = (tile_dims, valid_data) # TiledVideoArray always lays its axes out as # (frame, *tile_dims, *single_image_dims). For a plain acquisition # that already matches images_dims, but a timelapse export keeps each # timepoint in its own mp4, so timelapse_index is a tile dimension and # images_dims places it *before* frame_number. Retranspose the per-mp4 # frame axis into the position the metadata expects, mirroring how the # tiff branch above retransposes its stack dimensions. stack_dims_per_video = tuple( d for d in images_dims if d not in tile_dims + ('y', 'x', 'rgb', 'rgba') ) video_dims = stack_dims_per_video + tile_dims + single_image_dims transpose = tuple( video_dims.index(d) for d in images_dims ) final_stack = video_array.transpose(transpose) elif image_export_format == "tiledmp4": from ._mp4 import _load_tiledmp4 # TiledMP4 lists all data as valid valid_data = np.ones(N_cameras, dtype='bool') if update or 'valid_data' in metadata: metadata['valid_data'] = (('image_y', 'image_x'), valid_data) final_stack = _load_tiledmp4(directory, metadata) else: raise ValueError( f"Unknown value for exported format ({image_export_format})") # Don't pass the whole dataset as coordinates, just pass the coordinates # of the dataset # https://github.com/pydata/xarray/pull/2344#issuecomment-410564717 images = xr.DataArray( final_stack, dims=images_dims, ) mcam_dataset = metadata.assign({'images': images}) # mcam_dataset.images.set_close(final_stack.close) mcam_dataset.set_close( partial(_multi_file_closer, [mcam_dataset._close, final_stack.close]) ) return mcam_dataset def _adjust_image_dimensions_for_imagej(image, dims): dtype = image.dtype imagej_axes, axes_transposer = _get_imagej_axes(dims) # No transposing needed, just return the same array # When comparing, ensure you compare both as tuples # to avoid false due to changes in collection type if tuple(axes_transposer) == tuple(range(image.ndim)): return image, image.shape, dtype, imagej_axes # If we return a numpy array that is simply transposed, # it breaks our fast direct writer. read on to learn why # To avoid this, we return an "iterator" of the data # This helps tifffile be more "lazy" about how it it treats the # data iterating chunk by chunk. # # In fact that we "transpose" will cause a copy of the data # within the tiff writer, this is really problematic # since they don't use our aligned allocator # and thus it will cause the data to be copied to an unaligned buffer # Ultimately, it will cause the cache to fill for these big stacks that # we may be trying to write fast. # The specific code we are trying to avoid is # https://github.com/cgohlke/tifffile/blob/e61dfaeda97c2a337de50c5ffb4abd21b639363d/tifffile/tifffile.py#L2814 # where the following is done # dataarray = dataarray.reshape("flattened shape") # And thus this really breaks things transposed = np.transpose( image, axes=axes_transposer ) if 'rgb' in dims or 'rgba' in dims: d = 3 else: d = 2 data = [ transposed[i] for i in ndrange(transposed.shape[:-d]) ] return data, transposed.shape, dtype, imagej_axes # Jed 20231013 - currently we only target saving in imagej dimension order. # In the future we may decide to make this more user controlled? def _get_imagej_axes(dims): if 'timelapse_index' in dims and 'frame_number' in dims: converter = { 'frame_number': 'Z', 'timelapse_index': 'T', } else: converter = { 'frame_number': 'T', 'timelapse_index': 'T', } converter = converter | { 'channel': 'C', 'z_stage': 'Z', 'y': 'Y', 'x': 'X', 'rgb': 'S', 'rgba': 'S', } input_axes = '' for dim in dims: input_axes += converter[dim] imagej_axes = '' imagej_all_axes = 'TZCYXS' for a in imagej_all_axes: if a in input_axes: imagej_axes += a axes_transposer = tuple( imagej_axes.index(a) for a in input_axes ) return imagej_axes, axes_transposer def _adjust_images_dims_for_z_stack(images_dims, metadata): if _tiff_is_z_stack(images_dims, metadata): images_dims = tuple( dim if dim != 'frame_number' else 'z_stage' for dim in images_dims ) return images_dims def _tiff_is_z_stack(images_dims, metadata): return ( 'z_stage' in metadata and 'frame_number' in metadata.z_stage.dims and 'z_stage' not in images_dims ) def _get_pixel_width_from_metadata(metadata): if metadata is None: return 0.0 if 'pixel_width' not in metadata: return 0.0 pixel_width = metadata['pixel_width'] if hasattr(pixel_width, 'values'): pixel_width = pixel_width.values pixel_width = np.asarray(pixel_width) if pixel_width.shape == (): value = float(pixel_width) else: value = float(pixel_width.mean()) if not math.isfinite(value) or value <= 0: return 0.0 return value def _get_timestamp_from_metadata(metadata): if metadata is None: return None if 'software_timestamp' not in metadata: return None timestamp = metadata.software_timestamp if hasattr(timestamp, 'values'): timestamp = timestamp.values if np.asarray(timestamp).ndim > 0: timestamp = np.asarray(timestamp).flat[0] return timestamp def _get_serial_number_from_metadata(metadata): if metadata is None: return None if 'serial_number' not in metadata: return None return str(np.asarray(metadata['serial_number'])) def _get_description_from_metadata(metadata): if metadata is None: return None if 'description' not in metadata.attrs: return None return str(metadata.attrs['description']) def _get_exif_orientation_from_metadata(metadata): if metadata is None: return None if 'exif_orientation' not in metadata: return None value = int(np.asarray(metadata['exif_orientation'])) if value < 1 or value > 8: return None return value def _save_png_with_metadata( image, filename, metadata, include_timestamp=True, mode='w', ): from PIL.PngImagePlugin import PngInfo filename = Path(filename) if filename.suffix == '': filename = filename.with_suffix('.png') if include_timestamp: filename = filename.parent / (filename.stem + '_' + make_timestamp() + filename.suffix) if mode == 'x' and filename.exists(): raise FileExistsError(f"File already exists: {filename}") pixel_width = _get_pixel_width_from_metadata(metadata) timestamp = _get_timestamp_from_metadata(metadata) serial_number = _get_serial_number_from_metadata(metadata) description = _get_description_from_metadata(metadata) exif_orientation = _get_exif_orientation_from_metadata(metadata) if pixel_width > 0: dpi = int(round(0.0254 / pixel_width)) else: dpi = None png_info = PngInfo() png_info.add_itxt("Software", "Ramona Optics MCAM") if timestamp is not None: iso_timestamp = np.datetime_as_string(timestamp, unit='s') + 'Z' png_info.add_itxt("Timestamp", iso_timestamp) if serial_number is not None: png_info.add_itxt("SerialNumber", serial_number) if description is not None: png_info.add_itxt("Description", description) if exif_orientation is not None: png_info.add_itxt("Orientation", str(exif_orientation)) save_kwargs = { 'pnginfo': png_info, 'compress_level': 3, } if dpi is not None: save_kwargs['dpi'] = (dpi, dpi) imageio.imwrite(filename, image, plugin='pillow', **save_kwargs) return filename def _save_jpeg_with_metadata( image, filename, metadata, include_timestamp=True, mode='w', quality=95, ): from PIL.Image import Exif filename = Path(filename) if filename.suffix == '': filename = filename.with_suffix('.jpg') if include_timestamp: filename = filename.parent / (filename.stem + '_' + make_timestamp() + filename.suffix) if mode == 'x' and filename.exists(): raise FileExistsError(f"File already exists: {filename}") pixel_width = _get_pixel_width_from_metadata(metadata) timestamp = _get_timestamp_from_metadata(metadata) serial_number = _get_serial_number_from_metadata(metadata) description = _get_description_from_metadata(metadata) exif_orientation = _get_exif_orientation_from_metadata(metadata) # JPEG doesn't support alpha, strip it if image.ndim == 3 and image.shape[-1] == 4: image = image[..., :3] # Build EXIF data using numeric tag IDs # IFD0 tags: Software(305), Make(271), ImageDescription(270), # XResolution(282), YResolution(283), ResolutionUnit(296), # Orientation(274) # EXIF IFD (0x8769): DateTimeOriginal(36867), OffsetTimeOriginal(36881), # BodySerialNumber(42033) exif = Exif() exif[305] = "Ramona Optics MCAM" # Software exif[271] = "Ramona Optics" # Make if exif_orientation is not None: exif[274] = exif_orientation # Orientation if description is not None: exif[270] = description # ImageDescription if pixel_width > 0: resolution = 0.01 / pixel_width # pixels per centimeter exif[282] = resolution # XResolution exif[283] = resolution # YResolution exif[296] = 3 # ResolutionUnit (centimeters) exif_ifd = {} if timestamp is not None: exif_timestamp = ( np.datetime_as_string(timestamp, unit='s').replace('T', ' ').replace('-', ':') ) exif_ifd[36867] = exif_timestamp # DateTimeOriginal exif_ifd[36881] = "+00:00" # OffsetTimeOriginal if serial_number is not None: exif_ifd[42033] = serial_number # BodySerialNumber if exif_ifd: exif[0x8769] = exif_ifd # EXIF IFD imageio.imwrite(filename, image, plugin='pillow', quality=quality, exif=exif) return filename @kwarg_name_change('0.20', {'compressionargs': 'compression_args'}) def _save_tiff_with_metadata( image, filename, metadata, include_timestamp=True, mode='w', compression=None, compression_args=None, predictor=None, pyramid_levels=None, pyramid_downsample=None, rowsperstrip=None, **kwargs, ): """Save tiff with some of the given metadata Currently the resolution and orientation of the image as well as the date time of saving and the Ramona Optics software is included in the tiff file as metadata. Returns the saved filename as a Path object. Parameters ---------- image: numpy.array Array to be saved as a tiff file. filename: pathlike The filename where the data should be saved. metadata: xarray.dataset The metadata containing coords of the the desired calibration data to be added to the image. If `pixel width` is not present in the metadata then both the x and y resolution of the saved tiff file is set to '0 pixels per none units'. If `exif_orientation` is not present in the metadata then orientation is not set. include_timestamp: If set to `True`, this will append a timestamp to the provided directory. If set to false, no timestamp will be appended to the directory possibly overwriting any files currently within the previous directory. mode: str The mode to open the file in. Valid options are ``'w'`` (write), or ```x``` exclusive creation, failing if the file already exists. compression: str The compression to use when saving the tiff. Passed to tifffile.write. compression_args: dict[str, Any] The compression arguments to use when saving the tiff. Passed to tifffile.write. predictor: tifffile.PREDICTOR A predictor to applied with the compression algorithm. Not supported when the compression is unspecified. """ images_dims = metadata.images.dims stack_dims = get_stack_dims(metadata) if len(stack_dims) > 3: raise RuntimeError("We do not support more than 3 stack dimensions") if compression is None: # When unspecified, default to zstd (lossless) compression. # Pass an explicit 'none' to disable compression. compression = 'zstd' # This may become useful as an option in the future, but tifffile # just chokes if you specify this. if compression in ['None', 'NONE', 'none']: compression = None if predictor == 'horizontal': predictor = tifffile.PREDICTOR.HORIZONTAL if compression == 'zlib' and compression_args is None and predictor is None: # Some sensible default values for the 'zlib' compressor # We found these through some benchmarking and this gives a good mix # of compression ratio and compatibility with ImageJ. compression_args = dict(level=6) if image.dtype.kind != 'f': # Horizontal predictor is the best predictor for most images # However, it is not supported for floating point images predictor = tifffile.PREDICTOR.HORIZONTAL if compression == 'zstd' and compression_args is None and predictor is None: # Some sensible default values for the 'zstd' compressor # We found these through some benchmarking and this gives a good mix # of compression ratio and compatibility with ImageJ. compression_args = dict() if image.dtype.kind != 'f': # Horizontal predictor is the best predictor for most images # However, it is not supported for floating point images predictor = tifffile.PREDICTOR.HORIZONTAL if compression in ['jpeg', 'jpegxl', 'jpeg2000',] and predictor is not None: raise ValueError(f"{compression} does not support any predictor. got {predictor}") if compression == 'jpeg' and compression_args is None: # Default level for JPEG compression # Level ranges from 1 (worst) to 100 (best) compression_args = dict(level=85) if compression == 'jpegxl' and compression_args is None: # Default parameters for JPEG XL compression # Use high quality lossy compression by default (level 90, effort 1) # Level 90 provides excellent quality while still being lossy # I've found that level 90 on jpegxl gave similar compression ratios to level 85 on jpeg # and effort 1 is the fastest compression # Efforts lower than 5 would be faster, but honestly, they would essentially fall back # to jpeg quality, which was kinda bad and blocky. # See the following post # https://www.reddit.com/r/jpegxl/comments/10157np/comment/j2lr7nr/ # For a nice summary compression_args = dict(level=90, effort=5) if compression == 'jpeg2000-nr': # jpeg2000-nr is Ramona Short-hand for "jpeg2000 non-reversible" # Which turns on compression by default if compression_args is None: # By default jpeg2000 would use # 0.0039 for non-reversible (nr) compression # However, this is too argressive, I've found that # qstep 0.012 gives "acceptable results" for most images # It can probably be adjusted further as we learn to use # jpeg2000 compression_args = dict(qstep=0.012) if compression_args.get('reversible', False): raise ValueError("jpeg2000-nr requires reversible=False") compression_args['reversible'] = False compression = 'jpeg2000' if compression == 'jpeg2000' and compression_args is None: # Default parameters for JPEG2000 compression # Use lossless compression by default (reversible=True) # Users can specify reversible=False and qstep for lossy compression compression_args = dict(reversible=True) if compression == 'jpegxl' and compression_args is not None: # For JPEG XL, quality/level 100 means "distance 0", i.e. no lossy # compression. libjxl >= 0.12 only accepts distance 0 when lossless # encoding is requested explicitly, so route level >= 100 to the # dedicated lossless flag that imagecodecs understands. level = compression_args.get('level') if level is not None and level >= 100: compression_args = { key: value for key, value in compression_args.items() if key != 'level' } compression_args.setdefault('lossless', True) # This controls details of the zlib compression algorithm # It stores the number of given rows contiguously # It helps to set it to the entire y extent max_bytes_per_strip = 2 * 1024 ** 3 if rowsperstrip is not None: # allow us to experiment by specifying what rowsperstrip could be pass elif 'y' in metadata: if 'rgb' in metadata: colors_per_pixel = 3 elif 'rgba' in metadata: colors_per_pixel = 4 else: colors_per_pixel = 1 # By default, just compress the entire image all at once rowsperstrip = min( metadata.sizes['y'], max_bytes_per_strip // (metadata.sizes['x'] * image.dtype.itemsize * colors_per_pixel) ) elif 'rgb' in metadata: rowsperstrip = min( image.shape[-3], max_bytes_per_strip // (image.shape[-2] * image.dtype.itemsize * 3) ) elif 'rgba' in metadata: rowsperstrip = min( image.shape[-3], max_bytes_per_strip // (image.shape[-2] * image.dtype.itemsize * 4) ) else: rowsperstrip = min( image.shape[-2], max_bytes_per_strip // (image.shape[-2] * image.dtype.itemsize) ) # https://github.com/cgohlke/tifffile/blob/v2024.9.20/tifffile/tifffile.py#L2520-L2526 imagej = 'rgb' not in images_dims and 'rgba' not in images_dims if imagej: illumination_mode = np.atleast_1d(np.asarray(metadata.get('illumination_mode', []))) # Compression seems to negatively impact imageJ's ability to read the # LUT and Label metadata # https://gitlab.com/ramonaoptics/mcam/python-owl/-/merge_requests/6252#note_2171968748 luts = generate_luts(illumination_mode) labels = [ str(visual_illumination_modes.get(mode, mode)) for mode in np.asarray(illumination_mode) ] imagej_mode = 'composite' if len(illumination_mode) > 1 else None fps = float(metadata.frame_rate_setpoint) if 'frame_rate_setpoint' in metadata else None xorigin = int(metadata.get('x', np.array(0)).min()) yorigin = int(metadata.get('y', np.array(0)).min()) else: luts = None labels = None imagej_mode = None fps = None xorigin = None yorigin = None filename = Path(filename) if filename.suffix == '': filename = filename.with_suffix('.tif') if include_timestamp: filename = (filename.parent / (filename.name + '_' + make_timestamp() + filename.suffix)) extratags = [] if 'exif_orientation' in metadata: exif_orientation = np.asarray(metadata['exif_orientation']) if exif_orientation.shape == (): orientation = int(exif_orientation) # Stacked datasets may have different exif based on each image elif np.all(exif_orientation == exif_orientation.flat[0]): orientation = int(exif_orientation.flat[0]) else: orientation = None if orientation is not None: extratags.append(( tifffile.TIFF.TAGS['Orientation'], tifffile.DATATYPE.SHORT, 1, orientation, False )) else: extratags = [] pixel_width = _get_pixel_width_from_metadata(metadata) if pixel_width > 0.: # Try to detect if we have a z stack if ( 'z_stage' in metadata and metadata.z_stage.sizes.get('frame_number', 0) > 1 ): imagej_spacing = float(metadata.z_stage.diff('frame_number').mean()) else: imagej_spacing = None try: bin_mode = get_bin_mode(metadata) except AttributeError: # "graceful" fallback if y (and x) coordinates are not defined bin_mode = 1 image_pixel_width = pixel_width * bin_mode # I just don't think that units other than centimeter or inch are allowed # https://github.com/cgohlke/tifffile/commit/5a11183ec80a6f7073f6c595e9b3769892a3387e # https://github.com/cgohlke/tifffile/blob/master/tifffile/tifffile.py#L9151 # While some indications say that mm or um should be supported # libtiff just gives a huge warning that 5 is unknown imagej_unit = 'cm' imagej_unit_factor = 1E-2 resolutionunit = 'CENTIMETER' pixels_per_unit = float(imagej_unit_factor / (image_pixel_width)) resolution = (pixels_per_unit, pixels_per_unit) if imagej_spacing is not None: imagej_spacing = imagej_spacing / imagej_unit_factor else: resolution = None resolutionunit = None imagej_unit = None imagej_spacing = None description = _get_description_from_metadata(metadata) # Exif Tags can be complicated. Refer to # https://exiftool.org/TagNames/EXIF.html # https://www.loc.gov/preservation/digital/formats/content/tiff_tags.shtml # https://www.awaresystems.be/imaging/tiff/tifftags/privateifd/exif.html """ 0. code (int): Tag Id. 1. dtype (:py:class:`DATATYPE`): Data type of items in `value`. 2. count (int): Number of data values. Not used for string or bytes values. 3. value (Sequence[Any]): `count` values compatible with `dtype`. Bytes must contain count values of dtype packed as binary data. 4. writeonce (bool): If *True*, write tag to first page of a series only. """ timestamp = _get_timestamp_from_metadata(metadata) if timestamp is not None and len(stack_dims) == 0: # We specify things in UTC time, we don't have an offset date_time_original = np.datetime_as_string(timestamp, unit='s').replace('T', ' ') extratags.append(( tifffile.TIFF.TAGS['DateTimeOriginal'], tifffile.DATATYPE.ASCII, 1, date_time_original, False, )) extratags.append(( tifffile.TIFF.TAGS['OffsetTimeOriginal'], tifffile.DATATYPE.ASCII, 1, "Etc/UTC", False, )) local_timezone = get_localzone_name() if local_timezone is not None: extratags.append(( tifffile.TIFF.TAGS['OffsetTime'], tifffile.DATATYPE.ASCII, 1, local_timezone, False, )) serial_number = _get_serial_number_from_metadata(metadata) if serial_number is not None: extratags.append(( tifffile.TIFF.TAGS['BodySerialNumber'], tifffile.DATATYPE.ASCII, 1, serial_number, False, )) if 'rgb' in images_dims: if images_dims[-1] != 'rgb': raise ValueError( "We do not support writing images where RGB is not the last dimension") photometric = 'rgb' planarconfig = 'contig' extrasamples = None elif 'rgba' in images_dims: if images_dims[-1] != 'rgba': raise ValueError( "We do not support writing images where RGBA is not the last dimension") photometric = 'rgb' planarconfig = 'contig' extrasamples = 'assocalpha' else: photometric = 'minisblack' planarconfig = None extrasamples = None # I believe we can specify metadata['axes'] = 'TYX' or 'CYX' or 'ZYX' # but until we have a better usecase for this, I don't want to add # complexity # Without specifying it, tifffile infers this as 'QYX' # Necessary to save tiff stack >=4GB # we enable it for images that are larger than 3.5 GB to ensure that smaller # stacks continue to be supported by software that does not support bigtiffs bigtiff = image.nbytes > (3 * 1024 + 512) * (1024 ** 2) if ( compression == 'jpeg2000' and pyramid_downsample is None and pyramid_levels in (None, 'auto') ): if metadata.sizes['y'] % 3 == 0 and metadata.sizes['x'] % 3 == 0: pyramid_downsample = 3 pyramid_levels = 2 else: # jpeg2000 has native pyramids of ::2 pyramid_levels = 1 if pyramid_downsample is None: pyramid_downsample = 4 if pyramid_levels is None: # When unspecified, default to a 3-level pyramid. pyramid_levels = 3 if pyramid_levels == 'auto': _image_height = metadata.sizes['y'] _image_width = metadata.sizes['x'] pyramid_levels = 1 while ((_image_height % pyramid_downsample) == 0 and (_image_width % pyramid_downsample) == 0): # let's limit the smallest level to no smaller than 64x64 if _image_height < 64 or _image_width < 64: break if pyramid_levels >= 3: break _image_height //= pyramid_downsample _image_width //= pyramid_downsample pyramid_levels += 1 if pyramid_levels > 1: subifds = pyramid_levels - 1 def pyramid_reduce(image, has_rgb, downsample): # Also compute the output shape since it isn't obvious just from # the input size of the image due to rounding if isinstance(image, list): image = [ img[::downsample, ::downsample] for img in image ] elif has_rgb: # ND-Array with RGB image = image[..., ::downsample, ::downsample, :] else: # ND-Array monochrome image = image[..., ::downsample, ::downsample] return image else: subifds = None tiff_metadata = {} images_dims = _adjust_images_dims_for_z_stack(images_dims, metadata) image, shape, dtype, imagej_axes = _adjust_image_dimensions_for_imagej(image, images_dims) tiff_metadata['axes'] = imagej_axes if luts is not None: tiff_metadata['LUTs'] = luts if labels is not None: tiff_metadata['Labels'] = labels if imagej_mode is not None: tiff_metadata['mode'] = imagej_mode if fps is not None: tiff_metadata['fps'] = fps if xorigin is not None: tiff_metadata['xorigin'] = xorigin if yorigin is not None: tiff_metadata['yorigin'] = yorigin # It seems that since 2025 we need to start to add the 'unit' # to the metadata otherwise ImageJ doesn't display the units # correctly if imagej and imagej_unit is not None: tiff_metadata['unit'] = imagej_unit if imagej and imagej_spacing is not None: tiff_metadata['spacing'] = imagej_spacing has_rgb = 'rgb' in images_dims or 'rgba' in images_dims with FileHandle(filename, mode=mode) as f: with tifffile.TiffWriter(f, bigtiff=bigtiff, imagej=imagej) as tif_w: # Tiffile is really stubborn, and insists on re-openning the # file handle with their own class so we can't change how data # is written old_fh = tif_w._fh tif_w._fh = f tif_w.write( # For lists, provide it as an iterator # So that we trigger the lazy loading path in tifffile image if not isinstance(image, list) else iter(image), # Shape and dtype must be specified in the case # that image is a generator shape=shape, dtype=dtype, description=description, resolution=resolution, resolutionunit=resolutionunit, # Ensure they are aligned to 1 page since most of our images # are pretty huge and could benefit from the ability to skip # the cache align=4096, software='Ramona Optics MCAM', extratags=extratags, datetime=True, photometric=photometric, planarconfig=planarconfig, rowsperstrip=rowsperstrip, compression=compression, # tifffile uses no underscore for this # but we decided to add an underscore # https://gitlab.com/ramonaoptics/mcam/python-owl/-/merge_requests/8788#note_2985239520 compressionargs=compression_args, predictor=predictor, extrasamples=extrasamples, maxworkers=1, metadata=tiff_metadata, subifds=subifds, ) if subifds is not None: tif_w._imagej = False for i in range(1, pyramid_levels): image = pyramid_reduce( image, has_rgb=has_rgb, downsample=pyramid_downsample, ) if has_rgb: shape = shape[:-3] + ( (shape[-3] + pyramid_downsample - 1) // pyramid_downsample, (shape[-2] + pyramid_downsample - 1) // pyramid_downsample, shape[-1] ) else: shape = shape[:-2] + ( (shape[-2] + pyramid_downsample - 1) // pyramid_downsample, (shape[-1] + pyramid_downsample - 1) // pyramid_downsample, ) if resolution is not None: resolution = (resolution[0] * pyramid_downsample, resolution[1] * pyramid_downsample) if rowsperstrip is not None: rowsperstrip *= pyramid_downsample tif_w.write( # For lists, provide it as an iterator # So that we trigger the lazy loading path in tifffile image if not isinstance(image, list) else iter(image), shape=shape, dtype=dtype, description=description, resolution=resolution, resolutionunit=resolutionunit, align=4096, software='Ramona Optics MCAM', datetime=True, photometric=photometric, planarconfig=planarconfig, rowsperstrip=rowsperstrip, compression=compression, extrasamples=extrasamples, maxworkers=1, ) tif_w._fh = old_fh return filename def get_image_names(metadata, tile_dims, tile_shape, imagename_format): imagename_cache = metadata.get('_images_imagename_cache') if imagename_cache is None: return format_string_array( imagename_format, tile_dims, metadata, ) raw = np.asarray(imagename_cache) if raw.shape != tile_shape: return format_string_array( imagename_format, tile_dims, metadata, ) if raw.ndim == len(tile_dims) + 1 and raw.dtype.kind in ('S', 'U'): image_names = np.empty(tile_shape, dtype=object) for i in ndrange(tile_shape): part = raw[i] if part.dtype.kind == 'S': image_names[i] = part.tobytes().decode( 'utf-8', errors='replace' ).rstrip('\x00') else: image_names[i] = ''.join(part.flat).rstrip('\x00') return image_names return np.asarray(raw, dtype=object) def get_ref_ds(metadata): required = [ '_images_tiff_page_offsets', '_images_tiff_page_nbytes', '_images_tiff_page_compression', '_images_tiff_page_predictor', '_images_tiff_page_dtype', '_images_tiff_page_levels_downsampling', ] if any(v not in metadata.data_vars for v in required): return None return metadata[required] def _add_imagename_cache_to_metadata(metadata): tile_dims = metadata.attrs['image_export_tile_dims'] metadata_for_format = metadata.drop_vars('images', errors='ignore') names = np.asarray( format_string_array( metadata_for_format.attrs['imagename_format'], tile_dims, metadata_for_format, ), dtype=str, ) coords = { d: metadata.coords[d] for d in tile_dims if d in metadata.coords } imagename_ds = xr.Dataset( {'_images_imagename_cache': (tile_dims, names)}, coords=coords, ) metadata = metadata.drop_vars( list(imagename_ds.data_vars), errors="ignore", ) metadata = xr.merge( [metadata, imagename_ds], compat='override', join='right', ) return metadata def _add_tiff_refs_to_metadata( metadata, directory, ): directory = Path(directory) imagename_format = metadata.attrs['imagename_format'] tile_dims = metadata.attrs['image_export_tile_dims'] tile_shape = tuple(metadata.sizes[d] for d in tile_dims) filenames = get_image_names(metadata, tile_dims, tile_shape, imagename_format) # Determine array dimensions from first entry that successfully opens. n_pages = None n_levels = None for i in ndrange(tile_shape): filename = directory / str(filenames[i]) try: with tifffile.TiffFile(filename, mode='r') as tif: pages = len(tif.pages) if (sub_pages := tif.pages[0].pages) is None: levels = 1 else: levels = 1 + len(sub_pages) except Exception: continue if n_pages is None and n_levels is None: n_pages = pages n_levels = levels else: if pages != n_pages or levels != n_levels: raise ValueError( "All TIFF files used for reference export must have the same " f"number of pages and levels. Expected ({n_pages}, {n_levels}), " f"got ({pages}, {levels}) for {filename}." ) if n_pages is None or n_levels is None: raise RuntimeError( "Could not establish the correct number of pages or levels." ) offsets = np.zeros(tile_shape + (n_pages, n_levels), dtype=np.int64) nbytes = np.zeros(tile_shape + (n_pages, n_levels), dtype=np.int64) compression = np.zeros(tile_shape + (n_pages, n_levels), dtype=np.int64) predictor = np.zeros(tile_shape + (n_pages, n_levels), dtype=np.int64) dtype_array = np.zeros(tile_shape + (n_pages, n_levels), dtype='S8') levels_downsampling = np.zeros(tile_shape + (n_pages, n_levels), dtype=np.int64) for i in ndrange(tile_shape): filename = directory / str(filenames[i]) try: with tifffile.TiffFile(filename) as tif: if len(tif.pages) != n_pages: raise ValueError( "All TIFF files used for reference export must have the same " f"number of pages. Expected {n_pages}, got {len(tif.pages)} " f"for {filename}." ) page0 = tif.pages[0] level_downsampling = [1] if page0.pages is not None: for sub in page0.pages: level_downsampling.append(page0.shape[1] // sub.shape[1]) if len(level_downsampling) != n_levels: raise ValueError( "All TIFF files used for reference export must have the same " f"number of levels. Expected {n_levels}, got " f"{len(level_downsampling)} for {filename}." ) dtype_code = np.dtype(page0.dtype).str.encode('ascii') for page_index, page in enumerate(tif.pages): index = i + (page_index, 0) if len(page.dataoffsets) == 1: offsets[index] = int(page.dataoffsets[0]) nbytes[index] = int(page.databytecounts[0]) compression[index] = int(page.compression) predictor[index] = int(page.predictor) dtype_array[index] = dtype_code else: offsets[index] = 0 nbytes[index] = 0 compression[index] = 0 predictor[index] = 0 dtype_array[index] = dtype_code if page.pages is not None: for level_index, sub_page in enumerate(page.pages, start=1): index = i + (page_index, level_index) if len(sub_page.dataoffsets) == 1: offsets[index] = int(sub_page.dataoffsets[0]) nbytes[index] = int(sub_page.databytecounts[0]) compression[index] = int(sub_page.compression) predictor[index] = int(sub_page.predictor) dtype_array[index] = dtype_code else: offsets[index] = 0 nbytes[index] = 0 compression[index] = 0 predictor[index] = 0 dtype_array[index] = dtype_code for level_index, down in enumerate(level_downsampling): levels_downsampling[i + (page_index, level_index)] = int(down) except Exception: # The file might not be found. pass dims = tile_dims + ('_images_tiff_page', '_images_tiff_level') coords = { '_images_tiff_page': np.arange(n_pages), '_images_tiff_level': np.arange(n_levels), } tif_ref_ds = xr.Dataset( { '_images_tiff_page_offsets': (dims, offsets), '_images_tiff_page_nbytes': (dims, nbytes), '_images_tiff_page_compression': (dims, compression), '_images_tiff_page_predictor': (dims, predictor), '_images_tiff_page_dtype': (dims, dtype_array), '_images_tiff_page_levels_downsampling': (dims, levels_downsampling), }, coords=coords, ) metadata = metadata.drop_vars( list(tif_ref_ds.data_vars), errors="ignore", ) metadata = xr.merge( [metadata, tif_ref_ds], compat='override', join='right', ) return metadata def split_imagename_format(imagename_format): """Split the image name format into a base name and an extension. Parameters ---------- imagename_format: str The image name format to split. Returns ------- base: str The base name of the image name format. Preserving any slashes. extension: str The extension of the image name format. Examples -------- >>> split_image_name_format('cam{image_y}_{image_x}.tif') ('cam{image_y}_{image_x}', '.tif') >>> split_imagename_format('foo_{serial_number}/cam{image_y}_{image_x}.png') ('foo_{serial_number}/cam{image_y}_{image_x}', '.png') """ from pathlib import PurePosixPath p = PurePosixPath(imagename_format) base = str(p.with_suffix('')) extension = p.suffix return base, extension