DTI Metrics + FRF Error

Hi, tractoflow is running into an error when calculating the FRF and DTI metrics. Below is the command line along with the .command.err logs from both processes.

nextflow run /path/to/main.nf --input /path/to/data/ -with-singularity /path/to/scilus_1.5.0.sif -profile fully_reproducible -resume -with-report report.html --run_pft_tracking false --run_local_tracking true --local_seeding_mask_type fa --local_fa_seeding_mask_threshold 0.1 --local_tracking_mask_type fa --local_fa_tracking_mask_threshold 0.1

DTI_Metrics .command.err:
/work/86/f93c47147f953d97aacdeb3f776a03/.comand.err

Matplotlib created a temporary config/cache directory at /tmpdata/matplotlib-i8ca0s7s because the default path (/path/to/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.
/usr/local/lib/python3.10/dist-packages/dipy/reconst/dti.py:498: RuntimeWarning: invalid value encountered in divide
  return 3 * np.sqrt(6) * determinant((A_squiggle / A_s_norm))
Traceback (most recent call last):
  File "/usr/local/bin/scil_compute_dti_metrics.py", line 33, in <module>
    sys.exit(load_entry_point('scilpy', 'console_scripts', 'scil_compute_dti_metrics.py')())
  File "/scilpy/scripts/scil_compute_dti_metrics.py", line 321, in main
    pis_mask = np.max(S0 < DWI, axis=-1)
  File "<__array_function__ internals>", line 180, in amax
  File "/usr/local/lib/python3.10/dist-packages/numpy/core/fromnumeric.py", line 2793, in amax
    return _wrapreduction(a, np.maximum, 'max', axis, None, out,
  File "/usr/local/lib/python3.10/dist-packages/numpy/core/fromnumeric.py", line 86, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
ValueError: zero-size array to reduction operation maximum which has no identity

compute_FRF .command.err:
/work/d9/898ca05b3092d00a1f0ab2f92c3d11/.comand.err

Matplotlib created a temporary config/cache directory at /tmpdata/matplotlib-5zlpzs6t because the default path (/path/to/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.
WARNING:root:No white matter mask specified! Only mask will be used (if it has been supplied). 
Be *VERY* careful about the estimation of the fiber response function to ensure no invalid voxel was used.
Traceback (most recent call last):
  File "/usr/local/bin/scil_compute_ssst_frf.py", line 33, in <module>
    sys.exit(load_entry_point('scilpy', 'console_scripts', 'scil_compute_ssst_frf.py')())
  File "/scilpy/scripts/scil_compute_ssst_frf.py", line 117, in main
    full_response = compute_ssst_frf(
  File "/scilpy/scilpy/reconst/frf.py", line 116, in compute_ssst_frf
    raise ValueError(
ValueError: Could not find at least 300 voxels with sufficient FA to estimate the FRF!

Not sure how to proceed.

Hi @ameenq05,

I can see you are using the scilus_1.5.0.sif container but I don’t know which version of tractoflow you are using.

I highly encourage you to use the latest version of tractoflow 2.4.3 with the latest container scilus_1.6.0.sif.

Hope it helps,
If you still have some issues don’t hesitate to ask.
Arnaud

Hi Arnaud,

I can’t find the version of tractoflow in the log or main directory files (cloned from GitHub), but I can try using the latest container.

I also wanted to ask, what’s the correct way to include both full AP and PA series data in the tractoflow input structure, and not just the reverse-phase b0?

-Ameen

Hello @ameenq05 ,

Go inside your tractoflow folder and run this command: git log -1 it will give me the latest commit you have and I’ll be able to help you.
The only way to include both full AP/PA is to create a BIDS structure.

Arnaud

Hi Arnaud, here’s what I get:

git log -1

commit 7d69625136fbdb9e32fddbd23997f88d144e374c (**HEAD ->** **master**, **origin/master**, **origin/HEAD**)
Merge: d67e568 312e0c2
Author: Arnaud Bore <arnaud.bore@gmail.com>
Date: Fri Dec 8 09:36:55 2023 -0500

Merge pull request #88 from arnaudbore/enh_memory_gpu_processes

[WIP] Increase memory and cpu for eddy and local_tracking processes

I’ll try with a BIDS input, thanks!

Perfect, you’re using the latest version.
It means you need to use the latest container (scilus_1.6.0.sif).

For the BIDS conversion you can use dcm2bids, I’m also developing this tool if you have any issues don’t hesitate.

Arnaud

Hi Arnaud,

I downloaded the newest container and tried again with BIDS input, and got the following error:

Workflow execution completed unsuccessfully!
The exit status of the task that caused the workflow execution to fail was: 127.

The full error message was:

Error executing process > 'Read_BIDS (Read_BIDS)'

Caused by:
  Process `Read_BIDS (Read_BIDS)` terminated with an error exit status (127)

Command executed:

  scil_validate_bids.py sub-0122 tractoflow_bids_struct.json                --readout 0.062                                                 -v

Command exit status:
  127

Command output:
  (empty)

Command error:
  .command.sh: line 2: scil_validate_bids.py: command not found

Work dir:
  /path/to/work/

Check your command line, you may have forgotten the singularity. -with-singularity

Hi Arnaud,

The singularity was included. See below:

nextflow run /path/to/main.nf --bids /path/to/BIDS/ --with-singularity /path/to/scilus_1.6.0.sif -profile fully_reproducible -resume -with-report report.html --run_pft_tracking false --run_local_tracking true --local_seeding_mask_type fa --local_fa_seeding_mask_threshold 0.1 --local_tracking_mask_type fa --local_fa_tracking_mask_threshold 0.1

Hi !

You wrote --with-singularity instead of -with-singularity

Arnaud

Hi Arnaud,

Thanks, I’ve made the correction and it runs until reaching topup, where it fails. See the log and .command.err below:

.command.err in work directory:


usage: scil_image_math.py [-h] [--data_type DATA_TYPE] [--exclude_background]
                          [-f] [-v]
                          {lower_threshold,upper_threshold,lower_threshold_eq,upper_threshold_eq,lower_clip,upper_clip,absolute_value,round,ceil,floor,normaliz
e_sum,normalize_max,log_10,log_e,convert,invert,addition,subtraction,multiplication,division,mean,std,correlation,union,intersection,difference,concatenate,dil
ation,erosion,closing,opening,blur}
                          in_images [in_images ...] out_image
scil_image_math.py: error: Output file sub-0122_ses-02__rev_b0_mean.nii.gz exists. Use -f to force overwriting

Log file output:


executor >  local (14)
[24/f7b402] process > README (README)                [100%] 1 of 1 ✔
[-        ] process > Bet_Prelim_DWI                 -
[a9/8cbde7] process > Denoise_DWI (sub-0122_ses-02)  [100%] 2 of 2 ✔
[-        ] process > Gibbs_correction               -
[06/c89deb] process > Prepare_for_Topup (sub-0122... [100%] 2 of 2 ✔
[18/8b182d] process > Topup (sub-0122_ses-02)        [100%] 4 of 4, failed: 4...
[28/ac7270] process > Prepare_dwi_for_eddy (sub-0... [100%] 1 of 1 ✔

For your reference, here is my BIDS root directory. It is validated using the online validator:


.
├── code
│   └── tractoflow_bids
├── dataset_description.json
├── derivatives
├── participants.json
├── participants.tsv
├── README
└── sub-0122
    └── ses-02
        ├── anat
        │   ├── sub-0122_ses-02_T1w.json
        │   └── sub-0122_ses-02_T1w.nii.gz
        └── dwi
            ├── sub-0122_ses-02_dir-AP_dwi.bval
            ├── sub-0122_ses-02_dir-AP_dwi.bvec
            ├── sub-0122_ses-02_dir-AP_dwi.json
            ├── sub-0122_ses-02_dir-AP_dwi.nii.gz
            ├── sub-0122_ses-02_dir-PA_dwi.bval
            ├── sub-0122_ses-02_dir-PA_dwi.bvec
            ├── sub-0122_ses-02_dir-PA_dwi.json
            └── sub-0122_ses-02_dir-PA_dwi.nii.gz

Hello,

It’s a bug in the pipeline. I’m going to fix it.
Can you send me the .command.sh please ? I want to confirm where the issue comes from.

There are two ways to fix it quickly.

  1. Do you have any fmap (AP-PA b0) ? If you had these two images it would run topup on these image and then run eddy without any error.

  2. You can edit the main script:

nano ~/.nextflow/assets/scilus/tractoflow/main.nf

Change Line 656

scil_image_math.py mean ${sid}__concatenated_rev_b0.nii.gz ${sid}__rev_b0_mean.nii.gz

to

scil_image_math.py mean ${sid}__concatenated_rev_b0.nii.gz ${sid}__rev_b0_mean.nii.gz -f

It should fix your issue.

Thank you for reporting the bug
Arnaud

Thanks Arnaud. Below is the .command.sh found in the topup work directory:


#!/bin/bash -ue
export OMP_NUM_THREADS=4
export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
export OPENBLAS_NUM_THREADS=1
export ANTS_RANDOM_SEED=1234

scil_image_math.py concatenate sub-0122_ses-02__rev_b0_mean.nii.gz sub-0122_ses-02__rev_b0_mean.nii.gz sub-0122_ses-02__concatenated_rev_b0.nii.gz
scil_image_math.py mean sub-0122_ses-02__concatenated_rev_b0.nii.gz sub-0122_ses-02__rev_b0_mean.nii.gz
antsRegistrationSyNQuick.sh -d 3 -f sub-0122_ses-02__b0_mean.nii.gz -m sub-0122_ses-02__rev_b0_mean.nii.gz -o output -t r -e 1
mv outputWarped.nii.gz sub-0122_ses-02__rev_b0_warped.nii.gz
scil_prepare_topup_command.py sub-0122_ses-02__b0_mean.nii.gz sub-0122_ses-02__rev_b0_warped.nii.gz          --config b02b0.cnf          --encoding_direction y          --readout 0.0959097 --out_prefix topup_results          --out_script
sh topup.sh
cp corrected_b0s.nii.gz sub-0122_ses-02__corrected_b0s.nii.gz

Do I need to include a rev-b0 and change the script? Or should only one of these changes suffice?

Only one of these option would work. If you don’t want to change your input, modifying the line will do the trick.

Arnaud

Hi Arnaud,

I had to change lines 656 and 657 to get it to work:

656      scil_image_math.py concatenate $rev_b0 $rev_b0 ${sid}__concatenated_rev_b0.nii.gz -f
657      scil_image_math.py mean ${sid}__concatenated_rev_b0.nii.gz ${sid}__rev_b0_mean.nii.gz -f

The preprocessing steps run to completion, but now the script runs into an issue with calculating the DTI Metrics (we’ve come full circle!). Below are the log, .command.err, and .command.sh files:

log:

executor >  local (25)
[73/9294c4] process > README (README)                [100%] 1 of 1 ✔
[-        ] process > Bet_Prelim_DWI                 -
[a9/8cbde7] process > Denoise_DWI (sub-0122_ses-02)  [100%] 2 of 2 ✔
[-        ] process > Gibbs_correction               -
[b1/74308f] process > Prepare_for_Topup (sub-0122... [100%] 2 of 2 ✔
[c4/1fe81a] process > Topup (sub-0122_ses-02)        [100%] 1 of 1 ✔
[6c/bd8cc6] process > Prepare_dwi_for_eddy (sub-0... [100%] 1 of 1 ✔
[95/5dde27] process > Eddy_Topup (sub-0122_ses-02)   [100%] 1 of 1 ✔
[-        ] process > Eddy                           -
[f9/1d310b] process > Bet_DWI (sub-0122_ses-02)      [100%] 1 of 1 ✔
[13/0e289d] process > N4_DWI (sub-0122_ses-02)       [100%] 1 of 1 ✔
[23/fbd88c] process > Crop_DWI (sub-0122_ses-02)     [100%] 1 of 1 ✔
[-        ] process > Denoise_T1                     -
[9c/28aa39] process > N4_T1 (sub-0122_ses-02)        [100%] 1 of 1 ✔
[4f/b3dca6] process > Resample_T1 (sub-0122_ses-02)  [100%] 1 of 1 ✔
[9d/c0a7e0] process > Bet_T1 (sub-0122_ses-02)       [100%] 1 of 1 ✔
[16/c9355a] process > Crop_T1 (sub-0122_ses-02)      [100%] 1 of 1 ✔
[48/9dfb2d] process > Normalize_DWI (sub-0122_ses... [100%] 1 of 1 ✔
[7a/0c65ef] process > Resample_DWI (sub-0122_ses-02) [100%] 1 of 1 ✔
[20/938cd0] process > Extract_B0 (sub-0122_ses-02)   [100%] 1 of 1 ✔
[-        ] process > Extract_SH_Fitting_Shell       -
[-        ] process > SH_Fitting                     -
[c2/6454c9] process > Extract_DTI_Shell (sub-0122... [100%] 1 of 1 ✔
[8d/6f2711] process > DTI_Metrics (sub-0122_ses-02)  [100%] 4 of 4, failed: 4...
[9f/f48334] process > Extract_FODF_Shell (sub-012... [100%] 1 of 1 ✔
[-        ] process > Register_T1                    -
[-        ] process > Register_Freesurfer            -
[-        ] process > Segment_Freesurfer             -
[-        ] process > Segment_Tissues                -
[44/bba921] process > Compute_FRF (sub-0122_ses-02)  [100%] 1 of 1 ✔
[-        ] process > Mean_FRF                       -
[-        ] process > FODF_Metrics                   -
[-        ] process > PFT_Tracking_Maps              -
[-        ] process > PFT_Seeding_Mask               -
[-        ] process > PFT_Tracking                   -
[-        ] process > Local_Tracking_Mask            -
[-        ] process > Local_Seeding_Mask             -
[-        ] process > Local_Tracking                 -
Pipeline completed at: 2024-07-04T17:07:41.072+04:00
Execution status: OK
Execution duration: 1h 20m 33s
[8d/6f2711] NOTE: Process `DTI_Metrics (sub-0122_ses-02)` terminated with an error exit status (1) -- Error is ignored

Completed at: 04-Jul-2024 17:07:42
Duration    : 1h 20m 34s
CPU hours   : 6.9 (4.1% failed)
Succeeded   : 21
Ignored     : 1
Failed      : 4

.command.err in work/8d/6f2711a0da3bfa7a5303373cb706b0:


/usr/local/lib/python3.10/dist-packages/dipy/reconst/dti.py:499: RuntimeWarning: invalid value encountered in divide
  return 3 * np.sqrt(6) * determinant((A_squiggle / A_s_norm))
Traceback (most recent call last):
  File "/usr/local/bin/scil_compute_dti_metrics.py", line 33, in <module>
    sys.exit(load_entry_point('scilpy', 'console_scripts', 'scil_compute_dti_metrics.py')())
  File "/scilpy/scripts/scil_compute_dti_metrics.py", line 321, in main
    pis_mask = np.max(S0 < DWI, axis=-1)
  File "<__array_function__ internals>", line 180, in amax
  File "/usr/local/lib/python3.10/dist-packages/numpy/core/fromnumeric.py", line 2793, in amax
    return _wrapreduction(a, np.maximum, 'max', axis, None, out,
  File "/usr/local/lib/python3.10/dist-packages/numpy/core/fromnumeric.py", line 86, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
ValueError: zero-size array to reduction operation maximum which has no identity

.command.sh in work/8d/6f2711a0da3bfa7a5303373cb706b0:


#!/bin/bash -ue
export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
export OMP_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
scil_compute_dti_metrics.py sub-0122_ses-02__dwi_dti.nii.gz sub-0122_ses-02__bval_dti sub-0122_ses-02__bvec_dti --mask sub-0122_ses-02__b0_mask_resampled.nii.gz        --ad sub-0122_ses-02__ad.nii.gz --evecs sub-0122_ses-02__evecs.nii.gz        --evals sub-0122_ses-02__evals.nii.gz --fa sub-0122_ses-02__fa.nii.gz        --ga sub-0122_ses-02__ga.nii.gz --rgb sub-0122_ses-02__rgb.nii.gz        --md sub-0122_ses-02__md.nii.gz --mode sub-0122_ses-02__mode.nii.gz        --norm sub-0122_ses-02__norm.nii.gz --rd sub-0122_ses-02__rd.nii.gz        --tensor sub-0122_ses-02__tensor.nii.gz        --non-physical sub-0122_ses-02__nonphysical.nii.gz        --pulsation sub-0122_ses-02__pulsation.nii.gz        --residual sub-0122_ses-02__residual.nii.gz        -f --force_b0_threshold

Hello @ameenq05 ,

I suspect your acquisition have quite high b-values. Is that right ? What is your protocol (bvals and number of directions for each bval) ? If so, you need to add this option: --max_dti_shell_value in your command line and put the lowest b-value possible which gather at least 15 directions ish.
If it doesn’t fix your issue I will need to have access to the data at some point I guess.

Hope it helps,
Arnaud