[WIP] NF: Population specific bundle atlas creation #2425

drombas · 2021-07-26T10:44:50Z

Hi! This is part of GSoC 2021.

Goal

Given the segmented bundles of several subjects as input, compute a bundle atlas that is representative of that population.

Implementation

The input bundles are combined in pairs following a tree structure until a single bundle is obtained:

The main logic is implemented in atlasing.bundles.compute_bundle_atlas() and consists of these steps:

Check input arguments and retrieve the list of bundles and subjects to be processed
Preprocess the bundles (set the number of points per streamline, remove bundles with few streamlines...)
Create the tree structure with atlasing.bundles.get_pairwise_tree()
Build the atlas by following the tree. At each pair, bundles are aligned and then combined by atlasing.bundles.combine_bundles(), which implements several methods to match the streamlines and combine them into a single bundle.

Comments
This is WIP as we are still testing different approaches and need to reach a final decision on:

Which input file formats to support (currently it only supports trk)
Which bundle combination methods to support

List of tasks:

Create atlasing module and first tests
Extend the coverage of atlasing module tests
Add tests for the workflow
Polish docstrings and comments
Create an example/tutorial

pep8speaks · 2021-07-26T10:44:55Z

Hello @drombas, Thank you for updating !

Cheers ! There are no PEP8 issues in this Pull Request. 🍻

Comment last updated at 2021-08-17 10:47:06 UTC

codecov · 2021-07-26T11:38:58Z

Codecov Report

Merging #2425 (3e35297) into master (cf5f4f9) will increase coverage by 0.23%.
The diff coverage is 96.84%.

@@            Coverage Diff             @@
##           master    #2425      +/-   ##
==========================================
+ Coverage   84.49%   84.73%   +0.23%     
==========================================
  Files         125      127       +2     
  Lines       16624    16909     +285     
  Branches     2702     2769      +67     
==========================================
+ Hits        14046    14327     +281     
- Misses       1911     1913       +2     
- Partials      667      669       +2

Impacted Files	Coverage Δ
dipy/io/utils.py	`87.66% <78.78%> (+1.06%)`	⬆️
dipy/workflows/atlasing.py	`90.90% <90.90%> (ø)`
dipy/atlasing/bundles.py	`99.58% <99.58%> (ø)`
dipy/data/fetcher.py	`42.37% <100.00%> (+0.24%)`	⬆️

skoudoro · 2021-08-02T15:00:15Z

Thank you @drombas for this PR. Let us know when it is ready for a first review. Thank you

drombas · 2021-08-02T16:37:51Z

Hi @skoudoro ! Tests finally look mostly ok so yes, please, have a look at it!

arokem

This looks really interesting. I have a couple of small comments.

dipy/atlasing/bundles.py

arokem · 2021-08-03T21:57:23Z

dipy/io/utils.py

@@ -420,3 +425,71 @@ def read_img_arr_or_path(data, affine=None):
            affine = data.affine
        data = data.get_fdata()
    return data, affine
+
+
+def show_bundles(bundles, fname, colors=None, opacity=None, linewidth=None):


I wonder whether this shouldn't go into dipy.viz

I putted there cause I use it to save the output as png files but makes sense to move it to dipy.viz. @skoudoro what do you think?

dipy/io/utils.py

arokem · 2021-08-04T02:44:42Z

dipy/atlasing/bundles.py

+        raise ValueError("Subjects cannot be duplicated")
+
+    print(str(len(subjects)) + " subjects to be processed:")
+    print(subjects)


This can be pretty long, no?

Yes but I found it useful as a safety check. Using logger.debug() might be a good option?

BramshQamar

Hi @drombas,

Great work! Thank you for creating the PR. I have added my first review.

Bramsh

BramshQamar · 2021-08-17T20:11:53Z

dipy/atlasing/bundles.py

@@ -0,0 +1,506 @@
+"""Atlasing module: utilities to compute population specific bundle atlases.


population-specific streamline-based bundle atlases.

BramshQamar · 2021-08-17T20:13:00Z

dipy/atlasing/bundles.py

+"""Atlasing module: utilities to compute population specific bundle atlases.
+
+Available functions:
+    get_pairwise_tree: computes the matching pairs for a given number of items.


items or bundles?

BramshQamar · 2021-08-17T20:16:26Z

dipy/atlasing/bundles.py

+    combine_bundle: combines two bundles into a single one using different
+    streamline combination methods.
+
+    compute_bundle_atlas: given a list of input bundles computes the population


Does this create a whole-brain atlas from given bundles? eg: merging all bundles to create one whole brain atlas. Or does it create one atlas bundle per bundle type? I think it's the second one. Maybe just rephrase it.

BramshQamar · 2021-08-17T20:18:49Z

dipy/atlasing/bundles.py

+    """Pairwise tree structure calculation.
+
+    Constructs a pairwise tree by randomly matching the indexes of a given
+    number of items. The computed index structure is intended to be used for


Let's just use the word 'bundle' since it is a bundle-specific atlas or do you think this is more generic and could be used for other purposes? If it's the second, then maybe we need to move it to utils?

BramshQamar · 2021-08-17T20:23:07Z

dipy/atlasing/bundles.py

+        # If n_item is odd duplicate one of the others
+        if np.mod(n_item, 2) == 1:
+            # avoid removing again an item twice (index == 0)
+            lonely = np.random.randint(1, n_item)


How about the variable name single?

BramshQamar · 2021-08-17T21:03:50Z

dipy/workflows/atlasing.py

+from dipy.atlasing.bundles import compute_atlas_bundle
+
+
+class DiscreteBundleAtlasFlow(Workflow):


I think we should drop the word 'Discrete'.

BramshQamar · 2021-08-17T21:04:17Z

dipy/workflows/atlasing.py

+
+        Given several segmented bundles as input, compute the atlas by
+        combining the bundles pairwise.
+


Add some more details

BramshQamar · 2021-08-17T21:04:47Z

dipy/workflows/atlasing.py

+        distance : str, optional
+            Distance metric to be used to combine bundles. Default is 'mdf'.
+            The 'mdf_se' distance uses only start/end points of streamlines.
+        comb_method : str, optional


Same comment as in the other functions.

BramshQamar · 2021-08-17T21:08:02Z

dipy/atlasing/bundles.py

+                                            space=Space.RASMM)
+        save_trk(new_tractogram, file, bbox_valid_check=False)
+        return atlas, atlas_merged
+    return atlas


You can just return atlas, atlas_merged
atlas_merged can be None if there isn't one. This way the function output will be consistent (returns fixed number of objects).

Will solve the issue in the workflow as well where you have written the same function call twice with if and else conditions.

BramshQamar · 2021-08-17T21:10:21Z

dipy/workflows/atlasing.py

+                                     bundle_names=bundle_names,
+                                     model_bundle_dir=model_bundle_dir,
+                                     out_dir=out_dir,
+                                     merge_out=merge_out,


You can call it one time and just always have it return two outputs. See my comment above in compute_atlas_bundle function.

BramshQamar · 2021-09-03T19:06:43Z

dipy/atlasing/bundles.py

+                                                rm_small_clusters=1,
+                                                qbx_thr=[5])
+
+                points, offsets = unlist_streamlines(static)


Why unlist and relist here?

iPsych · 2022-11-16T13:56:00Z

When I run it and generate average tract using TRACULA training data (33 subjects), I found interesting error only happens in bilateral cab_PP (cingulum-angular bundle).

I guess it's due to the irregularity of each subject's tract shape, but cannot clearly figure out what happened.
The tracts also attached.

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [20], line 4
      2 for x in tracts:
      3   print(x)
----> 4   atlasing.bundles.compute_atlas_bundle('atlasing_test_flt/'+x+'/',model_bundle_dir='/Volumes/Jotunn/Suicide_Study/TRACULA_DIPY_Vis/atlasing_test_flt/rh.cab_PP/trc001',n_point=100,n_stream_min=5)
[rh.cab_PP.zip](https://github.com/dipy/dipy/files/10022300/rh.cab_PP.zip)

File ~/Library/jupyterlab-desktop/jlab_server/lib/python3.8/site-packages/dipy/atlasing/bundles.py:402, in compute_atlas_bundle(in_dir, subjects, group, mid_path, bundle_names, model_bundle_dir, out_dir, merge_out, save_temp, n_stream_min, n_stream_max, n_point, distance, comb_method, skip_pairs)
    399 streamlines = set_number_of_points(streamlines, n_point)
    401 if model_bundles is not None:
--> 402     streamlines, _, _, _ = slr_with_qbx(static=model_bundle,
    403                                         moving=streamlines,
    404                                         x0='affine',
    405                                         rm_small_clusters=1,
    406                                         qbx_thr=[5])
    407     header = header_model
    409 file = f'{step_dir}/bundle_{i}_prev_{sub}'

File ~/Library/jupyterlab-desktop/jlab_server/lib/python3.8/site-packages/dipy/align/streamlinear.py:980, in slr_with_qbx(static, moving, x0, rm_small_clusters, maxiter, select_random, verbose, greater_than, less_than, qbx_thr, nb_pts, progressive, rng, num_threads)
    977 rstreamlines2 = set_number_of_points(rstreamlines2, nb_pts)
    978 rstreamlines2._data.astype('f4')
--> 980 cluster_map2 = qbx_and_merge(rstreamlines2, thresholds=qbx_thr, rng=rng)
    982 qb_centroids2 = remove_clusters_by_size(cluster_map2, rm_small_clusters)
    984 if verbose:

File ~/Library/jupyterlab-desktop/jlab_server/lib/python3.8/site-packages/dipy/segment/clustering.py:746, in qbx_and_merge(streamlines, thresholds, nb_pts, select_randomly, rng, verbose)
    742 qbx_merge = QuickBundlesX([thresholds[-1]],
    743                           metric=AveragePointwiseEuclideanMetric())
    745 final_level = len(thresholds)
--> 746 len_qbx_fl = len(qbx_clusters.get_clusters(final_level))
    747 qbx_ordering_final = rng.choice(len_qbx_fl, len_qbx_fl, replace=False)
    749 qbx_merged_cluster_map = qbx_merge.cluster(
    750     qbx_clusters.get_clusters(final_level).centroids,
    751     ordering=qbx_ordering_final).get_clusters(1)

AttributeError: 'ClusterMapCentroid' object has no attribute 'get_clusters'

drombas · 2022-11-16T20:17:54Z

Hi @iPsych,
I suspect it is related with a bundle having very few streamlines (trc30 maybe).
Anyways, let me check it with my last version of this method (this PR is quite old)

drombas added 2 commits July 26, 2021 11:42

NF: add atlasing module and bundle atlas creation workflow

e322d37

TEST: add tests for bundle atlasing

6a9279f

drombas added 4 commits August 2, 2021 12:09

TEST: extend coverage of bundle atlasing

105c3c7

STYLE: improve readability of bundle atlasing

e9634d8

TEST: add pandas and fury as optional packages

8dd523d

RF: move function to save bundles in png to io.utils

db45ca5

drombas added 2 commits August 13, 2021 09:58

BF correct workflow launcher

6e52726

RF: optimize combine_bundles function

c9dda4c

arokem reviewed Aug 13, 2021

View reviewed changes

Improve exception messages and logging

3e35297

BramshQamar reviewed Aug 17, 2021

View reviewed changes

BramshQamar reviewed Sep 3, 2021

View reviewed changes

drombas mentioned this pull request Mar 7, 2022

NF: Unbiased groupwise linear bundle registration #2553

Merged

skoudoro force-pushed the master branch 7 times, most recently from 1419292 to ca6268a Compare December 8, 2023 22:25

skoudoro force-pushed the master branch 3 times, most recently from 5935e1e to 765963e Compare January 24, 2024 19:24

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] NF: Population specific bundle atlas creation #2425

[WIP] NF: Population specific bundle atlas creation #2425

drombas commented Jul 26, 2021 •

edited

pep8speaks commented Jul 26, 2021 •

edited

codecov bot commented Jul 26, 2021 •

edited

skoudoro commented Aug 2, 2021

drombas commented Aug 2, 2021

arokem left a comment

arokem Aug 3, 2021

drombas Aug 17, 2021

arokem Aug 4, 2021

drombas Aug 17, 2021

BramshQamar left a comment

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021 •

edited

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Aug 17, 2021

BramshQamar Sep 3, 2021

iPsych commented Nov 16, 2022

drombas commented Nov 16, 2022

		@@ -0,0 +1,506 @@
		"""Atlasing module: utilities to compute population specific bundle atlases.

		from dipy.atlasing.bundles import compute_atlas_bundle


		class DiscreteBundleAtlasFlow(Workflow):


		Given several segmented bundles as input, compute the atlas by
		combining the bundles pairwise.

[WIP] NF: Population specific bundle atlas creation #2425

Are you sure you want to change the base?

[WIP] NF: Population specific bundle atlas creation #2425

Conversation

drombas commented Jul 26, 2021 • edited

Goal

Implementation

List of tasks:

pep8speaks commented Jul 26, 2021 • edited

Comment last updated at 2021-08-17 10:47:06 UTC

codecov bot commented Jul 26, 2021 • edited

Codecov Report

skoudoro commented Aug 2, 2021

drombas commented Aug 2, 2021

arokem left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BramshQamar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BramshQamar Aug 17, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iPsych commented Nov 16, 2022

drombas commented Nov 16, 2022

drombas commented Jul 26, 2021 •

edited

pep8speaks commented Jul 26, 2021 •

edited

codecov bot commented Jul 26, 2021 •

edited

BramshQamar Aug 17, 2021 •

edited