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

Now build SimpleITK into static libs and use (auto)cxx.

Browse Source
This commit is contained in:
Wim Pomp
2025-09-12 18:03:24 +02:00
parent 0c6a38e4fa
commit 5750fd7f99
17 changed files with 8110 additions and 1426 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@@ -83,3 +83,6 @@ docs/_build/
*.nii *.nii
*.tif *.tif
TransformParameters* TransformParameters*
/cxx/CMakeLists.txt
/cxx/cmake-build-debug/
/src/autocxx-ffi-default-gen_unfmt.rs

3
.gitmodules vendored Normal file
View File

@@ -0,0 +1,3 @@
[submodule "SimpleITK"]
path = SimpleITK
url = https://github.com/SimpleITK/SimpleITK.git

29
Cargo.toml
View File

@@ -1,10 +1,10 @@
[package] [package]
name = "sitk-registration-sys" name = "sitk-registration-sys"
version = "2025.3.2" version = "0.2.0"
edition = "2024" edition = "2024"
license = "MIT OR Apache-2.0" license = "Apache-2.0"
description = "register and interpolate images" description = "register and interpolate images"
rust-version = "1.85.0" rust-version = "1.85.1"
authors = ["Wim Pomp <w.pomp@nki.nl>"] authors = ["Wim Pomp <w.pomp@nki.nl>"]
homepage = "https://github.com/wimpomp/sitk-registration-sys" homepage = "https://github.com/wimpomp/sitk-registration-sys"
repository = "https://github.com/wimpomp/sitk-registration-sys" repository = "https://github.com/wimpomp/sitk-registration-sys"
@@ -12,23 +12,32 @@ documentation = "https://docs.rs/sitk-registration-sys"
readme = "README.md" readme = "README.md"
keywords = ["registration", "affine", "bspline", "transform"] keywords = ["registration", "affine", "bspline", "transform"]
categories = ["multimedia::images", "science"] categories = ["multimedia::images", "science"]
exclude = ["SimpleITK/Testing/*", "SimpleITK/Examples/*", "SimpleITK/Wrapping/*"]
[lib] [lib]
name = "sitk_registration_sys" name = "sitk_registration_sys"
crate-type = ["cdylib", "rlib"] crate-type = ["cdylib", "rlib"]
[dependencies] [dependencies]
anyhow = "1.0.97" anyhow = "1.0.99"
libc = "0.2.170" autocxx = "0.30.0"
cxx = "1.0.183"
link-cplusplus = "1.0.12"
ndarray = "0.16.1" ndarray = "0.16.1"
num = "0.4.3" num = "0.4.3"
one_at_a_time_please = "1.0.1" serde = { version = "1.0.219", features = ["derive"] }
serde = { version = "1.0.218", features = ["derive"] }
serde_yaml = "0.9.33" serde_yaml = "0.9.33"
tempfile = "3.22.0"
[build-dependencies] [build-dependencies]
autocxx-build = "0.30.0"
cc = "1.2.36"
cmake = "0.1.54" cmake = "0.1.54"
git2 = "0.20.0" regex = "1.11.2"
[dev-dependencies] [profile.dev]
tempfile = "3.18.0" lto = true
strip = "symbols"
[features]
sitk_no_build = [] # do not build all the C++ stuff, just use the rust code to make docs etc.

0
LICENSE-APACHE → LICENSE
View File

23
LICENSE-MIT
View File

@@ -1,23 +0,0 @@
Permission is hereby granted, free of charge, to any
person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the
Software without restriction, including without
limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following
conditions:
The above copyright notice and this permission notice
shall be included in all copies or substantial portions
of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF
ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT
SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.

47
README.md
View File

@@ -2,30 +2,47 @@
This crate does two things: This crate does two things:
- find an affine transform or translation that transforms one image into the other - find an affine transform or translation that transforms one image into the other
- use bpline or nearest neighbor interpolation to apply a transformation to an image - use bspline or nearest neighbor interpolation to apply a transformation to an image
To do this, [SimpleITK](https://github.com/SimpleITK/SimpleITK.git), which is written in To do this, [SimpleITK](https://github.com/SimpleITK/SimpleITK.git), which is written in
C++, is used. An adapter library is created to expose the required functionality in SimpleITK C++, is used. An adapter library is created using [autocxx](https://crates.io/crates/autocxx)
in a shared library. Because of this, compilation of this crate requires quite some time, as to expose the required functionality in SimpleITK. Because of this, compilation of this crate
wel as cmake. requires quite some time, several GB of memory, up to 50 GB of hard disk space, as well as
cmake, a C++ compiler, llvm and git. Use at your own risk!
## Examples ## Examples
### Registration ### Registration
``` ```
let image_a = (some Array2); use ndarray::Array2;
let iameg_b = (some transformed Array2); use sitk_registration_sys::registration::{AffineTransform, julia_image};
let transform = Transform::register_affine(image_a.view(), image_b.view())?;
println!("transform: {:#?}", transform);
```
### Interpolation let j = julia_image(0f32, 0f32).unwrap();
``` let shape = j.shape();
let image = (Some Array2);
let shape = image.shape();
let origin = [ let origin = [
((shape[1] - 1) as f64) / 2f64, ((shape[1] - 1) as f64) / 2f64,
((shape[0] - 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 s = AffineTransform::new([1.2, 0., 0., 1., 5., 7.], origin, [shape[0], shape[1]]);
let transformed_image = transform.transform_image_bspline(image.view())?; let k: Array2<_> = s.transform_image_bspline(j.view()).unwrap().into();
let t = AffineTransform::register_affine(j.view(), k.view()).unwrap().inverse().unwrap();
let d = (t.matrix() - s.matrix()).powi(2).sum();
assert!(d < 0.025, "d: {}, t: {:?}", d, t.parameters);
```
### Interpolation
```
use ndarray::Array2;
use sitk_registration_sys::registration::{AffineTransform, julia_image};
let j = julia_image(-120f32, 10f32).unwrap();
let k = julia_image(0f32, 0f32).unwrap();
let shape = j.shape();
let origin = [
((shape[1] - 1) as f64) / 2f64,
((shape[0] - 1) as f64) / 2f64,
];
let transform = AffineTransform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
let n: Array2<_> = transform.transform_image_bspline(j.view()).unwrap().into();
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);
``` ```

1
SimpleITK Submodule

Submodule SimpleITK added at f0254cccec

124
build.rs
View File

@@ -1,34 +1,18 @@
use cmake::Config; use std::env;
use git2::Repository; use std::error::Error;
use std::ffi::OsStr; use std::fs;
use std::io::{Read, Write};
use std::path::PathBuf; use std::path::PathBuf;
use std::process::Command;
fn main() { fn main() -> Result<(), Box<dyn Error>> {
if std::env::var("DOCS_RS").is_err() { if env::var("DOCS_RS").is_err() & !cfg!(feature = "sitk_no_build") {
let out_dir = PathBuf::from(std::env::var("OUT_DIR").expect("OUT_DIR is undefined")); let sitk_dir = PathBuf::from("SimpleITK");
let mut target_dir = out_dir.clone();
while target_dir.file_name() != Some(OsStr::new("target")) {
if !target_dir.pop() {
panic!("Could not find target directory");
}
}
let sitk_dir = if let Some(d) = target_dir.parent() { // use cmake to compile the SimpleITK C++ code
d.join("sitk").to_path_buf() let dst = cmake::Config::new(sitk_dir.join("SuperBuild"))
} else {
target_dir.join("sitk")
};
if !sitk_dir.exists() {
Repository::clone("https://github.com/SimpleITK/SimpleITK.git", &sitk_dir)
.expect("unable to clone sitk");
}
let sitk_build_dir = sitk_dir.join("build");
if !sitk_build_dir.exists() {
println!("cargo::warning=Simple ITK; this will take a long time...");
Config::new(sitk_dir.join("SuperBuild"))
.out_dir(&sitk_dir)
.no_build_target(true) .no_build_target(true)
.define("BUILD_EXAMPLES", "OFF")
.define("BUILD_TESTING", "OFF") .define("BUILD_TESTING", "OFF")
.define("WRAP_CSHARP", "OFF") .define("WRAP_CSHARP", "OFF")
.define("WRAP_JAVA", "OFF") .define("WRAP_JAVA", "OFF")
@@ -40,22 +24,76 @@ fn main() {
.define("WRAP_DEFAULT", "OFF") .define("WRAP_DEFAULT", "OFF")
.define("SimpleITK_USE_ELASTIX", "ON") .define("SimpleITK_USE_ELASTIX", "ON")
.build(); .build();
}
println!( // there does not seem to be any way to tell cargo to group libraries together
"cargo::rustc-env=CMAKE_INSTALL_PREFIX={}", // when linking, so we group all objects together in three big static libraries
out_dir.display() let merged_lib = cmake::Config::new("merged_lib")
); .out_dir(dst.join("merged_lib"))
let path = Config::new("cpp") .no_build_target(true)
.very_verbose(true)
.define("Elastix_DIR", sitk_build_dir.join("Elastix-build"))
.define("ITK_DIR", sitk_build_dir.join("ITK-build"))
.define("SimpleITK_DIR", sitk_build_dir.join("SimpleITK-build"))
.define("CMAKE_INSTALL_PREFIX", out_dir)
.build(); .build();
println!("cargo::rustc-link-arg=-Wl,-rpath,{}", path.display());
println!("cargo::rustc-link-search={}", path.join("build").display()); let mut b = autocxx_build::Builder::new(
println!("cargo::rustc-link-lib=dylib=sitk_adapter"); "src/lib.rs",
println!("cargo::rerun-if-changed=build.rs"); ["cxx", dst.join("include/SimpleITK-3.0").to_str().unwrap()],
println!("cargo::rerun-if-changed=cpp"); )
.extra_clang_args(&["-std=c++17"])
.build()?;
b.flag("-std=c++17").compile("sitk_autocxx");
cc::Build::new()
.include(dst.join("include/SimpleITK-3.0"))
.file("cxx/ffi_extra.cxx")
.compile("sitk_ffi_extra");
let cxx_file = dst.join("autocxx-build-dir/rs/autocxx-ffi-default-gen.rs");
let mut cxx = String::new();
let _ = fs::File::open(&cxx_file)?.read_to_string(&mut cxx)?;
let ra = regex::Regex::new(r"\s*(#|\s(fn|use|mod|pub|struct|enum)\s)\s*")?;
let rb = regex::Regex::new(r"\s(;)\s")?;
let pointer = regex::Regex::new(r"pub type pointer = root\s+::\s+pointer;")?;
let cpp_string = regex::Regex::new(r"pub mod simple \{")?;
let rb_tree = regex::Regex::new(
r"pub type _Rb_tree_insert_return_type = root :: std :: _Node_insert_return < _Iterator , _NodeHandle >;",
)?;
{
let mut f = fs::OpenOptions::new()
.write(true)
.truncate(true)
.open(&cxx_file)?;
write!(
f,
"#[allow(unused_imports)]\n#[allow(unsafe_op_in_unsafe_fn)]\n#[allow(clippy::missing_safety_doc)]\n"
)?;
let cxx = ra.replace_all(&cxx, "\n$1");
let cxx = rb.replace_all(&cxx, "$1\n");
let cxx = pointer.replace_all(&cxx, r"// $0");
let cxx = rb_tree.replace_all(&cxx, r"// $0");
let cxx =
cpp_string.replace_all(&cxx, "pub mod simple {\nuse crate::ffi::ToCppString;");
write!(f, "{}", cxx)?;
} }
fs::copy(dst.join(&cxx_file), "src/autocxx-ffi-default-gen_unfmt.rs")?;
Command::new("rustfmt").arg(&cxx_file).status()?;
fs::copy(dst.join(&cxx_file), "src/autocxx-ffi-default-gen.rs")?;
println!("cargo:warning=merged_lib={}", merged_lib.join("build").display());
println!("cargo:rustc-link-search={}", merged_lib.join("build").display());
println!("cargo:rustc-link-lib=static=sitk_ffi_extra");
println!("cargo:rustc-link-lib=static=sitk");
println!("cargo:rustc-link-lib=static=elastix");
println!("cargo:rustc-link-lib=static=itk");
println!("cargo:rerun-if-changed=build.rs");
println!("cargo:rerun-if-changed=src/lib.rs");
println!("cargo:rerun-if-changed=cxx/ffi_extra.h");
println!("cargo:rerun-if-changed=cxx/ffi_extra.cxx");
} else {
let dst = PathBuf::from(env::var("OUT_DIR")?);
fs::create_dir_all(dst.join("autocxx-build-dir/rs"))?;
fs::copy(
"src/autocxx-ffi-default-gen.rs",
dst.join("autocxx-build-dir/rs/autocxx-ffi-default-gen.rs"),
)?;
}
Ok(())
} }

9
cpp/CMakeLists.txt
View File

@@ -1,9 +0,0 @@
cmake_minimum_required(VERSION 3.16.3)
project(sitk_adapter)
set(ENV{Elastix_DIR} "../sitk/build/Elastix-build" )
set(ENV{ITK_DIR} "../sitk/build/ITK-build" )
set(ENV{SimpleITK_DIR} "~../sitk/build/SimpleITK-build" )
find_package(SimpleITK)
add_library(sitk_adapter SHARED sitk_adapter.cxx)
target_link_libraries (sitk_adapter ${SimpleITK_LIBRARIES})
install(TARGETS sitk_adapter DESTINATION .)

476
cpp/sitk_adapter.cxx
View File

@@ -1,476 +0,0 @@
#include <SimpleITK.h>
#include <sitkImageOperators.h>
#include <cstring>
#include <filesystem>
namespace sitk = itk::simple;
using namespace std;
std::string gen_random(const int len) {
static const char alphanum[] =
"0123456789"
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"abcdefghijklmnopqrstuvwxyz";
std::string tmp_s;
tmp_s.reserve(len);
for (int i = 0; i < len; ++i) {
tmp_s += alphanum[rand() % (sizeof(alphanum) - 1)];
}
return tmp_s;
}
template <typename T>
sitk::Image make_image(
unsigned int width,
unsigned int height,
T* image,
sitk::PixelIDValueEnum id
) {
sitk::Image im(width, height, id);
if (id == sitk::PixelIDValueEnum::sitkUInt8) {
uint8_t* b = im.GetBufferAsUInt8();
memcpy(b, image, width * height);
} else if (id == sitk::PixelIDValueEnum::sitkInt8) {
int8_t* b = im.GetBufferAsInt8();
memcpy(b, image, width * height);
} else if (id == sitk::PixelIDValueEnum::sitkUInt16) {
uint16_t* b = im.GetBufferAsUInt16();
memcpy(b, image, width * height * 2);
} else if (id == sitk::PixelIDValueEnum::sitkInt16) {
int16_t* b = im.GetBufferAsInt16();
memcpy(b, image, width * height * 2);
} else if (id == sitk::PixelIDValueEnum::sitkUInt32) {
uint32_t* b = im.GetBufferAsUInt32();
memcpy(b, image, width * height * 4);
} else if (id == sitk::PixelIDValueEnum::sitkInt32) {
int32_t* b = im.GetBufferAsInt32();
memcpy(b, image, width * height * 4);
} else if (id == sitk::PixelIDValueEnum::sitkUInt64) {
uint64_t* b = im.GetBufferAsUInt64();
memcpy(b, image, width * height * 8);
} else if (id == sitk::PixelIDValueEnum::sitkInt64) {
int64_t* b = im.GetBufferAsInt64();
memcpy(b, image, width * height * 8);
} else if (id == sitk::PixelIDValueEnum::sitkFloat32) {
float* b = im.GetBufferAsFloat();
memcpy(b, image, width * height * 4);
} else if (id == sitk::PixelIDValueEnum::sitkFloat64) {
double* b = im.GetBufferAsDouble();
memcpy(b, image, width * height * 8);
}
return im;
}
sitk::Image
interp(
double* transform,
double* origin,
sitk::Image image,
bool bspline_or_nn
) {
try {
vector<double> matrix = {transform[0], transform[1], transform[2], transform[3]};
vector<double> translation = {transform[4], transform[5]};
vector<double> ori = {origin[0], origin[1]};
sitk::AffineTransform t(matrix, translation, ori);
sitk::InterpolatorEnum interpolator = (bspline_or_nn == false) ? sitk::sitkBSpline : sitk::sitkNearestNeighbor;
image = sitk::Resample(image, t, interpolator);
return image;
} catch (const std::exception &exc) {
cerr << exc.what();
return image;
}
}
extern "C" void
interp_u8(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
uint8_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkUInt8);
im = interp(transform, origin, im, bspline_or_nn);
uint8_t* c = im.GetBufferAsUInt8();
memcpy(*image, c, width * height);
}
extern "C" void
interp_i8(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
int8_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkInt8);
im = interp(transform, origin, im, bspline_or_nn);
int8_t* c = im.GetBufferAsInt8();
memcpy(*image, c, width * height);
}
extern "C" void
interp_u16(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
uint16_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkUInt16);
im = interp(transform, origin, im, bspline_or_nn);
uint16_t* c = im.GetBufferAsUInt16();
memcpy(*image, c, width * height * 2);
}
extern "C" void
interp_i16(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
int16_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkInt16);
im = interp(transform, origin, im, bspline_or_nn);
int16_t* c = im.GetBufferAsInt16();
memcpy(*image, c, width * height * 2);
}
extern "C" void
interp_u32(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
uint32_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkUInt32);
im = interp(transform, origin, im, bspline_or_nn);
uint32_t* c = im.GetBufferAsUInt32();
memcpy(*image, c, width * height * 4);
}
extern "C" void
interp_i32(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
int32_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkInt32);
im = interp(transform, origin, im, bspline_or_nn);
int32_t* c = im.GetBufferAsInt32();
memcpy(*image, c, width * height * 4);
}
extern "C" void
interp_u64(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
uint64_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkUInt64);
im = interp(transform, origin, im, bspline_or_nn);
uint64_t* c = im.GetBufferAsUInt64();
memcpy(*image, c, width * height * 8);
}
extern "C" void
interp_i64(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
int64_t** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkInt64);
im = interp(transform, origin, im, bspline_or_nn);
int64_t* c = im.GetBufferAsInt64();
memcpy(*image, c, width * height * 8);
}
extern "C" void
interp_f32(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
float** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkFloat32);
im = interp(transform, origin, im, bspline_or_nn);
float* c = im.GetBufferAsFloat();
memcpy(*image, c, width * height * 4);
}
extern "C" void
interp_f64(
unsigned int width,
unsigned int height,
double* transform,
double* origin,
double** image,
bool bspline_or_nn
) {
sitk::Image im = make_image(width, height, *image, sitk::PixelIDValueEnum::sitkFloat64);
im = interp(transform, origin, im, bspline_or_nn);
double* c = im.GetBufferAsDouble();
memcpy(*image, c, width * height * 8);
}
void
reg2(
sitk::Image fixed,
sitk::Image moving,
bool t_or_a,
double** transform
) {
try {
string kind = (t_or_a == false) ? "translation" : "affine";
sitk::ImageRegistrationMethod R;
R.SetMetricAsMattesMutualInformation();
const double maxStep = 4.0;
const double minStep = 0.01;
const unsigned int numberOfIterations = 200;
const double relaxationFactor = 0.5;
// R.SetOptimizerAsLBFGSB(maxStep, minStep, numberOfIterations, relaxationFactor);
R.SetOptimizerAsRegularStepGradientDescent(maxStep, minStep, numberOfIterations, relaxationFactor);
// R.SetOptimizerAsLBFGS2();
vector<double> matrix = {1.0, 0.0, 0.0, 1.0};
vector<double> translation = {0., 0.};
vector<double> origin = {399.5, 299.5};
R.SetInitialTransform(sitk::AffineTransform(matrix, translation, origin));
R.SetInterpolator(sitk::sitkBSpline);
sitk::Transform outTx = R.Execute(fixed, moving);
vector<double> t = outTx.GetParameters();
for (int i = 0; i < t.size(); i++) {
cout << t[i] << " ";
(*transform)[i] = t[i];
}
} catch (const std::exception &exc) {
cerr << exc.what();
}
}
void
reg(
sitk::Image fixed,
sitk::Image moving,
bool t_or_a,
double** transform
) {
try {
string kind = (t_or_a == false) ? "translation" : "affine";
// std::filesystem::path output_path = std::filesystem::temp_directory_path() / gen_random(12);
// std::filesystem::create_directory(output_path);
std::filesystem::path output_path = std::filesystem::temp_directory_path();
sitk::ElastixImageFilter tfilter = sitk::ElastixImageFilter();
tfilter.LogToConsoleOff();
tfilter.LogToFileOff();
tfilter.SetLogToFile(false);
tfilter.SetFixedImage(fixed);
tfilter.SetMovingImage(moving);
tfilter.SetParameterMap(sitk::GetDefaultParameterMap(kind));
tfilter.SetParameter("WriteResultImage", "false");
tfilter.SetOutputDirectory(output_path);
tfilter.Execute();
sitk::ElastixImageFilter::ParameterMapType parameter_map = tfilter.GetTransformParameterMap(0);
for (sitk::ElastixImageFilter::ParameterMapType::iterator parameter = parameter_map.begin(); parameter != parameter_map.end(); ++parameter) {
if (parameter->first == "TransformParameters") {
vector<string> tp = parameter->second;
if (t_or_a == true) {
for (int j = 0; j < tp.size(); j++) {
(*transform)[j] = stod(tp[j]);
}
} else {
(*transform)[0] = 1.0;
(*transform)[1] = 0.0;
(*transform)[2] = 0.0;
(*transform)[3] = 1.0;
for (int j = 0; j < tp.size(); j++) {
(*transform)[j + 4] = stod(tp[j]);
}
}
break;
}
}
} catch (const std::exception &exc) {
cerr << exc.what();
}
}
extern "C" void
register_u8(
unsigned int width,
unsigned int height,
uint8_t** fixed_arr,
uint8_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkUInt8;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_i8(
unsigned int width,
unsigned int height,
int8_t** fixed_arr,
int8_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkInt8;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_u16(
unsigned int width,
unsigned int height,
uint16_t** fixed_arr,
uint16_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkUInt16;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_i16(
unsigned int width,
unsigned int height,
int16_t** fixed_arr,
int16_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkInt16;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_u32(
unsigned int width,
unsigned int height,
uint32_t** fixed_arr,
uint32_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkUInt32;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_i32(
unsigned int width,
unsigned int height,
int32_t** fixed_arr,
int32_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkInt32;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_u64(
unsigned int width,
unsigned int height,
uint64_t** fixed_arr,
uint64_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkUInt64;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_i64(
unsigned int width,
unsigned int height,
int64_t** fixed_arr,
int64_t** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkInt64;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_f32(
unsigned int width,
unsigned int height,
float** fixed_arr,
float** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkFloat32;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}
extern "C" void
register_f64(
unsigned int width,
unsigned int height,
double** fixed_arr,
double** moving_arr,
bool t_or_a,
double** transform
) {
sitk::PixelIDValueEnum id = sitk::PixelIDValueEnum::sitkFloat64;
sitk::Image fixed = make_image(width, height, *fixed_arr, id);
sitk::Image moving = make_image(width, height, *moving_arr, id);
reg(fixed, moving, t_or_a, transform);
}

58
cxx/ffi_extra.cxx Normal file
View File

@@ -0,0 +1,58 @@
#include "ffi_extra.h"
#include <map>
#include <memory>
#include <utility>
#include <vector>
#include <string>
std::unique_ptr<ParameterMap> get_transform_parameter_map(itk::simple::ElastixImageFilter& tfilter, const unsigned int index) {
itk::simple::ElastixImageFilter::ParameterMapType parameter_map = tfilter.GetTransformParameterMap(index);
ParameterMap transform_parameter_map = ParameterMap(parameter_map);
auto ptr = std::make_unique<ParameterMap>(parameter_map);
return ptr;
}
std::unique_ptr<ParameterMap> get_default_parameter_map(const std::string& kind) {
itk::simple::ElastixImageFilter::ParameterMapType parameter_map = itk::simple::GetDefaultParameterMap(kind);
ParameterMap default_parameter_map = ParameterMap(parameter_map);
auto ptr = std::make_unique<ParameterMap>(default_parameter_map);
return ptr;
}
void set_parameter_map(itk::simple::ElastixImageFilter& tfilter, const std::unique_ptr<ParameterMap>& parameter_map) {
tfilter.SetParameterMap(parameter_map->parameters);
};
ParameterMap::ParameterMap() = default;
ParameterMap::~ParameterMap() = default;
ParameterMap::ParameterMap(itk::simple::ElastixImageFilter::ParameterMapType parameter_map) {
this->parameters = std::move(parameter_map);
}
std::unique_ptr<ParameterMap> new_parameter_map() {
ParameterMap parameter_map;
auto ptr = std::make_unique<ParameterMap>(parameter_map);
return ptr;
}
void ParameterMap::insert(const std::string &key, const std::vector<std::string> &value) {
this->parameters.emplace(key, value);
}
std::unique_ptr<std::vector<std::string>> ParameterMap::keys() const {
std::vector<std::string> keys;
for (auto &[key, value] : this->parameters) {
keys.push_back(key);
}
auto ptr = std::make_unique<std::vector<std::string>>(keys);
return ptr;
}
std::unique_ptr<std::vector<std::string>> ParameterMap::get(const std::string& key) const {
std::vector<std::string> value = this->parameters.at(key);
auto ptr = std::make_unique<std::vector<std::string>>(value);
return ptr;
}

20
cxx/ffi_extra.h Normal file
View File

@@ -0,0 +1,20 @@
#pragma once
#include <sitkElastixImageFilter.h>
#include <memory>
#include <vector>
class ParameterMap {
public:
itk::simple::ElastixImageFilter::ParameterMapType parameters;
ParameterMap();
~ParameterMap();
explicit ParameterMap(itk::simple::ElastixImageFilter::ParameterMapType parameter_map);
void insert(const std::string &key, const std::vector<std::string> &value);
[[nodiscard]] std::unique_ptr<std::vector<std::string>> get(const std::string& key) const;
[[nodiscard]] std::unique_ptr<std::vector<std::string>> keys() const;
};
std::unique_ptr<ParameterMap> new_parameter_map();
std::unique_ptr<ParameterMap> get_default_parameter_map(const std::string& kind);
std::unique_ptr<ParameterMap> get_transform_parameter_map(itk::simple::ElastixImageFilter& tfilter, unsigned int index);
void set_parameter_map(itk::simple::ElastixImageFilter& tfilter, const std::unique_ptr<ParameterMap>& parameter_map);

14
merged_lib/CMakeLists.txt Normal file
View File

@@ -0,0 +1,14 @@
cmake_minimum_required(VERSION 3.14.0)
project(merged_lib)
file(GLOB_RECURSE LIB_SITK ${CMAKE_BINARY_DIR}/../../build/SimpleITK-build/**/*.o)
file(GLOB_RECURSE LIB_ITK ${CMAKE_BINARY_DIR}/../../build/ITK-build/**/*.o)
file(GLOB_RECURSE LIB_ELASTIX ${CMAKE_BINARY_DIR}/../../build/Elastix-build/**/*.o)
add_library( sitk ${LIB_SITK} )
add_library( itk ${LIB_ITK} )
add_library( elastix ${LIB_ELASTIX} )
set_target_properties(sitk PROPERTIES LINKER_LANGUAGE CXX)
set_target_properties(itk PROPERTIES LINKER_LANGUAGE CXX)
set_target_properties(elastix PROPERTIES LINKER_LANGUAGE CXX)

7130
src/autocxx-ffi-default-gen.rs Normal file
View File

File diff suppressed because it is too large Load Diff

564
src/lib.rs
View File

@@ -1,485 +1,101 @@
mod sys; //! This crate does two things:
//! - find an affine transform or translation that transforms one image into the other
//! - use bspline or nearest neighbor interpolation to apply a transformation to an image
//!
//! To do this, [SimpleITK](https://github.com/SimpleITK/SimpleITK), which is written in
//! C++, is used. An adapter library is created using [autocxx](https://crates.io/crates/autocxx)
//! to expose the required functionality in SimpleITK. Because of this, compilation of this crate
//! requires quite some time, several GB of memory, up to 50 GB of hard disk space, as well as
//! cmake, a C++ compiler, llvm and git. Use at your own risk!
//!
//! # Examples
//! ## Registration
//! ```
//! use ndarray::Array2;
//! use sitk_registration_sys::registration::{AffineTransform, julia_image};
//!
//! let j = julia_image(0f32, 0f32).unwrap();
//! let shape = j.shape();
//! let origin = [
//! ((shape[1] - 1) as f64) / 2f64,
//! ((shape[0] - 1) as f64) / 2f64,
//! ];
//! let s = AffineTransform::new([1.2, 0., 0., 1., 5., 7.], origin, [shape[0], shape[1]]);
//! let k: Array2<_> = s.transform_image_bspline(j.view()).unwrap().into();
//! let t = AffineTransform::register_affine(j.view(), k.view()).unwrap().inverse().unwrap();
//! let d = (t.matrix() - s.matrix()).powi(2).sum();
//! assert!(d < 0.025, "d: {}, t: {:?}", d, t.parameters);
//! ```
//!
//! ## Interpolation
//! ```
//! use ndarray::Array2;
//! use sitk_registration_sys::registration::{AffineTransform, julia_image};
//!
//! let j = julia_image(-120f32, 10f32).unwrap();
//! let k = julia_image(0f32, 0f32).unwrap();
//! let shape = j.shape();
//! let origin = [
//! ((shape[1] - 1) as f64) / 2f64,
//! ((shape[0] - 1) as f64) / 2f64,
//! ];
//! let transform = AffineTransform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
//! let n: Array2<_> = transform.transform_image_bspline(j.view()).unwrap().into();
//! 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);
//! ```
use crate::sys::{interp, register}; extern crate link_cplusplus;
use anyhow::{Result, anyhow}; pub mod registration;
use ndarray::{Array2, ArrayView2, AsArray, Ix2, 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;
/// a trait marking number types that can be used in sitk: /// The bindings generated by [autocxx](https://google.github.io/autocxx).
/// (u/i)(8/16/32/64), (u/i)size, f(32/64) /// Everything in here is unsafe.
pub trait PixelType: Clone { pub use ffi::itk::simple;
const PT: u8;
use autocxx::prelude::*;
include_cpp! {
#include "sitkAdditionalProcedures.h"
#include "sitkAffineTransform.h"
#include "sitkElastixImageFilter.h"
#include "sitkInterpolator.h"
#include "sitkImage.h"
#include "sitkTransform.h"
safety!(unsafe)
generate!("itk::simple::AffineTransform")
generate!("itk::simple::ElastixImageFilter")
generate!("itk::simple::InterpolatorEnum")
generate!("itk::simple::Image")
generate!("itk::simple::Resample")
generate!("itk::simple::Transform")
opaque!("itk::simple::ElastixImageFilter") // deleted constructor
} }
macro_rules! sitk_impl { /// Manually generated bindings for some things that autocxx does not understand.
($($T:ty: $sitk:expr $(,)?)*) => { #[cxx::bridge]
$( pub mod ffi_extra {
impl PixelType for $T { unsafe extern "C++" {
const PT: u8 = $sitk; include!("ffi_extra.h");
}
)*
};
}
sitk_impl! { #[namespace = "itk::simple"]
u8: 1, type ElastixImageFilter = crate::ffi::itk::simple::ElastixImageFilter;
i8: 2, type ParameterMap;
u16: 3,
i16: 4,
u32: 5,
i32: 6,
u64: 7,
i64: 8,
f32: 9,
f64: 10,
}
#[cfg(target_pointer_width = "64")] fn new_parameter_map() -> UniquePtr<ParameterMap>;
sitk_impl!(usize: 7); fn insert(self: Pin<&mut ParameterMap>, key: &CxxString, value: &CxxVector<CxxString>);
#[cfg(target_pointer_width = "32")] fn keys(self: &ParameterMap) -> UniquePtr<CxxVector<CxxString>>;
sitk_impl!(usize: 5); fn get(self: &ParameterMap, key: &CxxString) -> UniquePtr<CxxVector<CxxString>>;
#[cfg(target_pointer_width = "64")]
sitk_impl!(isize: 8);
#[cfg(target_pointer_width = "32")]
sitk_impl!(isize: 6);
#[derive(Clone, Debug, Deserialize, Serialize)] fn get_transform_parameter_map(
pub struct Transform { tfilter: &mut ElastixImageFilter,
pub parameters: [f64; 6], index: u32,
pub dparameters: [f64; 6], ) -> UniquePtr<ParameterMap>;
pub origin: [f64; 2],
pub shape: [usize; 2],
}
impl Mul for Transform { fn get_default_parameter_map(kind: &CxxString) -> UniquePtr<ParameterMap>;
type Output = Transform;
#[allow(clippy::suspicious_arithmetic_impl)] fn set_parameter_map(
fn mul(self, other: Transform) -> Transform { tfilter: &mut ElastixImageFilter,
let m = self.matrix().dot(&other.matrix()); parameter_map: &UniquePtr<ParameterMap>,
let dm = self.dmatrix().dot(&other.matrix()) + self.matrix().dot(&other.dmatrix()); );
Transform {
parameters: [
m[[0, 0]],
m[[0, 1]],
m[[1, 0]],
m[[1, 1]],
m[[2, 0]],
m[[2, 1]],
],
dparameters: [
dm[[0, 0]],
dm[[0, 1]],
dm[[1, 0]],
dm[[1, 1]],
dm[[2, 0]],
dm[[2, 1]],
],
origin: self.origin,
shape: self.shape,
}
}
}
impl PartialEq<Self> for Transform {
fn eq(&self, other: &Self) -> bool {
self.parameters == other.parameters
&& self.dparameters == other.dparameters
&& self.origin == other.origin
&& self.shape == other.shape
}
}
impl Eq for Transform {}
impl Transform {
/// parameters: flat 2x2 part of matrix, translation; origin: center of rotation
pub fn new(parameters: [f64; 6], origin: [f64; 2], shape: [usize; 2]) -> Self {
Self {
parameters,
dparameters: [0f64; 6],
origin,
shape,
}
}
/// find the affine transform which transforms moving into fixed
pub fn register_affine<'a, A, T>(fixed: A, moving: A) -> Result<Transform>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
let (parameters, origin, shape) = register(fixed, moving, true)?;
Ok(Transform {
parameters,
dparameters: [0f64; 6],
origin,
shape,
})
}
/// find the translation which transforms moving into fixed
pub fn register_translation<'a, A, T>(fixed: A, moving: A) -> Result<Transform>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
let (parameters, origin, shape) = register(fixed, moving, false)?;
Ok(Transform {
parameters,
dparameters: [0f64; 6],
origin,
shape,
})
}
/// create a transform from a xy translation
pub fn from_translation(translation: [f64; 2]) -> Self {
Transform {
parameters: [1f64, 0f64, 0f64, 1f64, translation[0], translation[1]],
dparameters: [0f64; 6],
origin: [0f64; 2],
shape: [0usize; 2],
}
}
/// read a transform from a file
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, path: PathBuf) -> Result<()> {
let mut file = std::fs::OpenOptions::new()
.create(true)
.write(true)
.truncate(true)
.open(path)?;
to_writer(&mut file, self)?;
Ok(())
}
/// true if transform does nothing
pub fn is_unity(&self) -> bool {
self.parameters == [1f64, 0f64, 0f64, 1f64, 0f64, 0f64]
}
/// transform an image using nearest neighbor interpolation
pub fn transform_image_bspline<'a, A, T>(&self, image: A) -> Result<Array2<T>>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
interp(self.parameters, self.origin, image, false)
}
/// transform an image using bspline interpolation
pub fn transform_image_nearest_neighbor<'a, A, T>(&self, image: A) -> Result<Array2<T>>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
interp(self.parameters, self.origin, image, true)
}
/// get coordinates resulting from transforming input coordinates
pub fn transform_coordinates<'a, A, T>(&self, coordinates: A) -> Result<Array2<f64>>
where
T: 'a + Clone + Into<f64>,
A: AsArray<'a, T, Ix2>,
{
let coordinates = coordinates.into();
let s = coordinates.shape();
if s[1] != 2 {
return Err(anyhow!("coordinates must have two columns"));
}
let m = self.matrix();
let mut res = Array2::zeros([s[0], s[1]]);
for i in 0..s[0] {
let a = array![
coordinates[[i, 0]].clone().into(),
coordinates[[i, 1]].clone().into(),
1f64
]
.to_owned();
let b = m.dot(&a);
res.slice_mut(s![i, ..]).assign(&b.slice(s![..2]));
}
Ok(res)
}
/// get the matrix defining the transform
pub fn matrix(&self) -> Array2<f64> {
Array2::from_shape_vec(
(3, 3),
vec![
self.parameters[0],
self.parameters[1],
self.parameters[4],
self.parameters[2],
self.parameters[3],
self.parameters[5],
0f64,
0f64,
1f64,
],
)
.unwrap()
}
/// get the matrix describing the error of the transform
pub fn dmatrix(&self) -> Array2<f64> {
Array2::from_shape_vec(
(3, 3),
vec![
self.dparameters[0],
self.dparameters[1],
self.dparameters[4],
self.dparameters[2],
self.dparameters[3],
self.dparameters[5],
0f64,
0f64,
1f64,
],
)
.unwrap()
}
/// get the inverse transform
pub fn inverse(&self) -> Result<Transform> {
fn det(a: ArrayView2<f64>) -> f64 {
(a[[0, 0]] * a[[1, 1]]) - (a[[0, 1]] * a[[1, 0]])
}
let m = self.matrix();
let d = det(m.slice(s![..2, ..2]));
if d == 0f64 {
return Err(anyhow!("transform matrix is not invertible"));
}
let parameters = [
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 {
parameters,
dparameters: [0f64; 6],
origin: self.origin,
shape: self.shape,
})
}
/// adapt the transform to a new origin and shape
pub fn adapt(&mut self, origin: [f64; 2], shape: [usize; 2]) {
self.origin = [
origin[0] + (((self.shape[0] - shape[0]) as f64) / 2f64),
origin[1] + (((self.shape[1] - shape[1]) as f64) / 2f64),
];
self.shape = shape;
}
}
#[cfg(test)]
mod tests {
use super::*;
use anyhow::Result;
use ndarray::Array2;
use num::Complex;
use tempfile::NamedTempFile;
/// An example of generating julia fractals.
fn julia_image(shift_x: f32, shift_y: f32) -> Result<Array2<u8>> {
let imgx = 800;
let imgy = 600;
let scalex = 3.0 / imgx as f32;
let scaley = 3.0 / imgy as f32;
let mut im = Array2::<u8>::zeros((imgy, imgx));
for x in 0..imgx {
for y in 0..imgy {
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(cy, cx);
let mut i = 0;
while i < 255 && z.norm() <= 2.0 {
z = z * z + c;
i += 1;
}
im[[y, x]] = i as u8;
}
}
Ok(im)
}
#[test]
fn test_serialization() -> Result<()> {
let file = NamedTempFile::new()?;
let t = Transform::new([1.2, 0.3, -0.4, 0.9, 10.2, -9.5], [59.5, 49.5], [120, 100]);
t.to_file(file.path().to_path_buf())?;
let s = Transform::from_file(file.path().to_path_buf())?;
assert_eq!(s, t);
Ok(())
}
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(())
}
)*
}
}
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 j0 = j.clone();
let k0 = k.clone();
let transform = Transform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
// make sure j & k weren't mutated
assert!(j.iter().zip(j0.iter()).map(|(a, b)| a == b).all(|x| x));
assert!(k.iter().zip(k0.iter()).map(|(a, b)| a == b).all(|x| x));
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 j0 = j.clone();
let k0 = k.clone();
let t = Transform::register_translation(j.view(), k.view())?;
// make sure j & k weren't mutated
assert!(j.iter().zip(j0.iter()).map(|(a, b)| a == b).all(|x| x));
assert!(k.iter().zip(k0.iter()).map(|(a, b)| a == b).all(|x| x));
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,
} }
} }

645
src/registration.rs Normal file
View File

@@ -0,0 +1,645 @@
//! Some structs and methods to make working with registration and interpolation methods in
//! SimpleITK more Rust friendly.
use crate::simple;
use anyhow::{Result, anyhow};
use autocxx::prelude::*;
use cxx::{CxxVector, UniquePtr, let_cxx_string};
use ndarray::{Array2, ArrayView2, AsArray, Ix2, array, s};
use serde::{Deserialize, Serialize};
use serde_yaml::{from_reader, to_writer};
use std::fs::File;
use std::marker::PhantomData;
use std::ops::{Deref, DerefMut, Mul};
use std::path::PathBuf;
use num::Complex;
use tempfile::tempdir;
/// a trait marking number types that can be used in sitk:
/// (u/i)(8/16/32/64), (u/i)size, f(32/64)
pub trait PixelType: Clone {
const PT: simple::PixelIDValueEnum;
}
macro_rules! pixel_type_impl {
($($T:ty: $sitk:expr $(,)?)*) => {
$(
impl PixelType for $T {
const PT: simple::PixelIDValueEnum = $sitk;
}
)*
};
}
pixel_type_impl! {
u8: simple::PixelIDValueEnum::sitkUInt8,
i8: simple::PixelIDValueEnum::sitkInt8,
u16: simple::PixelIDValueEnum::sitkUInt16,
i16: simple::PixelIDValueEnum::sitkInt16,
u32: simple::PixelIDValueEnum::sitkUInt32,
i32: simple::PixelIDValueEnum::sitkInt32,
u64: simple::PixelIDValueEnum::sitkUInt64,
i64: simple::PixelIDValueEnum::sitkInt64,
f32: simple::PixelIDValueEnum::sitkFloat32,
f64: simple::PixelIDValueEnum::sitkFloat64,
}
#[cfg(target_pointer_width = "64")]
pixel_type_impl!(usize: simple::PixelIDValueEnum::sitkUInt64);
#[cfg(target_pointer_width = "32")]
pixel_type_impl!(usize: simple::PixelIDValueEnum::sitkUInt32);
#[cfg(target_pointer_width = "64")]
pixel_type_impl!(isize: simple::PixelIDValueEnum::sitkInt64);
#[cfg(target_pointer_width = "32")]
pixel_type_impl!(isize: simple::PixelIDValueEnum::sitkInt32);
/// Struct holding a pointer to an image
pub struct Image<T: PixelType> {
image: UniquePtr<simple::Image>,
pixel_type: PhantomData<T>,
}
impl<T: PixelType> Deref for Image<T> {
type Target = UniquePtr<simple::Image>;
fn deref(&self) -> &Self::Target {
&self.image
}
}
impl<T: PixelType> Image<T> {
/// encapsulate an itk::simple::Image
pub fn new(image: UniquePtr<simple::Image>) -> Self {
Self {
image,
pixel_type: PhantomData,
}
}
/// take an ndarray Array2 and turn it into a SimpleITK image
pub fn from_array<'a, A>(array: A) -> Self
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
let array = array.into();
let shape = array.shape();
let width = (shape[1] as u32).into();
let height = (shape[0] as u32).into();
let mut image = simple::Image::new3(width, height, T::PT).within_unique_ptr();
image.pin_mut().MakeUnique();
let buffer = image.pin_mut().GetBufferAsVoid();
unsafe { std::ptr::copy(array.as_ptr(), buffer as *mut T, shape[0] * shape[1]) };
Self {
image,
pixel_type: PhantomData,
}
}
/// return as an ndarray Array2
pub fn as_array(&self) -> Array2<T> {
let width = u32::from(self.image.GetWidth()) as usize;
let height = u32::from(self.image.GetHeight()) as usize;
let mut array = Array2::<T>::uninit((height, width));
let buffer = self.image.GetBufferAsVoid1();
unsafe {
std::ptr::copy(
buffer as *const T,
array.as_mut_ptr() as *mut T,
width * height,
);
array.assume_init()
}
}
}
impl<'a, A, T> From<A> for Image<T>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
fn from(value: A) -> Self {
Self::from_array(value.into())
}
}
impl<T> From<Image<T>> for Array2<T>
where
T: PixelType,
{
fn from(value: Image<T>) -> Self {
value.as_array()
}
}
/// a struct describing the transform
#[derive(Clone, Debug, Deserialize, Serialize, PartialEq)]
pub struct AffineTransform {
/// flattened 2x2 rotation matrix + translation
pub parameters: [f64; 6],
/// error / significance on parameters
pub dparameters: [f64; 6],
/// the point about which rotations are performed
pub origin: [f64; 2],
/// the shape of images for which this transform is meant
pub shape: [usize; 2],
}
impl Mul for AffineTransform {
type Output = AffineTransform;
#[allow(clippy::suspicious_arithmetic_impl)]
fn mul(self, other: AffineTransform) -> AffineTransform {
let m = self.matrix().dot(&other.matrix());
let dm = self.dmatrix().dot(&other.matrix()) + self.matrix().dot(&other.dmatrix());
AffineTransform {
parameters: [
m[[0, 0]],
m[[0, 1]],
m[[1, 0]],
m[[1, 1]],
m[[2, 0]],
m[[2, 1]],
],
dparameters: [
dm[[0, 0]],
dm[[0, 1]],
dm[[1, 0]],
dm[[1, 1]],
dm[[2, 0]],
dm[[2, 1]],
],
origin: self.origin,
shape: self.shape,
}
}
}
impl Eq for AffineTransform {}
impl AffineTransform {
/// parameters: flat 2x2 part of matrix, translation; origin: center of rotation
pub fn new(parameters: [f64; 6], origin: [f64; 2], shape: [usize; 2]) -> Self {
Self {
parameters,
dparameters: [0f64; 6],
origin,
shape,
}
}
/// find the affine transform which transforms moving into fixed
pub fn register_affine<F, M, T>(fixed: F, moving: M) -> Result<AffineTransform>
where
F: Into<Image<T>>,
M: Into<Image<T>>,
T: PixelType,
{
Self::register(fixed, moving, true)
}
/// find the translation which transforms moving into fixed
pub fn register_translation<F, M, T>(fixed: F, moving: M) -> Result<AffineTransform>
where
F: Into<Image<T>>,
M: Into<Image<T>>,
T: PixelType,
{
Self::register(fixed, moving, false)
}
/// find the transform which transforms moving into fixed
pub fn register<F, M, T>(fixed: F, moving: M, affine: bool) -> Result<AffineTransform>
where
F: Into<Image<T>>,
M: Into<Image<T>>,
T: PixelType,
{
let tmp_folder = tempdir()?;
let fixed = fixed.into();
let moving = moving.into();
let width = u32::from(fixed.GetWidth()) as usize;
let height = u32::from(fixed.GetHeight()) as usize;
let_cxx_string!(transform_name = if affine { "affine" } else { "translation" });
let parameter_map = crate::ffi_extra::get_default_parameter_map(&transform_name);
let mut tfilter = simple::ElastixImageFilter::new().within_box();
tfilter.as_mut().LogToConsoleOff();
tfilter.as_mut().LogToFileOff();
tfilter.as_mut().SetLogToFile(false);
tfilter.as_mut().SetFixedImage(&fixed);
tfilter.as_mut().SetMovingImage(&moving);
crate::ffi_extra::set_parameter_map(&mut tfilter, &parameter_map);
tfilter.as_mut().SetParameter("WriteResultImage", "False");
tfilter
.as_mut()
.SetOutputDirectory(tmp_folder.path().display().to_string());
let _ = tfilter.as_mut().Execute().within_unique_ptr();
let_cxx_string!(tp = "TransformParameters");
let p = crate::ffi_extra::get_transform_parameter_map(tfilter.as_mut().deref_mut(), 0)
.get(&tp)
.iter()
.map(|i| i.to_string_lossy().parse::<f64>())
.collect::<Result<Vec<f64>, _>>()?;
let parameters = if affine {
[p[0], p[1], p[2], p[3], p[4], p[5]]
} else {
[1.0, 0.0, 0.0, 1.0, p[0], p[1]]
};
let origin = [((height - 1) as f64) / 2.0, ((width - 1) as f64) / 2.0];
let shape = [height, width];
Ok(AffineTransform::new(parameters, origin, shape))
}
/// create a transform from a xy translation
pub fn from_translation(translation: [f64; 2]) -> Self {
AffineTransform {
parameters: [1f64, 0f64, 0f64, 1f64, translation[0], translation[1]],
dparameters: [0f64; 6],
origin: [0f64; 2],
shape: [0usize; 2],
}
}
/// read a transform from a file
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, path: PathBuf) -> Result<()> {
let mut file = std::fs::OpenOptions::new()
.create(true)
.write(true)
.truncate(true)
.open(path)?;
to_writer(&mut file, self)?;
Ok(())
}
/// true if transform does nothing
pub fn is_unity(&self) -> bool {
self.parameters == [1f64, 0f64, 0f64, 1f64, 0f64, 0f64]
}
/// transform an image using nearest neighbor interpolation
pub fn transform_image_bspline<I, T>(&self, image: I) -> Result<Image<T>>
where
I: Into<Image<T>>,
T: PixelType,
{
Self::interp(self, image, false)
}
/// transform an image using bspline interpolation
pub fn transform_image_nearest_neighbor<I, T>(&self, image: I) -> Result<Image<T>>
where
I: Into<Image<T>>,
T: PixelType,
{
Self::interp(self, image, true)
}
/// transform an image
pub fn interp<I, T>(&self, image: I, nearest_neighbor: bool) -> Result<Image<T>>
where
I: Into<Image<T>>,
T: PixelType,
{
let image = image.into();
let width = u32::from(image.GetWidth()) as usize;
let height = u32::from(image.GetHeight()) as usize;
let origin = [((width - 1) as f64) / 2f64, ((height - 1) as f64) / 2f64];
let p = self.parameters;
let matrix = cxx_vector([p[0], p[1], p[2], p[3]]);
let translation = cxx_vector([p[4], p[5]]);
let fixed_center = cxx_vector(origin);
let affine_transform =
simple::AffineTransform::new3(&matrix, &translation, &fixed_center).within_unique_ptr();
let transform = <_ as AsRef<simple::Transform>>::as_ref(affine_transform.as_ref().unwrap());
let interpolator = if nearest_neighbor {
simple::InterpolatorEnum::sitkBSpline
} else {
simple::InterpolatorEnum::sitkNearestNeighbor
};
Ok(Image::<T>::new(
simple::Resample(
&image,
transform,
interpolator,
0.0,
simple::PixelIDValueEnum::sitkUnknown,
nearest_neighbor,
)
.within_unique_ptr(),
))
}
/// get coordinates resulting from transforming input coordinates, coordinates must have two
/// columns: x & y
pub fn transform_coordinates<'a, A, T>(&self, coordinates: A) -> Result<Array2<f64>>
where
T: 'a + Clone + Into<f64>,
A: AsArray<'a, T, Ix2>,
{
let coordinates = coordinates.into();
let s = coordinates.shape();
if s[1] != 2 {
return Err(anyhow!("coordinates must have two columns"));
}
let m = self.matrix();
let mut res = Array2::zeros([s[0], s[1]]);
for i in 0..s[0] {
let a = array![
coordinates[[i, 0]].clone().into(),
coordinates[[i, 1]].clone().into(),
1f64
]
.to_owned();
let b = m.dot(&a);
res.slice_mut(s![i, ..]).assign(&b.slice(s![..2]));
}
Ok(res)
}
/// get the matrix defining the transform
pub fn matrix(&self) -> Array2<f64> {
Array2::from_shape_vec(
(3, 3),
vec![
self.parameters[0],
self.parameters[1],
self.parameters[4],
self.parameters[2],
self.parameters[3],
self.parameters[5],
0f64,
0f64,
1f64,
],
)
.unwrap()
}
/// get the matrix describing the error of the transform
pub fn dmatrix(&self) -> Array2<f64> {
Array2::from_shape_vec(
(3, 3),
vec![
self.dparameters[0],
self.dparameters[1],
self.dparameters[4],
self.dparameters[2],
self.dparameters[3],
self.dparameters[5],
0f64,
0f64,
1f64,
],
)
.unwrap()
}
/// get the inverse transform
pub fn inverse(&self) -> Result<AffineTransform> {
fn det(a: ArrayView2<f64>) -> f64 {
(a[[0, 0]] * a[[1, 1]]) - (a[[0, 1]] * a[[1, 0]])
}
let m = self.matrix();
let d = det(m.slice(s![..2, ..2]));
if d == 0f64 {
return Err(anyhow!("transform matrix is not invertible"));
}
let parameters = [
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(AffineTransform {
parameters,
dparameters: [0f64; 6],
origin: self.origin,
shape: self.shape,
})
}
/// adapt the transform to a new origin and shape
pub fn adapt(&mut self, origin: [f64; 2], shape: [usize; 2]) {
self.origin = [
origin[0] + (((self.shape[0] - shape[0]) as f64) / 2f64),
origin[1] + (((self.shape[1] - shape[1]) as f64) / 2f64),
];
self.shape = shape;
}
}
/// conveniently collect an iterator into a CxxVector
pub fn cxx_vector<T, I>(vec: I) -> UniquePtr<CxxVector<T>>
where
I: IntoIterator<Item = T>,
T: cxx::vector::VectorElement + cxx::ExternType<Kind = cxx::kind::Trivial>,
{
let mut v = CxxVector::new();
v.pin_mut().extend(vec);
v
}
/// An example of generating julia fractals, for testing purposes.
pub fn julia_image(shift_x: f32, shift_y: f32) -> Result<Array2<u8>> {
let imgx = 800;
let imgy = 600;
let scalex = 3.0 / imgx as f32;
let scaley = 3.0 / imgy as f32;
let mut im = Array2::<u8>::zeros((imgy, imgx));
for x in 0..imgx {
for y in 0..imgy {
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(cy, cx);
let mut i = 0;
while i < 255 && z.norm() <= 2.0 {
z = z * z + c;
i += 1;
}
im[[y, x]] = i as u8;
}
}
Ok(im)
}
#[cfg(test)]
mod tests {
use super::*;
use anyhow::Result;
use ndarray::Array2;
use tempfile::NamedTempFile;
#[test]
fn test_serialization() -> Result<()> {
let file = NamedTempFile::new()?;
let t = AffineTransform::new([1.2, 0.3, -0.4, 0.9, 10.2, -9.5], [59.5, 49.5], [120, 100]);
t.to_file(file.path().to_path_buf())?;
let s = AffineTransform::from_file(file.path().to_path_buf())?;
assert_eq!(s, t);
Ok(())
}
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 = AffineTransform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
let n: Array2<_> = transform.transform_image_bspline(j.view())?.into();
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_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 j0 = j.clone();
let k0 = k.clone();
let transform = AffineTransform::new([1., 0., 0., 1., 120., -10.], origin, [shape[0], shape[1]]);
// make sure j & k weren't mutated
assert!(j.iter().zip(j0.iter()).map(|(a, b)| a == b).all(|x| x));
assert!(k.iter().zip(k0.iter()).map(|(a, b)| a == b).all(|x| x));
let n: Array2<_> = transform.transform_image_nearest_neighbor(j.view())?.into();
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 j0 = j.clone();
let k0 = k.clone();
let t = AffineTransform::register_translation(j.view(), k.view())?;
// make sure j & k weren't mutated
assert!(j.iter().zip(j0.iter()).map(|(a, b)| a == b).all(|x| x));
assert!(k.iter().zip(k0.iter()).map(|(a, b)| a == b).all(|x| x));
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, "d: {}, t: {:?}", d, t.parameters);
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 = AffineTransform::new([1.2, 0., 0., 1., 5., 7.], origin, [shape[0], shape[1]]);
let k: Array2<_> = s.transform_image_bspline(j.view())?.into();
let t = AffineTransform::register_affine(j.view(), k.view())?.inverse()?;
let d = (t.matrix() - s.matrix()).powi(2).sum();
assert!(d < 0.025, "d: {}, t: {:?}", d, t.parameters);
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,
}
}

362
src/sys.rs
View File

@@ -1,362 +0,0 @@
use crate::PixelType;
use anyhow::Result;
use libc::{c_double, c_uint};
use ndarray::{Array2, AsArray, Ix2};
use one_at_a_time_please::one_at_a_time;
use std::ptr;
macro_rules! register_fn {
($($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 {
($($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! {
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(crate) fn interp<'a, A, T>(
parameters: [f64; 6],
origin: [f64; 2],
image: A,
bspline_or_nn: bool,
) -> Result<Array2<T>>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
let image = image.into();
let shape: Vec<usize> = image.shape().to_vec();
let width = shape[1] as c_uint;
let height = shape[0] as c_uint;
let mut im: Vec<_> = image.into_iter().cloned().collect();
let im_ptr: *mut T = ptr::from_mut(unsafe { &mut *im.as_mut_ptr() });
match T::PT {
1 => unsafe {
interp_u8(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut u8),
bspline_or_nn,
);
},
2 => unsafe {
interp_i8(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut i8),
bspline_or_nn,
);
},
3 => unsafe {
interp_u16(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut u16),
bspline_or_nn,
);
},
4 => unsafe {
interp_i16(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut i16),
bspline_or_nn,
);
},
5 => unsafe {
interp_u32(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut u32),
bspline_or_nn,
);
},
6 => unsafe {
interp_i32(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut i32),
bspline_or_nn,
);
},
7 => unsafe {
interp_u64(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut u64),
bspline_or_nn,
);
},
8 => unsafe {
interp_i64(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut i64),
bspline_or_nn,
);
},
9 => unsafe {
interp_f32(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut f32),
bspline_or_nn,
);
},
10 => unsafe {
interp_f64(
width,
height,
parameters.as_ptr(),
origin.as_ptr(),
&mut (im_ptr as *mut f64),
bspline_or_nn,
);
},
_ => {}
}
Ok(Array2::from_shape_vec(
(shape[0], shape[1]),
im.into_iter().collect(),
)?)
}
#[one_at_a_time]
pub(crate) fn register<'a, A, T>(
fixed: A,
moving: A,
translation_or_affine: bool,
) -> Result<([f64; 6], [f64; 2], [usize; 2])>
where
T: 'a + PixelType,
A: AsArray<'a, T, Ix2>,
{
let fixed = fixed.into();
let moving = moving.into();
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_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_ptr as *const u8,
moving_ptr as *const u8,
translation_or_affine,
&mut transform_ptr,
)
};
}
2 => {
unsafe {
register_i8(
width,
height,
fixed_ptr as *const i8,
moving_ptr as *const i8,
translation_or_affine,
&mut transform_ptr,
)
};
}
3 => {
unsafe {
register_u16(
width,
height,
fixed_ptr as *const u16,
moving_ptr as *const u16,
translation_or_affine,
&mut transform_ptr,
)
};
}
4 => {
unsafe {
register_i16(
width,
height,
fixed_ptr as *const i16,
moving_ptr as *const i16,
translation_or_affine,
&mut transform_ptr,
)
};
}
5 => {
unsafe {
register_u32(
width,
height,
fixed_ptr as *const u32,
moving_ptr as *const u32,
translation_or_affine,
&mut transform_ptr,
)
};
}
6 => {
unsafe {
register_i32(
width,
height,
fixed_ptr as *const i32,
moving_ptr as *const i32,
translation_or_affine,
&mut transform_ptr,
)
};
}
7 => {
unsafe {
register_u64(
width,
height,
fixed_ptr as *const u64,
moving_ptr as *const u64,
translation_or_affine,
&mut transform_ptr,
)
};
}
8 => {
unsafe {
register_i64(
width,
height,
fixed_ptr as *const i64,
moving_ptr as *const i64,
translation_or_affine,
&mut transform_ptr,
)
};
}
9 => {
unsafe {
register_f32(
width,
height,
fixed_ptr as *const f32,
moving_ptr as *const f32,
translation_or_affine,
&mut transform_ptr,
)
};
}
10 => {
unsafe {
register_f64(
width,
height,
fixed_ptr as *const f64,
moving_ptr as *const f64,
translation_or_affine,
&mut transform_ptr,
)
};
}
_ => {}
}
// 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,
transform[1] as f64,
transform[2] as f64,
transform[3] as f64,
transform[4] as f64,
transform[5] as f64,
],
[
((shape[0] - 1) as f64) / 2f64,
((shape[1] - 1) as f64) / 2f64,
],
[shape[0], shape[1]],
))
}