scilus / scilpy

The Sherbrooke Connectivity Imaging Lab (SCIL) Python dMRI processing toolbox

Scilpy

Scilpy is the main library supporting research and development at the Sherbrooke Connectivity Imaging Lab (SCIL).

Scilpy mainly comprises tools and utilities to quickly work with diffusion MRI. Most of the tools are based on or are wrappers of the DIPY library, and most of them will eventually be migrated to DIPY. Those tools implement the recommended workflows and parameters used in the lab.

The library is now built for Python 3.10 so be sure to create a virtual environnement for Python 3.10. If this version is not installed on your computer:

sudo add-apt-repository ppa:deadsnakes/ppa
sudo apt-get install python3.10 python3.10-dev python3.10-venv python3.10-minimal python3.10-tk

Make sure your pip is up-to-date before trying to install:

pip install --upgrade pip

The library's structure is mostly aligned on that of DIPY.

The library and scripts can be installed locally by using:

pip install -e .

If you don't want to install legacy scripts:

export SCILPY_LEGACY='False'
pip install -e .

(Then, without the legacy scripts, if you want to use pytest, use:)

pytest --ignore=scripts/legacy

On Linux, most likely you will have to install libraries for COMMIT/AMICO

sudo apt install libblas-dev liblapack-dev

While on MacOS you will have to use (most likely)

brew install openblas lapack

On Ubuntu >=20.04, you will have to install libraries for matplotlib

sudo apt install libfreetype6-dev

Note that using this technique will make it harder to remove the scripts when changing versions. We highly recommend working in a Python Virtual Environment.

Scilpy documentation is available: https://scilpy.readthedocs.io/en/latest/

scilpy's Issues

Generate documentation for scilpy scripts

#166 Is pretty great to find scripts. However, a more robust and user friendly option would be to generate documentation automatically from scripts.

Save seeds option take a lot of memory

Save seeds option in PFT and local tracking take a lot of memory comparing to the streamlines.

DTI tensor, coeff order, saved for "fiberNav"

In the "scil_compute_dti_metrics.py"
the script use the fiber Nav order for coef,
should we simply use the default dipy (lower triangle order)

        # Get the Tensor values and format them for visualisation
        # in the Fibernavigator.
        tensor_vals = lower_triangular(tenfit.quadratic_form)
        correct_order = [0, 1, 3, 2, 4, 5]
        tensor_vals_reordered = tensor_vals[..., correct_order]
        fiber_tensors = nib.Nifti1Image(
            tensor_vals_reordered.astype(np.float32), affine)
        nib.save(fiber_tensors, args.tensor)```

Rename assert_outputs_exists

Rename assert_outputs_exists -> assert_outputs_exist

visualize connectivity enhancement

If possible, it would be nice to merge visu options from @StongeEtienne in the what @frheault wrote. @StongeEtienne visualize connectivity matrix is called npy2png.
https://www.dropbox.com/s/uvgafmkyxo1ud37/npy2png.py?dl=0

Should not be too hard to implement. @StongeEtienne can you do it when you have a bit of time?

Thanks!

Connectivity: save as json

Add possibility to save as a .json file, for easier manipulation.

commit test fails

    mit.build_operator(build_dir=tmp_dir.name)
TypeError: build_operator() takes no keyword arguments

remove multiple_shell_energy.py redundancy with dipy code

Once https://github.com/nipy/dipy/pull/1711/files is merged in Dipy, we should use Dipy's version to have less code redundancy.

Originally posted by @jchoude in #32

pairwise bundle adjacency of two bundle always higher than one

Hi,

I have been checking the bundle adjacency of two bundles from test-retest, no matter how the two bundles are close to each other, none of the values fall between 0 and 1, it's always equal to or higher to 1.
Here is a result
https://drive.google.com/file/d/1VZL5w8Y66dn21P72n_sSeRLQd-1k8EKp/view?usp=sharing
the minimum value is exactly 1.0, which makes me wonder if there is a systematic reason that the optimal value of bundle_adjacency_voxel is one? Although I tried to compute two identical bundle and it gave me a value of zero.

FYI, when I'm computing the bundle_adjacency_voxel, the reference I passed to it is a T1 image with 1mm isotropic voxel size.

Maybe we should detect this with a heuristic or send a warning if 100+ images given?

Maybe we should detect this with a heuristic or send a warning if 100+ images given?

Originally posted by @mdesco in #233

Update VTK for python 3.8

Maybe we can set in the requirements.txt vtk >= 8.1.2 ?

@arnaudbore @jchoude What do you think ?

Update scipy for python 3.7

scipy==1.0.1 is not enough for python 3.7
Which version should we use ?
What needs to be modified ?

scil_decompose_connectivity problem when saving all streamlines

If I run the decompose command with the all the saving options, it works fine.

However, if I run this line, with or without invalid coordinates (ic) removed (there are 0 ic streamlines removed), I get the following error.

scil_decompose_connectivity.py Tracking_SET/final_tracks_noloops_ic.trk brainetomme_int16_2.nii.gz connectome_brainetomme --save_raw_connections --save_intermediate --save_discarded -f
/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py:208: DeprecationWarning: get_data() is deprecated in favor of get_fdata(), which has a more predictable return type. To obtain get_data() behavior going forward, use numpy.asanyarray(img.dataobj).

• deprecated from version: 3.0

• Will raise <class 'nibabel.deprecator.ExpiredDeprecationError'> as of version: 5.0
  data_labels = img_labels.get_data()
  /Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py:245: DeprecationWarning: 'ArraySequence.data' property is deprecated.
  Please use the 'ArraySequence.get_data()' method instead

• deprecated from version: 3.0

• Will raise <class 'nibabel.deprecator.ExpiredDeprecationError'> as of version: 4.0
  indices, points_to_idx = uncompress(sft.streamlines, return_mapping=True)
  /Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py:255: DeprecationWarning: get_data() is deprecated in favor of get_fdata(), which has a more predictable return type. To obtain get_data() behavior going forward, use numpy.asanyarray(img.dataobj).

• deprecated from version: 3.0

• Will raise <class 'nibabel.deprecator.ExpiredDeprecationError'> as of version: 5.0
  img_labels.get_data(), real_labels,
  /Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/scilpy/tractanalysis/tools.py:35: RuntimeWarning: invalid value encountered in true_divide
  ray = ray / np.linalg.norm(ray)
  /Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/dipy/io/stateful_tractogram.py:389: RuntimeWarning: invalid value encountered in less
  if np.any(bbox_corners < 0):
  2020-03-29 18:53:09 DINF-DESCOT-04J StatefulTractogram[36423] ERROR Voxel space values lower than 0.0
  2020-03-29 18:53:09 DINF-DESCOT-04J StatefulTractogram[36423] ERROR Voxel space values higher than dimensions
  Traceback (most recent call last):
  File "/Users/desm2239/my_env/scilpy_public/bin/scil_decompose_connectivity.py", line 4, in
  import('pkg_resources').run_script('scilpy==0.1.dev0', 'scil_decompose_connectivity.py')
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/pkg_resources/init.py", line 658, in run_script
  self.require(requires)[0].run_script(script_name, ns)
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/pkg_resources/init.py", line 1438, in run_script
  exec(code, namespace, namespace)
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py", line 386, in
  main()
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py", line 304, in main
  'raw', in_label, out_label)
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/scilpy-0.1.dev0-py3.7-macosx-10.14-x86_64.egg/EGG-INFO/scripts/scil_decompose_connectivity.py", line 105, in _save_if_needed
  save_tractogram(sft, out_name)
  File "/Users/desm2239/my_env/scilpy_public/lib/python3.7/site-packages/dipy/io/streamline.py", line 42, in save_tractogram
  raise ValueError('Bounding box is not valid in voxel space, cannot '
  ValueError: Bounding box is not valid in voxel space, cannot load a valid file if some coordinates are invalid. Please use the function remove_invalid_streamlines to discard invalid streamlines or set bbox_valid_check to False

Many connections are written and the last one is: 155_231.trk

Data is here if @frheault can test and correct problem.
https://www.dropbox.com/s/4m273c7y1bixj83/decompose_prob.zip?dl=0

Thanks!

assert_input_exist and assert_output_exist missing list

scilpy/scripts/scil_get_subset_streamlines.py

Line 37 in 151a2ff

assert_outputs_exists(parser, args, args.out_tractogram)

These functions always required a list, when a string is given the check is useless

Function gen_scheme has changed name

In the script scil_generate_smapling_scheme.py, there is an imported function whose name was not changed correctly.

scilpy/scripts/scil_generate_sampling_scheme.py

Line 10 in 47c8693

from scilpy.samplingscheme.gen_scheme import gen_scheme

In the module scilpy.samplingscheme.gen_scheme, the function is named generate_scheme

scil_reorder_connectivity better help

1. --labels_list LABELS_LIST
   List saved by the decomposition script,
   the json must contain labels rather than coordinates.

Please say this is a txt file and add (.txt) at the end.

2. Also, not a fan of --prefix why not just fix the output. Is it because the input can infact be *npy ?

3. would be cool that input is *npy because if I start working on a sub-connectome of the full connectome, I want to reorder my vol.npy, len.npy, fa.npy, track-count.noy, etc...

Thanks @frheault

--inverse sanity check in RBX script

Suggestion: maybe check the --inverse option before doing the streamline computation

Especially important in the single bundle version of the script

docstring look of scil_filter_connectivity.py

@frheault can you fix the first paragraph of scil_filter_connectivity docstring please. It is unreadable. I just see 1 big block of text.

Script to facilitate filtering of connectivity matrices. The same could be
achieved through a complex sequence of scil_connectivity_math.py. Can be used
with any connectivity matrix from scil_compute_connectivity.py. For example, a
simple filtering (Jasmeen style) would be: scil_filter_connectivity.py
out_mask.npy --greater_than */sc.npy 1 0.90 --lower_than */sim.npy 2 0.90
--greater_than */len.npy 40 0.90 -v; This will result in a binary mask where
each node with a value of 1 represents a node with at least 90% of the
population having at least 1 streamline, 90% of the population is similar to
the average (2mm) and 90% of the population having at least 40mm of average
streamlines length. --greater_than or --lower_than expect the same convention:
MATRICES_LIST VALUE_THR POPULATION_PERC It is strongly recommended (but not
enforced) that the same number of connectivity matrices is used for each
condition. This script performs an intersection of all conditions, meaning
that all conditions must be met in order not to be filtered. If the user wants
to manually handle the requirements, --keep_condition_count can be used and
manually binarized using scil_connectivity_math.py

thanks

Small typo

When running

scil_recognize_multi_bundles.py :

WARNING Model 486759/SLF_2_L_m.trk' simplified at threshod 3mm with 2783 centroids

Missing a "l" in threshold :) @frheault

Add script to print infos from tractogram/image

scil_tractogram_info.py ou scil_tractogram_header.py
scil_image_info.py ou scil_image_header.py

ou generic: scil_info.py or scil_header.py

Connectivity: manage self-connections

Currently, in the scil_compute_connectivity script, we allow for self connection. This could be parametrized.

scil_image_math.py can not concatenate 4D and 3D images

scil_compute_connectivity output extension

Please add to the help the extensions expected for --volume, --streamline_count, etc.....

@frheault

scil_visualize_seed --save option

@AntoineTheb is the behavior like we really want?
Right now, for example, I'm looping over 30 bundles and I scil_visualize_seed.py track.trk --save track.png to have the 30 PNG files.

But I need to close the VTK window 30 times. Is this wanted? Probably not?

Also, how hard to add the seed option to scil_visualize_bundle_mosaic?

Thanks!

Create Doc string for filter tractogram

Create doc string for the functions added in PR #24:
read_info_from_mb_bdo
target_line_based
filter_grid_roi
pre_filtering_for_geometrical_shape
filter_ellipsoid
filter_cuboid

Change call to add_verbose_arg

Suggestion: scil_image_math.py equivalent of mrcat

@frheault I'm only using mrtrix these days to run mrcat stuff.
Possible to include in scil_image_math?

Thanks!

Values loop around datatype when int16 DWI

Loop around the datatype in extract B0.

Moreover, when you generate the residual in dti metrics, the script crash because create a copy in float 64.

Must be patch and tested before release.

scil_recognize_multi_bundles.py crashes if --output does not exist

Could be fixed and maybe use --output_dir to be uniform with other scripts.
@frheault, easy one to fix.

Thanks!

scil_compute_seed_density_map.py bug

@AntoineTheb pense-tu que tu peux tester/checker ce bug.
Les bundles sortent de recobundle scilpy de Francois. Le bundle en question (AF) n'a aucune invalid coordinates. Il visualise bien avec tous les outils...

Data: https://www.dropbox.com/s/b443bzf5z9bmyvj/bug_antoine.zip?dl=0

Merci

scil_compute_seed_density_map.py AF_L_m.trk yo.nii.gz
Traceback (most recent call last):
File "/Users/desm2239/my_env/scilpy_public/bin/scil_compute_seed_density_map.py", line 4, in
import('pkg_resources').run_script('scilpy==0.1.dev0', 'scil_compute_seed_density_map.py')
File "/Users/desm2239/my_env/scilpy_public/lib/python3.6/site-packages/pkg_resources/init.py", line 658, in run_script
self.require(requires)[0].run_script(script_name, ns)
File "/Users/desm2239/my_env/scilpy_public/lib/python3.6/site-packages/pkg_resources/init.py", line 1438, in run_script
exec(code, namespace, namespace)
File "/Users/desm2239/my_env/scilpy_public/lib/python3.6/site-packages/scilpy-0.1.dev0-py3.6-macosx-10.9-x86_64.egg/EGG-INFO/scripts/scil_compute_seed_density_map.py", line 97, in
main()
File "/Users/desm2239/my_env/scilpy_public/lib/python3.6/site-packages/scilpy-0.1.dev0-py3.6-macosx-10.9-x86_64.egg/EGG-INFO/scripts/scil_compute_seed_density_map.py", line 88, in main
seed_density[tuple(seed_voxel)] += 1
IndexError: index 161 is out of bounds for axis 1 with size 144

argparse guideline

in_ for mandatory inputs
out_ for mandatory outputs

_build_arg_parser

etc...

Tests versus documentation

@mdesco about the data and example you mentioned earlier. I think we would need to fetch data like Dipy and do an example that can be run. Maybe every time there is a merge these examples could be run and if they crash we have a notification to modify them.

I say that because the problems with separate example is that after a month some will stop working (fixing parameters name or mandatory parameters order). It would act as tutorials as well as very high level 'tests' that do not actually verify output. But since we could run them in sequence the output of each is kind of tested as an input for the next.

I would prefer that to documentation or example since these 'tests' would force us to keep them up to date. They could be run on-demand or after every merge. The only problem with that idea is the size of the data required and the execution time (would be a couple Gb and likely that 2-3 hours to run)

@GuillaumeTh What do you think? Do you think it is possible, to download some data (wget from a lab cloud), run a bash script with a singularity (build using the master, and ants/mtrix/fsl I guess) and if there is no crash the 'test' pass. With basic comments, it would also be an example.

Originally posted by @frheault in #193 (comment)

A letter for the code.

Dear Sir,
I am reading your paper on SET.I want to reproduce the process in your paper, but the scilpy toolkit can't be downloaded for the time being. When will it be convenient for you to provide the code about scilpy?
I am looking forward to your reply.

scil_compute_dti_metrics takes forever

Need to deal with ext in scil_visualise multishell

We approved PR#127, but you still had a discussion going on about how to deal with the ext (with or without the point, how to deal with bval, bvecs).

Just opening the issue to help you remember.

check scil_extract_b0.py

No assert output exist

Use of Deprecated function

ren = window.Renderer()

In scil_visualize_bundles_mosaic

pycharm suggests: fury.windows.Scene instead

doc/tractogram_registration.md

Lets update the terms to correspond to the new scripts --help
@frheault merci!

label_list.txt info

We need an example or a description somewhere of what a label_list looks like for the scil_compute_connectivity.py.

The PRs had very rich info on example data and example scripts calling the functions which are now lost since the PR are merged. It would be amazing to find a way to bring into the doc or elsewhere...

I had to go back in old PRs to find an example how to call scil_compute_connectivity with examples as inputs.

@frheault any idea how to easy this? Otherwise, your mailbox and slack will explode

PFT script default params

I think we should change the default PFT script params, to the original in Gabriel Girard Paper
step size of 0.2
and @gabknight might want to check for others in "scil_compute_pft.py"
https://github.com/scilus/scilpy/blob/master/scripts/scil_compute_pft.py

Change the python versions in scripts

In the case that you don't have any virtualenv, we can not call the script because the script requiert python and not python3. We must to replace python by python3 in scripts.

scil_recognize_multi_bundles.py: Invalid example

The doc for scil_recognize_multi_bundles.py mentions

Example data and usage available at: https://zenodo.org/deposit/3613688

However, the link seems to be dead or something changed because when I try to access it, I get

Is it supposed to be restricted ? Or did I miss that I must ask for access ?

scil_compute_connectivity allclose 'affine' too strict

This command breaks:
scil_compute_connectivity.py decompose_sub-1144_ses-2.h5 atlas/sub-1144_ses-2/brainnetome.nii.gz --volume volume.npy --streamline_count track-count.npy --length length.npy --metrics *fa.nii.gz fa.npy --density_weighting --no_self_connection -f

On this data:
https://www.dropbox.com/s/ursgtsbn26iakx1/conn.zip?dl=0

As recommanded by @frheault, if I comment the allclose 'affine' lines of the code, things run and seem to be fine.

I don't really mind, but based on François's previous question about trying to have the same nomenclature everywhere:

I don't really mind, but based on François's previous question about trying to have the same nomenclature everywhere:

Everywhere else, it is called _mean_std, not _stats. I think I actually prefer stats but then maybe we should change it everywhere else?

Originally posted by @EmmaRenauld in #202

Deletes files with no warning

If there are files in --output_dir of scil_recognize_multi_bundles, they are deleted with no warnings!

scil_image_math.py float16 datatype

float16 is a typo, right? float32 or float64?
Anyways, it crashes if float16 is used.

@frheault

Should assert_outputs_exist() test directory by default

in assert_outputs_exist(...)
should "check_dir_exist" be True by default ?

dtype mandatory for getfdata

Need to fix this everywhere.

VTK not found for Python 3.7.3

pip install -r requirements.txt is not working under Python 3.7.3 (fresh venv)

Collecting vtk==8.1.0 (from -r requirements.txt (line 18))
  Could not find a version that satisfies the requirement vtk==8.1.0 (from -r requirements.txt (line 18)) (from versions: 8.1.2)
No matching distribution found for vtk==8.1.0 (from -r requirements.txt (line 18))

mac vtk-9.1.0 - black window - fury

catalina 10.15.6 with python 3.6

Black window - fury

scil_clean_qbx_clusters.py

2020-11-03 11:29:14.836 (  62.176s) [          7DF26A]vtkOpenGLPolyDataMapper:271   WARN| vtkOpenGLPolyDataMapper::GetVertexShaderCode was deprecated for VTK 9.0 and will be removed in a future version.  Use vtkOpenGLShaderProperty::GetVertexShaderCode instead.

2020-11-03 11:29:14.836 (  62.176s) [          7DF26A]vtkOpenGLPolyDataMapper:263   WARN| vtkOpenGLPolyDataMapper::SetVertexShaderCode was deprecated for VTK 9.0 and will be removed in a future version.  Use vtkOpenGLShaderProperty::SetVertexShaderCode instead.

Quick fix:
-> force vtk-8.1.2