Compare commits
10
Commits
aac61b981f
...
master
| Author | SHA1 | Date | |
|---|---|---|---|
 w.pomp
|
2b01f65cc9 | ||
|
w.pomp
|
e5136f00bb | ||
|
w.pomp
|
fccc3e6f1d | ||
|
w.pomp
|
0589b988a1 | ||
|
w.pomp
|
3ff4ffc6ea | ||
|
w.pomp
|
60a4f3a542 | ||
|
w.pomp
|
1c0bc4aa60 | ||
|
w.pomp
|
0049267997 | ||
|
Wim Pomp
|
9257dd1fcf | ||
|
Wim Pomp
|
b12982da48 |
No files matched your search
+2
-3
@@ -1,13 +1,12 @@
|
||||
[package]
|
||||
name = "image-registration"
|
||||
version = "0.1.0"
|
||||
version = "0.1.1"
|
||||
edition = "2024"
|
||||
rust-version = "1.88.0"
|
||||
authors = ["Wim Pomp <w.pomp@nki.nl>"]
|
||||
license = "MIT OR Apache-2.0"
|
||||
description = "Find and apply affine transformations."
|
||||
homepage = "https://git.pomppervova.nl/wim/image-registration/src/branch/rs"
|
||||
repository = "https://git.pomppervova.nl/wim/image-registration/src/branch/rs"
|
||||
repository = "https://git.pomppervova.nl/wim/image-registration"
|
||||
documentation = "https://docs.rs/image-registration"
|
||||
readme = "README.md"
|
||||
keywords = ["affine", "transformation", "ndarray"]
|
||||
|
||||
@@ -4,7 +4,7 @@ Find an ND affine transformation that registers an ND test image to a fixed imag
|
||||
other images.
|
||||
|
||||
`image-registration` is a pure-Rust, dependency-light implementation of intensity-based image registration in the spirit
|
||||
of [SimpleElastix]([https://elastix.lumc.nl](https://simpleelastix.github.io/)). It finds an affine (rigid + scale + shear + translation) transform
|
||||
of [SimpleElastix](https://simpleelastix.github.io/). It finds an affine (rigid + scale + shear + translation) transform
|
||||
that aligns a *moving* image to a *fixed* image by maximizing
|
||||
the [Mattes mutual-information](https://doi.org/10.1109/TMI.2003.809072) metric with a B-spline Parzen
|
||||
window, using a multi-resolution Gaussian image pyramid and a hybrid ASGD + L-BFGS optimization strategy.
|
||||
@@ -12,6 +12,10 @@ window, using a multi-resolution Gaussian image pyramid and a hybrid ASGD + L-BF
|
||||
Everything operates on [`ndarray`](https://docs.rs/ndarray) arrays, works in arbitrary dimensions (1D, 2D, 3D, …), and
|
||||
is parallelized with [rayon](https://docs.rs/rayon).
|
||||
|
||||
So far, 2D registration seems to work, but 3D and higher dimensions often do not. Also, this code is still a lot slower
|
||||
than (Simple)Elastix. The biggest benefit of this crate is that it's Rust and does not require a lengthy compilation of
|
||||
C++ code.
|
||||
|
||||
## Features
|
||||
|
||||
- **Affine registration** in any number of dimensions — rotation, scale, shear, and translation parameters
|
||||
@@ -36,13 +40,13 @@ Register a moving image to a fixed image and obtain the affine transform:
|
||||
```rust
|
||||
use image_registration::error::Error;
|
||||
use image_registration::transform::Transform;
|
||||
use ndarray::Array2;
|
||||
use ndarray::{Array2, Ix2};
|
||||
|
||||
fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform, Error> {
|
||||
fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform<Ix2>, Error> {
|
||||
// All 6 parameters of the 2D affine transform are free (None = optimize).
|
||||
// steps=None uses the default multi-resolution pyramid; initial_guess=None
|
||||
// initializes by aligning the geometric centers of the two images.
|
||||
Transform::<ndarray::Ix2>::register(
|
||||
Transform::<Ix2>::register(
|
||||
fixed.view(),
|
||||
moving.view(),
|
||||
vec![None; 6],
|
||||
@@ -55,9 +59,10 @@ fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform, Erro
|
||||
The returned `Transform` maps coordinates from the *fixed* image to the *moving* image, i.e.
|
||||
`moving(transform_point(x)) ≈ fixed(x)`. To warp the moving image into the fixed image's frame:
|
||||
|
||||
```rust
|
||||
```ignore
|
||||
let transform = register(&fixed, &moving)?;
|
||||
let registered = transform.interpolate_par::<3, _, _>(moving.view())?; // cubic B-spline, parallel
|
||||
// out-of-range pixels are filled with the default value (0.0 by default)
|
||||
let registered = transform.interpolate_par::<3, _, _>(moving.view(), &Default::default())?; // cubic B-spline, parallel
|
||||
```
|
||||
|
||||
## Concepts
|
||||
@@ -71,15 +76,17 @@ in total:
|
||||
- 2D: `[m00, m01, m10, m11, tx, ty]` (row-major linear part, then translation)
|
||||
- 3D: `[m00, m01, m02, m10, m11, m12, m20, m21, m22, tx, ty, tz]`
|
||||
|
||||
The transformation is applied about a center saved in ```Transform```:
|
||||
The transformation is applied about a center `c` that is supplied where the transform is used:
|
||||
|
||||
```text
|
||||
v = A·(p − c) + t + c
|
||||
```
|
||||
|
||||
where `A` is the linear part, `t` the translation, `c` the center of the transformation, and `p` the input coordinates.
|
||||
`Transform::new(parameters, shape)` sets the center to `(shape − 1) / 2` automatically; `Transform::new_with_center`
|
||||
lets you override it.
|
||||
`transform_point`/`transform_points` take the center as an argument, and resampling (`interpolate`, `interpolate_par`)
|
||||
uses the geometric center of the image, `(shape − 1) / 2`. `Transform::new_around_center(parameters, center)`
|
||||
re-expresses the parameters so that the transformation acts about a chosen center (this is how
|
||||
`from_rotation`/`with_rotation` rotate about a given point).
|
||||
|
||||
### Registration convention
|
||||
|
||||
@@ -101,19 +108,24 @@ fn main() -> Result<(), Error> {
|
||||
let rot = Transform::<Ix2>::from_rotation(std::f64::consts::FRAC_PI_4, &[99.5, 99.5]);
|
||||
|
||||
// Compose: translate, then scale, then rotate
|
||||
let t = rot.with_scaling(&[0.9, 1.1]).with_translation(&[10.0, -5.0]);
|
||||
let t = rot.add_scaling(&[0.9, 1.1]).with_translation(&[10.0, -5.0]);
|
||||
|
||||
// Transform a single point (the center is subtracted and re-added)
|
||||
let p = t.transform_point(&[200, 300]);
|
||||
let p = t.transform_point(&[200.0, 300.0], &[99.5, 99.5]);
|
||||
|
||||
// Transform many points at once (rows of an N×2 array)
|
||||
let points = ndarray::array![[0.0, 0.0], [199.0, 199.0]];
|
||||
let pts = t.transform_points(points.view())?;
|
||||
let pts = t.transform_points(points.view(), &[99.5, 99.5])?;
|
||||
|
||||
// Matrix form (3×3 homogeneous) and exact inverse
|
||||
let m = t.matrix();
|
||||
let inv = t.inverse()?;
|
||||
assert!(m.dot(&inv.matrix()).iter().all(|x| (x - 1.0).abs() < 1e-9));
|
||||
let identity = ndarray::Array2::<f64>::eye(3);
|
||||
assert!(m
|
||||
.dot(&inv.matrix())
|
||||
.iter()
|
||||
.zip(identity.iter())
|
||||
.all(|(x, y)| (x - y).abs() < 1e-9));
|
||||
Ok(())
|
||||
}
|
||||
```
|
||||
@@ -132,11 +144,11 @@ fn main() -> Result<(), Error> {
|
||||
let image = ndarray::Array2::<f64>::zeros((200, 150));
|
||||
let transform = Transform::<ndarray::Ix2>::from_rotation(0.3, &[99.5, 74.5]);
|
||||
|
||||
// Order-1 (linear) interpolation, single-threaded
|
||||
let warped_linear = transform.interpolate::<1, _, _>(image.view())?;
|
||||
// Order-1 (linear) interpolation, single-threaded (out-of-range pixels → 0.0)
|
||||
let warped_linear = transform.interpolate::<1, _, _>(image.view(), &Default::default())?;
|
||||
|
||||
// Order-3 (cubic) B-spline interpolation, parallel
|
||||
let warped_cubic = transform.interpolate_par::<3, _, _>(image.view())?;
|
||||
let warped_cubic = transform.interpolate_par::<3, _, _>(image.view(), &Default::default())?;
|
||||
Ok(())
|
||||
}
|
||||
```
|
||||
@@ -146,7 +158,7 @@ fn main() -> Result<(), Error> {
|
||||
Create a moving image by applying a known transform to a fixed image, register, and check that the recovered transform
|
||||
matches the inverse of the ground-truth transform (adapted from the crate's `register2_random_affine` test):
|
||||
|
||||
```rust
|
||||
```no_run
|
||||
use image_registration::error::Error;
|
||||
use image_registration::julia_image;
|
||||
use image_registration::transform::Transform;
|
||||
@@ -163,11 +175,10 @@ fn main() -> Result<(), Error> {
|
||||
let theta = 10.0f64.to_radians();
|
||||
let (s, c) = theta.sin_cos();
|
||||
let params = vec![c * 1.05, -s, s, c * 0.95, 12.0, -8.0];
|
||||
let ground_truth =
|
||||
Transform::<Ix2>::new_with_center(params.clone(), center.to_vec(), shape.to_vec());
|
||||
let ground_truth = Transform::<Ix2>::new_around_center(params.clone(), ¢er);
|
||||
|
||||
// Moving image = ground truth applied to the fixed image
|
||||
let moving = ground_truth.interpolate::<3, _, _>(fixed.view())?;
|
||||
let moving = ground_truth.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
|
||||
|
||||
// Register (all 6 parameters free, default steps)
|
||||
let t = Transform::<Ix2>::register(fixed.view(), moving.view(), vec![None; 6], None, None)?;
|
||||
@@ -185,9 +196,13 @@ fn main() -> Result<(), Error> {
|
||||
}
|
||||
```
|
||||
|
||||
> **Note:** registration uses random point sampling on its coarser levels, so whether `sse < 1.0` holds for a given
|
||||
> image pair and transform is stochastic; this example may need a few tries (or a simpler transform, e.g. pure
|
||||
> translation) to converge. The crate's `register2_random_affine` unit test is a verified case for this path.
|
||||
|
||||
### 4. Registration with inspection of every optimization step
|
||||
|
||||
```rust
|
||||
```ignore
|
||||
use image_registration::error::Error;
|
||||
use image_registration::transform::Transform;
|
||||
|
||||
@@ -212,13 +227,13 @@ fn main() -> Result<(), Error> {
|
||||
following uses two levels: a coarse level with strong Gaussian smoothing and the robust ASGD optimizer, then a fine
|
||||
level with L-BFGS for accuracy.
|
||||
|
||||
```rust
|
||||
```ignore
|
||||
use image_registration::error::Error;
|
||||
use image_registration::metric::{SamplingArg, Sigma};
|
||||
use image_registration::register::{Optimizer, Registration, RegistrationStep};
|
||||
|
||||
fn main() -> Result<(), Error> {
|
||||
let mut coarse = RegistrationStep::new(
|
||||
let coarse = RegistrationStep::new(
|
||||
Sigma::Absolute(vec![4.0, 4.0]), // Gaussian sigma per dimension
|
||||
SamplingArg::Fixed(3000), // 3000 sampled points (cached)
|
||||
32, // histogram bins
|
||||
@@ -226,10 +241,10 @@ fn main() -> Result<(), Error> {
|
||||
0.05, // edge fraction
|
||||
512, // max iterations
|
||||
1.0, // learning rate
|
||||
Optimizer::ASGD, // robust to noisy gradients
|
||||
);
|
||||
coarse.optimizer = Optimizer::ASGD; // robust to noisy gradients
|
||||
|
||||
let mut fine = RegistrationStep::new(
|
||||
let fine = RegistrationStep::new(
|
||||
Sigma::Absolute(vec![0.5, 0.5]),
|
||||
SamplingArg::Fixed(3000),
|
||||
32,
|
||||
@@ -237,9 +252,10 @@ fn main() -> Result<(), Error> {
|
||||
0.05,
|
||||
2048,
|
||||
1.0,
|
||||
); // optimizer stays LBFGS (the default)
|
||||
Optimizer::LBFGS, // the default optimizer
|
||||
);
|
||||
|
||||
let mut reg = Registration::new(vec![None; 6]);
|
||||
let mut reg = Registration::new(vec![None; 6])?;
|
||||
reg.set_steps(vec![coarse, fine]);
|
||||
let t = reg.register(fixed.view(), moving.view())?;
|
||||
Ok(())
|
||||
@@ -254,7 +270,7 @@ fn main() -> Result<(), Error> {
|
||||
Pass a `Vec<Option<f64>>` as `fixed_mu`: `None` optimizes a parameter, `Some(v)` keeps it fixed at `v`. For example,
|
||||
register for translation only:
|
||||
|
||||
```rust
|
||||
```ignore
|
||||
use image_registration::error::Error;
|
||||
use image_registration::transform::Transform;
|
||||
|
||||
@@ -284,18 +300,24 @@ Build a `MattesMetric` and evaluate the (negative) mutual information and its gr
|
||||
- `FixedAt(points)` — an explicit list of sample positions (fully deterministic; recommended for reproducible results).
|
||||
|
||||
```rust
|
||||
use algos::ObjectiveFunction;
|
||||
use image_registration::bspline::BSpline;
|
||||
use algos::ObjectiveFunction; // direct dependency of image-registration; add it to your crate to name the trait
|
||||
use image_registration::bspline::{BSpline, BSplineTrait};
|
||||
use image_registration::error::Error;
|
||||
use image_registration::julia_image;
|
||||
use image_registration::metric::{MattesMetric, SamplingArg, Sigma};
|
||||
use image_registration::transform::Transform;
|
||||
|
||||
fn main() -> Result<(), Error> {
|
||||
// `fixed` and `moving` are Array2<f64>
|
||||
let fixed = ndarray::Array2::<f64>::zeros((64, 64));
|
||||
let moving = ndarray::Array2::<f64>::zeros((64, 64));
|
||||
// `fixed` and `moving` are Array2<f64>; moving is a rotated + shifted copy of fixed
|
||||
let center = [31.5f64, 31.5];
|
||||
let fixed = julia_image(&[64, 64], &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0], ¢er, &[-0.8, 0.156])
|
||||
.mapv(|i| i as f64);
|
||||
let moving = Transform::<ndarray::Ix2>::from_rotation(10.0f64.to_radians(), ¢er)
|
||||
.interpolate_par::<3, _, _>(fixed.view(), &Default::default())?;
|
||||
|
||||
let f = Sigma::Absolute(vec![2.0, 2.0]).smooth(fixed.view())?;
|
||||
let m = Sigma::Absolute(vec![2.0, 2.0]).smooth(moving.view())?;
|
||||
let sigma = Sigma::Absolute(vec![2.0, 2.0]);
|
||||
let f = sigma.smooth(fixed.view())?;
|
||||
let m = sigma.smooth(moving.view())?;
|
||||
|
||||
let metric = MattesMetric::new(
|
||||
BSpline::<0, _>::new(f.view()), // fixed image B-spline (nearest-neighbor)
|
||||
@@ -336,11 +358,11 @@ fn main() -> Result<(), Error> {
|
||||
|
||||
| Module | Contents |
|
||||
|----------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
|
||||
| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_around_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `add_translation`, `add_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. |
|
||||
| `image_registration::register` | `Registration`, `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD). |
|
||||
| `image_registration::metric` | `MattesMetric` (Mattes mutual information with analytic gradient), `SamplingArg` (Fixed/FixedAt/Random), `Sigma` (None/Absolute/Relative), `FixedMu`. |
|
||||
| `image_registration::bspline` | `BSpline<N, D>` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, Parzen kernel helpers. |
|
||||
| `image_registration::filter` | Gaussian smoothing and FFT helpers: `gaussian_smooth`, `gaussian_kernel`, `fft`, `ifft`, `fft_freq`, `fft_shift`. |
|
||||
| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_around_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `add_translation`, `add_scaling`, `with_rotation`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`, `interpolate_order`, `interpolate_par_order`), registration (`register`, `register_debug`, `register_affine`, `register_translation`, `register_translation_fft`), serialization (`to_file`, `from_file`, and serde), plus the free function `transform_point`. |
|
||||
| `image_registration::register` | `Registration` (`new_affine`, `new_translation`, `with_steps`, `with_optimizer`, …), `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD). |
|
||||
| `image_registration::metric` | `MattesMetric` (Mattes mutual information with analytic gradient; `evaluate`/`gradient` via the `algos::ObjectiveFunction` trait), `new_from_arrays`, `residual_fit`, `SamplingArg` (Fixed/FixedAt/Random), `Sigma` (None/Absolute/Relative), `FixedMu`. |
|
||||
| `image_registration::bspline` | `BSpline<N, D>` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, `DefaultValue` (fill value for out-of-range pixels), Parzen kernel helpers. |
|
||||
| `image_registration::filter` | Gaussian smoothing and FFT helpers: `gaussian_smooth`, `gaussian_kernel`, `gaussian`, `fft`, `ifft`, `fft_freq`, `fft_shift`. |
|
||||
| `image_registration::optimize` | Optimizers: `lbfgs_minimize`, `asgd_minimize`, `AsgdConfig`. |
|
||||
| `image_registration::par_indexed_iter` | Parallel iterators over array indices. |
|
||||
| `image_registration::error` | `Error` type. |
|
||||
|
||||
+136
-2
@@ -7,8 +7,13 @@ use ndarray::{
|
||||
};
|
||||
use num::traits::FloatConst;
|
||||
use num::{Complex, cast::AsPrimitive, complex::ComplexFloat};
|
||||
use rand::RngExt;
|
||||
use rand::distr::StandardUniform;
|
||||
use rayon::iter::ParallelIterator;
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::hash::Hash;
|
||||
use std::ops::{Deref, MulAssign};
|
||||
use std::sync::{Arc, Mutex};
|
||||
|
||||
fn in_shape(shape_f: &[f64], v: &[f64]) -> bool {
|
||||
shape_f
|
||||
@@ -17,6 +22,123 @@ fn in_shape(shape_f: &[f64], v: &[f64]) -> bool {
|
||||
.all(|(s, x)| (-0.5 <= *x) && (x <= s))
|
||||
}
|
||||
|
||||
// TODO: sample nearby points
|
||||
#[derive(Debug, Serialize, Deserialize)]
|
||||
pub struct DefaultValue {
|
||||
#[serde(skip, default = "DefaultValue::make_rng")]
|
||||
rng: Arc<Mutex<rand::prelude::StdRng>>,
|
||||
values: Vec<f64>,
|
||||
weights: Vec<f64>,
|
||||
}
|
||||
|
||||
impl Clone for DefaultValue {
|
||||
fn clone(&self) -> Self {
|
||||
Self {
|
||||
rng: Self::make_rng(),
|
||||
values: self.values.clone(),
|
||||
weights: self.weights.clone(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl Default for DefaultValue {
|
||||
fn default() -> Self {
|
||||
Self {
|
||||
rng: Self::make_rng(),
|
||||
values: vec![0.0],
|
||||
weights: vec![1.0],
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl PartialEq for DefaultValue {
|
||||
fn eq(&self, other: &Self) -> bool {
|
||||
(self.values == other.values) && (self.weights == other.weights)
|
||||
}
|
||||
}
|
||||
|
||||
impl DefaultValue {
|
||||
pub fn new(mut values: Vec<f64>, mut weights: Vec<f64>) -> Self {
|
||||
if values.is_empty() {
|
||||
values.push(f64::NAN);
|
||||
}
|
||||
if weights.len() != values.len() {
|
||||
weights.resize(values.len(), 1.0);
|
||||
}
|
||||
weights.iter_mut().for_each(|w| {
|
||||
if !w.is_finite() {
|
||||
*w = 1.0;
|
||||
} else if *w < 0.0 {
|
||||
*w = 0.0
|
||||
}
|
||||
});
|
||||
let n = weights.iter().sum::<f64>();
|
||||
weights.iter_mut().fold(0.0, |mut acc, w| {
|
||||
acc += *w / n;
|
||||
*w = acc;
|
||||
acc
|
||||
});
|
||||
Self {
|
||||
rng: Self::make_rng(),
|
||||
values,
|
||||
weights,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn values(&self) -> &[f64] {
|
||||
&self.values
|
||||
}
|
||||
|
||||
pub fn weights(&self) -> &[f64] {
|
||||
&self.weights
|
||||
}
|
||||
|
||||
pub fn sample(&self) -> f64 {
|
||||
if self.values.len() == 1 {
|
||||
self.values[0]
|
||||
} else {
|
||||
let i = self.rng.lock().unwrap().sample::<f64, _>(StandardUniform);
|
||||
let n = self
|
||||
.weights
|
||||
.iter()
|
||||
.position(|&w| i <= w)
|
||||
.unwrap_or(self.weights.len() - 1);
|
||||
self.values[n]
|
||||
}
|
||||
}
|
||||
|
||||
fn make_rng() -> Arc<Mutex<rand::prelude::StdRng>> {
|
||||
Arc::new(Mutex::new(rand::make_rng()))
|
||||
}
|
||||
}
|
||||
|
||||
impl Eq for DefaultValue {}
|
||||
|
||||
impl Hash for DefaultValue {
|
||||
fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
|
||||
self.values.iter().for_each(|v| v.to_bits().hash(state));
|
||||
self.weights.iter().for_each(|v| v.to_bits().hash(state));
|
||||
}
|
||||
}
|
||||
|
||||
impl From<f64> for DefaultValue {
|
||||
fn from(value: f64) -> Self {
|
||||
Self::new(vec![value], vec![1.0])
|
||||
}
|
||||
}
|
||||
|
||||
impl From<Vec<f64>> for DefaultValue {
|
||||
fn from(value: Vec<f64>) -> Self {
|
||||
Self::new(value, Vec::new())
|
||||
}
|
||||
}
|
||||
|
||||
impl From<&[f64]> for DefaultValue {
|
||||
fn from(value: &[f64]) -> Self {
|
||||
Self::new(value.to_vec(), Vec::new())
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct BSplineMem {
|
||||
ndim: usize,
|
||||
@@ -74,7 +196,11 @@ where
|
||||
D: Dimension,
|
||||
BSpline<N, D>: BSplineTrait<D>,
|
||||
{
|
||||
pub fn interpolate(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
|
||||
pub fn interpolate(
|
||||
&self,
|
||||
transform: &Transform<D>,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error> {
|
||||
let mut mem = self.get_mem();
|
||||
let shape_f = self.shape_f();
|
||||
let center = shape_f.iter().map(|i| i / 2.0).collect::<Vec<_>>();
|
||||
@@ -91,13 +217,19 @@ where
|
||||
};
|
||||
if in_shape(&shape_f, &v) {
|
||||
*x = self.evaluate_at_continuous_index(&v, &mut mem)?;
|
||||
} else {
|
||||
*x = default_value.sample();
|
||||
}
|
||||
}
|
||||
|
||||
Ok(new)
|
||||
}
|
||||
|
||||
pub fn interpolate_par(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
|
||||
pub fn interpolate_par(
|
||||
&self,
|
||||
transform: &Transform<D>,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error> {
|
||||
let mem = self.get_mem();
|
||||
let shape_f = self.shape_f();
|
||||
let center = shape_f.iter().map(|i| i / 2.0).collect::<Vec<_>>();
|
||||
@@ -114,6 +246,8 @@ where
|
||||
.all(|(s, x)| (-0.5 <= *x) && (x <= s))
|
||||
{
|
||||
*x = self.evaluate_at_continuous_index(&v, &mut mem)?;
|
||||
} else {
|
||||
*x = default_value.sample();
|
||||
}
|
||||
Ok(mem)
|
||||
},
|
||||
|
||||
@@ -19,4 +19,6 @@ pub enum Error {
|
||||
InitialGuessLength(usize, usize),
|
||||
#[error("interpolation order {0} is not implemented")]
|
||||
InterpolationNotImplemented(usize),
|
||||
#[error("matrix is singular")]
|
||||
SingularMatrix,
|
||||
}
|
||||
+7
-5
@@ -1,3 +1,6 @@
|
||||
//! Module containing gaussian smooth functions as well as helper functions to calculate ND fast
|
||||
//! fourier transforms
|
||||
|
||||
use crate::error::Error;
|
||||
use itertools::Itertools;
|
||||
use ndarray::{Array, Array1, ArrayD, AsArray, Axis, Dimension, RemoveAxis, concatenate, s};
|
||||
@@ -113,7 +116,6 @@ mod tests {
|
||||
use super::*;
|
||||
use crate::julia_image;
|
||||
use ndarray::{Ix1, array};
|
||||
use tiffwrite::IJTiffFile;
|
||||
|
||||
#[test]
|
||||
fn smooth() -> Result<(), Box<dyn std::error::Error>> {
|
||||
@@ -135,10 +137,10 @@ mod tests {
|
||||
&[29.5, 39.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let im_b = gaussian_smooth(im_a.mapv(|i| i as f64).view(), &[4.0, 4.0])?;
|
||||
let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
||||
t.save(im_a.mapv(|i| i as u32).view(), 0, 0, 0)?;
|
||||
t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
|
||||
let _im_b = gaussian_smooth(im_a.mapv(|i| i as f64).view(), &[4.0, 4.0])?;
|
||||
// let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
||||
// t.save(im_a.mapv(|i| i as u32).view(), 0, 0, 0)?;
|
||||
// t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
|
||||
Ok(())
|
||||
}
|
||||
|
||||
|
||||
+11
-10
@@ -1,3 +1,5 @@
|
||||
#![cfg_attr(docsrs, feature(doc_cfg))]
|
||||
#![doc = include_str!("../README.md")]
|
||||
use rayon::iter::ParallelIterator;
|
||||
pub mod bspline;
|
||||
pub mod error;
|
||||
@@ -48,7 +50,6 @@ mod tests {
|
||||
use super::*;
|
||||
use crate::transform::Transform;
|
||||
use tempfile::NamedTempFile;
|
||||
use tiffwrite::IJTiffFile;
|
||||
|
||||
#[test]
|
||||
fn serialization() -> Result<(), Error> {
|
||||
@@ -79,12 +80,12 @@ mod tests {
|
||||
.collect::<Vec<_>>()
|
||||
.try_into()
|
||||
.unwrap();
|
||||
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, ¢er);
|
||||
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &[0.0, 0.0]);
|
||||
let transform_j = transform.parameters.clone().try_into().unwrap();
|
||||
let c = [-0.8, 0.156];
|
||||
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
||||
let k = julia_image(&shape, &transform_k, ¢er, &c).mapv(|x| x as f64);
|
||||
let n = transform.interpolate::<3, _, _>(k.view())?;
|
||||
let n = transform.interpolate::<3, _, _>(k.view(), &Default::default())?;
|
||||
let sj = j.iter().sum::<f64>();
|
||||
// let sk = k.iter().sum::<f64>();
|
||||
let sn = n.iter().sum::<f64>();
|
||||
@@ -119,12 +120,12 @@ mod tests {
|
||||
.collect::<Vec<_>>()
|
||||
.try_into()
|
||||
.unwrap();
|
||||
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, ¢er);
|
||||
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &[0.0, 0.0]);
|
||||
let transform_j = transform.parameters.clone().try_into().unwrap();
|
||||
let c = [-0.8, 0.156];
|
||||
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
||||
let k = julia_image(&shape, &transform_k, ¢er, &c).mapv(|x| x as f64);
|
||||
let n = transform.interpolate_par::<3, _, _>(k.view())?;
|
||||
let n = transform.interpolate_par::<3, _, _>(k.view(), &Default::default())?;
|
||||
let sj = j.iter().sum::<f64>();
|
||||
let sn = n.iter().sum::<f64>();
|
||||
|
||||
@@ -148,22 +149,22 @@ mod tests {
|
||||
&[3999.5, 2999.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view())?;
|
||||
let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
|
||||
assert_ne!(image_b[[4000, 3000]], 0.0);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn julia() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let image_a = julia_image(
|
||||
let _image_a = julia_image(
|
||||
&[60, 80],
|
||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||
&[29.5, 39.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
||||
t.save(image_a.view(), 0, 0, 0)?;
|
||||
// let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
||||
// t.save(image_a.view(), 0, 0, 0)?;
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
+271
-2
@@ -9,8 +9,7 @@ use itertools::Itertools;
|
||||
use ndarray::{Array, Array1, Array2, Array3, AsArray, Axis, Dimension, s};
|
||||
use num::cast::AsPrimitive;
|
||||
use rand::RngExt;
|
||||
use rayon::iter::IntoParallelIterator;
|
||||
use rayon::iter::ParallelIterator;
|
||||
use rayon::iter::{IntoParallelIterator, ParallelIterator};
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::cell::RefCell;
|
||||
use std::ops::Deref;
|
||||
@@ -65,6 +64,22 @@ impl FixedMu {
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// combine the dparameters of the variable parameters into a complete
|
||||
/// parameter set, with no uncertainty (0) on the fixed parameters
|
||||
pub fn combine_dparameters(&self, dparameters: &[f64]) -> Vec<f64> {
|
||||
debug_assert_eq!(dparameters.len(), self.len_variable);
|
||||
let mut dparameters_iter = dparameters.iter();
|
||||
self.iter()
|
||||
.map(|i| {
|
||||
if i.is_some() {
|
||||
0.0
|
||||
} else {
|
||||
*dparameters_iter.next().unwrap()
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// extract mu from complete parameter set
|
||||
pub fn extract_variable(&self, mu: &[f64]) -> Vec<f64> {
|
||||
self.iter()
|
||||
@@ -204,6 +219,24 @@ fn image_jacobian(point: &[f64], center: &[f64], image_derivative: &[f64]) -> Ve
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Local Gauss-Newton statistics of the intensity fit
|
||||
/// `fixed(x) = moving(T(x)) + noise` at a parameter point.
|
||||
///
|
||||
/// `hessian` is the weighted normal matrix `JᵀWJ` and `variance` is the
|
||||
/// weighted residual variance `Σ w r² / (Σ w − q − 1)`, with `q` the number of
|
||||
/// variable parameters. The extra degree of freedom accounts for the weighted
|
||||
/// mean intensity offset that is removed from the residuals. The parameter
|
||||
/// covariance of the fit is `variance * hessian⁻¹`.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct ResidualFit {
|
||||
/// weighted normal matrix `JᵀWJ`
|
||||
pub hessian: Vec<Vec<f64>>,
|
||||
/// residual variance (mean squared residual, corrected for the dof)
|
||||
pub variance: f64,
|
||||
/// number of sample points that contributed
|
||||
pub n_samples: usize,
|
||||
}
|
||||
|
||||
pub struct MattesMetric<D>
|
||||
where
|
||||
D: Dimension,
|
||||
@@ -440,6 +473,93 @@ where
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Compute the Gauss-Newton normal matrix and the residual variance of the
|
||||
/// intensity fit at `point`, using the same sample points and edge
|
||||
/// weighting as [`Self::eval`].
|
||||
///
|
||||
/// Residuals are `fixed(x) - moving(T(x))`, so a well-registered image pair
|
||||
/// gives a small variance. Returns `None` when interpolation fails or when
|
||||
/// too few samples contribute to estimate the residual variance.
|
||||
pub fn residual_fit(&self, point: &[f64]) -> Option<ResidualFit> {
|
||||
let complete_mu = self.fixed_mu.combine(point);
|
||||
let idx = self.sampling.index();
|
||||
let shape = self.shape.as_slice();
|
||||
let shape_f = self.shape.iter().map(|s| s - 0.5).collect::<Vec<_>>();
|
||||
let center = self.center.as_slice();
|
||||
let e = 1.0 / self.edge;
|
||||
let n_params = point.len();
|
||||
let mut hessian = vec![vec![0.0; n_params]; n_params];
|
||||
let mut sum_w = 0.0;
|
||||
let mut sum_wr = 0.0;
|
||||
let mut sum_wr2 = 0.0;
|
||||
let mut n_samples = 0;
|
||||
let mut mem_f = self.fixed.get_mem();
|
||||
let mut mem_m = self.moving.get_mem();
|
||||
for u in idx {
|
||||
let v = transform_point(&u, center, &complete_mu);
|
||||
if !shape_f
|
||||
.iter()
|
||||
.zip_eq(v.iter())
|
||||
.all(|(s, x)| (-0.5 <= *x) && (x <= s))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
let f = self
|
||||
.fixed
|
||||
.evaluate_at_continuous_index(&u, &mut mem_f)
|
||||
.ok()?;
|
||||
let (m, d) = self
|
||||
.moving
|
||||
.evaluate_value_and_derivative_at_continuous_index(&v, &mut mem_m)
|
||||
.ok()?;
|
||||
// same smooth edge weighting as in `eval`
|
||||
let mut w = 1.0;
|
||||
for (s, x) in shape.iter().zip_eq(&v) {
|
||||
if e * (x + 0.5) < *s {
|
||||
let k = 3.0 * e * (x + 0.5) / s - 1.5;
|
||||
w *= square_bspline_integral(k);
|
||||
} else if e * (x + 0.5) > (e - 1.0) * s {
|
||||
let k = 3.0 * e * (s - x - 0.5) / s - 1.5;
|
||||
w *= square_bspline_integral(k);
|
||||
}
|
||||
}
|
||||
if w <= 0.0 || !w.is_finite() {
|
||||
continue;
|
||||
}
|
||||
let r = f - m;
|
||||
sum_w += w;
|
||||
sum_wr += w * r;
|
||||
sum_wr2 += w * r * r;
|
||||
n_samples += 1;
|
||||
let j = self
|
||||
.fixed_mu
|
||||
.extract_variable(&image_jacobian(&u, center, &d));
|
||||
for (a, ja) in j.iter().enumerate() {
|
||||
for (b, jb) in j.iter().enumerate() {
|
||||
hessian[a][b] += w * ja * jb;
|
||||
}
|
||||
}
|
||||
}
|
||||
// the residual variance needs at least one degree of freedom
|
||||
if n_samples == 0 || sum_w <= n_params as f64 + 1.0 {
|
||||
return None;
|
||||
}
|
||||
// subtract the weighted mean intensity mismatch (e.g. a constant
|
||||
// offset between the images) before estimating the variance
|
||||
let mean_r = sum_wr / sum_w;
|
||||
let sse = (sum_wr2 - sum_w * mean_r * mean_r).max(0.0);
|
||||
// q parameters plus the mean intensity offset that was removed
|
||||
let variance = sse / (sum_w - n_params as f64 - 1.0);
|
||||
if !variance.is_finite() {
|
||||
return None;
|
||||
}
|
||||
Some(ResidualFit {
|
||||
hessian,
|
||||
variance,
|
||||
n_samples,
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
impl<D> ObjectiveFunction<f64> for MattesMetric<D>
|
||||
@@ -455,6 +575,39 @@ where
|
||||
self.eval(point).unwrap();
|
||||
Some(self.metric.borrow().derivative.clone())
|
||||
}
|
||||
|
||||
/// The hessian is computed numerically: central finite differences of
|
||||
/// `gradient` with respect to each parameter (2 * n full evaluations).
|
||||
/// `MattesMetric` only accumulates first-order terms, so an analytic hessian
|
||||
/// is not available. The result is symmetrized, which also damps the small
|
||||
/// asymmetry introduced by finite differencing.
|
||||
fn hessian(&self, point: &[f64]) -> Option<Vec<Vec<f64>>> {
|
||||
let n = point.len();
|
||||
let mut hessian = vec![vec![0.0; n]; n];
|
||||
let mut p = point.to_vec();
|
||||
for j in 0..n {
|
||||
let h = 1e-5 * f64::max(1.0, point[j].abs());
|
||||
p[j] += h;
|
||||
self.eval(&p).unwrap();
|
||||
let g_plus = self.metric.borrow().derivative.clone();
|
||||
p[j] -= 2.0 * h;
|
||||
self.eval(&p).unwrap();
|
||||
let g_minus = self.metric.borrow().derivative.clone();
|
||||
p[j] += h;
|
||||
for i in 0..n {
|
||||
hessian[i][j] = (g_plus[i] - g_minus[i]) / (2.0 * h);
|
||||
}
|
||||
}
|
||||
#[expect(clippy::needless_range_loop)]
|
||||
for i in 0..n {
|
||||
for j in (i + 1)..n {
|
||||
let w = (hessian[i][j] + hessian[j][i]) / 2.0;
|
||||
hessian[i][j] = w;
|
||||
hessian[j][i] = w;
|
||||
}
|
||||
}
|
||||
Some(hessian)
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
|
||||
@@ -588,6 +741,122 @@ mod tests {
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn residual_fit_1d() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let a = array![1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0, 2.0, 0.0, 1.0];
|
||||
// b is a shifted right by 3 (padded on the left with the first value)
|
||||
let b = array![1.0, 1.0, 1.0, 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0];
|
||||
let points = (0..11).map(|i| vec![i as f64]).collect::<Vec<_>>();
|
||||
let m = MattesMetric::new_from_arrays(
|
||||
a.view(),
|
||||
b.view(),
|
||||
SamplingArg::FixedAt(points),
|
||||
16,
|
||||
0.05,
|
||||
)?
|
||||
.with_fixed_mu(vec![Some(1.0), None]);
|
||||
|
||||
// fixed(x) = moving(x + 3): no residual
|
||||
let fit = m.residual_fit(&[3.0]).expect("residual fit");
|
||||
assert!(fit.n_samples > 0);
|
||||
assert!(fit.hessian[0][0] > 0.0);
|
||||
assert!(fit.variance < 1e-12);
|
||||
|
||||
// away from the optimum the residual variance must grow
|
||||
let worse = m.residual_fit(&[0.0]).expect("residual fit");
|
||||
assert!(worse.variance > fit.variance);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hessian_1d() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let im_a = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0];
|
||||
let im_b = array![1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0];
|
||||
let a = gaussian_smooth(im_a.view(), &[2.0])?;
|
||||
let b = gaussian_smooth(im_b.view(), &[2.0])?;
|
||||
|
||||
let points = vec![
|
||||
vec![3.0],
|
||||
vec![4.0],
|
||||
vec![5.0],
|
||||
vec![6.0],
|
||||
vec![7.0],
|
||||
vec![8.0],
|
||||
];
|
||||
|
||||
let m = MattesMetric::new_from_arrays(
|
||||
a.view(),
|
||||
b.view(),
|
||||
SamplingArg::FixedAt(points),
|
||||
5,
|
||||
0.05,
|
||||
)?
|
||||
.with_fixed_mu(vec![Some(1.0), None]);
|
||||
|
||||
let p = -3.49;
|
||||
let h = m.hessian(&[p]).unwrap();
|
||||
assert_eq!(h.len(), 1);
|
||||
assert_eq!(h[0].len(), 1);
|
||||
|
||||
// cross-check against the second central difference of `evaluate`
|
||||
let step = 1e-4;
|
||||
let fd = (m.evaluate(&[p + step]) - 2.0 * m.evaluate(&[p]) + m.evaluate(&[p - step]))
|
||||
/ (step * step);
|
||||
println!("hessian = {:?}\nfinite diff = {:?}", h[0][0], fd);
|
||||
assert!((h[0][0] - fd).abs() < 1e-3 + 1e-2 * fd.abs());
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hessian_2d() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let a = Array2::from_shape_fn((12, 12), |(i, j)| i as f64 + 0.7 * j as f64);
|
||||
let b = Array2::from_shape_fn((12, 12), |(i, j)| 0.5 * i as f64 + 24.0 - 0.7 * j as f64);
|
||||
let a = gaussian_smooth(a.view(), &[2.0, 2.0])?;
|
||||
let b = gaussian_smooth(b.view(), &[2.0, 2.0])?;
|
||||
|
||||
let points = (0..12)
|
||||
.flat_map(|i| (0..12).map(move |j| vec![i as f64, j as f64]))
|
||||
.step_by(2)
|
||||
.collect::<Vec<_>>();
|
||||
|
||||
let m = MattesMetric::new_from_arrays(
|
||||
a.view(),
|
||||
b.view(),
|
||||
SamplingArg::FixedAt(points),
|
||||
7,
|
||||
0.05,
|
||||
)?;
|
||||
let n = m.fixed_mu().len_variable();
|
||||
assert_eq!(n, 2 * 2 + 2);
|
||||
|
||||
let p = vec![0.5, -0.3, 0.2, 0.7, -0.4, 0.3];
|
||||
let h = m.hessian(&p).unwrap();
|
||||
assert_eq!(h.len(), n);
|
||||
assert!(h.iter().all(|row| row.len() == n));
|
||||
// symmetrized by construction
|
||||
for (i, row) in h.iter().enumerate() {
|
||||
for (j, v) in row.iter().enumerate() {
|
||||
assert_eq!(*v, h[j][i]);
|
||||
}
|
||||
}
|
||||
|
||||
// cross-check each diagonal entry against the 2nd central difference of `evaluate`
|
||||
let step = 1e-4;
|
||||
for j in 0..n {
|
||||
let mut pp = p.clone();
|
||||
let mut pm = p.clone();
|
||||
pp[j] += step;
|
||||
pm[j] -= step;
|
||||
let fd = (m.evaluate(&pp) - 2.0 * m.evaluate(&p) + m.evaluate(&pm)) / (step * step);
|
||||
println!("diagonal[{j}] = {:?}\nfinite diff = {:?}", h[j][j], fd);
|
||||
assert!(
|
||||
(h[j][j] - fd).abs() < 1e-3 + 1e-2 * fd.abs(),
|
||||
"diagonal entry {j}"
|
||||
);
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_square_bspline_integral() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let x = Array1::linspace(-2.0, 2.0, 500);
|
||||
|
||||
+288
-11
@@ -1,10 +1,10 @@
|
||||
use crate::bspline::{BSpline, BSplineTrait};
|
||||
use crate::error::Error;
|
||||
use crate::metric::{FixedMu, MattesMetric, SamplingArg, Sigma};
|
||||
use crate::metric::{FixedMu, MattesMetric, ResidualFit, SamplingArg, Sigma};
|
||||
use crate::optimize::{AsgdConfig, asgd_minimize, lbfgs_minimize};
|
||||
use crate::transform::Transform;
|
||||
use crate::transform::{Transform, try_matrix_inverse};
|
||||
use algos::{ObjectiveFunction, OptimizationConfig};
|
||||
use ndarray::{AsArray, Dimension};
|
||||
use ndarray::{Array2, AsArray, Dimension};
|
||||
use num::cast::AsPrimitive;
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::marker::PhantomData;
|
||||
@@ -175,7 +175,7 @@ impl<D: Dimension> Registration<D> {
|
||||
} else {
|
||||
return Err(Error::NumberOfDimensionsNotDefined);
|
||||
};
|
||||
let mut fixed_mu = vec![Some(0.0); ndim * ndim];
|
||||
let mut fixed_mu = vec![Some(0.0); ndim * ndim + ndim];
|
||||
for i in 0..ndim {
|
||||
fixed_mu[i * (ndim + 1)] = Some(1.0);
|
||||
fixed_mu[ndim * ndim + i] = None;
|
||||
@@ -286,8 +286,11 @@ impl<D: Dimension> Registration<D> {
|
||||
.unwrap_or_else(|| RegistrationStep::default_steps(ndim, n));
|
||||
|
||||
let mut p_var = self.initial_guess.as_ref().cloned().unwrap_or_else(|| {
|
||||
Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim)
|
||||
let full = Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim);
|
||||
self.fixed_mu.extract_variable(&full)
|
||||
});
|
||||
let mut last_metric: Option<MattesMetric<D>> = None;
|
||||
let mut last_result: Option<algos::OptimizationResult<f64>> = None;
|
||||
|
||||
for RegistrationStep {
|
||||
sigma,
|
||||
@@ -324,13 +327,23 @@ impl<D: Dimension> Registration<D> {
|
||||
&scales,
|
||||
);
|
||||
|
||||
if let Some(result) = optimization_result
|
||||
last_metric = Some(metric);
|
||||
last_result = optimization_result;
|
||||
if let Some(result) = last_result.as_ref()
|
||||
&& result.optimal_point.iter().all(|i| i.is_finite())
|
||||
{
|
||||
p_var = result.optimal_point;
|
||||
p_var = result.optimal_point.clone();
|
||||
}
|
||||
}
|
||||
Ok(Transform::<D>::new(self.fixed_mu.combine(p_var.as_slice())))
|
||||
let mut transform = Transform::<D>::new(self.fixed_mu.combine(p_var.as_slice()));
|
||||
if let (Some(metric), Some(result)) = (&last_metric, last_result.as_ref())
|
||||
&& result.converged
|
||||
&& let Some(fit) = metric.residual_fit(&p_var)
|
||||
{
|
||||
let dparameters = dparameters_from_residual_fit(&fit);
|
||||
transform = transform.with_dparameters(self.fixed_mu.combine_dparameters(&dparameters));
|
||||
}
|
||||
Ok(transform)
|
||||
}
|
||||
|
||||
fn optimize_metric<M: ObjectiveFunction<f64>>(
|
||||
@@ -389,9 +402,12 @@ impl<D: Dimension> Registration<D> {
|
||||
.unwrap_or_else(|| RegistrationStep::default_steps(ndim, n));
|
||||
|
||||
let mut p_var = self.initial_guess.as_ref().cloned().unwrap_or_else(|| {
|
||||
Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim)
|
||||
let full = Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim);
|
||||
self.fixed_mu.extract_variable(&full)
|
||||
});
|
||||
let mut p = vec![0.0; n * n + n];
|
||||
let mut p = vec![0.0; ndim * ndim + ndim];
|
||||
let mut last_metric: Option<MattesMetric<D>> = None;
|
||||
let mut last_result: Option<algos::OptimizationResult<f64>> = None;
|
||||
|
||||
let mut registration_results = Vec::new();
|
||||
for RegistrationStep {
|
||||
@@ -461,6 +477,8 @@ impl<D: Dimension> Registration<D> {
|
||||
.fixed_mu()
|
||||
.combine(&optimization_result.optimal_point);
|
||||
}
|
||||
last_metric = Some(metric);
|
||||
last_result = Some(optimization_result.clone());
|
||||
|
||||
registration_results.push(RegistrationResult {
|
||||
sigma_fixed: sigma.sigma(f_orig_shape.as_slice()),
|
||||
@@ -477,6 +495,265 @@ impl<D: Dimension> Registration<D> {
|
||||
converged: optimization_result.converged,
|
||||
});
|
||||
}
|
||||
Ok((Transform::<D>::new(p), registration_results))
|
||||
let mut transform = Transform::<D>::new(p);
|
||||
if let (Some(metric), Some(result)) = (&last_metric, last_result.as_ref())
|
||||
&& result.converged
|
||||
&& let Some(fit) = metric.residual_fit(&p_var)
|
||||
{
|
||||
let dparameters = dparameters_from_residual_fit(&fit);
|
||||
transform = transform.with_dparameters(self.fixed_mu.combine_dparameters(&dparameters));
|
||||
}
|
||||
Ok((transform, registration_results))
|
||||
}
|
||||
}
|
||||
|
||||
/// Diagonal of the inverse of a symmetric positive-definite matrix.
|
||||
///
|
||||
/// Non-finite entries are treated as zero. Returns `None` if the matrix is
|
||||
/// singular, because a singular normal matrix means the corresponding
|
||||
/// directions are not determined by the data. Diagonal entries of the inverse
|
||||
/// that are not finite and positive are reported as 0.
|
||||
fn variance_from_hessian(hessian: &[Vec<f64>]) -> Option<Vec<f64>> {
|
||||
let n = hessian.len();
|
||||
if hessian.iter().any(|row| row.len() != n) {
|
||||
return None;
|
||||
}
|
||||
let mut data: Vec<f64> = Vec::with_capacity(n * n);
|
||||
data.extend(
|
||||
hessian
|
||||
.iter()
|
||||
.flatten()
|
||||
.map(|v| if v.is_finite() { *v } else { 0.0 }),
|
||||
);
|
||||
let matrix = Array2::from_shape_vec((n, n), data).ok()?;
|
||||
let inv = try_matrix_inverse(matrix)?;
|
||||
Some(
|
||||
(0..n)
|
||||
.map(|i| {
|
||||
let v = inv[[i, i]];
|
||||
if v.is_finite() && v > 0.0 { v } else { 0.0 }
|
||||
})
|
||||
.collect(),
|
||||
)
|
||||
}
|
||||
|
||||
/// Compute parameter uncertainty (standard error, the `dparameters` field) from
|
||||
/// the intensity fit at the optimum.
|
||||
///
|
||||
/// The moving image is modelled as `moving(T(x)) = fixed(x) + noise` with
|
||||
/// independent residuals. Around the optimum this is a linear least-squares
|
||||
/// problem with Jacobian `J`, so the parameter covariance is
|
||||
///
|
||||
/// ```text
|
||||
/// Cov = σ² (JᵀWJ)⁻¹
|
||||
/// ```
|
||||
///
|
||||
/// with `σ²` the residual variance and `JᵀWJ` the weighted Gauss-Newton normal
|
||||
/// matrix, both computed by [`MattesMetric::residual_fit`]. This is the usual
|
||||
/// standard error of a nonlinear least-squares fit: it is calibrated in
|
||||
/// parameter units and grows when the images fit poorly or when the image
|
||||
/// gradients are small.
|
||||
///
|
||||
/// This deliberately does *not* use the inverse Hessian of `-MI`: mutual
|
||||
/// information is not a likelihood, so its curvature has an arbitrary scale
|
||||
/// (it changes with the number of histogram bins and the Parzen width) and
|
||||
/// `H⁻¹` is not a standard error.
|
||||
///
|
||||
/// Each uncertainty is the square root of the corresponding diagonal. If the
|
||||
/// normal matrix is singular or the residuals cannot be estimated, the result
|
||||
/// is 0 for the affected parameters, meaning "not estimated".
|
||||
pub(crate) fn dparameters_from_residual_fit(fit: &ResidualFit) -> Vec<f64> {
|
||||
let n = fit.hessian.len();
|
||||
if !fit.variance.is_finite() || fit.variance <= 0.0 {
|
||||
return vec![0.0; n];
|
||||
}
|
||||
variance_from_hessian(&fit.hessian)
|
||||
.map(|variance| {
|
||||
variance
|
||||
.into_iter()
|
||||
.map(|v| {
|
||||
let v = fit.variance * v;
|
||||
if v.is_finite() && v > 0.0 {
|
||||
v.sqrt()
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
})
|
||||
.unwrap_or_else(|| vec![0.0; n])
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use ndarray::{Array1, Ix1};
|
||||
|
||||
#[test]
|
||||
fn dparameters_from_residual_fit_1d() {
|
||||
// Cov = variance / (JᵀWJ), so with JᵀWJ = 4 and variance = 1/4 the
|
||||
// standard error is sqrt(1/16) = 1/4.
|
||||
let fit = ResidualFit {
|
||||
hessian: vec![vec![4.0]],
|
||||
variance: 0.25,
|
||||
n_samples: 10,
|
||||
};
|
||||
let d = dparameters_from_residual_fit(&fit);
|
||||
assert!((d[0] - 0.25).abs() < 1e-12);
|
||||
|
||||
// the covariance scales linearly with the residual variance
|
||||
let fit = ResidualFit {
|
||||
variance: 1.0,
|
||||
..fit
|
||||
};
|
||||
let d = dparameters_from_residual_fit(&fit);
|
||||
assert!((d[0] - 0.5).abs() < 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn dparameters_from_residual_fit_2d() {
|
||||
// inverse of [[2, 1], [1, 2]] has diagonal 2/3
|
||||
let fit = ResidualFit {
|
||||
hessian: vec![vec![2.0, 1.0], vec![1.0, 2.0]],
|
||||
variance: 0.5,
|
||||
n_samples: 20,
|
||||
};
|
||||
let d = dparameters_from_residual_fit(&fit);
|
||||
let expected: f64 = (1.0f64 / 3.0).sqrt();
|
||||
assert!((d[0] - expected).abs() < 1e-9);
|
||||
assert!((d[1] - expected).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn dparameters_from_residual_fit_not_estimated() {
|
||||
// a singular normal matrix leaves the parameters undetermined; there
|
||||
// is no diagonal fallback
|
||||
let singular = ResidualFit {
|
||||
hessian: vec![vec![2.0, -2.0], vec![-2.0, 2.0]],
|
||||
variance: 1.0,
|
||||
n_samples: 20,
|
||||
};
|
||||
assert_eq!(dparameters_from_residual_fit(&singular), vec![0.0, 0.0]);
|
||||
|
||||
// a non-positive-definite or non-finite normal matrix is not
|
||||
// informative either
|
||||
for hessian in [vec![vec![0.0]], vec![vec![-4.0]], vec![vec![f64::NAN]]] {
|
||||
let fit = ResidualFit {
|
||||
hessian,
|
||||
variance: 1.0,
|
||||
n_samples: 10,
|
||||
};
|
||||
assert_eq!(dparameters_from_residual_fit(&fit), vec![0.0]);
|
||||
}
|
||||
|
||||
// no residual variance means no standard error is available
|
||||
for variance in [0.0, -1.0, f64::NAN, f64::INFINITY] {
|
||||
let fit = ResidualFit {
|
||||
hessian: vec![vec![4.0]],
|
||||
variance,
|
||||
n_samples: 10,
|
||||
};
|
||||
assert_eq!(dparameters_from_residual_fit(&fit), vec![0.0]);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn combine_dparameters() {
|
||||
let fixed_mu: FixedMu = vec![Some(1.0), None, None].into();
|
||||
assert_eq!(
|
||||
fixed_mu.combine_dparameters(&[0.1, 0.2]),
|
||||
vec![0.0, 0.1, 0.2]
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn register_deterministic_dparameters() -> Result<(), Box<dyn std::error::Error>> {
|
||||
// dense, non-periodic 1D signal; b is a, shifted right by 2: the
|
||||
// optimal translation is +2
|
||||
// incommensurate frequencies keep the signal effectively aperiodic,
|
||||
// so the mutual information peaks at the true shift
|
||||
let n = 256;
|
||||
let a: Array1<f64> = (0..n)
|
||||
.map(|i| {
|
||||
let x = i as f64;
|
||||
1.5 + (0.131 * x).sin()
|
||||
+ 0.7 * (0.397 * x + 1.3).sin()
|
||||
+ 0.3 * (0.883 * x + 0.7).sin()
|
||||
})
|
||||
.collect();
|
||||
let b: Array1<f64> = (0..n)
|
||||
.map(|i| if i < 2 { 0.0 } else { a[[i - 2]] })
|
||||
.collect();
|
||||
let points: Vec<Vec<f64>> = (0..n).map(|i| vec![i as f64]).collect();
|
||||
let steps = vec![
|
||||
RegistrationStep::new(
|
||||
Sigma::Absolute(vec![2.0]),
|
||||
SamplingArg::FixedAt(points.clone()),
|
||||
24,
|
||||
1e-6,
|
||||
0.05,
|
||||
1024,
|
||||
1.0,
|
||||
Optimizer::LBFGS,
|
||||
),
|
||||
RegistrationStep::new(
|
||||
Sigma::None,
|
||||
SamplingArg::FixedAt(points),
|
||||
48,
|
||||
1e-6,
|
||||
0.05,
|
||||
2048,
|
||||
1.0,
|
||||
Optimizer::LBFGS,
|
||||
),
|
||||
];
|
||||
let registration = Registration::<Ix1>::new_translation()?.with_steps(steps);
|
||||
let transform = registration.register(a.view(), b.view())?;
|
||||
assert_eq!(transform.parameters.len(), 2);
|
||||
assert_eq!(transform.dparameters.len(), 2);
|
||||
// the fixed scale parameter has no uncertainty
|
||||
assert_eq!(transform.dparameters[0], 0.0);
|
||||
// the data are a perfect shifted copy, so the residual variance is
|
||||
// (numerically) zero and so is the translation uncertainty
|
||||
assert!(transform.dparameters[1].is_finite());
|
||||
assert!(transform.dparameters[1] < 0.05);
|
||||
// the found translation should be close to the applied one
|
||||
assert!((transform.parameters[1] - 2.0).abs() < 1.0);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noisy_fit_has_positive_dparameters() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let n = 128;
|
||||
let a: Array1<f64> = (0..n)
|
||||
.map(|i| {
|
||||
let x = i as f64;
|
||||
1.5 + (0.131 * x).sin()
|
||||
+ 0.7 * (0.397 * x + 1.3).sin()
|
||||
+ 0.3 * (0.883 * x + 0.7).sin()
|
||||
})
|
||||
.collect();
|
||||
let b: Array1<f64> = (0..n)
|
||||
.map(|i| {
|
||||
let base = if i < 2 { 0.0 } else { a[[i - 2]] };
|
||||
// deterministic broadband perturbation
|
||||
base + 0.05 * ((i * 37 % 19) as f64 / 18.0 - 0.5)
|
||||
})
|
||||
.collect();
|
||||
let points: Vec<Vec<f64>> = (0..n).map(|i| vec![i as f64]).collect();
|
||||
let metric = MattesMetric::new_from_arrays(
|
||||
a.view(),
|
||||
b.view(),
|
||||
SamplingArg::FixedAt(points),
|
||||
16,
|
||||
0.05,
|
||||
)?
|
||||
.with_fixed_mu(vec![Some(1.0), None]);
|
||||
|
||||
let fit = metric.residual_fit(&[2.0]).expect("residual fit");
|
||||
assert!(fit.variance > 0.0);
|
||||
let d = dparameters_from_residual_fit(&fit);
|
||||
assert!(d[0] > 0.0);
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
+129
-69
@@ -1,10 +1,14 @@
|
||||
use crate::bspline::{BSpline, BSplineTrait};
|
||||
use crate::bspline::{BSpline, BSplineTrait, DefaultValue};
|
||||
use crate::error::Error;
|
||||
use crate::filter::{fft, ifft};
|
||||
use crate::metric::FixedMu;
|
||||
use crate::register::{Registration, RegistrationResult, RegistrationStep};
|
||||
use itertools::Itertools;
|
||||
use ndarray::{Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, Ix2, IxDyn, s};
|
||||
use ndarray::{
|
||||
Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, IntoDimension, Ix2, IxDyn, s,
|
||||
};
|
||||
use num::cast::AsPrimitive;
|
||||
use num::complex::ComplexFloat;
|
||||
use serde::{Deserialize, Serialize};
|
||||
use serde_yaml::{from_reader, to_writer};
|
||||
use std::fs::File;
|
||||
@@ -44,11 +48,18 @@ where
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Calculate the inverse of a square matrix using LU factorization with partial pivoting.
|
||||
fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
|
||||
/// Calculate the inverse of a square matrix using LU factorization with partial
|
||||
/// pivoting.
|
||||
///
|
||||
/// Returns `None` if the matrix is singular, or numerically rank-deficient.
|
||||
pub(crate) fn try_matrix_inverse(matrix: Array2<f64>) -> Option<Array2<f64>> {
|
||||
let n = matrix.nrows();
|
||||
debug_assert_eq!(n, matrix.ncols(), "matrix must be square");
|
||||
|
||||
// pivot magnitudes must exceed this relative threshold to count as significant
|
||||
let max_abs = matrix.iter().fold(0.0f64, |m, &x| m.max(x.abs()));
|
||||
let threshold = max_abs * (n as f64) * f64::EPSILON;
|
||||
|
||||
// PA = LU factorization (Doolittle with partial pivoting)
|
||||
let mut lu = matrix;
|
||||
let mut perm: Vec<usize> = (0..n).collect();
|
||||
@@ -60,7 +71,9 @@ fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
|
||||
pivot_row = i;
|
||||
}
|
||||
}
|
||||
assert_ne!(lu[[pivot_row, k]], 0.0, "matrix is singular");
|
||||
if lu[[pivot_row, k]].abs() <= threshold || !lu[[pivot_row, k]].is_finite() {
|
||||
return None;
|
||||
}
|
||||
if pivot_row != k {
|
||||
for j in 0..n {
|
||||
lu.swap([k, j], [pivot_row, j]);
|
||||
@@ -98,7 +111,7 @@ fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
|
||||
}
|
||||
}
|
||||
|
||||
inv
|
||||
Some(inv)
|
||||
}
|
||||
|
||||
/// a struct describing the transform
|
||||
@@ -153,10 +166,10 @@ impl<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
|
||||
.cloned()
|
||||
.collect(),
|
||||
dparameters: dm
|
||||
.slice(s![..self.ndim(), ..])
|
||||
.slice(s![..self.ndim(), ..self.ndim()])
|
||||
.flatten()
|
||||
.iter()
|
||||
.chain(dm.slice(s![..self.ndim(), ..]).iter())
|
||||
.chain(dm.slice(s![..self.ndim(), self.ndim()]).iter())
|
||||
.cloned()
|
||||
.collect(),
|
||||
dimension: self.dimension,
|
||||
@@ -399,8 +412,8 @@ impl<D: Dimension> Transform<D> {
|
||||
}
|
||||
|
||||
pub fn new_around_center(parameters: Vec<f64>, center: &[f64]) -> Self {
|
||||
(Self::new(vec![1.0, 0.0, 0.0, 1.0, -center[0], -center[1]]) * Self::new(parameters))
|
||||
* Self::new(vec![1.0, 0.0, 0.0, 1.0, center[0], center[1]])
|
||||
Self::new(parameters).with_translation(center)
|
||||
* Self::from_translation(¢er.iter().map(|i| -i).collect::<Vec<_>>())
|
||||
}
|
||||
|
||||
/// find the affine transform which transforms moving into fixed
|
||||
@@ -420,7 +433,46 @@ impl<D: Dimension> Transform<D> {
|
||||
M: AsArray<'a, T, D>,
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
{
|
||||
Registration::new_translation()?.register(fixed, moving)
|
||||
let fixed = fixed.into();
|
||||
let moving = moving.into();
|
||||
let initial_guess = Self::register_affine(&fixed, &moving)?;
|
||||
let initial_guess =
|
||||
initial_guess.parameters[initial_guess.ndim() * initial_guess.ndim()..].to_vec();
|
||||
Registration::new_translation()?
|
||||
.with_initial_guess(initial_guess)?
|
||||
.register(fixed, moving)
|
||||
}
|
||||
|
||||
/// find the translation which transforms moving into fixed
|
||||
/// using the cross-correlation theorem and fft
|
||||
pub fn register_translation_fft<'a, F, M, T>(fixed: F, moving: M) -> Result<Self, Error>
|
||||
where
|
||||
F: AsArray<'a, T, D>,
|
||||
M: AsArray<'a, T, D>,
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
{
|
||||
let f = fft(fixed.into().mapv(|i| i.as_()).view())?.mapv(|i| i.conj());
|
||||
let m = fft(moving.into().mapv(|i| i.as_()).view())?;
|
||||
let c = ifft(&(f * m))?.mapv(|i| i.re);
|
||||
let (i, _) = c
|
||||
.indexed_iter()
|
||||
.max_by(|(_, a), (_, b)| a.total_cmp(b))
|
||||
.expect("no max because array is empty");
|
||||
let d = i.into_dimension();
|
||||
let j = d
|
||||
.as_array_view()
|
||||
.iter()
|
||||
.zip(c.shape())
|
||||
.map(|(&i, &s)| {
|
||||
let mut fi = i as f64;
|
||||
let fs = s as f64;
|
||||
if 2.0 * fi > fs {
|
||||
fi -= fs;
|
||||
}
|
||||
fi
|
||||
})
|
||||
.collect::<Vec<_>>();
|
||||
Ok(Self::from_translation(&j))
|
||||
}
|
||||
|
||||
/// find the transform which transforms moving into fixed, using fixed_mu to specify which
|
||||
@@ -753,7 +805,7 @@ impl<D: Dimension> Transform<D> {
|
||||
/// get the inverse transform
|
||||
pub fn inverse(&self) -> Result<Self, Error> {
|
||||
let matrix = self.matrix();
|
||||
let inverse = matrix_inverse(matrix);
|
||||
let inverse = try_matrix_inverse(matrix).ok_or(Error::SingularMatrix)?;
|
||||
let parameters = inverse
|
||||
.slice(s![..self.ndim(), ..self.ndim()])
|
||||
.iter()
|
||||
@@ -768,21 +820,10 @@ impl<D: Dimension> Transform<D> {
|
||||
})
|
||||
}
|
||||
|
||||
pub fn interpolate<'a, const B: usize, T, A>(&self, image: A) -> Result<Array<f64, D>, Error>
|
||||
where
|
||||
A: AsArray<'a, T, D>,
|
||||
D: Dimension,
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
BSpline<B, D>: BSplineTrait<D>,
|
||||
{
|
||||
let image = image.into();
|
||||
let bspline = BSpline::<B, _>::new(&image);
|
||||
bspline.interpolate(self)
|
||||
}
|
||||
|
||||
pub fn interpolate_par<'a, const B: usize, T, A>(
|
||||
pub fn interpolate<'a, const B: usize, T, A>(
|
||||
&self,
|
||||
image: A,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error>
|
||||
where
|
||||
A: AsArray<'a, T, D>,
|
||||
@@ -792,13 +833,30 @@ impl<D: Dimension> Transform<D> {
|
||||
{
|
||||
let image = image.into();
|
||||
let bspline = BSpline::<B, _>::new(&image);
|
||||
bspline.interpolate_par(self)
|
||||
bspline.interpolate(self, default_value)
|
||||
}
|
||||
|
||||
pub fn interpolate_par<'a, const B: usize, T, A>(
|
||||
&self,
|
||||
image: A,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error>
|
||||
where
|
||||
A: AsArray<'a, T, D>,
|
||||
D: Dimension,
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
BSpline<B, D>: BSplineTrait<D>,
|
||||
{
|
||||
let image = image.into();
|
||||
let bspline = BSpline::<B, _>::new(&image);
|
||||
bspline.interpolate_par(self, default_value)
|
||||
}
|
||||
|
||||
pub fn interpolate_order<'a, T, A>(
|
||||
&self,
|
||||
image: A,
|
||||
order: usize,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error>
|
||||
where
|
||||
A: AsArray<'a, T, D>,
|
||||
@@ -806,12 +864,12 @@ impl<D: Dimension> Transform<D> {
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
{
|
||||
match order {
|
||||
0 => Ok(self.interpolate::<0, _, _>(image)?),
|
||||
1 => Ok(self.interpolate::<1, _, _>(image)?),
|
||||
2 => Ok(self.interpolate::<2, _, _>(image)?),
|
||||
3 => Ok(self.interpolate::<3, _, _>(image)?),
|
||||
4 => Ok(self.interpolate::<4, _, _>(image)?),
|
||||
5 => Ok(self.interpolate::<5, _, _>(image)?),
|
||||
0 => Ok(self.interpolate::<0, _, _>(image, default_value)?),
|
||||
1 => Ok(self.interpolate::<1, _, _>(image, default_value)?),
|
||||
2 => Ok(self.interpolate::<2, _, _>(image, default_value)?),
|
||||
3 => Ok(self.interpolate::<3, _, _>(image, default_value)?),
|
||||
4 => Ok(self.interpolate::<4, _, _>(image, default_value)?),
|
||||
5 => Ok(self.interpolate::<5, _, _>(image, default_value)?),
|
||||
_ => Err(Error::InterpolationNotImplemented(order)),
|
||||
}
|
||||
}
|
||||
@@ -820,6 +878,7 @@ impl<D: Dimension> Transform<D> {
|
||||
&self,
|
||||
image: A,
|
||||
order: usize,
|
||||
default_value: &DefaultValue,
|
||||
) -> Result<Array<f64, D>, Error>
|
||||
where
|
||||
A: AsArray<'a, T, D>,
|
||||
@@ -827,12 +886,12 @@ impl<D: Dimension> Transform<D> {
|
||||
T: 'a + Clone + AsPrimitive<f64>,
|
||||
{
|
||||
match order {
|
||||
0 => Ok(self.interpolate_par::<0, _, _>(image)?),
|
||||
1 => Ok(self.interpolate_par::<1, _, _>(image)?),
|
||||
2 => Ok(self.interpolate_par::<2, _, _>(image)?),
|
||||
3 => Ok(self.interpolate_par::<3, _, _>(image)?),
|
||||
4 => Ok(self.interpolate_par::<4, _, _>(image)?),
|
||||
5 => Ok(self.interpolate_par::<5, _, _>(image)?),
|
||||
0 => Ok(self.interpolate_par::<0, _, _>(image, default_value)?),
|
||||
1 => Ok(self.interpolate_par::<1, _, _>(image, default_value)?),
|
||||
2 => Ok(self.interpolate_par::<2, _, _>(image, default_value)?),
|
||||
3 => Ok(self.interpolate_par::<3, _, _>(image, default_value)?),
|
||||
4 => Ok(self.interpolate_par::<4, _, _>(image, default_value)?),
|
||||
5 => Ok(self.interpolate_par::<5, _, _>(image, default_value)?),
|
||||
_ => Err(Error::InterpolationNotImplemented(order)),
|
||||
}
|
||||
}
|
||||
@@ -948,7 +1007,7 @@ impl Transform<Ix2> {
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::matrix_inverse;
|
||||
use super::try_matrix_inverse;
|
||||
use crate::julia_image;
|
||||
use crate::transform::File;
|
||||
use crate::transform::Transform;
|
||||
@@ -965,7 +1024,7 @@ mod tests {
|
||||
|
||||
let a = array![[4.0, 7.0], [2.0, 6.0]];
|
||||
let expected = array![[0.6, -0.7], [-0.2, 0.4]];
|
||||
let got = matrix_inverse(a.clone());
|
||||
let got = try_matrix_inverse(a.clone()).expect("matrix should not be singular");
|
||||
assert!(
|
||||
got.iter()
|
||||
.zip_eq(expected.iter())
|
||||
@@ -974,7 +1033,7 @@ mod tests {
|
||||
|
||||
let b = array![[0.0, 2.0], [3.0, 4.0]];
|
||||
let expected = array![[-2.0 / 3.0, 1.0 / 3.0], [0.5, 0.0]];
|
||||
let got = matrix_inverse(b.clone());
|
||||
let got = try_matrix_inverse(b.clone()).expect("matrix should not be singular");
|
||||
assert!(
|
||||
got.iter()
|
||||
.zip_eq(expected.iter())
|
||||
@@ -983,13 +1042,16 @@ mod tests {
|
||||
|
||||
let c = array![[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]];
|
||||
let reference = array![[-24.0, 18.0, 5.0], [20.0, -15.0, -4.0], [-5.0, 4.0, 1.0]];
|
||||
let got = matrix_inverse(c);
|
||||
let got = try_matrix_inverse(c).expect("matrix should not be singular");
|
||||
assert!(
|
||||
got.iter()
|
||||
.zip_eq(reference.iter())
|
||||
.all(|(x, y)| (x - y).abs() < 1e-10)
|
||||
);
|
||||
|
||||
let singular = array![[1.0, 2.0], [2.0, 4.0]];
|
||||
assert!(try_matrix_inverse(singular).is_none());
|
||||
|
||||
Ok(())
|
||||
}
|
||||
|
||||
@@ -1001,9 +1063,9 @@ mod tests {
|
||||
&[3999.5, 2999.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view())?;
|
||||
let image_c = transform.interpolate::<3, _, _>(image_a.view())?;
|
||||
let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
|
||||
let image_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
||||
|
||||
assert!(
|
||||
image_b
|
||||
@@ -1022,8 +1084,8 @@ mod tests {
|
||||
&[3999.5, 2999.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
||||
let image_b = transform.interpolate::<3, _, _>(image_a.view())?;
|
||||
let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
|
||||
let image_b = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
||||
assert_ne!(image_b[[4000, 3000]], 0.0);
|
||||
Ok(())
|
||||
}
|
||||
@@ -1036,8 +1098,8 @@ mod tests {
|
||||
&[3999.5, 2999.5],
|
||||
&[-0.8, 0.156],
|
||||
);
|
||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view())?;
|
||||
let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
|
||||
let image_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
|
||||
assert_ne!(image_b[[4000, 3000]], 0.0);
|
||||
Ok(())
|
||||
}
|
||||
@@ -1056,7 +1118,7 @@ mod tests {
|
||||
.slice(s![.., 0])
|
||||
.mapv(|i| i as f64);
|
||||
let q = vec![0.85, 4.0];
|
||||
let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a)?;
|
||||
let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a, &Default::default())?;
|
||||
|
||||
// Test with ALL points (like grad_check uses)
|
||||
let metric_all = MattesMetric::<ndarray::Ix1>::new(
|
||||
@@ -1140,7 +1202,8 @@ mod tests {
|
||||
)
|
||||
.slice(s![.., 0])
|
||||
.mapv(|i| i as f64);
|
||||
let im_b = Transform::new(vec![0.85, 4.0]).interpolate::<1, _, _>(&im_a)?;
|
||||
let im_b =
|
||||
Transform::new(vec![0.85, 4.0]).interpolate::<1, _, _>(&im_a, &Default::default())?;
|
||||
|
||||
// Smooth both images with sigma=8 like level 0
|
||||
let sigma = Sigma::Absolute(vec![8.0]);
|
||||
@@ -1249,7 +1312,8 @@ mod tests {
|
||||
];
|
||||
|
||||
let q = vec![0.85, 2.0];
|
||||
let moving = Transform::new(q.clone()).interpolate::<3, _, _>(&fixed)?;
|
||||
let moving =
|
||||
Transform::new(q.clone()).interpolate::<3, _, _>(&fixed, &Default::default())?;
|
||||
|
||||
// The registration finds T such that im_b(T(x)) = im_a(x), which is the inverse of q
|
||||
let q_inv = Transform::<ndarray::Ix1>::new(q.clone()).inverse()?;
|
||||
@@ -1305,7 +1369,7 @@ mod tests {
|
||||
.slice(s![.., 0])
|
||||
.mapv(|i| i as f64);
|
||||
let q = vec![0.85, 2.0];
|
||||
let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a)?;
|
||||
let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a, &Default::default())?;
|
||||
|
||||
let f = Sigma::Absolute(vec![1.0]).smooth(im_a.view())?;
|
||||
let m = Sigma::Absolute(vec![1.0]).smooth(im_b.view())?;
|
||||
@@ -1425,7 +1489,7 @@ mod tests {
|
||||
let params = vec![c * sx, -s * sy + shx, s * sx + shy, c * sy, tx, ty];
|
||||
|
||||
let transform = Transform::<Ix2>::new_around_center(params.clone(), ¢er);
|
||||
let moving = transform.interpolate::<3, _, _>(fixed.view())?;
|
||||
let moving = transform.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
|
||||
|
||||
let q_inv = transform.inverse()?.parameters;
|
||||
|
||||
@@ -1457,28 +1521,26 @@ mod tests {
|
||||
use crate::metric::{MattesMetric, SamplingArg};
|
||||
use algos::ObjectiveFunction;
|
||||
|
||||
let shape = [200, 200];
|
||||
let center = [99.5, 99.5];
|
||||
let im_a = julia_image(
|
||||
&shape,
|
||||
&[200, 200],
|
||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||
¢er,
|
||||
&[99.5, 99.5],
|
||||
&[-0.8, 0.156],
|
||||
)
|
||||
.mapv(|i| i as f64);
|
||||
let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, ¢er);
|
||||
let im_b = rotation.interpolate::<3, _, _>(im_a.view())?;
|
||||
let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, &[0.0, 0.0]);
|
||||
let im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
|
||||
|
||||
let q_inv = rotation.inverse()?.parameters;
|
||||
let p = rotation.parameters.clone();
|
||||
let identity = vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0];
|
||||
let edge = 0.01;
|
||||
|
||||
// Test at full resolution with different sigma levels (matching elastix approach)
|
||||
// Test at full resolution with different sigma levels (matching Elastix approach)
|
||||
// Elastix does NOT downsample — it smooths at full resolution with
|
||||
// sigma = 0.5 * factor where factor is from the pyramid schedule [8, 4, 2, 1]
|
||||
println!(
|
||||
"=== Full resolution (no downsampling, matching elastix FixedSmoothingImagePyramid) ==="
|
||||
"=== Full resolution (no downsampling, matching Elastix FixedSmoothingImagePyramid) ==="
|
||||
);
|
||||
for (level, sigma_val) in [4.0, 2.0, 1.0, 0.5].iter().enumerate() {
|
||||
let f = gaussian_smooth(im_a.view(), &[*sigma_val; 2])?;
|
||||
@@ -1536,20 +1598,18 @@ mod tests {
|
||||
use crate::filter::gaussian_smooth;
|
||||
use crate::metric::{MattesMetric, SamplingArg};
|
||||
use algos::ObjectiveFunction;
|
||||
let shape = [200, 200];
|
||||
let center = [99.5, 99.5];
|
||||
let im_a = julia_image(
|
||||
&shape,
|
||||
&[200, 200],
|
||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||
¢er,
|
||||
&[99.5, 99.5],
|
||||
&[-0.8, 0.156],
|
||||
)
|
||||
.mapv(|i| i as f64);
|
||||
let edge = 0.01;
|
||||
// Full resolution, sigma=4.0 — matching elastix FixedSmoothingImagePyramid (no downsampling)
|
||||
for angle_deg in [5.0f64, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0] {
|
||||
let rotation = Transform::<Ix2>::from_rotation(angle_deg.to_radians(), ¢er);
|
||||
let im_b = rotation.interpolate::<3, _, _>(im_a.view())?;
|
||||
let rotation = Transform::<Ix2>::from_rotation(angle_deg.to_radians(), &[0.0, 0.0]);
|
||||
let im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
|
||||
let q_inv = rotation.inverse()?.parameters;
|
||||
let p = rotation.parameters.clone();
|
||||
let f = gaussian_smooth(im_a.view(), &[4.0, 4.0])?;
|
||||
@@ -1684,7 +1744,7 @@ mod tests {
|
||||
fn register_real_images2() -> Result<(), Box<dyn std::error::Error>> {
|
||||
let fixed = read_tiff("test_files/fixed.tif")?;
|
||||
let e = Transform::<Ix2>::new(vec![0.8, 0.0, 0.0, 1.0, 0.0, 0.0]).inverse()?;
|
||||
let moving = e.interpolate::<3, _, _>(fixed.view())?;
|
||||
let moving = e.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
|
||||
|
||||
let t = Transform::<Ix2>::register(
|
||||
fixed.view(),
|
||||
|
||||
Reference in new issue
Block a user