diff --git a/.appveyor.yml b/.appveyor.yml index 3743e67..a1503a6 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -17,18 +17,27 @@ build: matrix: fast_finish: true -before_build: - - cmd: if "%PLATFORM%"=="x86" call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvars32.bat" - - cmd: if "%PLATFORM%"=="x64" call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvars64.bat" - - mkdir build - - cd build +environment: + matrix: + - PYTHON: "C:\\Python36-x64" + +install: +- git submodule update --init --recursive +- set PATH=%PYTHON%;%PYTHON%\Scripts;%PATH% +- pip install -r requirements.txt - - cmd: if "%PLATFORM%"=="x86" set GENERATOR="Visual Studio 15 2017" - - cmd: if "%PLATFORM%"=="x64" set GENERATOR="Visual Studio 15 2017 Win64" - - cmd: cmake .. -DSTK_BUILD_TESTS=1 -DSTK_WARNINGS_ARE_ERRORS=1 -DCMAKE_BUILD_TYPE=%CONFIGURATION% -G%GENERATOR% +before_build: +- mkdir build +- cd build +- cmake .. -DCMAKE_BUILD_TYPE=%CONFIGURATION% -DCMAKE_GENERATOR_PLATFORM=%PLATFORM% -DSTK_BUILD_TESTS=1 -DSTK_WARNINGS_ARE_ERRORS=1 build_script: - - cmd: cmake --build . --config %CONFIGURATION% +- cmake --build . --config %CONFIGURATION% +- cd .. +- python setup.py install -DCMAKE_GENERATOR_PLATFORM=%PLATFORM% test_script: - - ctest --output-on-failure \ No newline at end of file +- cd build +- ctest --output-on-failure +- cd .. +- python -m unittest discover test diff --git a/.gitignore b/.gitignore index 8234407..5276635 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,9 @@ build_debug/ build_vs/ build/ +dist/ +*.egg-info +__pycache__ +.env/ .vscode/ -.ycm_extra_conf.py +.ycm_extra_conf.py \ No newline at end of file diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..9bb357f --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "third_party/pybind11"] + path = third_party/pybind11 + url = https://github.com/pybind/pybind11.git diff --git a/.travis.yml b/.travis.yml index 6b82c07..dbf23e1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,37 +1,56 @@ -language: cpp +language: python +python: +- 3.6 matrix: - fast_finish: true include: - - os: linux + - name: trusty + os: linux dist: trusty - - os: linux + env: + - PYTHON_EXECUTABLE=~/virtualenv/python3.6/bin/python + - name: xenial + os: linux dist: xenial - - os: linux + env: + - PYTHON_EXECUTABLE=~/virtualenv/python3.6/bin/python + - name: g++8 + os: linux addons: apt: sources: - - ubuntu-toolchain-r-test + - ubuntu-toolchain-r-test packages: - - g++-6 - env: - - MATRIX_EVAL="CC=gcc-6 && CXX=g++-6" - - os: linux - addons: - apt: - sources: - - ubuntu-toolchain-r-test - packages: - - g++-7 - env: - - MATRIX_EVAL="CC=gcc-7 && CXX=g++-7" + - g++-8 + before_install: + - sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-8 90 + - sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-8 90 + env: + - PYTHON_EXECUTABLE=~/virtualenv/python3.6/bin/python - name: osx os: osx language: generic + before_install: + - brew install llvm libomp + - python3 -m pip install virtualenv + - virtualenv venv -p python3 + - source venv/bin/activate + env: + - PYTHON_EXECUTABLE=../venv/bin/python + - HOMEBREW_NO_AUTO_UPDATE=1 + - CC=/usr/local/opt/llvm/bin/clang + - CXX=/usr/local/opt/llvm/bin/clang++ + +install: +- pip install -r requirements.txt script: - mkdir build - cd build -- cmake .. -DCMAKE_BUILD_TYPE=Release -DSTK_BUILD_TESTS=1 -DSTK_WARNINGS_ARE_ERRORS=1 +- cmake .. -DCMAKE_BUILD_TYPE=Release -DSTK_BUILD_TESTS=1 + -DSTK_BUILD_PYTHON_WRAPPER=ON -DPYTHON_EXECUTABLE=$PYTHON_EXECUTABLE - cmake --build . --config Release -- ctest --output-on-failure \ No newline at end of file +- ctest --output-on-failure +- cd .. +- python setup.py install +- python -m unittest discover test diff --git a/CMakeLists.txt b/CMakeLists.txt index 2626c45..9987a62 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.8) option(STK_BUILD_EXAMPLES "Build examples" ON) option(STK_BUILD_TESTS "Build unit tests" ON) +option(STK_BUILD_PYTHON_WRAPPER "Build Python wrapper" OFF) option(STK_USE_CUDA "Enables CUDA support" OFF) option(STK_WARNINGS_ARE_ERRORS "Warnings are treated as errors" OFF) option(STK_BUILD_WITH_DEBUG_INFO "Includes debug info in release builds" OFF) @@ -56,7 +57,9 @@ if (MSVC) set(EXTRA_FLAGS "${extra_flags} /fp:fast") endif() elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU") - set(EXTRA_FLAGS "-Wall -Wextra -pedantic -Wno-missing-field-initializers") + set(EXTRA_FLAGS "-Wall -Wextra -pedantic -Wno-missing-field-initializers -Wno-attributes") + + # -Wno-attributes thanks to pybind11 if (STK_WARNINGS_ARE_ERRORS) set(EXTRA_FLAGS "${extra_flags} -Werror") @@ -88,6 +91,7 @@ endif() if (STK_USE_CUDA) add_definitions(-DSTK_USE_CUDA) enable_language(CUDA) + set(CMAKE_CUDA_STANDARD 11) if (STK_ENABLE_FAST_MATH) set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --use_fast_math") @@ -109,10 +113,11 @@ set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} ${EXTRA_FLAGS_RELEASE}") set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${EXTRA_LINK_FLAGS}") set(CMAKE_EXE_LINKER_FLAGS_RELEASE "${CMAKE_EXE_LINKER_FLAGS_RELEASE} ${EXTRA_LINK_FLAGS_RELEASE}") +include_directories(src) + # TODO: Determine QNANHIBIT for NrrdIO set(QNANHIBIT 1) - if(NOT WIN32) find_package(ZLIB QUIET) find_package(NIFTI QUIET) diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..6e9c09d --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,4 @@ +recursive-include src * +recursive-include third_party * +include cmake/* +include CMakeLists.txt diff --git a/README.md b/README.md index 835e66c..f61b17d 100644 --- a/README.md +++ b/README.md @@ -23,3 +23,37 @@ CMake was not compatible with the one required by CUDA, it is possible to specify a different executable with `-DCMAKE_CUDA_FLAGS="-ccbin gcc-XX"`, where `gcc-XX` is a version of `gcc` compatible with your CUDA version. +# Python API + +A minimalistic Python API is also provided. + +## Install + +```bash +python setup.py install +``` + +## Example usage + +```python +import stk +import numpy as np + +# Create volume directly from numpy +vol = stk.Volume(np.zeros((5,5,5)).astype(np.float32), spacing=(2,2,2)) +# or read volume from file +vol = stk.read_volume('test.nrrd') + +# Modify data (numpy array points to the volume data) +data = np.array(vol, copy=False) +data[0:10] = 0.0 + +# Access meta data +vol.origin = (2, 2, 2) +vol.spacing = (3, 3, 3) + +# Write volume +stk.write_volume('test-out.nrrd', vol) +``` + + diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..296d654 --- /dev/null +++ b/requirements.txt @@ -0,0 +1 @@ +numpy \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..31dd519 --- /dev/null +++ b/setup.py @@ -0,0 +1,91 @@ +import os +import sys +import platform +import subprocess + +from pprint import pprint + +from setuptools import setup, Extension +from setuptools.command.build_ext import build_ext + +with open("README.md", encoding="utf-8") as f: + readme = f.read() + +# Parse command line flags +flags = {k: 'OFF' for k in ['--debug', '--use-cuda', '--use-itk']} +for flag in flags.keys(): + if flag in sys.argv: + flags[flag] = 'ON' + sys.argv.remove(flag) + +# Command line flags forwarded to CMake +cmake_cmd_args = [] +for f in sys.argv: + if f.startswith('-D'): + cmake_cmd_args.append(f) + +for f in cmake_cmd_args: + sys.argv.remove(f) + + +class CMakeExtension(Extension): + def __init__(self, name, cmake_lists_dir=''): + Extension.__init__(self, name, sources=[]) + self.cmake_lists_dir = os.path.abspath(cmake_lists_dir) + + +class CMakeBuild(build_ext): + def run(self): + try: + subprocess.check_output(['cmake', '--version']) + except OSError: + raise RuntimeError('Cannot find CMake executable') + + for ext in self.extensions: + self.build_extension(ext) + + def build_extension(self, ext): + + extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) + cfg = 'Debug' if flags['--debug'] == 'ON' else 'Release' + build_args = ['--config', cfg] + + cmake_args = [ + '-DSTK_BUILD_PYTHON_WRAPPER=ON', + '-DSTK_BUILD_TESTS=OFF', + '-DSTK_BUILD_WITH_DEBUG_INFO=%s' % ('ON' if cfg == 'Debug' else 'OFF'), + '-DCMAKE_BUILD_TYPE=%s' % cfg, + '-DSTK_USE_CUDA=%s' % flags['--use-cuda'], + '-DSTK_ITK_BRIDGE=%s' % flags['--use-itk'], + '-DPYTHON_EXECUTABLE=' + sys.executable, + '-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, + ] + + if platform.system() == "Windows": + cmake_args += ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}'.format(cfg.upper(), extdir)] + + cmake_args += cmake_cmd_args + pprint(cmake_args) + + if not os.path.exists(self.build_temp): + os.makedirs(self.build_temp) + + subprocess.check_call(['cmake', ext.cmake_lists_dir] + cmake_args, cwd=self.build_temp) + subprocess.check_call(['cmake', '--build', '.'] + build_args, cwd=self.build_temp) + + +setup( + name='python-stk', + version='0.4', + author='Simon Ekström', + author_email='', + description='', + long_description=readme, + long_description_content_type='text/markdown', + install_requires=['numpy'], + packages=['stk'], + ext_modules=[CMakeExtension('_stk', '.')], + cmdclass={'build_ext': CMakeBuild}, + zip_safe=False, +) + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 65fdd29..d7e040f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,105 +1,5 @@ -set(STK_SRCS - "stk/common/error.cpp" - "stk/common/error.h" - "stk/common/log.cpp" - "stk/common/log.h" - "stk/common/stack_trace.cpp" - "stk/common/stack_trace.h" - "stk/cuda/volume.h" - "stk/filters/decomposable_filter.h" - "stk/filters/decomposable_filter.inl" - "stk/filters/gaussian_filter.cpp" - "stk/filters/gaussian_filter.h" - "stk/filters/map.h" - "stk/filters/sobel.cpp" - "stk/filters/sobel.h" - "stk/filters/vector_calculus.h" - "stk/filters/vector_calculus.inl" - "stk/io/io.cpp" - "stk/io/io.h" - "stk/io/nifti.cpp" - "stk/io/nifti.h" - "stk/io/nrrd.cpp" - "stk/io/nrrd.h" - "stk/io/vtk.cpp" - "stk/io/vtk.h" - "stk/math/float2.h" - "stk/math/float3.h" - "stk/math/float4.h" - "stk/math/math.cpp" - "stk/math/math.h" - "stk/math/types.h" - "stk/image/dim3.h" - "stk/image/types.cpp" - "stk/image/types.h" - "stk/image/volume.cpp" - "stk/image/volume.h" - "stk/image/volume.inl" -) +add_subdirectory(stk) -if (STK_USE_CUDA) - set(STK_SRCS - ${STK_SRCS} - "stk/cuda/cuda.cpp" - "stk/cuda/cuda.h" - "stk/cuda/stream.cpp" - "stk/cuda/stream.h" - "stk/filters/gpu/normalize.cu" - "stk/image/gpu_volume.cpp" - "stk/image/gpu_volume.cu" - "stk/image/gpu_volume.h" - ) +if(STK_BUILD_PYTHON_WRAPPER) + add_subdirectory(python_wrapper) endif() - -if(STK_ITK_BRIDGE) - set(STK_SRCS - ${STK_SRCS} - "stk/itk_bridge/itk_convert.h" - "stk/itk_bridge/itk_io.cpp" - "stk/itk_bridge/itk_io.h" - ) -endif() - -add_library(stk STATIC ${STK_SRCS}) - -target_include_directories(stk PUBLIC - ${CMAKE_CURRENT_SOURCE_DIR} - ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES} - ) - -if(NIFTI_FOUND) - target_include_directories(stk PUBLIC ${NIFTI_INCLUDE_DIR}) - target_link_libraries(stk PUBLIC ${NIFTI_LIBRARIES}) -else() - target_link_libraries(stk PUBLIC niftiio) -endif() - -if(NrrdIO_FOUND) - target_include_directories(stk PUBLIC ${NrrdIO_INCLUDE_DIR}) - target_link_libraries(stk PUBLIC ${NrrdIO_LIBRARIES}) -else() - target_link_libraries(stk PUBLIC NrrdIO) -endif() - -if(STK_ITK_BRIDGE) - add_definitions(-DSTK_ITK_BRIDGE) - find_package(ITK REQUIRED) - include(${ITK_USE_FILE}) - - target_include_directories(stk PRIVATE ${ITK_INCLUDE_DIRS}) - target_link_libraries(stk PRIVATE ${ITK_LIBRARIES}) -endif() - -if(STK_STACK_TRACE) - add_definitions(-DSTK_STACK_TRACE) - - if(UNIX) - find_package(Backward) - if(Backward_FOUND) - add_backward(stk) - target_link_libraries(stk PUBLIC ${BACKWARD_LIBRARIES}) - add_definitions(-DSTK_USE_BACKWARD) - endif() - endif() -endif() - diff --git a/src/python_wrapper/CMakeLists.txt b/src/python_wrapper/CMakeLists.txt new file mode 100644 index 0000000..39a33a7 --- /dev/null +++ b/src/python_wrapper/CMakeLists.txt @@ -0,0 +1,8 @@ +pybind11_add_module(_stk _stk.cpp) + +target_include_directories(_stk PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} + ${PYTHON_INCLUDE_DIRS} ${PYBIND11_INCLUDE_DIR} ${pybind11_INCLUDE_DIR}) +target_link_libraries(_stk PUBLIC stk) + +set(STK_PYTHON_PATH "$ENV{PYTHONPATH}${SEP}${CMAKE_CURRENT_BINARY_DIR}" CACHE STRING "Python path for Sphinx") + diff --git a/src/python_wrapper/_stk.cpp b/src/python_wrapper/_stk.cpp new file mode 100644 index 0000000..529b18f --- /dev/null +++ b/src/python_wrapper/_stk.cpp @@ -0,0 +1,491 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace py = pybind11; + + +namespace pybind11 { namespace detail { + template <> struct type_caster { + public: + PYBIND11_TYPE_CASTER(float3, _("float3")); + + bool load(handle src, bool) { + try { + tuple t = py::cast(src); + if (t.size() != 3) { + return false; + } + value = float3{ + t[0].cast(), + t[1].cast(), + t[2].cast() + }; + } catch (std::exception& e) { + return false; + } + return true; + } + static handle cast(float3 src, return_value_policy /* policy */, handle /* parent */) { + return make_tuple( + src.x, + src.y, + src.z + ).release(); + } + }; + template <> struct type_caster { + public: + PYBIND11_TYPE_CASTER(dim3, _("dim3")); + + bool load(handle src, bool) { + try { + tuple t = py::cast(src); + if (t.size() != 3) { + return false; + } + value = dim3{ + t[0].cast(), + t[1].cast(), + t[2].cast() + }; + } catch (std::exception& e) { + return false; + } + return true; + } + static handle cast(dim3 src, return_value_policy /* policy */, handle /* parent */) { + return make_tuple( + src.x, + src.y, + src.z + ).release(); + } + }; + template <> struct type_caster { + public: + PYBIND11_TYPE_CASTER(int3, _("int3")); + + bool load(handle src, bool) { + try { + tuple t = py::cast(src); + if (t.size() != 3) { + return false; + } + value = int3{ + t[0].cast(), + t[1].cast(), + t[2].cast() + }; + } catch (std::exception& e) { + return false; + } + return true; + } + static handle cast(int3 src, return_value_policy /* policy */, handle /* parent */) { + return make_tuple( + src.x, + src.y, + src.z + ).release(); + } + }; + template <> struct type_caster { + public: + PYBIND11_TYPE_CASTER(Matrix3x3f, _("Matrix3x3f")); + + bool load(handle src, bool) { + try { + array_t a = py::cast>(src); + if (a.ndim() != 2 + || a.shape(0) != 3 + || a.shape(1) != 3) { + return false; + } + value = Matrix3x3f{ + float3{a.at(0, 0), a.at(0, 1), a.at(0, 2)}, + float3{a.at(1, 0), a.at(1, 1), a.at(1, 2)}, + float3{a.at(2, 0), a.at(2, 1), a.at(2, 2)} + }; + } catch (std::exception& e) { + return false; + } + return true; + } + static handle cast(Matrix3x3f src, return_value_policy /* policy */, handle /* parent */) { + return py::array_t( + {3, 3}, + &src._rows[0].x + ).release(); + } + }; + +}} // namespace pybind11::detail + +std::vector get_shape(const stk::Volume& vol) +{ + dim3 size = vol.size(); + size_t ncomp = stk::num_components(vol.voxel_type()); + if (ncomp == 1) { + // Flip directions since size is in (width, height, depth), + // numpy expects (depth, height, width) + return { + size.z, + size.y, + size.x + }; + } + else { + return { + size.z, + size.y, + size.x, + ncomp + }; + } +} + +std::vector get_strides(const stk::Volume& vol) +{ + const size_t* strides = vol.strides(); + size_t ncomp = stk::num_components(vol.voxel_type()); + + // Reverse, because of numpy ordering + if (ncomp == 1) { + return { + strides[2], + strides[1], + strides[0] + }; + } + else { + return { + strides[2], + strides[1], + strides[0], + strides[0]/ncomp + }; + } +} + +const char* format_descriptor(stk::Type base_type) +{ + if (base_type == stk::Type_Char) return "b"; + else if (base_type == stk::Type_UChar) return "B"; + else if (base_type == stk::Type_Short) return "h"; + else if (base_type == stk::Type_UShort) return "H"; + else if (base_type == stk::Type_Int) return "i"; + else if (base_type == stk::Type_UInt) return "I"; + else if (base_type == stk::Type_Float) return "f"; + else if (base_type == stk::Type_Double) return "d"; + return ""; +} + +py::buffer_info get_buffer_info(stk::Volume& v) +{ + stk::Type base_type = stk::base_type(v.voxel_type()); + int ncomp = stk::num_components(v.voxel_type()); + + return py::buffer_info( + v.ptr(), /* Pointer to buffer */ + stk::type_size(base_type), /* Size of one scalar */ + format_descriptor(base_type), /* Python struct-style format descriptor */ + ncomp == 1 ? 3 : 4, /* Number of dimensions */ + + get_shape(v), /* Buffer dimensions */ + get_strides(v) /* Strides (in bytes) for each index */ + ); +} + +stk::Volume numpy_to_volume(const py::array arr) +{ + // Ensure ordering + py::array buf = py::array::ensure(arr, py::array::c_style | py::array::forcecast); + + if (!(buf.ndim() == 3 || buf.ndim() == 4)) { + throw std::invalid_argument("Expected a 3D or a 4D array"); + } + + dim3 size { + std::uint32_t(buf.shape(2)), + std::uint32_t(buf.shape(1)), + std::uint32_t(buf.shape(0)), + }; + + stk::Type base_type = stk::Type_Unknown; + if (py::isinstance>(arr)) { + base_type = stk::Type_Char; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_Char; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_UChar; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_Short; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_UShort; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_Int; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_UInt; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_Float; + } + else if (py::isinstance>(arr)) { + base_type = stk::Type_Double; + } + else { + throw std::invalid_argument("Unsupported type"); + } + + int ncomp = 1; + if (buf.ndim() == 4) { + ncomp = buf.shape(3); + } + + if (ncomp < 1 || ncomp > 4) { + throw std::invalid_argument("Unsupported number of channels"); + } + + // NOTE: the value of ndim can be ambiguous, e.g. + // ndim == 3 may be a scalar volume or a vector 2D image... + + return stk::Volume( + size, + stk::build_type(base_type, ncomp), + buf.data() + ); +} + +stk::Volume make_volume( + py::array arr, + const float3& origin, + const float3& spacing, + const Matrix3x3f& direction +) { + stk::Volume vol = numpy_to_volume(arr); + vol.set_origin(origin); + vol.set_spacing(spacing); + vol.set_direction(direction); + return vol; +} + + +std::string divergence_docstring = +R"(Compute the divergence of a displacement field. + +The divergence of a vector field + +.. math:: + f(\boldsymbol{x}) = + (f_1(\boldsymbol{x}), + f_2(\boldsymbol{x}), + f_3(\boldsymbol{x})) + +with :math:`\boldsymbol{x} = (x_1, x_2, x_3)` +is defined as + +.. math:: + \nabla \cdot f (\boldsymbol{x}) = + \sum_{i=1}^3 + \frac{\partial f_i}{\partial x_i} (\boldsymbol{x}) + +.. note:: + All the arrays must be C-contiguous. + +Parameters +---------- +displacement: np.ndarray + Displacement field used to resample the image. + +origin: np.ndarray + Origin of the displacement field. + +spacing: np.ndarray + Spacing of the displacement field. + +direction: Tuple[Int] + Cosine direction matrix of the displacement field. + +Returns +------- +np.ndarray + Scalar volume image containing the divergence of + the input displacement. +)"; + + + +std::string rotor_docstring = +R"(Compute the rotor of a displacement field. + +The rotor of a 3D 3-vector field + +.. math:: + f(\boldsymbol{x}) = + (f_1(\boldsymbol{x}), + f_2(\boldsymbol{x}), + f_3(\boldsymbol{x})) + +with :math:`\boldsymbol{x} = (x_1, x_2, x_3)` +is defined as + +.. math:: + \nabla \times f(\boldsymbol{x}) = + \left( + \frac{\partial f_3}{\partial x_2} - + \frac{\partial f_2}{\partial x_3}, + \frac{\partial f_1}{\partial x_3} - + \frac{\partial f_3}{\partial x_1}, + \frac{\partial f_2}{\partial x_1} - + \frac{\partial f_1}{\partial x_2} + \right) + +.. note:: + All the arrays must be C-contiguous. + +Parameters +---------- +displacement: stk.Volume + Displacement field used to resample the image. + +Returns +------- +stk.Volume + Vector volume image containing the rotor of + the input displacement. +)"; + + + +std::string circulation_density_docstring = +R"(Compute the circulation density of a displacement field. + +The circulation density for a 3D 3-vector field is defined as the +norm of the rotor. + +.. note:: + All the arrays must be C-contiguous. + +Parameters +---------- +displacement: stk.Volume + Displacement field used to resample the image. + +Returns +------- +stk.Volume + Vector volume image containing the circulation + density of the input displacement. +)"; + +PYBIND11_MODULE(_stk, m) +{ + py::enum_(m, "Type") + .value("Unknown", stk::Type_Unknown) + .value("Char", stk::Type_Char) + .value("Char2", stk::Type_Char2) + .value("Char3", stk::Type_Char3) + .value("Char4", stk::Type_Char4) + .value("UChar", stk::Type_UChar) + .value("UChar2", stk::Type_UChar2) + .value("UChar3", stk::Type_UChar3) + .value("UChar4", stk::Type_UChar4) + .value("Short", stk::Type_Short) + .value("Short2", stk::Type_Short2) + .value("Short3", stk::Type_Short3) + .value("Short4", stk::Type_Short4) + .value("UShort", stk::Type_UShort) + .value("UShort2", stk::Type_UShort2) + .value("UShort3", stk::Type_UShort3) + .value("UShort4", stk::Type_UShort4) + .value("Int", stk::Type_Int) + .value("Int2", stk::Type_Int2) + .value("Int3", stk::Type_Int3) + .value("Int4", stk::Type_Int4) + .value("UInt", stk::Type_UInt) + .value("UInt2", stk::Type_UInt2) + .value("UInt3", stk::Type_UInt3) + .value("UInt4", stk::Type_UInt4) + .value("Float", stk::Type_Float) + .value("Float2", stk::Type_Float2) + .value("Float3", stk::Type_Float3) + .value("Float4", stk::Type_Float4) + .value("Double", stk::Type_Double) + .value("Double2", stk::Type_Double2) + .value("Double3", stk::Type_Double3) + .value("Double4", stk::Type_Double4) + .export_values(); + + py::class_(m, "Volume", py::buffer_protocol()) + .def_buffer(get_buffer_info) + .def(py::init<>()) + .def(py::init(&make_volume), + py::arg("array"), + py::arg("origin") = float3{0.0, 0.0, 0.0}, + py::arg("spacing") = float3{1.0, 1.0, 1.0}, + py::arg("direction") = Matrix3x3f::Identity + ) + .def("valid", &stk::Volume::valid) + .def("copy_meta_from", &stk::Volume::copy_meta_from) + .def_property_readonly("size", &stk::Volume::size) + .def_property_readonly("type", &stk::Volume::voxel_type) + .def_property("origin", &stk::Volume::origin, &stk::Volume::set_origin) + .def_property("spacing", &stk::Volume::spacing, &stk::Volume::set_spacing) + .def_property("direction", &stk::Volume::direction, &stk::Volume::set_direction) + ; + + py::implicitly_convertible(); + + m.def("read_volume", &stk::read_volume, ""); + m.def("write_volume", &stk::write_volume, ""); + + m.def("divergence", + &stk::divergence, + divergence_docstring.c_str() + ); + + m.def("rotor", + &stk::rotor, + rotor_docstring.c_str() + ); + + m.def("circulation_density", + &stk::circulation_density, + circulation_density_docstring.c_str() + ); + + // Translate relevant exception types. The exceptions not handled + // here will be translated autmatically according to pybind11's rules. + py::register_exception_translator([](std::exception_ptr p) { + try { + if (p) { + std::rethrow_exception(p); + } + } + catch (const stk::FatalException &e) { + // Map stk::FatalException to Python RutimeError + // +13 to remove the "Fatal error: " from the message + PyErr_SetString(PyExc_RuntimeError, e.what() + 13); + } + }); +} + diff --git a/src/stk/CMakeLists.txt b/src/stk/CMakeLists.txt new file mode 100644 index 0000000..7e0bac3 --- /dev/null +++ b/src/stk/CMakeLists.txt @@ -0,0 +1,105 @@ +set(STK_SRCS + "common/error.cpp" + "common/error.h" + "common/log.cpp" + "common/log.h" + "common/stack_trace.cpp" + "common/stack_trace.h" + "cuda/volume.h" + "filters/decomposable_filter.h" + "filters/decomposable_filter.inl" + "filters/gaussian_filter.cpp" + "filters/gaussian_filter.h" + "filters/map.h" + "filters/sobel.cpp" + "filters/sobel.h" + "filters/vector_calculus.h" + "filters/vector_calculus.inl" + "io/io.cpp" + "io/io.h" + "io/nifti.cpp" + "io/nifti.h" + "io/nrrd.cpp" + "io/nrrd.h" + "io/vtk.cpp" + "io/vtk.h" + "math/float2.h" + "math/float3.h" + "math/float4.h" + "math/math.cpp" + "math/math.h" + "math/types.h" + "image/dim3.h" + "image/types.cpp" + "image/types.h" + "image/volume.cpp" + "image/volume.h" + "image/volume.inl" +) + +if (STK_USE_CUDA) + set(STK_SRCS + ${STK_SRCS} + "cuda/cuda.cpp" + "cuda/cuda.h" + "cuda/stream.cpp" + "cuda/stream.h" + "filters/gpu/normalize.cu" + "image/gpu_volume.cpp" + "image/gpu_volume.cu" + "image/gpu_volume.h" + ) +endif() + +if(STK_ITK_BRIDGE) + set(STK_SRCS + ${STK_SRCS} + "itk_bridge/itk_convert.h" + "itk_bridge/itk_io.cpp" + "itk_bridge/itk_io.h" + ) +endif() + +add_library(stk STATIC ${STK_SRCS}) + +target_include_directories(stk PUBLIC + "${CMAKE_CURRENT_SOURCE_DIR}/.." + ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES} + ) + +if(NIFTI_FOUND) + target_include_directories(stk PUBLIC ${NIFTI_INCLUDE_DIR}) + target_link_libraries(stk PUBLIC ${NIFTI_LIBRARIES}) +else() + target_link_libraries(stk PUBLIC niftiio) +endif() + +if(NrrdIO_FOUND) + target_include_directories(stk PUBLIC ${NrrdIO_INCLUDE_DIR}) + target_link_libraries(stk PUBLIC ${NrrdIO_LIBRARIES}) +else() + target_link_libraries(stk PUBLIC NrrdIO) +endif() + +if(STK_ITK_BRIDGE) + add_definitions(-DSTK_ITK_BRIDGE) + find_package(ITK REQUIRED) + include(${ITK_USE_FILE}) + + target_include_directories(stk PRIVATE ${ITK_INCLUDE_DIRS}) + target_link_libraries(stk PRIVATE ${ITK_LIBRARIES}) +endif() + +if(STK_STACK_TRACE) + add_definitions(-DSTK_STACK_TRACE) + + if(UNIX) + find_package(Backward) + if(Backward_FOUND) + add_backward(stk) + target_link_libraries(stk PUBLIC ${BACKWARD_LIBRARIES}) + add_definitions(-DSTK_USE_BACKWARD) + endif() + endif() +endif() + diff --git a/src/stk/image/gpu_volume.cpp b/src/stk/image/gpu_volume.cpp index 5bc1074..54a0ec0 100644 --- a/src/stk/image/gpu_volume.cpp +++ b/src/stk/image/gpu_volume.cpp @@ -449,11 +449,6 @@ void GpuVolume::set_direction(const Matrix3x3f& direction) _direction = direction; _inverse_direction = _direction.inverse(); } -void GpuVolume::set_direction(const std::initializer_list direction) -{ - _direction.set(direction); - _inverse_direction = _direction.inverse(); -} const float3& GpuVolume::origin() const { diff --git a/src/stk/image/gpu_volume.h b/src/stk/image/gpu_volume.h index 857ddab..45d63bf 100644 --- a/src/stk/image/gpu_volume.h +++ b/src/stk/image/gpu_volume.h @@ -110,7 +110,6 @@ namespace stk void set_origin(const float3& origin); void set_spacing(const float3& spacing); void set_direction(const Matrix3x3f& direction); - void set_direction(const std::initializer_list direction); const float3& origin() const; const float3& spacing() const; diff --git a/src/stk/image/volume.cpp b/src/stk/image/volume.cpp index 28085cf..a22a534 100644 --- a/src/stk/image/volume.cpp +++ b/src/stk/image/volume.cpp @@ -255,11 +255,6 @@ void Volume::set_direction(const Matrix3x3f& direction) _direction = direction; _inverse_direction = _direction.inverse(); } -void Volume::set_direction(const std::initializer_list direction) -{ - _direction.set(direction); - _inverse_direction = _direction.inverse(); -} const float3& Volume::origin() const { return _origin; diff --git a/src/stk/image/volume.h b/src/stk/image/volume.h index 9950136..a09f2ce 100644 --- a/src/stk/image/volume.h +++ b/src/stk/image/volume.h @@ -150,7 +150,6 @@ namespace stk void set_origin(const float3& origin); void set_spacing(const float3& spacing); void set_direction(const Matrix3x3f& direction); - void set_direction(const std::initializer_list direction); const float3& origin() const; const float3& spacing() const; diff --git a/src/stk/math/math.cpp b/src/stk/math/math.cpp index 66b41cd..e75cd22 100644 --- a/src/stk/math/math.cpp +++ b/src/stk/math/math.cpp @@ -1,5 +1,10 @@ #include "math.h" +Matrix3x3f Matrix3x3f::Identity { + float3{1, 0, 0}, + float3{0, 1, 0}, + float3{0, 0, 1} +}; namespace stk { diff --git a/src/stk/math/types.h b/src/stk/math/types.h index 0d91dbe..8a12fff 100644 --- a/src/stk/math/types.h +++ b/src/stk/math/types.h @@ -136,6 +136,8 @@ struct STK_ALIGN(16) double4 struct Matrix3x3f { + static Matrix3x3f Identity; + enum : unsigned int { rows = 3, cols = 3 diff --git a/stk/__init__.py b/stk/__init__.py new file mode 100644 index 0000000..a5a299c --- /dev/null +++ b/stk/__init__.py @@ -0,0 +1 @@ +from _stk import * \ No newline at end of file diff --git a/test/test_volume.py b/test/test_volume.py new file mode 100644 index 0000000..3f8bade --- /dev/null +++ b/test/test_volume.py @@ -0,0 +1,85 @@ +import os +import random +import unittest +import shutil, tempfile + +import numpy as np + +import stk + +class Test_Volume(unittest.TestCase): + def setUp(self): + self.temp_dir = tempfile.mkdtemp() + + def tearDown(self): + shutil.rmtree(self.temp_dir) + + def test_constructor(self): + arr = np.zeros((4, 5, 6)).astype(np.int32) + vol = stk.Volume(arr) + + np.testing.assert_equal(vol.size, (6, 5, 4)) + np.testing.assert_equal(vol.origin, (0, 0, 0)) + np.testing.assert_equal(vol.spacing, (1, 1, 1)) + np.testing.assert_equal(vol.direction, ((1, 0, 0), (0, 1, 0), (0, 0, 1))) + np.testing.assert_equal(np.array(vol), arr) + + d = ((10, 0, 0), (0, 11, 0), (0, 0, 12)) + vol = stk.Volume(arr, origin=(1,2,3), spacing=(4,5,6), direction=d) + + np.testing.assert_equal(vol.size, (6, 5, 4)) + np.testing.assert_equal(vol.origin, (1, 2, 3)) + np.testing.assert_equal(vol.spacing, (4, 5, 6)) + np.testing.assert_equal(vol.direction, d) + np.testing.assert_equal(np.array(vol), arr) + + def test_meta(self): + vol = stk.Volume() + vol.origin = (2, 3, 4) + np.testing.assert_equal(vol.origin, (2, 3, 4)) + + vol = stk.Volume() + vol.spacing = (2, 3, 4) + np.testing.assert_equal(vol.spacing, (2, 3, 4)) + + d = ((10, 0, 0), (0, 11, 0), (0, 0, 12)) + vol = stk.Volume() + vol.direction = d + np.testing.assert_equal(vol.direction, d) + + def test_copy_meta_from(self): + vol1 = stk.Volume() + vol2 = stk.Volume() + + vol1.origin = (2,3,4) + vol1.spacing = (5,6,7) + d = ((10, 0, 0), (0, 11, 0), (0, 0, 12)) + vol1.direction = d + + vol2.copy_meta_from(vol1) + + np.testing.assert_equal(vol2.origin, (2, 3, 4)) + np.testing.assert_equal(vol2.spacing, (5, 6, 7)) + np.testing.assert_equal(vol2.direction, d) + + def test_type(self): + data = np.zeros((5,5,5)).astype(np.float32) + vol = stk.Volume(data) + self.assertEqual(vol.type, stk.Type.Float) + + data = np.zeros((5,5,5,3)).astype(np.float32) + vol = stk.Volume(data) + self.assertEqual(vol.type, stk.Type.Float3) + + def test_io(self): + data = np.zeros((5,5,5)).astype(np.float32) + vol1 = stk.Volume(data, origin=(2,3,4), spacing=(5,6,7)) + + f = os.path.join(self.temp_dir, 'test.nrrd') + stk.write_volume(f, vol1) + + vol2 = stk.read_volume(f) + np.testing.assert_equal(vol2.origin, vol1.origin) + np.testing.assert_equal(vol2.spacing, vol1.spacing) + np.testing.assert_equal(vol2.direction, vol1.direction) + np.testing.assert_equal(np.array(vol2), np.array(vol1)) diff --git a/third_party/CMakeLists.txt b/third_party/CMakeLists.txt index 18daeee..62d76b8 100644 --- a/third_party/CMakeLists.txt +++ b/third_party/CMakeLists.txt @@ -61,3 +61,7 @@ if(NOT NrrdIO_FOUND) endif() endif() +if(STK_BUILD_PYTHON_WRAPPER) + add_subdirectory(pybind11) +endif() + diff --git a/third_party/pybind11 b/third_party/pybind11 new file mode 160000 index 0000000..7ec2ddf --- /dev/null +++ b/third_party/pybind11 @@ -0,0 +1 @@ +Subproject commit 7ec2ddfc95f65d1e986d359466a6c254aa514ef3