Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add order arg to shift and zoom #24

Open
wants to merge 13 commits into
base: master
Choose a base branch
from
Open

Add order arg to shift and zoom #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
f6fdd94
Add order arg to zoom and shift
nilgoyette Sep 17, 2023
d25a408
Order 0
nilgoyette Sep 17, 2023
ec05f83
Add tests for order 0
nilgoyette Sep 17, 2023
7113b81
Add a Nearest failing test
nilgoyette Sep 17, 2023
ac27075
Refactor run_zoom_shift to fix Nearest
nilgoyette Sep 17, 2023
3fed594
Add failing tests for order=1
nilgoyette Sep 17, 2023
be1a30b
Refactor interpolate arrays for various orders
nilgoyette Sep 17, 2023
a832918
Add doc for the order=0 case
nilgoyette Sep 17, 2023
9cdc973
Add failing tests for order=2
nilgoyette Sep 17, 2023
0181f2e
Fix order=2
nilgoyette Sep 17, 2023
fe049c5
Add tests for order=4
nilgoyette Sep 17, 2023
b8dc55d
Add failing tests for order=5
nilgoyette Sep 17, 2023
9c3ef83
Fix order=5
nilgoyette Sep 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
163 changes: 98 additions & 65 deletions src/interpolation/zoom_shift.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
use std::ops::{Add, Sub};

use ndarray::{Array, Array2, ArrayBase, Data, Ix3, Zip};
use ndarray::{s, Array, Array2, ArrayBase, ArrayViewMut1, Data, Ix3, Zip};
use num_traits::{FromPrimitive, Num, ToPrimitive};

use crate::{array_like, pad, round_ties_even, spline_filter, BorderMode, PadMode};
Expand All @@ -12,6 +12,7 @@ use crate::{array_like, pad, round_ties_even, spline_filter, BorderMode, PadMode
///
/// * `data` - A 3D array of the data to shift.
/// * `shift` - The shift along the axes.
/// * `order` - The order of the spline.
/// * `mode` - The mode parameter determines how the input array is extended beyond its boundaries.
/// * `prefilter` - Determines if the input array is prefiltered with spline_filter before
/// interpolation. The default is `true`, which will create a temporary `f64` array of filtered
Expand All @@ -20,20 +21,17 @@ use crate::{array_like, pad, round_ties_even, spline_filter, BorderMode, PadMode
pub fn shift<S, A>(
data: &ArrayBase<S, Ix3>,
shift: [f64; 3],
order: usize,
mode: BorderMode<A>,
prefilter: bool,
) -> Array<A, Ix3>
where
S: Data<Elem = A>,
A: Copy + Num + FromPrimitive + PartialOrd + ToPrimitive,
{
let order = 3;
let dim = [data.dim().0, data.dim().1, data.dim().2];
let shift = shift.map(|s| -s);
let reslicer = ZoomShiftReslicer::new(dim, dim, [1.0, 1.0, 1.0], shift, order, mode);

let out = array_like(&data, data.raw_dim(), A::zero());
run_zoom_shift(data, order, mode, prefilter, &reslicer, out)
run_zoom_shift(data, dim, [1.0, 1.0, 1.0], shift, order, mode, prefilter)
}

/// Zoom an array.
Expand All @@ -42,6 +40,7 @@ where
///
/// * `data` - A 3D array of the data to zoom
/// * `zoom` - The zoom factor along the axes.
/// * `order` - The order of the spline.
/// * `mode` - The mode parameter determines how the input array is extended beyond its boundaries.
/// * `prefilter` - Determines if the input array is prefiltered with spline_filter before
/// interpolation. The default is `true`, which will create a temporary `f64` array of filtered
Expand All @@ -50,15 +49,14 @@ where
pub fn zoom<S, A>(
data: &ArrayBase<S, Ix3>,
zoom: [f64; 3],
order: usize,
mode: BorderMode<A>,
prefilter: bool,
) -> Array<A, Ix3>
where
S: Data<Elem = A>,
A: Copy + Num + FromPrimitive + PartialOrd + ToPrimitive,
{
let order = 3;
let i_dim = [data.dim().0, data.dim().1, data.dim().2];
let mut o_dim = data.raw_dim();
for (ax, (&ax_len, zoom)) in data.shape().iter().zip(zoom.iter()).enumerate() {
o_dim[ax] = round_ties_even(ax_len as f64 * zoom) as usize;
Expand All @@ -76,47 +74,51 @@ where
nom[1] as f64 / div[1] as f64,
nom[2] as f64 / div[2] as f64,
];
let reslicer = ZoomShiftReslicer::new(i_dim, o_dim, zoom, [0.0, 0.0, 0.0], order, mode);

let out = array_like(&data, o_dim, A::zero());
run_zoom_shift(data, order, mode, prefilter, &reslicer, out)
run_zoom_shift(data, o_dim, zoom, [0.0, 0.0, 0.0], order, mode, prefilter)
}

fn run_zoom_shift<S, A>(
data: &ArrayBase<S, Ix3>,
odim: [usize; 3],
zooms: [f64; 3],
shifts: [f64; 3],
order: usize,
mode: BorderMode<A>,
prefilter: bool,
reslicer: &ZoomShiftReslicer,
mut out: Array<A, Ix3>,
) -> Array<A, Ix3>
where
S: Data<Elem = A>,
A: Copy + Num + FromPrimitive + PartialOrd + ToPrimitive,
{
let idim = [data.dim().0, data.dim().1, data.dim().2];
let mut out = array_like(&data, odim, A::zero());
if prefilter && order > 1 {
let data = match mode {
// We need to allocate and work on filtered data
let (data, nb_prepad) = match mode {
BorderMode::Nearest => {
let padded = pad(data, &[[12, 12]], PadMode::Edge);
spline_filter(&padded, order, mode)
(spline_filter(&padded, order, mode), 12)
}
_ => spline_filter(data, order, mode),
_ => (spline_filter(data, order, mode), 0),
};
let reslicer = ZoomShiftReslicer::new(idim, odim, zooms, shifts, order, mode, nb_prepad);
Zip::indexed(&mut out).for_each(|idx, o| {
*o = A::from_f64(reslicer.interpolate(&data, idx)).unwrap();
});
} else {
// We can use the &data as-is
let reslicer = ZoomShiftReslicer::new(idim, odim, zooms, shifts, order, mode, 0);
Zip::indexed(&mut out).for_each(|idx, o| {
*o = A::from_f64(reslicer.interpolate(&data, idx)).unwrap();
*o = A::from_f64(reslicer.interpolate(data, idx)).unwrap();
});
}
out
}

/// Zoom shift transformation (only scaling and translation).
///
/// Hardcoded to use spline interpolation of order 3 and mirror mode.
struct ZoomShiftReslicer {
order: usize,
offsets: [Vec<isize>; 3],
edge_offsets: [Array2<isize>; 3],
is_edge_case: [Vec<bool>; 3],
Expand All @@ -134,6 +136,7 @@ impl ZoomShiftReslicer {
shifts: [f64; 3],
order: usize,
mode: BorderMode<A>,
nb_prepad: isize,
) -> ZoomShiftReslicer
where
A: Copy + ToPrimitive,
Expand All @@ -160,8 +163,8 @@ impl ZoomShiftReslicer {
};

let mut reslicer =
ZoomShiftReslicer { offsets, edge_offsets, is_edge_case, splvals, zeros, cval };
reslicer.build_arrays(idim, odim, zooms, shifts, order, mode);
ZoomShiftReslicer { order, offsets, edge_offsets, is_edge_case, splvals, zeros, cval };
reslicer.build_arrays(idim, odim, zooms, shifts, order, mode, nb_prepad);
reslicer
}

Expand All @@ -173,6 +176,7 @@ impl ZoomShiftReslicer {
shifts: [f64; 3],
order: usize,
mode: BorderMode<A>,
nb_prepad: isize,
) where
A: Copy,
{
Expand All @@ -181,10 +185,6 @@ impl ZoomShiftReslicer {
BorderMode::Constant(_) | BorderMode::Wrap => BorderMode::Mirror,
_ => mode,
};
let nb_prepad = match mode {
BorderMode::Nearest => 12,
_ => 0,
};
let iorder = order as isize;
let idim = [
idim[0] as isize + 2 * nb_prepad,
Expand All @@ -207,12 +207,13 @@ impl ZoomShiftReslicer {
_ => to = map_coordinates(to, idim[axis] as f64, mode),
};
if to > -1.0 {
if order > 0 {
build_splines(to, &mut splvals.row_mut(from), order);
}
if order & 1 == 0 {
to += 0.5;
}

build_splines(from, to, splvals, order);

let start = to.floor() as isize - iorder / 2;
offsets[from] = start;
if start < 0 || start + iorder >= idim[axis] {
Expand Down Expand Up @@ -240,46 +241,46 @@ impl ZoomShiftReslicer {
return self.cval;
}

// Linear interpolation use a 4x4x4 block. This is simple enough, but we must adjust this
// Order = 0
// We do not want to go further because
// - it would be uselessly slower
// - self.splvals is empty so it would crash (although we could fill it with 1.0)
if self.edge_offsets[0].is_empty() {
let x = self.offsets[0][start.0] as usize;
let y = self.offsets[1][start.1] as usize;
let z = self.offsets[2][start.2] as usize;
return data[(x, y, z)].to_f64().unwrap();
}

// Linear interpolation use a nxnxn block. This is simple enough, but we must adjust this
// block when the `start` is near the edges.
let n = self.order + 1;
let valid_index = |original_offset, is_edge, start, d: usize, v| {
(original_offset + if is_edge { self.edge_offsets[d][(start, v)] } else { v as isize })
as usize
};

let original_offset = self.offsets[0][start.0];
let is_edge = self.is_edge_case[0][start.0];
let xs = [
valid_index(original_offset, is_edge, start.0, 0, 0),
valid_index(original_offset, is_edge, start.0, 0, 1),
valid_index(original_offset, is_edge, start.0, 0, 2),
valid_index(original_offset, is_edge, start.0, 0, 3),
];

let original_offset = self.offsets[1][start.1];
let is_edge = self.is_edge_case[1][start.1];
let ys = [
valid_index(original_offset, is_edge, start.1, 1, 0),
valid_index(original_offset, is_edge, start.1, 1, 1),
valid_index(original_offset, is_edge, start.1, 1, 2),
valid_index(original_offset, is_edge, start.1, 1, 3),
];

let original_offset = self.offsets[2][start.2];
let is_edge = self.is_edge_case[2][start.2];
let zs = [
valid_index(original_offset, is_edge, start.2, 2, 0),
valid_index(original_offset, is_edge, start.2, 2, 1),
valid_index(original_offset, is_edge, start.2, 2, 2),
valid_index(original_offset, is_edge, start.2, 2, 3),
];
let original_offset_x = self.offsets[0][start.0];
let is_edge_x = self.is_edge_case[0][start.0];
let mut xs = [0; 6];
let original_offset_y = self.offsets[1][start.1];
let is_edge_y = self.is_edge_case[1][start.1];
let mut ys = [0; 6];
let original_offset_z = self.offsets[2][start.2];
let is_edge_z = self.is_edge_case[2][start.2];
let mut zs = [0; 6];
for i in 0..n {
xs[i] = valid_index(original_offset_x, is_edge_x, start.0, 0, i);
ys[i] = valid_index(original_offset_y, is_edge_y, start.1, 1, i);
zs[i] = valid_index(original_offset_z, is_edge_z, start.2, 2, i);
}

let mut t = 0.0;
for (z, &idx_z) in zs.iter().enumerate() {
for (z, &idx_z) in zs[..n].iter().enumerate() {
let spline_z = self.splvals[2][(start.2, z)];
for (y, &idx_y) in ys.iter().enumerate() {
for (y, &idx_y) in ys[..n].iter().enumerate() {
let spline_yz = self.splvals[1][(start.1, y)] * spline_z;
for (x, &idx_x) in xs.iter().enumerate() {
for (x, &idx_x) in xs[..n].iter().enumerate() {
let spline_xyz = self.splvals[0][(start.0, x)] * spline_yz;
t += data[(idx_x, idx_y, idx_z)].to_f64().unwrap() * spline_xyz;
}
Expand All @@ -289,14 +290,46 @@ impl ZoomShiftReslicer {
}
}

fn build_splines(from: usize, to: f64, spline: &mut Array2<f64>, order: usize) {
let x = to - to.floor();
let y = x;
let z = 1.0 - x;
spline[(from, 0)] = z * z * z / 6.0;
spline[(from, 1)] = (y * y * (y - 2.0) * 3.0 + 4.0) / 6.0;
spline[(from, 2)] = (z * z * (z - 2.0) * 3.0 + 4.0) / 6.0;
spline[(from, order)] = 1.0 - (spline[(from, 0)] + spline[(from, 1)] + spline[(from, 2)]);
fn build_splines(to: f64, spline: &mut ArrayViewMut1<f64>, order: usize) {
let x = to - if order & 1 == 1 { to } else { to + 0.5 }.floor();
match order {
1 => spline[0] = 1.0 - x,
2 => {
spline[0] = 0.5 * (0.5 - x).powi(2);
spline[1] = 0.75 - x * x;
}
3 => {
let z = 1.0 - x;
spline[0] = z * z * z / 6.0;
spline[1] = (x * x * (x - 2.0) * 3.0 + 4.0) / 6.0;
spline[2] = (z * z * (z - 2.0) * 3.0 + 4.0) / 6.0;
}
4 => {
let t = x * x;
let y = 1.0 + x;
let z = 1.0 - x;
spline[0] = (0.5 - x).powi(4) / 24.0;
spline[1] = y * (y * (y * (5.0 - y) / 6.0 - 1.25) + 5.0 / 24.0) + 55.0 / 96.0;
spline[2] = t * (t * 0.25 - 0.625) + 115.0 / 192.0;
spline[3] = z * (z * (z * (5.0 - z) / 6.0 - 1.25) + 5.0 / 24.0) + 55.0 / 96.0;
}
5 => {
let y = 1.0 - x;
let t = y * y;
spline[0] = y * t * t / 120.0;
let y = x + 1.0;
spline[1] = y * (y * (y * (y * (y / 24.0 - 0.375) + 1.25) - 1.75) + 0.625) + 0.425;
let t = x * x;
spline[2] = t * (t * (0.25 - x / 12.0) - 0.5) + 0.55;
let z = 1.0 - x;
let t = z * z;
spline[3] = t * (t * (0.25 - z / 12.0) - 0.5) + 0.55;
let z = z + 1.0;
spline[4] = z * (z * (z * (z * (z / 24.0 - 0.375) + 1.25) - 1.75) + 0.625) + 0.425;
}
_ => panic!("order must be between 1 and 5"),
}
spline[order] = 1.0 - spline.slice(s![..order]).sum();
}

fn map_coordinates<A>(mut idx: f64, len: f64, mode: BorderMode<A>) -> f64 {
Expand Down
Loading