From f65d256a396124b10221b0a8a03b2214e2e073b8 Mon Sep 17 00:00:00 2001 From: 36000 Date: Thu, 18 Jun 2026 16:18:30 -0700 Subject: [PATCH 1/4] First Draft of TRX support --- CMakeLists.txt | 89 ++++++++++++++ commit/trk2dictionary/trk2dictionary.pyx | 28 ++++- commit/trk2dictionary/trk2dictionary_c.cpp | 136 ++++++++++++++------- pyproject.toml | 24 ++-- requirements.txt | 2 + setup.py | 84 ------------- 6 files changed, 211 insertions(+), 152 deletions(-) create mode 100644 CMakeLists.txt delete mode 100644 setup.py 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..9999c9d 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,30 @@ 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 != ".trx" and extension != ".trk" and extension != ".tck": + logger.error( 'Invalid input file: only .trx, .trk, and .tck are supported') + if extension == ".trx": + logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) + 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[1] + Py = voxel_sizes[2] + Pz = voxel_sizes[3] + + 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 if extension == ".trk": logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) + 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 +386,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 ) diff --git a/commit/trk2dictionary/trk2dictionary_c.cpp b/commit/trk2dictionary/trk2dictionary_c.cpp index db0928f..5f945e8 100644 --- a/commit/trk2dictionary/trk2dictionary_c.cpp +++ b/commit/trk2dictionary/trk2dictionary_c.cpp @@ -11,6 +11,8 @@ #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, 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,17 +168,23 @@ 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 // ----------------------------------------------------------------- unsigned long long int current; unsigned long long int *OffsetArr = new unsigned long long int[threads_count](); @@ -183,51 +192,53 @@ int trk2dictionary( 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 fpTractogram1; + + if ( isTRX ) { + fpTractogram1 = trx::AnyTrxFile::load(str_filename); + } else { + FILE* fp = fopen( str_filename, "rb" ); + if ( fp == NULL ) return 0; + fseek(fp, offset, SEEK_SET); + fpTractogram1 = fp; + } tempTotFibers = 0; temp_totICSegments = 0; @@ -499,11 +517,12 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp for(int f=startpos; f(fpTractogram1), fiber, idx, ptrToVOXMM ); + else if ( isTRK ) + N = read_fiberTRK( get(fpTractogram1), fiber, n_scalars, n_properties ); else - N = read_fiberTCK( fpTractogram1, fiber , ptrToVOXMM ); - + N = read_fiberTCK( get(fpTractogram1), fiber, ptrToVOXMM ); fiberForwardModel( fiber, N, nReplicas, ptrBlurRho, ptrBlurAngle, ptrBlurWeights, ptrBlurApplyTo[f], ptrHashTable, P ); @@ -564,7 +583,10 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp } } } - fclose( fpTractogram1 ); + if ( isTRX ) + get(fpTractogram1).close(); + else + fclose( get(fpTractogram1) ); fclose( pDict_TRK_norm ); fclose( pDict_IC_f ); fclose( pDict_IC_v ); @@ -869,6 +891,26 @@ bool rayBoxIntersection( Vector& origin, Vector& direction, Vect } +unsigned int read_fiberTRX(const trx::AnyTrxFile &fp, float fiber[3][MAX_FIB_LEN], int idx, float* ptrToVOXMM) +{ + const uint64_t start = fp.offsets_u64[idx]; + const uint64_t end = fp.offsets_u64[idx + 1]; + const unsigned int N = static_cast(end - start); + + std::vector> streamline = fp.get_streamline(idx); + + for (uint64_t i = start; i < end; ++i) { + const size_t local = static_cast(i - start); + const float x = static_cast(streamline[local][0]); + const float y = static_cast(streamline[local][1]); + const float z = static_cast(streamline[local][2]); + fiber[0][local] = x * ptrToVOXMM[0] + y * ptrToVOXMM[1] + z * ptrToVOXMM[2] + ptrToVOXMM[3]; + fiber[1][local] = x * ptrToVOXMM[4] + y * ptrToVOXMM[5] + z * ptrToVOXMM[6] + ptrToVOXMM[7]; + fiber[2][local] = 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() -) From b02b3fc82bb7b9ecf5aae5be2f1f51d53b63ebab Mon Sep 17 00:00:00 2001 From: 36000 Date: Thu, 18 Jun 2026 16:30:03 -0700 Subject: [PATCH 2/4] typo fixes in trx --- commit/trk2dictionary/trk2dictionary.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index 9999c9d..3a92fc2 100644 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/commit/trk2dictionary/trk2dictionary.pyx @@ -353,15 +353,15 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam Nx = int(hdr['DIMENSIONS'][0]) Ny = int(hdr['DIMENSIONS'][1]) Nz = int(hdr['DIMENSIONS'][2]) - Px = voxel_sizes[1] - Py = voxel_sizes[2] - Pz = voxel_sizes[3] + 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 - if extension == ".trk": + elif extension == ".trk": logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) hdr = nibabel.streamlines.load( filename_tractogram, lazy_load=True ).header From 4f7d90667164a569e7b48829594c9257182a9bf8 Mon Sep 17 00:00:00 2001 From: 36000 Date: Mon, 22 Jun 2026 22:28:21 -0700 Subject: [PATCH 3/4] fix trx affine error, load TRX file outside of parallelization --- commit/trk2dictionary/trk2dictionary.pyx | 7 ++- commit/trk2dictionary/trk2dictionary_c.cpp | 64 +++++++++++----------- 2 files changed, 38 insertions(+), 33 deletions(-) diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index 3a92fc2..9ad0a44 100644 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/commit/trk2dictionary/trk2dictionary.pyx @@ -426,8 +426,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 #include #define _FILE_OFFSET_BITS 64 @@ -90,13 +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_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 isTRX, 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 ); @@ -185,15 +185,19 @@ int trk2dictionary( // 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; - if (!isTRX) - { + if (isTRX) { + trxFile = trx::AnyTrxFile::load(str_filename); + } else { FILE* fpTractogram = fopen(str_filename,"rb"); if (fpTractogram == NULL) return 0; fseek( fpTractogram, data_offset, SEEK_SET ); // skip the header @@ -254,7 +258,7 @@ int trk2dictionary( } // ---- Original ------ for( int i = 0; i fpTractogram1; + FILE* fpTractogram1; - if ( isTRX ) { - fpTractogram1 = trx::AnyTrxFile::load(str_filename); - } else { + if ( !isTRX ) { FILE* fp = fopen( str_filename, "rb" ); if ( fp == NULL ) return 0; fseek(fp, offset, SEEK_SET); @@ -518,11 +523,11 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp { if ( isTRX ) - N = read_fiberTRX( get(fpTractogram1), fiber, idx, ptrToVOXMM ); + N = read_fiberTRX( trxFile, fiber, idx, ptrToVOXMM ); else if ( isTRK ) - N = read_fiberTRK( get(fpTractogram1), fiber, n_scalars, n_properties ); + N = read_fiberTRK( fpTractogram1, fiber, n_scalars, n_properties ); else - N = read_fiberTCK( get(fpTractogram1), fiber, ptrToVOXMM ); + N = read_fiberTCK( fpTractogram1, fiber, ptrToVOXMM ); fiberForwardModel( fiber, N, nReplicas, ptrBlurRho, ptrBlurAngle, ptrBlurWeights, ptrBlurApplyTo[f], ptrHashTable, P ); @@ -583,10 +588,8 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp } } } - if ( isTRX ) - get(fpTractogram1).close(); - else - fclose( get(fpTractogram1) ); + if ( !isTRX ) + fclose( fpTractogram1 ); fclose( pDict_TRK_norm ); fclose( pDict_IC_f ); fclose( pDict_IC_v ); @@ -891,22 +894,21 @@ bool rayBoxIntersection( Vector& origin, Vector& direction, Vect } -unsigned int read_fiberTRX(const trx::AnyTrxFile &fp, float fiber[3][MAX_FIB_LEN], int idx, float* ptrToVOXMM) +unsigned int read_fiberTRX(const trx::AnyTrxFile& fp, float fiber[3][MAX_FIB_LEN], int idx, float* ptrToVOXMM) { - const uint64_t start = fp.offsets_u64[idx]; - const uint64_t end = fp.offsets_u64[idx + 1]; - const unsigned int N = static_cast(end - start); - 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 = start; i < end; ++i) { - const size_t local = static_cast(i - start); - const float x = static_cast(streamline[local][0]); - const float y = static_cast(streamline[local][1]); - const float z = static_cast(streamline[local][2]); - fiber[0][local] = x * ptrToVOXMM[0] + y * ptrToVOXMM[1] + z * ptrToVOXMM[2] + ptrToVOXMM[3]; - fiber[1][local] = x * ptrToVOXMM[4] + y * ptrToVOXMM[5] + z * ptrToVOXMM[6] + ptrToVOXMM[7]; - fiber[2][local] = x * ptrToVOXMM[8] + y * ptrToVOXMM[9] + z * ptrToVOXMM[10] + ptrToVOXMM[11]; + 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; } From 07b68bf819553dacf5364f53a5b6c685759433a0 Mon Sep 17 00:00:00 2001 From: 36000 Date: Mon, 22 Jun 2026 22:33:41 -0700 Subject: [PATCH 4/4] only mention reference file for tck, not trx --- commit/trk2dictionary/trk2dictionary.pyx | 2 -- 1 file changed, 2 deletions(-) diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index 9ad0a44..e335727 100644 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/commit/trk2dictionary/trk2dictionary.pyx @@ -345,7 +345,6 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam logger.error( 'Invalid input file: only .trx, .trk, and .tck are supported') if extension == ".trx": - logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) hdr = load_trx(filename_tractogram, "same").header affine = hdr["VOXEL_TO_RASMM"] voxel_sizes = nibabel.affines.voxel_sizes(affine) @@ -362,7 +361,6 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam n_scalars = 0 # stored separately in .trx n_properties = 0 # stored separately in .trx elif extension == ".trk": - logger.subinfo ( f'geometry taken from "{filename_tractogram}"', indent_lvl=3, indent_char='-' ) hdr = nibabel.streamlines.load( filename_tractogram, lazy_load=True ).header Nx = int(hdr['dimensions'][0])