dHCP fMRI without SBREF or Fieldmap

[Johann_Drayne]

I have run the structural pipeline on our neonatal with relatively high success and would like to also use the fMRI pipeline. Unfortunately, I do not believe we have SBREF or spin-echo EPI images. I noticed a solution to having no Fieldmap and am hacking the script to either

1. Use the scripts motion correction methods to register the func to anat
2. Slot into the script our own motion corrected func in anat space

Or maybe it just be easier to write our own script to use the trained FIX classifier and not use the fMRI pipeline (unless below question 2) applies)?

My two main questions

1. Has anyone faced this issue of not having SBREF. Did you edit the script to go straight from func to anat. Or, did you skip the registration process and input your own motion corrected func to anat images.
2. I can get my hands on ‘DTI field maps’, can I insert these directly into the pipeline (in place of the spin-echo EPI images) or do I need to pre-process them in any way?

Thanking in advance for any advice offered!

[seanfitz]

Hi Johann

@Johann_Drayne

Has anyone faced this issue of not having SBREF. Did you edit the script to go straight from func to anat. Or, did you skip the registration process and input your own motion corrected func to anat images.

My advice, in the absence of SBREF, would be to take the temporal mean of the fMRI and call it your SBREF. The SBREF is just used an intermediary registration target in the func-to-T2 alignment because it has better tissue contrast than the func. The temporal mean of the fMRI would have better SNR than the func and potentially also have better contrast. This way the registration would go from the 1st-volume of the func, to the mean of the func, then to the T2w. It is important that the 1st volume of the func is used because the motion correction step (with EDDY) expects the field map to be aligned with the first volume.

@Johann_Drayne

I can get my hands on ‘DTI field maps’, can I insert these directly into the pipeline (in place of the spin-echo EPI images) or do I need to pre-process them in any way?

You can use the DTI field maps if there has not been a shim between the fMRI and the DTI, however keep in mind the field map will be less appropriate the further away in time it is from the fMRI acquisition. IHMO I think it would likely still be better than no field map.

In terms of pre-processing, the field map will need to be in rad/s as per FSL Fugue: FUGUE/Guide - FslWiki

Alternatively, if you have B0s with two different phase encode directions, you could use FSL topup to generate a field map for you. topup - FslWiki

I hope this helps.

Cheers, Sean

[Johann_Drayne]

Hi Sean,

Many thanks for the detailed reply! That is a good pointer/workaround in the absence of SBREF.

I have a couple more question about the DTI which I was hoping you could iron out?
*We unfortunately do not have B0s with two different phase encode directions so will have to stick with the DTI ‘method’.

In terms of shim settings we seem to have a range

• Small differences of: 0 -1 -1 -1 -1 1 6 1
• Average differences of: 1 5 -19 -82 9 12 8 9 && 7 10 -7 -60 40 -48 24 -15
• Large differences of: -85 174 -30 100 -423 -107 -315 396
  I assume the large differences are too large, but would the average distances be acceptable? Maybe it is something we just have to test out.

Furthermore, our DTI images are 3D. I assume it is adequate to create 4D images at TR intervals with the same number of time points as the fMRI images. But then our spinecho_pedir will not have alternating directions.
So maybe this method cannot work in our case?

On a final point, when running the pipeline I receive a

• Traceback (most recent call last): File "/home/johann.drayne/.local/bin/dhcp_neofmri.py", line 609, in <module> p.import_data(**args) File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 147, in import_data importdata.import_func( File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/util/util.py", line 51, in timed result = f(*args, **kw) TypeError: import_func() missing 2 required positional arguments: 'func_echospacing' and 'func_pedir'

• Traceback (most recent call last): File "/home/johann.drayne/.local/bin/dhcp_neofmri.py", line 623, in <module> p.run_all(**args) File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 930, in run_all self.prepare_fieldmap( File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 204, in prepare_fieldmap self.defaults.on_disk(['spinecho'], error=True) File "/home/johann.drayne/.local/lib/python3.8/site-packages/fsl/utils/filetree/filetree.py", line 363, in on_disk raise IOError("Failed to find any files on disk for {}".format(missing)) OSError: Failed to find any files on disk for ('spinecho',)

But I will post about these in a different topic, seeing as they are probably unrelated.

[seanfitz]

Hi Johann

Unfortunately I do not have the expertise to comment in more detail on the shim. I will ask one of my colleagues to take a look at this thread and chime in.

@Johann_Drayne

Furthermore, our DTI images are 3D. I assume it is adequate to create 4D images at TR intervals with the same number of time points as the fMRI images. But then our spinecho_pedir will not have alternating directions.
So maybe this method cannot work in our case?

By DTI image, I am assuming you are talking about a dual-echo-time fieldmap acquired along-side your diffusion weighted imaging? If so, this should be 3D. The pipeline will align it with the first volume of the fMRI and then it will be used as the initial estimate for the calculation of a 4D motion-by-susceptibility field.

@Johann_Drayne

On a final point, when running the pipeline I receive a

• Traceback (most recent call last): File "/home/johann.drayne/.local/bin/dhcp_neofmri.py", line 609, in <module> p.import_data(**args) File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 147, in import_data importdata.import_func( File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/util/util.py", line 51, in timed result = f(*args, **kw) TypeError: import_func() missing 2 required positional arguments: 'func_echospacing' and 'func_pedir'
• Traceback (most recent call last): File "/home/johann.drayne/.local/bin/dhcp_neofmri.py", line 623, in <module> p.run_all(**args) File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 930, in run_all self.prepare_fieldmap( File "/home/johann.drayne/.local/lib/python3.8/site-packages/dhcp/func/pipeline.py", line 204, in prepare_fieldmap self.defaults.on_disk(['spinecho'], error=True) File "/home/johann.drayne/.local/lib/python3.8/site-packages/fsl/utils/filetree/filetree.py", line 363, in on_disk raise IOError("Failed to find any files on disk for {}".format(missing)) OSError: Failed to find any files on disk for ('spinecho',)

But I will post about these in a different topic, seeing as they are probably unrelated.

Firstly, are you using the master branch? or the release v1.1 branch?

The master branch is for development, and it is entirely possible that low level changes might have broken high level code such as dhcp_neofmri.py. However, it should work in the release branch.

Irrespective, I don’t believe the high-level interface currently supports the dual-echo-time fieldmap. You would need to use the lower-level python interface. I could add an example to the examples folder early next week if that is helpful.

Cheers, Sean

[seanfitz]

Hi Johann

I have discussed the shim further with a colleague, and we think that, as you suggest, trial and error will be your best option.

I hope this helps.

Cheers, Sean

[Johann_Drayne]

Hi Sean,

I appreciate the swift replies!

That’s sounds good, I will trial and error the subjects with varying difference in shim settings. If I notice any ‘trend’ I will note my findings in this thread for future readers.

@seanfitz

By DTI image, I am assuming you are talking about a dual-echo-time fieldmap acquired along-side your diffusion weighted imaging? If so, this should be 3D. The pipeline will align it with the first volume of the fMRI and then it will be used as the initial estimate for the calculation of a 4D motion-by-susceptibility field.

Yes the image is 3D, however when I pass a 3D image into the pipeline (under --spinecho) I receive an error.
ValueError: incorrect shape for 4D nifti: (84, 84, 43)
The pipeline seems to prefer when it is 4D
But I could be doing something wrong.

I wrote my script based off this example
But maybe I will try this one in the meantime.

@seanfitz

Firstly, are you using the master branch? or the release v1.1 branch?

The master branch is for development, and it is entirely possible that low level changes might have broken high level code such as dhcp_neofmri.py. However, it should work in the release branch.

Switching to the v1.1 branch cleared up those previous errors. So thank-you for that advice! To note, I had to run with --no-deps flag (i.e. pip install --no-deps ...) and install the python dependencies after, as pprint was causing an issue.

@seanfitz

Irrespective, I don’t believe the high-level interface currently supports the dual-echo-time fieldmap. You would need to use the lower-level python interface. I could add an example to the examples folder early next week if that is helpful.

The latest error is that when the pipeline runs eddy, I get a returned non-zero exit status 1. This could be another silly issue I am making, so if you could write an example and if time permits, mock images as well (however just an example would be excellent) I would highly appreciate that! That would probably clear a lot of issues where I could just be making a silly mistake.

All the best,

Johann

[seanfitz]

Hi Johann

I have created an example script for running the pipeline (at a low-level) using a dual-echo-time-derived fieldmap in the master branch:

dhcp/examples/dhcp-dual-echo-time-fieldmap.py · master · Sean Fitzgibbon /...

Hopefully, this will help you to progress. I have added some comments highlighting the main differences to the default dHCP pipeline which uses a spin-echo-EPI-derived fieldmap.

@Johann_Drayne

Yes the image is 3D, however when I pass a 3D image into the pipeline (under --spinecho) I receive an error.
ValueError: incorrect shape for 4D nifti: (84, 84, 43)
The pipeline seems to prefer when it is 4D
But I could be doing something wrong.

So the problem here is that the --spinecho argument is intended for a spinecho EPI comprised of multiple volumes with opposing phase-encode directions. A fieldmap is then derived from this spinecho. However, as you already have a fieldmap, you can pass it directly which will hopefully be clearer in the example.

Unfortunately this functionality is only currently exposed at the lower python level, and not at the higher dhcp_neofmri.py command line level. I am in the process of updating dhcp_neofmri.py to support some of these alternative use cases, hence why it is currently broken on the master branch. I will attempt to raise the priority in getting it fixed.

Cheers, Sean

[Johann_Drayne]

Hi Sean,

Apologies for the late reply.

After installing the master branch, and running the example script you wrote, the pipeline runs well until eddy.

I made the wrong assumption that I could use eddy_openmp, we cannot use cuda as we have a Radeon PRO WX 2100. Running eddy_openmp with the same inputs as the pipeline, I get EddyInputError: Slice-to-vol is not yet implemented for the CPU version. It looks like we will need GPU support ( MRtrix3 issue similar FSL issue). So the only workaround is buy an NVIDIA GPU or modify the eddy command so that it can run on the CPU. I may look into how I can still use eddy but without slice-to-volume correction, however, that seemed like one of the key benefits; to correct for intra slice movement.
For the moment, I have passed the use_mcflirt = True at the start of the file. And in
mcdc.mcdc(... do_dc=False, do_s2v=False, do_mbs=False, use_mcflirt=use_mcflirt,)

Unfortunate I cannot use eddy for the time being, but at least I can use more of the pipeline.

The next issue I had was in FIX. When the FIX directory is set up, the mask.nii.gz seems to be created; 1. By taking the T2 mask
2. Warping it into fMRI space
3. Binarising it
4. Thresholding out lowest 5% of values.
Since the original T2 mask I provided was already binarised, after the thresholding, the image turned to all zeros.

I also had a hard time finding the trained data using the pipeline, so I manually called fix -f and fix -c, to iron out these issues I inelegantly:

1. Commented out line 136 in func/denoise.py
2. Commented out lines 146-149 in func/denoise.py
3. Commented out lines 304-306 in examples/ dhcp-dual-echo-time-fieldmap.py
4. Then I removed the mask, added another one from earlier in the pipeline and ran fix extract and classify with some ugly code which I added in at line 303 in examples/ dhcp-dual-echo-time-fieldmap.py

fixDir = f"/home/johann.drayne/fmridhcpnew/output/sub-{subid}/ses-session1/denoise/"
textinFile = f"/home/johann.drayne/fmridhcpnew/output/sub-{subid}/ses-session1/denoise/fix4melview_train35_thr5.txt "
textoutFile = f"/home/johann.drayne/fmridhcpnew/output/sub-{subid}/ses-session1/denoise/fix_labels.txt"

bashCmd1 = "rm -r /home/johann.drayne/fmridhcpnew/output/sub-MS040053/ses-session1/denoise/mask.nii.gz"
bashCmd2 = "cp /home/johann.drayne/fmridhcpnew/output/sub-MS040053/ses-session1/mcdc/func_mcdc_brainmask.nii.gz /home/johann.drayne/fmr>
bashCmd3 = "/usr/local/fix/fix -f " + fixDir
bashCmd4 = "/usr/local/fix/fix -c " + fixDir + " /home/johann.drayne/dhcpfmri/data/train35.RData 5"
bashCmd5 = "mv " + textinFile + textoutFile

process = subprocess.Popen(bashCmd1.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
print(output)
process = subprocess.Popen(bashCmd2.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
print(output)
process = subprocess.Popen(bashCmd3.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
print(output)
process = subprocess.Popen(bashCmd4.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
print(output)
process = subprocess.Popen(bashCmd5.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
print(output)

That ran FIX well, however there many of the ICA maps are unclassified so we are going to create our own training file and will update if there is much difference.

We had an issue in the STANDARD section with the pipeline not being able to find the atlas tree.

So we:

1. Commented out line 320 in examples/ dhcp-dual-echo-time-fieldmap.py
2. Added underneath atlas_tree = Path('/home/johann.drayne/dhcpfmri/data/dhcp_volumetric_atlas_extended/atlas/atlas.tree')

We again faced a similar issue in with QC

1. We did the same on line 650 in func/pipeline.py as above
2. Similarly on line 500 in func/pipeline.py we changed
   group_map = Path(util.get_setting('dhcp_group_maps'))
   →
   group_map = Path(‘/home/johann.drayne/dhcpfmri/data/group_maps.nii.gz’)

This brings me to the most recent error which I am struggling with.
The function I have gotten to in the main .py file is p.qc() on line 369.

The last function that the logger outputs is Add func_raw to QC standard this seems to be the function on line 702 in func/pipeline.py that is causing the issue. I have uploaded the last error output texterror.txt (6.2 KB)

It is true qc/raw_dr2.nii.gz is not created, however, there are a couple things that suggest something else is at play here.

1. The output qc/raw_dseg.nii.gz is almost a fully empty image (). To me this suggests there is something going on with the warp file reg/func-mcdc_to_standard/func-mcdc_to_standard_invwarp.nii.gz

   

2. The output qc/raw_standard_tsnr.nii.gz seems to have been warped, however it is weirdly clipped. Maybe suggesting the qc/raw_tsnr.nii.gz was not aligned with the warped file?

   

3. The output qc/tmpkmepy1sl/template_native.nii.gz seems to be a purely grey image. I think this is to do with the warp reg/func-mcdc_to_standard/func-mcdc_to_standard_invwarp.nii.gz again.

   

4. The output qc/raw_dr1.txt is a text file of nans’ but I think this is because of the previous reasons.
   

   

5. Because dr1.txt is full of nans’ this is probably stopping qc/raw_dr2.nii.gz from being created.

To conclude, do you know a fix to run eddy, but only the CPU version?
and
Do you have any ideas as to how to combat this last issue? I will create our own training file and will also try the pipeline with another subject and see if that brings around any differences.

I would like to thank you for all your previous help so far, this pipeline is very exciting!

All the best,

Johann

[seanfitz]

Hi Johann

@Johann_Drayne

do you know a fix to run eddy, but only the CPU version?

The current version of EDDY will only run on the GPU. However, I am reliably informed that a CPU version is currently in testing with the intention of being in the next FSL release (6.0.5). I am not sure of the exact ETA for the next FSL release but hopefully next month…

Alternatively the pipeline can use MCFLIRT (without slice-to-volume nor motion-by-susceptility) for motion correction, but then you would need to do a static distortion correction of the functional image with the sub-{subbed}/ses-{seaside}/reg/func_to_struct/func_dc_warp.nii.gz transform.

@Johann_Drayne

The next issue I had was in FIX. When the FIX directory is set up, the mask.nii.gz seems to be created; 1. By taking the T2 mask
2. Warping it into fMRI space
3. Binarising it
4. Thresholding out lowest 5% of values.
Since the original T2 mask I provided was already binarised, after the thresholding, the image turned to all zeros.

The T2 mask is expected to be binarised. After warping to the fMRI space it will no longer be binary because of interpolation, therefore it is thresholded at 0.5 and binarised. So your binary T2 mask should should still be binary in the fMRI space. Are you sure that the T2w mask is good? and did the registration of func-to-struct work well?

@Johann_Drayne

I also had a hard time finding the trained data using the pipeline
…
We had an issue in the STANDARD section with the pipeline not being able to find the atlas tree
…
We again faced a similar issue in with QC

The location of these resource files is defined in a settings file dhcp-neonatal-fmri-pipeline/dhcp/resources/settings.json. Did you use dhcp_fetch_resources.py to download these resources? if so, the settings.json file should have automatically been created.

You can create it manually. Mine looks like this:

{
    "dhcp_trained_fix": "/vols/Data/baby/seantemp/dhcp-neonatal-fmri-pipeline/dhcp/resources/train35.RData",
    "dhcp_trained_fix_threshold": 10,
    "dhcp_volumetric_atlas_tree": "/vols/Data/baby/seantemp/dhcp-neonatal-fmri-pipeline/dhcp/resources/dhcp_volumetric_atlas_extended/atlas/atlas.tree",
    "dhcp_group_maps": "/vols/Data/baby/seantemp/dhcp-neonatal-fmri-pipeline/dhcp/resources/group_maps.nii.gz",
    "dhcp_group_qc": "/vols/Data/baby/seantemp/dhcp-neonatal-fmri-pipeline/dhcp/resources/grp_qc_512_anon.json"
}

@Johann_Drayne

That ran FIX well, however there many of the ICA maps are unclassified so we are going to create our own training file and will update if there is much difference.

I would not be surprised if you ultimately had to retrain FIX, however I would wait until you had the upstream parts of the pipeline working such as slice-to-volume and motion-by-susceptibiity in case the differences created by their absence has impacted the classification performance.

@Johann_Drayne

The output qc/raw_dseg.nii.gz is almost a fully empty image (). To me this suggests there is something going on with the warp file reg/func-mcdc_to_standard/func-mcdc_to_standard_invwarp.nii.gz

I agree. The problem likely relates to registration, which might also explain your earlier issue with the FIX mask. Each registration directory contains a *_img.nii.gz file which is the source/moving image transformed to the ref/fixed space. I would check these to see if the registrations look ok.

Good luck.

Cheers, Sean

[Johann_Drayne]

Hi Sean,

@seanfitz

However, I am reliably informed that a CPU version is currently in testing with the intention of being in the next FSL release (6.0.5)

That is excellent news, very much looking forward to that.
It turns out I was able to connect to another lab computer that has a GeForce RTX 2060, Driver Version: 460.73.01 GPU. So we installed cuda 9.1 and ran the original pipeline. Eddy runs… however, after so many iterations we get EDDY::: EddyInternalGpuUtils::make_XtX_cuBLAS: cublasSsyrk error (similar issue in NVIDIA forum)
From what I can tell it seems that cuda has a minimum and maximum compute capability and it seems like our GPU would need cuda version 10.x or 11.x to work optimally (also talked about in this FSL forum).
Running eddy with a specific cuda version seems possible and if I successfully run it, I will comment my steps.
I mainly wanted to write the above in case anyone else was having issues with eddy.

The following points were noted using the use_mcflirt flag in the mcdc step

@seanfitz

Alternatively the pipeline can use MCFLIRT (without slice-to-volume nor motion-by-susceptility) for motion correction, but then you would need to do a static distortion correction of the functional image with the sub-{subbed}/ses-{seaside}/reg/func_to_struct/func_dc_warp.nii.gz transform.

I assume this is setting the
do_dc = True
do_s2v = True
do_mbs = True
at the start of the script? I just set these equal to False in the mcdc function.

@seanfitz

The T2 mask is expected to be binarised. After warping to the fMRI space it will no longer be binary because of interpolation, therefore it is thresholded at 0.5 and binarised. So your binary T2 mask should should still be binary in the fMRI space. Are you sure that the T2w mask is good? and did the registration of func-to-struct work well?

Ah ok, makes sense.
I think I found the culprit. I chose the wrong dseg file, from the anat dHCP file I am now using .
If the pipeline didn’t fully run
${dhcpanat}workdir/${subid}-${sesid}/segmentations/${subid}-${sesid}_tissue_labels

If the pipeline fully ran
${dhcpanat}derivatives/sub-${subid}/ses-${sesid}/anat/sub-${subid}_ses-${sesid}_drawem_tissue_labels

After this fix, the registration of func-to-struct looks good.

As a side note, I am using the restored T1 and T2 images from the anat dHCP pipeline. Should I be using the skull stripped version or the defaced versions, or are the restored versions (I assume N4 correction) of the images enough?

@seanfitz

The location of these resource files is defined in a settings file dhcp-neonatal-fmri-pipeline/dhcp/resources/settings.json . Did you use dhcp_fetch_resources.py to download these resources? if so, the settings.json file should have automatically been created.

You are completely right, what happened was I installed the pipeline. Downloaded the resources. Which created the file, but then deleting and re-installing to a different pipeline version, so the settings.sh file was empty.

So after a fresh install and subsequent re-downloading of the data (so the settings.json was created) I only had to make one edit to the scripts for the pipeline to run through all the stages (this is for the mcflirt version, instead of eddy).

During the FIX stage the fix_classify() function in func/denoise.py puts out the error.
mix = np.loadtxt(op.join(fixdir, 'filtered_func_data.ica/melodic_mix'))[:, noise_idx]
IndexError: arrays used as indices must be of integer (or boolean) type

So I changed line 199 to
noise_idx = (np.array(noise_idx)).astype(int)-1
This let the pipeline to run successfully. However, I am not sure if this is correct.

@seanfitz

I would not be surprised if you ultimately had to retrain FIX, however I would wait until you had the upstream parts of the pipeline working such as slice-to-volume and motion-by-susceptibiity in case the differences created by their absence has impacted the classification performance.

Unless we get eddy working and so FIX will have a better chance at classifying maps it looks like we will have to create our own training file.
What image would you recommend we use to create our training file?
Raw fMRI images that we input into the pipeline or perhaps the motion corrected fMRI images that are input at the ICA step?

As before, many thanks for your comments, they have been highly valued!

All the best,

Johann

[seanfitz]

Hi Johann

@Johann_Drayne

I assume this is setting the
do_dc = True
do_s2v = True
do_mbs = True
at the start of the script? I just set these equal to False in the mcdc function.

Yes that would work, alternatively could just set use_mcflirt=True

@Johann_Drayne

As a side note, I am using the restored T1 and T2 images from the anat dHCP pipeline. Should I be using the skull stripped version or the defaced versions, or are the restored versions (I assume N4 correction) of the images enough?

You should use the not-skull-stripped restored images. The pipeline will do skull stripping when and where it needs to.

@Johann_Drayne

So I changed line 199 to
noise_idx = (np.array(noise_idx)).astype(int)-1
This let the pipeline to run successfully. However, I am not sure if this is correct.

Yes that will work. Strange that they were no already integers… not sure why.

@Johann_Drayne

Unless we get eddy working and so FIX will have a better chance at classifying maps it looks like we will have to create our own training file.
What image would you recommend we use to create our training file?
Raw fMRI images that we input into the pipeline or perhaps the motion corrected fMRI images that are input at the ICA step?

The procedure for training fix would be:

1. Run the pipeline for 15-20 subjects, stopping prior to fix.fix_classify (line 304 in the example I sent earlier)
2. Manually inspect the ICA components and create a hand_labels_noise.txt file per subject which lists the artefact IC indices
3. Run FIX training as per, FIX/UserGuide - FslWiki
4. Run the remainder of the pipeline with the *.RData that you have generated in step 3.

If you do need to go down this path I can put together an example.

Cheers, Sean

[Johann_Drayne]

Hi Sean,

We got a new computer for the lab and decided to go with a V2 NVIDIA Quadro P1000 4Gb GPU. The exciting news is that eddy slice to volume correction will run on this GPU.
I have set,
do_dc = True
do_s2v = True
do_mbs = True
Looking at the eddy output I have a couple of questions that I wonder if you would have any insight.

This is the output from the eddy function.

_
_
The input functional image and fmap_to_func_img_Hz.nii.gz are in good alignment.
The below image is the fmap_to_func_img_Hz.nii.gz image.

_
_

The below image is the imported func.nii.gz image

_
_
I have tried looking on FSL forums as to why our output looks like the first image. But I cannot seem to find much information.
I ran the same eddy function in the pipeline while varying the parameters (namely, fwhm, mporder, s2v_lambdha, s2v_niter and mbs_niter), but did not find any differences.
I assume it is to do with our inputs.

Do you have any insight, would you require the full input files or had you any issue before and found there was no solution?

As always, I appreciate your input!

Many thanks, Johann

[seanfitz]

Hi Johann

So that clearly does not look right. It looks more like a temporal standard deviation image than the mcdc image. Can you clarify that the first “output” screenshot you present is the mcdc/func_mcdc.nii.gz file? and does it have the correct number of timepoints/volumes? and are you showing only the first timepoint/volume?

If so, it would be helpful for me to see the the actual EDDY call so we can clarify exactly which args were provided. The EDDY call should be the log files. You can get it quickly with:

cat logs/* | grep "eddy --imain"

It should look something like:

[2020-10-21 13:22:16,031 - dhcp.util.shellops.run ] INFO : /rds/general/user/sfitzgib/home/bin/eddy --imain=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/import/func.nii.gz --mask=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/import/func_brainmask.nii.gz --index=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/index.txt --bvals=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/bvals --bvecs=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/bvecs --acqp=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/eddy.acqp --out=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/func_mcdc --very_verbose --estimate_move_by_susceptibility --data_is_shelled --mbs_niter=20 --mbs_lambda=5 --mbs_ksp=5 --niter=10 --fwhm=10,10,5,5,0,0,0,0,0,0 --s2v_fwhm=0 --s2v_niter=10 --s2v_interp=trilinear --mporder=4 --nvoxhp=1000 --slspec=/rds/general/user/sfitzgib/home/dhcp-neonatal-fmri-pipeline/dhcp/resources/default_func.slorder --field=/rds/general/user/sfitzgib/home/workdir/sub-CC00108XX09/ses-36800/mcdc/fmap_to_func_img_Hz --b0_only --field_mat=/rds/general/user/sfitzgib/home/bin/fsl-6.0.3/etc/flirtsch/ident.mat --dont_mask_output --s2v_lambda=1

Cheers, Sean

[Johann_Drayne]

Hi Sean,

Yes, I can confirm that the “output” screenshot is the mcdc/func_mcdc.nii.gz image.
It also has the correct number of timepoints (compared to the imported func image).
Yes, I believe that was the first volume, however, the other volumes are also very similar.

Here is the eddy function

/usr/local/fsl/bin/eddy_cuda 
--imain=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/import/func.nii.gz 
--mask=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/import/func_brainmask.nii.gz
--index=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/index.txt
--bvals=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/bvals
--bvecs=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/bvecs
--acqp=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/eddy.acqp
--out=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/func_mcdc
--very_verbose
--estimate_move_by_susceptibility
--data_is_shelled
--mbs_niter=20
--mbs_lambda=5
--mbs_ksp=5
--niter=10
--fwhm=10,10,5,5,0,0,0,0,0,0
--s2v_fwhm=0
--s2v_niter=10
--s2v_interp=trilinear
--mporder=16
--nvoxhp=1000
--slspec=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/import/func.slorder
--field=/home/johann.drayne/dhcp_func/dhcp_func_output/v02/sub-MS040064/ses-session1/mcdc/fmap_to_func_img_Hz
--b0_only
--field_mat=/usr/local/fsl/etc/flirtsch/ident.mat
--dont_mask_output
--s2v_lambda=1

I should note that I changed line 111 in util/fsl.py
cmd = op.join(path, 'eddy') —> cmd = op.join(path, 'eddy_cuda’)
However, I assume this should not make that much of a difference.

Cheers,

Johann

[seanfitz]

Hi Johann

That looks sensible. The easiest way for me to debug this would be if I could play with the files directly. Is there any chance you could share the input files with me? If so, drop me an email (sean.fitzgibbon[at]ndcn.ox.ac.uk) and we can arrange a transfer.

Cheers, Sean

[seanfitz]

Hi Johann

The problem seems to be that there is a very tight brainmask on the fieldmap. Thus edge voxels in the fMRI are not covered by the fieldmap and that is messing EDDY up. I have modified your code to extrapolate the fieldmap in the fMRI phase-encode direction and it seems to work fine now. I will DM you the modified code but if useful to others the extrapolation was done with FSL fugue:

fmap = defaults.update(src_space='fmap', ref_space='func').get('resampled_image')
fmap_brainmask = defaults.update(src_space='fmap', ref_space='func').get('resampled_brainmask')

unwarpdir = get_axis(defaults.get('func'), func_pedir.name)[0]

fsl.fugue(
    loadfmap=fmap,
    mask=fmap_brainmask,
    unmaskfmap=True,
    savefmap=fmap,
    unwarpdir=unwarpdir,
)

The extrapolated fieldmap (in func space) looks like this:

And now the motion and distortion corrected fMRI looks like:

I hope this helps.

Cheers, Sean