wim/sitk-registration-sys
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- add more tests

Browse Source 
- do not allow multiple registrations to run at once because it causes a memory error in the sitk code
This commit is contained in:
Wim Pomp
2025-03-03 09:40:46 +01:00
parent 0aeb3c8c7e
commit 934f038ea1
8 changed files with 378 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/lib.rs
+157 -37
View File
@@ -3,10 +3,13 @@ mod sys;
use crate::sys::{PixelType, interp, register};
use anyhow::{Result, anyhow};
use ndarray::{Array2, ArrayView2, array, s};
use serde::{Deserialize, Serialize};
use serde_yaml::{from_reader, to_writer};
use std::fs::File;
use std::ops::Mul;
use std::path::PathBuf;
#[derive(Clone, Debug)]
#[derive(Clone, Debug, Deserialize, Serialize)]
pub struct Transform {
pub parameters: [f64; 6],
pub dparameters: [f64; 6],
@@ -94,13 +97,16 @@ impl Transform {
}
/// read a transform from a file
pub fn from_file(file: PathBuf) -> Result<Self> {
todo!()
pub fn from_file(path: PathBuf) -> Result<Self> {
let file = File::open(path)?;
Ok(from_reader(file)?)
}
/// write a transform to a file
pub fn to_file(&self, file: PathBuf) -> Result<()> {
todo!()
pub fn to_file(&self, path: PathBuf) -> Result<()> {
let mut file = File::open(path)?;
to_writer(&mut file, self)?;
Ok(())
}
/// true if transform does nothing
@@ -190,18 +196,17 @@ impl Transform {
}
let m = self.matrix();
let d0 = det(m.slice(s![1.., 1..]));
if d0 == 0f64 {
let d = det(m.slice(s![..2, ..2]));
if d == 0f64 {
return Err(anyhow!("transform matrix is not invertible"));
}
let d2 = det(m.slice(s![..2, ..2]));
let parameters = [
d0 / d2,
-det(m.slice(s![..;2, 1..])) / d2,
-det(m.slice(s![1.., ..;2])) / d2,
det(m.slice(s![..;2, ..;2])) / d2,
det(m.slice(s![..2, 1..])) / d2,
-det(m.slice(s![..2, ..;2])) / d2,
det(m.slice(s![1.., 1..])) / d,
-det(m.slice(s![..;2, 1..])) / d,
-det(m.slice(s![1.., ..;2])) / d,
det(m.slice(s![..;2, ..;2])) / d,
det(m.slice(s![..2, 1..])) / d,
-det(m.slice(s![..2, ..;2])) / d,
];
Ok(Transform {
@@ -228,10 +233,9 @@ mod tests {
use anyhow::Result;
use ndarray::Array2;
use num::Complex;
use tiffwrite::IJTiffFile;
/// An example of generating julia fractals.
fn julia_image() -> Result<Array2<u8>> {
fn julia_image(shift_x: f32, shift_y: f32) -> Result<Array2<u8>> {
let imgx = 800;
let imgy = 600;
@@ -241,11 +245,11 @@ mod tests {
let mut im = Array2::<u8>::zeros((imgy, imgx));
for x in 0..imgx {
for y in 0..imgy {
let cx = y as f32 * scalex - 1.5;
let cy = x as f32 * scaley - 1.5;
let cy = (y as f32 + shift_y) * scalex - 1.5;
let cx = (x as f32 + shift_x) * scaley - 1.5;
let c = Complex::new(-0.4, 0.6);
let mut z = Complex::new(cx, cy);
let mut z = Complex::new(cy, cx);
let mut i = 0;
while i < 255 && z.norm() <= 2.0 {
@@ -259,24 +263,140 @@ mod tests {
Ok(im)
}
#[test]
fn test_interp() -> Result<()> {
let j = julia_image()?;
let mut tif = IJTiffFile::new("interp_test.tif")?;
tif.save(&j, 0, 0, 0)?;
let shape = j.shape();
let origin = [
((shape[1] - 1) as f64) / 2f64,
((shape[0] - 1) as f64) / 2f64,
];
let transform = Transform::new([1.2, 0., 0., 1., 10., 0.], origin, [shape[0], shape[1]]);
let k = transform.transform_image_bspline(j.view())?;
tif.save(&k, 1, 0, 0)?;
macro_rules! interp_tests_bspline {
($($name:ident: $t:ty $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let j = julia_image(-120f32, 10f32)?.mapv(|x| x as $t);
let k = julia_image(0f32, 0f32)?.mapv(|x| x as $t);
let shape = j.shape();
let origin = [
((shape[1] - 1) as f64) / 2f64,
((shape[0] - 1) as f64) / 2f64,
];
let transform = Transform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
let n = transform.transform_image_bspline(j.view())?;
let d = (k.mapv(|x| x as f64) - n.mapv(|x| x as f64)).powi(2).sum();
assert!(d <= (shape[0] * shape[1]) as f64);
Ok(())
}
)*
}
}
let t = Transform::register_affine(k.view(), j.view())?;
println!("t: {:#?}", t);
println!("m: {:#?}", t.matrix());
println!("i: {:#?}", t.inverse()?.matrix());
Ok(())
interp_tests_bspline! {
interpbs_u8: u8,
interpbs_i8: i8,
interpbs_u16: u16,
interpbs_i16: i16,
interpbs_u32: u32,
interpbs_i32: i32,
interpbs_u64: u64,
interpbs_i64: i64,
interpbs_f32: f32,
interpbs_f64: f64,
}
macro_rules! interp_tests_nearest_neighbor {
($($name:ident: $t:ty $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let j = julia_image(-120f32, 10f32)?.mapv(|x| x as $t);
let k = julia_image(0f32, 0f32)?.mapv(|x| x as $t);
let shape = j.shape();
let origin = [
((shape[1] - 1) as f64) / 2f64,
((shape[0] - 1) as f64) / 2f64,
];
let transform = Transform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
let n = transform.transform_image_nearest_neighbor(j.view())?;
let d = (k.mapv(|x| x as f64) - n.mapv(|x| x as f64)).powi(2).sum();
assert!(d <= (shape[0] * shape[1]) as f64);
Ok(())
}
)*
}
}
interp_tests_nearest_neighbor! {
interpnn_u8: u8,
interpnn_i8: i8,
interpnn_u16: u16,
interpnn_i16: i16,
interpnn_u32: u32,
interpnn_i32: i32,
interpnn_u64: u64,
interpnn_i64: i64,
interpnn_f32: f32,
interpnn_f64: f64,
}
macro_rules! registration_tests_translation {
($($name:ident: $t:ty $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let j = julia_image(0f32, 0f32)?.mapv(|x| x as $t);
let k = julia_image(10f32, 20f32)?.mapv(|x| x as $t);
let t = Transform::register_translation(j.view(), k.view())?;
let mut m = Array2::eye(3);
m[[0, 2]] = -10f64;
m[[1, 2]] = -20f64;
let d = (t.matrix() - m).powi(2).sum();
assert!(d < 0.01);
Ok(())
}
)*
}
}
registration_tests_translation! {
registration_translation_u8: u8,
registration_translation_i8: i8,
registration_translation_u16: u16,
registration_translation_i16: i16,
registration_translation_u32: u32,
registration_translation_i32: i32,
registration_translation_u64: u64,
registration_translation_i64: i64,
registration_translation_f32: f32,
registration_translation_f64: f64,
}
macro_rules! registration_tests_affine {
($($name:ident: $t:ty $(,)?)*) => {
$(
#[test]
fn $name() -> Result<()> {
let j = julia_image(0f32, 0f32)?.mapv(|x| x as $t);
let shape = j.shape();
let origin = [
((shape[1] - 1) as f64) / 2f64,
((shape[0] - 1) as f64) / 2f64,
];
let s = Transform::new([1.2, 0., 0., 1., 5., 7.], origin, [shape[0], shape[1]]);
let k = s.transform_image_bspline(j.view())?;
let t = Transform::register_affine(j.view(), k.view())?.inverse()?;
let d = (t.matrix() - s.matrix()).powi(2).sum();
assert!(d < 0.01);
Ok(())
}
)*
}
}
registration_tests_affine! {
registration_tests_affine_u8: u8,
registration_tests_affine_i8: i8,
registration_tests_affine_u16: u16,
registration_tests_affine_i16: i16,
registration_tests_affine_u32: u32,
registration_tests_affine_i32: i32,
registration_tests_affine_u64: u64,
registration_tests_affine_i64: i64,
registration_tests_affine_f32: f32,
registration_tests_affine_f64: f64,
}
}
src/sys.rs
+105 -78
View File
@@ -1,55 +1,65 @@
use anyhow::Result;
use libc::{c_double, c_uint};
use ndarray::{Array2, ArrayView2};
use one_at_a_time_please::one_at_a_time;
use std::ptr;
macro_rules! register_fn {
($T:ty, $t:ident) => {
fn $t(
width: c_uint,
height: c_uint,
fixed_arr: *const $T,
moving_arr: *const $T,
translation_or_affine: bool,
transform: &mut *mut c_double,
);
($($name:ident: $T:ty $(,)?)*) => {
$(
fn $name(
width: c_uint,
height: c_uint,
fixed_arr: *const $T,
moving_arr: *const $T,
translation_or_affine: bool,
transform: &mut *mut c_double,
);
)*
};
}
macro_rules! interp_fn {
($T:ty, $t:ident) => {
fn $t(
width: c_uint,
height: c_uint,
transform: *const c_double,
origin: *const c_double,
image: &mut *mut $T,
bspline_or_nn: bool,
);
($($name:ident: $T:ty $(,)?)*) => {
$(
fn $name(
width: c_uint,
height: c_uint,
transform: *const c_double,
origin: *const c_double,
image: &mut *mut $T,
bspline_or_nn: bool,
);
)*
};
}
unsafe extern "C" {
register_fn!(u8, register_u8);
register_fn!(i8, register_i8);
register_fn!(u16, register_u16);
register_fn!(i16, register_i16);
register_fn!(u32, register_u32);
register_fn!(i32, register_i32);
register_fn!(u64, register_u64);
register_fn!(i64, register_i64);
register_fn!(f32, register_f32);
register_fn!(f64, register_f64);
interp_fn!(u8, interp_u8);
interp_fn!(i8, interp_i8);
interp_fn!(u16, interp_u16);
interp_fn!(i16, interp_i16);
interp_fn!(u32, interp_u32);
interp_fn!(i32, interp_i32);
interp_fn!(u64, interp_u64);
interp_fn!(i64, interp_i64);
interp_fn!(f32, interp_f32);
interp_fn!(f64, interp_f64);
register_fn! {
register_u8: u8,
register_i8: i8,
register_u16: u16,
register_i16: i16,
register_u32: u32,
register_i32: i32,
register_u64: u64,
register_i64: i64,
register_f32: f32,
register_f64: f64,
}
interp_fn! {
interp_u8: u8,
interp_i8: i8,
interp_u16: u16,
interp_i16: i16,
interp_u32: u32,
interp_i32: i32,
interp_u64: u64,
interp_i64: i64,
interp_f32: f32,
interp_f64: f64,
}
}
pub trait PixelType: Clone {
@@ -57,31 +67,36 @@ pub trait PixelType: Clone {
}
macro_rules! sitk_impl {
($T:ty, $sitk:expr) => {
impl PixelType for $T {
const PT: u8 = $sitk;
}
($($T:ty: $sitk:expr $(,)?)*) => {
$(
impl PixelType for $T {
const PT: u8 = $sitk;
}
)*
};
}
sitk_impl!(u8, 1);
sitk_impl!(i8, 2);
sitk_impl!(u16, 3);
sitk_impl!(i16, 4);
sitk_impl!(u32, 5);
sitk_impl!(i32, 6);
sitk_impl!(u64, 7);
sitk_impl!(i64, 8);
sitk_impl! {
u8: 1,
i8: 2,
u16: 3,
i16: 4,
u32: 5,
i32: 6,
u64: 7,
i64: 8,
f32: 9,
f64: 10,
}
#[cfg(target_pointer_width = "64")]
sitk_impl!(usize, 7);
sitk_impl!(usize: 7);
#[cfg(target_pointer_width = "32")]
sitk_impl!(usize, 5);
sitk_impl!(usize: 5);
#[cfg(target_pointer_width = "64")]
sitk_impl!(isize, 8);
sitk_impl!(isize: 8);
#[cfg(target_pointer_width = "32")]
sitk_impl!(isize, 6);
sitk_impl!(f32, 9);
sitk_impl!(f64, 10);
sitk_impl!(isize: 6);
pub(crate) fn interp<T: PixelType>(
parameters: [f64; 6],
@@ -204,6 +219,7 @@ pub(crate) fn interp<T: PixelType>(
)?)
}
#[one_at_a_time]
pub(crate) fn register<T: PixelType>(
fixed: ArrayView2<T>,
moving: ArrayView2<T>,
@@ -212,19 +228,24 @@ pub(crate) fn register<T: PixelType>(
let shape: Vec<usize> = fixed.shape().to_vec();
let width = shape[1] as c_uint;
let height = shape[0] as c_uint;
let fixed: Vec<_> = fixed.into_iter().collect();
let moving: Vec<_> = moving.into_iter().collect();
let fixed: Vec<_> = fixed.into_iter().cloned().collect();
let moving: Vec<_> = moving.into_iter().cloned().collect();
let fixed_ptr = fixed.as_ptr();
let moving_ptr = moving.as_ptr();
let mut transform: Vec<c_double> = vec![0.0; 6];
let mut transform_ptr: *mut c_double = ptr::from_mut(unsafe { &mut *transform.as_mut_ptr() });
// let ma0 = &mut moving as *mut Vec<T> as usize;
// println!("ma0: {:#x}", ma0);
match T::PT {
1 => {
unsafe {
register_u8(
width,
height,
fixed.as_ptr() as *const u8,
moving.as_ptr() as *const u8,
fixed_ptr as *const u8,
moving_ptr as *const u8,
translation_or_affine,
&mut transform_ptr,
)
@@ -235,8 +256,8 @@ pub(crate) fn register<T: PixelType>(
register_i8(
width,
height,
fixed.as_ptr() as *const i8,
moving.as_ptr() as *const i8,
fixed_ptr as *const i8,
moving_ptr as *const i8,
translation_or_affine,
&mut transform_ptr,
)
@@ -247,8 +268,8 @@ pub(crate) fn register<T: PixelType>(
register_u16(
width,
height,
fixed.as_ptr() as *const u16,
moving.as_ptr() as *const u16,
fixed_ptr as *const u16,
moving_ptr as *const u16,
translation_or_affine,
&mut transform_ptr,
)
@@ -259,8 +280,8 @@ pub(crate) fn register<T: PixelType>(
register_i16(
width,
height,
fixed.as_ptr() as *const i16,
moving.as_ptr() as *const i16,
fixed_ptr as *const i16,
moving_ptr as *const i16,
translation_or_affine,
&mut transform_ptr,
)
@@ -271,8 +292,8 @@ pub(crate) fn register<T: PixelType>(
register_u32(
width,
height,
fixed.as_ptr() as *const u32,
moving.as_ptr() as *const u32,
fixed_ptr as *const u32,
moving_ptr as *const u32,
translation_or_affine,
&mut transform_ptr,
)
@@ -283,8 +304,8 @@ pub(crate) fn register<T: PixelType>(
register_i32(
width,
height,
fixed.as_ptr() as *const i32,
moving.as_ptr() as *const i32,
fixed_ptr as *const i32,
moving_ptr as *const i32,
translation_or_affine,
&mut transform_ptr,
)
@@ -295,8 +316,8 @@ pub(crate) fn register<T: PixelType>(
register_u64(
width,
height,
fixed.as_ptr() as *const u64,
moving.as_ptr() as *const u64,
fixed_ptr as *const u64,
moving_ptr as *const u64,
translation_or_affine,
&mut transform_ptr,
)
@@ -307,8 +328,8 @@ pub(crate) fn register<T: PixelType>(
register_i64(
width,
height,
fixed.as_ptr() as *const i64,
moving.as_ptr() as *const i64,
fixed_ptr as *const i64,
moving_ptr as *const i64,
translation_or_affine,
&mut transform_ptr,
)
@@ -319,8 +340,8 @@ pub(crate) fn register<T: PixelType>(
register_f32(
width,
height,
fixed.as_ptr() as *const f32,
moving.as_ptr() as *const f32,
fixed_ptr as *const f32,
moving_ptr as *const f32,
translation_or_affine,
&mut transform_ptr,
)
@@ -331,8 +352,8 @@ pub(crate) fn register<T: PixelType>(
register_f64(
width,
height,
fixed.as_ptr() as *const f64,
moving.as_ptr() as *const f64,
fixed_ptr as *const f64,
moving_ptr as *const f64,
translation_or_affine,
&mut transform_ptr,
)
@@ -341,6 +362,12 @@ pub(crate) fn register<T: PixelType>(
_ => {}
}
// let ma1 = &mut moving as *mut Vec<T> as usize;
// println!("ma1: {:#x}", ma1);
// println!("{}", fixed.len());
// println!("{}", moving.len());
Ok((
[
transform[0] as f64,