NODDI recon producing weird values for ICVF and ISOVF

weberam2 commented 1 year ago

Hello!

I have the following DTI scheme:

-> Loading data:
    * DWI signal
        - dim    = 108 x 129 x 109 x 96
        - pixdim = 1.800 x 1.800 x 1.800
    * Acquisition scheme
        - 96 samples, 2 shells
        - 6 @ b=0 , 60 @ b=2000.0 , 30 @ b=1000.0 
    * Binary mask
        - dim    = 108 x 129 x 109
        - pixdim = 1.800 x 1.800 x 1.800
        - voxels = 206359

Images were acquired on a GE scanner Two separate sequences were acquired:

b=1000 with 30 directions, in one phase direction,
then b=2000 with 60 directions in opposite phase direction

QSIprep was used to preprocess the images and they are concatenated into one .nii.gz file (which is default behaviour for QSIprep)

running AMICO recon runs no problem

However, the ISOVF and ICVF look like they are being miscalculated in some regions over others, in a very obvious way. Best way to visualize is with an image:

What am I doing wrong?

mattcieslak commented 1 year ago

Hi @weberam2 how does the distortion correction look? There are b=0 images from both phase encoding directions, right? I also wonder if there's something odd with the colormap in your image plotting tool. What does the histogram of ICVF and ISOVF values look like inside the brain mask?

weberam2 commented 1 year ago

Hi @mattcieslak

Thanks for getting back to me!

"how does the distortion correction look?" It looks great!

"here are b=0 images from both phase encoding directions, right?" Correct: The two scans are as follows:

b=0 x 3, then b=1000 x 30, all in AP direction
then a separate scan: b=0 x 3, then b=2000 x 60, all in PA direction

" I also wonder if there's something odd with the colormap in your image plotting tool. What does the histogram of ICVF and ISOVF values look like inside the brain mask?" No, this does not seem to be the case. Here is the histogram from ICVF:

I do have an update: I managed to do the preprocessing with MRtrix3:

mrcat sub-Rett01_dir-AP_dwi.mif sub-Rett01_dir-PA_dwi.mif all_dwis.mif -axis 3
dwidenoise all_dwis.mif all_dwis_denoise.mif -noise noise.mif
dwifslpreproc all_dwis.mif dwi_preprocessed.mif -rpe_header -eddy_options " --slm=linear "

And when I run AMICO NODDI, my images look like this (and the histogram) ICVF:

So I do believe it is something happening with QSIprep preprocessing...

Any suggestions of what we could try? I am happy to share the dwi .nii.gz / .json and the anat .nii.gz / .json files

mattcieslak commented 1 year ago

with default settings the qsiprep workflow is extremely similar to the mrtrix preprocessing pipeline. I wonder if it is the b=0 image harmonization, which isn't present in the mrtrix version of the workflow. Could you try rerunning preprocessing with --no-b0-harmonization?

weberam2 commented 1 year ago

I'm running that now, and will let you know the results But in the mean-time, I also realize the preprocessing doesn't finish properly

230217-21:39:41,928 nipype.workflow ERROR:
could not run node: qsiprep_wf.single_subject_Rett01_wf.dwi_preproc_wf.pre_hmc_wf.raw_rpe_concat

But it still outputs the preprocessed files:

sub-Rett01_space-T1w_desc-eddy_cnr.nii.gz
sub-Rett01_space-T1w_desc-preproc_dwi.bvec
sub-Rett01_desc-SliceQC_dwi.json
sub-Rett01_space-T1w_desc-preproc_dwi.b
sub-Rett01_space-T1w_desc-preproc_dwi.nii.gz
sub-Rett01_space-T1w_desc-brain_mask.nii.gz
sub-Rett01_space-T1w_desc-preproc_dwi.bval
sub-Rett01_space-T1w_dwiref.nii.gz

weberam2 commented 1 year ago

more info of the error:

Traceback:                                                                                                                                                                                                 
        Traceback (most recent call last):                                                                                                                                                                 
          File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/base/core.py", line 398, in run                                                                                         
            runtime = self._run_interface(runtime)                                                                                                                                                         
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/interfaces/nilearn.py", line 134, in _run_interface                                                                               
            new_nii = concat_imgs(self.inputs.in_files, dtype=self.inputs.dtype)                                                                                                                           
          File "/usr/local/miniconda/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py", line 487, in concat_niimgs                                                                          
            for index, (size, niimg) in enumerate(zip(lengths, _iter_check_niimg(                                                                                                                          
          File "/usr/local/miniconda/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py", line 160, in _iter_check_niimg                                                                      
            raise ValueError(                                                                                                                                                                              
        ValueError: Field of view of image #1 is different from reference FOV.                                                                                                                             
        Reference affine:                                                                                                                                                                                  
        array([[-1.79477441e+00, -6.63911998e-02,  1.19901352e-01,                                                                                                                                         
                 1.28985382e+02],                                                                                                                                                                          
               [-8.74731839e-02,  1.76716912e+00, -3.30843806e-01,                                                                                                                                         
                -9.69735107e+01],                                                                                                                                                                          
               [ 1.05513543e-01,  3.35716337e-01,  1.76523221e+00,                                                                                                                                         
                -1.38833572e+02],                                                                                                                                                                          
               [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,                                                                                                                                         
                 1.00000000e+00]])                                                                                                                                                                         
        Image affine:                                                                                                                                                                                      
        array([[-1.79477453e+00, -6.63908795e-02,  1.19907245e-01,                                                                                                                                         
                 1.28985336e+02],                                                                                                                                                                          
               [-8.74728560e-02,  1.76716900e+00, -3.30862254e-01,                                                                                                                                         
                -9.69734955e+01],                                                                                                                                                                          
               [ 1.05512999e-01,  3.35716873e-01,  1.76532757e+00,                                                                                                                                         
                -1.38833649e+02],                                                                                                                                                                          
               [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,                                                                                                                                         
                 1.00000000e+00]])                                                                                                                                                                         
        Reference shape:                                                                                                                                                                                   
        (140, 140, 80)                                                                                                                                                                                     
        Image shape:                                                                                                                                                                                       
        (140, 140, 80, 33)

mattcieslak commented 1 year ago

What version are you using? There was a similar error a few versions ago

weberam2 commented 1 year ago

What version are you using? There was a similar error a few versions ago

pennbbl/qsiprep:0.16.1

cookpa commented 1 year ago

@mattcieslak is there a solution to this when the source data has different orientation?

@weberam2 I've seen this in some data as well. If you have the DICOM headers, you can compare ImagePositionPatient and ShimSetting - in affected data I've seen, there's been a re-shim and change in alignment of the image FOV.

weberam2 commented 1 year ago

Just a quick update, I ran the QSIprep with the option: --no-b0-harmonization and still get the bad NODDI images

weberam2 commented 1 year ago

@mattcieslak is there a solution to this when the source data has different orientation?

@weberam2 I've seen this in some data as well. If you have the DICOM headers, you can compare ImagePositionPatient and ShimSetting - in affected data I've seen, there's been a re-shim and change in alignment of the image FOV.

OK, so the first series: REL Image Position Patient//-121.315-144.362-69.2205 The second: REL Image Position Patient//-127.311-127.82\19.0426

For ShimSetting, I can't find anything with that in the dicom files...

cookpa commented 1 year ago

Sorry, I got mixed up - ShimSetting is in the BIDS JSON sidecar

weberam2 commented 1 year ago

I don't see the ShimSetting in there either.

Here are both .json files:

{
    "Modality": "MR",
    "MagneticFieldStrength": 3,
    "ImagingFrequency": 127.773,
    "Manufacturer": "GE",
    "InternalPulseSequenceName": "EPI2",
    "ManufacturersModelName": "DISCOVERY MR750",
    "InstitutionName": "BCCH MRI Research Facility",
    "DeviceSerialNumber": "00000604875MR750",
    "StationName": "MR01OC01",
    "BodyPartExamined": "BRAIN",
    "PatientPosition": "HFS",
    "SoftwareVersions": "27\\LX\\MR Software release:DV26.0_R06_2132.a",
    "MRAcquisitionType": "2D",
    "SeriesDescription": "DTI B1000_30dir_AP",
    "ProtocolName": "RESEARCH - AMRETT",
    "ScanningSequence": "EP\\SE",
    "SequenceVariant": "NONE",
    "ScanOptions": "CL_GEMS\\SAT_GEMS\\EDR_GEMS\\EPI_GEMS\\HYPERBAND_GEMS\\ACC_GEMS\\PFF\\FS",
    "ImageType": [
        "ORIGINAL",
        "PRIMARY",
        "OTHER"
    ],
    "SeriesNumber": 7,
    "AcquisitionTime": "10:07:41.000000",
    "AcquisitionNumber": 1,
    "SliceThickness": 1.8,
    "SpacingBetweenSlices": 1.8,
    "SAR": 0.219551,
    "EchoTime": 0.082,
    "RepetitionTime": 4.1,
    "FlipAngle": 90,
    "PhaseEncodingPolarityGE": "Unflipped",
    "CoilString": "Nova32ch",
    "MultibandAccelerationFactor": 3,
    "PercentPhaseFOV": 100,
    "PercentSampling": 100,
    "AcquisitionMatrixPE": 140,
    "ReconMatrixPE": 140,
    "ParallelReductionFactorInPlane": 2,
    "EffectiveEchoSpacing": 0.000382072,
    "TotalReadoutTime": 0.053108,
    "PixelBandwidth": 3571.43,
    "PhaseEncodingDirection": "j-",
    "ImageOrientationPatientDICOM": [
        0.997097,
        0.048596,
        0.0586183,
        -0.0368838,
        0.981761,
        -0.186509
    ],
    "InPlanePhaseEncodingDirectionDICOM": "COL",
    "ConversionSoftware": "dcm2niix",
    "ConversionSoftwareVersion": "v1.0.20211006",
    "Dcm2bidsVersion": "2.1.6"
}

{
    "Modality": "MR",
    "MagneticFieldStrength": 3,
    "ImagingFrequency": 127.773,
    "Manufacturer": "GE",
    "InternalPulseSequenceName": "EPI2",
    "ManufacturersModelName": "DISCOVERY MR750",
    "InstitutionName": "BCCH MRI Research Facility",
    "DeviceSerialNumber": "00000604875MR750",
    "StationName": "MR01OC01",
    "BodyPartExamined": "BRAIN",
    "PatientPosition": "HFS",
    "SoftwareVersions": "27\\LX\\MR Software release:DV26.0_R06_2132.a",
    "MRAcquisitionType": "2D",
    "SeriesDescription": "DTI B2000_60dir_PA",
    "ProtocolName": "RESEARCH - AMRETT",
    "ScanningSequence": "EP\\SE",
    "SequenceVariant": "NONE",
    "ScanOptions": "CL_GEMS\\SAT_GEMS\\EDR_GEMS\\EPI_GEMS\\HYPERBAND_GEMS\\ACC_GEMS\\PFF\\FS",
    "ImageType": [
        "ORIGINAL",
        "PRIMARY",
        "OTHER"
    ],
    "SeriesNumber": 8,
    "AcquisitionTime": "10:10:26.000000",
    "AcquisitionNumber": 1,
    "SliceThickness": 1.8,
    "SpacingBetweenSlices": 1.8,
    "SAR": 0.219551,
    "EchoTime": 0.082,
    "RepetitionTime": 4.1,
    "FlipAngle": 90,
    "PhaseEncodingPolarityGE": "Flipped",
    "CoilString": "Nova32ch",
    "MultibandAccelerationFactor": 3,
    "PercentPhaseFOV": 100,
    "PercentSampling": 100,
    "AcquisitionMatrixPE": 140,
    "ReconMatrixPE": 140,
    "ParallelReductionFactorInPlane": 2,
    "EffectiveEchoSpacing": 0.000382072,
    "TotalReadoutTime": 0.053108,
    "PixelBandwidth": 3571.43,
    "PhaseEncodingDirection": "j",
    "ImageOrientationPatientDICOM": [
        0.997097,
        0.0485962,
        0.0586186,
        -0.036884,
        0.981761,
        -0.186509
    ],
    "InPlanePhaseEncodingDirectionDICOM": "COL",
    "ConversionSoftware": "dcm2niix",
    "ConversionSoftwareVersion": "v1.0.20211006",
    "Dcm2bidsVersion": "2.1.6"
}

weberam2 commented 1 year ago

If I run diff sub-Rett01_dir-AP_dwi.json sub-Rett01_dir-PA_dwi.json on these I get:

15c15
<     "SeriesDescription": "DTI B1000_30dir_AP",
---
>     "SeriesDescription": "DTI B2000_60dir_PA",
25,26c25,26
<     "SeriesNumber": 7,
<     "AcquisitionTime": "10:07:41.000000",
---
>     "SeriesNumber": 8,
>     "AcquisitionTime": "10:10:26.000000",
34c34
<     "PhaseEncodingPolarityGE": "Unflipped",
---
>     "PhaseEncodingPolarityGE": "Flipped",
45c45
<     "PhaseEncodingDirection": "j-",
---
>     "PhaseEncodingDirection": "j",
48,50c48,50
<         0.048596,
<         0.0586183,
<         -0.0368838,
---
>         0.0485962,
>         0.0586186,
>         -0.036884,

cookpa commented 1 year ago

Ah, I see you have a GE scanner, I'm looking at Siemens data here. I guess the ShimSetting might be a Siemens-specific output from dcm2niix.

Anyway, it looks like your patient orientation is identical between series - this is good because it means the bvecs should have the same coordinate system. Sorry to derail your thread!

weberam2 commented 1 year ago

Ah, I see you have a GE scanner, I'm looking at Siemens data here. I guess the ShimSetting might be a Siemens-specific output from dcm2niix.

Anyway, it looks like your patient orientation is identical between series - this is good because it means the bvecs should have the same coordinate system. Sorry to derail your thread!

Don't be sorry! I'll take any help I can get. I like the idea of QSIprep and would like to see it work...

mattcieslak commented 1 year ago

@weberam2 do you know if the mrtrix pipeline is bias correcting the images individually before running eddy? I wonder if there is something happening in the bias correction that is making the b>0 images have very different intensities in them across the two scans.

QSIPrep will refuse to concatenate images with different phase encoding directions before denoising (which for us means unringing, mp-pca and b1 bias correction). Maybe --dwi-no-biascorr in addition to --no-b0-harmonization will make this less of a problem?

weberam2 commented 1 year ago

That worked!!

So, this is what I ran:

docker run --rm -v /mnt/WeberLab/license.txt:/opt/freesurfer/license.txt:ro -v $PWD:/data:ro -v $PWD/derivatives/FastSurfer:/sngl/freesurfer-input:ro -v $PWD/derivatives:/out -v $HOME/work:/scratch pennbbl/qsiprep:0.16.1 /data /out participant --freesurfer-input /sngl/freesurfer-input --output-resolution 1.8 --participant-label sub-Rett01 -w /scratch --no-b0-harmonization --dwi-no-biascorr

and then

docker run --rm -v /mnt/WeberLab/license.txt:/opt/freesurfer/license.txt:ro -v $PWD:/data:ro -v $PWD/derivatives/qsiprep:/qsiprep-output:ro -v $PWD/derivatives:/out -v $HOME/work:/scratch pennbbl/qsiprep:0.16.1 /data /out participant --recon-input /qsiprep-output --recon-spec amico_noddi --participant-label sub-Rett01 --recon-only -w /scratch

and now I get nice NODDI images (here is the ISOVF):

I still have some questions though:

What are the negative consequences of using --dwi-no-biascorr and --no-b0-harmonization for my preprocessing?
My qsiprep preprocessing still doesn't spit out a boilerplate, even though it does seem to finish the preprocessing (attaching the crash report). Why is that? crash-20230219-013059-root-raw_rpe_concat-42fbb6b6-65ca-4d00-be47-71d3bd6bfdbe.txt

cookpa commented 1 year ago

What are the negative consequences of using --dwi-no-biascorr and --no-b0-harmonization for my preprocessing?

B0 harmonization is designed to reduce signal drift over time. I believe this was inspired by the HCP pipelines, where there are typically four DWI series taking a total of 20-30+ minutes, that are separated in time (ie, not acquired back to back). It scales each series by a constant factor, and with many b=0 volumes in each series, should be fairly reliable.

The DWI bias correction estimates a spatial multiplicative correction caused by B1 inhomogeneity. The bias field affects SIFT results. In mrtrix, this is applied after running topup and eddy

https://mrtrix.readthedocs.io/en/latest/fixel_based_analysis/mt_fibre_density_cross-section.html

QSIPrep will refuse to concatenate images with different phase encoding directions before denoising (which for us means unringing, mp-pca and b1 bias correction). Maybe --dwi-no-biascorr in addition to --no-b0-harmonization will make this less of a problem?

This seems different to what mrtrix does, and a bit worrying. If the bias field is estimated independently between series, it can provide a different correction both temporally and spatially. Subsequent b=0 harmonization could account for a global scaling between the series, but it won't affect spatial differences in the bias field. There's also the possibility that brain masking will be different between series, which will also change the bias field.

mattcieslak commented 1 year ago

Should we move the bias correction to after the resampling? Thinking more about it now, I'm leaning towards moving it. In addition to the other downsides you mention I wonder if it may also confuse uncorrected signal pileup and dropout for B1 inhomogeneity, which is obviously also bad.

cookpa commented 1 year ago

I wonder if it may also confuse uncorrected signal pileup and dropout for B1 inhomogeneity, which is obviously also bad.

It definitely will, this is a great point. For optimal performance, outliers should be winsorized. However, if the bias field is less than perfect, the impact of any errors will be mitigated if it is applied once to the whole data set.

PennLINC / qsiprep

NODDI recon producing weird values for ICVF and ISOVF #517