Compare commits
1
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
e5136f00bb |
No files matched your search
@@ -59,7 +59,7 @@ fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform<Ix2>,
|
|||||||
The returned `Transform` maps coordinates from the *fixed* image to the *moving* image, i.e.
|
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:
|
`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 transform = register(&fixed, &moving)?;
|
||||||
// out-of-range pixels are filled with the default value (0.0 by default)
|
// 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
|
let registered = transform.interpolate_par::<3, _, _>(moving.view(), &Default::default())?; // cubic B-spline, parallel
|
||||||
@@ -158,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
|
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):
|
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::error::Error;
|
||||||
use image_registration::julia_image;
|
use image_registration::julia_image;
|
||||||
use image_registration::transform::Transform;
|
use image_registration::transform::Transform;
|
||||||
@@ -196,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
|
### 4. Registration with inspection of every optimization step
|
||||||
|
|
||||||
```rust
|
```ignore
|
||||||
use image_registration::error::Error;
|
use image_registration::error::Error;
|
||||||
use image_registration::transform::Transform;
|
use image_registration::transform::Transform;
|
||||||
|
|
||||||
@@ -223,7 +227,7 @@ fn main() -> Result<(), Error> {
|
|||||||
following uses two levels: a coarse level with strong Gaussian smoothing and the robust ASGD optimizer, then a fine
|
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.
|
level with L-BFGS for accuracy.
|
||||||
|
|
||||||
```rust
|
```ignore
|
||||||
use image_registration::error::Error;
|
use image_registration::error::Error;
|
||||||
use image_registration::metric::{SamplingArg, Sigma};
|
use image_registration::metric::{SamplingArg, Sigma};
|
||||||
use image_registration::register::{Optimizer, Registration, RegistrationStep};
|
use image_registration::register::{Optimizer, Registration, RegistrationStep};
|
||||||
@@ -266,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,
|
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:
|
register for translation only:
|
||||||
|
|
||||||
```rust
|
```ignore
|
||||||
use image_registration::error::Error;
|
use image_registration::error::Error;
|
||||||
use image_registration::transform::Transform;
|
use image_registration::transform::Transform;
|
||||||
|
|
||||||
@@ -299,6 +303,7 @@ Build a `MattesMetric` and evaluate the (negative) mutual information and its gr
|
|||||||
use algos::ObjectiveFunction; // direct dependency of image-registration; add it to your crate to name the trait
|
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::bspline::{BSpline, BSplineTrait};
|
||||||
use image_registration::error::Error;
|
use image_registration::error::Error;
|
||||||
|
use image_registration::julia_image;
|
||||||
use image_registration::metric::{MattesMetric, SamplingArg, Sigma};
|
use image_registration::metric::{MattesMetric, SamplingArg, Sigma};
|
||||||
use image_registration::transform::Transform;
|
use image_registration::transform::Transform;
|
||||||
|
|
||||||
|
|||||||
+4
-5
@@ -113,7 +113,6 @@ mod tests {
|
|||||||
use super::*;
|
use super::*;
|
||||||
use crate::julia_image;
|
use crate::julia_image;
|
||||||
use ndarray::{Ix1, array};
|
use ndarray::{Ix1, array};
|
||||||
use tiffwrite::IJTiffFile;
|
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn smooth() -> Result<(), Box<dyn std::error::Error>> {
|
fn smooth() -> Result<(), Box<dyn std::error::Error>> {
|
||||||
@@ -135,10 +134,10 @@ mod tests {
|
|||||||
&[29.5, 39.5],
|
&[29.5, 39.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
);
|
);
|
||||||
let im_b = gaussian_smooth(im_a.mapv(|i| i as f64).view(), &[4.0, 4.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"))?;
|
// 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_a.mapv(|i| i as u32).view(), 0, 0, 0)?;
|
||||||
t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
|
// t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
|
||||||
Ok(())
|
Ok(())
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
+5
-6
@@ -50,7 +50,6 @@ mod tests {
|
|||||||
use super::*;
|
use super::*;
|
||||||
use crate::transform::Transform;
|
use crate::transform::Transform;
|
||||||
use tempfile::NamedTempFile;
|
use tempfile::NamedTempFile;
|
||||||
use tiffwrite::IJTiffFile;
|
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn serialization() -> Result<(), Error> {
|
fn serialization() -> Result<(), Error> {
|
||||||
@@ -81,7 +80,7 @@ mod tests {
|
|||||||
.collect::<Vec<_>>()
|
.collect::<Vec<_>>()
|
||||||
.try_into()
|
.try_into()
|
||||||
.unwrap();
|
.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 transform_j = transform.parameters.clone().try_into().unwrap();
|
||||||
let c = [-0.8, 0.156];
|
let c = [-0.8, 0.156];
|
||||||
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
||||||
@@ -121,7 +120,7 @@ mod tests {
|
|||||||
.collect::<Vec<_>>()
|
.collect::<Vec<_>>()
|
||||||
.try_into()
|
.try_into()
|
||||||
.unwrap();
|
.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 transform_j = transform.parameters.clone().try_into().unwrap();
|
||||||
let c = [-0.8, 0.156];
|
let c = [-0.8, 0.156];
|
||||||
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64);
|
||||||
@@ -158,14 +157,14 @@ mod tests {
|
|||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn julia() -> Result<(), Box<dyn std::error::Error>> {
|
fn julia() -> Result<(), Box<dyn std::error::Error>> {
|
||||||
let image_a = julia_image(
|
let _image_a = julia_image(
|
||||||
&[60, 80],
|
&[60, 80],
|
||||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||||
&[29.5, 39.5],
|
&[29.5, 39.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
);
|
);
|
||||||
let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
// let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
|
||||||
t.save(image_a.view(), 0, 0, 0)?;
|
// t.save(image_a.view(), 0, 0, 0)?;
|
||||||
Ok(())
|
Ok(())
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
+9
-13
@@ -1063,7 +1063,7 @@ mod tests {
|
|||||||
&[3999.5, 2999.5],
|
&[3999.5, 2999.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
);
|
);
|
||||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
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_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
|
||||||
let image_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
let image_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
||||||
|
|
||||||
@@ -1084,7 +1084,7 @@ mod tests {
|
|||||||
&[3999.5, 2999.5],
|
&[3999.5, 2999.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
);
|
);
|
||||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
|
||||||
let image_b = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
let image_b = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
|
||||||
assert_ne!(image_b[[4000, 3000]], 0.0);
|
assert_ne!(image_b[[4000, 3000]], 0.0);
|
||||||
Ok(())
|
Ok(())
|
||||||
@@ -1098,7 +1098,7 @@ mod tests {
|
|||||||
&[3999.5, 2999.5],
|
&[3999.5, 2999.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
);
|
);
|
||||||
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
|
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_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
|
||||||
assert_ne!(image_b[[4000, 3000]], 0.0);
|
assert_ne!(image_b[[4000, 3000]], 0.0);
|
||||||
Ok(())
|
Ok(())
|
||||||
@@ -1521,16 +1521,14 @@ mod tests {
|
|||||||
use crate::metric::{MattesMetric, SamplingArg};
|
use crate::metric::{MattesMetric, SamplingArg};
|
||||||
use algos::ObjectiveFunction;
|
use algos::ObjectiveFunction;
|
||||||
|
|
||||||
let shape = [200, 200];
|
|
||||||
let center = [99.5, 99.5];
|
|
||||||
let im_a = julia_image(
|
let im_a = julia_image(
|
||||||
&shape,
|
&[200, 200],
|
||||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||||
¢er,
|
&[99.5, 99.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
)
|
)
|
||||||
.mapv(|i| i as f64);
|
.mapv(|i| i as f64);
|
||||||
let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, ¢er);
|
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 im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
|
||||||
|
|
||||||
let q_inv = rotation.inverse()?.parameters;
|
let q_inv = rotation.inverse()?.parameters;
|
||||||
@@ -1600,19 +1598,17 @@ mod tests {
|
|||||||
use crate::filter::gaussian_smooth;
|
use crate::filter::gaussian_smooth;
|
||||||
use crate::metric::{MattesMetric, SamplingArg};
|
use crate::metric::{MattesMetric, SamplingArg};
|
||||||
use algos::ObjectiveFunction;
|
use algos::ObjectiveFunction;
|
||||||
let shape = [200, 200];
|
|
||||||
let center = [99.5, 99.5];
|
|
||||||
let im_a = julia_image(
|
let im_a = julia_image(
|
||||||
&shape,
|
&[200, 200],
|
||||||
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
&[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
|
||||||
¢er,
|
&[99.5, 99.5],
|
||||||
&[-0.8, 0.156],
|
&[-0.8, 0.156],
|
||||||
)
|
)
|
||||||
.mapv(|i| i as f64);
|
.mapv(|i| i as f64);
|
||||||
let edge = 0.01;
|
let edge = 0.01;
|
||||||
// Full resolution, sigma=4.0 — matching elastix FixedSmoothingImagePyramid (no downsampling)
|
// 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] {
|
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 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 im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
|
||||||
let q_inv = rotation.inverse()?.parameters;
|
let q_inv = rotation.inverse()?.parameters;
|
||||||
let p = rotation.parameters.clone();
|
let p = rotation.parameters.clone();
|
||||||
|
|||||||
Reference in new issue
Block a user