commit e2184b2d331cc85b6fcf038893133c2e201ddb67 Author: Wim Pomp Date: Sat Mar 1 19:40:37 2025 +0100 - First commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a1a09d9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,84 @@ +/target +/Cargo.lock + +# Byte-compiled / optimized / DLL files +__pycache__/ +.pytest_cache/ +*.py[cod] + +# C extensions +*.so +*.dylib +*.dll + +# Distribution / packaging +.Python +.venv/ +env/ +bin/ +build/ +develop-eggs/ +dist/ +eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +include/ +man/ +venv/ +*.egg-info/ +.installed.cfg +*.egg + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt +pip-selfcheck.json + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.cache +nosetests.xml +coverage.xml + +# Translations +*.mo + +# Mr Developer +.mr.developer.cfg +.project +.pydevproject + +# Rope +.ropeproject + +# Django stuff: +*.log +*.pot + +.DS_Store + +# Sphinx documentation +docs/_build/ + +# PyCharm +.idea/ + +# VSCode +.vscode/ + +# Pyenv +.python-version + +/tests/files/* +/cpp/CMakeFiles/ +/cpp/cmake_install.cmake +/cpp/CMakeCache.txt +/cpp/Makefile +/sitk +*.nii +TransformParameters* diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..90927b5 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,31 @@ +[package] +name = "sitk-registration-sys" +version = "2025.3.0" +edition = "2024" +license = "MIT OR Apache-2.0" +description = "register and interpolate images" +rust-version = "1.85.0" +authors = ["Wim Pomp "] +homepage = "https://github.com/wimpomp/sitk-registration-sys" +repository = "https://github.com/wimpomp/sitk-registration-sys" +documentation = "https://docs.rs/sitk-registration-sys" +readme = "README.md" +keywords = ["registration", "affine", "bspline", "transform"] +categories = ["multimedia::images", "science"] + +[lib] +name = "sitk_regsitration_sys" +crate-type = ["cdylib", "rlib"] + +[dependencies] +anyhow = "1.0.96" +libc = "0.2.170" +ndarray = "0.16.1" +num = "0.4.3" + +[dev-dependencies] +tiffwrite = "2025.2.0" + +[build-dependencies] +cmake = "0.1.54" +git2 = "0.20.0" \ No newline at end of file diff --git a/LICENSE-APACHE b/LICENSE-APACHE new file mode 100644 index 0000000..1b5ec8b --- /dev/null +++ b/LICENSE-APACHE @@ -0,0 +1,176 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + +TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + +1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + +2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + +3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + +4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + +5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + +6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + +7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + +8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + +9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + +END OF TERMS AND CONDITIONS diff --git a/LICENSE-MIT b/LICENSE-MIT new file mode 100644 index 0000000..31aa793 --- /dev/null +++ b/LICENSE-MIT @@ -0,0 +1,23 @@ +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. diff --git a/README.md b/README.md new file mode 100644 index 0000000..e8969d4 --- /dev/null +++ b/README.md @@ -0,0 +1,31 @@ +# sitk-registration-sys + +This crate does two things: +- 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 + +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 +in a shared library. Because of this, compilation of this crate requires quite some time, as +wel as cmake. + +## Examples +### Registration +``` + let image_a = (some Array2); + let iameg_b = (some transformed Array2); + let transform = Transform::register_affine(image_a.view(), image_b.view())?; + println!("transform: {:#?}", transform); +``` + +### Interpolation +``` + let image = (Some Array2); + let shape = image.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 transformed_image = transform.transform_image_bspline(image.view())?; +``` \ No newline at end of file diff --git a/build.rs b/build.rs new file mode 100644 index 0000000..7ce670d --- /dev/null +++ b/build.rs @@ -0,0 +1,62 @@ +use cmake::Config; +use git2::Repository; +use std::ffi::OsStr; +use std::path::PathBuf; + +fn main() { + if std::env::var("DOCS_RS").is_err() { + let out_dir = PathBuf::from(std::env::var("OUT_DIR").expect("OUT_DIR is undefined")); + 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() { + d.join("sitk").to_path_buf() + } 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) + .define("BUILD_TESTING", "OFF") + .define("WRAP_CSHARP", "OFF") + .define("WRAP_JAVA", "OFF") + .define("WRAP_LUA", "OFF") + .define("WRAP_R", "OFF") + .define("WRAP_RUBY", "OFF") + .define("WRAP_TCL", "OFF") + .define("WRAP_PYTHON", "OFF") + .define("WRAP_DEFAULT", "OFF") + .define("SimpleITK_USE_ELASTIX", "ON") + .build(); + } + // println!("cargo::rustc-env=CMAKE_INSTALL_PREFIX=/home/wim/code/rust/sitk-sys/cpp"); + println!( + "cargo::rustc-env=CMAKE_INSTALL_PREFIX={}", + out_dir.display() + ); + let path = Config::new("cpp") + .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(); + println!("cargo::rustc-link-arg=-Wl,-rpath,{}", path.display()); + println!("cargo::rustc-link-search={}", path.join("build").display()); + println!("cargo::rustc-link-lib=dylib=sitk_adapter"); + println!("cargo::rerun-if-changed=build.rs"); + println!("cargo::rerun-if-changed=cpp/*.cxx"); + } +} diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt new file mode 100644 index 0000000..bf9a83b --- /dev/null +++ b/cpp/CMakeLists.txt @@ -0,0 +1,11 @@ +cmake_minimum_required(VERSION 3.16.3) + +project(sitk_adapter) + +set(ENV{Elastix_DIR} "../sitk/build/Elastix-build" ) +set(ENV{ITK_DIR} "~/code/c/SimpleITK/build/ITK-build" ) +set(ENV{SimpleITK_DIR} "~/code/c/SimpleITK/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 .) diff --git a/cpp/sitk_adapter.cxx b/cpp/sitk_adapter.cxx new file mode 100644 index 0000000..50dccf6 --- /dev/null +++ b/cpp/sitk_adapter.cxx @@ -0,0 +1,415 @@ +#include +#include +#include + +namespace sitk = itk::simple; + +using namespace std; + + +template +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 matrix = {transform[0], transform[1], transform[2], transform[3]}; + vector translation = {transform[4], transform[5]}; + vector 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); +} + +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); +} + +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); +} + +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); +} + +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); +} + +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); +} + +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); +} + +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); +} + + +void +reg( + sitk::Image fixed, + sitk::Image moving, + bool t_or_a, + double** transform +) { + try { + string kind = (t_or_a == false) ? "translation" : "affine"; + sitk::ElastixImageFilter tfilter = sitk::ElastixImageFilter(); + tfilter.LogToConsoleOff(); + tfilter.SetFixedImage(fixed); + tfilter.SetMovingImage(moving); + tfilter.SetParameterMap(sitk::GetDefaultParameterMap(kind)); + 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 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]); + } + } + } + } + } 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); +} \ No newline at end of file diff --git a/interp_test.tif b/interp_test.tif new file mode 100644 index 0000000..089d6b3 Binary files /dev/null and b/interp_test.tif differ diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..fc5bacc --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,282 @@ +mod sys; + +use crate::sys::{PixelType, interp, register}; +use anyhow::{Result, anyhow}; +use ndarray::{Array2, ArrayView2, array, s}; +use std::ops::Mul; +use std::path::PathBuf; + +#[derive(Clone, Debug)] +pub struct Transform { + pub parameters: [f64; 6], + pub dparameters: [f64; 6], + pub origin: [f64; 2], + pub shape: [usize; 2], +} + +impl Mul for Transform { + type Output = Transform; + + #[allow(clippy::suspicious_arithmetic_impl)] + fn mul(self, other: Transform) -> Transform { + let m = self.matrix().dot(&other.matrix()); + 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 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( + fixed: ArrayView2, + moving: ArrayView2, + ) -> Result { + 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( + fixed: ArrayView2, + moving: ArrayView2, + ) -> Result { + 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(file: PathBuf) -> Result { + todo!() + } + + /// write a transform to a file + pub fn to_file(&self, file: PathBuf) -> Result<()> { + todo!() + } + + /// 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(&self, image: ArrayView2) -> Result> { + interp(self.parameters, self.origin, image, false) + } + + /// transform an image using bspline interpolation + pub fn transform_image_nearest_neighbor( + &self, + image: ArrayView2, + ) -> Result> { + interp(self.parameters, self.origin, image, true) + } + + /// get coordinates resulting from transforming input coordinates + pub fn transform_coordinates(&self, coordinates: ArrayView2) -> Result> + where + T: Clone + 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 { + 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 { + 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 { + fn det(a: ArrayView2) -> f64 { + (a[[0, 0]] * a[[1, 1]]) - (a[[0, 1]] * a[[1, 0]]) + } + + let m = self.matrix(); + let d0 = det(m.slice(s![1.., 1..])); + if d0 == 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, + ]; + + 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 tiffwrite::IJTiffFile; + + /// An example of generating julia fractals. + fn julia_image() -> Result> { + let imgx = 800; + let imgy = 600; + + let scalex = 3.0 / imgx as f32; + let scaley = 3.0 / imgy as f32; + + let mut im = Array2::::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 c = Complex::new(-0.4, 0.6); + let mut z = Complex::new(cx, cy); + + 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_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)?; + + let t = Transform::register_affine(k.view(), j.view())?; + println!("t: {:#?}", t); + println!("m: {:#?}", t.matrix()); + println!("i: {:#?}", t.inverse()?.matrix()); + Ok(()) + } +} diff --git a/src/sys.rs b/src/sys.rs new file mode 100644 index 0000000..42a2059 --- /dev/null +++ b/src/sys.rs @@ -0,0 +1,359 @@ +use anyhow::Result; +use libc::{c_double, c_uint}; +use ndarray::{Array2, ArrayView2}; +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, + ); + }; +} + +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, + ); + }; +} + +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); +} + +pub trait PixelType: Clone { + const PT: u8; +} + +macro_rules! sitk_impl { + ($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); +#[cfg(target_pointer_width = "64")] +sitk_impl!(usize, 7); +#[cfg(target_pointer_width = "32")] +sitk_impl!(usize, 5); +#[cfg(target_pointer_width = "64")] +sitk_impl!(isize, 8); +#[cfg(target_pointer_width = "32")] +sitk_impl!(isize, 6); +sitk_impl!(f32, 9); +sitk_impl!(f64, 10); + +pub(crate) fn interp( + parameters: [f64; 6], + origin: [f64; 2], + image: ArrayView2, + bspline_or_nn: bool, +) -> Result> { + let shape: Vec = 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(), + )?) +} + +pub(crate) fn register( + fixed: ArrayView2, + moving: ArrayView2, + translation_or_affine: bool, +) -> Result<([f64; 6], [f64; 2], [usize; 2])> { + let shape: Vec = 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 mut transform: Vec = vec![0.0; 6]; + let mut transform_ptr: *mut c_double = ptr::from_mut(unsafe { &mut *transform.as_mut_ptr() }); + + match T::PT { + 1 => { + unsafe { + register_u8( + width, + height, + fixed.as_ptr() as *const u8, + moving.as_ptr() as *const u8, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 2 => { + unsafe { + register_i8( + width, + height, + fixed.as_ptr() as *const i8, + moving.as_ptr() as *const i8, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 3 => { + unsafe { + register_u16( + width, + height, + fixed.as_ptr() as *const u16, + moving.as_ptr() as *const u16, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 4 => { + unsafe { + register_i16( + width, + height, + fixed.as_ptr() as *const i16, + moving.as_ptr() as *const i16, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 5 => { + unsafe { + register_u32( + width, + height, + fixed.as_ptr() as *const u32, + moving.as_ptr() as *const u32, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 6 => { + unsafe { + register_i32( + width, + height, + fixed.as_ptr() as *const i32, + moving.as_ptr() as *const i32, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 7 => { + unsafe { + register_u64( + width, + height, + fixed.as_ptr() as *const u64, + moving.as_ptr() as *const u64, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 8 => { + unsafe { + register_i64( + width, + height, + fixed.as_ptr() as *const i64, + moving.as_ptr() as *const i64, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 9 => { + unsafe { + register_f32( + width, + height, + fixed.as_ptr() as *const f32, + moving.as_ptr() as *const f32, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + 10 => { + unsafe { + register_f64( + width, + height, + fixed.as_ptr() as *const f64, + moving.as_ptr() as *const f64, + translation_or_affine, + &mut transform_ptr, + ) + }; + } + _ => {} + } + + 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]], + )) +}