"""mcam_data functions.
The mcam_data is fundamentally a numpy array with metadata that describes
how the data was acquired.
Ramona Optics strives to ensure that data saved using the functions below is:
1. Forward compatible. If you save your data today, we will help ensure you
can read it tomorrow.
2. Round trippable. Data saved with these methods is ensured to be
recoverable bit by bit. If this isn't possible with the chosen format,
a warning will be provided.
3. Metadata completeness: We attempt to store as much information about the
experiment to ensure that experimental parameters can be recovered
when analyzing the data.
We classify storage methods in two general categories:
* Saved: Data using this method is meant to be used by Ramona Optics
applications. We strive to document how the data is stored, to ensure
that you have the necessary information to recovery what you need.
* Export: Data saved using these methods is meant to be used to standard
image processing software that are not compatible with the images above.
We also provide a third methodology called:
* Archival: This data format has default settings chosen such that the data
is compressed in terms of final storage space.
"""
from datetime import UTC, datetime
from pathlib import Path
from typing import Optional
from warnings import warn
import numpy as np
import xarray as xr
# Not exactly part of the public API, but they seem to be open to it?
# https://github.com/pydata/xarray/issues/5081
from xarray.core.indexing import LazilyIndexedArray
from owl.util import get_localzone_name
from ..array._bayer2rgb_array import BayerToRGBAArray, BayerToRGBArray
from ._export import load as _load_exported_data
from ._netcdf import load as _load_netcdf
from ._netcdf import save, save_video
from ._new import new_dataset, new_rgb_dataset, new_rgba_dataset
from ._properties import get_bayer_pattern, get_stack_dimension
__all__ = [
'save',
'save_video',
'load',
'to_rgb',
'to_rgba',
'get_photometric_sensor_corrections',
'get_photometric_corrections',
]
[docs]
def to_rgb(mcam_dataarray, bayer_pattern=None):
"""Convert raw data obtained from the mcam to RGB data.
Parameters
----------
mcam_dataarray:
The raw data acquired from the MCAM. Should be an xarray DataArray object
that was ideally created with `mcam_data.new`.
bayer_pattern: [None, 'rggb', 'bggr', 'grbg', 'gbrg']
If `None` is provided, then the bayer_pattern is inferred from the
`bayer_pattern` coordinate of the `raw_data`. In recent versions,
this coordinate should be automatically populated if the data was
acquired with an MCAM. If `bayer_pattern` is specified, then the
coordinate of the `raw_data` is ignored and this parameter is
passed to `bayer2rgb` as the assumed color ordering.
Returns
-------
rgb_dataarray: xarray DataArray
Deep copy of `mcam_dataarray` with rgb color `mcam_data`.
"""
if 'rgb' in mcam_dataarray.dims or 'rgba' in mcam_dataarray.dims:
raise ValueError('Can only convert raw bayer data to rgb')
if isinstance(mcam_dataarray, xr.DataArray):
return _to_rgb_dataarray(mcam_dataarray, bayer_pattern=bayer_pattern)
# isinstance(mcam_dataarray, xr.Dataset):
if bayer_pattern is None:
bayer_pattern = get_bayer_pattern(mcam_dataarray)
images = _to_rgb_dataarray(
mcam_dataarray.images, bayer_pattern=bayer_pattern)
# Shallow copy to not corrupt the user's dataset
dataset = mcam_dataarray.copy(deep=False)
dataset['images'] = images
if 'bayer_pattern' in dataset:
# Removing the bayer pattern indicates that it is a color image???
del dataset['bayer_pattern']
return dataset
def _to_rgb_dataarray(mcam_dataarray, bayer_pattern):
rgb_data = BayerToRGBArray(mcam_dataarray.variable, bayer_pattern)
rgb_data = LazilyIndexedArray(rgb_data)
rgb_dataarray = new_rgb_dataset(
array=rgb_data,
dims=(*mcam_dataarray.dims, 'rgb'),
coords={**mcam_dataarray.coords, 'rgb': ['r', 'g', 'b']},
)["images"]
rgb_dataarray.attrs = mcam_dataarray.attrs.copy()
return rgb_dataarray
[docs]
def to_rgba(mcam_dataarray, bayer_pattern=None):
"""Converts raw data obtained from the mcam to RGBA data.
Parameters
----------
raw_data:
The raw data acquired from the MCAM. Should be an xarray object
that was ideally created with `mcam_data.new`.
bayer_pattern: [None, 'rggb', 'bggr', 'grbg', 'gbrg']
If `None` is provided, then the bayer_pattern is inferred from the
`bayer_pattern` coordinate of the `raw_data`. In recent versions,
this coordinate should be automatically populated if the data was
acquired with an MCAM. If `bayer_pattern` is specified, then the
coordinate of the `raw_data` is ignored and this parameter is
passed to `bayer2rgba` as the assumed color ordering.
Returns
-------
rgba_dataarray: xarray DataArray
Deep copy of `mcam_dataarray` with rgba color `mcam_data`.
"""
if 'rgb' in mcam_dataarray.dims or 'rgba' in mcam_dataarray.dims:
raise ValueError('Can only convert raw bayer data to rgb.')
if isinstance(mcam_dataarray, xr.DataArray):
return _to_rgba_dataarray(mcam_dataarray, bayer_pattern=bayer_pattern)
# isinstance(mcam_dataarray, xr.Dataset):
if bayer_pattern is None:
bayer_pattern = get_bayer_pattern(mcam_dataarray)
images = _to_rgba_dataarray(
mcam_dataarray.images, bayer_pattern=bayer_pattern)
# Shallow copy to not corrupt the user's dataset
dataset = mcam_dataarray.copy(deep=False)
dataset['images'] = images
return dataset
def _to_rgba_dataarray(mcam_dataarray, bayer_pattern):
rgba_data = BayerToRGBAArray(mcam_dataarray.variable, bayer_pattern)
rgba_data = LazilyIndexedArray(rgba_data)
rgba_dataarray = new_rgba_dataset(
array=rgba_data,
dims=(*mcam_dataarray.dims, 'rgba'),
coords={**mcam_dataarray.coords, 'rgba': ['r', 'g', 'b', 'a']},
)["images"]
rgba_dataarray.attrs = mcam_dataarray.attrs.copy()
return rgba_dataarray
[docs]
def load(
directory: Path | str,
delayed: Optional[bool] = None,
progress: bool = True,
scheduler: str = 'threading',
update: bool = True,
**_kwargs,
) -> xr.Dataset:
"""Load mcam_data and convert it to whatever we want.
Parameters
---------
directory:
The path where the data is stored as bmps.
metadata_filename:
The name of the metadata file. If not provided, a file name
``metadata.json`` or ``metadata.nc`` will be looked for in that order.
delayed:
If True, the computation will return a lazy object. As the user of this
library, you will have to explicitly force the computation of the lazy
object. If None, this function will attempt to automatically determine
if the data should be loaded lazily or eagerly based on the size of the
data. If a particular behavior is desired in your application, the
delayed parameter should be specified. of the lazy object. Metadata is
loaded eagerly.
progress:
If True, then a progress bar is shown during loading operations. This
is only valid if delayed=True.
scheduler:
Parameter passed to ``dask.compute`` to select which schedule is used
to load the data in parallel.
update:
Set to False if you wish to load the raw un-updated data. Setting this
parameter to False keeps the raw metadata as is in the `.nc` file but
means that the loaded dataset is unsupported by the remainder of the
``owl`` analysis functions.
Returns
-------
mcam_data: xarray DataArray or Dataset
Return the mcam_data in an xarray DataArray with all the metadata.
"""
# Resolve the directory now so that delayed operations
# and directory changes don't affect the loading of files
filename_or_directory = Path(directory).resolve()
if filename_or_directory.is_dir():
if filename_or_directory.suffix == '.zarr':
_kwargs.setdefault('metadata_filename', filename_or_directory.name)
directory = filename_or_directory.parent
else:
directory = filename_or_directory
return _load_exported_data(
directory=directory,
delayed=delayed,
progress=progress,
scheduler=scheduler,
update=update,
_stacklevel_increment=1,
**_kwargs
)
return _load_netcdf(
filename=filename_or_directory,
delayed=delayed,
scheduler=scheduler,
progress=progress,
update=update,
_stacklevel_increment=1,
**_kwargs
)
[docs]
def dataset_from_large_image(
image,
N_cameras=(6, 4),
*,
tile_shape=None,
contiguous=True,
trim_edge=False,
):
"""Creates an mcam_data like Dataset from a large image.
Parameters
----------
image: array-like, grayscale [N, M], or rgb [N, M, 3], or rgba [N, M, 4]
The image you wish to chunk up into smaller pieces.
N_cameras: (int, int)
Tuple of int for the number of image tiles to create in in each dimension.
tile_shape: (int, int)
Tuple of int for the shape of each tile. If not provided, the tile shape
is inferred from the image shape and N_cameras.
contiguous: bool, optional
If true, will ensure that the resulting array is contiguous in memory.
trim_edge: bool, optional
If true, will trim away edges to make a smaller image in the case that
the chunks do not line up with the original image size.
Return
------
m: mcam_data
mcam_data like xarray.Dataset structure that holds the image data as
the variable ``images``.
Example
-------
>>> from owl import mcam_data
>>> from imageio import imread
>>> image = imread('large_image_example.tif')
>>> dataset = mcam_data.from_large_image(image, trim_edge=True)
>>> dataset
<xarray.Dataset>
Dimensions: (image_y: 6, image_x: 4, y: 256, x: 128)
Coordinates:
* image_y (image_y) int64 0 1 2 3 4 5
* image_x (image_x) int64 0 1 2 3
* y (y) int64 0 1 2 3 4 5 6 7 ... 249 250 251 252 253 254 255
* x (x) int64 0 1 2 3 4 5 6 7 ... 121 122 123 124 125 126 127
__owl_version__ <U29 '0.18.123.dev6+ga658ddd8.dirty'
__sys_version__ <U80 '3.9.13 | packaged by Ramona Optics | (main, Aug 3...
__owl_sys_info__ <U674 "{'python': '3.9.13 | packaged by Ramona Optics |...
Data variables:
images (image_y, image_x, y, x) float64 0.0 0.0 0.0 ... 0.0 0.0
>>> mcam_data.save(mcam_data, "large_image_chunked", include_timestamp=False)
"""
if image.ndim == 2:
chromaticity = 'gray'
elif image.ndim != 3:
raise ValueError("We only support grayscale (2D) or color images (3D).")
elif image.shape[-1] == 3:
chromaticity = 'rgb'
elif image.shape[-1] == 4:
chromaticity = 'rgba'
else:
raise ValueError(f"Unknown image shape, got shape={image.shape}")
if tile_shape is None:
tile_shape = (image.shape[0] // N_cameras[0], image.shape[1] // N_cameras[1])
new_image_shape = (tile_shape[0] * N_cameras[0], tile_shape[1] * N_cameras[1])
# Any RGB or RGBA channels
new_image_shape += image.shape[2:]
tile_shape += image.shape[2:]
if not trim_edge and new_image_shape != image.shape:
raise ValueError(
"Image shape is not divisible by chunk size. "
f"Input image has shape: {image.shape}, but output image would "
f"have shape {new_image_shape}. To trim edges so that the image "
"fits, set the parameter `trim_edge=True`.")
image = image[:new_image_shape[0], :new_image_shape[1], ...]
new_image = image.reshape(
(N_cameras[0], tile_shape[0], N_cameras[1], tile_shape[1]) + tile_shape[2:]
)
new_image = new_image.transpose(
(0, 2, 1, 3) + tuple(range(4, new_image.ndim))
)
if contiguous:
new_image = np.ascontiguousarray(new_image)
if chromaticity == 'gray':
return new_dataset(
array=new_image,
dims=('image_y', 'image_x', 'y', 'x'))
elif chromaticity == 'rgb':
return new_rgb_dataset(
array=new_image,
dims=('image_y', 'image_x', 'y', 'x', 'rgb'))
elif chromaticity == 'rgba':
return new_rgba_dataset(
array=new_image,
dims=('image_y', 'image_x', 'y', 'x', 'rgba'))
def _get_photometric_pixel_corrections(mcam_dataset, illumination_type):
from owl.calibration.color_correction import create_pixel_corrections
needed_params = ['_pixel_response_coefficient',
'_pixel_response_offset']
if not all(illumination_type + param in mcam_dataset for param in needed_params):
raise ValueError('given mcam_data does not contain the required transform parameters')
else:
warn("`get_photometric_corrections` is deprecated as of version 0.16.22. "
"This function will be removed in version 0.17.0. "
"pixel_corrections can be accessed through\n"
"get_crop_binned_pixel_polynomial().",
stacklevel=2)
image_shape = (len(mcam_dataset['y']), len(mcam_dataset['x']))
pixel_polyco_coefficient = mcam_dataset[illumination_type +
'_pixel_response_coefficient'].data
pixel_polyco_offset = mcam_dataset[illumination_type +
'_pixel_response_offset'].data
(coefficient_pixel_corrections,
offset_pixel_corrections) = create_pixel_corrections(pixel_polyco_coefficient,
pixel_polyco_offset,
image_shape=image_shape)
return coefficient_pixel_corrections, offset_pixel_corrections
[docs]
def get_photometric_sensor_corrections(mcam_dataset, illumination_type):
"""Get the matrix used to create the average value of the individual images.
Parameters
----------
mcam_dataset : xarray.Dataset
An ``xarray.Dataset`` of mcam_data.
illumination_type : string ('reflection', 'transmission')
Which illumination board to use to illuminate the sensors.
Raises
------
ValueError
Unable to get correction if photometric calibration data is not in `mcam_dataset`.
Returns
-------
response_matrix : numpy array
The matrix that describes the sensor's response to the leds. It is a
M x N x 4 x 4 array of float where M and N are the shape of the sensor array.
"""
param = illumination_type + '_photometric_response'
if param not in mcam_dataset:
raise ValueError('given mcam_data does not contain the required transform parameters')
response_matrix = mcam_dataset[param].data
return response_matrix
[docs]
def get_photometric_corrections(mcam_dataset, illumination_type):
"""Get all all color corrections available for a calibrated system.
Parameters
----------
mcam_dataset : xarray.Dataset
An ``xarray.Dataset`` of mcam_data.
illumination_type : string ('reflection', 'transmission')
Which illumination board to use to illuminate the sensors.
Returns
-------
response_matrix : numpy array
The matrix that describes the sensor's response to the leds. It is a
M x N x 4 x 4 array of float where M and N are the shape of the sensor array.
coefficient_corrections : numpy array
Array of pixel correction coefficients of shape M x N x image_shape.
offset_corrections : numpy array
Array of pixel correction offsets of shape M x N x image_shape.
"""
response_matrix = get_photometric_sensor_corrections(mcam_dataset, illumination_type)
(coefficient_pixel_corrections,
offset_pixel_corrections) = _get_photometric_pixel_corrections(mcam_dataset, illumination_type)
return response_matrix, coefficient_pixel_corrections, offset_pixel_corrections
def add_software_timestamp(dataset, *, dims=None, since=None):
"""Add a variable to store the timestamp of each frame to the dataset.
The provided dataset is mutated to add (or overwwrite) a variable called
``'software_timestamp'``
Parameters
----------
dataset: mcam_data
The dataset which will get a the software timestamp.
dims: Tuple[str]
The dimensions of the software_timestamp variable. If not provided
this is inferred from the stack dimensions.
since: datetime
The time from which to store the information in the xarray. For best
results, this time should be close to the times that will be stored
in the array. This object should be timezone aware. Otherwise, the
it will be assume that it refers to localtime.
Returns
-------
dataset:
dataset with software_timestamp added.
"""
if since is None:
# Specifying a timezone aware timestamp will add
# +00:00 at the end of the isoformat
# which will make it explicit that it is UTC time
since = datetime.now(UTC)
tzname = since.tzname()
if tzname is None:
raise RuntimeError(
"An timezone naive since object was passed. "
"Please provide a timezone aware object. "
"See "
"https://docs.python.org/3/library/datetime.html#determining-if-an-object-is-aware-or-naive " # noqa
"for more information",
)
if dims is None:
dims = get_stack_dimension(dataset)
if isinstance(dims, str):
dims = (dims,)
shape = tuple(
dataset.sizes[dim]
for dim in dims
)
variables = {}
local_timezone_name = get_localzone_name()
if local_timezone_name is not None:
# tzinfo and tzname may return things like EDT
# which.... to my surprise, just isn't in th IANA....
# There are a few others like that
# https://github.com/stub42/pytz/issues/71
# instead, get_localzone_name will return fuller qualifiers
# like
# "America/New_York"
# Which is a timezone
variables['local_timezone'] = ((), local_timezone_name)
variables['software_timestamp_timezone'] = ((), tzname)
# Mark - 2023/12/23
# This is a bit of a hack, but I don't want to test the ability to save
# time as nanoseconds with all versions of pandas and xarray
# instead this happened to be the versions installed on my machine
# that seemed to support it.
# if we do not have at least these versions installed, it will fall back
# to the code path we had before Dec 23, 2023 and we will save time as
# it used to.
encoding = {
'units': 'nanoseconds since ' + since.isoformat(),
'calendar': 'proleptic_gregorian',
'dtype': 'int64',
}
variables['software_timestamp'] = (
dims,
np.zeros(shape, dtype='datetime64[ns]'),
{},
encoding
)
return dataset.merge(variables, compat='override')
def _update_timestamp(dataset, frame_rate, time_end, *, dim='frame_number'):
# Update the software timestamp
N_frames = len(dataset.software_timestamp)
delta_time = 1 / frame_rate
dt_ns = (np.arange(N_frames) * delta_time * 1E9).astype(int)
dt_ns = dt_ns - dt_ns[-1]
index = dataset.software_timestamp.dims.index(dim)
ndim = dataset.software_timestamp.ndim
broadcast_slice = (np.newaxis,) * index + (slice(None),) + (np.newaxis,) * (ndim - index - 1)
dt_ns = dt_ns[broadcast_slice]
# Mark: 2021/11/03
# This is quite hacky, but keeps xarray happy
# 1. We convert datetime object to an xarray object
dataset.software_timestamp.data[...] = time_end.replace(tzinfo=None)
# 2. We let xarray do the math between ints and datetime objects
dataset.software_timestamp.data[...] = (
dataset.software_timestamp.data[...] + dt_ns
)