diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..75e2e17 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,89 @@ +cmake_minimum_required(VERSION 3.21...4.0) +project(dmri-commit CXX C) +find_package(Python COMPONENTS Interpreter Development.Module NumPy REQUIRED) +include(UseCython) + +set(CYTHON_ARGS "-3") + +execute_process( + COMMAND "${Python_EXECUTABLE}" -c + "import sys; sys.path.insert(0,'${CMAKE_CURRENT_SOURCE_DIR}'); from setup_operator import write_operator_c_file; write_operator_c_file()" + WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}" + COMMAND_ERROR_IS_FATAL ANY +) + +# trx-cpp via FetchContent +set(CMAKE_POSITION_INDEPENDENT_CODE ON) +set(CMAKE_POLICY_DEFAULT_CMP0077 NEW) +include(FetchContent) + +set(TRX_BUILD_TESTS OFF) +set(TRX_BUILD_EXAMPLES OFF) +set(TRX_BUILD_DOCS OFF) +set(TRX_BUILD_BENCHMARKS OFF) +set(TRX_ENABLE_NIFTI OFF) +set(TRX_ENABLE_INSTALL OFF) + +FetchContent_Declare( + trx-cpp + GIT_REPOSITORY https://github.com/tee-ar-ex/trx-cpp.git + GIT_TAG main + GIT_SHALLOW TRUE +) +FetchContent_MakeAvailable(trx-cpp) + +# json11.cpp is missing '#include '; force-include it +if(MSVC) + target_compile_options(trx PRIVATE /FIcstdint) +else() + target_compile_options(trx PRIVATE -include cstdint) +endif() + +cython_transpile(commit/trk2dictionary/trk2dictionary.pyx LANGUAGE CXX OUTPUT_VARIABLE trk2dict_src) +Python_add_library(trk2dictionary MODULE "${trk2dict_src}" WITH_SOABI) +target_compile_features(trk2dictionary PRIVATE cxx_std_17) +target_compile_options(trk2dictionary PRIVATE $<$>:-w>) +target_include_directories(trk2dictionary PRIVATE + "${CMAKE_CURRENT_SOURCE_DIR}/commit/trk2dictionary") +target_link_libraries(trk2dictionary PRIVATE trx-cpp::trx Python::NumPy) +install(TARGETS trk2dictionary DESTINATION commit) + +cython_transpile(commit/core.pyx LANGUAGE CXX OUTPUT_VARIABLE core_src) +Python_add_library(core MODULE "${core_src}" WITH_SOABI) +target_compile_features(core PRIVATE cxx_std_11) +target_compile_options(core PRIVATE $<$>:-w>) +target_link_libraries(core PRIVATE Python::NumPy) +install(TARGETS core DESTINATION commit) + +cython_transpile(commit/proximals.pyx LANGUAGE CXX OUTPUT_VARIABLE proximals_src) +Python_add_library(proximals MODULE "${proximals_src}" WITH_SOABI) +target_compile_features(proximals PRIVATE cxx_std_11) +target_compile_options(proximals PRIVATE $<$>:-w>) +target_link_libraries(proximals PRIVATE Python::NumPy) +install(TARGETS proximals DESTINATION commit) + +cython_transpile(commit/models.pyx LANGUAGE CXX OUTPUT_VARIABLE models_src) +Python_add_library(models MODULE "${models_src}" WITH_SOABI) +target_compile_features(models PRIVATE cxx_std_11) +target_compile_options(models PRIVATE $<$>:-w>) +target_link_libraries(models PRIVATE Python::NumPy) +install(TARGETS models DESTINATION commit) + +cython_transpile(commit/operator/operator.pyx LANGUAGE C OUTPUT_VARIABLE operator_src) +Python_add_library(commit_operator MODULE + "${operator_src}" + "${CMAKE_CURRENT_SOURCE_DIR}/commit/operator/operator_c.c" + WITH_SOABI +) +set_target_properties(commit_operator PROPERTIES OUTPUT_NAME operator) +target_compile_options(commit_operator PRIVATE + $<$>:-w -O3 -Ofast> + $<$:/fp:fast /DHAVE_STRUCT_TIMESPEC> +) +target_link_libraries(commit_operator PRIVATE Python::NumPy) +if(WIN32) + target_include_directories(commit_operator PRIVATE "$ENV{PTHREAD_WIN_INCLUDE}") + target_link_directories(commit_operator PRIVATE "$ENV{PTHREAD_WIN_LIB}") + target_link_libraries(commit_operator PRIVATE pthread) +endif() +install(TARGETS commit_operator DESTINATION commit/operator) diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index 072968d..e335727 100644 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/commit/trk2dictionary/trk2dictionary.pyx @@ -18,6 +18,7 @@ import pickle from importlib import metadata import shutil import time +from trx.io import load as load_trx logger = setup_logger('trk2dictionary') @@ -340,14 +341,28 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam if not exists(filename_tractogram): logger.error( f'Tractogram file not found: {filename_tractogram}' ) extension = splitext(filename_tractogram)[1] - if extension != ".trk" and extension != ".tck": - logger.error( 'Invalid input file: only .trk and .tck are supported' ) - - hdr = nibabel.streamlines.load( filename_tractogram, lazy_load=True ).header - - - if extension == ".trk": - logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) + if extension != ".trx" and extension != ".trk" and extension != ".tck": + logger.error( 'Invalid input file: only .trx, .trk, and .tck are supported') + + if extension == ".trx": + hdr = load_trx(filename_tractogram, "same").header + affine = hdr["VOXEL_TO_RASMM"] + voxel_sizes = nibabel.affines.voxel_sizes(affine) + + Nx = int(hdr['DIMENSIONS'][0]) + Ny = int(hdr['DIMENSIONS'][1]) + Nz = int(hdr['DIMENSIONS'][2]) + Px = voxel_sizes[0] + Py = voxel_sizes[1] + Pz = voxel_sizes[2] + + data_offset = 0 # stored separately in .trx + n_count = hdr['NB_STREAMLINES'] + n_scalars = 0 # stored separately in .trx + n_properties = 0 # stored separately in .trx + elif extension == ".trk": + hdr = nibabel.streamlines.load( filename_tractogram, lazy_load=True ).header + Nx = int(hdr['dimensions'][0]) Ny = int(hdr['dimensions'][1]) Nz = int(hdr['dimensions'][2]) @@ -369,6 +384,7 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam else: logger.error( 'TCK files do not contain information about the geometry. Use "TCK_ref_image" for that' ) logger.subinfo ( f'geometry taken from "{TCK_ref_image}"', indent_lvl=3, indent_char='-' ) + hdr = nibabel.streamlines.load( filename_tractogram, lazy_load=True ).header niiREF = nibabel.load( TCK_ref_image ) niiREF_hdr = _get_header( niiREF ) @@ -408,8 +424,11 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam # get toVOXMM matrix (remove voxel scaling from affine) in case of TCK cdef float [:] toVOXMM cdef float* ptrToVOXMM - if extension == ".tck": - M = _get_affine( niiREF ).copy() + if extension == ".tck" or extension == ".trx": + if extension == ".tck": + M = _get_affine( niiREF ).copy() + else: + M = np.asarray(affine, dtype=np.float64).copy() # float64 conversion added to comply with the new cast policy of numpy v2 M[:3, :3] = M[:3, :3].dot( np.diag([np.float64(1)/Px,np.float64(1)/Py,np.float64(1)/Pz]) ) toVOXMM = np.ravel(np.linalg.inv(M)).astype(' #include #include +#include #include #include +#include #define _FILE_OFFSET_BITS 64 #define MAX_FIB_LEN 10000 @@ -88,12 +90,13 @@ int verbosity = 0; bool rayBoxIntersection( Vector& origin, Vector& direction, Vector& vmin, Vector& vmax, double & t); void fiberForwardModel( float fiber[3][MAX_FIB_LEN], unsigned int pts, int nReplicas, double* ptrBlurRho, double* ptrBlurAngle, double* ptrBlurWeights, bool doApplyBlur, short* ptrHashTable, vector>& P ); void segmentForwardModel( const Vector& P1, const Vector& P2, int k, double w, short* ptrHashTable); +unsigned int read_fiberTRX(const trx::AnyTrxFile& fp, float fiber[3][MAX_FIB_LEN], int idx, float* ptrToVOXMM); unsigned int read_fiberTRK( FILE* fp, float fiber[3][MAX_FIB_LEN], int ns, int np ); unsigned int read_fiberTCK( FILE* fp, float fiber[3][MAX_FIB_LEN] , float* toVOXMM ); // ---------- Parallel fuction -------------- -int ICSegments( char* str_filename, int isTRK, int n_count, int nReplicas, int n_scalars, int n_properties, float* ptrToVOXMM, +int ICSegments( char* str_filename, const trx::AnyTrxFile& trxFile, int isTRX, int isTRK, int n_count, int nReplicas, int n_scalars, int n_properties, float* ptrToVOXMM, double* ptrTDI , double* ptrBlurRho, double* ptrBlurAngle, double* ptrBlurWeights, bool* ptrBlurApplyTo, short* ptrHashTable, char* path_out, unsigned long long int offset, int idx, unsigned int startpos, unsigned int endpos ); @@ -165,69 +168,81 @@ int trk2dictionary( // Check the file extension // ------------------------------------- int isTRK; + int isTRX; char *ext = strrchr(str_filename, '.'); - if (strcmp(ext,".trk")==0) // for .trk file + if (strcmp(ext,".trx")==0) { // for .trx file + isTRK = 0; + isTRX = 1; + } else if (strcmp(ext,".trk")==0) { // for .trk file isTRK = 1; - else if (strcmp(ext,".tck")==0) // for .tck file + isTRX = 0; + } else if (strcmp(ext,".tck")==0) { // for .tck file isTRK = 0; - else - return 0; + isTRX = 0; + } else return 0; // Open tractogram file and compute the offset for each thread + // This is only needed for .trk and .tck files + // For trx, we can use the trx::AnyTrxFile class to read + // the streamlines in parallel // ----------------------------------------------------------------- + trx::AnyTrxFile trxFile; unsigned long long int current; unsigned long long int *OffsetArr = new unsigned long long int[threads_count](); int f = 0; float *Buff = new float[3](); int N; - FILE* fpTractogram = fopen(str_filename,"rb"); - if (fpTractogram == NULL) return 0; - fseek( fpTractogram, data_offset, SEEK_SET ); // skip the header - - OffsetArr[0] = ftell( fpTractogram ); - - if(isTRK) { - while( f != n_count) { - fread( (char*)&N, 1, 4, fpTractogram ); // read the number of points in each streamline - if( N >= MAX_FIB_LEN || N <= 0 ) return 0; // check - for( int k=0; k= MAX_FIB_LEN || N <= 0 ) return 0; // check + for( int k=0; k& origin, Vector& direction, Vect } +unsigned int read_fiberTRX(const trx::AnyTrxFile& fp, float fiber[3][MAX_FIB_LEN], int idx, float* ptrToVOXMM) +{ + std::vector> streamline = fp.get_streamline(idx); + const unsigned int N = streamline.size(); + + if ( N >= MAX_FIB_LEN || N <= 0 ) + return 0; + + for (uint64_t i = 0; i < N; i++) { + const float x = static_cast(streamline[i][0]); + const float y = static_cast(streamline[i][1]); + const float z = static_cast(streamline[i][2]); + fiber[0][i] = x * ptrToVOXMM[0] + y * ptrToVOXMM[1] + z * ptrToVOXMM[2] + ptrToVOXMM[3]; + fiber[1][i] = x * ptrToVOXMM[4] + y * ptrToVOXMM[5] + z * ptrToVOXMM[6] + ptrToVOXMM[7]; + fiber[2][i] = x * ptrToVOXMM[8] + y * ptrToVOXMM[9] + z * ptrToVOXMM[10] + ptrToVOXMM[11]; + } + return N; +} + // Read a fiber from file .trk unsigned int read_fiberTRK( FILE* fp, float fiber[3][MAX_FIB_LEN], int ns, int np ) { diff --git a/pyproject.toml b/pyproject.toml index 4d6f5cb..b2ea870 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,10 +1,12 @@ [build-system] requires = [ + "scikit-build-core>=0.5", + "cmake>=3.21", "cython>=3.0.10", + "cython-cmake", "numpy>=1.24.4", - "setuptools>=69.2.0" ] -build-backend = "setuptools.build_meta" +build-backend = "scikit_build_core.build" [project] name = "dmri-commit" @@ -13,7 +15,8 @@ dependencies = [ "dmri-amico>=2.0.1", "dmri-dicelib>=1.1.0", "numpy>=1.24.4", - "scipy>=1.10.1" + "scipy>=1.10.1", + "trx-python>=0.1" ] requires-python = ">=3.8" authors = [ @@ -49,19 +52,8 @@ Repository = "https://github.com/daducci/COMMIT.git" Issues = "https://github.com/daducci/COMMIT/issues" Changelog = "https://github.com/daducci/COMMIT/blob/master/CHANGELOG.md" -[tool.setuptools] -packages = [ - "commit", - "commit.operator" -] -include-package-data = false - -[tool.setuptools.package-data] -"commit.operator" = [ - "*.c", - "*.pyx", - "*.pyxbld" -] +[tool.scikit-build] +wheel.packages = ["commit"] [tool.cibuildwheel] build-verbosity = 3 diff --git a/requirements.txt b/requirements.txt index 1bad000..9972e72 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1,10 @@ wheel setuptools>=46.1 Cython>=0.29 +cmake>=3.16 numpy>=1.12 scipy>=1.0 dipy>=1.0 dmri-dicelib==1.1.0 dmri-amico>=2.0.1 +trx-python>=0.1 diff --git a/setup.py b/setup.py deleted file mode 100644 index 587d86d..0000000 --- a/setup.py +++ /dev/null @@ -1,84 +0,0 @@ -import os -import sys -from setuptools import Extension, setup -from setuptools.command.build_ext import build_ext - -# name of the package -package_name = 'commit' - -def get_extensions(): - # Cython extension to create the sparse data structure from a tractogram - # for the computation of matrix-vector multiplications - trk2dictionary = Extension(name=f'{package_name}.trk2dictionary', - sources=[f'{package_name}/trk2dictionary/trk2dictionary.pyx'], - libraries=[] if sys.platform == 'win32' else ['stdc++'], - extra_compile_args=[] if sys.platform == 'win32' else ['-w', '-std=c++11'], - language='c++') - core = Extension(name=f'{package_name}.core', - sources=[f'{package_name}/core.pyx'], - extra_compile_args=[] if sys.platform == 'win32' else ['-w', '-std=c++11'], - language='c++') - proximals = Extension(name=f'{package_name}.proximals', - sources=[f'{package_name}/proximals.pyx'], - extra_compile_args=[] if sys.platform == 'win32' else ['-w', '-std=c++11'], - language='c++') - models = Extension(name=f'{package_name}.models', - sources=[f'{package_name}/models.pyx'], - extra_compile_args=[] if sys.platform == 'win32' else ['-w', '-std=c++11'], - language='c++') - # NOTE: Windows requires the pthread-win32 static library to compile the operator extension - # The library can be downloaded from https://github.com/GerHobbelt/pthread-win32 - # The PTHREAD_WIN_INCLUDE and PTHREAD_WIN_LIB environment variables must be set to the include and lib directories - if sys.platform == 'win32': - try: - pthread_win_include = os.environ['PTHREAD_WIN_INCLUDE'] - pthread_win_lib = os.environ['PTHREAD_WIN_LIB'] - except KeyError: - raise RuntimeError('PTHREAD_WIN_INCLUDE and PTHREAD_WIN_LIB must be set') - operator = Extension(name=f'{package_name}.operator.operator', - sources=[f'{package_name}/operator/operator.pyx', f'{package_name}/operator/operator_c.c'], - include_dirs=[pthread_win_include] if sys.platform == 'win32' else [], - libraries=['pthread'] if sys.platform == 'win32' else [], - library_dirs=[pthread_win_lib] if sys.platform == 'win32' else [], - extra_compile_args=['/fp:fast', '/DHAVE_STRUCT_TIMESPEC'] if sys.platform == 'win32' else ['-w', '-O3', '-Ofast'], - language='c') - - return [trk2dictionary, core, proximals, operator, models] - -class CustomBuildExtCommand(build_ext): - """ build_ext command to use when numpy headers are needed. """ - def run(self): - # Now that the requirements are installed, get everything from numpy - from multiprocessing import cpu_count - - from Cython.Build import cythonize - from numpy import get_include - - # Add everything requires for build - self.swig_opts = None - self.include_dirs.extend([get_include()]) - self.distribution.ext_modules[:] = cythonize( self.distribution.ext_modules, build_dir='build' ) - - # if not specified via '-j N' option, set compilation using max number of cores - if self.parallel is None: - self.parallel = cpu_count() - print( f'Parallel compilation using {self.parallel} threads' ) - - # Call original build_ext command - build_ext.finalize_options(self) - build_ext.run(self) - -# generate the operator_c.c file -sys.path.insert(0, os.path.dirname(__file__)) -from setup_operator import write_operator_c_file -write_operator_c_file() - -# create the 'build' directory -if not os.path.exists('build'): - os.makedirs('build') - -# install the package -setup( - cmdclass={'build_ext': CustomBuildExtCommand}, - ext_modules=get_extensions() -)