Skip to content

Comments

[ENH] PyAFQ in the superficial white matter#103

Merged
arokem merged 167 commits intotractometry:mainfrom
36000:dam
Dec 18, 2025
Merged

[ENH] PyAFQ in the superficial white matter#103
arokem merged 167 commits intotractometry:mainfrom
36000:dam

Conversation

@36000
Copy link
Collaborator

@36000 36000 commented Jun 23, 2025

This is a few things I have been trying out to modernize pyAFQ on multi-shell images, as well as some random things.

These are used:

  1. Require T1; use Brainchop to make brain mask, synthseg to segment into CSF, GM, WM; generates WM/GM interface highly reliably using some basic scikit-image functions and the 3T segmentations.
  2. I am adding an option to use IsolationForest from scikit-learn for cleaning and making that the default for VOF and callosal bundles over the old system. This makes scikit-learn a dependency. It does not assume a tube-like shape for a given bundle, instead finding outlier nodes and removing streamlines with too many outlier nodes.
  3. I made it so mahalanobis cleaning by default only considers the middle 60% of the bundle. This means streamlines are allowed to deviate in the starting and ending 20% of the bundle. This is useful for allowing more diverse endpoints.
  4. Made some changes to VOF segmentation. There is still more work to be done in a separate PR.
  5. Streamlined API for parallelization. Made parallelization with Ray turned on by default.
  6. Made PFT, ACT, AODF CSD and WMGMI seeding the default over local tracking and uniform seeding in the white matter. More work needs to be done on smarter WMGMI seeding, but this is a good start.
  7. I am vendorizing here some code from scilpy for doing unified filtering to get asymmetric ODFs. I have modified this to also use numba and ray, and it is efficient enough, ~20 min. could be an option for anyone interested in endpoints.
  8. Use osqpy directly to fit MSMT

@36000 36000 changed the title [WIP/ENH] Retooling Multishell pipeline [WIP/ENH] PyAFQ in the superficial white matter Jun 30, 2025
@36000
Copy link
Collaborator Author

36000 commented Jun 30, 2025

@36000
Copy link
Collaborator Author

36000 commented Jul 8, 2025

Some interesting ORs found using seeding by endpoints on wmgmi using 2M seeds
Screenshot 2025-07-08 at 2 08 44 PM
Screenshot 2025-07-08 at 2 05 49 PM

@36000
Copy link
Collaborator Author

36000 commented Jul 22, 2025

Added an interesting output. For each bundle, you can see the distance to the nearest endpoint for every voxel in the gray matter. The distances are in millimeters. This could be thresholded to get a boolean endpoint mask for each bundle.

endpoint_distance_visualization.mov

@36000 36000 changed the title [WIP/ENH] PyAFQ in the superficial white matter [ENH] PyAFQ in the superficial white matter Jul 29, 2025
@36000 36000 requested a review from Copilot July 30, 2025 21:42
Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

This PR introduces modern multi-shell DWI processing capabilities to pyAFQ, along with improvements to tractography, segmentation, and bundle recognition workflows. The primary focus is enhancing analysis of superficial white matter through better seeding, tracking, and cleaning methods.

Key Changes:

  • Require T1-weighted images: All workflows now require T1w data for tissue segmentation and interface generation
  • Enhanced tracking defaults: Switch to PFT tracking with WM/GM interface seeding and ACT stopping criteria by default
  • Improved bundle cleaning: Add IsolationForest cleaning as default for VOF and callosal bundles, plus enhanced Mahalanobis cleaning

Reviewed Changes

Copilot reviewed 43 out of 46 changed files in this pull request and generated 6 comments.

Show a summary per file
File Description
setup.cfg Add scikit-learn, numba, and osqp dependencies
AFQ/tasks/data.py Add T1w tissue segmentation, DAM fitting, MSMT CSD, and asymmetric ODF processing
AFQ/tasks/tractography.py Update default tracking to PFT with WM/GM interface seeding
AFQ/recognition/cleaning.py Add IsolationForest cleaning method and orientation-based Mahalanobis cleaning
AFQ/api/bundle_dict.py Update VOF and callosal bundle definitions with new cleaning methods
AFQ/models/ Add new models for MSMT, DAM fitting, asymmetric filtering, and WM/GM interface generation
Comments suppressed due to low confidence (2)

AFQ/models/msmt.py:328

  • [nitpick] The parameter name 'use_osqppy' is unclear. Consider renaming to 'use_osqp' to better reflect that it uses the OSQP solver.
         use_osqppy=True,

@36000
Copy link
Collaborator Author

36000 commented Aug 5, 2025

Now with brainchop and isolation forest cleaning, here are the callosal bundles from 1million seeds on example HBN subject
callosal

@36000
Copy link
Collaborator Author

36000 commented Aug 7, 2025

@arokem right now I am just trying to get the documentation to finish without running out of memory, I think due to increases in parallelization. Anyways, this is otherwise ready for review!

@arokem
Copy link
Member

arokem commented Aug 7, 2025

I'll get right to it! Shouldn't take me more than...

Screenshot 2025-08-07 at 11 26 34 AM

... 6-8 months 😄

Copy link
Member

@arokem arokem left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First batch of comments

@arokem
Copy link
Member

arokem commented Aug 19, 2025

I am experiencing an interesting regression with this PR, where calling:


my_afq = GroupAFQ(
    bids_path=study_dir,
    preproc_pipeline="qsiprep",
    brain_mask_definition=brain_mask_definition,
    tracking_params={"n_seeds": 4,
                     "directions": "prob",
                     "odf_model": "CSD",
                     "seed_mask": RoiImage()},
    bundle_info=bundles)

proceeds to use only 4 seeds for tracking (as though I am passing "random_seeds": True as part of the tracking_params dictionary. Is that an intentional consequence of these changes?

@36000
Copy link
Collaborator Author

36000 commented Aug 19, 2025

This is not a regression, its a change in defaults. The old defaults were:

"random_seeds": False
"n_seeds": 1

New defaults are:

"random_seeds": True
"n_seeds": 2000000

I think this is more intuitive but old users might get confused initially.

@arokem
Copy link
Member

arokem commented Aug 19, 2025

This old user certainly did 😆

@arokem
Copy link
Member

arokem commented Sep 19, 2025

Hey @36000 : when you get a chance, could you please rebase this on main? Thanks!

@arokem
Copy link
Member

arokem commented Sep 19, 2025

Interestingly, with this branch, we found that installing pyAFQ into a clean environment doesn't work (under some circumstances?). Instead, both @johndromero and I see this:

Obtaining file:///Users/arokem/tmp/afq_dam_test/pyAFQ
  Installing build dependencies ... done
  Checking if build backend supports build_editable ... done
  Getting requirements to build editable ... done
  Preparing editable metadata (pyproject.toml) ... error
  error: subprocess-exited-with-error
  
  × Preparing editable metadata (pyproject.toml) did not run successfully.
  │ exit code: 1
  ╰─> [71 lines of output]
      /private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools_scm/_integration/version_inference.py:51: UserWarning: version of pyAFQ already set
        warnings.warn(self.message)
      Traceback (most recent call last):
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/expand.py", line 71, in __getattr__
          return next(
                 ^^^^^
      StopIteration
      
      The above exception was the direct cause of the following exception:
      
      Traceback (most recent call last):
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/expand.py", line 185, in read_attr
          value = getattr(StaticModule(module_name, spec), attr_name)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/expand.py", line 77, in __getattr__
          raise AttributeError(f"{self.name} has no attribute {attr}") from e
      AttributeError: AFQ has no attribute __version__
      
      During handling of the above exception, another exception occurred:
      
      Traceback (most recent call last):
        File "/Users/arokem/miniforge3/envs/afq_dam_test/lib/python3.11/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 389, in <module>
          main()
        File "/Users/arokem/miniforge3/envs/afq_dam_test/lib/python3.11/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 373, in main
          json_out["return_val"] = hook(**hook_input["kwargs"])
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/Users/arokem/miniforge3/envs/afq_dam_test/lib/python3.11/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 209, in prepare_metadata_for_build_editable
          return hook(metadata_directory, config_settings)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/build_meta.py", line 478, in prepare_metadata_for_build_editable
          return self.prepare_metadata_for_build_wheel(
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/build_meta.py", line 374, in prepare_metadata_for_build_wheel
          self.run_setup()
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/build_meta.py", line 317, in run_setup
          exec(code, locals())
        File "<string>", line 33, in <module>
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/__init__.py", line 115, in setup
          return distutils.core.setup(**attrs)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/_distutils/core.py", line 160, in setup
          dist.parse_config_files()
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/dist.py", line 752, in parse_config_files
          setupcfg.parse_configuration(
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 188, in parse_configuration
          meta.parse()
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 502, in parse
          section_parser_method(section_options)
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 477, in parse_section
          self[name] = value
          ~~~~^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 294, in __setitem__
          parsed = self.parsers.get(option_name, lambda x: x)(value)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 600, in _parse_version
          return expand.version(self._parse_attr(value, self.package_dir, self.root_dir))
                                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/setupcfg.py", line 419, in _parse_attr
          return expand.read_attr(attr_desc, package_dir, root_dir)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/expand.py", line 190, in read_attr
          module = _load_spec(spec, module_name)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        File "/private/var/folders/px/t5ym05pn5k16nz97mzcpqqkr0000gn/T/pip-build-env-2rgpqld6/overlay/lib/python3.11/site-packages/setuptools/config/expand.py", line 211, in _load_spec
          spec.loader.exec_module(module)
        File "<frozen importlib._bootstrap_external>", line 940, in exec_module
        File "<frozen importlib._bootstrap>", line 241, in _call_with_frames_removed
        File "/Users/arokem/tmp/afq_dam_test/pyAFQ/AFQ/__init__.py", line 1, in <module>
          from .version import version as __version__  # noqa
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      ModuleNotFoundError: No module named 'AFQ.version'
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
error: metadata-generation-failed

× Encountered error while generating package metadata.
╰─> See above for output.

note: This is an issue with the package mentioned above, not pip.
hint: See above for details.

I think that it's related to changes in setup/pyproject, but not sure.

If I try to install this on top of an already-existing installation (i.e., in an environment where pyAFQ was previously installed) then no such issue. Not sure why this works in the CI and not locally, but thought I'd raise this here so we can debug.

@36000 36000 force-pushed the dam branch 4 times, most recently from ea5ed84 to a756209 Compare September 19, 2025 22:56
num_chunks=True)
rng_seed=2025,
odf_model="csd_aodf",
trx=False)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is trx switched off? Do we want to add some text explaining these settings (particularly the "aodf" part)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should be on, and odf_model="csd_aodf" is no longer necessary (it is the default). I must have changed this during testing at some point

@arokem
Copy link
Member

arokem commented Dec 16, 2025

The first tutorial example now runs on my primitive mac (no onnxruntime), but I end up with this:

INFO:AFQ:Calculating _desc-slCount_tractography.csv
INFO:AFQ:_desc-slCount_tractography.csv completed. Saving to /Users/arokem/AFQ_data/HBN/derivatives/afq/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/stats/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_desc-slCount_tractography.csv
Traceback (most recent call last):
  File "/Users/arokem/source/pyAFQ/examples/tutorial_examples/plot_001_group_afq_api.py", line 268, in <module>
    raise ValueError((
ValueError: Small number of streamlines found for bundle(s):
                                 n_streamlines
Left Anterior Thalamic                       0
Right Anterior Thalamic                      0
Left Cingulum Cingulate                      0
Right Cingulum Cingulate                     0
Left Corticospinal                           0
Right Corticospinal                          0
Left Inferior Fronto-occipital               0
Right Inferior Fronto-occipital              0
Left Inferior Longitudinal                   0
Right Inferior Longitudinal                  0
Left Superior Longitudinal                  20
Right Superior Longitudinal                  0
Left Arcuate                                 0
Right Arcuate                                0
Left Uncinate                                0
Right Uncinate                               0
Left Posterior Arcuate                       0
Right Posterior Arcuate                      0
Left Vertical Occipital                     20
Right Vertical Occipital                     0
Callosum Anterior Frontal                    0
Callosum Motor                               0
Callosum Occipital                           0
Callosum Orbital                             0
Callosum Posterior Parietal                  0
Callosum Superior Frontal                    0
Callosum Superior Parietal                   0
Callosum Temporal                            0
Total Recognized                            40

@36000
Copy link
Collaborator Author

36000 commented Dec 16, 2025

@arokem try again. I introduced a bug in the last commit while I was trying to fix other things with PVE resampling. Now it should be working

sh_coeff = nib.load(csd_params).get_fdata()

logger.info("Applying unified filtering to generate "
"asymmetric CSD CSD ODFs...")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"asymmetric CSD CSD ODFs...")
"asymmetric CSD ODFs...")

@arokem
Copy link
Member

arokem commented Dec 16, 2025

Yep, I can run the first example! I did notice that the html output no longer has an anatomical image embedded in it. I think that we had the T1 in there at some point?

@arokem
Copy link
Member

arokem commented Dec 16, 2025

I guess we can sort that out on a separate PR after we merge this one. Looks like there are still some test failures on the CI, though.

@36000
Copy link
Collaborator Author

36000 commented Dec 17, 2025

@arokem I think I fixed the T1 not showing up

@arokem
Copy link
Member

arokem commented Dec 18, 2025

Yep! It works now.

I think we might be ready to merge this PR. WDYT?

@arokem
Copy link
Member

arokem commented Dec 18, 2025

Oh sorry, I think I am still running into some issues. In particular, I tried to run the pyafq2 example and running into some errors. I'll work to figure out what these errors and report back.

@arokem
Copy link
Member

arokem commented Dec 18, 2025

OK - looks like when I run things in the "pyAFQ2" mode I get the onnxruntime errors, unless I also provide a brain_mask_definition. I think that's probably fine, but I might suggest to use the HBN data, and to add the brain_mask_definition in the example. Something like this:


diff --git a/examples/howto_examples/pyafq_2.py b/examples/howto_examples/pyafq_2.py
index b3f9b80d..e22e0b85 100644
--- a/examples/howto_examples/pyafq_2.py
+++ b/examples/howto_examples/pyafq_2.py
@@ -64,7 +64,7 @@ pve = afm.PVEImages(
 # 3. it must be lateral to the inferior fronto-occipital fasciculus
 #    instead of the inferior longitudinal fasciculus;
 # 4. cleaning has been changed: there is now mahalanobis cleaning on
-#    orientation, and isolation forest cleaning instead of mahalanobis for 
+#    orientation, and isolation forest cleaning instead of mahalanobis for
 #    distance.
 # Additionally, in the new version, the inferior endpoints of the
 # corticospinal tracts (CST) were removed.
@@ -214,12 +214,22 @@ bundle_info = abd.default18_bd() + \
 # Bundle definitions for VOF, pAF, and CST to use the old definitions;
 # Callosal bundles to use mahalanobis cleaning.
 
+bids_path = afd.fetch_hbn_preproc(
+    ["NDARAA948VFH"],
+    clear_previous_afq="all")[1]
+
+brain_mask_definition = afm.ImageFile(
+    suffix="mask", filters={"desc": "brain", "scope": "qsiprep"})
+
+
 myafq = GroupAFQ(
-    bids_path=op.join(afd.afq_home, 'stanford_hardi'),
-    dwi_preproc_pipeline='vistasoft',
-    t1_preproc_pipeline='freesurfer',
+    bids_path=bids_path,
+    dwi_preproc_pipeline='qsiprep',
+    t1_preproc_pipeline='qsiprep',
     tracking_params=tracking_params,
     pve=pve,
+    brain_mask_definition=brain_mask_definition,
+    participant_labels=["NDARAA948VFH"],
     bundle_info=bundle_info)
 
 myafq.export_all()

@arokem
Copy link
Member

arokem commented Dec 18, 2025

Or potentially at least as a commented-out block of code:

diff --git a/examples/howto_examples/pyafq_2.py b/examples/howto_examples/pyafq_
2.py
index b3f9b80d..0d7397ca 100644
--- a/examples/howto_examples/pyafq_2.py
+++ b/examples/howto_examples/pyafq_2.py
@@ -64,7 +64,7 @@ pve = afm.PVEImages(
 # 3. it must be lateral to the inferior fronto-occipital fasciculus
 #    instead of the inferior longitudinal fasciculus;
 # 4. cleaning has been changed: there is now mahalanobis cleaning on
-#    orientation, and isolation forest cleaning instead of mahalanobis for 
+#    orientation, and isolation forest cleaning instead of mahalanobis for
 #    distance.
 # Additionally, in the new version, the inferior endpoints of the
 # corticospinal tracts (CST) were removed.
@@ -214,12 +214,32 @@ bundle_info = abd.default18_bd() + \
 # Bundle definitions for VOF, pAF, and CST to use the old definitions;
 # Callosal bundles to use mahalanobis cleaning.
 
+bids_path = afd.fetch_hbn_preproc(
+    ["NDARAA948VFH"],
+    clear_previous_afq="all")[1]
+
+################################################################
+# .. note::
+#
+# In pyAFQ 2.x, if a brain mask was not provided, it was automatically calculated
+# using a median Otsu algorithm on the b0 image. In pyAFQ 3.x, a brain mask
+# would instead be calculated using the mindgrab neural network model. However,
+# this may not work on some older computers/os. In which case, the user should fall
+# back to provide a brain mask. For example, by specifying::
+#
+#    brain_mask_definition = afm.ImageFile(
+#         suffix="mask", filters={"desc": "brain", "scope": "qsiprep"})
+
+
 myafq = GroupAFQ(
-    bids_path=op.join(afd.afq_home, 'stanford_hardi'),
-    dwi_preproc_pipeline='vistasoft',
-    t1_preproc_pipeline='freesurfer',
+    bids_path=bids_path,
+    dwi_preproc_pipeline='qsiprep',
+    t1_preproc_pipeline='qsiprep',
     tracking_params=tracking_params,
     pve=pve,
+    # Only pass a brain mask here if you defined it above.
+    # brain_mask_definition=brain_mask_definition,

@36000
Copy link
Collaborator Author

36000 commented Dec 18, 2025

@arokem I made it so the pyafq 2 example uses a brain mask based on the freesurfer segmentation. This excludes the CSF, but seems to work.

@arokem
Copy link
Member

arokem commented Dec 18, 2025

This is now ready to be merged, I believe. We'll continue to polish and work on remaining items for version 3.0 on separate PRs.

@arokem arokem merged commit 6fa4987 into tractometry:main Dec 18, 2025
5 of 9 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants