PennLINC / qsiprep

Preprocessing of diffusion MRI
http://qsiprep.readthedocs.io
BSD 3-Clause "New" or "Revised" License
138 stars 55 forks source link

dwi2response: [ERROR] Error trying to calculate statistic 'median' from image 'safe_sdm.mif' #337

Closed m-petersen closed 2 years ago

m-petersen commented 2 years ago

Hi Matt,

I am currently working on a qsiprep implementation with the datalad fairly big processing workflow. I use the containerized qsiprep version 0.14.2 from dockerhub.

I face an error with dwi2response that I don't exactly know how to diagnose and address.

QSIRecon failed: Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node
    result["result"] = node.run(updatehash=updatehash)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 516, in run
    result = self._run_interface(execute=True)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 635, in _run_interface
    return self._run_command(execute)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 741, in _run_command
    result = self._interface.run(cwd=outdir)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 428, in run
    runtime = self._run_interface(runtime)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 822, in _run_interface
    self.raise_exception(runtime)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 753, in raise_exception
    ).format(**runtime.dictcopy())
RuntimeError: Command:
dwi2response dhollander -mask /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/resample_mask/sub-1332847_desc-brain_mask_resample.nii.gz -nthreads 16 /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/create_mif/sub-1332847_ses-2_space-T1w_desc-preproc_dwi.mif sub-1332847_ses-2_space-T1w_desc-preproc_dwi_wm.txt sub-1332847_ses-2_space-T1w_desc-preproc_dwi_gm.txt sub-1332847_ses-2_space-T1w_desc-preproc_dwi_csf.txt
Standard output:

Standard error:
mrstats: [ERROR] Cannot output statistic of interest; no values read (empty mask?)
dwi2response:
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response:
dwi2response: Generated scratch directory: /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/estimate_response/dwi2response-tmp-BW9OCD/
dwi2response: Importing DWI data (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/create_mif/sub-1332847_ses-2_space-T1w_desc-preproc_dwi.mif)...
dwi2response: Importing mask (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/resample_mask/sub-1332847_desc-brain_mask_resample.nii.gz)...
dwi2response: Changing to scratch directory (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/estimate_response/dwi2response-tmp-BW9OCD/)
dwi2response: -------
dwi2response: 2 unique b-value(s) detected: 5,2003 with 3,3 volumes
dwi2response: -------
dwi2response: Preparation:
dwi2response: * Eroding brain mask by 3 pass(es)...
dwi2response:   [ mask: 189527 -> 144840 ]
dwi2response: * Computing signal decay metric (SDM):
dwi2response:  * b=5...
dwi2response:  * b=2003...
dwi2response: * Removing erroneous voxels from mask and correcting SDM...
dwi2response:   [ mask: 144840 -> 144817 ]
dwi2response: -------
dwi2response: Crude segmentation:
dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 144817 -> 144817 (WM) & 0 (GM-CSF) ]
dwi2response: * Crude GM versus CSF separation...

dwi2response: [ERROR] Error trying to calculate statistic 'median' from image 'safe_sdm.mif'
Return code: 1

Traceback (most recent call last):
  File "/usr/local/miniconda/bin/qsiprep", line 8, in <module>
    sys.exit(main())
  File "/usr/local/miniconda/lib/python3.7/site-packages/qsiprep/cli/run.py", line 721, in main
    qsirecon_post_wf.run(**plugin_settings)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/workflows.py", line 638, in run
    runner.run(execgraph, updatehash=updatehash, config=self.config)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/plugins/base.py", line 164, in run
    self._clean_queue(jobid, graph, result=result)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/plugins/base.py", line 227, in _clean_queue
    raise RuntimeError("".join(result["traceback"]))
RuntimeError: Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node
    result["result"] = node.run(updatehash=updatehash)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 516, in run
    result = self._run_interface(execute=True)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 635, in _run_interface
    return self._run_command(execute)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/pipeline/engine/nodes.py", line 741, in _run_command
    result = self._interface.run(cwd=outdir)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 428, in run
    runtime = self._run_interface(runtime)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 822, in _run_interface
    self.raise_exception(runtime)
  File "/usr/local/miniconda/lib/python3.7/site-packages/nipype/interfaces/base/core.py", line 753, in raise_exception
    ).format(**runtime.dictcopy())
RuntimeError: Command:
dwi2response dhollander -mask /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/resample_mask/sub-1332847_desc-brain_mask_resample.nii.gz -nthreads 16 /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/create_mif/sub-1332847_ses-2_space-T1w_desc-preproc_dwi.mif sub-1332847_ses-2_space-T1w_desc-preproc_dwi_wm.txt sub-1332847_ses-2_space-T1w_desc-preproc_dwi_gm.txt sub-1332847_ses-2_space-T1w_desc-preproc_dwi_csf.txt
Standard output:

Standard error:
mrstats: [ERROR] Cannot output statistic of interest; no values read (empty mask?)
dwi2response:
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response:
dwi2response: Generated scratch directory: /dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/estimate_response/dwi2response-tmp-BW9OCD/
dwi2response: Importing DWI data (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/create_mif/sub-1332847_ses-2_space-T1w_desc-preproc_dwi.mif)...
dwi2response: Importing mask (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/resample_mask/sub-1332847_desc-brain_mask_resample.nii.gz)...
dwi2response: Changing to scratch directory (/dev/shm/1332847.10143549/ds/qsirecon_wf/sub-1332847_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwi_file_..dev..shm..1332847.10143549..ds..outputs..qsiprep..sub-1332847..ses-2..dwi..sub-1332847_ses-2_space-T1w_desc-preproc_dwi.nii.gz/estimate_response/dwi2response-tmp-BW9OCD/)
dwi2response: -------
dwi2response: 2 unique b-value(s) detected: 5,2003 with 3,3 volumes
dwi2response: -------
dwi2response: Preparation:
dwi2response: * Eroding brain mask by 3 pass(es)...
dwi2response:   [ mask: 189527 -> 144840 ]
dwi2response: * Computing signal decay metric (SDM):
dwi2response:  * b=5...
dwi2response:  * b=2003...
dwi2response: * Removing erroneous voxels from mask and correcting SDM...
dwi2response:   [ mask: 144840 -> 144817 ]
dwi2response: -------
dwi2response: Crude segmentation:
dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 144817 -> 144817 (WM) & 0 (GM-CSF) ]
dwi2response: * Crude GM versus CSF separation...

dwi2response: [ERROR] Error trying to calculate statistic 'median' from image 'safe_sdm.mif'
Return code: 1

Apparently the mask isn't empty but GM-CSF is. I get the same error for all 16 test subjects I run in parallel. Inputs (uk biobank data) look fine to me.

Some notes to my approach. I use a slightly customized version from https://github.com/psychoinformatics-de/fairly-big-processing-workflow and I don't find an obvious reason why dwi2response is failing and while the remaining non-dependent commands succeeds. Max RAM usage of the jobs is 256gb of 1 TB available. Suspecting an issue with /tmp I mounted the current directory /dev/shm/1332847.10143549/ds/ also as /tmp and set it explicitly with SINGULARITYENV_TMPFILE_DIR=/tmp. If I run qsiprep without any datalad involvement on the same system it runs without issues pointing me away from qsiprep and towards datalad / the fairly big workflow. However, to get an idea of what might be wrong I hope for some insights from your site. I'll also post on the mrtrix forum.

Happy to provide further details and grateful for any input.

Cheers, Marvin

call_qsiprep.sh

#!/bin/bash

subid=$1

# pybids (inside qsiprep) will try to read all JSON files in a dataset. In case
# of a recomputation, JSON files of other subjects can be dangling symlinks.
# We prevent pybids from crashing the qsiprep run when it can't read those, by
# wiping them out temporarily via renaming.
# We spare only those that belong to the participant we want to process.
# After job completion, the jsons will be restored.
# See https://github.com/bids-standard/pybids/issues/631 for more information.
# NO PROBLEM HERE, as only jsons of the respective subject are available
export SINGULARITYENV_MRTRIX_TMPFILE_DIR=/tmp

# execute qsiprep. Its runscript is available as /singularity within the
# container. Custom qsiprep parametrization can be done here.
/singularity inputs/ukb outputs participant    --work-dir .    --participant-label $subid    --recon_input outputs/qsiprep/${subid}    --recon_spec mrtrix_multishell_msmt    --denoise-method dwidenoise    --unringing-method mrdegibbs    --skull-strip-template OASIS    --use-syn-sdc    --force-syn    --output-space T1w template    --template MNI152NLin2009cAsym    --output-resolution 2    --nthreads 16    --omp-nthreads 16    --mem_mb 16000    --stop-on-first-crash    --fs-license-file code/license.txt    --skull-strip-fixed-seed --skip_bids_validation    --notrack
m-petersen commented 2 years ago

https://community.mrtrix.org/t/dwi2response-error-error-trying-to-calculate-statistic-median-from-image-safe-sdm-mif/5508

jdtournier commented 2 years ago

To investigate further, I would definitely urge the addition of the -debug flag to the dwi2response call. That should provide a lot more diagnostic information for us to look through.

In the meantime, I've had a brief go at following the code path that would lead to this, and I'm struggling a bit. There's an initial error right at the beginning:

mrstats: [ERROR] Cannot output statistic of interest; no values read (empty mask?)

which seems highly relevant... Since this occurs before the invocation of the dhollander algorithm proper (which starts by printing out a line of hyphens), I suspect this occurs at this point where we try to count the number of voxels in the mask. This essentially invokes the mrstats command, which would indeed print the message seen if the mask is empty. But it would also throw an exception, which should result in a non-zero exit code for that command. This should in turn be caught in the Python code and raise an exception there, along with the message "Error trying to calculate statistics from image blah" - which we don't see in your output.

So in short, I can't make sense of this. Note that all this is looking at the code in the latest release of MRtrix 3.0.3 - things might be different using older versions. What version was used here?

But clearly, the fact that things work as expected when run without datalad suggests that the issue is not necessarily within the dwi2response script per se - though I really fail to understand how.

You state that the mask provided to the command is non-empty, yet when the command runs, the first thing it does is invoke mrstats on that mask, and the error message clearly indicates that there were no errors reading the file as such (I would expect to see error messages related to I/O if that were the case), but that the mask genuinely contains nothing but zeros. I don't see how we can reconcile these two observations. If there were filesystem or I/O related issues introduced by the use of datalad, I would expect errors to occur much earlier, when reading the image header...? :confused:

mattcieslak commented 2 years ago

Our group also uses the fairly-big workflow, and I think you're on the right track with looking into how /tmp gets mounted in the container. We've seen some truly bizarre read/write errors with the mixture of datalad and singularity. You're mounting your datalad dataset at /tmp, but also using this as the mrtrix tmpdir? That makes me nervous. Could you try using a directory external to your datalad dataset as the mrtrix tempdir? It's also a good idea to use a working directory somewhere that's not tracked by datalad. Our group's implementation is here with

If you're just getting started processing this data, I think you should upgrade to qsiprep-0.14.3. This release includes mrtrix 3.0.3, which will be helpful to check @jdtournier's suggestions.

unrelated notes: Doesn't UKB have EPI fieldmaps? I'd strongly recommend using actual fieldmap acquisitions instead of SyN-SDC. In our testing the SyN displacements have had relatively low correlation with actual fieldmap displacements.

It looks like you're doing reconstruction and preprocessing in a single commandline run. I'd suggest removing the --recon-input argument. It currently points to a subject directory, which would cause it to crash if you tried using a similar command with --recon-only.

m-petersen commented 2 years ago

Thanks so much for your thorough answers.

@mattcieslak

Thanks for sharing your workflow. I'll try to adapt my approach accordingly. I will follow your recommendations and set /tmp somewhere outside the datalad dataset, unset `--recon-input and upgrade to qsiprep 0.14.3.

UKB does indeed have EPI field maps. However I do not have access to the field maps because I am working on a collaborators datalad dataset missing them.

@jdtournier

The mask and the image aren't empty. The according screenshots are attached. However, while looking into them again I noticed a problem. The DWI input into dwi2response only has 6 volumes. I checked it visually and I am quite sure that qsiprep took the PA acquired DWI (6 volumes) instead of the AP acquired one (105 volumes) and preprocessed it. This might cause the error. After some research in the documentation I am not quite sure why qsiprep takes the PA instead of the AP as they follow the BIDS standard as far as I can tell (according filenames in the below code cells).

mrinfo outputs of one subject:

Raw AP dwi

************************************************
Image name:          "sub-1333562_ses-2_acq-AP_dwi.nii.gz"
************************************************
  Dimensions:        104 x 104 x 72 x 105
  Voxel size:        2.01923 x 2.01923 x 2 x 3.6
  Data strides:      [ -1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         unsigned 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9987   -0.003835     0.05051      -104.8
                          0.02193      0.9316     -0.3629      -91.83
                         -0.04566      0.3635      0.9305      -133.5
  comments:          TE=92;Time=174255.383;phase=1;dwell=0.670

Raw PA dwi

************************************************
Image name:          "sub-1333562_ses-2_acq-PA_dwi.nii.gz"
************************************************
  Dimensions:        104 x 104 x 72 x 6
  Voxel size:        2.01923 x 2.01923 x 2 x 3.6
  Data strides:      [ -1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         unsigned 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9987   -0.003835     0.05051      -104.8
                          0.02193      0.9316     -0.3629      -91.83
                         -0.04566      0.3635      0.9305      -133.5
  comments:          TE=92;Time=174215.358;phase=0;dwell=0.670

Preprocessed DWI (aligned with ACPC?)

************************************************
Image name:          "sub-1333562_ses-2_space-T1w_desc-preproc_dwi.nii.gz"
************************************************
  Dimensions:        97 x 116 x 98 x 6
  Voxel size:        2 x 2 x 2 x 0
  Data strides:      [ -1 -2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1          -0           0         -97
                               -0           1           0      -127.5
                               -0          -0           1         -85

Before that I tried to reproduce the error locally with the -debug flag. For the sake of completeness I attached the error outputs. Interestingly a different error occurs when I use my locally installed mrtrix3 version 3.0.3. Please find the error output in the attached error_mrtrix303_local.txt. However, when I use the qsiprep container I can reproduce the error (error_qsiprep_local.txt). The initial error is in error_qsiprep_hpc.txt.

Attachments

sub-1000201_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwifile..dev..shm..1000201.10149155..ds..outputs..qsiprep..sub-1000201..ses-2..dwi..sub-1000201_ses-2_space-T1w_desc-preproc_dwi.nii.gz/resample_mask/sub-1000201_desc-brain_mask_resample.nii.gz dwi2response0000

sub-1000201_mrtrix_multishell_msmt/recon_wf/msmt_csd/_dwifile..dev..shm..1000201.10149155..ds..outputs..qsiprep..sub-1000201..ses-2..dwi..sub-1000201_ses-2_space-T1w_desc-preproc_dwi.nii.gz/create_mif/sub-1000201_ses-2_space-T1w_desc-preproc_dwi.mif dwi2response0001

error_mrtrix303_local.txt error_qsiprep_hpc.txt error_qsiprep_local.txt

mattcieslak commented 2 years ago

This sounds like there's something strange in the BIDS inputs. Maybe they acquired a scan that was intended to be an opposite blip b=0 series for distortion correction and put it in the dwi/ directory?

m-petersen commented 2 years ago

Yes there are the DWI (blip up, sub-1333562_ses-2_acq-AP_dwi.nii.gz) and a reverse phase encoded image (blip down, sub-1333562_ses-2_acq-PA_dwi.nii.gz) in the dwi/ directory.

mattcieslak commented 2 years ago

And one only has 6 images in it? If so, the 6-image series should probably be moved into the fmap/ directory, have its extension changed to _epi and have the >6 image series listed in the "IntendedFor" field of its sidecar. This will get you distortion correction and may solve the reconstruction issue as well.

To circumvent an issue with BIDS, you can add the .bval and .bvec file to the fmap/ directory as well and add **/fmap/*.bv* to ./bidsignore.

m-petersen commented 2 years ago

Following your recommendations resolved the issue. Thanks again for your help!