845 lines
32 KiB
Rust
845 lines
32 KiB
Rust
use crate::axes::Axis;
|
|
use crate::error::Error;
|
|
use crate::readers::{DynReader, Frame, PixelType, Reader};
|
|
use crate::view::View;
|
|
use image_registration::{bspline::DefaultValue, filter, transform::Transform};
|
|
use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, array, concatenate, s};
|
|
use rayon::iter::{ParallelIterator, IntoParallelIterator};
|
|
use serde::{Deserialize, Serialize};
|
|
use std::borrow::Cow;
|
|
use std::collections::{HashMap, HashSet};
|
|
use std::fmt::Debug;
|
|
use std::hash::{Hash, Hasher};
|
|
use std::io::Seek;
|
|
use std::path::Path;
|
|
use std::sync::Arc;
|
|
use tiffwrite::IJTiffFile;
|
|
|
|
#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
|
|
struct TransformOld {
|
|
#[serde(alias = "CenterOfRotationPoint")]
|
|
center_of_rotation_point: Vec<f64>,
|
|
#[serde(alias = "Size")]
|
|
size: Vec<f64>,
|
|
#[serde(alias = "TransformParameters")]
|
|
transform_parameters: Vec<f64>,
|
|
#[serde(alias = "dTransformParameters")]
|
|
dtransform_parameters: Vec<f64>,
|
|
}
|
|
|
|
impl TryFrom<HashMap<String, TransformOld>> for Transforms {
|
|
type Error = Error;
|
|
|
|
fn try_from(t: HashMap<String, TransformOld>) -> Result<Self, Error> {
|
|
Ok(Self {
|
|
default_value: Default::default(),
|
|
interpolation_order: None,
|
|
channel: t
|
|
.into_iter()
|
|
.map(|(n, i)| {
|
|
let center = i
|
|
.size
|
|
.iter()
|
|
.zip(i.center_of_rotation_point.iter())
|
|
.map(|(s, c)| c - s / 2.0 + 0.5)
|
|
.collect::<Vec<_>>();
|
|
Ok((
|
|
n,
|
|
Transform::new_around_center(i.transform_parameters, ¢er)
|
|
.with_dparameters(i.dtransform_parameters)
|
|
.inverse()?,
|
|
))
|
|
})
|
|
.collect::<Result<_, Error>>()?,
|
|
drift: Vec::new(),
|
|
})
|
|
}
|
|
}
|
|
|
|
pub type DefaultValues = HashMap<String, DefaultValue>;
|
|
|
|
#[derive(Clone, Debug, Default, Eq, PartialEq, Deserialize)]
|
|
pub struct Transforms {
|
|
#[serde(default)]
|
|
pub(crate) interpolation_order: Option<usize>,
|
|
#[serde(default)]
|
|
pub(crate) default_value: DefaultValues,
|
|
#[serde(default)]
|
|
pub channel: HashMap<String, Transform<IxDyn>>,
|
|
#[serde(default)]
|
|
pub drift: Vec<Transform<IxDyn>>,
|
|
}
|
|
|
|
impl Hash for Transforms {
|
|
fn hash<H: Hasher>(&self, state: &mut H) {
|
|
let mut keys: Vec<_> = self.channel.keys().collect();
|
|
keys.sort();
|
|
for key in keys {
|
|
key.hash(state);
|
|
self.channel[key].hash(state);
|
|
}
|
|
self.drift.hash(state);
|
|
}
|
|
}
|
|
|
|
impl PartialEq<Transforms> for &Transforms {
|
|
fn eq(&self, other: &Transforms) -> bool {
|
|
**self == *other
|
|
}
|
|
}
|
|
|
|
impl Serialize for Transforms {
|
|
fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
|
|
where
|
|
S: serde::Serializer,
|
|
{
|
|
#[derive(Serialize)]
|
|
struct HumanReadable<'a> {
|
|
#[serde(skip_serializing_if = "HashMap::is_empty")]
|
|
channel: &'a HashMap<String, Transform<IxDyn>>,
|
|
#[serde(skip_serializing_if = "Vec::is_empty")]
|
|
drift: &'a Vec<Transform<IxDyn>>,
|
|
}
|
|
|
|
#[derive(Serialize)]
|
|
struct Binary<'a> {
|
|
interpolation_order: Option<usize>,
|
|
channel: &'a HashMap<String, Transform<IxDyn>>,
|
|
drift: &'a Vec<Transform<IxDyn>>,
|
|
}
|
|
|
|
if serializer.is_human_readable() {
|
|
HumanReadable {
|
|
channel: &self.channel,
|
|
drift: &self.drift,
|
|
}
|
|
.serialize(serializer)
|
|
} else {
|
|
Binary {
|
|
interpolation_order: self.interpolation_order,
|
|
channel: &self.channel,
|
|
drift: &self.drift,
|
|
}
|
|
.serialize(serializer)
|
|
}
|
|
}
|
|
}
|
|
|
|
#[derive(Clone, Debug, PartialEq, Copy)]
|
|
pub enum TransformKind {
|
|
None,
|
|
Channel,
|
|
Drift,
|
|
Both,
|
|
}
|
|
|
|
impl Transforms {
|
|
pub fn kind(&self) -> TransformKind {
|
|
match (self.channel.is_empty(), self.drift.is_empty()) {
|
|
(true, true) => TransformKind::None,
|
|
(false, true) => TransformKind::Channel,
|
|
(true, false) => TransformKind::Drift,
|
|
(false, false) => TransformKind::Both,
|
|
}
|
|
}
|
|
|
|
pub fn load<P: AsRef<Path>>(path: P) -> Result<Self, Error> {
|
|
let mut file = std::fs::File::open(&path)?;
|
|
Ok(
|
|
if let Ok(t) = yaml_serde::from_reader::<_, Self>(&file)
|
|
&& !(t.channel.is_empty() && t.drift.is_empty())
|
|
{
|
|
t
|
|
} else {
|
|
file.seek(std::io::SeekFrom::Start(0))?;
|
|
yaml_serde::from_reader::<_, HashMap<String, TransformOld>>(&file)?.try_into()?
|
|
},
|
|
)
|
|
}
|
|
|
|
pub fn save<P: AsRef<Path>>(&self, path: P) -> Result<(), Error> {
|
|
let file = std::fs::OpenOptions::new()
|
|
.write(true)
|
|
.create(true)
|
|
.truncate(true)
|
|
.open(path)?;
|
|
Ok(yaml_serde::to_writer(file, self)?)
|
|
}
|
|
|
|
pub fn calculate_transforms_from<P, Q, R>(
|
|
transform_file: Option<P>,
|
|
bead_files: Option<&[Q]>,
|
|
drift_file: Option<R>,
|
|
main_channel: Option<usize>,
|
|
default_transform: Option<Vec<f64>>,
|
|
) -> Result<Self, Error>
|
|
where
|
|
P: AsRef<Path>,
|
|
Q: AsRef<Path>,
|
|
R: AsRef<Path>,
|
|
{
|
|
if let Some(transform_file) = transform_file.as_ref()
|
|
&& let Ok(transforms) = Self::load(transform_file.as_ref().with_extension("yml"))
|
|
{
|
|
Ok(transforms)
|
|
} else {
|
|
let default_transform = default_transform.map(Transform::new);
|
|
let mut transforms = Self::default();
|
|
if let Some(bead_files) = bead_files {
|
|
let main_channel = main_channel.unwrap_or(0);
|
|
transforms.channel = Self::calculate_channel_transforms_2d(
|
|
bead_files,
|
|
main_channel,
|
|
default_transform,
|
|
)?
|
|
.into_iter()
|
|
.map(|(k, t)| (k, t.into_dyn()))
|
|
.collect();
|
|
};
|
|
if let Some(drift_file) = drift_file {
|
|
transforms.drift = Self::calculate_drift_transforms_2d(drift_file)?
|
|
.into_iter()
|
|
.map(|t| t.into_dyn())
|
|
.collect();
|
|
};
|
|
if let Some(transform_file) = transform_file {
|
|
let transform_file = transform_file.as_ref();
|
|
transforms.save(transform_file.with_extension("yml"))?;
|
|
if let Some(bead_files) = bead_files {
|
|
transforms.save_channel_transform_tiff(
|
|
transform_file.with_extension("tif"),
|
|
bead_files,
|
|
)?;
|
|
}
|
|
}
|
|
Ok(transforms)
|
|
}
|
|
}
|
|
|
|
/// save a tiff file with the uncorrected and corrected frames of the bead files side by side
|
|
pub fn save_channel_transform_tiff<P, Q>(&self, file: P, bead_files: &[Q]) -> Result<(), Error>
|
|
where
|
|
P: AsRef<Path>,
|
|
Q: AsRef<Path>,
|
|
{
|
|
let mut transforms = self.clone();
|
|
transforms.drift = Vec::new();
|
|
let mut tif = IJTiffFile::new(file)?;
|
|
for (t, bead_file) in bead_files.iter().enumerate() {
|
|
let mut view = View::<_, DynReader>::from_path(bead_file)?.into_dyn();
|
|
let axes = view.axes();
|
|
if axes.contains(&Axis::Z) {
|
|
view = view.max_proj(Axis::Z)?;
|
|
}
|
|
if axes.contains(&Axis::T) {
|
|
view = view.max_proj(Axis::T)?;
|
|
}
|
|
let view_u = view.into_dimensionality::<Ix3>()?;
|
|
let view_t = view_u.with_transforms(transforms.clone());
|
|
for c in 0..view_u.shape().c {
|
|
let a = view_u.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::<u16>()?;
|
|
let b = view_t.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::<u16>()?;
|
|
let frame = concatenate(ndarray::Axis(1), &[a.view(), b.view()])?;
|
|
tif.save(frame, c, 0, t)?;
|
|
}
|
|
}
|
|
Ok(())
|
|
}
|
|
|
|
pub fn calculate_drift_transforms_2d<P: AsRef<Path>>(
|
|
file: P,
|
|
) -> Result<Vec<Transform<Ix2>>, Error> {
|
|
View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_2d()
|
|
}
|
|
|
|
pub fn calculate_drift_transforms_3d<P: AsRef<Path>>(
|
|
file: P,
|
|
) -> Result<Vec<Transform<Ix3>>, Error> {
|
|
View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_3d()
|
|
}
|
|
|
|
pub fn calculate_channel_transforms_2d<P: AsRef<Path>>(
|
|
bead_files: &[P],
|
|
main_channel: usize,
|
|
default_transform: Option<Transform<Ix2>>,
|
|
) -> Result<HashMap<String, Transform<Ix2>>, Error> {
|
|
// TODO: autodetect main channel
|
|
let mut transforms = Vec::new();
|
|
let default_transform = default_transform.unwrap_or_default();
|
|
for file in bead_files {
|
|
let view = View::<_, DynReader>::from_path(file)?;
|
|
transforms.push(view.calculate_channel_transforms_2d(main_channel)?)
|
|
}
|
|
let channels = transforms
|
|
.iter()
|
|
.flat_map(|t| t.keys())
|
|
.collect::<HashSet<_>>();
|
|
let mut average_transforms = HashMap::new();
|
|
for channel in channels {
|
|
let (s, n) = transforms
|
|
.iter()
|
|
.filter_map(|t| {
|
|
if t.contains_key(channel) {
|
|
Some((&t[channel] * &default_transform).matrix())
|
|
} else {
|
|
None
|
|
}
|
|
})
|
|
.fold((Array2::zeros((3, 3)), 0), |(s, n), t| (s + t, n + 1));
|
|
let matrix = s / (n as f64);
|
|
let (s, n) = transforms
|
|
.iter()
|
|
.filter_map(|t| {
|
|
if t.contains_key(channel) {
|
|
Some(((&t[channel] * &default_transform).matrix() - &matrix).powi(2))
|
|
} else {
|
|
None
|
|
}
|
|
})
|
|
.fold((Array2::zeros((3, 3)), 0), |(s, n), t| (s + t, n + 1));
|
|
let dmatrix = (s / n as f64).sqrt();
|
|
average_transforms.insert(
|
|
channel.clone(),
|
|
Transform::default()
|
|
.with_matrix(matrix.view())
|
|
.with_dmatrix(dmatrix.view()),
|
|
);
|
|
}
|
|
Ok(average_transforms)
|
|
}
|
|
|
|
pub fn calculate_channel_transforms_3d<P: AsRef<Path>>(
|
|
bead_files: &[P],
|
|
main_channel: usize,
|
|
default_transform: Option<Transform<Ix3>>,
|
|
) -> Result<HashMap<String, Transform<Ix3>>, Error> {
|
|
let mut transforms = Vec::new();
|
|
let default_transform = default_transform.unwrap_or_default();
|
|
for file in bead_files {
|
|
let view = View::<_, DynReader>::from_path(file)?;
|
|
transforms.push(view.calculate_channel_transforms_3d(main_channel)?)
|
|
}
|
|
let channels = transforms
|
|
.iter()
|
|
.flat_map(|t| t.keys())
|
|
.collect::<HashSet<_>>();
|
|
let mut average_transforms = HashMap::new();
|
|
for channel in channels {
|
|
let (s, n) = transforms
|
|
.iter()
|
|
.filter_map(|t| {
|
|
if t.contains_key(channel) {
|
|
Some((&t[channel] * &default_transform).matrix())
|
|
} else {
|
|
None
|
|
}
|
|
})
|
|
.fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
|
|
let matrix = s / (n as f64);
|
|
let (s, n) = transforms
|
|
.iter()
|
|
.filter_map(|t| {
|
|
if t.contains_key(channel) {
|
|
Some(((&t[channel] * &default_transform).matrix() - &matrix).powi(2))
|
|
} else {
|
|
None
|
|
}
|
|
})
|
|
.fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
|
|
let dmatrix = (s / n as f64).sqrt();
|
|
average_transforms.insert(
|
|
channel.clone(),
|
|
Transform::default()
|
|
.with_matrix(matrix.view())
|
|
.with_dmatrix(dmatrix.view()),
|
|
);
|
|
}
|
|
Ok(average_transforms)
|
|
}
|
|
|
|
pub fn get_transform(
|
|
&self,
|
|
channel: &str,
|
|
t: usize,
|
|
) -> Result<Option<Cow<'_, Transform<IxDyn>>>, Error> {
|
|
match (self.channel.get(channel), self.drift.get(t)) {
|
|
(Some(channel), Some(drift)) => match (channel.ndim(), drift.ndim()) {
|
|
(2, 2) | (3, 3) => Ok(Some(Cow::Owned(drift * channel))),
|
|
(2, 3) => Ok(Some(Cow::Owned(drift.add_dimension(0) * channel))),
|
|
(3, 2) => Ok(Some(Cow::Owned(drift * channel.add_dimension(0)))),
|
|
(2 | 3, n) | (n, _) => Err(Error::TransformationNotImplemented(n)),
|
|
},
|
|
(Some(channel), None) => Ok(Some(Cow::Borrowed(channel))),
|
|
(None, Some(drift)) => Ok(Some(Cow::Borrowed(drift))),
|
|
(None, None) => Ok(None),
|
|
}
|
|
}
|
|
}
|
|
|
|
impl<D: Dimension, R: Reader> View<D, R> {
|
|
pub fn get_transforms(&self) -> &Transforms {
|
|
&self.transforms
|
|
}
|
|
|
|
pub fn set_transforms(&mut self, transforms: Transforms) {
|
|
self.transforms = transforms;
|
|
}
|
|
|
|
pub fn with_transforms(&self, transforms: Transforms) -> Self {
|
|
let mut new = self.clone();
|
|
new.transforms = transforms;
|
|
new
|
|
}
|
|
|
|
pub fn with_transforms_from_yaml(mut self, path: &Path) -> Result<Self, Error> {
|
|
self.set_transforms(Transforms::load(path)?);
|
|
Ok(self)
|
|
}
|
|
|
|
pub fn load_transforms_from_yaml(&mut self, path: &Path) -> Result<(), Error> {
|
|
self.set_transforms(Transforms::load(path)?);
|
|
Ok(())
|
|
}
|
|
|
|
pub fn calculate_channel_transforms_2d(
|
|
&self,
|
|
main_channel: usize,
|
|
) -> Result<HashMap<String, Transform<Ix2>>, Error> {
|
|
let main_max = self
|
|
.slice_cztyx(s![main_channel, .., 0, .., ..])?
|
|
.max_proj(Axis::Z)?
|
|
.as_array::<f64>()?;
|
|
let mut transforms = HashMap::new();
|
|
for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
|
|
let name = metadata
|
|
.name
|
|
.clone()
|
|
.unwrap_or_else(|| format!("channel{}", channel));
|
|
if channel == main_channel {
|
|
transforms.insert(name, Transform::default());
|
|
} else {
|
|
let max = self
|
|
.slice_cztyx(s![channel, .., 0, .., ..])?
|
|
.max_proj(Axis::Z)?
|
|
.as_array::<f64>()?;
|
|
transforms.insert(
|
|
name,
|
|
Transform::register_affine(main_max.view(), max.view())?,
|
|
);
|
|
}
|
|
}
|
|
Ok(transforms)
|
|
}
|
|
|
|
pub fn calculate_channel_transforms_3d(
|
|
&self,
|
|
main_channel: usize,
|
|
) -> Result<HashMap<String, Transform<Ix3>>, Error> {
|
|
let main_max = self
|
|
.slice_cztyx(s![main_channel, .., 0, .., ..])?
|
|
.as_array::<f64>()?;
|
|
let mut transforms = HashMap::new();
|
|
for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
|
|
let name = metadata
|
|
.name
|
|
.clone()
|
|
.unwrap_or_else(|| format!("channel{}", channel));
|
|
if channel == main_channel {
|
|
transforms.insert(name, Transform::default());
|
|
} else {
|
|
let max = self
|
|
.slice_cztyx(s![channel, .., 0, .., ..])?
|
|
.as_array::<f64>()?;
|
|
transforms.insert(
|
|
name,
|
|
Transform::register_affine(main_max.view(), max.view())?,
|
|
);
|
|
}
|
|
}
|
|
Ok(transforms)
|
|
}
|
|
|
|
pub fn calculate_drift_transforms_2d(&self) -> Result<Vec<Transform<Ix2>>, Error> {
|
|
let mut view = self.clone().into_dyn();
|
|
let axes = view.axes();
|
|
if axes.contains(&Axis::Z) {
|
|
view = view.max_proj(Axis::Z)?;
|
|
}
|
|
if axes.contains(&Axis::C) {
|
|
view = view.max_proj(Axis::C)?;
|
|
}
|
|
view = view.permute_axes(&[Axis::T, Axis::Y, Axis::X])?;
|
|
let n = view.shape().t;
|
|
let m = (n as f64).sqrt().round() as usize;
|
|
let t_groups = (0..=m)
|
|
.map(|i| (((i * n) as f64) / (m as f64)).round() as usize)
|
|
.collect::<Vec<_>>();
|
|
let t_keys = t_groups
|
|
.array_windows()
|
|
.map(|[start, end]| (start + end) / 2)
|
|
.collect::<Vec<_>>();
|
|
let t_pairs = t_keys
|
|
.iter()
|
|
.zip(t_groups.array_windows())
|
|
.flat_map(|(t_key, [start, end])| (*start..*end).map(|i| (*t_key, i)))
|
|
.chain(t_keys.array_windows().map(|[start, end]| (*start, *end)))
|
|
.collect::<Vec<_>>();
|
|
let fmaxz_keys = t_keys
|
|
.iter()
|
|
.map(|t| {
|
|
Ok((
|
|
t,
|
|
filter::gaussian_smooth(
|
|
&view.slice(s![*t, .., ..])?.as_array::<f64>()?,
|
|
&[5.0; 2],
|
|
)?,
|
|
))
|
|
})
|
|
.collect::<Result<HashMap<_, _>, Error>>()?;
|
|
let shifts = t_pairs
|
|
.into_par_iter()
|
|
.map(|(t_key, t)| {
|
|
if t_key == t {
|
|
Ok(Transform::<Ix2>::from_translation(&vec![0.0; view.ndim()]))
|
|
} else {
|
|
if let Some(fmaxz) = fmaxz_keys.get(&t_key) {
|
|
Ok(Transform::<Ix2>::register_translation(
|
|
fmaxz_keys[&t_key].view(),
|
|
fmaxz,
|
|
)?)
|
|
} else {
|
|
let fmaxz = filter::gaussian_smooth(
|
|
&view.slice(s![t, .., ..])?.as_array::<f64>()?,
|
|
&[5.0; 2],
|
|
)?;
|
|
Ok(Transform::<Ix2>::register_translation(
|
|
fmaxz_keys[&t_key].view(),
|
|
fmaxz.view(),
|
|
)?)
|
|
}
|
|
}
|
|
})
|
|
.collect::<Result<Vec<_>, Error>>()?;
|
|
let mut t = Transform::<Ix2>::from_translation(&[0.0; 2]);
|
|
let mut key_shifts = vec![t.clone()];
|
|
key_shifts.extend(shifts.iter().skip(n).map(|i| {
|
|
t = i * &t;
|
|
t.clone()
|
|
}));
|
|
Ok(key_shifts
|
|
.iter()
|
|
.zip(t_groups.array_windows())
|
|
.flat_map(|(s, [start, end])| {
|
|
let mut t = s.clone();
|
|
(*start..*end)
|
|
.map(|i| {
|
|
t = &shifts[i] * &t;
|
|
t.clone()
|
|
})
|
|
.collect::<Vec<_>>()
|
|
})
|
|
.collect())
|
|
}
|
|
|
|
pub fn calculate_drift_transforms_3d(&self) -> Result<Vec<Transform<Ix3>>, Error> {
|
|
todo!()
|
|
}
|
|
|
|
pub(crate) fn get_transformed_frame(
|
|
&self,
|
|
c: usize,
|
|
z: usize,
|
|
t: usize,
|
|
transform: &Transform<IxDyn>,
|
|
order: usize,
|
|
default_value: &DefaultValue,
|
|
) -> Result<Arc<Frame>, Error> {
|
|
match transform.ndim() {
|
|
2 => {
|
|
let tr = transform.clone().into_dimensionality::<Ix2>()?;
|
|
|
|
macro_rules! interpolate {
|
|
($array:expr ,$var:ident, bool) => {
|
|
Ok(Arc::from(Frame::$var(
|
|
tr.interpolate_par_order(
|
|
&$array.mapv(|i| if i { 1.0 } else { 0.0 }),
|
|
order,
|
|
default_value,
|
|
)?
|
|
.mapv(|i| i > 0.5),
|
|
)))
|
|
};
|
|
($array:expr ,$var:ident, f64) => {
|
|
Ok(Arc::from(Frame::$var(tr.interpolate_par_order(
|
|
$array,
|
|
order,
|
|
default_value,
|
|
)?)))
|
|
};
|
|
($array:expr ,$var:ident, $t:ty) => {
|
|
Ok(Arc::from(Frame::$var(
|
|
tr.interpolate_par_order($array, order, default_value)?
|
|
.mapv(|i| i as $t),
|
|
)))
|
|
};
|
|
}
|
|
|
|
match self.get_cached_frame(c, z, t, None)?.as_ref() {
|
|
Frame::Bool(frame) => interpolate!(frame, Bool, bool),
|
|
Frame::I8(frame) => interpolate!(frame, I8, i8),
|
|
Frame::U8(frame) => interpolate!(frame, U8, u8),
|
|
Frame::I16(frame) => interpolate!(frame, I16, i16),
|
|
Frame::U16(frame) => interpolate!(frame, U16, u16),
|
|
Frame::I32(frame) => interpolate!(frame, I32, i32),
|
|
Frame::U32(frame) => interpolate!(frame, U32, u32),
|
|
Frame::U64(frame) => interpolate!(frame, U64, u64),
|
|
Frame::I64(frame) => interpolate!(frame, I64, i64),
|
|
Frame::I128(frame) => interpolate!(frame, I128, i128),
|
|
Frame::U128(frame) => interpolate!(frame, U128, u128),
|
|
Frame::F32(frame) => interpolate!(frame, F32, f32),
|
|
Frame::F64(frame) => interpolate!(frame, F64, f64),
|
|
Frame::F128(frame) => interpolate!(frame, F128, f64),
|
|
}
|
|
}
|
|
3 => {
|
|
let tr = transform.clone().into_dimensionality::<Ix3>()?;
|
|
let shape = self.shape();
|
|
let center = shape.iter().map(|&i| (i as f64) / 2.0).collect::<Vec<_>>();
|
|
let zf = z as f64;
|
|
let yf = shape.y as f64 - 1.0;
|
|
let xf = shape.x as f64 - 1.0;
|
|
let corners = array![[zf, 0.0, 0.0], [zf, 0.0, xf], [zf, yf, 0.0], [zf, yf, xf]];
|
|
let extent = tr.transform_points(&corners, ¢er)?;
|
|
let z_extent = extent.slice(s![.., 2]);
|
|
let z_min = (z_extent
|
|
.iter()
|
|
.min_by(|a, b| a.total_cmp(b))
|
|
.unwrap()
|
|
.floor() as isize
|
|
- 2)
|
|
.max(0) as usize;
|
|
let z_max = (z_extent
|
|
.iter()
|
|
.max_by(|a, b| a.total_cmp(b))
|
|
.unwrap()
|
|
.ceil() as isize
|
|
+ 2)
|
|
.min(shape.z as isize - 1)
|
|
.max(0) as usize;
|
|
let frames = (z_min..=z_max)
|
|
.map(|zi| self.get_cached_frame(c, zi, t, None))
|
|
.collect::<Result<Vec<_>, Error>>()?;
|
|
|
|
// TODO; check these translations
|
|
let u0 = tr.as_translation(&[z_min as f64, 0.0, 0.0]);
|
|
let u1 = tr.as_translation(&[-(z_min as f64), 0.0, 0.0]);
|
|
let tr = u0 * tr * u1;
|
|
|
|
match frames.first().map(|f| f.as_ref()) {
|
|
Some(first) => {
|
|
macro_rules! interpolate {
|
|
($var:ident, bool) => {{
|
|
let views = frames
|
|
.iter()
|
|
.map(|f| match f.as_ref() {
|
|
Frame::$var(a) => Ok(a.view()),
|
|
_ => unreachable!(),
|
|
})
|
|
.collect::<Result<Vec<_>, Error>>()?;
|
|
let stack = ndarray::stack(ndarray::Axis(0), &views)?;
|
|
let transformed = tr.interpolate_par_order(
|
|
&stack.mapv(|i| if i { 1.0 } else { 0.0 }),
|
|
order,
|
|
default_value,
|
|
)?;
|
|
Ok(Arc::from(Frame::$var(
|
|
transformed.slice(s![0, .., ..]).mapv(|i| i > 0.5),
|
|
)))
|
|
}};
|
|
($var:ident, f64) => {{
|
|
let views = frames
|
|
.iter()
|
|
.map(|f| match f.as_ref() {
|
|
Frame::$var(a) => Ok(a.view()),
|
|
_ => unreachable!(),
|
|
})
|
|
.collect::<Result<Vec<_>, Error>>()?;
|
|
let stack = ndarray::stack(ndarray::Axis(0), &views)?;
|
|
let transformed =
|
|
tr.interpolate_par_order(&stack, order, default_value)?;
|
|
Ok(Arc::from(Frame::$var(
|
|
transformed.slice(s![0, .., ..]).to_owned(),
|
|
)))
|
|
}};
|
|
($var:ident, $t:ty) => {{
|
|
let views = frames
|
|
.iter()
|
|
.map(|f| match f.as_ref() {
|
|
Frame::$var(a) => Ok(a.view()),
|
|
_ => unreachable!(),
|
|
})
|
|
.collect::<Result<Vec<_>, Error>>()?;
|
|
let stack = ndarray::stack(ndarray::Axis(0), &views)?;
|
|
let transformed =
|
|
tr.interpolate_par_order(&stack, order, default_value)?;
|
|
Ok(Arc::from(Frame::$var(
|
|
transformed.slice(s![0, .., ..]).mapv(|i| i as $t),
|
|
)))
|
|
}};
|
|
}
|
|
match first {
|
|
Frame::Bool(_) => interpolate!(Bool, bool),
|
|
Frame::I8(_) => interpolate!(I8, i8),
|
|
Frame::U8(_) => interpolate!(U8, u8),
|
|
Frame::I16(_) => interpolate!(I16, i16),
|
|
Frame::U16(_) => interpolate!(U16, u16),
|
|
Frame::I32(_) => interpolate!(I32, i32),
|
|
Frame::U32(_) => interpolate!(U32, u32),
|
|
Frame::F32(_) => interpolate!(F32, f32),
|
|
Frame::F64(_) => interpolate!(F64, f64),
|
|
Frame::I64(_) => interpolate!(I64, i64),
|
|
Frame::U64(_) => interpolate!(U64, u64),
|
|
Frame::I128(_) => interpolate!(I128, i128),
|
|
Frame::U128(_) => interpolate!(U128, u128),
|
|
Frame::F128(_) => interpolate!(F128, f64),
|
|
}
|
|
}
|
|
None => {
|
|
let shape = self.shape();
|
|
Ok(Arc::new(match self.pixel_type() {
|
|
PixelType::Bool => Frame::Bool(Array2::default((shape.y, shape.x))),
|
|
PixelType::I8 => Frame::I8(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::U8 => Frame::U8(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::I16 => Frame::I16(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::U16 => Frame::U16(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::I32 => Frame::I32(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::U32 => Frame::U32(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::F32 => Frame::F32(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::F64 => Frame::F64(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::I64 => Frame::I64(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::U64 => Frame::U64(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::I128 => Frame::I128(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::U128 => Frame::U128(Array2::zeros((shape.y, shape.x))),
|
|
PixelType::F128 => Frame::F128(Array2::zeros((shape.y, shape.x))),
|
|
}))
|
|
}
|
|
}
|
|
}
|
|
d => Err(Error::TransformationNotImplemented(d)),
|
|
}
|
|
}
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use crate::axes::Axis;
|
|
use crate::error::Error;
|
|
use crate::readers::{DynReader, Reader};
|
|
use crate::transforms::Transforms;
|
|
use image_registration::transform::Transform;
|
|
use ndarray::s;
|
|
use std::collections::HashMap;
|
|
|
|
fn open(file: &str) -> Result<DynReader, Error> {
|
|
let path = std::env::current_dir()?
|
|
.join("tests")
|
|
.join("files")
|
|
.join(file);
|
|
DynReader::new(&path, 0, 0)
|
|
}
|
|
|
|
#[test]
|
|
fn test_transforms() -> Result<(), Box<dyn std::error::Error>> {
|
|
let _t = Transforms::default();
|
|
|
|
Ok(())
|
|
}
|
|
|
|
#[test]
|
|
fn load_transforms() -> Result<(), Box<dyn std::error::Error>> {
|
|
let t = Transforms::load("tests/files/transforms.yml")?;
|
|
println!("t: {:?}", t);
|
|
assert_eq!(t.channel["TV1"].ndim(), 2);
|
|
assert_eq!(t.channel["TV2"].ndim(), 2);
|
|
Ok(())
|
|
}
|
|
|
|
#[test]
|
|
fn load_old_transforms() -> Result<(), Box<dyn std::error::Error>> {
|
|
let t = Transforms::load("tests/files/transform_red.yml")?;
|
|
println!("t: {:?}", t);
|
|
assert_eq!(t.channel["TV1"].ndim(), 2);
|
|
assert_eq!(t.channel["TV2"].ndim(), 2);
|
|
Ok(())
|
|
}
|
|
|
|
#[test]
|
|
fn transforms_as_array() -> Result<(), Box<dyn std::error::Error>> {
|
|
let file = "czi/beads_2020_08_05__11_43_07.czi";
|
|
let reader = open(file)?;
|
|
let view = reader.view();
|
|
let mut channel = HashMap::new();
|
|
channel.insert(
|
|
"TV1".to_string(),
|
|
Transform::new(vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0]),
|
|
);
|
|
channel.insert(
|
|
"TV2".to_string(),
|
|
Transform::new(vec![
|
|
0.9995868322142533,
|
|
0.012140955373934502,
|
|
0.012592501668218145,
|
|
0.8571690967587099,
|
|
-5.270859902118885,
|
|
-5.81118547158953,
|
|
]),
|
|
);
|
|
let view_t = view.with_transforms(Transforms {
|
|
default_value: Default::default(),
|
|
interpolation_order: None,
|
|
channel,
|
|
drift: Vec::new(),
|
|
});
|
|
|
|
let s_a = view.max_proj(Axis::Z)?.slice_cztyx(s![1, 0, .., .., ..])?;
|
|
|
|
println!(
|
|
"s_a: {:?}, {:?}, {:?}",
|
|
s_a.get_slice(),
|
|
s_a.get_axes(),
|
|
s_a.get_operations()
|
|
);
|
|
|
|
let s_b = view_t
|
|
.max_proj(Axis::Z)?
|
|
.slice_cztyx(s![1, 0, .., .., ..])?;
|
|
|
|
println!(
|
|
"s_b: {:?}, {:?}, {:?}",
|
|
s_b.get_slice(),
|
|
s_b.get_axes(),
|
|
s_b.get_operations()
|
|
);
|
|
|
|
let a = s_a.as_array::<u16>()?;
|
|
let b = s_b.as_array::<u16>()?;
|
|
assert_ne!(a, b);
|
|
|
|
// let a = view.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
|
|
// let b = view_t.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
|
|
// assert_ne!(a, b);
|
|
//
|
|
// let a = view
|
|
// .max_proj(Axis::Z)?
|
|
// .slice(s![1, 0, .., ..])?
|
|
// .as_array::<u16>()?;
|
|
// let b = view_t
|
|
// .max_proj(Axis::Z)?
|
|
// .slice(s![1, 0, .., ..])?
|
|
// .as_array::<u16>()?;
|
|
// assert_ne!(a, b);
|
|
|
|
Ok(())
|
|
}
|
|
}
|