Skip to content

Commit

Permalink
Begin work on DEMO
Browse files Browse the repository at this point in the history
  • Loading branch information
VincentBeaud committed Oct 2, 2024
1 parent db99c84 commit a7c9e94
Show file tree
Hide file tree
Showing 13 changed files with 108 additions and 76 deletions.
86 changes: 44 additions & 42 deletions docs/fibertube/DEMO.md
Original file line number Diff line number Diff line change
@@ -1,15 +1,19 @@
# Demo Workshop
In this demo, you will be introduced to the main scripts of this project as you apply them on simple data.
<br><br>
Our main objective is better understand and quantify the fundamental limitations of tractography algorithms, and how they might evolve as we introduce histology as *a-priori* information for fiber orientation. To do so, we will be evaluating tractography's ability to reconstruct individual white matter fiber strands at various simulated extreme resolutions.
Our main objective is better understand and quantify the fundamental limitations of tractography algorithms, and how they might evolve as we approach microscopy resolution where individual axons can be seen, tracked or segmented. To do so, we will be evaluating tractography's ability to reconstruct individual white matter fiber strands at various simulated extreme resolutions.
## Terminology
Here is a list of special terms and definitions used in this project:
- Axon: Physical object. Portion of the nerve cell that carries out the electrical impulse to other neurons.
- Streamline: Virtual object. Series of coordinates propagated through a stream of directional data using a tractography algorithm.
- Fibertube: Virtual representation of an axon. It is composed of a centerline and a diameter.
- Fibertube Tractography: The application of a tractography algorithm directly on fibertubes (to reconstruct them), without using a discretized grid of fODFs or peaks.
Here is a list of terms and definitions used in this project.

In the context of this project, centerlines of fibertubes will come on the form of a tractogram, and we will provide each of them with an artificial axonal diameter coming from a text file. The resulting set of fibertubes will sometimes be referred to as "ground-truth data".
General:
- Axon: Bio-physical object. Portion of the nerve cell that carries out the electrical impulse to other neurons. (On the order of 0.1 to 1um)
- Streamline: Virtual object. Series of equidistant 3D coordinates approximating an underlying fiberous structure.

Fibertube Tracking:
- Centerline: Virtual object. Series of equidistant 3D coordinates representing the directional information of a fibertube.
- Fibertube: Virtual representation of an axon. It is composed of a centerline and a single diameter for its whole length.
- Fibertube segment: Because centerlines are made of discrete coordinates, fibertubes end up being composed of a series of equally lengthed adjacent cylinders. A fibertube segment is any single one of those cylinders.
- Fibertube Tractography: The application of a tractography algorithm directly on fibertubes to reconstruct them. Contrary to traditional white matter fiber tractography, fibertube tractography does not rely on a discretized grid of fODFs or peaks.

![Fibertube visualized in Blender](https://github.com/VincentBeaud/fibertube_tracking/assets/77688542/25494d10-a8d5-46fa-93d9-0072287d0105)

Expand All @@ -21,102 +25,100 @@ This project can be split into 3 major steps:
- Tracking and experimentation <br>
We will perform 'Fibertube Tracking' on our newly formed set of fibertubes with a variety of parameter combinations.
- Calculation of metrics <br>
By passing the resulting tractogram through different analytic pipelines and scripts, we will acquire connectivity and fiber reconstruction scores for each of the parameter combinations.

Of course, some new questions will arise from results and we may branch out of this 3-step frame temporarily.
By passing the resulting tractogram through different analytic pipelines and scripts (like Tractometer), we will acquire connectivity and fiber reconstruction scores for each of the parameter combinations.

## Preparing the data
> [!IMPORTANT]
> All commands written down below assume that your console is positioned in the `demo/` folder.
> All commands written down below assume that your console is positioned in the folder containing your data.
The data required to perform fibertube tractography comes in two files:
- `./disco_centerlines.tck` is a 2000 streamlines subset of the DISCO dataset.
- `./disco_diameters.txt` contains the diameters.
- `./disco_mask.nii.gz` is a white matter mask for spatial reference.
- `./centerlines.trk` contains the entire ground-truth of the DISCO dataset.
- `./diameters.txt` contains the diameters.

![DISCO subset visualized in MI-Brain](https://github.com/VincentBeaud/fibertube_tracking/assets/77688542/197b3f1f-2f57-41d0-af0a-5f7377bab274)

The first thing to do is resample `disco_centerlines.trk` so that each centerline is formed of
The first thing to do is resample `centerlines.trk` so that each centerline is formed of
segments no longer than 0.2 mm.

> [!NOTE]
> This is because the next script will rely on a KDTree to find all neighboring fibertube segments of any given point. Because the search radius is set at the length of the longest fibertube segment, the performance drops significantly if they are not shortened to ~0.2mm.
To resample a tractogram, we can use this script from scilpy:
```
scil_tractogram_resample_nb_points.py disco_centerlines.tck disco_centerlines_resampled.tck --step_size 0.2 --reference disco_mask.nii.gz
scil_tractogram_resample_nb_points.py centerlines.trk centerlines_resampled.trk --step_size 0.2
```

Next, we want to filter out intersecting fibertubes, to make the data anatomically plausible and remove any partial volume effect. This step is crucial to ensure perfect fiber reconstruction at lower scale.
Next, we want to filter out intersecting fibertubes, to make the data anatomically plausible and remove any partial volume effect.

![Fibertube intersection visualized in Blender](https://github.com/VincentBeaud/perfect_tracking/assets/77688542/ede5d949-d7a5-4619-b75b-72fd41d65b38)

This is accomplished using `ft_filter_collisions.py`. <br>
For this demo, let's go all-in and turn on every option.
This is accomplished using `scil_tractogram_filter_collisions.py`. <br>

```
python ../ft_filter_collisions.py disco_centerlines_resampled.tck disco_diameters.txt disco_centerlines_clean.tck -cr -cd -v -f --reference disco_mask.nii.gz
scil_tractogram_filter_collisions.py centerlines_resampled.trk diameters.txt fibertubes.trk --save_colliding --out_metrics metrics.txt -v
```
> [!IMPORTANT]
> Because this is a script from the project (and not in our environment), we need to call it with "python".

After a short wait, you should get something like:
```
...
├── disco_centerlines_clean_obstacle.tck
├── disco_centerlines_clean_invalid.tck
├── disco_centerlines_clean_diameters.tck
├── disco_centerlines_clean.tck
├── centerlines_resampled_obstacle.trk
├── centerlines_resampled_invalid.trk
├── fibertubes.trk
...
```

## Visualising collisions
Throughout the entire flow, you will be provided with scripts to visualize the data at different steps.
As you may have guessed from the output name, this script automatically combines the diameter to the centerlines as data_per_streamline in the output tractogram. This is why we named it "fibertubes.trk".

## Visualising collisions
By calling:
```
python ../ft_visualize_collisions.py disco_centerlines_clean_invalid.trk --obstacle disco_centerlines_clean_obstacle.trk --ref_tractogram disco_centerlines.tck --reference disco_mask.nii.gz
scil_viz_tractogram_collisions.py centerlines_resampled_invalid.trk --obstacle centerlines_resampled_obstacle.trk --ref_tractogram centerlines.trk
```
You are able to see exactly which streamline has been filtered ("invalid" - In red) as well as the streamlines they collided with ("obstacle" - In green).
In white and lower opacity is the original tractogram passed as `--ref_tractogram`.

![Filtered intersections visualized in 3D](https://github.com/VincentBeaud/fibertube_tracking/assets/77688542/4bc75029-0d43-4664-8502-fd528e9d93f4)

### Fibertube metrics
Before we get into tracking. Here is an overview of the metrics that can be computed from the clean data, using `ft_fibers_metrics.py`.
Before we get into tracking. Here is an overview of the metrics that we saved in `metrics.txt`:

- `min_external_distance`: Smallest distance separating two fibers in the entire set.
- `max_voxel_anisotropic`: Diagonal vector of the largest possible anisotropic voxel that would not intersect two fibers.
- `min_external_distance`: Smallest distance separating two fibertubes, outside their diameter.
- `max_voxel_anisotropic`: Diagonal vector of the largest possible anisotropic voxel that would not intersect two fibertubes.
- `max_voxel_isotropic`: Isotropic version of max_voxel_anisotropic made by using the smallest component. <br>
Ex: max_voxel_anisotropic: (3, 5, 5) => max_voxel_isotropic: (3, 3, 3)
- `max_voxel_rotated`: Rotated version of max_voxel_anisotropic to align it will (1, 1, 1). This makes it an isotropic voxel, but is only valid if the entire tractogram is rotated the same way.
- `max_voxel_rotated`: Largest possible isotropic voxel obtainable if the tractogram is rotated. It is only usable if the entire tractogram is rotated according to [rotation_matrix].
Ex: max_voxel_anisotropic: (1, 0, 0) => max_voxel_rotated: (0.5774, 0.5774, 0.5774)
- `rotation_matrix`: 4D transformation matrix representing the rotation to be applied on [in_centerlines] for transforming `max_voxel_anisotropic` into `max_voxel_rotated` (see scil_tractogram_apply_transform.py)
- `rotation_matrix`: 4D transformation matrix representing the rotation to be applied on the tractogram to align max_voxel_rotated with the coordinate system (see scil_tractogram_apply_transform.py).

![Metrics (without max_voxel_rotated) visualized in Blender](https://github.com/VincentBeaud/perfect_tracking/assets/77688542/95cd4e50-1a36-49af-ac11-0d5f33d3f32e)
<br>
![max_voxel_rotated visualized in Blender](https://github.com/VincentBeaud/perfect_tracking/assets/77688542/72812e47-371f-4005-b289-1de0d70d2f33)

> [!NOTE]
> This information can be useful for analyzing the reconstruction obtained through tracking, as well as for performing track density imaging. The latter will however require a more aggressive fibertube filtering using
> the `--min_distance` argument in the filtering script.
> This information can be useful for analyzing the reconstruction obtained through tracking, as well as for performing track density imaging.
## Performing fibertube tracking
We're finally at the tracking phase! Using the script `ft_fibertube_tracking.py`, you are able to track without relying on a discretized grid of directions. Instead, you will be propagating a streamline through fibertubes and degrading the resolution during the process by using a sphere of "blur".
We're finally at the tracking phase! Using the script `scil_fibertube_tracking.py`, you are able to track without relying on a discretized grid of directions or fODFs. Instead, you will be propagating a streamline through fibertubes and degrading the resolution by using a `blur_radius`. The way it works is as follows:

### Tracking
When the tracking algorithm is about to select a new direction to propagate the current streamline, it will build a sphere of radius `blur_radius` and pick randomely from all the fibertube segments intersecting with it. The larger the intersection volume, the more likely a fibertube segment is to be picked and used as a tracking direction. This makes fibertube tracking inherently probabilistic.
Theoretically, with a `blur_radius` of 0, any given set of coordinates has either a single tracking direction because it is within a fibertube, or no direction at all from being out of one. In fact, this behavior won't change until the diameter of the sphere is larger than the smallest distance separating two fibertubes. When this happens, more than one fibertubes will intersect the `blur_radius` sphere and introduce partial volume effect.

With a sphere of radius 0, any given set of coordinates has either a single tracking direction because it is within a fibertube, or no direction at all from being out of one. In fact, this behavior won't change until the diameter of the sphere is larger than the smallest distance separating two fibers.

### Seeding
For now, a number of seeds is set randomly within the first segment of every fibertube. We can however change how many fibers will be tracked, as well as the amount of seeds within each. (See Seeding options in the help menu).

The interface of the script is very similar to `scil_tracking_local_dev.py`, but simplified and with a `sampling_radius` mandatory option. This will be the radius of our "blurring" sphere. Let us do:
<br>
The interface of the script is very similar to `scil_tracking_local_dev.py`, but simplified and with a `blur_radius` option. Let us do:

```
python ../ft_tracking.py disco_centerlines_clean.tck disco_centerlines_clean_diameters.txt disco_mask.nii.gz reconstruction.tck 0.01 0.01 -v -f --nb_seeds_per_fiber 2 --nb_fibers 2 --save_seeds --save_config --processes 4 --reference disco_mask.nii.gz
scil_fibertube_tracking.py fibertubes.trk tracking.trk 0.01 0.01 --nb_fibertubes 3 --out_config tracking_config.txt --processes 4 -v -f
```
This should take a few minutes at most. However, if you don't mind waiting a little bit, feel free to play with the parameters and explore the resulting tractogram.

> [!NOTE]
> Given the time required for each streamline, the `--processes` parameter will become a good friend of yours.
> Given the time required for each streamline, the `--processes` parameter will be very useful.
HERE
## Visualizing fibertube coverage of tracked streamlines
Another useful script is:
```
Expand Down
Binary file added docs/fibertube/centerlines_resampled_invalid.trk
Binary file not shown.
Binary file not shown.
Binary file added docs/fibertube/fibertubes.trk
Binary file not shown.
20 changes: 20 additions & 0 deletions docs/fibertube/metrics.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
{
"min_external_distance": [
1.2642728071303542e-06
],
"max_voxel_anisotropic": [
9.599388306469336e-08,
9.921478749674861e-07,
7.776976076456776e-07
],
"max_voxel_isotropic": [
9.599388306469336e-08,
9.599388306469336e-08,
9.599388306469336e-08
],
"max_voxel_rotated": [
7.299282455258108e-07,
7.299282455258108e-07,
7.299282455258108e-07
]
}
Binary file added docs/fibertube/tracking.trk
Binary file not shown.
6 changes: 6 additions & 0 deletions docs/fibertube/tracking_config.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
{
"step_size": 0.01,
"blur_radius": 0.01,
"nb_fibertubes": 3,
"nb_seeds_per_fibertube": 5
}
10 changes: 5 additions & 5 deletions scilpy/tracking/seed.py
Original file line number Diff line number Diff line change
Expand Up @@ -273,22 +273,22 @@ class FibertubeSeedGenerator(SeedGenerator):
fibertube tracking. Generates a given number of seed within the first
segment of a given number of fibertubes.
"""
def __init__(self, centerlines, diameters, nb_seeds_per_fiber):
def __init__(self, centerlines, diameters, nb_seeds_per_fibertube):
"""
Parameters
----------
centerlines: list
Tractogram containing the fibertube centerlines
diameters: list
Diameters of each fibertube
nb_seeds_per_fiber: int
nb_seeds_per_fibertube: int
"""
self.space = Space.VOXMM
self.origin = Origin.NIFTI

self.centerlines = centerlines
self.diameters = diameters
self.nb_seeds_per_fiber = nb_seeds_per_fiber
self.nb_seeds_per_fibertube = nb_seeds_per_fibertube

def init_generator(self, rng_seed, numbers_to_skip):
"""
Expand Down Expand Up @@ -344,7 +344,7 @@ def init_generator(self, rng_seed, numbers_to_skip):
def get_next_pos(self, random_generator: np.random.Generator,
shuffled_indices, which_seed):

which_fi = which_seed // self.nb_seeds_per_fiber
which_fi = which_seed // self.nb_seeds_per_fibertube

fiber = self.centerlines[shuffled_indices[which_fi]]
radius = self.diameters[shuffled_indices[which_fi]] / 2
Expand All @@ -357,7 +357,7 @@ def get_next_pos(self, random_generator: np.random.Generator,
def get_next_n_pos(self, random_generator, shuffled_indices,
which_seed_start, n):

which_fi = which_seed_start // self.nb_seeds_per_fiber
which_fi = which_seed_start // self.nb_seeds_per_fibertube

fiber = self.centerlines[shuffled_indices[which_fi]]
radius = self.diameters[shuffled_indices[which_fi]] / 2
Expand Down
28 changes: 16 additions & 12 deletions scripts/scil_fibertube_score_tractogram.py
Original file line number Diff line number Diff line change
Expand Up @@ -72,15 +72,19 @@ def _build_arg_parser():
formatter_class=argparse.RawTextHelpFormatter)

p.add_argument('in_fibertubes',
help='Path to the tractogram file containing the \n'
'fibertubes with their respective diameter saved \n'
'data_per_streamline (must be .trk). \n'
'The fibertubes must be void of any collision \n'
'(see scil_filter_intersections.py). \n')

p.add_argument('in_tractogram',
help='Path to a text file containing the ground-truth \n'
'reconstruction (must be .trk) with seeds saved as \n'
help='Path to the tractogram file (must be .trk) \n'
'containing ground-truth fibertubes. They must be: \n'
'1- Void of any collision. \n'
'2- With their respective diameter saved \n'
'as data_per_streamline. \n'
'For both of these requirements, see \n'
'scil_tractogram_filter_collisions.')

p.add_argument('in_tracking',
help='Path to the tractogram file (must be .trk) \n'
'containing the reconstruction of ground-truth \n'
'fibertubes made from fibertube tracking. Seeds \n'
'used for tracking must be saved as \n'
'data_per_streamline.')

p.add_argument('in_config',
Expand Down Expand Up @@ -127,7 +131,7 @@ def main():
.format(args.out_metrics))

assert_inputs_exist(parser, [args.in_fibertubes, args.in_config,
args.in_tractogram])
args.in_tracking])
assert_outputs_exist(parser, args, [args.out_metrics])

our_space = Space.VOXMM
Expand All @@ -147,7 +151,7 @@ def main():
len(centerlines))

logging.debug('Loading reconstructed tractogram')
in_sft = load_tractogram(args.in_tractogram, 'same', our_space,
in_sft = load_tractogram(args.in_tracking, 'same', our_space,
our_origin)
streamlines = in_sft.get_streamlines_copy()
streamlines, streamlines_length = get_streamlines_as_fixed_array(
Expand All @@ -156,7 +160,7 @@ def main():
logging.debug("Loading seeds")
if "seeds" not in in_sft.data_per_streamline:
parser.error('No seeds found as data per streamline on ' +
args.in_tractogram)
args.in_tracking)

seeds = in_sft.data_per_streamline['seeds']
seeds_fiber = resolve_origin_seeding(seeds, centerlines, diameters)
Expand Down
14 changes: 9 additions & 5 deletions scripts/scil_fibertube_tracking.py
Original file line number Diff line number Diff line change
Expand Up @@ -168,9 +168,8 @@ def main():
parser = _build_arg_parser()
args = parser.parse_args()

if args.verbose:
logging.getLogger().setLevel(logging.DEBUG)
logging.getLogger('numba').setLevel(logging.WARNING)
logging.getLogger().setLevel(logging.getLevelName(args.verbose))
logging.getLogger('numba').setLevel(logging.WARNING)

if not nib.streamlines.is_supported(args.in_fibertubes):
parser.error('Invalid input streamline file format (must be trk ' +
Expand All @@ -180,11 +179,16 @@ def main():
parser.error('Invalid output streamline file format (must be trk ' +
'or tck): {0}'.format(args.out_tractogram))

out_tractogram_no_ext, ext = os.path.splitext(args.out_tractogram)
_, out_config_ext = os.path.splitext(args.out_config)
out_tractogram_no_ext, out_tractogram_ext = os.path.splitext(args.out_tractogram)

if out_config_ext != '.txt':
parser.error('Invalid output file format (must be txt): {0}'
.format(args.out_config))

outputs = [args.out_tractogram]
if not args.do_not_save_seeds:
outputs.append(out_tractogram_no_ext + '_seeds' + ext)
outputs.append(out_tractogram_no_ext + '_seeds' + out_tractogram_ext)

assert_inputs_exist(parser, [args.in_fibertubes])
assert_outputs_exist(parser, args, outputs, [args.out_config])
Expand Down
Loading

0 comments on commit a7c9e94

Please sign in to comment.