- a comment and a TODO

Cargo.toml

@@ -1,6 +1,6 @@
 [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>"]

README.md

@@ -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.
 `moving(transform_point(x)) ≈ fixed(x)`. To warp the moving image into the fixed image's frame:
 
-```rust
+```ignore
 let transform = register(&fixed, &moving)?;
 // 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
@@ -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
 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;
@@ -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
 
-```rust
+```ignore
 use image_registration::error::Error;
 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
 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};
@@ -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,
 register for translation only:
 
-```rust
+```ignore
 use image_registration::error::Error;
 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 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;
 

src/bspline.rs

@@ -22,6 +22,7 @@ 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")]

src/filter.rs

@@ -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(())
     }
 

src/lib.rs

@@ -1,4 +1,4 @@
-#![cfg_attr(docsrs, feature(doc_cfg, doc_auto_cfg))]
+#![cfg_attr(docsrs, feature(doc_cfg))]
 #![doc = include_str!("../README.md")]
 use rayon::iter::ParallelIterator;
 pub mod bspline;
@@ -50,7 +50,6 @@ mod tests {
     use super::*;
     use crate::transform::Transform;
     use tempfile::NamedTempFile;
-    use tiffwrite::IJTiffFile;
 
     #[test]
     fn serialization() -> Result<(), Error> {
@@ -81,7 +80,7 @@ mod tests {
             .collect::<Vec<_>>()
             .try_into()
             .unwrap();
-        let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &center);
+        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, &center, &c).mapv(|x| x as f64);
@@ -121,7 +120,7 @@ mod tests {
             .collect::<Vec<_>>()
             .try_into()
             .unwrap();
-        let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &center);
+        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, &center, &c).mapv(|x| x as f64);
@@ -158,14 +157,14 @@ mod tests {
 
     #[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(())
     }
 }

src/transform.rs

@@ -1063,7 +1063,7 @@ mod tests {
             &[3999.5, 2999.5],
             &[-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_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
 
@@ -1084,7 +1084,7 @@ mod tests {
             &[3999.5, 2999.5],
             &[-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())?;
         assert_ne!(image_b[[4000, 3000]], 0.0);
         Ok(())
@@ -1098,7 +1098,7 @@ mod tests {
             &[3999.5, 2999.5],
             &[-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())?;
         assert_ne!(image_b[[4000, 3000]], 0.0);
         Ok(())
@@ -1521,16 +1521,14 @@ 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],
-            &center,
+            &[99.5, 99.5],
             &[-0.8, 0.156],
         )
         .mapv(|i| i as f64);
-        let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, &center);
+        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;
@@ -1600,19 +1598,17 @@ 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],
-            &center,
+            &[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(), &center);
+            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();