From b362b631eaa63b6031acf413b4030279dbc23795 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Thu, 1 Feb 2024 09:26:00 -0500 Subject: [PATCH 01/19] [DOC] fix some link [ci skip] --- doc/conf.py | 2 +- doc/index.rst | 21 --------------------- doc/links_names.inc | 10 ++++------ 3 files changed, 5 insertions(+), 28 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index 196ead3594e..8f7eabad72f 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -402,7 +402,7 @@ "link_type": "external" }, { - "name": "Google supported DIPY through the Google Summer of Code Program during Summer 2015, 2016, 2018", + "name": "Google supported DIPY through the Google Summer of Code Program during multiple Summer (2015-2024)", "link": "https://summerofcode.withgoogle.com/", "link_type": "external" }, diff --git a/doc/index.rst b/doc/index.rst index 54f7dea1c91..7e884443f76 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -128,29 +128,8 @@ Announcements - :doc:`DIPY 1.6.0 ` released January 16, 2023. - - - See some of our :ref:`Past Announcements ` -******* -Support -******* - -We acknowledge support from the following organizations: - -- The department of Intelligent Systems Engineering of Indiana University. - -- The National Institute of Biomedical Imaging and Bioengineering, NIH. - -- The Gordon and Betty Moore Foundation and the Alfred P. Sloan Foundation, through the - University of Washington eScience Institute Data Science Environment. - -- Google supported DIPY through the Google Summer of Code Program during - Summer 2015, 2016 and 2018. - -- The International Neuroinformatics Coordination Facility. - .. This tree is helping populate the side navigation panel .. toctree:: diff --git a/doc/links_names.inc b/doc/links_names.inc index 740252e8050..18dc59e9088 100644 --- a/doc/links_names.inc +++ b/doc/links_names.inc @@ -22,7 +22,6 @@ .. _nibabel: https://nipy.org/nibabel .. _nibabel pypi: https://pypi.python.org/pypi/nibabel .. _nipy development guidelines: https://nipy.org/nipy/devel/guidelines/index.html -.. _buildbots: https://nipy.bic.berkeley.edu/builders .. _`dipy gitter`: https://gitter.im/dipy/dipy .. _neurostars: https://neurostars.org/ .. _h5py: https://www.h5py.org/ @@ -47,7 +46,6 @@ .. _LGPL: https://www.gnu.org/copyleft/lesser.html .. Working process -.. _pynifti: https://niftilib.sourceforge.net/pynifti/ .. _nifticlibs: https://nifti.nimh.nih.gov .. _nifti: https://nifti.nimh.nih.gov .. _`nipy launchpad`: https://launchpad.net/nipy @@ -57,13 +55,13 @@ .. _`dipy mailing list`: https://mail.python.org/mailman3/lists/dipy.python.org/ .. _`nipy bugs`: https://bugs.launchpad.net/nipy .. _pep8: https://www.python.org/dev/peps/pep-0008/ -.. _`numpy coding style`: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt +.. _`numpy coding style`: https://numpydoc.readthedocs.io/en/latest/format.html#docstring-standard .. _`python module path`: https://docs.python.org/tutorial/modules.html#the-module-search-path .. Code support stuff .. _pychecker: https://pychecker.sourceforge.net/ -.. _pylint: https://www.logilab.org/project/pylint -.. _pyflakes: https://divmod.org/trac/wiki/DivmodPyflakes +.. _pylint: https://pylint.readthedocs.io/en/latest/ +.. _pyflakes: https://github.com/PyCQA/pyflakes .. _virtualenv: https://pypi.python.org/pypi/virtualenv .. _git: https://git.or.cz/ .. _github: https://github.com @@ -160,7 +158,7 @@ .. _powershell: https://www.microsoft.com/powershell .. _msysgit: https://code.google.com/p/msysgit .. _putty: https://www.chiark.greenend.org.uk/~sgtatham/putty -.. _visualstudiobuildtools: https://landinghub.visualstudio.com/visual-cpp-build-tools +.. _visualstudiobuildtools: https://visualstudio.microsoft.com/downloads/?q=build+tools#build-tools-for-visual-studio-2022 .. Functional imaging labs .. _`functional imaging laboratory`: https://www.fil.ion.ucl.ac.uk From 916757ccfa99fb4506f689681d5de6a9064295e9 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Thu, 1 Feb 2024 09:44:33 -0500 Subject: [PATCH 02/19] apply arokem suggestion Co-authored-by: Ariel Rokem --- doc/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/conf.py b/doc/conf.py index 8f7eabad72f..850f449695d 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -402,7 +402,7 @@ "link_type": "external" }, { - "name": "Google supported DIPY through the Google Summer of Code Program during multiple Summer (2015-2024)", + "name": "Google supported DIPY through the Google Summer of Code Program (2015-2024)", "link": "https://summerofcode.withgoogle.com/", "link_type": "external" }, From c2a6c2fddb470baff178283fba7c87f72c844313 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 5 Feb 2024 14:41:16 +0000 Subject: [PATCH 03/19] Bump codecov/codecov-action from 3 to 4 Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 3 to 4. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v3...v4) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/test_template.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test_template.yml b/.github/workflows/test_template.yml index 8dfd1eab8cf..fe7b7797715 100644 --- a/.github/workflows/test_template.yml +++ b/.github/workflows/test_template.yml @@ -126,7 +126,7 @@ jobs: tools/ci/run_tests.sh fi - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 + uses: codecov/codecov-action@v4 if: ${{ fromJSON(env.COVERAGE) }} with: token: ${{ secrets.CODECOV_TOKEN }} From 1c1974d101ad9ab6bed8e3c6becb35b81c5c4ef6 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Tue, 6 Feb 2024 14:58:15 -0500 Subject: [PATCH 04/19] [OPT] enable openmp for macOS --- .github/workflows/nightly.yml | 2 ++ .github/workflows/test_template.yml | 15 ++++++--------- dipy/meson.build | 16 ++++++++++++++++ doc/devel/installation_from_source.rst | 12 +++++++++++- 4 files changed, 35 insertions(+), 10 deletions(-) diff --git a/.github/workflows/nightly.yml b/.github/workflows/nightly.yml index 07b1cd86dcf..826c08b9387 100644 --- a/.github/workflows/nightly.yml +++ b/.github/workflows/nightly.yml @@ -95,11 +95,13 @@ jobs: - name: Build the wheel run: python -m cibuildwheel --output-dir dist env: + CIBW_BEFORE_ALL_MACOS: "brew install llvm libomp" CIBW_BUILD: ${{ matrix.cibw_python }} CIBW_ARCHS_MACOS: ${{ matrix.cibw_arch }} CIBW_TEST_SKIP: "*_aarch64 *-macosx_arm64" CIBW_MANYLINUX_X86_64_IMAGE: ${{ matrix.cibw_manylinux }} CIBW_MANYLINUX_I686_IMAGE: ${{ matrix.cibw_manylinux }} + CC: clang - name: Rename Python version run: echo "PY_VERSION=$(echo ${{ matrix.cibw_python }} | cut -d- -f1)" >> $GITHUB_ENV - uses: actions/upload-artifact@v4 diff --git a/.github/workflows/test_template.yml b/.github/workflows/test_template.yml index 8dfd1eab8cf..5accae5e141 100644 --- a/.github/workflows/test_template.yml +++ b/.github/workflows/test_template.yml @@ -99,15 +99,12 @@ jobs: else tools/ci/install_dependencies.sh fi - # No need to update mingw-w64, we use msvc - # mingw-w64 does not manage well openmp threads so we avoid it for now - # Note that compilation works with mingw-w64 but 2-3 tests fail due to openmp - # - name: Install rtools (mingw-w64) - # if runner.os == 'Windows' - # run: | - # choco install rtools -y --no-progress --force --version=4.0.0.20220206 - # echo "/c/rtools40/ucrt64/bin;" >> $GITHUB_PATH - # echo "PKG_CONFIG_PATH=/c/opt/64/lib/pkgconfig;" >> $GITHUB_ENV + - name: Install OpenMP on macOS + if: runner.os == 'macOS' + run: | + brew install llvm + brew install libomp + echo "CC=clang" >> $GITHUB_ENV - name: Install DIPY run: | if [ "${{ inputs.use-pre }}" == "true" ]; then diff --git a/dipy/meson.build b/dipy/meson.build index 1bfc1b02ab3..3f14a196c8a 100644 --- a/dipy/meson.build +++ b/dipy/meson.build @@ -139,6 +139,22 @@ np_dep = declare_dependency(include_directories: inc_np) # Define Optimisation for cython extensions # ------------------------------------------------------------------------ omp = dependency('openmp', required: false) +if not omp.found() and meson.get_compiler('c').get_id() == 'clang' + # Check for libomp (OpenMP) using Homebrew + brew = find_program('brew', required : false) + if brew.found() + output = run_command(brew, 'list', 'libomp', check: true) + output = output.stdout().strip() + if output.contains('/libomp/') + omp_prefix = fs.parent(output.split('\n')[0]) + message('OpenMP Found: YES (Manual search) - ', omp_prefix) + omp = declare_dependency(compile_args : ['-Xpreprocessor', '-fopenmp'], + link_args : ['-L' + omp_prefix + '/lib', '-lomp'], + include_directories : include_directories(omp_prefix / 'include') + ) + endif + endif +endif # SSE intrinsics sse2_cflags = [] diff --git a/doc/devel/installation_from_source.rst b/doc/devel/installation_from_source.rst index ebb2e9c7d1f..deeb6ae89d4 100644 --- a/doc/devel/installation_from_source.rst +++ b/doc/devel/installation_from_source.rst @@ -198,7 +198,17 @@ If you are already using the Homebrew Python, or the standard python.org Python, you will need to use the llvm compiler with OMP. Run:: brew install llvm - export CC=/usr/local/bin/clang-omp + brew install libomp + export CC="/opt/homebrew/opt/llvm/bin/clang" + +In case the compiler or OpenMP are not detected , you can specify some environment +variables. For example, you can add the following lines to your ``~/.bash_profile``:: + + export CC="/opt/homebrew/opt/llvm/bin/clang" + export CXX="/opt/homebrew/opt/llvm/bin/clang++" + export LDFLAGS="-L/opt/homebrew/opt/libomp/lib -L/opt/homebrew/opt/llvm/lib" + export CPPFLAGS="-I/opt/homebrew/opt/libomp/include -I/opt/homebrew/opt/llvm/include" + export CFLAGS="-I/opt/homebrew/opt/libomp/include -I/opt/homebrew/opt/llvm/include" Building and installing ~~~~~~~~~~~~~~~~~~~~~~~ From f84003424b382999bc8822a3646ec6c8618ae58d Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Thu, 9 Nov 2023 13:29:28 -0500 Subject: [PATCH 05/19] WIP - adding fast tracking module --- dipy/tracking/fast_tracking.pxd | 14 ++++++++++++++ dipy/tracking/fast_tracking.pyx | 20 ++++++++++++++++++++ 2 files changed, 34 insertions(+) create mode 100644 dipy/tracking/fast_tracking.pxd create mode 100644 dipy/tracking/fast_tracking.pyx diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd new file mode 100644 index 00000000000..448e3a6d9d9 --- /dev/null +++ b/dipy/tracking/fast_tracking.pxd @@ -0,0 +1,14 @@ +cimport numpy as cnp + +from dipy.tracking.stopping_criterion cimport StoppingCriterion +from dipy.direction.pmf cimport PmfGen + +cpdef generate_tractogram(cnp.ndarray[cnp.float_t, ndim=2] seed_positons, + cnp.ndarray[cnp.float_t, ndim=2] seed_directions, + StoppingCriterion sc, + PmfGen pmf_gen) + +cdef double[:, :] generate_tractogram_c(double[:, :] seed_positons, + double[:, :] seed_directions, + StoppingCriterion sc, + PmfGen pmf_gen) \ No newline at end of file diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx new file mode 100644 index 00000000000..f560131f645 --- /dev/null +++ b/dipy/tracking/fast_tracking.pyx @@ -0,0 +1,20 @@ +cimport numpy as cnp + +from dipy.tracking.stopping_criterion cimport StoppingCriterion +from dipy.direction.pmf cimport PmfGen + +cpdef generate_tractogram(cnp.ndarray[cnp.float_t, ndim=2] seed_positons, + cnp.ndarray[cnp.float_t, ndim=2] seed_directions, + StoppingCriterion sc, + PmfGen pmf_gen): + cdef: + cnp.npy_intp _len=seed_positons.shape[0] + + return generate_tractogram_c( seed_positons, + seed_directions, + sc, pmf_gen) + +cdef double[:, :] generate_tractogram_c(double[:, :] seed_positons, + double[:, :] seed_directions, + StoppingCriterion sc, + PmfGen pmf_gen) \ No newline at end of file From ebbde750532dc439a4c5d679ced5bb23083c6bc4 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 17 Nov 2023 14:08:38 -0500 Subject: [PATCH 06/19] WIP --- dipy/tracking/fast_tracking.pxd | 82 ++++++++++++++++-- dipy/tracking/fast_tracking.pyx | 144 +++++++++++++++++++++++++++++--- 2 files changed, 205 insertions(+), 21 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 448e3a6d9d9..46e0f3a32e5 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -3,12 +3,76 @@ cimport numpy as cnp from dipy.tracking.stopping_criterion cimport StoppingCriterion from dipy.direction.pmf cimport PmfGen -cpdef generate_tractogram(cnp.ndarray[cnp.float_t, ndim=2] seed_positons, - cnp.ndarray[cnp.float_t, ndim=2] seed_directions, - StoppingCriterion sc, - PmfGen pmf_gen) - -cdef double[:, :] generate_tractogram_c(double[:, :] seed_positons, - double[:, :] seed_directions, - StoppingCriterion sc, - PmfGen pmf_gen) \ No newline at end of file + +cpdef list generate_tractogram(double[:,::1] seed_positons, + double[:,::1] seed_directions, + StoppingCriterion sc, + TrackingParameters params) + + +cdef int generate_tractogram_c(double[:, ::1] seed_positons, + double[:, ::1] seed_directions, + int nbr_seeds, + StoppingCriterion sc, + TrackingParameters params, + double[:,:,:] streamlines) nogil + +cdef double* get_pmf(double* point, + PmfGen pmf_gen, + double pmf_threshold, + int pmf_len) nogil + +cdef class TrackingParameters(): + pass + +cdef class ProbabilisticTrackingParameters(TrackingParameters): + cdef: + double cos_similarity + double pmf_threshold + PmfGen pmf_gen + int pmf_len + double step_size + double[:, :] vertices + + +cdef int probabilistic_tracker(double* point, + double* direction, + ProbabilisticTrackingParameters params) nogil + +cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): + pass + + +cdef int deterministic_maximum_tracker(double* point, + double* direction, + DeterministicTrackingParameters params) nogil + +cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): + cdef: + double angular_separation + double data_support_exponent + double[3][3] frame + double k1 + double k2 + double k_small + double last_val + double last_val_cand + double max_angle + double max_curvature + double[3] position + int probe_count + double probe_length + double probe_normalizer + int probe_quality + double probe_radius + double probe_step_size + double[9] propagator + int rejection_sampling_max_try + int rejection_sampling_nbr_sample + double[3] voxel_size + double[3] inv_voxel_size + + +cdef int paralle_transport_tracker(double* point, + double* direction, + ParalleTransportTrackingParameters params) nogil \ No newline at end of file diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index f560131f645..48dc5205d66 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -1,20 +1,140 @@ +import numpy as np cimport numpy as cnp -from dipy.tracking.stopping_criterion cimport StoppingCriterion from dipy.direction.pmf cimport PmfGen +from dipy.tracking.stopping_criterion cimport StoppingCriterion +from dipy.utils.fast_numpy cimport (copy_point, cumsum, norm, normalize, + where_to_insert, random) + -cpdef generate_tractogram(cnp.ndarray[cnp.float_t, ndim=2] seed_positons, - cnp.ndarray[cnp.float_t, ndim=2] seed_directions, - StoppingCriterion sc, - PmfGen pmf_gen): +cpdef list generate_tractogram(double[:,::1] seed_positons, + double[:,::1] seed_directions, + StoppingCriterion sc, + TrackingParameters params): cdef: cnp.npy_intp _len=seed_positons.shape[0] + double[:,:,:] streamlines + + streamlines = np.array(np.zeros((_len, 500, 3)), order='F') + + generate_tractogram_c(seed_positons, + seed_directions, + _len, sc, params, streamlines) + + return [] + + +cdef int generate_tractogram_c(double[:,::1] seed_positons, + double[:,::1] seed_directions, + int nbr_seeds, + StoppingCriterion sc, + TrackingParameters params, + double[:,:,:] streamlines) nogil: + + cdef: + cnp.npy_intp i + + + #for i in range(nbr_seeds): + + # for loop over all seed position and directions <> + + #initialize an empty streamline array + #do while stream_status is valid forward and backward + + #call the tracker to get the next direction + #probabilistic_tracker(point, dir, params) + + #stream_status = sc.check_point_c(new_pos) + + # copy the streamline to the results array + + return 1 + + +cdef double* get_pmf(double* point, + PmfGen pmf_gen, + double pmf_threshold, + int pmf_len) nogil: + cdef: + cnp.npy_intp i + double* pmf + double absolute_pmf_threshold + double max_pmf=0 + + pmf = pmf_gen.get_pmf_c(point) + for i in range(pmf_len): + if pmf[i] > max_pmf: + max_pmf = pmf[i] + absolute_pmf_threshold = pmf_threshold * max_pmf + + for i in range(pmf_len): + if pmf[i] < absolute_pmf_threshold: + pmf[i] = 0.0 + return pmf + + +cdef int probabilistic_tracker(double* point, + double* direction, + ProbabilisticTrackingParameters params) nogil: + cdef: + cnp.npy_intp i, idx, + double[:] newdir + double* pmf + double last_cdf, cos_sim + + pmf = get_pmf(point, params.pmf_gen, params.pmf_threshold, params.pmf_len) + + if norm(direction) == 0: + return 1 + normalize(direction) + + for i in range(params.pmf_len): + cos_sim = params.vertices[i][0] * direction[0] \ + + params.vertices[i][1] * direction[1] \ + + params.vertices[i][2] * direction[2] + if cos_sim < 0: + cos_sim = cos_sim * -1 + if cos_sim < params.cos_similarity: + pmf[i] = 0 + + cumsum(pmf, pmf, params.pmf_len) + last_cdf = pmf[params.pmf_len - 1] + if last_cdf == 0: + return 1 + + idx = where_to_insert(pmf, random() * last_cdf, params.pmf_len) + + newdir = params.vertices[idx] + # Update direction and return 0 for error + if (direction[0] * newdir[0] + + direction[1] * newdir[1] + + direction[2] * newdir[2] > 0): + copy_point(&newdir[0], direction) + else: + newdir[0] = newdir[0] * -1 + newdir[1] = newdir[1] * -1 + newdir[2] = newdir[2] * -1 + copy_point(&newdir[0], direction) + return 0 + +#get_direction_c of the DG +cdef int deterministic_maximum_tracker(double* point, + double* direction, + DeterministicTrackingParameters params) nogil: + # update point and dir with new position and direction + + # return 1 if the propagation failed. + + return 1 + +#get_direction_c of the DG +cdef int paralle_transport_tracker(double* point, + double* direction, + ParalleTransportTrackingParameters params) nogil: + # update point and dir with new position and direction + + # return 1 if the propagation failed. - return generate_tractogram_c( seed_positons, - seed_directions, - sc, pmf_gen) + return 1 -cdef double[:, :] generate_tractogram_c(double[:, :] seed_positons, - double[:, :] seed_directions, - StoppingCriterion sc, - PmfGen pmf_gen) \ No newline at end of file From 59f88eba68eaaee4bb4026f482f23a450cc39c11 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Wed, 22 Nov 2023 09:06:27 -0500 Subject: [PATCH 07/19] wip - cython generate tractogram/streamlines --- dipy/tracking/fast_tracking.pxd | 29 +++++--- dipy/tracking/fast_tracking.pyx | 115 +++++++++++++++++++++++++------- 2 files changed, 113 insertions(+), 31 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 46e0f3a32e5..92d8fe73990 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -6,16 +6,24 @@ from dipy.direction.pmf cimport PmfGen cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, - StoppingCriterion sc, TrackingParameters params) -cdef int generate_tractogram_c(double[:, ::1] seed_positons, - double[:, ::1] seed_directions, +cdef int generate_tractogram_c(double[:,::1] seed_positons, + double[:,::1] seed_directions, int nbr_seeds, - StoppingCriterion sc, TrackingParameters params, - double[:,:,:] streamlines) nogil + double[:,:,:] streamlines, + double[:] status) nogil + + +cdef int generate_local_streamline(double[::1] seed, + double[::1] position, + double[::1] stream_x, + double[::1] stream_y, + double[::1] stream_z, + TrackingParameters params) nogil + cdef double* get_pmf(double* point, PmfGen pmf_gen, @@ -23,6 +31,11 @@ cdef double* get_pmf(double* point, int pmf_len) nogil cdef class TrackingParameters(): + cdef: + StoppingCriterion sc + int max_len + double step_size + double[3] voxel_size pass cdef class ProbabilisticTrackingParameters(TrackingParameters): @@ -31,9 +44,9 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): double pmf_threshold PmfGen pmf_gen int pmf_len - double step_size - double[:, :] vertices + double[:, :] vertices + pass cdef int probabilistic_tracker(double* point, double* direction, @@ -69,8 +82,8 @@ cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): double[9] propagator int rejection_sampling_max_try int rejection_sampling_nbr_sample - double[3] voxel_size double[3] inv_voxel_size + pass cdef int paralle_transport_tracker(double* point, diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index 48dc5205d66..31b090b2271 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -6,50 +6,119 @@ from dipy.tracking.stopping_criterion cimport StoppingCriterion from dipy.utils.fast_numpy cimport (copy_point, cumsum, norm, normalize, where_to_insert, random) +from dipy.tracking.stopping_criterion cimport (StreamlineStatus, + StoppingCriterion, + TRACKPOINT, + ENDPOINT, + OUTSIDEIMAGE, + INVALIDPOINT) + +from libc.stdlib cimport malloc, free + cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, - StoppingCriterion sc, TrackingParameters params): cdef: cnp.npy_intp _len=seed_positons.shape[0] - double[:,:,:] streamlines + double[:,:,:] streamlines_arr + double[:] status_arr - streamlines = np.array(np.zeros((_len, 500, 3)), order='F') + # temporary array to store generated streamlines + streamlines_arr = np.array(np.zeros((_len, params.max_len, 3)), order='C') + status_arr = np.array(np.zeros((_len)) - 2, order='C') generate_tractogram_c(seed_positons, seed_directions, - _len, sc, params, streamlines) + _len, params, streamlines_arr, status_arr) - return [] + streamlines = [] + for i in range(_len): + # array to list + pass + + return streamlines cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,::1] seed_directions, int nbr_seeds, - StoppingCriterion sc, TrackingParameters params, - double[:,:,:] streamlines) nogil: - + double[:,:,:] streamlines, + double[:] status) nogil: cdef: - cnp.npy_intp i - - - #for i in range(nbr_seeds): + cnp.npy_intp i, j, k + double *stream + + # <> + for i in range(nbr_seeds): + stream_x = malloc(params.max_len * sizeof(double)) + stream_y = malloc(params.max_len * sizeof(double)) + stream_z = malloc(params.max_len * sizeof(double)) + + status[i] = generate_local_streamline(seed_positons[i], + seed_directions[i], + stream_x, + stream_y, + stream_z, + params) + for j in range(params.max_len): + streamlines[i,j,0] = stream_x[j] + streamlines[i,j,1] = stream_y[j] + streamlines[i,j,2] = stream_z[j] + free(stream_x) + free(stream_y) + free(stream_z) - # for loop over all seed position and directions <> - - #initialize an empty streamline array - #do while stream_status is valid forward and backward - - #call the tracker to get the next direction - #probabilistic_tracker(point, dir, params) + return 0 - #stream_status = sc.check_point_c(new_pos) - # copy the streamline to the results array +cdef int generate_local_streamline(double[::1] seed, + double[::1] direction, + double[::1] stream_x, + double[::1] stream_y, + double[::1] stream_z, + TrackingParameters params) nogil: + cdef: + cnp.npy_intp i, j + double point[3] + double voxdir[3] + StreamlineStatus stream_status + + # set the initial position + copy_point(&seed[0], point) + #copy_point(&seed[0], &streamline[0,0]) + stream_x[i] = seed[0] + stream_y[i] = seed[1] + stream_z[i] = seed[2] + + stream_status = TRACKPOINT + for i in range(1, params.max_len): + if probabilistic_tracker(point, direction, params): + break + + # update position + for j in range(3): + point[j] += direction[j] / params.voxel_size[j] * params.step_size + + #copy_point(point, &streamline[i, 0]) + stream_x[i] = point[0] + stream_y[i] = point[1] + stream_z[i] = point[2] + + stream_status = params.sc.check_point_c(point) + if stream_status == TRACKPOINT: + continue + elif (stream_status == ENDPOINT or + stream_status == INVALIDPOINT or + stream_status == OUTSIDEIMAGE): + break + else: + # maximum length of streamline has been reached, return everything + i = params.max_len - return 1 + # i should be return to know the length of the streamline + return stream_status cdef double* get_pmf(double* point, @@ -78,7 +147,7 @@ cdef int probabilistic_tracker(double* point, double* direction, ProbabilisticTrackingParameters params) nogil: cdef: - cnp.npy_intp i, idx, + cnp.npy_intp i, idx double[:] newdir double* pmf double last_cdf, cos_sim From 2e2b180a1c0340bd549200bb4a3442e2be34d01a Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Wed, 22 Nov 2023 09:27:44 -0500 Subject: [PATCH 08/19] wip - temporarly removed nogil --- dipy/tracking/fast_tracking.pxd | 18 +++++++++--------- dipy/tracking/fast_tracking.pyx | 23 +++++++++++------------ 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 92d8fe73990..9b4a0274809 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -14,21 +14,21 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, int nbr_seeds, TrackingParameters params, double[:,:,:] streamlines, - double[:] status) nogil + double[:] status) cdef int generate_local_streamline(double[::1] seed, double[::1] position, - double[::1] stream_x, - double[::1] stream_y, - double[::1] stream_z, - TrackingParameters params) nogil + double* stream_x, + double* stream_y, + double* stream_z, + TrackingParameters params) cdef double* get_pmf(double* point, PmfGen pmf_gen, double pmf_threshold, - int pmf_len) nogil + int pmf_len) cdef class TrackingParameters(): cdef: @@ -50,7 +50,7 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): cdef int probabilistic_tracker(double* point, double* direction, - ProbabilisticTrackingParameters params) nogil + ProbabilisticTrackingParameters params) cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): pass @@ -58,7 +58,7 @@ cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): cdef int deterministic_maximum_tracker(double* point, double* direction, - DeterministicTrackingParameters params) nogil + DeterministicTrackingParameters params) cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): cdef: @@ -88,4 +88,4 @@ cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): cdef int paralle_transport_tracker(double* point, double* direction, - ParalleTransportTrackingParameters params) nogil \ No newline at end of file + ParalleTransportTrackingParameters params) \ No newline at end of file diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index 31b090b2271..47dd1d5491d 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -45,10 +45,9 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, int nbr_seeds, TrackingParameters params, double[:,:,:] streamlines, - double[:] status) nogil: + double[:] status): cdef: - cnp.npy_intp i, j, k - double *stream + cnp.npy_intp i, j # <> for i in range(nbr_seeds): @@ -75,10 +74,10 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, cdef int generate_local_streamline(double[::1] seed, double[::1] direction, - double[::1] stream_x, - double[::1] stream_y, - double[::1] stream_z, - TrackingParameters params) nogil: + double* stream_x, + double* stream_y, + double* stream_z, + TrackingParameters params): cdef: cnp.npy_intp i, j double point[3] @@ -94,7 +93,7 @@ cdef int generate_local_streamline(double[::1] seed, stream_status = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(point, direction, params): + if probabilistic_tracker(&point[0], &direction[0], params): break # update position @@ -124,7 +123,7 @@ cdef int generate_local_streamline(double[::1] seed, cdef double* get_pmf(double* point, PmfGen pmf_gen, double pmf_threshold, - int pmf_len) nogil: + int pmf_len): cdef: cnp.npy_intp i double* pmf @@ -145,7 +144,7 @@ cdef double* get_pmf(double* point, cdef int probabilistic_tracker(double* point, double* direction, - ProbabilisticTrackingParameters params) nogil: + ProbabilisticTrackingParameters params): cdef: cnp.npy_intp i, idx double[:] newdir @@ -190,7 +189,7 @@ cdef int probabilistic_tracker(double* point, #get_direction_c of the DG cdef int deterministic_maximum_tracker(double* point, double* direction, - DeterministicTrackingParameters params) nogil: + DeterministicTrackingParameters params): # update point and dir with new position and direction # return 1 if the propagation failed. @@ -200,7 +199,7 @@ cdef int deterministic_maximum_tracker(double* point, #get_direction_c of the DG cdef int paralle_transport_tracker(double* point, double* direction, - ParalleTransportTrackingParameters params) nogil: + ParalleTransportTrackingParameters params): # update point and dir with new position and direction # return 1 if the propagation failed. From ecd2c6108b09ab92129727a635c24d14baa3db65 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 8 Dec 2023 16:41:36 -0500 Subject: [PATCH 09/19] WIP - working prototype --- dipy/tracking/fast_tracking.pxd | 11 +-- dipy/tracking/fast_tracking.pyx | 158 +++++++++++++++++++++----------- 2 files changed, 108 insertions(+), 61 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 9b4a0274809..138494dcf5c 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -17,11 +17,9 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:] status) -cdef int generate_local_streamline(double[::1] seed, - double[::1] position, - double* stream_x, - double* stream_y, - double* stream_z, +cdef int generate_local_streamline(double* seed, + double* position, + double* stream, TrackingParameters params) @@ -44,8 +42,9 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): double pmf_threshold PmfGen pmf_gen int pmf_len - double[:, :] vertices + + pass cdef int probabilistic_tracker(double* point, diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index 47dd1d5491d..c3ef77a8223 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -16,6 +16,7 @@ from dipy.tracking.stopping_criterion cimport (StreamlineStatus, from libc.stdlib cimport malloc, free +#need cpdef to access cdef functions? cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, TrackingParameters params): @@ -23,19 +24,29 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, cnp.npy_intp _len=seed_positons.shape[0] double[:,:,:] streamlines_arr double[:] status_arr + cnp.npy_intp i, j + + #print(_len, flush=True) # temporary array to store generated streamlines - streamlines_arr = np.array(np.zeros((_len, params.max_len, 3)), order='C') - status_arr = np.array(np.zeros((_len)) - 2, order='C') + streamlines_arr = np.array(np.zeros((_len, params.max_len * 2 + 1, 3)), order='C') + status_arr = np.array(np.zeros((_len)), order='C') + #stream = np.empty((self.max_length + 1, 3), dtype=float) generate_tractogram_c(seed_positons, seed_directions, _len, params, streamlines_arr, status_arr) - streamlines = [] for i in range(_len): # array to list - pass + s = [] + for j in range(params.max_len): + if norm(&streamlines_arr[i,j,0]) > 0: + s.append(list(np.copy(streamlines_arr[i,j]))) + else: + break + if len(s) > 1: + streamlines.append(np.array(s)) return streamlines @@ -47,77 +58,99 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,:,:] streamlines, double[:] status): cdef: - cnp.npy_intp i, j + cnp.npy_intp i, j, k + # <> for i in range(nbr_seeds): - stream_x = malloc(params.max_len * sizeof(double)) - stream_y = malloc(params.max_len * sizeof(double)) - stream_z = malloc(params.max_len * sizeof(double)) - - status[i] = generate_local_streamline(seed_positons[i], - seed_directions[i], - stream_x, - stream_y, - stream_z, + #stream_x = malloc((params.max_len * 2 + 1) * sizeof(double)) + #stream_y = malloc((params.max_len * 2 + 1) * sizeof(double)) + #stream_z = malloc((params.max_len * 2 + 1) * sizeof(double)) + stream = malloc((params.max_len * 3 * 2 + 1) * sizeof(double)) + + #for j in range(params.max_len * 2 + 1): + # stream_x[j] = 0 + # stream_y[j] = 0 + # stream_z[j] = 0 + + for j in range(params.max_len * 3 * 2 + 1): + stream[j] = 0 + + status[i] = generate_local_streamline(&seed_positons[i][0], + &seed_directions[i][0], + stream, params) - for j in range(params.max_len): - streamlines[i,j,0] = stream_x[j] - streamlines[i,j,1] = stream_y[j] - streamlines[i,j,2] = stream_z[j] - free(stream_x) - free(stream_y) - free(stream_z) + k = 0 + for j in range(params.max_len * 2 + 1): + if (stream[j * 3] != 0 + and stream[j * 3 + 1] !=0 + and stream[j * 3 + 2] != 0): + streamlines[i,k,0] = stream[j * 3] + streamlines[i,k,1] = stream[j * 3 + 1] + streamlines[i,k,2] = stream[j * 3 + 2] + k = k + 1 + #free(stream_x) + #free(stream_y) + #free(stream_z) + free(stream) return 0 -cdef int generate_local_streamline(double[::1] seed, - double[::1] direction, - double* stream_x, - double* stream_y, - double* stream_z, +cdef int generate_local_streamline(double* seed, + double* direction, + double* stream, TrackingParameters params): cdef: cnp.npy_intp i, j double point[3] double voxdir[3] - StreamlineStatus stream_status - + StreamlineStatus stream_status_forward, stream_status_backward # set the initial position - copy_point(&seed[0], point) - #copy_point(&seed[0], &streamline[0,0]) - stream_x[i] = seed[0] - stream_y[i] = seed[1] - stream_z[i] = seed[2] + copy_point(seed, point) + copy_point(direction, voxdir) + copy_point(seed, &stream[params.max_len * 3]) - stream_status = TRACKPOINT + # forward tracking + stream_status_forward = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(&point[0], &direction[0], params): + if probabilistic_tracker(&point[0], &voxdir[0], params): break - # update position for j in range(3): - point[j] += direction[j] / params.voxel_size[j] * params.step_size - - #copy_point(point, &streamline[i, 0]) - stream_x[i] = point[0] - stream_y[i] = point[1] - stream_z[i] = point[2] - - stream_status = params.sc.check_point_c(point) - if stream_status == TRACKPOINT: - continue - elif (stream_status == ENDPOINT or - stream_status == INVALIDPOINT or - stream_status == OUTSIDEIMAGE): + point[j] += voxdir[j] / params.voxel_size[j] * params.step_size + + copy_point(point, &stream[(params.max_len + i )* 3]) + + stream_status_forward = params.sc.check_point_c(point) + if (stream_status_forward == ENDPOINT or + stream_status_forward == INVALIDPOINT or + stream_status_forward == OUTSIDEIMAGE): break - else: - # maximum length of streamline has been reached, return everything - i = params.max_len - # i should be return to know the length of the streamline - return stream_status + # backward tracking + copy_point(seed, point) + copy_point(direction, voxdir) + for j in range(3): + voxdir[j] = voxdir[j] * -1 + stream_status_backward = TRACKPOINT + for i in range(1, params.max_len): + if probabilistic_tracker(&point[0], &voxdir[0], params): + break + # update position + for j in range(3): + point[j] += voxdir[j] / params.voxel_size[j] * params.step_size + + copy_point(point, &stream[(params.max_len + i )* 3]) + + + stream_status_backward = params.sc.check_point_c(point) + if (stream_status_backward == ENDPOINT or + stream_status_backward == INVALIDPOINT or + stream_status_backward == OUTSIDEIMAGE): + break + # need to handle stream status + return 0 #stream_status cdef double* get_pmf(double* point, @@ -174,7 +207,7 @@ cdef int probabilistic_tracker(double* point, idx = where_to_insert(pmf, random() * last_cdf, params.pmf_len) newdir = params.vertices[idx] - # Update direction and return 0 for error + # Update direction if (direction[0] * newdir[0] + direction[1] * newdir[1] + direction[2] * newdir[2] > 0): @@ -206,3 +239,18 @@ cdef int paralle_transport_tracker(double* point, return 1 + + +cdef class ProbabilisticTrackingParameters(TrackingParameters): + + def __cinit__(self, sc, max_len, step_size, voxel_size, cos_similarity, + pmf_threshold, pmf_gen, pmf_len, vertices): + self.sc = sc + self.max_len = max_len + self.step_size = step_size + self.voxel_size = voxel_size + self.cos_similarity = cos_similarity + self.pmf_threshold = pmf_threshold + self.pmf_gen = pmf_gen + self.pmf_len = pmf_len + self.vertices = vertices From 1aa30cce7361f00ce4b76e0f043898504b9bc0ff Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Sat, 9 Dec 2023 08:39:47 -0500 Subject: [PATCH 10/19] OPT - cython code --- dipy/tracking/fast_tracking.pyx | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index c3ef77a8223..d3708dd67ed 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -32,7 +32,6 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, streamlines_arr = np.array(np.zeros((_len, params.max_len * 2 + 1, 3)), order='C') status_arr = np.array(np.zeros((_len)), order='C') - #stream = np.empty((self.max_length + 1, 3), dtype=float) generate_tractogram_c(seed_positons, seed_directions, _len, params, streamlines_arr, status_arr) @@ -63,16 +62,10 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, # <> for i in range(nbr_seeds): - #stream_x = malloc((params.max_len * 2 + 1) * sizeof(double)) - #stream_y = malloc((params.max_len * 2 + 1) * sizeof(double)) - #stream_z = malloc((params.max_len * 2 + 1) * sizeof(double)) stream = malloc((params.max_len * 3 * 2 + 1) * sizeof(double)) - #for j in range(params.max_len * 2 + 1): - # stream_x[j] = 0 - # stream_y[j] = 0 - # stream_z[j] = 0 - + # initialize to 0. It will be replaced when better handling various + # streamline lengtyh for j in range(params.max_len * 3 * 2 + 1): stream[j] = 0 @@ -80,6 +73,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, &seed_directions[i][0], stream, params) + # copy the v k = 0 for j in range(params.max_len * 2 + 1): if (stream[j * 3] != 0 @@ -89,9 +83,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, streamlines[i,k,1] = stream[j * 3 + 1] streamlines[i,k,2] = stream[j * 3 + 2] k = k + 1 - #free(stream_x) - #free(stream_y) - #free(stream_z) + free(stream) return 0 From ddaa2bc03d18c8c1806e1fc84b49a44adfeabc32 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Mon, 4 Dec 2023 23:18:36 -0500 Subject: [PATCH 11/19] functor example --- dipy/tracking/fast_tracking.pxd | 21 ++++++--- dipy/tracking/fast_tracking.pyx | 4 +- dipy/tracking/tracker.py | 81 +++++++++++++++++++++++++++++++++ 3 files changed, 98 insertions(+), 8 deletions(-) create mode 100644 dipy/tracking/tracker.py diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 138494dcf5c..abddda95181 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -4,6 +4,17 @@ from dipy.tracking.stopping_criterion cimport StoppingCriterion from dipy.direction.pmf cimport PmfGen +cdef class TrackingParameters(): + cdef: + StoppingCriterion sc + int max_len + double step_size + double[3] voxel_size + pass + + +ctypedef int (*func_ptr)(double* point, double* direction, ProbabilisticTrackingParameters) + cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, TrackingParameters params) @@ -20,6 +31,9 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, cdef int generate_local_streamline(double* seed, double* position, double* stream, + func_ptr tracker, + # sc_ptr stopping_criterion, + # pmf_ptr pmf_gen, TrackingParameters params) @@ -28,13 +42,6 @@ cdef double* get_pmf(double* point, double pmf_threshold, int pmf_len) -cdef class TrackingParameters(): - cdef: - StoppingCriterion sc - int max_len - double step_size - double[3] voxel_size - pass cdef class ProbabilisticTrackingParameters(TrackingParameters): cdef: diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index d3708dd67ed..e10d2595564 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -72,6 +72,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, status[i] = generate_local_streamline(&seed_positons[i][0], &seed_directions[i][0], stream, + &probabilistic_tracker, params) # copy the v k = 0 @@ -92,6 +93,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, cdef int generate_local_streamline(double* seed, double* direction, double* stream, + func_ptr tracker, TrackingParameters params): cdef: cnp.npy_intp i, j @@ -106,7 +108,7 @@ cdef int generate_local_streamline(double* seed, # forward tracking stream_status_forward = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(&point[0], &voxdir[0], params): + if tracker(&point[0], &voxdir[0], params): # , pmf_gen_func): break # update position for j in range(3): diff --git a/dipy/tracking/tracker.py b/dipy/tracking/tracker.py new file mode 100644 index 00000000000..44815b640d1 --- /dev/null +++ b/dipy/tracking/tracker.py @@ -0,0 +1,81 @@ +from dipy.tracking.fast_tracking import generate_tractogram + + + +def probabilistic_directions(): + pass + + +def deterministic_directions(): + pass + + +def eudx_directions(): + pass + +def build_stopping_criterion(mask=None, check_point=None, include=None, exclude=None): + sc_element = { + "mask": mask, + "wm_map": None, + "gm_map": None, + "csf_map": None, + "checkpoint": check_point, + "include": include, + "exclude": exclude + } + return sc_element + + + +def generate_streamlines(): + pass + + +def generate_tractogram_py(seeds, streamline_generator=generate_streamlines, postprocess_streamline=None): + for s in seeds: + streamline = streamline_generator(s, stopping_criterion, direction_getter, max_cross, step_size, maxlen, fixedstep) + if postprocess_streamline is not None: + streamline = postprocess_streamline(streamlines) + yield streamline + + +def build_tractogram_algorithm(streamline_generator, tracker_func, pmf_gen_func, + stopping_criterion_func, postprocess_streamline=None): + + def generate_tractogram(seed_positions, seed_directions, tracking_parameters): + for s in len(seed_positions): + streamline = streamline_generator(seed_position[i], seed_directions[i], + tracker_func, pmf_gen_func, + stopping_criterion_func) + if postprocess_streamline is not None: + streamline = postprocess_streamline(streamlines) + yield streamline + + return generate_tractogram + + +def local_tracking(seed_positons, seed_directions, *, tracking_parameters=None, tractogram_func=None, postprocess_streamline=None): + tractogram_func = tractogram_func or generate_tractogram + if tractogram_func and 'cython' in type(tractogram_func): + return tractogram_func(seed_positons, seed_directions, tracking_parameters) + + return generate_tractogram_py(seed_positons, seed_directions, tracking_parameters) + + +probabilistic_tracking = partial(local_tracking, tractogram_func=generate_tractogram) +probabilistic_tracking_python = partial(local_tracking, tractogram_func=generate_tractogram_py) + +deterministic_tracking = partial(local_tracking, direction_getter='deterministic') +eudx_tracking = partial(local_tracking, direction_getter='eudx') +fact_tracking = partial(local_tracking, direction_getter='fact') +ptt_tracking = partial(local_tracking, direction_getter='ptt') +pft_tracking = partial(local_tracking, direction_getter='pft') + + +def test_probabilistic_tracking(): + pass + + +if __name__ == "__main__": + my_generate_tractogram = build_tractogram_algorithm() + streamlines = local_tracking(p_seed, d_seed, params, tractogram_func=my_generate_tractogram) From d9ee878ff7880455149c55640f32ec4e1ca6d652 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Thu, 11 Jan 2024 10:55:53 -0500 Subject: [PATCH 12/19] add fast tracking to meson --- dipy/tracking/meson.build | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dipy/tracking/meson.build b/dipy/tracking/meson.build index 14aa621796a..1ab20a02f09 100644 --- a/dipy/tracking/meson.build +++ b/dipy/tracking/meson.build @@ -1,6 +1,7 @@ cython_sources = [ 'direction_getter', 'distances', + 'fast_tracking', 'fbcmeasures', 'localtrack', 'propspeed', @@ -11,6 +12,7 @@ cython_sources = [ cython_headers = [ 'direction_getter.pxd', + 'fast_tracking.pxd', 'fbcmeasures.pxd', 'propspeed.pxd', 'stopping_criterion.pxd', From 9b719c2710e87507d928e066f92106ed1ccc1319 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Thu, 11 Jan 2024 14:14:29 -0500 Subject: [PATCH 13/19] RF - move stopping criteria out of tracking params --- dipy/tracking/fast_tracking.pxd | 8 ++++++-- dipy/tracking/fast_tracking.pyx | 18 +++++++++--------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index abddda95181..595d5463b8d 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -6,23 +6,26 @@ from dipy.direction.pmf cimport PmfGen cdef class TrackingParameters(): cdef: - StoppingCriterion sc int max_len double step_size double[3] voxel_size pass -ctypedef int (*func_ptr)(double* point, double* direction, ProbabilisticTrackingParameters) +ctypedef int (*func_ptr)(double* point, + double* direction, + ProbabilisticTrackingParameters) cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, + StoppingCriterion sc, TrackingParameters params) cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,::1] seed_directions, int nbr_seeds, + StoppingCriterion sc, TrackingParameters params, double[:,:,:] streamlines, double[:] status) @@ -34,6 +37,7 @@ cdef int generate_local_streamline(double* seed, func_ptr tracker, # sc_ptr stopping_criterion, # pmf_ptr pmf_gen, + StoppingCriterion sc, TrackingParameters params) diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index e10d2595564..faefaef25a6 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -19,6 +19,7 @@ from libc.stdlib cimport malloc, free #need cpdef to access cdef functions? cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, + StoppingCriterion sc, TrackingParameters params): cdef: cnp.npy_intp _len=seed_positons.shape[0] @@ -26,15 +27,12 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, double[:] status_arr cnp.npy_intp i, j - #print(_len, flush=True) - # temporary array to store generated streamlines streamlines_arr = np.array(np.zeros((_len, params.max_len * 2 + 1, 3)), order='C') status_arr = np.array(np.zeros((_len)), order='C') - generate_tractogram_c(seed_positons, - seed_directions, - _len, params, streamlines_arr, status_arr) + generate_tractogram_c(seed_positons, seed_directions, _len, sc, params, + streamlines_arr, status_arr) streamlines = [] for i in range(_len): # array to list @@ -53,6 +51,7 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,::1] seed_directions, int nbr_seeds, + StoppingCriterion sc, TrackingParameters params, double[:,:,:] streamlines, double[:] status): @@ -73,6 +72,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, &seed_directions[i][0], stream, &probabilistic_tracker, + sc, params) # copy the v k = 0 @@ -94,6 +94,7 @@ cdef int generate_local_streamline(double* seed, double* direction, double* stream, func_ptr tracker, + StoppingCriterion sc, TrackingParameters params): cdef: cnp.npy_intp i, j @@ -116,7 +117,7 @@ cdef int generate_local_streamline(double* seed, copy_point(point, &stream[(params.max_len + i )* 3]) - stream_status_forward = params.sc.check_point_c(point) + stream_status_forward = sc.check_point_c(point) if (stream_status_forward == ENDPOINT or stream_status_forward == INVALIDPOINT or stream_status_forward == OUTSIDEIMAGE): @@ -138,7 +139,7 @@ cdef int generate_local_streamline(double* seed, copy_point(point, &stream[(params.max_len + i )* 3]) - stream_status_backward = params.sc.check_point_c(point) + stream_status_backward = sc.check_point_c(point) if (stream_status_backward == ENDPOINT or stream_status_backward == INVALIDPOINT or stream_status_backward == OUTSIDEIMAGE): @@ -237,9 +238,8 @@ cdef int paralle_transport_tracker(double* point, cdef class ProbabilisticTrackingParameters(TrackingParameters): - def __cinit__(self, sc, max_len, step_size, voxel_size, cos_similarity, + def __cinit__(self, max_len, step_size, voxel_size, cos_similarity, pmf_threshold, pmf_gen, pmf_len, vertices): - self.sc = sc self.max_len = max_len self.step_size = step_size self.voxel_size = voxel_size From 2c4c6ffa8c221e6b637b5bf238dfbddead431cd8 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 12 Jan 2024 11:56:04 -0500 Subject: [PATCH 14/19] OPT - removed python objects --- dipy/tracking/fast_tracking.pxd | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 595d5463b8d..063c258a3dc 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -9,7 +9,7 @@ cdef class TrackingParameters(): int max_len double step_size double[3] voxel_size - pass + double[3] inv_voxel_size ctypedef int (*func_ptr)(double* point, @@ -27,6 +27,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, int nbr_seeds, StoppingCriterion sc, TrackingParameters params, + func_ptr traker, double[:,:,:] streamlines, double[:] status) @@ -56,8 +57,6 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): double[:, :] vertices - pass - cdef int probabilistic_tracker(double* point, double* direction, ProbabilisticTrackingParameters params) @@ -92,8 +91,6 @@ cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): double[9] propagator int rejection_sampling_max_try int rejection_sampling_nbr_sample - double[3] inv_voxel_size - pass cdef int paralle_transport_tracker(double* point, From ea8242f618d7c9e56eadcfe0ca06418eef0bf9c0 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 12 Jan 2024 12:27:16 -0500 Subject: [PATCH 15/19] OPT - move pmf_gen out of params --- dipy/tracking/fast_tracking.pxd | 23 +++++++---- dipy/tracking/fast_tracking.pyx | 67 +++++++++++++++++++++++++-------- 2 files changed, 66 insertions(+), 24 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 063c258a3dc..ec3ee2c73f3 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -8,18 +8,20 @@ cdef class TrackingParameters(): cdef: int max_len double step_size - double[3] voxel_size + double[:] voxel_size double[3] inv_voxel_size ctypedef int (*func_ptr)(double* point, double* direction, - ProbabilisticTrackingParameters) + ProbabilisticTrackingParameters, + PmfGen) cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, StoppingCriterion sc, - TrackingParameters params) + TrackingParameters params, + PmfGen pmf_gen) cdef int generate_tractogram_c(double[:,::1] seed_positons, @@ -27,6 +29,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, int nbr_seeds, StoppingCriterion sc, TrackingParameters params, + PmfGen pmf_gen, func_ptr traker, double[:,:,:] streamlines, double[:] status) @@ -39,7 +42,8 @@ cdef int generate_local_streamline(double* seed, # sc_ptr stopping_criterion, # pmf_ptr pmf_gen, StoppingCriterion sc, - TrackingParameters params) + TrackingParameters params, + PmfGen pmf_gen) cdef double* get_pmf(double* point, @@ -52,14 +56,15 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): cdef: double cos_similarity double pmf_threshold - PmfGen pmf_gen + #PmfGen pmf_gen int pmf_len double[:, :] vertices cdef int probabilistic_tracker(double* point, double* direction, - ProbabilisticTrackingParameters params) + ProbabilisticTrackingParameters params, + PmfGen pmf_gen) cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): pass @@ -67,7 +72,8 @@ cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): cdef int deterministic_maximum_tracker(double* point, double* direction, - DeterministicTrackingParameters params) + DeterministicTrackingParameters params, + PmfGen pmf_gen,) cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): cdef: @@ -95,4 +101,5 @@ cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): cdef int paralle_transport_tracker(double* point, double* direction, - ParalleTransportTrackingParameters params) \ No newline at end of file + ParalleTransportTrackingParameters params, + PmfGen pmf_gen) \ No newline at end of file diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index faefaef25a6..a38c7fccfdc 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -1,3 +1,7 @@ +# cython: boundscheck=False +# cython: initializedcheck=False +# cython: wraparound=False + import numpy as np cimport numpy as cnp @@ -20,7 +24,8 @@ from libc.stdlib cimport malloc, free cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, StoppingCriterion sc, - TrackingParameters params): + TrackingParameters params, + PmfGen pmf_gen): cdef: cnp.npy_intp _len=seed_positons.shape[0] double[:,:,:] streamlines_arr @@ -32,7 +37,7 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, status_arr = np.array(np.zeros((_len)), order='C') generate_tractogram_c(seed_positons, seed_directions, _len, sc, params, - streamlines_arr, status_arr) + pmf_gen, &probabilistic_tracker, streamlines_arr, status_arr) streamlines = [] for i in range(_len): # array to list @@ -53,6 +58,8 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, int nbr_seeds, StoppingCriterion sc, TrackingParameters params, + PmfGen pmf_gen, + func_ptr tracker, double[:,:,:] streamlines, double[:] status): cdef: @@ -71,9 +78,10 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, status[i] = generate_local_streamline(&seed_positons[i][0], &seed_directions[i][0], stream, - &probabilistic_tracker, + tracker, sc, - params) + params, + pmf_gen) # copy the v k = 0 for j in range(params.max_len * 2 + 1): @@ -95,7 +103,8 @@ cdef int generate_local_streamline(double* seed, double* stream, func_ptr tracker, StoppingCriterion sc, - TrackingParameters params): + TrackingParameters params, + PmfGen pmf_gen): cdef: cnp.npy_intp i, j double point[3] @@ -109,11 +118,11 @@ cdef int generate_local_streamline(double* seed, # forward tracking stream_status_forward = TRACKPOINT for i in range(1, params.max_len): - if tracker(&point[0], &voxdir[0], params): # , pmf_gen_func): + if tracker(&point[0], &voxdir[0], params, pmf_gen): # , pmf_gen_func): break # update position for j in range(3): - point[j] += voxdir[j] / params.voxel_size[j] * params.step_size + point[j] += voxdir[j] * params.inv_voxel_size[j] * params.step_size copy_point(point, &stream[(params.max_len + i )* 3]) @@ -130,11 +139,11 @@ cdef int generate_local_streamline(double* seed, voxdir[j] = voxdir[j] * -1 stream_status_backward = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(&point[0], &voxdir[0], params): + if tracker(&point[0], &voxdir[0], params, pmf_gen): break # update position for j in range(3): - point[j] += voxdir[j] / params.voxel_size[j] * params.step_size + point[j] += voxdir[j] * params.inv_voxel_size[j] * params.step_size copy_point(point, &stream[(params.max_len + i )* 3]) @@ -172,14 +181,32 @@ cdef double* get_pmf(double* point, cdef int probabilistic_tracker(double* point, double* direction, - ProbabilisticTrackingParameters params): + ProbabilisticTrackingParameters params, + PmfGen pmf_gen): cdef: cnp.npy_intp i, idx double[:] newdir double* pmf double last_cdf, cos_sim + double max_pmf=0 + double absolute_pmf_threshold + + #pmf = get_pmf(point, params.pmf_gen, params.pmf_threshold, params.pmf_len) + + + pmf = pmf_gen.get_pmf_c(point) + for i in range(params.pmf_len): + if pmf[i] > max_pmf: + max_pmf = pmf[i] + absolute_pmf_threshold = params.pmf_threshold * max_pmf + + for i in range(params.pmf_len): + if pmf[i] < absolute_pmf_threshold: + pmf[i] = 0.0 + + + - pmf = get_pmf(point, params.pmf_gen, params.pmf_threshold, params.pmf_len) if norm(direction) == 0: return 1 @@ -217,7 +244,8 @@ cdef int probabilistic_tracker(double* point, #get_direction_c of the DG cdef int deterministic_maximum_tracker(double* point, double* direction, - DeterministicTrackingParameters params): + DeterministicTrackingParameters params, + PmfGen pmf_gen): # update point and dir with new position and direction # return 1 if the propagation failed. @@ -227,7 +255,8 @@ cdef int deterministic_maximum_tracker(double* point, #get_direction_c of the DG cdef int paralle_transport_tracker(double* point, double* direction, - ParalleTransportTrackingParameters params): + ParalleTransportTrackingParameters params, + PmfGen pmf_gen): # update point and dir with new position and direction # return 1 if the propagation failed. @@ -238,13 +267,19 @@ cdef int paralle_transport_tracker(double* point, cdef class ProbabilisticTrackingParameters(TrackingParameters): - def __cinit__(self, max_len, step_size, voxel_size, cos_similarity, - pmf_threshold, pmf_gen, pmf_len, vertices): + def __cinit__(self, int max_len, double step_size, double[:] voxel_size, + double cos_similarity, double pmf_threshold, #PmfGen pmf_gen, + int pmf_len, double[:,:] vertices): + cdef: + cnp.npy_intp i self.max_len = max_len self.step_size = step_size self.voxel_size = voxel_size self.cos_similarity = cos_similarity self.pmf_threshold = pmf_threshold - self.pmf_gen = pmf_gen + #self.pmf_gen = pmf_gen self.pmf_len = pmf_len self.vertices = vertices + + for i in range(3): + self.inv_voxel_size[i] = 1. / voxel_size[i] From 7bf71c399166af27f28e77615193a7f675b9d8fd Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 12 Jan 2024 13:17:09 -0500 Subject: [PATCH 16/19] RF - removed unused PmfGen class member --- dipy/direction/pmf.pxd | 1 - dipy/direction/pmf.pyx | 1 - 2 files changed, 2 deletions(-) diff --git a/dipy/direction/pmf.pxd b/dipy/direction/pmf.pxd index e49364f37d0..0f3a06aa654 100644 --- a/dipy/direction/pmf.pxd +++ b/dipy/direction/pmf.pxd @@ -5,7 +5,6 @@ cdef class PmfGen: double[:] pmf double[:, :, :, :] data double[:, :] vertices - object sphere cpdef double[:] get_pmf(self, double[::1] point) cdef double* get_pmf_c(self, double* point) noexcept nogil diff --git a/dipy/direction/pmf.pyx b/dipy/direction/pmf.pyx index cbee0fc94c2..1a7ecdfc3d0 100644 --- a/dipy/direction/pmf.pyx +++ b/dipy/direction/pmf.pyx @@ -16,7 +16,6 @@ cdef class PmfGen: double[:, :, :, :] data, object sphere): self.data = np.asarray(data, dtype=float, order='C') - self.sphere = sphere self.vertices = np.asarray(sphere.vertices, dtype=float) cpdef double[:] get_pmf(self, double[::1] point): From 146c3d8bf7880a4cee8991c5f75cbf90cd7e02eb Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Thu, 18 Jan 2024 15:28:19 -0500 Subject: [PATCH 17/19] [ENH] allow noexcept nogil --- dipy/tracking/fast_tracking.pxd | 14 ++++----- dipy/tracking/fast_tracking.pyx | 43 ++++++++++++---------------- dipy/tracking/stopping_criterion.pxd | 2 +- dipy/tracking/stopping_criterion.pyx | 21 ++++++++------ dipy/utils/fast_numpy.pxd | 2 +- 5 files changed, 40 insertions(+), 42 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index ec3ee2c73f3..7f07ddeaaf1 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -15,12 +15,12 @@ cdef class TrackingParameters(): ctypedef int (*func_ptr)(double* point, double* direction, ProbabilisticTrackingParameters, - PmfGen) + PmfGen) nogil cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, StoppingCriterion sc, - TrackingParameters params, + ProbabilisticTrackingParameters params, PmfGen pmf_gen) @@ -28,7 +28,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,::1] seed_directions, int nbr_seeds, StoppingCriterion sc, - TrackingParameters params, + ProbabilisticTrackingParameters params, PmfGen pmf_gen, func_ptr traker, double[:,:,:] streamlines, @@ -42,14 +42,14 @@ cdef int generate_local_streamline(double* seed, # sc_ptr stopping_criterion, # pmf_ptr pmf_gen, StoppingCriterion sc, - TrackingParameters params, - PmfGen pmf_gen) + ProbabilisticTrackingParameters params, + PmfGen pmf_gen) noexcept nogil cdef double* get_pmf(double* point, PmfGen pmf_gen, double pmf_threshold, - int pmf_len) + int pmf_len) noexcept nogil cdef class ProbabilisticTrackingParameters(TrackingParameters): @@ -64,7 +64,7 @@ cdef class ProbabilisticTrackingParameters(TrackingParameters): cdef int probabilistic_tracker(double* point, double* direction, ProbabilisticTrackingParameters params, - PmfGen pmf_gen) + PmfGen pmf_gen) noexcept nogil cdef class DeterministicTrackingParameters(ProbabilisticTrackingParameters): pass diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index a38c7fccfdc..ca2b0872e67 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -1,7 +1,9 @@ # cython: boundscheck=False # cython: initializedcheck=False # cython: wraparound=False +# cython: Nonecheck=False +from cython.parallel import prange import numpy as np cimport numpy as cnp @@ -24,7 +26,7 @@ from libc.stdlib cimport malloc, free cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,::1] seed_directions, StoppingCriterion sc, - TrackingParameters params, + ProbabilisticTrackingParameters params, PmfGen pmf_gen): cdef: cnp.npy_intp _len=seed_positons.shape[0] @@ -57,7 +59,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, double[:,::1] seed_directions, int nbr_seeds, StoppingCriterion sc, - TrackingParameters params, + ProbabilisticTrackingParameters params, PmfGen pmf_gen, func_ptr tracker, double[:,:,:] streamlines, @@ -66,8 +68,8 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, cnp.npy_intp i, j, k - # <> - for i in range(nbr_seeds): + for i in prange(nbr_seeds, nogil=True): + # for i in range(nbr_seeds): stream = malloc((params.max_len * 3 * 2 + 1) * sizeof(double)) # initialize to 0. It will be replaced when better handling various @@ -103,13 +105,14 @@ cdef int generate_local_streamline(double* seed, double* stream, func_ptr tracker, StoppingCriterion sc, - TrackingParameters params, - PmfGen pmf_gen): + ProbabilisticTrackingParameters params, + PmfGen pmf_gen) noexcept nogil: cdef: cnp.npy_intp i, j double point[3] double voxdir[3] StreamlineStatus stream_status_forward, stream_status_backward + # set the initial position copy_point(seed, point) copy_point(direction, voxdir) @@ -118,7 +121,7 @@ cdef int generate_local_streamline(double* seed, # forward tracking stream_status_forward = TRACKPOINT for i in range(1, params.max_len): - if tracker(&point[0], &voxdir[0], params, pmf_gen): # , pmf_gen_func): + if probabilistic_tracker(&point[0], &voxdir[0], params, pmf_gen): # , pmf_gen_func): break # update position for j in range(3): @@ -139,7 +142,7 @@ cdef int generate_local_streamline(double* seed, voxdir[j] = voxdir[j] * -1 stream_status_backward = TRACKPOINT for i in range(1, params.max_len): - if tracker(&point[0], &voxdir[0], params, pmf_gen): + if probabilistic_tracker(&point[0], &voxdir[0], params, pmf_gen): break # update position for j in range(3): @@ -160,7 +163,7 @@ cdef int generate_local_streamline(double* seed, cdef double* get_pmf(double* point, PmfGen pmf_gen, double pmf_threshold, - int pmf_len): + int pmf_len) noexcept nogil: cdef: cnp.npy_intp i double* pmf @@ -182,27 +185,17 @@ cdef double* get_pmf(double* point, cdef int probabilistic_tracker(double* point, double* direction, ProbabilisticTrackingParameters params, - PmfGen pmf_gen): + PmfGen pmf_gen) noexcept nogil: cdef: cnp.npy_intp i, idx - double[:] newdir + double* newdir double* pmf double last_cdf, cos_sim double max_pmf=0 double absolute_pmf_threshold - #pmf = get_pmf(point, params.pmf_gen, params.pmf_threshold, params.pmf_len) - - - pmf = pmf_gen.get_pmf_c(point) - for i in range(params.pmf_len): - if pmf[i] > max_pmf: - max_pmf = pmf[i] - absolute_pmf_threshold = params.pmf_threshold * max_pmf + pmf = get_pmf(point, pmf_gen, params.pmf_threshold, params.pmf_len) - for i in range(params.pmf_len): - if pmf[i] < absolute_pmf_threshold: - pmf[i] = 0.0 @@ -228,17 +221,17 @@ cdef int probabilistic_tracker(double* point, idx = where_to_insert(pmf, random() * last_cdf, params.pmf_len) - newdir = params.vertices[idx] + newdir = ¶ms.vertices[idx][0] # Update direction if (direction[0] * newdir[0] + direction[1] * newdir[1] + direction[2] * newdir[2] > 0): - copy_point(&newdir[0], direction) + copy_point(newdir, direction) else: newdir[0] = newdir[0] * -1 newdir[1] = newdir[1] * -1 newdir[2] = newdir[2] * -1 - copy_point(&newdir[0], direction) + copy_point(newdir, direction) return 0 #get_direction_c of the DG diff --git a/dipy/tracking/stopping_criterion.pxd b/dipy/tracking/stopping_criterion.pxd index 044c8fe46a6..74b9a023d62 100644 --- a/dipy/tracking/stopping_criterion.pxd +++ b/dipy/tracking/stopping_criterion.pxd @@ -12,7 +12,7 @@ cdef class StoppingCriterion: double interp_out_double[1] double[::1] interp_out_view cpdef StreamlineStatus check_point(self, double[::1] point) - cdef StreamlineStatus check_point_c(self, double* point) + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil cdef class BinaryStoppingCriterion(StoppingCriterion): diff --git a/dipy/tracking/stopping_criterion.pyx b/dipy/tracking/stopping_criterion.pyx index f08114ba002..3cd94e4278d 100644 --- a/dipy/tracking/stopping_criterion.pyx +++ b/dipy/tracking/stopping_criterion.pyx @@ -6,6 +6,7 @@ cdef extern from "dpy_math.h" nogil: int dpy_rint(double) +from dipy.utils.fast_numpy cimport random from dipy.core.interpolation cimport trilinear_interpolate4d_c import numpy as np @@ -17,7 +18,7 @@ cdef class StoppingCriterion: return self.check_point_c(&point[0]) - cdef StreamlineStatus check_point_c(self, double* point): + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil: pass @@ -31,7 +32,7 @@ cdef class BinaryStoppingCriterion(StoppingCriterion): self.interp_out_view = self.interp_out_double self.mask = (mask > 0).astype('uint8') - cdef StreamlineStatus check_point_c(self, double* point): + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil: cdef: unsigned char result int err @@ -68,7 +69,7 @@ cdef class ThresholdStoppingCriterion(StoppingCriterion): self.metric_map = np.asarray(metric_map, 'float64') self.threshold = threshold - cdef StreamlineStatus check_point_c(self, double* point): + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil: cdef: double result int err @@ -189,7 +190,7 @@ cdef class ActStoppingCriterion(AnatomicalStoppingCriterion): self.include_map = np.asarray(include_map, 'float64') self.exclude_map = np.asarray(exclude_map, 'float64') - cdef StreamlineStatus check_point_c(self, double* point): + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil: cdef: double include_result, exclude_result int include_err, exclude_err @@ -251,10 +252,11 @@ cdef class CmcStoppingCriterion(AnatomicalStoppingCriterion): self.average_voxel_size = average_voxel_size self.correction_factor = step_size / average_voxel_size - cdef StreamlineStatus check_point_c(self, double* point): + cdef StreamlineStatus check_point_c(self, double* point) noexcept nogil: cdef: double include_result, exclude_result int include_err, exclude_err + double p include_err = trilinear_interpolate4d_c(self.include_map[..., None], point, self.interp_out_view) @@ -277,19 +279,22 @@ cdef class CmcStoppingCriterion(AnatomicalStoppingCriterion): raise RuntimeError("Unexpected interpolation error " + "(exclude_map - code:%i)" % exclude_err) - rng = np.random.default_rng() + # rng = np.random.default_rng() + # test if the tracking continues if include_result + exclude_result <= 0: return TRACKPOINT num = max(0, (1 - include_result - exclude_result)) den = num + include_result + exclude_result p = (num / den) ** self.correction_factor - if rng.random() < p: + # if rng.random() < p: + if random() < p: return TRACKPOINT # test if the tracking stopped in the include tissue map p = (include_result / (include_result + exclude_result)) - if rng.random() < p: + # if rng.random() < p: + if random() < p: return ENDPOINT # the tracking stopped in the exclude tissue map diff --git a/dipy/utils/fast_numpy.pxd b/dipy/utils/fast_numpy.pxd index 26f58793f2a..0e9f9ab3e29 100644 --- a/dipy/utils/fast_numpy.pxd +++ b/dipy/utils/fast_numpy.pxd @@ -52,7 +52,7 @@ cdef void random_perpendicular_vector( cpdef (double, double) random_point_within_circle( double r) noexcept nogil -cpdef double random() nogil +cpdef double random() noexcept nogil cpdef void seed(cnp.npy_uint32 s) nogil From a84fd917b0ba0254baa223618392e3f3d3a7ebca Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 2 Feb 2024 14:23:24 -0500 Subject: [PATCH 18/19] WIP - alpha version of paralle tracking --- dipy/tracking/fast_tracking.pxd | 14 ++-- dipy/tracking/fast_tracking.pyx | 120 +++++++++++++++++++++++++------- 2 files changed, 104 insertions(+), 30 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 7f07ddeaaf1..6b9d02a3baf 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -45,11 +45,15 @@ cdef int generate_local_streamline(double* seed, ProbabilisticTrackingParameters params, PmfGen pmf_gen) noexcept nogil - -cdef double* get_pmf(double* point, - PmfGen pmf_gen, - double pmf_threshold, - int pmf_len) noexcept nogil +cdef int trilinear_interpolate4d_c(double[:, :, :, :] data, + double* point, + double* result) noexcept nogil + +cdef int get_pmf(double* pmf, + double* point, + PmfGen pmf_gen, + double pmf_threshold, + int pmf_len) noexcept nogil cdef class ProbabilisticTrackingParameters(TrackingParameters): diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index ca2b0872e67..4223f247f5b 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -3,10 +3,12 @@ # cython: wraparound=False # cython: Nonecheck=False +cimport cython from cython.parallel import prange import numpy as np cimport numpy as cnp +#from dipy.core.interpolation cimport trilinear_interpolate4d_c from dipy.direction.pmf cimport PmfGen from dipy.tracking.stopping_criterion cimport StoppingCriterion from dipy.utils.fast_numpy cimport (copy_point, cumsum, norm, normalize, @@ -20,6 +22,8 @@ from dipy.tracking.stopping_criterion cimport (StreamlineStatus, INVALIDPOINT) from libc.stdlib cimport malloc, free +from libc.math cimport floor +from libc.time cimport time, time_t #need cpdef to access cdef functions? @@ -33,13 +37,18 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, double[:,:,:] streamlines_arr double[:] status_arr cnp.npy_intp i, j + time_t t1, t2 # temporary array to store generated streamlines streamlines_arr = np.array(np.zeros((_len, params.max_len * 2 + 1, 3)), order='C') status_arr = np.array(np.zeros((_len)), order='C') - + print("1") + t1 = time(NULL) generate_tractogram_c(seed_positons, seed_directions, _len, sc, params, pmf_gen, &probabilistic_tracker, streamlines_arr, status_arr) + print("2") + t2 = time(NULL) + print(t2-t1) streamlines = [] for i in range(_len): # array to list @@ -69,7 +78,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, for i in prange(nbr_seeds, nogil=True): - # for i in range(nbr_seeds): + #for i in range(nbr_seeds): stream = malloc((params.max_len * 3 * 2 + 1) * sizeof(double)) # initialize to 0. It will be replaced when better handling various @@ -126,7 +135,6 @@ cdef int generate_local_streamline(double* seed, # update position for j in range(3): point[j] += voxdir[j] * params.inv_voxel_size[j] * params.step_size - copy_point(point, &stream[(params.max_len + i )* 3]) stream_status_forward = sc.check_point_c(point) @@ -135,7 +143,7 @@ cdef int generate_local_streamline(double* seed, stream_status_forward == OUTSIDEIMAGE): break - # backward tracking + # # backward tracking copy_point(seed, point) copy_point(direction, voxdir) for j in range(3): @@ -147,30 +155,93 @@ cdef int generate_local_streamline(double* seed, # update position for j in range(3): point[j] += voxdir[j] * params.inv_voxel_size[j] * params.step_size - copy_point(point, &stream[(params.max_len + i )* 3]) - stream_status_backward = sc.check_point_c(point) if (stream_status_backward == ENDPOINT or stream_status_backward == INVALIDPOINT or stream_status_backward == OUTSIDEIMAGE): break - # need to handle stream status + # # need to handle stream status return 0 #stream_status -cdef double* get_pmf(double* point, - PmfGen pmf_gen, - double pmf_threshold, - int pmf_len) noexcept nogil: +@cython.boundscheck(False) +@cython.wraparound(False) +cdef int trilinear_interpolate4d_c( + double[:, :, :, :] data, + double* point, + double* result) noexcept nogil: + """Tri-linear interpolation along the last dimension of a 4d array + + Parameters + ---------- + point : 1d array (3,) + 3 doubles representing a 3d point in space. If point has integer values + ``[i, j, k]``, the result will be the same as ``data[i, j, k]``. + data : 4d array + Data to be interpolated. + result : 1d array + The result of interpolation. Should have length equal to the + ``data.shape[3]``. + Returns + ------- + err : int + 0 : successful interpolation. + -1 : point is outside the data area, meaning round(point) is not a + valid index to data. + -2 : point has the wrong shape + -3 : shape of data and result do not match + + """ + cdef: + cnp.npy_intp flr, N + double w, rem + cnp.npy_intp index[3][2] + double weight[3][2] + + #if data.shape[3] != result.shape[0]: + # return -3 + + for i in range(3): + if point[i] < -.5 or point[i] >= (data.shape[i] - .5): + return -1 + + flr = floor(point[i]) + rem = point[i] - flr + + index[i][0] = flr + (flr == -1) + index[i][1] = flr + (flr != (data.shape[i] - 1)) + weight[i][0] = 1 - rem + weight[i][1] = rem + + N = data.shape[3] + for i in range(N): + result[i] = 0 + + for i in range(2): + for j in range(2): + for k in range(2): + w = weight[0][i] * weight[1][j] * weight[2][k] + for L in range(N): + result[L] += w * data[index[0][i], index[1][j], + index[2][k], L] + return 0 + + +cdef int get_pmf(double* pmf, + double* point, + PmfGen pmf_gen, + double pmf_threshold, + int pmf_len) noexcept nogil: cdef: cnp.npy_intp i - double* pmf double absolute_pmf_threshold double max_pmf=0 - pmf = pmf_gen.get_pmf_c(point) + if trilinear_interpolate4d_c(pmf_gen.data, point, pmf): + return 1 + for i in range(pmf_len): if pmf[i] > max_pmf: max_pmf = pmf[i] @@ -179,7 +250,8 @@ cdef double* get_pmf(double* point, for i in range(pmf_len): if pmf[i] < absolute_pmf_threshold: pmf[i] = 0.0 - return pmf + + return 0 cdef int probabilistic_tracker(double* point, @@ -187,21 +259,17 @@ cdef int probabilistic_tracker(double* point, ProbabilisticTrackingParameters params, PmfGen pmf_gen) noexcept nogil: cdef: - cnp.npy_intp i, idx + cnp.npy_intp i, idx, res double* newdir double* pmf double last_cdf, cos_sim - double max_pmf=0 - double absolute_pmf_threshold - - pmf = get_pmf(point, pmf_gen, params.pmf_threshold, params.pmf_len) - - - - - + pmf = malloc(params.pmf_len * sizeof(double)) + if get_pmf(pmf, point, pmf_gen, params.pmf_threshold, params.pmf_len): + free(pmf) + return 1 if norm(direction) == 0: + free(pmf) return 1 normalize(direction) @@ -217,10 +285,11 @@ cdef int probabilistic_tracker(double* point, cumsum(pmf, pmf, params.pmf_len) last_cdf = pmf[params.pmf_len - 1] if last_cdf == 0: + free(pmf) return 1 idx = where_to_insert(pmf, random() * last_cdf, params.pmf_len) - + #idx=0 ######################################################################### REMOVE newdir = ¶ms.vertices[idx][0] # Update direction if (direction[0] * newdir[0] @@ -232,6 +301,7 @@ cdef int probabilistic_tracker(double* point, newdir[1] = newdir[1] * -1 newdir[2] = newdir[2] * -1 copy_point(newdir, direction) + free(pmf) return 0 #get_direction_c of the DG From 4fd1b4299b0d18e87a252f3a4d6f033981ad9378 Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Thu, 8 Feb 2024 09:55:01 -0500 Subject: [PATCH 19/19] add functor management --- dipy/tracking/fast_tracking.pxd | 16 ++++++++-------- dipy/tracking/fast_tracking.pyx | 23 +++++++++++------------ 2 files changed, 19 insertions(+), 20 deletions(-) diff --git a/dipy/tracking/fast_tracking.pxd b/dipy/tracking/fast_tracking.pxd index 6b9d02a3baf..e62aa5d9d7a 100644 --- a/dipy/tracking/fast_tracking.pxd +++ b/dipy/tracking/fast_tracking.pxd @@ -15,16 +15,16 @@ cdef class TrackingParameters(): ctypedef int (*func_ptr)(double* point, double* direction, ProbabilisticTrackingParameters, - PmfGen) nogil + PmfGen) noexcept nogil -cpdef list generate_tractogram(double[:,::1] seed_positons, +cpdef list generate_tractogram(double[:,::1] seed_positions, double[:,::1] seed_directions, StoppingCriterion sc, ProbabilisticTrackingParameters params, PmfGen pmf_gen) -cdef int generate_tractogram_c(double[:,::1] seed_positons, +cdef int generate_tractogram_c(double[:,::1] seed_positions, double[:,::1] seed_directions, int nbr_seeds, StoppingCriterion sc, @@ -79,7 +79,7 @@ cdef int deterministic_maximum_tracker(double* point, DeterministicTrackingParameters params, PmfGen pmf_gen,) -cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): +cdef class ParallelTransportTrackingParameters(ProbabilisticTrackingParameters): cdef: double angular_separation double data_support_exponent @@ -103,7 +103,7 @@ cdef class ParalleTransportTrackingParameters(ProbabilisticTrackingParameters): int rejection_sampling_nbr_sample -cdef int paralle_transport_tracker(double* point, - double* direction, - ParalleTransportTrackingParameters params, - PmfGen pmf_gen) \ No newline at end of file +cdef int parallel_transport_tracker(double* point, + double* direction, + ParallelTransportTrackingParameters params, + PmfGen pmf_gen) \ No newline at end of file diff --git a/dipy/tracking/fast_tracking.pyx b/dipy/tracking/fast_tracking.pyx index 4223f247f5b..0e2295070cc 100644 --- a/dipy/tracking/fast_tracking.pyx +++ b/dipy/tracking/fast_tracking.pyx @@ -27,13 +27,13 @@ from libc.time cimport time, time_t #need cpdef to access cdef functions? -cpdef list generate_tractogram(double[:,::1] seed_positons, +cpdef list generate_tractogram(double[:,::1] seed_positions, double[:,::1] seed_directions, StoppingCriterion sc, ProbabilisticTrackingParameters params, PmfGen pmf_gen): cdef: - cnp.npy_intp _len=seed_positons.shape[0] + cnp.npy_intp _len=seed_positions.shape[0] double[:,:,:] streamlines_arr double[:] status_arr cnp.npy_intp i, j @@ -44,7 +44,7 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, status_arr = np.array(np.zeros((_len)), order='C') print("1") t1 = time(NULL) - generate_tractogram_c(seed_positons, seed_directions, _len, sc, params, + generate_tractogram_c(seed_positions, seed_directions, _len, sc, params, pmf_gen, &probabilistic_tracker, streamlines_arr, status_arr) print("2") t2 = time(NULL) @@ -64,7 +64,7 @@ cpdef list generate_tractogram(double[:,::1] seed_positons, return streamlines -cdef int generate_tractogram_c(double[:,::1] seed_positons, +cdef int generate_tractogram_c(double[:,::1] seed_positions, double[:,::1] seed_directions, int nbr_seeds, StoppingCriterion sc, @@ -76,7 +76,6 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, cdef: cnp.npy_intp i, j, k - for i in prange(nbr_seeds, nogil=True): #for i in range(nbr_seeds): stream = malloc((params.max_len * 3 * 2 + 1) * sizeof(double)) @@ -86,7 +85,7 @@ cdef int generate_tractogram_c(double[:,::1] seed_positons, for j in range(params.max_len * 3 * 2 + 1): stream[j] = 0 - status[i] = generate_local_streamline(&seed_positons[i][0], + status[i] = generate_local_streamline(&seed_positions[i][0], &seed_directions[i][0], stream, tracker, @@ -130,7 +129,7 @@ cdef int generate_local_streamline(double* seed, # forward tracking stream_status_forward = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(&point[0], &voxdir[0], params, pmf_gen): # , pmf_gen_func): + if tracker(&point[0], &voxdir[0], params, pmf_gen): # , pmf_gen_func): break # update position for j in range(3): @@ -150,7 +149,7 @@ cdef int generate_local_streamline(double* seed, voxdir[j] = voxdir[j] * -1 stream_status_backward = TRACKPOINT for i in range(1, params.max_len): - if probabilistic_tracker(&point[0], &voxdir[0], params, pmf_gen): + if tracker(&point[0], &voxdir[0], params, pmf_gen): break # update position for j in range(3): @@ -316,10 +315,10 @@ cdef int deterministic_maximum_tracker(double* point, return 1 #get_direction_c of the DG -cdef int paralle_transport_tracker(double* point, - double* direction, - ParalleTransportTrackingParameters params, - PmfGen pmf_gen): +cdef int parallel_transport_tracker(double* point, + double* direction, + ParallelTransportTrackingParameters params, + PmfGen pmf_gen): # update point and dir with new position and direction # return 1 if the propagation failed.