wim/ndbioimage
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- implement sliced views, including min, max, sum and mean operations

Browse Source
This commit is contained in:
Wim Pomp
2025-04-27 20:07:49 +02:00
parent 87e9715f97
commit 5195ccfcb5
15 changed files with 3566 additions and 1068 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
Cargo.toml
View File

@@ -1,6 +1,6 @@
[package]
name = "ndbioimage"
version = "2025.2.3"
version = "2025.4.0"
edition = "2021"
rust-version = "1.78.0"
authors = ["Wim Pomp <w.pomp@nki.nl>"]
@@ -19,15 +19,20 @@ name = "ndbioimage"
crate-type = ["cdylib", "rlib"]
[dependencies]
anyhow = "1.0.95"
anyhow = "1.0.98"
itertools = "0.14.0"
indexmap = { version = "2.9.0", features = ["serde"] }
j4rs = "0.22.0"
ndarray = "0.16.1"
ndarray = { version = "0.16.1", features = ["serde"] }
num = "0.4.3"
numpy = { version = "0.23.0", optional = true }
numpy = { version = "0.24.0", optional = true }
serde = { version = "1.0.219", features = ["rc"] }
serde_json = { version = "1.0.140", optional = true }
serde_with = "3.12.0"
thread_local = "1.1.8"
[dependencies.pyo3]
version = "0.23.4"
version = "0.24.2"
features = ["extension-module", "abi3-py310", "generate-import-lib", "anyhow"]
optional = true
@@ -35,10 +40,10 @@ optional = true
rayon = "1.10.0"
[build-dependencies]
anyhow = "1.0.95"
anyhow = "1.0.98"
j4rs = "0.22.0"
retry = "2.0.0"
retry = "2.1.0"
[features]
python = ["dep:pyo3", "dep:numpy"]
python = ["dep:pyo3", "dep:numpy", "dep:serde_json"]
gpl-formats = []

32
README.md
View File

@@ -74,7 +74,7 @@ use ndarray::Array2;
use ndbioimage::Reader;
let path = "/path/to/file";
let reader = Reader::new(&path, 0).unwrap();
let reader = Reader::new(&path, 0)?;
println!("size: {}, {}", reader.size_y, reader.size_y);
let frame = reader.get_frame(0, 0, 0).unwrap();
if let Ok(arr) = <Frame as TryInto<Array2<i8>>>::try_into(frame) {
@@ -86,30 +86,22 @@ let xml = reader.get_ome_xml().unwrap();
println!("{}", xml);
```
```
use ndarray::Array2;
use ndbioimage::Reader;
let path = "/path/to/file";
let reader = Reader::new(&path, 0)?;
let view = reader.view();
let view = view.max_proj(3)?;
let array = view.as_array::<u16>()?
```
### Command line
```ndbioimage --help```: show help
```ndbioimage image```: show metadata about image
```ndbioimage image -w {name}.tif -r```: copy image into image.tif (replacing {name} with image), while registering channels
```ndbioimage image -w image.mp4 -C cyan lime red``` copy image into image.mp4 (z will be max projected), make channel colors cyan lime and red
## Adding more formats
Readers for image formats subclass AbstractReader. When an image reader is imported, Imread will
automatically recognize it and use it to open the appropriate file format. Image readers
are required to implement the following methods:
- staticmethod _can_open(path): return True if path can be opened by this reader
- \_\_frame__(self, c, z, t): return the frame at channel=c, z-slice=z, time=t from the file
Optional methods:
- get_ome: reads metadata from file and adds them to an OME object imported
from the ome-types library
- open(self): maybe open some file handle
- close(self): close any file handles
Optional fields:
- priority (int): Imread will try readers with a lower number first, default: 99
- do_not_pickle (strings): any attributes that should not be included when the object is pickled,
for example: any file handles
# TODO
- more image formats

2
build.rs
View File

@@ -8,7 +8,7 @@ use retry::{delay, delay::Exponential, retry};
use j4rs::Jvm;
fn main() -> anyhow::Result<()> {
println!("cargo:rerun-if-changed=build.rs");
println!("cargo::rerun-if-changed=build.rs");
#[cfg(not(feature = "python"))]
if std::env::var("DOCS_RS").is_err() {

1414
py/ndbioimage/__init__.py
View File

File diff suppressed because it is too large Load Diff

296
py/ndbioimage/transforms.py
View File

@@ -21,7 +21,7 @@ except ImportError:
DataFrame, Series, concat = None, None, None
if hasattr(yaml, 'full_load'):
if hasattr(yaml, "full_load"):
yamlload = yaml.full_load
else:
yamlload = yaml.load
@@ -34,7 +34,7 @@ class Transforms(dict):
@classmethod
def from_file(cls, file, C=True, T=True):
with open(Path(file).with_suffix('.yml')) as f:
with open(Path(file).with_suffix(".yml")) as f:
return cls.from_dict(yamlload(f), C, T)
@classmethod
@@ -42,7 +42,9 @@ class Transforms(dict):
new = cls()
for key, value in d.items():
if isinstance(key, str) and C:
new[key.replace(r'\:', ':').replace('\\\\', '\\')] = Transform.from_dict(value)
new[key.replace(r"\:", ":").replace("\\\\", "\\")] = (
Transform.from_dict(value)
)
elif T:
new[key] = Transform.from_dict(value)
return new
@@ -69,11 +71,19 @@ class Transforms(dict):
return new
def asdict(self):
return {key.replace('\\', '\\\\').replace(':', r'\:') if isinstance(key, str) else key: value.asdict()
for key, value in self.items()}
return {
key.replace("\\", "\\\\").replace(":", r"\:")
if isinstance(key, str)
else key: value.asdict()
for key, value in self.items()
}
def __getitem__(self, item):
return np.prod([self[i] for i in item[::-1]]) if isinstance(item, tuple) else super().__getitem__(item)
return (
np.prod([self[i] for i in item[::-1]])
if isinstance(item, tuple)
else super().__getitem__(item)
)
def __missing__(self, key):
return self.default
@@ -88,7 +98,7 @@ class Transforms(dict):
return hash(frozenset((*self.__dict__.items(), *self.items())))
def save(self, file):
with open(Path(file).with_suffix('.yml'), 'w') as f:
with open(Path(file).with_suffix(".yml"), "w") as f:
yaml.safe_dump(self.asdict(), f, default_flow_style=None)
def copy(self):
@@ -109,8 +119,10 @@ class Transforms(dict):
transform_channels = {key for key in self.keys() if isinstance(key, str)}
if set(channel_names) - transform_channels:
mapping = key_map(channel_names, transform_channels)
warnings.warn(f'The image file and the transform do not have the same channels,'
f' creating a mapping: {mapping}')
warnings.warn(
f"The image file and the transform do not have the same channels,"
f" creating a mapping: {mapping}"
)
for key_im, key_t in mapping.items():
self[key_im] = self[key_t]
@@ -124,37 +136,54 @@ class Transforms(dict):
def coords_pandas(self, array, channel_names, columns=None):
if isinstance(array, DataFrame):
return concat([self.coords_pandas(row, channel_names, columns) for _, row in array.iterrows()], axis=1).T
return concat(
[
self.coords_pandas(row, channel_names, columns)
for _, row in array.iterrows()
],
axis=1,
).T
elif isinstance(array, Series):
key = []
if 'C' in array:
key.append(channel_names[int(array['C'])])
if 'T' in array:
key.append(int(array['T']))
if "C" in array:
key.append(channel_names[int(array["C"])])
if "T" in array:
key.append(int(array["T"]))
return self[tuple(key)].coords(array, columns)
else:
raise TypeError('Not a pandas DataFrame or Series.')
raise TypeError("Not a pandas DataFrame or Series.")
def with_beads(self, cyllens, bead_files):
assert len(bead_files) > 0, 'At least one file is needed to calculate the registration.'
transforms = [self.calculate_channel_transforms(file, cyllens) for file in bead_files]
assert len(bead_files) > 0, (
"At least one file is needed to calculate the registration."
)
transforms = [
self.calculate_channel_transforms(file, cyllens) for file in bead_files
]
for key in {key for transform in transforms for key in transform.keys()}:
new_transforms = [transform[key] for transform in transforms if key in transform]
new_transforms = [
transform[key] for transform in transforms if key in transform
]
if len(new_transforms) == 1:
self[key] = new_transforms[0]
else:
self[key] = Transform()
self[key].parameters = np.mean([t.parameters for t in new_transforms], 0)
self[key].dparameters = (np.std([t.parameters for t in new_transforms], 0) /
np.sqrt(len(new_transforms))).tolist()
self[key].parameters = np.mean(
[t.parameters for t in new_transforms], 0
)
self[key].dparameters = (
np.std([t.parameters for t in new_transforms], 0)
/ np.sqrt(len(new_transforms))
).tolist()
return self
@staticmethod
def get_bead_files(path):
from . import Imread
files = []
for file in path.iterdir():
if file.name.lower().startswith('beads'):
if file.name.lower().startswith("beads"):
try:
with Imread(file):
files.append(file)
@@ -162,32 +191,36 @@ class Transforms(dict):
pass
files = sorted(files)
if not files:
raise Exception('No bead file found!')
raise Exception("No bead file found!")
checked_files = []
for file in files:
try:
if file.is_dir():
file /= 'Pos0'
file /= "Pos0"
with Imread(file): # check for errors opening the file
checked_files.append(file)
except (Exception,):
continue
if not checked_files:
raise Exception('No bead file found!')
raise Exception("No bead file found!")
return checked_files
@staticmethod
def calculate_channel_transforms(bead_file, cyllens):
""" When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162
in the horizontal direction """
"""When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162
in the horizontal direction"""
from . import Imread
with Imread(bead_file, axes='zcyx') as im: # noqa
max_ims = im.max('z')
with Imread(bead_file, axes="zcyx") as im: # noqa
max_ims = im.max("z")
goodch = [c for c, max_im in enumerate(max_ims) if not im.is_noise(max_im)]
if not goodch:
goodch = list(range(len(max_ims)))
untransformed = [c for c in range(im.shape['c']) if cyllens[im.detector[c]].lower() == 'none']
untransformed = [
c
for c in range(im.shape["c"])
if cyllens[im.detector[c]].lower() == "none"
]
good_and_untrans = sorted(set(goodch) & set(untransformed))
if good_and_untrans:
@@ -200,54 +233,81 @@ class Transforms(dict):
matrix[0, 0] = 0.86
transform.matrix = matrix
transforms = Transforms()
for c in tqdm(goodch, desc='Calculating channel transforms'): # noqa
for c in tqdm(goodch, desc="Calculating channel transforms"): # noqa
if c == masterch:
transforms[im.channel_names[c]] = transform
else:
transforms[im.channel_names[c]] = Transform.register(max_ims[masterch], max_ims[c]) * transform
transforms[im.channel_names[c]] = (
Transform.register(max_ims[masterch], max_ims[c]) * transform
)
return transforms
@staticmethod
def save_channel_transform_tiff(bead_files, tiffile):
from . import Imread
n_channels = 0
for file in bead_files:
with Imread(file) as im:
n_channels = max(n_channels, im.shape['c'])
n_channels = max(n_channels, im.shape["c"])
with IJTiffFile(tiffile) as tif:
for t, file in enumerate(bead_files):
with Imread(file) as im:
with Imread(file).with_transform() as jm:
for c in range(im.shape['c']):
tif.save(np.hstack((im(c=c, t=0).max('z'), jm(c=c, t=0).max('z'))), c, 0, t)
for c in range(im.shape["c"]):
tif.save(
np.hstack(
(im(c=c, t=0).max("z"), jm(c=c, t=0).max("z"))
),
c,
0,
t,
)
def with_drift(self, im):
""" Calculate shifts relative to the first frame
divide the sequence into groups,
compare each frame to the frame in the middle of the group and compare these middle frames to each other
"""Calculate shifts relative to the first frame
divide the sequence into groups,
compare each frame to the frame in the middle of the group and compare these middle frames to each other
"""
im = im.transpose('tzycx')
t_groups = [list(chunk) for chunk in Chunks(range(im.shape['t']), size=round(np.sqrt(im.shape['t'])))]
im = im.transpose("tzycx")
t_groups = [
list(chunk)
for chunk in Chunks(
range(im.shape["t"]), size=round(np.sqrt(im.shape["t"]))
)
]
t_keys = [int(np.round(np.mean(t_group))) for t_group in t_groups]
t_pairs = [(int(np.round(np.mean(t_group))), frame) for t_group in t_groups for frame in t_group]
t_pairs = [
(int(np.round(np.mean(t_group))), frame)
for t_group in t_groups
for frame in t_group
]
t_pairs.extend(zip(t_keys, t_keys[1:]))
fmaxz_keys = {t_key: filters.gaussian(im[t_key].max('z'), 5) for t_key in t_keys}
fmaxz_keys = {
t_key: filters.gaussian(im[t_key].max("z"), 5) for t_key in t_keys
}
def fun(t_key_t, im, fmaxz_keys):
t_key, t = t_key_t
if t_key == t:
return 0, 0
else:
fmaxz = filters.gaussian(im[t].max('z'), 5)
return Transform.register(fmaxz_keys[t_key], fmaxz, 'translation').parameters[4:]
fmaxz = filters.gaussian(im[t].max("z"), 5)
return Transform.register(
fmaxz_keys[t_key], fmaxz, "translation"
).parameters[4:]
shifts = np.array(pmap(fun, t_pairs, (im, fmaxz_keys), desc='Calculating image shifts.'))
shifts = np.array(
pmap(fun, t_pairs, (im, fmaxz_keys), desc="Calculating image shifts.")
)
shift_keys_cum = np.zeros(2)
for shift_keys, t_group in zip(np.vstack((-shifts[0], shifts[im.shape['t']:])), t_groups):
for shift_keys, t_group in zip(
np.vstack((-shifts[0], shifts[im.shape["t"] :])), t_groups
):
shift_keys_cum += shift_keys
shifts[t_group] += shift_keys_cum
for i, shift in enumerate(shifts[:im.shape['t']]):
for i, shift in enumerate(shifts[: im.shape["t"]]):
self[i] = Transform.from_shift(shift)
return self
@@ -257,9 +317,11 @@ class Transform:
if sitk is None:
self.transform = None
else:
self.transform = sitk.ReadTransform(str(Path(__file__).parent / 'transform.txt'))
self.dparameters = [0., 0., 0., 0., 0., 0.]
self.shape = [512., 512.]
self.transform = sitk.ReadTransform(
str(Path(__file__).parent / "transform.txt")
)
self.dparameters = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
self.shape = [512.0, 512.0]
self.origin = [255.5, 255.5]
self._last, self._inverse = None, None
@@ -274,12 +336,14 @@ class Transform:
@classmethod
def register(cls, fix, mov, kind=None):
""" kind: 'affine', 'translation', 'rigid' """
"""kind: 'affine', 'translation', 'rigid'"""
if sitk is None:
raise ImportError('SimpleElastix is not installed: '
'https://simpleelastix.readthedocs.io/GettingStarted.html')
raise ImportError(
"SimpleElastix is not installed: "
"https://simpleelastix.readthedocs.io/GettingStarted.html"
)
new = cls()
kind = kind or 'affine'
kind = kind or "affine"
new.shape = fix.shape
fix, mov = new.cast_image(fix), new.cast_image(mov)
# TODO: implement RigidTransform
@@ -290,16 +354,18 @@ class Transform:
tfilter.SetParameterMap(sitk.GetDefaultParameterMap(kind))
tfilter.Execute()
transform = tfilter.GetTransformParameterMap()[0]
if kind == 'affine':
new.parameters = [float(t) for t in transform['TransformParameters']]
new.shape = [float(t) for t in transform['Size']]
new.origin = [float(t) for t in transform['CenterOfRotationPoint']]
elif kind == 'translation':
new.parameters = [1.0, 0.0, 0.0, 1.0] + [float(t) for t in transform['TransformParameters']]
new.shape = [float(t) for t in transform['Size']]
if kind == "affine":
new.parameters = [float(t) for t in transform["TransformParameters"]]
new.shape = [float(t) for t in transform["Size"]]
new.origin = [float(t) for t in transform["CenterOfRotationPoint"]]
elif kind == "translation":
new.parameters = [1.0, 0.0, 0.0, 1.0] + [
float(t) for t in transform["TransformParameters"]
]
new.shape = [float(t) for t in transform["Size"]]
new.origin = [(t - 1) / 2 for t in new.shape]
else:
raise NotImplementedError(f'{kind} tranforms not implemented (yet)')
raise NotImplementedError(f"{kind} tranforms not implemented (yet)")
new.dparameters = 6 * [np.nan]
return new
@@ -315,18 +381,35 @@ class Transform:
@classmethod
def from_file(cls, file):
with open(Path(file).with_suffix('.yml')) as f:
with open(Path(file).with_suffix(".yml")) as f:
return cls.from_dict(yamlload(f))
@classmethod
def from_dict(cls, d):
new = cls()
new.origin = None if d['CenterOfRotationPoint'] is None else [float(i) for i in d['CenterOfRotationPoint']]
new.parameters = ((1., 0., 0., 1., 0., 0.) if d['TransformParameters'] is None else
[float(i) for i in d['TransformParameters']])
new.dparameters = ([(0., 0., 0., 0., 0., 0.) if i is None else float(i) for i in d['dTransformParameters']]
if 'dTransformParameters' in d else 6 * [np.nan] and d['dTransformParameters'] is not None)
new.shape = None if d['Size'] is None else [None if i is None else float(i) for i in d['Size']]
new.origin = (
None
if d["CenterOfRotationPoint"] is None
else [float(i) for i in d["CenterOfRotationPoint"]]
)
new.parameters = (
(1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
if d["TransformParameters"] is None
else [float(i) for i in d["TransformParameters"]]
)
new.dparameters = (
[
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0) if i is None else float(i)
for i in d["dTransformParameters"]
]
if "dTransformParameters" in d
else 6 * [np.nan] and d["dTransformParameters"] is not None
)
new.shape = (
None
if d["Size"] is None
else [None if i is None else float(i) for i in d["Size"]]
)
return new
def __mul__(self, other): # TODO: take care of dmatrix
@@ -359,9 +442,13 @@ class Transform:
@property
def matrix(self):
return np.array(((*self.parameters[:2], self.parameters[4]),
(*self.parameters[2:4], self.parameters[5]),
(0, 0, 1)))
return np.array(
(
(*self.parameters[:2], self.parameters[4]),
(*self.parameters[2:4], self.parameters[5]),
(0, 0, 1),
)
)
@matrix.setter
def matrix(self, value):
@@ -370,9 +457,13 @@ class Transform:
@property
def dmatrix(self):
return np.array(((*self.dparameters[:2], self.dparameters[4]),
(*self.dparameters[2:4], self.dparameters[5]),
(0, 0, 0)))
return np.array(
(
(*self.dparameters[:2], self.dparameters[4]),
(*self.dparameters[2:4], self.dparameters[5]),
(0, 0, 0),
)
)
@dmatrix.setter
def dmatrix(self, value):
@@ -384,7 +475,7 @@ class Transform:
if self.transform is not None:
return list(self.transform.GetParameters())
else:
return [1., 0., 0., 1., 0., 0.]
return [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]
@parameters.setter
def parameters(self, value):
@@ -420,43 +511,62 @@ class Transform:
self.shape = shape[:2]
def asdict(self):
return {'CenterOfRotationPoint': self.origin, 'Size': self.shape, 'TransformParameters': self.parameters,
'dTransformParameters': np.nan_to_num(self.dparameters, nan=1e99).tolist()}
return {
"CenterOfRotationPoint": self.origin,
"Size": self.shape,
"TransformParameters": self.parameters,
"dTransformParameters": np.nan_to_num(self.dparameters, nan=1e99).tolist(),
}
def frame(self, im, default=0):
if self.is_unity():
return im
else:
if sitk is None:
raise ImportError('SimpleElastix is not installed: '
'https://simpleelastix.readthedocs.io/GettingStarted.html')
raise ImportError(
"SimpleElastix is not installed: "
"https://simpleelastix.readthedocs.io/GettingStarted.html"
)
dtype = im.dtype
im = im.astype('float')
intp = sitk.sitkBSpline if np.issubdtype(dtype, np.floating) else sitk.sitkNearestNeighbor
return self.cast_array(sitk.Resample(self.cast_image(im), self.transform, intp, default)).astype(dtype)
im = im.astype("float")
intp = (
sitk.sitkBSpline
if np.issubdtype(dtype, np.floating)
else sitk.sitkNearestNeighbor
)
return self.cast_array(
sitk.Resample(self.cast_image(im), self.transform, intp, default)
).astype(dtype)
def coords(self, array, columns=None):
""" Transform coordinates in 2 column numpy array,
or in pandas DataFrame or Series objects in columns ['x', 'y']
"""Transform coordinates in 2 column numpy array,
or in pandas DataFrame or Series objects in columns ['x', 'y']
"""
if self.is_unity():
return array.copy()
elif DataFrame is not None and isinstance(array, (DataFrame, Series)):
columns = columns or ['x', 'y']
columns = columns or ["x", "y"]
array = array.copy()
if isinstance(array, DataFrame):
array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))
elif isinstance(array, Series):
array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))[0]
array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))[
0
]
return array
else: # somehow we need to use the inverse here to get the same effect as when using self.frame
return np.array([self.inverse.transform.TransformPoint(i.tolist()) for i in np.asarray(array)])
return np.array(
[
self.inverse.transform.TransformPoint(i.tolist())
for i in np.asarray(array)
]
)
def save(self, file):
""" save the parameters of the transform calculated
with affine_registration to a yaml file
"""save the parameters of the transform calculated
with affine_registration to a yaml file
"""
if not file[-3:] == 'yml':
file += '.yml'
with open(file, 'w') as f:
if not file[-3:] == "yml":
file += ".yml"
with open(file, "w") as f:
yaml.safe_dump(self.asdict(), f, default_flow_style=None)

218
src/axes.rs Normal file
View File

@@ -0,0 +1,218 @@
use crate::stats::MinMax;
use anyhow::{anyhow, Error, Result};
use ndarray::{Array, Dimension, Ix2, SliceInfo, SliceInfoElem};
use serde::{Deserialize, Deserializer, Serialize, Serializer};
use serde_with::{DeserializeAs, SerializeAs};
use std::hash::{Hash, Hasher};
use std::str::FromStr;
/// a trait to find axis indices from any object
pub trait Ax {
/// C: 0, Z: 1, T: 2, Y: 3, X: 4
fn n(&self) -> usize;
/// the indices of axes in self.axes, which always has all of CZTYX
fn pos(&self, axes: &[Axis], slice: &[SliceInfoElem]) -> Result<usize>;
/// the indices of axes in self.axes, which always has all of CZTYX, but skip axes with an operation
fn pos_op(&self, axes: &[Axis], slice: &[SliceInfoElem], op_axes: &[Axis]) -> Result<usize>;
}
/// Enum for CZTYX axes or a new axis
#[derive(Clone, Copy, Debug, Eq, Ord, PartialOrd, Serialize, Deserialize)]
pub enum Axis {
C,
Z,
T,
Y,
X,
New,
}
impl Hash for Axis {
fn hash<H: Hasher>(&self, state: &mut H) {
(*self as usize).hash(state);
}
}
impl FromStr for Axis {
type Err = Error;
fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
match s.to_uppercase().as_str() {
"C" => Ok(Axis::C),
"Z" => Ok(Axis::Z),
"T" => Ok(Axis::T),
"Y" => Ok(Axis::Y),
"X" => Ok(Axis::X),
"NEW" => Ok(Axis::New),
_ => Err(anyhow!("invalid axis: {}", s)),
}
}
}
impl Ax for Axis {
fn n(&self) -> usize {
*self as usize
}
fn pos(&self, axes: &[Axis], _slice: &[SliceInfoElem]) -> Result<usize> {
if let Some(pos) = axes.iter().position(|a| a == self) {
Ok(pos)
} else {
Err(Error::msg("Axis not found in axes"))
}
}
fn pos_op(&self, axes: &[Axis], _slice: &[SliceInfoElem], _op_axes: &[Axis]) -> Result<usize> {
self.pos(axes, _slice)
}
}
impl Ax for usize {
fn n(&self) -> usize {
*self
}
fn pos(&self, _axes: &[Axis], slice: &[SliceInfoElem]) -> Result<usize> {
let idx: Vec<_> = slice
.iter()
.enumerate()
.filter_map(|(i, s)| if s.is_index() { None } else { Some(i) })
.collect();
Ok(idx[*self])
}
fn pos_op(&self, axes: &[Axis], slice: &[SliceInfoElem], op_axes: &[Axis]) -> Result<usize> {
let idx: Vec<_> = axes
.iter()
.zip(slice.iter())
.enumerate()
.filter_map(|(i, (ax, s))| {
if s.is_index() | op_axes.contains(ax) {
None
} else {
Some(i)
}
})
.collect();
assert!(*self < idx.len(), "self: {}, idx: {:?}", self, idx);
Ok(idx[*self])
}
}
#[derive(Clone, Debug, Serialize, Deserialize)]
pub(crate) enum Operation {
Max,
Min,
Sum,
Mean,
}
impl Operation {
pub(crate) fn operate<T, D>(
&self,
array: Array<T, D>,
axis: usize,
) -> Result<<Array<T, D> as MinMax>::Output>
where
D: Dimension,
Array<T, D>: MinMax,
{
match self {
Operation::Max => array.max(axis),
Operation::Min => array.min(axis),
Operation::Sum => array.sum(axis),
Operation::Mean => array.mean(axis),
}
}
}
impl PartialEq for Axis {
fn eq(&self, other: &Self) -> bool {
(*self as u8) == (*other as u8)
}
}
pub(crate) fn slice_info<D: Dimension>(
info: &[SliceInfoElem],
) -> Result<SliceInfo<&[SliceInfoElem], Ix2, D>> {
Ok(info.try_into()?)
}
#[derive(Serialize, Deserialize)]
#[serde(remote = "SliceInfoElem")]
pub(crate) enum SliceInfoElemDef {
Slice {
start: isize,
end: Option<isize>,
step: isize,
},
Index(isize),
NewAxis,
}
impl SerializeAs<SliceInfoElem> for SliceInfoElemDef {
fn serialize_as<S>(
source: &SliceInfoElem,
serializer: S,
) -> std::result::Result<S::Ok, S::Error>
where
S: Serializer,
{
SliceInfoElemDef::serialize(source, serializer)
}
}
impl<'de> DeserializeAs<'de, SliceInfoElem> for SliceInfoElemDef {
fn deserialize_as<D>(deserializer: D) -> std::result::Result<SliceInfoElem, D::Error>
where
D: Deserializer<'de>,
{
SliceInfoElemDef::deserialize(deserializer)
}
}
#[derive(Clone, Debug)]
pub(crate) struct Slice {
start: isize,
end: isize,
step: isize,
}
impl Slice {
pub(crate) fn new(start: isize, end: isize, step: isize) -> Self {
Self { start, end, step }
}
pub(crate) fn empty() -> Self {
Self {
start: 0,
end: 0,
step: 1,
}
}
}
impl Iterator for Slice {
type Item = isize;
fn next(&mut self) -> Option<Self::Item> {
if self.end - self.start >= self.step {
let r = self.start;
self.start += self.step;
Some(r)
} else {
None
}
}
}
impl IntoIterator for &Slice {
type Item = isize;
type IntoIter = Slice;
fn into_iter(self) -> Self::IntoIter {
self.clone()
}
}

18
src/bioformats.rs
View File

@@ -96,6 +96,18 @@ macro_rules! method {
};
}
fn transmute_vec<T, U>(vec: Vec<T>) -> Vec<U> {
unsafe {
// Ensure the original vector is not dropped.
let mut v_clone = std::mem::ManuallyDrop::new(vec);
Vec::from_raw_parts(
v_clone.as_mut_ptr() as *mut U,
v_clone.len(),
v_clone.capacity(),
)
}
}
/// Wrapper around bioformats java class loci.common.DebugTools
pub struct DebugTools;
@@ -125,8 +137,7 @@ impl ChannelSeparator {
}
pub(crate) fn open_bytes(&self, index: i32) -> Result<Vec<u8>> {
let bi8 = self.open_bi8(index)?;
Ok(unsafe { std::mem::transmute::<Vec<i8>, Vec<u8>>(bi8) })
Ok(transmute_vec(self.open_bi8(index)?))
}
method!(open_bi8, "openBytes", [index: i32|p] => Vec<i8>|c);
@@ -149,8 +160,7 @@ impl ImageReader {
}
pub(crate) fn open_bytes(&self, index: i32) -> Result<Vec<u8>> {
let bi8 = self.open_bi8(index)?;
Ok(unsafe { std::mem::transmute::<Vec<i8>, Vec<u8>>(bi8) })
Ok(transmute_vec(self.open_bi8(index)?))
}
pub(crate) fn ome_xml(&self) -> Result<String> {

488
src/lib.rs
View File

@@ -1,333 +1,23 @@
mod bioformats;
mod axes;
#[cfg(feature = "python")]
mod py;
use anyhow::{anyhow, Result};
use bioformats::{DebugTools, ImageReader, MetadataTools};
use ndarray::Array2;
use num::{FromPrimitive, Zero};
use std::any::type_name;
use std::fmt::Debug;
use std::ops::Deref;
use std::path::{Path, PathBuf};
use thread_local::ThreadLocal;
/// Pixel types (u)int(8/16/32) or float(32/64)
#[derive(Clone, Debug)]
pub enum PixelType {
INT8 = 0,
UINT8 = 1,
INT16 = 2,
UINT16 = 3,
INT32 = 4,
UINT32 = 5,
FLOAT = 6,
DOUBLE = 7,
}
impl TryFrom<i32> for PixelType {
type Error = anyhow::Error;
fn try_from(value: i32) -> Result<Self, Self::Error> {
match value {
0 => Ok(PixelType::INT8),
1 => Ok(PixelType::UINT8),
2 => Ok(PixelType::INT16),
3 => Ok(PixelType::UINT16),
4 => Ok(PixelType::INT32),
5 => Ok(PixelType::UINT32),
6 => Ok(PixelType::FLOAT),
7 => Ok(PixelType::DOUBLE),
_ => Err(anyhow::anyhow!("Unknown pixel type {}", value)),
}
}
}
/// Struct containing frame data in one of eight pixel types. Cast to `Array2<T>` using try_into.
#[derive(Clone, Debug)]
pub enum Frame {
INT8(Array2<i8>),
UINT8(Array2<u8>),
INT16(Array2<i16>),
UINT16(Array2<u16>),
INT32(Array2<i32>),
UINT32(Array2<u32>),
FLOAT(Array2<f32>),
DOUBLE(Array2<f64>),
}
macro_rules! impl_frame_cast {
($t:tt, $s:ident) => {
impl From<Array2<$t>> for Frame {
fn from(value: Array2<$t>) -> Self {
Frame::$s(value)
}
}
};
}
impl_frame_cast!(i8, INT8);
impl_frame_cast!(u8, UINT8);
impl_frame_cast!(i16, INT16);
impl_frame_cast!(u16, UINT16);
impl_frame_cast!(i32, INT32);
impl_frame_cast!(u32, UINT32);
impl_frame_cast!(f32, FLOAT);
impl_frame_cast!(f64, DOUBLE);
impl<T> TryInto<Array2<T>> for Frame
where
T: FromPrimitive + Zero + 'static,
{
type Error = anyhow::Error;
fn try_into(self) -> std::result::Result<Array2<T>, Self::Error> {
let mut err = Ok(());
let arr = match self {
Frame::INT8(v) => v.mapv_into_any(|x| {
T::from_i8(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::UINT8(v) => v.mapv_into_any(|x| {
T::from_u8(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::INT16(v) => v.mapv_into_any(|x| {
T::from_i16(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::UINT16(v) => v.mapv_into_any(|x| {
T::from_u16(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::INT32(v) => v.mapv_into_any(|x| {
T::from_i32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::UINT32(v) => v.mapv_into_any(|x| {
T::from_u32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::FLOAT(v) => v.mapv_into_any(|x| {
T::from_f32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::DOUBLE(v) => v.mapv_into_any(|x| {
T::from_f64(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
};
match err {
Err(err) => Err(err),
Ok(()) => Ok(arr),
}
}
}
/// Reader interface to file. Use get_frame to get data.
pub struct Reader {
image_reader: ThreadLocal<ImageReader>,
/// path to file
pub path: PathBuf,
/// which (if more than 1) of the series in the file to open
pub series: i32,
/// size x (horizontal)
pub size_x: usize,
/// size y (vertical)
pub size_y: usize,
/// size c (# channels)
pub size_c: usize,
/// size z (# slices)
pub size_z: usize,
/// size t (# time/frames)
pub size_t: usize,
/// pixel type ((u)int(8/16/32) or float(32/64))
pub pixel_type: PixelType,
little_endian: bool,
}
impl Deref for Reader {
type Target = ImageReader;
fn deref(&self) -> &Self::Target {
self.image_reader.get_or(|| {
let reader = ImageReader::new().unwrap();
let meta_data_tools = MetadataTools::new().unwrap();
let ome_meta = meta_data_tools.create_ome_xml_metadata().unwrap();
reader.set_metadata_store(ome_meta).unwrap();
reader.set_id(self.path.to_str().unwrap()).unwrap();
reader.set_series(self.series).unwrap();
reader
})
}
}
impl Clone for Reader {
fn clone(&self) -> Self {
Reader::new(&self.path, self.series).unwrap()
}
}
impl Debug for Reader {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Reader")
.field("path", &self.path)
.field("series", &self.series)
.field("size_x", &self.size_x)
.field("size_y", &self.size_y)
.field("size_c", &self.size_c)
.field("size_z", &self.size_z)
.field("size_t", &self.size_t)
.field("pixel_type", &self.pixel_type)
.field("little_endian", &self.little_endian)
.finish()
}
}
impl Reader {
/// Create new reader for image file at path.
pub fn new(path: &Path, series: i32) -> Result<Self> {
DebugTools::set_root_level("ERROR")?;
let mut reader = Reader {
image_reader: ThreadLocal::new(),
path: PathBuf::from(path),
series,
size_x: 0,
size_y: 0,
size_c: 0,
size_z: 0,
size_t: 0,
pixel_type: PixelType::INT8,
little_endian: false,
};
reader.size_x = reader.get_size_x()? as usize;
reader.size_y = reader.get_size_y()? as usize;
reader.size_c = reader.get_size_c()? as usize;
reader.size_z = reader.get_size_z()? as usize;
reader.size_t = reader.get_size_t()? as usize;
reader.pixel_type = PixelType::try_from(reader.get_pixel_type()?)?;
reader.little_endian = reader.is_little_endian()?;
Ok(reader)
}
/// Get ome metadata as xml string
pub fn get_ome_xml(&self) -> Result<String> {
self.ome_xml()
}
fn deinterleave(&self, bytes: Vec<u8>, channel: usize) -> Result<Vec<u8>> {
let chunk_size = match self.pixel_type {
PixelType::INT8 => 1,
PixelType::UINT8 => 1,
PixelType::INT16 => 2,
PixelType::UINT16 => 2,
PixelType::INT32 => 4,
PixelType::UINT32 => 4,
PixelType::FLOAT => 4,
PixelType::DOUBLE => 8,
};
Ok(bytes
.chunks(chunk_size)
.skip(channel)
.step_by(self.size_c)
.flat_map(|a| a.to_vec())
.collect())
}
/// Retrieve fame at channel c, slize z and time t.
pub fn get_frame(&self, c: usize, z: usize, t: usize) -> Result<Frame> {
let bytes = if self.is_rgb()? & self.is_interleaved()? {
let index = self.get_index(z as i32, 0, t as i32)?;
self.deinterleave(self.open_bytes(index)?, c)?
} else if self.get_rgb_channel_count()? > 1 {
let channel_separator = bioformats::ChannelSeparator::new(self)?;
let index = channel_separator.get_index(z as i32, c as i32, t as i32)?;
channel_separator.open_bytes(index)?
} else if self.is_indexed()? {
let index = self.get_index(z as i32, 0, t as i32)?;
self.open_bytes(index)?
// TODO: apply LUT
// let _bytes_lut = match self.pixel_type {
// PixelType::INT8 | PixelType::UINT8 => {
// let _lut = self.image_reader.get_8bit_lookup_table()?;
// }
// PixelType::INT16 | PixelType::UINT16 => {
// let _lut = self.image_reader.get_16bit_lookup_table()?;
// }
// _ => {}
// };
} else {
let index = self.get_index(z as i32, c as i32, t as i32)?;
self.open_bytes(index)?
};
self.bytes_to_frame(bytes)
}
fn bytes_to_frame(&self, bytes: Vec<u8>) -> Result<Frame> {
macro_rules! get_frame {
($t:tt, <$n:expr) => {
Ok(Frame::from(Array2::from_shape_vec(
(self.size_y, self.size_x),
bytes
.chunks($n)
.map(|x| $t::from_le_bytes(x.try_into().unwrap()))
.collect(),
)?))
};
($t:tt, >$n:expr) => {
Ok(Frame::from(Array2::from_shape_vec(
(self.size_y, self.size_x),
bytes
.chunks($n)
.map(|x| $t::from_be_bytes(x.try_into().unwrap()))
.collect(),
)?))
};
}
match (&self.pixel_type, self.little_endian) {
(PixelType::INT8, true) => get_frame!(i8, <1),
(PixelType::UINT8, true) => get_frame!(u8, <1),
(PixelType::INT16, true) => get_frame!(i16, <2),
(PixelType::UINT16, true) => get_frame!(u16, <2),
(PixelType::INT32, true) => get_frame!(i32, <4),
(PixelType::UINT32, true) => get_frame!(u32, <4),
(PixelType::FLOAT, true) => get_frame!(f32, <4),
(PixelType::DOUBLE, true) => get_frame!(f64, <8),
(PixelType::INT8, false) => get_frame!(i8, >1),
(PixelType::UINT8, false) => get_frame!(u8, >1),
(PixelType::INT16, false) => get_frame!(i16, >2),
(PixelType::UINT16, false) => get_frame!(u16, >2),
(PixelType::INT32, false) => get_frame!(i32, >4),
(PixelType::UINT32, false) => get_frame!(u32, >4),
(PixelType::FLOAT, false) => get_frame!(f32, >4),
(PixelType::DOUBLE, false) => get_frame!(f64, >8),
}
}
}
mod reader;
mod stats;
mod view;
#[cfg(test)]
mod tests {
use super::*;
use crate::stats::MinMax;
use ndarray::{Array, Array4, Array5, NewAxis};
use rayon::prelude::*;
use crate::axes::Axis;
use crate::reader::{Frame, Reader};
use anyhow::Result;
use ndarray::{s, Array2};
fn open(file: &str) -> Result<Reader> {
let path = std::env::current_dir()?
.join("tests")
@@ -413,4 +103,160 @@ mod tests {
println!("{}", xml);
Ok(())
}
#[test]
fn view() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let a = view.slice(s![0, 5, 0, .., ..])?;
let b = reader.get_frame(0, 5, 0)?;
let c: Array2<isize> = a.try_into()?;
let d: Array2<isize> = b.try_into()?;
assert_eq!(c, d);
Ok(())
}
#[test]
fn view_shape() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let a = view.slice(s![0, ..5, 0, .., 100..200])?;
let shape = a.shape();
assert_eq!(shape, vec![5, 1024, 100]);
Ok(())
}
#[test]
fn view_new_axis() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let a = Array5::<u8>::zeros((1, 9, 1, 1024, 1024));
let a = a.slice(s![0, ..5, 0, NewAxis, 100..200, ..]);
let v = view.slice(s![0, ..5, 0, NewAxis, 100..200, ..])?;
assert_eq!(v.shape(), a.shape());
println!("\nshape: {:?}", v.shape());
let a = a.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis]);
let v = v.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis])?;
assert_eq!(v.shape(), a.shape());
Ok(())
}
#[test]
fn view_permute_axes() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let s = view.shape();
let mut a = Array5::<u8>::zeros((s[0], s[1], s[2], s[3], s[4]));
assert_eq!(view.shape(), a.shape());
let b: Array5<usize> = view.clone().try_into()?;
assert_eq!(b.shape(), a.shape());
let view = view.swap_axes(Axis::C, Axis::Z)?;
a.swap_axes(0, 1);
assert_eq!(view.shape(), a.shape());
let b: Array5<usize> = view.clone().try_into()?;
assert_eq!(b.shape(), a.shape());
println!("{:?}", view.axes());
let view = view.permute_axes(&[Axis::X, Axis::Z, Axis::Y])?;
println!("{:?}", view.axes());
let a = a.permuted_axes([4, 1, 2, 0, 3]);
assert_eq!(view.shape(), a.shape());
let b: Array5<usize> = view.clone().try_into()?;
assert_eq!(b.shape(), a.shape());
Ok(())
}
macro_rules! test_max {
($($name:ident: $b:expr $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let array: Array5<usize> = view.clone().try_into()?;
let view = view.max_proj($b)?;
let a: Array4<usize> = view.clone().try_into()?;
let b = array.max($b)?;
assert_eq!(a.shape(), b.shape());
assert_eq!(a, b);
Ok(())
}
)*
};
}
test_max! {
max_c: 0
max_z: 1
max_t: 2
max_y: 3
max_x: 4
}
macro_rules! test_index {
($($name:ident: $b:expr $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let view = reader.view();
let v4: Array<usize, _> = view.slice($b)?.try_into()?;
let a5: Array5<usize> = reader.view().try_into()?;
let a4 = a5.slice($b).to_owned();
assert_eq!(a4, v4);
Ok(())
}
)*
};
}
test_index! {
index_0: s![.., .., .., .., ..]
index_1: s![0, .., .., .., ..]
index_2: s![.., 0, .., .., ..]
index_3: s![.., .., 0, .., ..]
index_4: s![.., .., .., 0, ..]
index_5: s![.., .., .., .., 0]
index_6: s![0, 0, .., .., ..]
index_7: s![0, .., 0, .., ..]
index_8: s![0, .., .., 0, ..]
index_9: s![0, .., .., .., 0]
index_a: s![.., 0, 0, .., ..]
index_b: s![.., 0, .., 0, ..]
index_c: s![.., 0, .., .., 0]
index_d: s![.., .., 0, 0, ..]
index_e: s![.., .., 0, .., 0]
index_f: s![.., .., .., 0, 0]
index_g: s![0, 0, 0, .., ..]
index_h: s![0, 0, .., 0, ..]
index_i: s![0, 0, .., .., 0]
index_j: s![0, .., 0, 0, ..]
index_k: s![0, .., 0, .., 0]
index_l: s![0, .., .., 0, 0]
index_m: s![0, 0, 0, 0, ..]
index_n: s![0, 0, 0, .., 0]
index_o: s![0, 0, .., 0, 0]
index_p: s![0, .., 0, 0, 0]
index_q: s![.., 0, 0, 0, 0]
index_r: s![0, 0, 0, 0, 0]
}
#[test]
fn dyn_view() -> Result<()> {
let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif";
let reader = open(file)?;
let a = reader.view().into_dyn();
let b = a.max_proj(1)?;
let c = b.slice(s![0, 0, .., ..])?;
let d = c.as_array::<usize>()?;
assert_eq!(d.shape(), [1024, 1024]);
Ok(())
}
}

605
src/py.rs
View File

@@ -1,64 +1,600 @@
use crate::axes::Axis;
use crate::bioformats::download_bioformats;
use crate::{Frame, Reader};
use numpy::ToPyArray;
use crate::reader::{PixelType, Reader};
use crate::view::{Item, View};
use anyhow::{anyhow, Result};
use ndarray::{Ix0, IxDyn, SliceInfoElem};
use numpy::IntoPyArray;
use pyo3::prelude::*;
use pyo3::types::{PyEllipsis, PyInt, PyList, PySlice, PySliceMethods, PyString, PyTuple};
use pyo3::IntoPyObjectExt;
use serde::{Deserialize, Serialize};
use serde_json::{from_str, to_string};
use std::path::PathBuf;
#[pyclass(subclass)]
#[pyo3(name = "Reader")]
#[derive(Debug)]
struct PyReader {
reader: Reader,
#[pyclass(module = "ndbioimage.ndbioimage_rs")]
struct ViewConstructor;
#[pymethods]
impl ViewConstructor {
#[new]
fn new() -> Self {
Self
}
fn __getstate__(&self) -> (u8,) {
(0,)
}
fn __setstate__(&self, _state: (u8,)) {}
#[staticmethod]
fn __call__(state: String) -> PyResult<PyView> {
if let Ok(new) = from_str(&state) {
Ok(new)
} else {
Err(anyhow!("cannot parse state").into())
}
}
}
#[pyclass(subclass, module = "ndbioimage.ndbioimage_rs")]
#[pyo3(name = "View")]
#[derive(Debug, Serialize, Deserialize)]
struct PyView {
view: View<IxDyn>,
dtype: PixelType,
}
#[pymethods]
impl PyReader {
impl PyView {
#[new]
fn new(path: &str, series: usize) -> PyResult<Self> {
#[pyo3(signature = (path, series = 0, dtype = "uint16"))]
/// new view on a file at path, open series #, open as dtype: (u)int(8/16/32) or float(32/64)
fn new(path: &str, series: usize, dtype: &str) -> PyResult<Self> {
let mut path = PathBuf::from(path);
if path.is_dir() {
for file in path.read_dir()? {
if let Ok(f) = file {
let p = f.path();
if f.path().is_file() & (p.extension() == Some("tif".as_ref())) {
path = p;
break;
}
for file in path.read_dir()?.flatten() {
let p = file.path();
if file.path().is_file() & (p.extension() == Some("tif".as_ref())) {
path = p;
break;
}
}
}
Ok(PyReader {
reader: Reader::new(&path, series as i32)?,
Ok(Self {
view: Reader::new(&path, series as i32)?.view().into_dyn(),
dtype: dtype.parse()?,
})
}
/// close the file: does nothing as this is handled automatically
fn close(&self) -> PyResult<()> {
Ok(())
}
fn copy(&self) -> PyView {
PyView {
view: self.view.clone(),
dtype: self.dtype.clone(),
}
}
/// slice the view and return a new view or a single number
fn __getitem__<'py>(
&self,
py: Python<'py>,
n: Bound<'py, PyAny>,
) -> PyResult<Bound<'py, PyAny>> {
let slice: Vec<_> = if n.is_instance_of::<PyTuple>() {
n.downcast_into::<PyTuple>()?.into_iter().collect()
} else if n.is_instance_of::<PyList>() {
n.downcast_into::<PyList>()?.into_iter().collect()
} else {
vec![n]
};
let mut new_slice = Vec::new();
let mut ellipsis = None;
let shape = self.view.shape();
for (i, (s, t)) in slice.iter().zip(shape.iter()).enumerate() {
if s.is_instance_of::<PyInt>() {
new_slice.push(SliceInfoElem::Index(
s.downcast::<PyInt>()?.extract::<isize>()?,
));
} else if s.is_instance_of::<PySlice>() {
let u = s.downcast::<PySlice>()?.indices(*t as isize)?;
new_slice.push(SliceInfoElem::Slice {
start: u.start,
end: Some(u.stop),
step: u.step,
});
} else if s.is_instance_of::<PyEllipsis>() {
if ellipsis.is_some() {
return Err(anyhow!("cannot have more than one ellipsis").into());
}
let _ = ellipsis.insert(i);
} else {
return Err(anyhow!("cannot convert {:?} to slice", s).into());
}
}
if new_slice.len() > shape.len() {
return Err(anyhow!(
"got more indices ({}) than dimensions ({})",
new_slice.len(),
shape.len()
)
.into());
}
while new_slice.len() < shape.len() {
if let Some(i) = ellipsis {
new_slice.insert(
i,
SliceInfoElem::Slice {
start: 0,
end: None,
step: 1,
},
)
} else {
new_slice.push(SliceInfoElem::Slice {
start: 0,
end: None,
step: 1,
})
}
}
let view = self.view.slice(new_slice.as_slice())?;
if view.ndim() == 0 {
Ok(match self.dtype {
PixelType::I8 => view
.into_dimensionality::<Ix0>()?
.item::<i8>()?
.into_pyobject(py)?
.into_any(),
PixelType::U8 => view
.into_dimensionality::<Ix0>()?
.item::<u8>()?
.into_pyobject(py)?
.into_any(),
PixelType::I16 => view
.into_dimensionality::<Ix0>()?
.item::<i16>()?
.into_pyobject(py)?
.into_any(),
PixelType::U16 => view
.into_dimensionality::<Ix0>()?
.item::<u16>()?
.into_pyobject(py)?
.into_any(),
PixelType::I32 => view
.into_dimensionality::<Ix0>()?
.item::<i32>()?
.into_pyobject(py)?
.into_any(),
PixelType::U32 => view
.into_dimensionality::<Ix0>()?
.item::<u32>()?
.into_pyobject(py)?
.into_any(),
PixelType::F32 => view
.into_dimensionality::<Ix0>()?
.item::<f32>()?
.into_pyobject(py)?
.into_any(),
PixelType::F64 => view
.into_dimensionality::<Ix0>()?
.item::<f64>()?
.into_pyobject(py)?
.into_any(),
PixelType::I64 => view
.into_dimensionality::<Ix0>()?
.item::<i64>()?
.into_pyobject(py)?
.into_any(),
PixelType::U64 => view
.into_dimensionality::<Ix0>()?
.item::<u64>()?
.into_pyobject(py)?
.into_any(),
PixelType::I128 => view
.into_dimensionality::<Ix0>()?
.item::<i128>()?
.into_pyobject(py)?
.into_any(),
PixelType::U128 => view
.into_dimensionality::<Ix0>()?
.item::<u128>()?
.into_pyobject(py)?
.into_any(),
PixelType::F128 => view
.into_dimensionality::<Ix0>()?
.item::<f64>()?
.into_pyobject(py)?
.into_any(),
})
} else {
PyView {
view,
dtype: self.dtype.clone(),
}
.into_bound_py_any(py)
}
}
fn __reduce__(&self) -> PyResult<(ViewConstructor, (String,))> {
if let Ok(s) = to_string(self) {
Ok((ViewConstructor, (s,)))
} else {
Err(anyhow!("cannot get state").into())
}
}
/// retrieve a single frame at czt, sliced accordingly
fn get_frame<'py>(
&self,
py: Python<'py>,
c: usize,
z: usize,
t: usize,
c: isize,
z: isize,
t: isize,
) -> PyResult<Bound<'py, PyAny>> {
Ok(match self.reader.get_frame(c, z, t)? {
Frame::INT8(arr) => arr.to_pyarray(py).into_any(),
Frame::UINT8(arr) => arr.to_pyarray(py).into_any(),
Frame::INT16(arr) => arr.to_pyarray(py).into_any(),
Frame::UINT16(arr) => arr.to_pyarray(py).into_any(),
Frame::INT32(arr) => arr.to_pyarray(py).into_any(),
Frame::UINT32(arr) => arr.to_pyarray(py).into_any(),
Frame::FLOAT(arr) => arr.to_pyarray(py).into_any(),
Frame::DOUBLE(arr) => arr.to_pyarray(py).into_any(),
Ok(match self.dtype {
PixelType::I8 => self
.view
.get_frame::<i8>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::U8 => self
.view
.get_frame::<u8>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::I16 => self
.view
.get_frame::<i16>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::U16 => self
.view
.get_frame::<u16>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::I32 => self
.view
.get_frame::<i32>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::U32 => self
.view
.get_frame::<u32>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::F32 => self
.view
.get_frame::<f32>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::F64 => self
.view
.get_frame::<f64>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::I64 => self
.view
.get_frame::<i64>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::U64 => self
.view
.get_frame::<u64>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::I128 => self
.view
.get_frame::<i64>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::U128 => self
.view
.get_frame::<u64>(c, z, t)?
.into_pyarray(py)
.into_any(),
PixelType::F128 => self
.view
.get_frame::<f64>(c, z, t)?
.into_pyarray(py)
.into_any(),
})
}
/// retrieve the ome metadata as an xml string
fn get_ome_xml(&self) -> PyResult<String> {
Ok(self.reader.get_ome_xml()?)
Ok(self.view.get_ome_xml()?)
}
fn close(&mut self) -> PyResult<()> {
self.reader.close()?;
/// the file path
#[getter]
fn path(&self) -> PyResult<String> {
Ok(self.view.path.display().to_string())
}
/// the series in the file
#[getter]
fn series(&self) -> PyResult<i32> {
Ok(self.view.series)
}
/// the axes in the view
#[getter]
fn axes(&self) -> Vec<String> {
self.view
.axes()
.iter()
.map(|a| format!("{:?}", a))
.collect()
}
/// the shape of the view
#[getter]
fn shape(&self) -> Vec<usize> {
self.view.shape()
}
#[getter]
fn slice(&self) -> PyResult<Vec<String>> {
Ok(self
.view
.get_slice()
.iter()
.map(|s| format!("{:#?}", s))
.collect())
}
/// the number of pixels in the view
#[getter]
fn size(&self) -> usize {
self.view.size()
}
/// the number of dimensions in the view
#[getter]
fn ndim(&self) -> usize {
self.view.ndim()
}
/// find the position of an axis
#[pyo3(text_signature = "axis: str | int")]
fn get_ax(&self, axis: Bound<'_, PyAny>) -> PyResult<usize> {
if axis.is_instance_of::<PyString>() {
let axis = axis
.downcast_into::<PyString>()?
.extract::<String>()?
.parse::<Axis>()?;
Ok(self
.view
.axes()
.iter()
.position(|a| *a == axis)
.ok_or_else(|| anyhow!("cannot find axis {:?}", axis))?)
} else if axis.is_instance_of::<PyInt>() {
Ok(axis.downcast_into::<PyInt>()?.extract::<usize>()?)
} else {
Err(anyhow!("cannot convert to axis").into())
}
}
/// swap two axes
#[pyo3(text_signature = "ax0: str | int, ax1: str | int")]
fn swap_axes(&self, ax0: Bound<'_, PyAny>, ax1: Bound<'_, PyAny>) -> PyResult<Self> {
let ax0 = self.get_ax(ax0)?;
let ax1 = self.get_ax(ax1)?;
let view = self.view.swap_axes(ax0, ax1)?;
Ok(PyView {
view,
dtype: self.dtype.clone(),
})
}
/// permute the order of the axes
#[pyo3(signature = (axes = None), text_signature = "axes: list[str | int] = None")]
fn transpose(&self, axes: Option<Vec<Bound<'_, PyAny>>>) -> PyResult<Self> {
let view = if let Some(axes) = axes {
let ax = axes
.into_iter()
.map(|a| self.get_ax(a))
.collect::<Result<Vec<_>, _>>()?;
self.view.permute_axes(&ax)?
} else {
self.view.transpose()?
};
Ok(PyView {
view,
dtype: self.dtype.clone(),
})
}
/// collect data into a numpy array
fn as_array<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
Ok(match self.dtype {
PixelType::I8 => self.view.as_array_dyn::<i8>()?.into_pyarray(py).into_any(),
PixelType::U8 => self.view.as_array_dyn::<u8>()?.into_pyarray(py).into_any(),
PixelType::I16 => self.view.as_array_dyn::<i16>()?.into_pyarray(py).into_any(),
PixelType::U16 => self.view.as_array_dyn::<u16>()?.into_pyarray(py).into_any(),
PixelType::I32 => self.view.as_array_dyn::<i32>()?.into_pyarray(py).into_any(),
PixelType::U32 => self.view.as_array_dyn::<u32>()?.into_pyarray(py).into_any(),
PixelType::F32 => self.view.as_array_dyn::<f32>()?.into_pyarray(py).into_any(),
PixelType::F64 => self.view.as_array_dyn::<f64>()?.into_pyarray(py).into_any(),
PixelType::I64 => self.view.as_array_dyn::<i64>()?.into_pyarray(py).into_any(),
PixelType::U64 => self.view.as_array_dyn::<u64>()?.into_pyarray(py).into_any(),
PixelType::I128 => self.view.as_array_dyn::<i64>()?.into_pyarray(py).into_any(),
PixelType::U128 => self.view.as_array_dyn::<u64>()?.into_pyarray(py).into_any(),
PixelType::F128 => self.view.as_array_dyn::<f64>()?.into_pyarray(py).into_any(),
})
}
/// change the data type of the view: (u)int(8/16/32) or float(32/64)
fn as_type(&self, dtype: String) -> PyResult<Self> {
Ok(PyView {
view: self.view.clone(),
dtype: dtype.parse()?,
})
}
#[getter]
fn get_dtype(&self) -> PyResult<&str> {
Ok(match self.dtype {
PixelType::I8 => "int8",
PixelType::U8 => "uint8",
PixelType::I16 => "int16",
PixelType::U16 => "uint16",
PixelType::I32 => "int32",
PixelType::U32 => "uint32",
PixelType::F32 => "float32",
PixelType::F64 => "float64",
PixelType::I64 => "int64",
PixelType::U64 => "uint64",
PixelType::I128 => "int128",
PixelType::U128 => "uint128",
PixelType::F128 => "float128",
})
}
#[setter]
fn set_dtype(&mut self, dtype: String) -> PyResult<()> {
self.dtype = dtype.parse()?;
Ok(())
}
/// get the maximum overall or along a given axis
#[pyo3(signature = (axis = None), text_signature = "axis: str | int")]
fn max<'py>(
&self,
py: Python<'py>,
axis: Option<Bound<'py, PyAny>>,
) -> PyResult<Bound<'py, PyAny>> {
if let Some(axis) = axis {
PyView {
dtype: self.dtype.clone(),
view: self.view.max_proj(self.get_ax(axis)?)?,
}
.into_bound_py_any(py)
} else {
Ok(match self.dtype {
PixelType::I8 => self.view.max::<i8>()?.into_pyobject(py)?.into_any(),
PixelType::U8 => self.view.max::<u8>()?.into_pyobject(py)?.into_any(),
PixelType::I16 => self.view.max::<i16>()?.into_pyobject(py)?.into_any(),
PixelType::U16 => self.view.max::<u16>()?.into_pyobject(py)?.into_any(),
PixelType::I32 => self.view.max::<i32>()?.into_pyobject(py)?.into_any(),
PixelType::U32 => self.view.max::<u32>()?.into_pyobject(py)?.into_any(),
PixelType::F32 => self.view.max::<f32>()?.into_pyobject(py)?.into_any(),
PixelType::F64 => self.view.max::<f64>()?.into_pyobject(py)?.into_any(),
PixelType::I64 => self.view.max::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U64 => self.view.max::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::I128 => self.view.max::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U128 => self.view.max::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::F128 => self.view.max::<f64>()?.into_pyobject(py)?.into_any(),
})
}
}
/// get the minimum overall or along a given axis
#[pyo3(signature = (axis = None), text_signature = "axis: str | int")]
fn min<'py>(
&self,
py: Python<'py>,
axis: Option<Bound<'py, PyAny>>,
) -> PyResult<Bound<'py, PyAny>> {
if let Some(axis) = axis {
PyView {
dtype: self.dtype.clone(),
view: self.view.min_proj(self.get_ax(axis)?)?,
}
.into_bound_py_any(py)
} else {
Ok(match self.dtype {
PixelType::I8 => self.view.min::<i8>()?.into_pyobject(py)?.into_any(),
PixelType::U8 => self.view.min::<u8>()?.into_pyobject(py)?.into_any(),
PixelType::I16 => self.view.min::<i16>()?.into_pyobject(py)?.into_any(),
PixelType::U16 => self.view.min::<u16>()?.into_pyobject(py)?.into_any(),
PixelType::I32 => self.view.min::<i32>()?.into_pyobject(py)?.into_any(),
PixelType::U32 => self.view.min::<u32>()?.into_pyobject(py)?.into_any(),
PixelType::F32 => self.view.min::<f32>()?.into_pyobject(py)?.into_any(),
PixelType::F64 => self.view.min::<f64>()?.into_pyobject(py)?.into_any(),
PixelType::I64 => self.view.min::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U64 => self.view.min::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::I128 => self.view.min::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U128 => self.view.min::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::F128 => self.view.min::<f64>()?.into_pyobject(py)?.into_any(),
})
}
}
/// get the mean overall or along a given axis
#[pyo3(signature = (axis = None), text_signature = "axis: str | int")]
fn mean<'py>(
&self,
py: Python<'py>,
axis: Option<Bound<'py, PyAny>>,
) -> PyResult<Bound<'py, PyAny>> {
if let Some(axis) = axis {
let dtype = if let PixelType::F32 = self.dtype {
PixelType::F32
} else {
PixelType::F64
};
PyView {
dtype,
view: self.view.mean_proj(self.get_ax(axis)?)?,
}
.into_bound_py_any(py)
} else {
Ok(match self.dtype {
PixelType::F32 => self.view.mean::<f32>()?.into_pyobject(py)?.into_any(),
_ => self.view.mean::<f64>()?.into_pyobject(py)?.into_any(),
})
}
}
/// get the sum overall or along a given axis
#[pyo3(signature = (axis = None), text_signature = "axis: str | int")]
fn sum<'py>(
&self,
py: Python<'py>,
axis: Option<Bound<'py, PyAny>>,
) -> PyResult<Bound<'py, PyAny>> {
let dtype = match self.dtype {
PixelType::I8 => PixelType::I16,
PixelType::U8 => PixelType::U16,
PixelType::I16 => PixelType::I32,
PixelType::U16 => PixelType::U32,
PixelType::I32 => PixelType::I64,
PixelType::U32 => PixelType::U64,
PixelType::F32 => PixelType::F32,
PixelType::F64 => PixelType::F64,
PixelType::I64 => PixelType::I128,
PixelType::U64 => PixelType::U128,
PixelType::I128 => PixelType::I128,
PixelType::U128 => PixelType::U128,
PixelType::F128 => PixelType::F128,
};
if let Some(axis) = axis {
PyView {
dtype,
view: self.view.sum_proj(self.get_ax(axis)?)?,
}
.into_bound_py_any(py)
} else {
Ok(match self.dtype {
PixelType::F32 => self.view.sum::<f32>()?.into_pyobject(py)?.into_any(),
PixelType::F64 => self.view.sum::<f64>()?.into_pyobject(py)?.into_any(),
PixelType::I64 => self.view.sum::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U64 => self.view.sum::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::I128 => self.view.sum::<i64>()?.into_pyobject(py)?.into_any(),
PixelType::U128 => self.view.sum::<u64>()?.into_pyobject(py)?.into_any(),
PixelType::F128 => self.view.sum::<f64>()?.into_pyobject(py)?.into_any(),
_ => self.view.sum::<i64>()?.into_pyobject(py)?.into_any(),
})
}
}
}
pub(crate) fn ndbioimage_file() -> anyhow::Result<PathBuf> {
@@ -75,7 +611,8 @@ pub(crate) fn ndbioimage_file() -> anyhow::Result<PathBuf> {
#[pymodule]
#[pyo3(name = "ndbioimage_rs")]
fn ndbioimage_rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
m.add_class::<PyReader>()?;
m.add_class::<PyView>()?;
m.add_class::<ViewConstructor>()?;
#[pyfn(m)]
#[pyo3(name = "download_bioformats")]

442
src/reader.rs Normal file
View File

@@ -0,0 +1,442 @@
use crate::axes::Axis;
use crate::bioformats;
use crate::bioformats::{DebugTools, ImageReader, MetadataTools};
use crate::view::View;
use anyhow::{anyhow, Error, Result};
use ndarray::{s, Array2, Ix5};
use num::{FromPrimitive, Zero};
use serde::{Deserialize, Serialize};
use std::any::type_name;
use std::fmt::Debug;
use std::ops::Deref;
use std::path::{Path, PathBuf};
use std::str::FromStr;
use std::sync::Arc;
use thread_local::ThreadLocal;
/// Pixel types (u)int(8/16/32) or float(32/64), (u/i)(64/128) are not included in bioformats
#[allow(clippy::upper_case_acronyms)]
#[derive(Clone, Debug, Serialize, Deserialize)]
pub enum PixelType {
I8,
U8,
I16,
U16,
I32,
U32,
F32,
F64,
I64,
U64,
I128,
U128,
F128,
}
impl TryFrom<i32> for PixelType {
type Error = Error;
fn try_from(value: i32) -> Result<Self, Self::Error> {
match value {
0 => Ok(PixelType::I8),
1 => Ok(PixelType::U8),
2 => Ok(PixelType::I16),
3 => Ok(PixelType::U16),
4 => Ok(PixelType::I32),
5 => Ok(PixelType::U32),
6 => Ok(PixelType::F32),
7 => Ok(PixelType::F64),
8 => Ok(PixelType::I64),
9 => Ok(PixelType::U64),
10 => Ok(PixelType::I128),
11 => Ok(PixelType::U128),
12 => Ok(PixelType::F128),
_ => Err(anyhow::anyhow!("Unknown pixel type {}", value)),
}
}
}
impl FromStr for PixelType {
type Err = Error;
fn from_str(s: &str) -> Result<Self, Self::Err> {
match s.to_lowercase().as_str() {
"int8" | "i8" => Ok(PixelType::I8),
"uint8" | "u8" => Ok(PixelType::U8),
"int16" | "i16" => Ok(PixelType::I16),
"uint16" | "u16" => Ok(PixelType::U16),
"int32" | "i32" => Ok(PixelType::I32),
"uint32" | "u32" => Ok(PixelType::U32),
"float" | "f32" | "float32" => Ok(PixelType::F32),
"double" | "f64" | "float64" => Ok(PixelType::F64),
"int64" | "i64" => Ok(PixelType::I64),
"uint64" | "u64" => Ok(PixelType::U64),
"int128" | "i128" => Ok(PixelType::I128),
"uint128" | "u128" => Ok(PixelType::U128),
"extended" | "f128" => Ok(PixelType::F128),
_ => Err(anyhow::anyhow!("Unknown pixel type {}", s)),
}
}
}
/// Struct containing frame data in one of eight pixel types. Cast to `Array2<T>` using try_into.
#[allow(clippy::upper_case_acronyms)]
#[derive(Clone, Debug)]
pub enum Frame {
I8(Array2<i8>),
U8(Array2<u8>),
I16(Array2<i16>),
U16(Array2<u16>),
I32(Array2<i32>),
U32(Array2<u32>),
F32(Array2<f32>),
F64(Array2<f64>),
I64(Array2<i64>),
U64(Array2<u64>),
I128(Array2<i128>),
U128(Array2<u128>),
F128(Array2<f64>), // f128 is nightly
}
macro_rules! impl_frame_cast {
($($t:tt: $s:ident $(,)?)*) => {
$(
impl From<Array2<$t>> for Frame {
fn from(value: Array2<$t>) -> Self {
Frame::$s(value)
}
}
)*
};
}
impl_frame_cast! {
u8: U8
i8: I8
i16: I16
u16: U16
i32: I32
u32: U32
f32: F32
f64: F64
i64: I64
u64: U64
i128: I128
u128: U128
}
#[cfg(target_pointer_width = "32")]
impl_frame_cast! {
usize: UINT32
isize: INT32
}
impl<T> TryInto<Array2<T>> for Frame
where
T: FromPrimitive + Zero + 'static,
{
type Error = Error;
fn try_into(self) -> Result<Array2<T>, Self::Error> {
let mut err = Ok(());
let arr = match self {
Frame::I8(v) => v.mapv_into_any(|x| {
T::from_i8(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::U8(v) => v.mapv_into_any(|x| {
T::from_u8(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::I16(v) => v.mapv_into_any(|x| {
T::from_i16(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::U16(v) => v.mapv_into_any(|x| {
T::from_u16(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::I32(v) => v.mapv_into_any(|x| {
T::from_i32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::U32(v) => v.mapv_into_any(|x| {
T::from_u32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::F32(v) => v.mapv_into_any(|x| {
T::from_f32(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::F64(v) | Frame::F128(v) => v.mapv_into_any(|x| {
T::from_f64(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::I64(v) => v.mapv_into_any(|x| {
T::from_i64(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::U64(v) => v.mapv_into_any(|x| {
T::from_u64(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::I128(v) => v.mapv_into_any(|x| {
T::from_i128(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
Frame::U128(v) => v.mapv_into_any(|x| {
T::from_u128(x).unwrap_or_else(|| {
err = Err(anyhow!("cannot convert {} into {}", x, type_name::<T>()));
T::zero()
})
}),
};
match err {
Err(err) => Err(err),
Ok(()) => Ok(arr),
}
}
}
/// Reader interface to file. Use get_frame to get data.
#[derive(Serialize, Deserialize)]
pub struct Reader {
#[serde(skip)]
image_reader: ThreadLocal<ImageReader>,
/// path to file
pub path: PathBuf,
/// which (if more than 1) of the series in the file to open
pub series: i32,
/// size x (horizontal)
pub size_x: usize,
/// size y (vertical)
pub size_y: usize,
/// size c (# channels)
pub size_c: usize,
/// size z (# slices)
pub size_z: usize,
/// size t (# time/frames)
pub size_t: usize,
/// pixel type ((u)int(8/16/32) or float(32/64))
pub pixel_type: PixelType,
little_endian: bool,
}
impl Deref for Reader {
type Target = ImageReader;
fn deref(&self) -> &Self::Target {
self.image_reader.get_or(|| {
let reader = ImageReader::new().unwrap();
let meta_data_tools = MetadataTools::new().unwrap();
let ome_meta = meta_data_tools.create_ome_xml_metadata().unwrap();
reader.set_metadata_store(ome_meta).unwrap();
reader.set_id(self.path.to_str().unwrap()).unwrap();
reader.set_series(self.series).unwrap();
reader
})
}
}
impl Clone for Reader {
fn clone(&self) -> Self {
Reader::new(&self.path, self.series).unwrap()
}
}
impl Debug for Reader {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Reader")
.field("path", &self.path)
.field("series", &self.series)
.field("size_x", &self.size_x)
.field("size_y", &self.size_y)
.field("size_c", &self.size_c)
.field("size_z", &self.size_z)
.field("size_t", &self.size_t)
.field("pixel_type", &self.pixel_type)
.field("little_endian", &self.little_endian)
.finish()
}
}
impl Reader {
/// Create a new reader for the image file at a path, and open series #.
pub fn new(path: &Path, series: i32) -> Result<Self> {
DebugTools::set_root_level("ERROR")?;
let mut reader = Reader {
image_reader: ThreadLocal::default(),
path: PathBuf::from(path),
series,
size_x: 0,
size_y: 0,
size_c: 0,
size_z: 0,
size_t: 0,
pixel_type: PixelType::I8,
little_endian: false,
};
reader.size_x = reader.get_size_x()? as usize;
reader.size_y = reader.get_size_y()? as usize;
reader.size_c = reader.get_size_c()? as usize;
reader.size_z = reader.get_size_z()? as usize;
reader.size_t = reader.get_size_t()? as usize;
reader.pixel_type = PixelType::try_from(reader.get_pixel_type()?)?;
reader.little_endian = reader.is_little_endian()?;
Ok(reader)
}
/// Get ome metadata as xml string
pub fn get_ome_xml(&self) -> Result<String> {
self.ome_xml()
}
fn deinterleave(&self, bytes: Vec<u8>, channel: usize) -> Result<Vec<u8>> {
let chunk_size = match self.pixel_type {
PixelType::I8 => 1,
PixelType::U8 => 1,
PixelType::I16 => 2,
PixelType::U16 => 2,
PixelType::I32 => 4,
PixelType::U32 => 4,
PixelType::F32 => 4,
PixelType::F64 => 8,
PixelType::I64 => 8,
PixelType::U64 => 8,
PixelType::I128 => 16,
PixelType::U128 => 16,
PixelType::F128 => 8,
};
Ok(bytes
.chunks(chunk_size)
.skip(channel)
.step_by(self.size_c)
.flat_map(|a| a.to_vec())
.collect())
}
/// Retrieve fame at channel c, slize z and time t.
pub fn get_frame(&self, c: usize, z: usize, t: usize) -> Result<Frame> {
let bytes = if self.is_rgb()? & self.is_interleaved()? {
let index = self.get_index(z as i32, 0, t as i32)?;
self.deinterleave(self.open_bytes(index)?, c)?
} else if self.get_rgb_channel_count()? > 1 {
let channel_separator = bioformats::ChannelSeparator::new(self)?;
let index = channel_separator.get_index(z as i32, c as i32, t as i32)?;
channel_separator.open_bytes(index)?
} else if self.is_indexed()? {
let index = self.get_index(z as i32, 0, t as i32)?;
self.open_bytes(index)?
// TODO: apply LUT
// let _bytes_lut = match self.pixel_type {
// PixelType::INT8 | PixelType::UINT8 => {
// let _lut = self.image_reader.get_8bit_lookup_table()?;
// }
// PixelType::INT16 | PixelType::UINT16 => {
// let _lut = self.image_reader.get_16bit_lookup_table()?;
// }
// _ => {}
// };
} else {
let index = self.get_index(z as i32, c as i32, t as i32)?;
self.open_bytes(index)?
};
self.bytes_to_frame(bytes)
}
fn bytes_to_frame(&self, bytes: Vec<u8>) -> Result<Frame> {
macro_rules! get_frame {
($t:tt, <$n:expr) => {
Ok(Frame::from(Array2::from_shape_vec(
(self.size_y, self.size_x),
bytes
.chunks($n)
.map(|x| $t::from_le_bytes(x.try_into().unwrap()))
.collect(),
)?))
};
($t:tt, >$n:expr) => {
Ok(Frame::from(Array2::from_shape_vec(
(self.size_y, self.size_x),
bytes
.chunks($n)
.map(|x| $t::from_be_bytes(x.try_into().unwrap()))
.collect(),
)?))
};
}
match (&self.pixel_type, self.little_endian) {
(PixelType::I8, true) => get_frame!(i8, <1),
(PixelType::U8, true) => get_frame!(u8, <1),
(PixelType::I16, true) => get_frame!(i16, <2),
(PixelType::U16, true) => get_frame!(u16, <2),
(PixelType::I32, true) => get_frame!(i32, <4),
(PixelType::U32, true) => get_frame!(u32, <4),
(PixelType::F32, true) => get_frame!(f32, <4),
(PixelType::F64, true) => get_frame!(f64, <8),
(PixelType::I64, true) => get_frame!(i64, <8),
(PixelType::U64, true) => get_frame!(u64, <8),
(PixelType::I128, true) => get_frame!(i128, <16),
(PixelType::U128, true) => get_frame!(u128, <16),
(PixelType::F128, true) => get_frame!(f64, <8),
(PixelType::I8, false) => get_frame!(i8, >1),
(PixelType::U8, false) => get_frame!(u8, >1),
(PixelType::I16, false) => get_frame!(i16, >2),
(PixelType::U16, false) => get_frame!(u16, >2),
(PixelType::I32, false) => get_frame!(i32, >4),
(PixelType::U32, false) => get_frame!(u32, >4),
(PixelType::F32, false) => get_frame!(f32, >4),
(PixelType::F64, false) => get_frame!(f64, >8),
(PixelType::I64, false) => get_frame!(i64, >8),
(PixelType::U64, false) => get_frame!(u64, >8),
(PixelType::I128, false) => get_frame!(i128, >16),
(PixelType::U128, false) => get_frame!(u128, >16),
(PixelType::F128, false) => get_frame!(f64, >8),
}
}
/// get a slicable view on the image file
pub fn view(&self) -> View<Ix5> {
let slice = s![
0..self.size_c,
0..self.size_z,
0..self.size_t,
0..self.size_y,
0..self.size_x
];
View::new(
Arc::new(self.clone()),
slice.as_ref().to_vec(),
vec![Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X],
)
}
}
impl Drop for Reader {
fn drop(&mut self) {
let _ = self.close();
}
}

201
src/stats.rs Normal file
View File

@@ -0,0 +1,201 @@
use anyhow::{anyhow, Result};
use ndarray::{Array, ArrayD, ArrayView, Axis, Dimension, RemoveAxis};
/// a trait to define the min, max, sum and mean operations along an axis
pub trait MinMax {
type Output;
fn max(self, axis: usize) -> Result<Self::Output>;
fn min(self, axis: usize) -> Result<Self::Output>;
fn sum(self, axis: usize) -> Result<Self::Output>;
fn mean(self, axis: usize) -> Result<Self::Output>;
}
macro_rules! impl_frame_stats_float_view {
($($t:tt),+ $(,)?) => {
$(
impl<D> MinMax for ArrayView<'_, $t, D>
where
D: Dimension + RemoveAxis,
{
type Output = Array<$t, D::Smaller>;
fn max(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| {
x.iter()
.fold($t::NEG_INFINITY, |prev, curr| prev.max(*curr))
})
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn min(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| {
x.iter()
.fold($t::INFINITY, |prev, curr| prev.min(*curr))
})
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn sum(self, axis: usize) -> Result<Self::Output> {
Ok(self.sum_axis(Axis(axis)))
}
fn mean(self, axis: usize) -> Result<Self::Output> {
self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean"))
}
}
)*
};
}
macro_rules! impl_frame_stats_int_view {
($($t:tt),+ $(,)?) => {
$(
impl<D> MinMax for ArrayView<'_, $t, D>
where
D: Dimension + RemoveAxis,
{
type Output = Array<$t, D::Smaller>;
fn max(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| *x.iter().max().unwrap())
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn min(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| *x.iter().min().unwrap())
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn sum(self, axis: usize) -> Result<Self::Output> {
Ok(self.sum_axis(Axis(axis)))
}
fn mean(self, axis: usize) -> Result<Self::Output> {
self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean"))
}
}
)*
};
}
macro_rules! impl_frame_stats_float {
($($t:tt),+ $(,)?) => {
$(
impl<D> MinMax for Array<$t, D>
where
D: Dimension + RemoveAxis,
{
type Output = Array<$t, D::Smaller>;
fn max(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| {
x.iter()
.fold($t::NEG_INFINITY, |prev, curr| prev.max(*curr))
})
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn min(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| {
x.iter()
.fold($t::INFINITY, |prev, curr| prev.min(*curr))
})
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn sum(self, axis: usize) -> Result<Self::Output> {
Ok(self.sum_axis(Axis(axis)))
}
fn mean(self, axis: usize) -> Result<Self::Output> {
self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean"))
}
}
)*
};
}
macro_rules! impl_frame_stats_int {
($($t:tt),+ $(,)?) => {
$(
impl<D> MinMax for Array<$t, D>
where
D: Dimension + RemoveAxis,
{
type Output = Array<$t, D::Smaller>;
fn max(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| *x.iter().max().unwrap())
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn min(self, axis: usize) -> Result<Self::Output> {
let a: Vec<_> = self
.lanes(Axis(axis))
.into_iter()
.map(|x| *x.iter().min().unwrap())
.collect();
let mut shape = self.shape().to_vec();
shape.remove(axis);
Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?)
}
fn sum(self, axis: usize) -> Result<Self::Output> {
Ok(self.sum_axis(Axis(axis)))
}
fn mean(self, axis: usize) -> Result<Self::Output> {
self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean"))
}
}
)*
};
}
impl_frame_stats_float_view!(f32, f64);
impl_frame_stats_int_view!(u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize);
impl_frame_stats_float!(f32, f64);
impl_frame_stats_int!(u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize);

821
src/view.rs Normal file
View File

@@ -0,0 +1,821 @@
use crate::axes::{slice_info, Ax, Axis, Operation, Slice, SliceInfoElemDef};
use crate::reader::Reader;
use crate::stats::MinMax;
use anyhow::{anyhow, Error, Result};
use indexmap::IndexMap;
use itertools::iproduct;
use ndarray::{
s, Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn,
SliceArg, SliceInfoElem,
};
use num::{Bounded, FromPrimitive, Zero};
use serde::{Deserialize, Serialize};
use serde_with::serde_as;
use std::any::type_name;
use std::collections::HashMap;
use std::iter::Sum;
use std::marker::PhantomData;
use std::mem::{transmute, transmute_copy};
use std::ops::{AddAssign, Deref, Div};
use std::path::PathBuf;
use std::sync::Arc;
fn idx_bnd(idx: isize, bnd: isize) -> Result<isize> {
if idx < -bnd {
Err(anyhow!("Index {} out of bounds {}", idx, bnd))
} else if idx < 0 {
Ok(bnd - idx)
} else if idx < bnd {
Ok(idx)
} else {
Err(anyhow!("Index {} out of bounds {}", idx, bnd))
}
}
fn slc_bnd(idx: isize, bnd: isize) -> Result<isize> {
if idx < -bnd {
Err(anyhow!("Index {} out of bounds {}", idx, bnd))
} else if idx < 0 {
Ok(bnd - idx)
} else if idx <= bnd {
Ok(idx)
} else {
Err(anyhow!("Index {} out of bounds {}", idx, bnd))
}
}
pub trait Number:
'static + AddAssign + Bounded + Clone + Div<Self, Output = Self> + FromPrimitive + PartialOrd + Zero
{
}
impl<T> Number for T where
T: 'static
+ AddAssign
+ Bounded
+ Clone
+ Div<Self, Output = Self>
+ FromPrimitive
+ PartialOrd
+ Zero
{
}
/// sliceable view on an image file
#[serde_as]
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct View<D: Dimension> {
reader: Arc<Reader>,
#[serde_as(as = "Vec<SliceInfoElemDef>")]
slice: Vec<SliceInfoElem>,
axes: Vec<Axis>,
operations: IndexMap<Axis, Operation>,
dimensionality: PhantomData<D>,
}
impl<D: Dimension> View<D> {
pub(crate) fn new(reader: Arc<Reader>, slice: Vec<SliceInfoElem>, axes: Vec<Axis>) -> Self {
Self {
reader,
slice,
axes,
operations: IndexMap::new(),
dimensionality: PhantomData,
}
}
/// the file path
pub fn path(&self) -> &PathBuf {
&self.reader.path
}
/// the series in the file
pub fn series(&self) -> i32 {
self.reader.series
}
fn with_operations(mut self, operations: IndexMap<Axis, Operation>) -> Self {
self.operations = operations;
self
}
/// change the dimension into a dynamic dimension
pub fn into_dyn(self) -> View<IxDyn> {
View {
reader: self.reader,
slice: self.slice,
axes: self.axes,
operations: self.operations,
dimensionality: PhantomData,
}
}
/// change the dimension into a concrete dimension
pub fn into_dimensionality<D2: Dimension>(self) -> Result<View<D2>> {
if let Some(d) = D2::NDIM {
if d == self.ndim() {
Ok(View {
reader: self.reader,
slice: self.slice,
axes: self.axes,
operations: self.operations,
dimensionality: PhantomData,
})
} else {
Err(anyhow!("Dimensionality mismatch: {} != {}", d, self.ndim()))
}
} else {
Ok(View {
reader: self.reader,
slice: self.slice,
axes: self.axes,
operations: self.operations,
dimensionality: PhantomData,
})
}
}
/// the order of the axes, including axes sliced out
pub fn get_axes(&self) -> &[Axis] {
&self.axes
}
/// the slice defining the view
pub fn get_slice(&self) -> &[SliceInfoElem] {
&self.slice
}
/// the axes in the view
pub fn axes(&self) -> Vec<Axis> {
self.axes
.iter()
.zip(self.slice.iter())
.filter_map(|(ax, s)| {
if s.is_index() || self.operations.contains_key(ax) {
None
} else {
Some(*ax)
}
})
.collect()
}
pub(crate) fn op_axes(&self) -> Vec<Axis> {
self.operations.keys().cloned().collect()
}
/// the number of dimensions in the view
pub fn ndim(&self) -> usize {
if let Some(d) = D::NDIM {
d
} else {
self.shape().len()
}
}
/// the number of pixels in the view
pub fn size(&self) -> usize {
self.shape().into_iter().product()
}
/// the shape of the view
pub fn shape(&self) -> Vec<usize> {
let mut shape = Vec::<usize>::new();
for (ax, s) in self.axes.iter().zip(self.slice.iter()) {
match s {
SliceInfoElem::Slice { start, end, step } => {
if !self.operations.contains_key(ax) {
if let Some(e) = end {
shape.push(((e - start).max(0) / step) as usize);
} else {
panic!("slice has no end")
}
}
}
SliceInfoElem::Index(_) => {}
SliceInfoElem::NewAxis => {
if !self.operations.contains_key(ax) {
shape.push(1);
}
}
}
}
shape
}
fn shape_all(&self) -> Vec<usize> {
let mut shape = Vec::<usize>::new();
for s in self.slice.iter() {
match s {
SliceInfoElem::Slice { start, end, step } => {
if let Some(e) = end {
shape.push(((e - start).max(0) / step) as usize);
} else {
panic!("slice has no end")
}
}
_ => shape.push(1),
}
}
shape
}
/// swap two axes
pub fn swap_axes<A: Ax>(&self, axis0: A, axis1: A) -> Result<Self> {
let idx0 = axis0.pos_op(&self.axes, &self.slice, &self.op_axes())?;
let idx1 = axis1.pos_op(&self.axes, &self.slice, &self.op_axes())?;
let mut slice = self.slice.to_vec();
slice.swap(idx0, idx1);
let mut axes = self.axes.clone();
axes.swap(idx0, idx1);
Ok(View::new(self.reader.clone(), slice, axes))
}
/// subset of gives axes will be reordered in given order
pub fn permute_axes<A: Ax>(&self, axes: &[A]) -> Result<Self> {
let idx: Vec<usize> = axes
.iter()
.map(|a| a.pos_op(&self.axes, &self.slice, &self.op_axes()).unwrap())
.collect();
let mut jdx = idx.clone();
jdx.sort();
let mut slice = self.slice.to_vec();
let mut axes = self.axes.clone();
for (&i, j) in idx.iter().zip(jdx) {
slice[j] = self.slice[i];
axes[j] = self.axes[i];
}
Ok(View::new(self.reader.clone(), slice, axes))
}
/// reverse the order of the axes
pub fn transpose(&self) -> Result<Self> {
Ok(View::new(
self.reader.clone(),
self.slice.iter().rev().cloned().collect(),
self.axes.iter().rev().cloned().collect(),
))
}
fn operate<A: Ax>(&self, axis: A, operation: Operation) -> Result<View<D::Smaller>> {
let pos = axis.pos_op(&self.axes, &self.slice, &self.op_axes())?;
let mut operations = self.operations.clone();
operations.insert(self.axes[pos], operation);
Ok(
View::new(self.reader.clone(), self.slice.clone(), self.axes.clone())
.with_operations(operations),
)
}
/// maximum along axis
pub fn max_proj<A: Ax>(&self, axis: A) -> Result<View<D::Smaller>> {
self.operate(axis, Operation::Max)
}
/// minimum along axis
pub fn min_proj<A: Ax>(&self, axis: A) -> Result<View<D::Smaller>> {
self.operate(axis, Operation::Min)
}
/// sum along axis
pub fn sum_proj<A: Ax>(&self, axis: A) -> Result<View<D::Smaller>> {
self.operate(axis, Operation::Sum)
}
/// mean along axis
pub fn mean_proj<A: Ax>(&self, axis: A) -> Result<View<D::Smaller>> {
self.operate(axis, Operation::Mean)
}
/// created a new sliced view
pub fn slice<I>(&self, info: I) -> Result<View<I::OutDim>>
where
I: SliceArg<D>,
{
if self.slice.out_ndim() < info.in_ndim() {
return Err(Error::msg("not enough free dimensions"));
}
let info = info.as_ref();
let mut n_idx = 0;
let mut r_idx = 0;
let mut new_slice = Vec::new();
let mut new_axes = Vec::new();
let reader_slice = self.slice.as_slice();
while (r_idx < reader_slice.len()) & (n_idx < info.len()) {
let n = &info[n_idx];
let r = &reader_slice[r_idx];
let a = &self.axes[r_idx];
if self.operations.contains_key(a) {
new_slice.push(*r);
new_axes.push(*a);
r_idx += 1;
} else {
match (n, r) {
(
SliceInfoElem::Slice {
start: info_start,
end: info_end,
step: info_step,
},
SliceInfoElem::Slice { start, end, step },
) => {
let new_start = start + info_start;
let new_end = match (info_end, end) {
(Some(m), Some(n)) => *n.min(&(start + info_step * m)),
(None, Some(n)) => *n,
_ => panic!("slice has no end"),
};
let new_step = (step * info_step).abs();
new_slice.push(SliceInfoElem::Slice {
start: new_start,
end: Some(new_end),
step: new_step,
});
new_axes.push(*a);
n_idx += 1;
r_idx += 1;
}
(SliceInfoElem::Index(k), SliceInfoElem::Slice { start, end, step }) => {
if *k < 0 {
new_slice.push(SliceInfoElem::Index(end.unwrap_or(0) + step.abs() * k))
} else {
new_slice.push(SliceInfoElem::Index(start + step.abs() * k));
}
new_axes.push(*a);
n_idx += 1;
r_idx += 1;
}
(SliceInfoElem::NewAxis, SliceInfoElem::NewAxis) => {
new_slice.push(SliceInfoElem::NewAxis);
new_slice.push(SliceInfoElem::NewAxis);
new_axes.push(Axis::New);
new_axes.push(Axis::New);
n_idx += 1;
r_idx += 1;
}
(_, SliceInfoElem::NewAxis) => {
new_slice.push(SliceInfoElem::NewAxis);
new_axes.push(Axis::New);
n_idx += 1;
r_idx += 1;
}
(SliceInfoElem::NewAxis, _) => {
new_slice.push(SliceInfoElem::NewAxis);
new_axes.push(Axis::New);
n_idx += 1;
}
(_, SliceInfoElem::Index(k)) => {
new_slice.push(SliceInfoElem::Index(*k));
new_axes.push(*a);
r_idx += 1;
}
}
}
}
debug_assert_eq!(r_idx, reader_slice.len());
while n_idx < info.len() {
debug_assert!(info[n_idx].is_new_axis());
new_slice.push(SliceInfoElem::NewAxis);
new_axes.push(Axis::New);
n_idx += 1;
}
Ok(View::new(self.reader.clone(), new_slice, new_axes)
.with_operations(self.operations.clone()))
}
/// slice, but slice elements are in cztyx order, all cztyx must be given,
/// but axes not present in view will be ignored, view axes are reordered in cztyx order
pub fn slice_cztyx<I>(&self, info: I) -> Result<View<I::OutDim>>
where
I: SliceArg<Ix5>,
{
let axes = self.axes();
let slice: Vec<_> = info
.as_ref()
.iter()
.zip(self.axes.iter())
.filter_map(|(&s, ax)| if axes.contains(ax) { Some(s) } else { None })
.collect();
let new_axes: Vec<_> = [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X]
.into_iter()
.filter(|ax| axes.contains(ax))
.collect();
self.clone()
.into_dyn()
.permute_axes(&new_axes)?
.slice(slice.as_slice())?
.into_dimensionality()
}
/// the pixel intensity at a given index
pub fn item_at<T>(&self, index: &[isize]) -> Result<T>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
let slice: Vec<_> = index.iter().map(|s| SliceInfoElem::Index(*s)).collect();
let view = self.clone().into_dyn().slice(slice.as_slice())?;
let arr = view.as_array()?;
Ok(arr.first().unwrap().clone())
}
/// collect the view into an ndarray
pub fn as_array<T>(&self) -> Result<Array<T, D>>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
Ok(self.as_array_dyn()?.into_dimensionality()?)
}
/// collect the view into a dynamic-dimension ndarray
pub fn as_array_dyn<T>(&self) -> Result<ArrayD<T>>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
let mut op_xy = IndexMap::new();
if let Some((&ax, op)) = self.operations.first() {
if (ax == Axis::X) || (ax == Axis::Y) {
op_xy.insert(ax, op.clone());
if let Some((&ax2, op2)) = self.operations.get_index(1) {
if (ax2 == Axis::X) || (ax2 == Axis::Y) {
op_xy.insert(ax2, op2.clone());
}
}
}
}
let op_czt = if let Some((&ax, op)) = self.operations.get_index(op_xy.len()) {
IndexMap::from([(ax, op.clone())])
} else {
IndexMap::new()
};
let mut shape_out = Vec::new();
let mut slice = Vec::new();
let mut ax_out = Vec::new();
for (s, a) in self.slice.iter().zip(&self.axes) {
match s {
SliceInfoElem::Slice { start, end, step } => {
if let Some(e) = end {
if !op_xy.contains_key(a) && !op_czt.contains_key(a) {
shape_out.push(((e - start).max(0) / step) as usize);
slice.push(SliceInfoElem::Slice {
start: 0,
end: None,
step: 1,
});
ax_out.push(*a);
}
} else {
panic!("slice has no end")
}
}
SliceInfoElem::Index(_) => {}
SliceInfoElem::NewAxis => {
shape_out.push(1);
slice.push(SliceInfoElem::Index(0));
ax_out.push(*a);
}
}
}
let mut slice_reader = vec![Slice::empty(); 5];
let mut xy_dim = 0usize;
let shape = [
self.size_c as isize,
self.size_z as isize,
self.size_t as isize,
self.size_y as isize,
self.size_x as isize,
];
for (s, &axis) in self.slice.iter().zip(&self.axes) {
match axis {
Axis::New => {}
_ => match s {
SliceInfoElem::Slice { start, end, step } => {
if let Axis::X | Axis::Y = axis {
if !op_xy.contains_key(&axis) {
xy_dim += 1;
}
}
slice_reader[axis as usize] = Slice::new(
idx_bnd(*start, shape[axis as usize])?,
slc_bnd(end.unwrap(), shape[axis as usize])?,
*step,
);
}
SliceInfoElem::Index(j) => {
slice_reader[axis as usize] = Slice::new(
idx_bnd(*j, shape[axis as usize])?,
slc_bnd(*j + 1, shape[axis as usize])?,
1,
);
}
SliceInfoElem::NewAxis => panic!("axis cannot be a new axis"),
},
}
}
let xy = [
self.slice[Axis::Y.pos(&self.axes, &self.slice)?],
self.slice[Axis::X.pos(&self.axes, &self.slice)?],
];
let mut array = if let Some((_, op)) = op_czt.first() {
match op {
Operation::Max => {
ArrayD::<T>::from_elem(shape_out.into_dimension(), T::min_value())
}
Operation::Min => {
ArrayD::<T>::from_elem(shape_out.into_dimension(), T::max_value())
}
_ => ArrayD::<T>::zeros(shape_out.into_dimension()),
}
} else {
ArrayD::<T>::zeros(shape_out.into_dimension())
};
let size_c = self.reader.size_c as isize;
let size_z = self.reader.size_z as isize;
let size_t = self.reader.size_t as isize;
let mut axes_out_idx = [None; 5];
for (i, ax) in ax_out.iter().enumerate() {
axes_out_idx[*ax as usize] = Some(i);
}
for (c, z, t) in iproduct!(&slice_reader[0], &slice_reader[1], &slice_reader[2]) {
if let Some(i) = axes_out_idx[0] {
slice[i] = SliceInfoElem::Index(c)
};
if let Some(i) = axes_out_idx[1] {
slice[i] = SliceInfoElem::Index(z)
};
if let Some(i) = axes_out_idx[2] {
slice[i] = SliceInfoElem::Index(t)
};
let frame = self.reader.get_frame(
(c % size_c) as usize,
(z % size_z) as usize,
(t % size_t) as usize,
)?;
let arr_frame: Array2<T> = frame.try_into()?;
let arr_frame = match xy_dim {
0 => {
if op_xy.contains_key(&Axis::X) && op_xy.contains_key(&Axis::Y) {
let xys = slice_info::<Ix2>(&xy)?;
let (&ax0, op0) = op_xy.first().unwrap();
let (&ax1, op1) = op_xy.get_index(1).unwrap();
let a = arr_frame.slice(xys).to_owned();
let b = op0.operate(a, ax0 as usize - 3)?;
let c: &Array1<T> = unsafe { transmute(&b) };
let d = op1.operate(c.to_owned(), ax1 as usize - 3)?;
let e: &Array0<T> = unsafe { transmute(&d) };
e.to_owned().into_dyn()
} else if op_xy.contains_key(&Axis::X) || op_xy.contains_key(&Axis::Y) {
let xys = slice_info::<Ix1>(&xy)?;
let (&ax, op) = op_xy.first().unwrap();
let a = arr_frame.slice(xys).to_owned();
let b = op.operate(a, ax as usize - 3)?;
let c: &Array0<T> = unsafe { transmute(&b) };
c.to_owned().into_dyn()
} else {
let xys = slice_info::<Ix0>(&xy)?;
arr_frame.slice(xys).to_owned().into_dyn()
}
}
1 => {
if op_xy.contains_key(&Axis::X) || op_xy.contains_key(&Axis::Y) {
let xys = slice_info::<Ix2>(&xy)?;
let (&ax, op) = op_xy.first().unwrap();
let a = arr_frame.slice(xys).to_owned();
let b = op.operate(a, ax as usize - 3)?;
let c: &Array1<T> = unsafe { transmute(&b) };
c.to_owned().into_dyn()
} else {
let xys = slice_info::<Ix1>(&xy)?;
arr_frame.slice(xys).to_owned().into_dyn()
}
}
2 => {
let xys = slice_info::<Ix2>(&xy)?;
if axes_out_idx[4] < axes_out_idx[3] {
arr_frame.t().slice(xys).to_owned().into_dyn()
} else {
arr_frame.slice(xys).to_owned().into_dyn()
}
}
_ => {
panic!("xy cannot be 3d or more");
}
};
if let Some((_, op)) = op_czt.first() {
match op {
Operation::Max => {
array
.slice_mut(slice.as_slice())
.zip_mut_with(&arr_frame, |x, y| {
*x = if *x >= *y { x.clone() } else { y.clone() }
});
}
Operation::Min => {
array
.slice_mut(slice.as_slice())
.zip_mut_with(&arr_frame, |x, y| {
*x = if *x < *y { x.clone() } else { y.clone() }
});
}
Operation::Sum => {
array
.slice_mut(slice.as_slice())
.zip_mut_with(&arr_frame, |x, y| *x += y.clone());
}
Operation::Mean => {
array
.slice_mut(slice.as_slice())
.zip_mut_with(&arr_frame, |x, y| *x += y.clone());
}
}
} else {
array.slice_mut(slice.as_slice()).assign(&arr_frame)
}
}
let mut out = Some(array);
let ax_out: HashMap<Axis, usize> = ax_out
.into_iter()
.enumerate()
.map(|(i, a)| (a, i))
.collect();
for (ax, op) in self.operations.iter().skip(op_xy.len() + op_czt.len()) {
if let Some(&idx) = ax_out.get(ax) {
let arr = out.take().unwrap();
let _ = out.insert(unsafe { transmute_copy(&op.operate(arr, idx)?) });
}
}
let mut n = 1;
for (&ax, size) in self.axes.iter().zip(self.shape_all().iter()) {
if let Some(Operation::Mean) = self.operations.get(&ax) {
if (ax == Axis::C) || (ax == Axis::Z) || (ax == Axis::T) {
n *= size;
}
}
}
let array = if n == 1 {
out.take().unwrap()
} else {
let m = T::from_usize(n).unwrap_or_else(|| T::zero());
out.take().unwrap().mapv(|x| x / m.clone())
};
Ok(array)
}
/// retrieve a single frame at czt, sliced accordingly
pub fn get_frame<T>(&self, c: isize, z: isize, t: isize) -> Result<Array2<T>>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
self.slice_cztyx(s![c, z, t, .., ..])?.as_array()
}
fn get_stat<T>(&self, operation: Operation) -> Result<T>
where
T: Number + Sum,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
let arr: ArrayD<T> = self.as_array_dyn()?;
Ok(match operation {
Operation::Max => arr
.flatten()
.into_iter()
.reduce(|a, b| if a > b { a } else { b })
.unwrap_or_else(|| T::min_value()),
Operation::Min => arr
.flatten()
.into_iter()
.reduce(|a, b| if a < b { a } else { b })
.unwrap_or_else(|| T::max_value()),
Operation::Sum => arr.flatten().into_iter().sum(),
Operation::Mean => {
arr.flatten().into_iter().sum::<T>()
/ T::from_usize(arr.len()).ok_or_else(|| {
anyhow!("cannot convert {} into {}", arr.len(), type_name::<T>())
})?
}
})
}
/// maximum intensity
pub fn max<T>(&self) -> Result<T>
where
T: Number + Sum,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
self.get_stat(Operation::Max)
}
/// minimum intensity
pub fn min<T>(&self) -> Result<T>
where
T: Number + Sum,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
self.get_stat(Operation::Min)
}
/// sum intensity
pub fn sum<T>(&self) -> Result<T>
where
T: Number + Sum,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
self.get_stat(Operation::Sum)
}
/// mean intensity
pub fn mean<T>(&self) -> Result<T>
where
T: Number + Sum,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
self.get_stat(Operation::Mean)
}
}
impl<D: Dimension> Deref for View<D> {
type Target = Reader;
fn deref(&self) -> &Self::Target {
self.reader.as_ref()
}
}
impl<T, D> TryFrom<View<D>> for Array<T, D>
where
T: Number,
D: Dimension,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
type Error = Error;
fn try_from(view: View<D>) -> Result<Self, Self::Error> {
view.as_array()
}
}
impl<T, D> TryFrom<&View<D>> for Array<T, D>
where
T: Number,
D: Dimension,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
type Error = Error;
fn try_from(view: &View<D>) -> Result<Self, Self::Error> {
view.as_array()
}
}
/// trait to define a function to retrieve the only item in a 0d array
pub trait Item {
fn item<T>(&self) -> Result<T>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax;
}
impl Item for View<Ix0> {
fn item<T>(&self) -> Result<T>
where
T: Number,
ArrayD<T>: MinMax,
Array1<T>: MinMax,
Array2<T>: MinMax,
{
Ok(self
.as_array()?
.first()
.ok_or_else(|| anyhow!("Empty view"))?
.clone())
}
}

24
tests/test_open.py
View File

@@ -6,16 +6,20 @@ import pytest
from ndbioimage import Imread
@pytest.mark.parametrize('file',
[file for file in (Path(__file__).parent / 'files').iterdir() if not file.suffix == '.pzl'])
@pytest.mark.parametrize(
"file",
[
file
for file in (Path(__file__).parent / "files").iterdir()
if not file.suffix == ".pzl"
],
)
def test_open(file):
with Imread(file) as im:
mean = im[dict(c=0, z=0, t=0)].mean()
with Imread(file, axes="cztyx") as im:
mean = im[0, 0, 0].mean()
b = pickle.dumps(im)
jm = pickle.loads(b)
assert jm[dict(c=0, z=0, t=0)].mean() == mean
v = im.view()
assert v[dict(c=0, z=0, t=0)].mean() == mean
b = pickle.dumps(v)
w = pickle.loads(b)
assert w[dict(c=0, z=0, t=0)].mean() == mean
assert jm.get_frame(0, 0, 0).mean() == mean
b = pickle.dumps(im)
jm = pickle.loads(b)
assert jm[0, 0, 0].mean() == mean

15
tests/test_slicing.py
View File

@@ -12,19 +12,24 @@ from ndbioimage import Imread
@pytest.fixture
def array():
return np.random.randint(0, 255, (64, 64, 2, 3, 4), 'uint8')
return np.random.randint(0, 255, (64, 64, 2, 3, 4), "uint8")
@pytest.fixture()
def image(array):
with tempfile.TemporaryDirectory() as folder:
file = Path(folder) / "test.tif"
tiffwrite(file, array, 'yxczt')
with Imread(file, axes='yxczt') as im:
tiffwrite(file, array, "yxczt")
with Imread(file, axes="yxczt") as im:
yield im
@pytest.mark.parametrize('s', combinations_with_replacement(
(0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1)), 5))
@pytest.mark.parametrize(
"s",
combinations_with_replacement(
(0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1)), 5
),
)
def test_slicing(s, image, array):
s_im, s_a = image[s], array[s]
if isinstance(s_a, Number):

37
tests/test_ufuncs.py
View File

@@ -11,21 +11,42 @@ from ndbioimage import Imread
@pytest.fixture
def array():
return np.random.randint(0, 255, (64, 64, 2, 3, 4), 'uint16')
return np.random.randint(0, 255, (64, 64, 2, 3, 4), "uint16")
@pytest.fixture()
def image(array):
with tempfile.TemporaryDirectory() as folder:
file = Path(folder) / "test.tif"
tiffwrite(file, array, 'yxczt')
with Imread(file, axes='yxczt') as im:
tiffwrite(file, array, "yxczt")
with Imread(file, axes="yxczt") as im:
yield im
@pytest.mark.parametrize('fun_and_axis', product(
(np.sum, np.nansum, np.min, np.nanmin, np.max, np.nanmax, np.argmin, np.argmax,
np.mean, np.nanmean, np.var, np.nanvar, np.std, np.nanstd), (None, 0, 1, 2, 3, 4)))
@pytest.mark.parametrize(
"fun_and_axis",
product(
(
np.sum,
np.nansum,
np.min,
np.nanmin,
np.max,
np.nanmax,
np.argmin,
np.argmax,
np.mean,
np.nanmean,
np.var,
np.nanvar,
np.std,
np.nanstd,
),
(None, 0, 1, 2, 3, 4),
),
)
def test_ufuncs(fun_and_axis, image, array):
fun, axis = fun_and_axis
assert np.all(np.isclose(fun(image, axis), fun(array, axis))), \
f'function {fun.__name__} over axis {axis} does not give the correct result'
assert np.all(np.isclose(np.asarray(fun(image, axis)), fun(array, axis))), (
f"function {fun.__name__} over axis {axis} does not give the correct result"
)