use crate::axes::{Ax, Axis, Operation, Shape, Slice, SliceInfoElemDef, slice_info}; use crate::cache::{ArrayCache, ArrayKeyRef, FrameCache, FrameKeyRef, ReaderKeyRef}; use crate::error::Error; use crate::metadata::Metadata; use crate::readers::{Dimensions, DynReader, Frame, Reader}; use crate::stats::MinMax; #[cfg(feature = "transforms")] use crate::transforms::Transforms; use indexmap::IndexMap; use itertools::{Itertools, iproduct}; use ndarray::{ Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn, SliceArg, SliceInfoElem, s, }; use num::traits::ToBytes; use num::{Bounded, FromPrimitive, ToPrimitive, Zero}; use serde::{Deserialize, Serialize}; use serde_with::serde_as; use std::any::type_name; use std::collections::{HashMap, HashSet}; use std::fmt::{Debug, Display, Formatter}; use std::hash::{Hash, Hasher}; use std::iter::Sum; use std::marker::PhantomData; use std::ops::{AddAssign, Deref, Div}; use std::path::Path; use std::sync::Arc; fn idx_bnd(idx: isize, bnd: isize) -> Result { if idx < -bnd { Err(Error::OutOfBounds(idx, bnd)) } else if idx < 0 { Ok(bnd - idx) } else if idx < bnd { Ok(idx) } else { Err(Error::OutOfBounds(idx, bnd)) } } fn slc_bnd(idx: isize, bnd: isize) -> Result { if idx < -bnd { Err(Error::OutOfBounds(idx, bnd)) } else if idx < 0 { Ok(bnd - idx) } else if idx <= bnd { Ok(idx) } else { Err(Error::OutOfBounds(idx, bnd)) } } /// a trait for numeric types that can be used as pixel values pub trait Number: 'static + Send + Sync + AddAssign + Bounded + Clone + Div + FromPrimitive + PartialOrd + Zero { } impl Number for T where T: 'static + Send + Sync + AddAssign + Bounded + Clone + Div + FromPrimitive + PartialOrd + Zero { } /// sliceable view on an image file #[serde_as] #[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] pub struct View { reader: R, /// same order as axes #[serde_as(as = "Vec")] slice: Vec, /// always has all of cztyx with possibly some new axes added axes: Vec, operations: IndexMap, dimensionality: PhantomData, #[cfg(feature = "transforms")] pub(crate) transforms: Transforms, } impl Hash for View where D: Dimension, R: Reader, { fn hash(&self, state: &mut H) { self.reader.hash(state); self.slice.hash(state); self.axes.hash(state); for (ax, op) in self.operations.iter() { ax.hash(state); op.hash(state); } } } impl View { pub(crate) fn new(reader: R, slice: Vec, axes: Vec) -> Self { Self { reader, slice, axes, operations: IndexMap::new(), dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: Transforms::default(), } } #[expect(dead_code)] pub(crate) fn new_with_axes(reader: R, axes: Vec) -> Result { let mut slice = Vec::new(); let shape = reader.shape(); for axis in axes.iter() { match axis { Axis::C => slice.push(SliceInfoElem::Slice { start: 0, end: Some(shape.c as isize), step: 1, }), Axis::Z => slice.push(SliceInfoElem::Slice { start: 0, end: Some(shape.z as isize), step: 1, }), Axis::T => slice.push(SliceInfoElem::Slice { start: 0, end: Some(shape.t as isize), step: 1, }), Axis::Y => slice.push(SliceInfoElem::Slice { start: 0, end: Some(shape.y as isize), step: 1, }), Axis::X => slice.push(SliceInfoElem::Slice { start: 0, end: Some(shape.x as isize), step: 1, }), Axis::New => { slice.push(SliceInfoElem::NewAxis); } } } let mut axes = axes.clone(); for axis in [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] { if !axes.contains(&axis) { let size = match axis { Axis::C => shape.c, Axis::Z => shape.z, Axis::T => shape.t, Axis::Y => shape.y, Axis::X => shape.x, Axis::New => 1, }; if size > 1 { return Err(Error::OutOfBoundsAxis(format!("{:?}", axis), size)); } slice.push(SliceInfoElem::Index(0)); axes.push(axis); } } Ok(Self { reader, slice, axes, operations: IndexMap::new(), dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: Transforms::default(), }) } /// the file path pub fn path(&self) -> &Path { self.reader.path() } /// the series in the file pub fn series(&self) -> usize { self.reader.series() } /// set the maximum number of frames the process-wide cache may hold, /// evicting the least recently used frames pub fn with_cache_capacity(self, capacity: usize) -> Self { FrameCache::global().set_capacity(capacity); self } /// the maximum number of frames the process-wide cache may hold pub fn cache_capacity(&self) -> usize { FrameCache::global().capacity() } /// set the maximum number of materialized arrays the process-wide cache /// may hold, evicting the least recently used arrays pub fn with_array_cache_capacity(self, capacity: usize) -> Self { ArrayCache::global().set_capacity(capacity); self } /// the maximum number of materialized arrays the process-wide cache may hold pub fn array_cache_capacity(&self) -> usize { ArrayCache::global().capacity() } fn with_operations(mut self, operations: IndexMap) -> Self { self.operations = operations; self } /// change the dimension into a dynamic dimension pub fn into_dyn(self) -> View { View { reader: self.reader, slice: self.slice, axes: self.axes, operations: self.operations, dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: Transforms::default(), } } /// change the dimension into a concrete dimension pub fn into_dimensionality(self) -> Result, Error> { 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, #[cfg(feature = "transforms")] transforms: self.transforms, }) } else { Err(Error::DimensionalityMismatch(d, self.ndim())) } } else { Ok(View { reader: self.reader, slice: self.slice, axes: self.axes, operations: self.operations, dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: self.transforms, }) } } /// the order of the axes, including axes sliced out pub fn get_axes(&self) -> &[Axis] { &self.axes } #[expect(dead_code)] pub(crate) fn get_operations(&self) -> &IndexMap { &self.operations } /// the slice defining the view pub fn get_slice(&self) -> &[SliceInfoElem] { &self.slice } /// the axes in the view pub fn axes(&self) -> Vec { 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() } /// remove axes of size 1 pub fn squeeze(&self) -> Result, Error> { let view = self.clone().into_dyn(); let slice: Vec<_> = self .shape() .into_iter() .map(|s| { if s == 1 { SliceInfoElem::Index(0) } else { SliceInfoElem::Slice { start: 0, end: None, step: 1, } } }) .collect(); view.slice(slice.as_slice()) } pub(crate) fn op_axes(&self) -> Vec { 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 first dimension pub fn len(&self) -> usize { self.shape()[0] } /// whether the view has no pixels in the first dimension pub fn is_empty(&self) -> bool { self.shape()[0] == 0 } /// the number of pixels in the view pub fn size(&self) -> usize { self.shape().into_iter().product() } /// the size of the view along an axis, or `None` if the axis is not present pub fn size_ax(&self, ax: Axis) -> Option { self.axes() .iter() .position(|a| *a == ax) .map(|i| self.shape()[i]) } /// the shape of the view pub fn shape(&self) -> Shape { let mut shape = Shape::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.order.push(*ax); shape.set_axis(ax, ((e - start).max(0) / step) as usize); } else { panic!("slice has no end") } } } SliceInfoElem::Index(_) => {} SliceInfoElem::NewAxis => { if !self.operations.contains_key(ax) { shape.order.push(*ax); } } } } shape } /// swap two axes pub fn swap_axes(&self, axis0: A, axis1: A) -> Result { 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).with_operations(self.operations.clone())) } /// subset of gives axes will be reordered in given order pub fn permute_axes(&self, axes: &[A]) -> Result { let idx: Vec = 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).with_operations(self.operations.clone())) } /// subset of gives axes will be reordered in given order, only axes in axes will be returned pub fn permute_axes_dyn(&self, axes: &[A]) -> Result, Error> { let idx: Vec = 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 new_slice = self.slice.to_vec(); let mut new_axes = self.axes.clone(); let axes = axes.iter().map(|a| a.n()).collect::>(); for (&i, j) in idx.iter().zip(jdx) { new_slice[j] = self.slice[i]; new_axes[j] = self.axes[i]; } for (ax, s) in new_axes.iter().zip(new_slice.iter_mut()) { if let SliceInfoElem::Slice { start, end, step } = s && !axes.contains(&ax.n()) { let size = ((end.expect("slice has no end") - *start).max(0) / *step) as usize; if size != 1 { return Err(Error::SizeMismatch(ax.to_string(), size)); } *s = SliceInfoElem::Index(*start) } } Ok(View::new(self.reader.clone(), new_slice, new_axes) .with_operations(self.operations.clone())) } /// reverse the order of the axes pub fn transpose(&self) -> Result { Ok(View::new( self.reader.clone(), self.slice.iter().rev().cloned().collect(), self.axes.iter().rev().cloned().collect(), ) .with_operations(self.operations.clone())) } /// reduce an axis with an operation (max, min, sum or mean), returning a view of lower dimensionality pub fn operate( &self, axis: A, operation: Operation, ) -> Result, Error> { let pos = axis.pos_op(&self.axes, &self.slice, &self.op_axes())?; let ax = self.axes[pos]; let (axes, slice, operations) = if Axis::New == ax { let mut axes = self.axes.clone(); let mut slice = self.slice.clone(); axes.remove(pos); slice.remove(pos); (axes, slice, self.operations.clone()) } else if self.operations.contains_key(&ax) { if D::NDIM.is_none() { ( self.axes.clone(), self.slice.clone(), self.operations.clone(), ) } else { return Err(Error::AxisAlreadyOperated(pos, ax.to_string())); } } else { let mut operations = self.operations.clone(); operations.insert(ax, operation); (self.axes.clone(), self.slice.clone(), operations) }; Ok(View::new(self.reader.clone(), slice, axes).with_operations(operations)) } /// maximum along axis pub fn max_proj(&self, axis: A) -> Result, Error> { self.operate(axis, Operation::Max) } /// minimum along axis pub fn min_proj(&self, axis: A) -> Result, Error> { self.operate(axis, Operation::Min) } /// sum along axis pub fn sum_proj(&self, axis: A) -> Result, Error> { self.operate(axis, Operation::Sum) } /// mean along axis pub fn mean_proj(&self, axis: A) -> Result, Error> { self.operate(axis, Operation::Mean) } /// created a new sliced view pub fn slice(&self, info: I) -> Result, Error> where I: SliceArg, { if self.slice.out_ndim() < info.in_ndim() { return Err(Error::NotEnoughFreeDimensions); } 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.get(n_idx); let r = reader_slice.get(r_idx); let a = self.axes.get(r_idx); match a { Some(i) if self.operations.contains_key(i) => { new_slice.push(*r.expect("slice should exist for axes under operation")); new_axes.push(*i); r_idx += 1; } _ => match (n, r) { ( Some(SliceInfoElem::Slice { start: info_start, end: info_end, step: info_step, }), Some(SliceInfoElem::Slice { start, end, step }), ) => { let new_start = start + info_start; let end = end.expect("slice has no end"); let new_end = if let Some(m) = info_end { end.min(start + info_step * m) } else { end }; let new_step = (step * info_step).abs(); if new_start > end { return Err(Error::OutOfBounds(*info_start, (end - start) / step)); } new_slice.push(SliceInfoElem::Slice { start: new_start, end: Some(new_end), step: new_step, }); new_axes.push(*a.expect("axis should exist when slice exists")); n_idx += 1; r_idx += 1; } ( Some(SliceInfoElem::Index(k)), Some(SliceInfoElem::Slice { start, end, step }), ) => { let i = if *k < 0 { end.unwrap_or(0) + step.abs() * k } else { start + step.abs() * k }; let end = end.expect("slice has no end"); if i >= end { return Err(Error::OutOfBounds(i, (end - start) / step)); } new_slice.push(SliceInfoElem::Index(i)); new_axes.push(*a.expect("axis should exist when slice exists")); n_idx += 1; r_idx += 1; } (Some(SliceInfoElem::Slice { start, .. }), Some(SliceInfoElem::NewAxis)) => { if *start != 0 { return Err(Error::OutOfBounds(*start, 1)); } new_slice.push(SliceInfoElem::NewAxis); new_axes.push(Axis::New); n_idx += 1; r_idx += 1; } (Some(SliceInfoElem::Index(k)), Some(SliceInfoElem::NewAxis)) => { if *k != 0 { return Err(Error::OutOfBounds(*k, 1)); } n_idx += 1; r_idx += 1; } (Some(SliceInfoElem::NewAxis), Some(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; } (Some(SliceInfoElem::NewAxis), _) => { new_slice.push(SliceInfoElem::NewAxis); new_axes.push(Axis::New); n_idx += 1; } (_, Some(SliceInfoElem::Index(k))) => { new_slice.push(SliceInfoElem::Index(*k)); new_axes.push(*a.expect("axis should exist when slice exists")); r_idx += 1; } _ => unreachable!(), }, } } 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())) } /// resets axes to cztyx order, with all 5 axes present, /// inserts new axes in place of axes under operation (max_proj etc.) pub fn reset_axes(&self) -> Result, Error> { let mut axes = Vec::new(); let mut slice = Vec::new(); for ax in [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] { axes.push(ax); let s = self.slice[ax.pos(&self.axes, &self.slice)?]; match s { SliceInfoElem::Slice { .. } => slice.push(s), SliceInfoElem::Index(i) => slice.push(SliceInfoElem::Slice { start: i, end: Some(i + 1), step: 1, }), SliceInfoElem::NewAxis => { panic!("slice should not be NewAxis when axis is one of cztyx") } } if self.operations.contains_key(&ax) { axes.push(Axis::New); slice.push(SliceInfoElem::NewAxis) } } Ok(View::new(self.reader.clone(), slice, 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(&self, info: I) -> Result, Error> where I: SliceArg, { self.reset_axes()?.slice(info)?.into_dimensionality() } /// the pixel intensity at a given index pub fn item_at(&self, index: &[isize]) -> Result where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: 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(&self) -> Result, Error> where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { Ok(self.as_array_dyn()?.into_dimensionality()?) } /// collect the view into a dynamic-dimension ndarray pub fn as_array_dyn(&self) -> Result, Error> where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { let key = self.array_key::(); if let Some(arr) = ArrayCache::global().get(&key) { return Ok(arr.as_ref().clone()); } let mut op_xy = IndexMap::new(); if let Some((&ax, op)) = self.operations.first() && ((ax == Axis::X) || (ax == Axis::Y)) { op_xy.insert(ax, op.clone()); if let Some((&ax2, op2)) = self.operations.get_index(1) && ((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 } => { let end = end.expect("slice has no end"); if !op_xy.contains_key(a) && !op_czt.contains_key(a) { shape_out.push(((end - start).max(0) / step) as usize); slice.push(SliceInfoElem::Slice { start: 0, end: None, step: 1, }); ax_out.push(*a); } } 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_reader = self.reader.shape(); 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 && !op_xy.contains_key(&axis) { xy_dim += 1; } slice_reader[axis as usize] = Slice::new( idx_bnd(*start, shape_reader[axis] as isize)?, slc_bnd(end.unwrap(), shape_reader[axis] as isize)?, *step, ); } SliceInfoElem::Index(j) => { slice_reader[axis as usize] = Slice::new( idx_bnd(*j, shape_reader[axis] as isize)?, slc_bnd(*j + 1, shape_reader[axis] as isize)?, 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::::from_elem(shape_out.into_dimension(), T::min_value()) } Operation::Min => { ArrayD::::from_elem(shape_out.into_dimension(), T::max_value()) } _ => ArrayD::::zeros(shape_out.into_dimension()), } } else { ArrayD::::zeros(shape_out.into_dimension()) }; let shape = self.reader.shape(); let mut axes_out_idx = [None; 5]; for (i, ax) in ax_out.iter().enumerate() { if *ax < Axis::New { 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.get_cached_frame( (c % shape.c as isize) as usize, (z % shape.z as isize) as usize, (t % shape.t as isize) as usize, )?; let arr_frame: Array2 = frame.as_ref().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::(&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 = op1.operate(b.to_owned(), ax1 as usize - 3)?; c.to_owned().into_dyn() } else if op_xy.contains_key(&Axis::X) || op_xy.contains_key(&Axis::Y) { let xys = slice_info::(&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)?; b.to_owned().into_dyn() } else { let xys = slice_info::(&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::(&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)?; b.to_owned().into_dyn() } else { let xys = slice_info::(&xy)?; arr_frame.slice(xys).to_owned().into_dyn() } } 2 => { let xys = slice_info::(&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() } } _ => { unreachable!("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 mut ax_out: HashMap = 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.remove(ax) { for i in ax_out.values_mut() { if *i > idx { *i -= 1; } } let arr = out.take().unwrap(); let a = op.operate(arr, idx)?; let _ = out.insert(a); } } let n = if let Some((&ax, op)) = op_czt.first() && *op == Operation::Mean { self.axes .iter() .zip(self.slice.iter()) .find(|(a, _)| **a == ax) .and_then(|(_, s)| match s { SliceInfoElem::Slice { start, end, step } => { end.map(|e| (((e - start).max(0) / step) as usize).max(1)) } _ => Some(1), }) .unwrap_or(1) } else { 1 }; 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()) }; ArrayCache::global().insert(key.to_owned(), array.clone()); Ok(array) } /// turn the view into a 1d array pub fn flatten(&self) -> Result, Error> where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { Ok(Array1::from_iter(self.as_array()?.iter().cloned())) } /// turn the data into a byte vector pub fn to_bytes(&self) -> Result, Error> where T: Number + ToBytesVec, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { Ok(self .as_array()? .iter() .flat_map(|i| i.to_bytes_vec()) .collect()) } fn frame_key(&self, c: usize, z: usize, t: usize) -> FrameKeyRef<'_> { FrameKeyRef { reader: ReaderKeyRef { name: self.reader.reader_name(), path: self.reader.path(), series: self.reader.series(), position: self.reader.position(), }, c, z, t, } } fn array_key(&self) -> ArrayKeyRef<'_> { ArrayKeyRef { reader: ReaderKeyRef { name: self.reader.reader_name(), path: self.reader.path(), series: self.reader.series(), position: self.reader.position(), }, dtype: type_name::(), slice: &self.slice, axes: &self.axes, operations: &self.operations, } } fn get_cached_frame(&self, c: usize, z: usize, t: usize) -> Result, Error> { let key = self.frame_key(c, z, t); if let Some(frame) = FrameCache::global().get(&key) { return Ok(frame); } let frame = Arc::new(self.reader.get_frame(c, z, t)?); FrameCache::global().insert(key.to_owned(), frame.clone()); Ok(frame) } /// retrieve a single frame at czt, sliced accordingly pub fn get_frame(&self, c: N, z: N, t: N) -> Result, Error> where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, N: Display + ToPrimitive, { let c = c .to_isize() .ok_or_else(|| Error::Cast(c.to_string(), "isize".to_string()))?; let z = z .to_isize() .ok_or_else(|| Error::Cast(z.to_string(), "isize".to_string()))?; let t = t .to_isize() .ok_or_else(|| Error::Cast(t.to_string(), "isize".to_string()))?; self.slice_cztyx(s![c, z, t, .., ..])?.as_array() } fn get_stat(&self, operation: Operation) -> Result where T: Number + Sum, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { let arr: ArrayD = 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::from_usize(arr.len()).ok_or_else(|| { Error::Cast(arr.len().to_string(), type_name::().to_string()) })? } }) } /// maximum intensity pub fn max(&self) -> Result where T: Number + Sum, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { self.get_stat(Operation::Max) } /// minimum intensity pub fn min(&self) -> Result where T: Number + Sum, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { self.get_stat(Operation::Min) } /// sum intensity pub fn sum(&self) -> Result where T: Number + Sum, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { self.get_stat(Operation::Sum) } /// mean intensity pub fn mean(&self) -> Result where T: Number + Sum, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { self.get_stat(Operation::Mean) } /// gives a helpful summary of the recorded experiment pub fn summary(&self) -> Result { let mut s = "".to_string(); s.push_str(&format!("path/filename: {}\n", self.path().display())); s.push_str(&format!("series/pos: {}\n", self.series())); s.push_str(&format!("reader: {}\n", self.reader_name())); s.push_str(&format!("dtype: {:?}\n", self.pixel_type())); let axes = self .axes() .into_iter() .map(|ax| format!("{}", ax)) .join("") .to_lowercase(); let shape = self .shape() .into_iter() .map(|s| format!("{}", s)) .join(" x "); let space = " ".repeat(6usize.saturating_sub(axes.len())); s.push_str(&format!("shape ({}):{}{}\n", axes, space, shape)); s.push_str(&self.metadata()?.summary()?); Ok(s) } } impl Deref for View { type Target = R; fn deref(&self) -> &Self::Target { &self.reader } } impl TryFrom> for Array where T: Number, D: Dimension, R: Reader, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { type Error = Error; fn try_from(view: View) -> Result { view.as_array() } } impl TryFrom<&View> for Array where T: Number, D: Dimension, R: Reader, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { type Error = Error; fn try_from(view: &View) -> Result { view.as_array() } } /// trait to define a function to retrieve the only item in a 0d array pub trait Item { /// the single item in the 0d array, cast to `T` fn item(&self) -> Result where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>; } impl View { /// create a view on the image at `path`, parsing series and position from the file name pub fn from_path

(path: P) -> Result where P: AsRef, { let (path, dimensions) = Dimensions::parse_path(path)?; Ok(R::new(path, dimensions.s.unwrap_or(0), dimensions.p.unwrap_or(0))?.view()) } } impl Item for View { fn item(&self) -> Result where T: Number, ArrayD: MinMax>, Array1: MinMax>, Array2: MinMax>, { Ok(self.as_array()?.first().ok_or(Error::EmptyView)?.clone()) } } impl Display for View { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { if let Ok(summary) = self.summary() { write!(f, "{}", summary) } else { write!(f, "{}", self.path().display()) } } } /// trait to convert numbers to bytes pub trait ToBytesVec { /// the number as a vector of bytes in native endianness fn to_bytes_vec(&self) -> Vec; } macro_rules! to_bytes_vec_impl { ($($t:ty $(,)?)*) => { $( impl ToBytesVec for $t { #[inline] fn to_bytes_vec(&self) -> Vec { self.to_ne_bytes().to_vec() } } )* }; } to_bytes_vec_impl!( u8, u16, u32, u64, u128, usize, i8, i16, i32, i64, i128, isize, f32, f64 ); #[cfg(test)] mod tests { use crate::axes::{Axis, Operation}; use crate::cache::{ArrayCache, ArrayKey, ArrayKeyRef, FrameCache, ReaderKey, ReaderKeyRef}; use crate::error::Error; use crate::readers::{DynReader, Frame, Reader}; use crate::stats::MinMax; use crate::view::Item; use indexmap::IndexMap; use ndarray::{Array, Array4, Array5, NewAxis}; use ndarray::{Array2, ArrayD, IxDyn, SliceInfoElem, s}; use std::path::{Path, PathBuf}; use std::sync::Arc; fn open(file: &str) -> Result { let path = std::env::current_dir()? .join("tests") .join("files") .join(file); DynReader::new(&path, 0, 0) } #[test] fn view() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); println!("view: {:?}", view); let a = view.slice(s![0, 5, 0, .., ..])?; println!("a: {:?}", a); let b = reader.get_frame(0, 5, 0)?; let c: Array2 = a.try_into()?; println!("c {:?}", c); let d: Array2 = b.try_into()?; println!("d {:?}", d); assert_eq!(c, d); Ok(()) } #[test] fn view_shape() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; 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.to_vec(), vec![5, 1024, 100]); Ok(()) } #[test] fn view_new_axis() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); let a = Array5::::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().to_vec(), a.shape()); let a = a.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis]); let v = v.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis])?; assert_eq!(v.shape().to_vec(), a.shape()); Ok(()) } #[test] fn view_permute_axes() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); let s = view.shape(); let mut a = Array5::::zeros((s[0], s[1], s[2], s[3], s[4])); assert_eq!(view.shape().to_vec(), a.shape()); let b: Array5 = 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().to_vec(), a.shape()); let b: Array5 = view.clone().try_into()?; assert_eq!(b.shape(), a.shape()); let view = view.permute_axes(&[Axis::X, Axis::Z, Axis::Y])?; let a = a.permuted_axes([4, 1, 2, 0, 3]); assert_eq!(view.shape().to_vec(), a.shape()); let b: Array5 = view.clone().try_into()?; assert_eq!(b.shape(), a.shape()); Ok(()) } macro_rules! test_max { ($($name:ident: $b:expr $(,)?)*) => { $( #[test] fn $name() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); let array: Array5 = view.clone().try_into()?; let view = view.max_proj($b)?; let a: Array4 = 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<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); let v4: Array = view.slice($b)?.try_into()?; let a5: Array5 = 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<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; 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::()?; assert_eq!(d.shape(), [1024, 1024]); Ok(()) } #[test] fn item() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; let reader = open(file)?; let view = reader.view(); let a = view.slice(s![.., 0, 0, 0, 0])?; let b = a.slice(s![0])?; let item = b.item::()?; assert_eq!(item, 2); Ok(()) } #[test] fn slice_cztyx() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; let reader = open(file)?; let view = reader.view().max_proj(Axis::Z)?.into_dyn(); println!("view.axes: {:?}", view.get_axes()); println!("view.slice: {:?}", view.get_slice()); let r = view.reset_axes()?; println!("r.axes: {:?}", r.get_axes()); println!("r.slice: {:?}", r.get_slice()); let a = view.slice_cztyx(s![0, 0, 0, .., ..])?; println!("a.axes: {:?}", a.get_axes()); println!("a.slice: {:?}", a.get_slice()); assert_eq!(a.axes(), [Axis::Y, Axis::X]); Ok(()) } #[test] fn reset_axes() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; let reader = open(file)?; let view = reader.view().max_proj(Axis::Z)?; let view = view.reset_axes()?; assert_eq!(view.axes(), [Axis::C, Axis::New, Axis::T, Axis::Y, Axis::X]); let a = view.as_array::()?; assert_eq!(a.ndim(), 5); Ok(()) } #[test] fn reset_axes2() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; let reader = open(file)?; let view = reader.view().squeeze()?; let a = view.reset_axes()?; assert_eq!(a.axes(), [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X]); Ok(()) } #[test] fn reset_axes3() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; let reader = open(file)?; let view4 = reader.view().squeeze()?; let view = view4.max_proj(Axis::Z)?.into_dyn(); let slice = view.slice_cztyx(s![0, .., .., .., ..])?.into_dyn(); let a = slice.as_array::()?; assert_eq!(slice.shape().to_vec(), [1, 10, 1280, 1280]); assert_eq!(a.shape(), [1, 10, 1280, 1280]); let r = slice.reset_axes()?; let b = r.as_array::()?; assert_eq!(r.shape().to_vec(), [1, 1, 10, 1280, 1280]); assert_eq!(b.shape(), [1, 1, 10, 1280, 1280]); let q = slice.max_proj(Axis::C)?.max_proj(Axis::T)?; let c = q.as_array::()?; assert_eq!(q.shape().to_vec(), [1, 1280, 1280]); assert_eq!(c.shape().to_vec(), [1, 1280, 1280]); let p = q.reset_axes()?; let d = p.as_array::()?; println!("axes: {:?}", p.get_axes()); println!("operations: {:?}", p.get_operations()); println!("slice: {:?}", p.get_slice()); assert_eq!(p.shape().to_vec(), [1, 1, 1, 1280, 1280]); assert_eq!(d.shape(), [1, 1, 1, 1280, 1280]); Ok(()) } #[test] fn max() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; let reader = open(file)?; let view = reader.view(); let m = view.max_proj(Axis::T)?; let a = m.as_array::()?; assert_eq!(m.shape().to_vec(), [2, 1, 1280, 1280]); assert_eq!(a.shape(), [2, 1, 1280, 1280]); let mc = view.max_proj(Axis::C)?; let a = mc.as_array::()?; assert_eq!(mc.shape().to_vec(), [1, 10, 1280, 1280]); assert_eq!(a.shape(), [1, 10, 1280, 1280]); let mz = mc.max_proj(Axis::Z)?; let a = mz.as_array::()?; assert_eq!(mz.shape().to_vec(), [10, 1280, 1280]); assert_eq!(a.shape(), [10, 1280, 1280]); let mt = mz.max_proj(Axis::T)?; let a = mt.as_array::()?; assert_eq!(mt.shape().to_vec(), [1280, 1280]); assert_eq!(a.shape(), [1280, 1280]); Ok(()) } #[test] fn frame_cache() -> Result<(), Error> { let cache = FrameCache::new(2); let key = |c: usize, z: usize, t: usize| { ( ReaderKey { name: "test".to_string(), path: PathBuf::from("test.tif"), series: 0, position: 0, }, c, z, t, ) }; cache.insert( key(0, 0, 0), Arc::new(Frame::from(Array2::from_elem((2, 2), 1u16))), ); cache.insert( key(0, 0, 1), Arc::new(Frame::from(Array2::from_elem((2, 2), 2u16))), ); assert!(cache.get(&key(0, 0, 0)).is_some()); cache.insert( key(0, 0, 2), Arc::new(Frame::from(Array2::from_elem((2, 2), 3u16))), ); assert!(cache.get(&key(0, 0, 0)).is_some()); assert!(cache.get(&key(0, 0, 1)).is_none()); assert!(cache.get(&key(0, 0, 2)).is_some()); assert_eq!(cache.len(), 2); Ok(()) } #[test] fn array_cache() -> Result<(), Error> { let cache = ArrayCache::new(2); let key = |dtype: &'static str, index: isize, axes: &[Axis]| ArrayKey { reader: ReaderKey { name: "test".to_string(), path: PathBuf::from("test.tif"), series: 0, position: 0, }, dtype, slice: vec![SliceInfoElem::Index(index)], axes: axes.to_vec(), operations: vec![], }; let k0 = key("u16", 0, &[Axis::T]); let k1 = key("u16", 1, &[Axis::T]); cache.insert(k0.clone(), ArrayD::::zeros(IxDyn(&[2, 2]))); cache.insert(k1.clone(), ArrayD::::zeros(IxDyn(&[2, 2]))); assert!(cache.get::(&k0).is_some()); cache.insert( key("u16", 2, &[Axis::T]), ArrayD::::zeros(IxDyn(&[2, 2])), ); assert!(cache.get::(&k0).is_some()); assert!(cache.get::(&k1).is_none()); assert_eq!(cache.len(), 2); let ops: IndexMap = IndexMap::new(); let borrowed = ArrayKeyRef { reader: ReaderKeyRef { name: "test", path: Path::new("test.tif"), series: 0, position: 0, }, dtype: "u16", slice: &[SliceInfoElem::Index(0)], axes: &[Axis::T], operations: &ops, }; assert!(cache.get::(&borrowed).is_none()); assert!(cache.get::(&borrowed).is_some()); Ok(()) } #[test] fn as_array_caches_frames() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; let reader = open(file)?; let view = reader.view(); let shape = view.shape(); let (c, z, t) = (shape[0] - 1, shape[1] - 1, shape[2] - 1); let cache = FrameCache::global(); let a: Array5 = view.clone().try_into()?; assert!(cache.get(&view.frame_key(c, z, t)).is_some()); let b: Array5 = view.clone().try_into()?; assert_eq!(a, b); // a fresh reader on the same file shares the same cache key let view2 = open(file)?.view(); assert_eq!(view.frame_key(c, z, t), view2.frame_key(c, z, t)); let c2: Array5 = view2.clone().try_into()?; assert_eq!(a, c2); Ok(()) } }