Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: generate CIFTI derivatives #1001

Merged
merged 66 commits into from
Apr 20, 2018
Merged

ENH: generate CIFTI derivatives #1001

merged 66 commits into from
Apr 20, 2018

Conversation

mgxd
Copy link
Collaborator

@mgxd mgxd commented Feb 21, 2018

ref. #794

  • converts bold/gii data to CIFTI dtseries file
    • fsaverage5/fsaverage6 surfaces + MNI2009NLin2009cAsym subcortical volumes
  • adds compress input to DerivativesDataSink
  • increases required niworkflows to 0.3.7
TODOs
  • complete interface
  • add to bold workflow
  • add testing
To be added in another PR
  • create HCP 91282 greyordinate space CIFTIs (Support HCP-compatible CIFTI outputs #1025)
  • support fsnative surface + subject space aparc+aseg
    • number of indices in fsaverage5 (~10k) vs fsnative (+100k) will make these files fairly large

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks pretty reasonable from my limited knowledge of CIFTI. Could you add the --cifti-output flag to the ds005 tests?

https://github.com/poldracklab/fmriprep/blob/d13d39ef97a613c1901c0cc916834aca681f540e/.circleci/config.yml#L283-L289

subjects_dir = Directory(mandatory=True, desc="FreeSurfer SUBJECTS_DIR")
TR = traits.Float(mandatory=True, desc="repetition time")
gifti_files = traits.List(File(exists=True), mandatory=True,
desc="surface geometry files")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably worth specifying that this should be length 2 with order (left, right).

'label',
'*h.aparc.annot')))
# TODO: fetch label_file
return annotation_files
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it make sense to use a FreeSurferSource to find the aparc.annot files?

Also, do we want to consider the possibility people will want to use DKTatlas or a2009s (at least to the extent of being aware of whether we're making a move in that direction easy or hard)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could but we'd still have to parse the results - I can change if you feel strongly.

Could you elaborate on your second point? Use how?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FreeSurfer now automatically produces

% ls /data/out/ds000114/derivatives/freesurfer/sub-07/label/lh.aparc*.annot 
/data/out/ds000114/derivatives/freesurfer/sub-07/label/lh.aparc.a2009s.annot
/data/out/ds000114/derivatives/freesurfer/sub-07/label/lh.aparc.annot
/data/out/ds000114/derivatives/freesurfer/sub-07/label/lh.aparc.DKTatlas.annot

It seems plausible that someone might want to generate their CIFTI files with labels from another parcellation. Again, not necessarily that we need to code it up right now, but just to consider how easy it would be to augment.

(It may also be that I'm horribly wrong and CIFTI pretty much demands vanilla aparc.)

self.gifti_files,
self.volume_target,
self.surface_target,
self.TR)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self.inputs.*


def _fetch_data(self):
"""Converts inputspec to files"""
if self.surface_target == "fsnative":
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self.inputs.surface_target

raise NotImplementedError

annotation_files = sorted(glob(os.path.join(self.subjects_dir,
self.surface_target,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self.inputs.*

('targets.out', 'surface_target'),
('outputnode.surfaces', 'gifti_files')]),
(inputnode, cifti, [('subjects_dir', 'subjects_dir')]),
(bold_mni_trans_wf, outputnode, [('outputnode.bold_mni', 'bold_mni')]),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suspect this is meant to be:

(cifti, outputnode, [('out_file', 'bold_cifti')])

model_type = "CIFTI_MODEL_TYPE_SURFACE"
# use the corresponding annotation
hemi = structure.split('_')[-1][0]
annots = [fl for fl in annotation_files if ".aparc.annot" in fl]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While probably fine, this could be made a little more precise:

annots = [fl for fl in annotation_files if fl.endswith('.aparc.annot')]

else nb.freesurfer.read_annot(annots[1]))
# currently only supports L/R cortex
gii = (nb.load(gii_files[0]) if hemi == "LEFT"
else nb.load(gii_files[1]))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gii = nb.load(gii_files[hemi == 'RIGHT'])

output_spec = GenerateCiftiOutputSpec

def _run_interface(self, runtime):
self.annotation_files, self.label_file = self._fetch_data()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No clear reason to assign these to instance variables.

hemi = structure.split('_')[-1][0]
annots = [fl for fl in annotation_files if ".aparc.annot" in fl]
annot = (nb.freesurfer.read_annot(annots[0]) if hemi == "LEFT"
else nb.freesurfer.read_annot(annots[1]))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

annot = nb.freesurfer.read_annot(annots[hemi == 'RIGHT'])

@oesteban
Copy link
Member

we should probably store the label_file template to make this process faster - perhaps to the OSF project used in niworkflows?

👍 I'll give you access to the OSF repository. What is your username/orcid/handle?

@mgxd
Copy link
Collaborator Author

mgxd commented Feb 22, 2018

@oesteban orcid # 0000-0002-7252-7771

'CIFTI_STRUCTURE_CAUDATE_LEFT': [11],
'CIFTI_STRUCTURE_CAUDATE_RIGHT': [50],
'CIFTI_STRUCTURE_CEREBELLAR_WHITE_MATTER_LEFT': [7],
'CIFTI_STRUCTURE_CEREBELLAR_WHITE_MATTER_RIGHT': [46],
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i don't think we need the white matter voxels included.

'CIFTI_STRUCTURE_CEREBELLUM_LEFT': [6],
'CIFTI_STRUCTURE_CEREBELLUM_RIGHT': [45],
'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_LEFT': [28],
'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_RIGHT': [60],
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

or the diencephalon at present. i'm not sure these are labeled by the atlas.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I used the labels provided on mindboggle's page - I'll remove these though

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok you can leave these in place. do you have an example file on openmind that i can take a look at?

subjects_dir = Directory(mandatory=True, desc="FreeSurfer SUBJECTS_DIR")
TR = traits.Float(mandatory=True, desc="repetition time")
gifti_files = traits.List(File(exists=True), mandatory=True,
desc="list of surface geometry files (length 2 with order [L,R])")
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@effigies - does the gifti file contain any info about the surface target, in which case surface_target should be removed and just derived from the gifti file.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So these are functional time series, and it does not appear so.

% mris_info /data/out/ds000114/derivatives/fmriprep/sub-07/ses-test/func/sub-07_ses-test_task-fingerfootlips_bold_space-fsaverage5.L.func.gii
gifti_image struct
    version    = 1.0
    numDA      = 184
gim->meta nvpairs struct, len = 4 :
    nvpair: 'AnatomicalStructurePrimary' = 'CortexLeft'
    nvpair: 'UserName' = 'root'
    nvpair: 'Date' = 'Mon Oct 16 18:34:03 2017'
    nvpair: 'gifticlib-version' = 'gifti library version 1.09, 28 June, 2010'

gim->labeltable giiLabelTable struct, len = 0 :
--------------------------------------------------
gim->darray[0] giiDataArray struct
    intent   2001 = NIFTI_INTENT_TIME_SERIES
    datatype   16 = NIFTI_TYPE_FLOAT32
    ind_ord     1 = RowMajorOrder
    num_dim       = 1
    dims          = 10242, 0, 0, 0, 0, 0
    encoding    3 = GZipBase64Binary
    endian      2 = LittleEndian
    ext_fname     = 
    ext_offset    = 0
darray->meta nvpairs struct, len = 1 :
    nvpair: 'TimeStep' = '2500.000000'

darray->coordsys giiCoordSystem struct
    dataspace  = NIFTI_XFORM_UNKNOWN
    xformspace = NIFTI_XFORM_UNKNOWN
    xform[0] :  1.000000  0.000000  0.000000  0.000000
    xform[1] :  0.000000  1.000000  0.000000  0.000000
    xform[2] :  0.000000  0.000000  1.000000  0.000000
    xform[3] :  0.000000  0.000000  0.000000  1.000000
    data       = <set>
    nvals      = 10242
    nbyper     = 4
    numCS      = 1
darray->ex_atrs nvpairs struct, len = 0 :
--------------------------------------------------
--------------------------------------------------
gim->darray[1] giiDataArray struct
    intent   2001 = NIFTI_INTENT_TIME_SERIES
    datatype   16 = NIFTI_TYPE_FLOAT32
    ind_ord     1 = RowMajorOrder
    num_dim       = 1
    dims          = 10242, 0, 0, 0, 0, 0
    encoding    3 = GZipBase64Binary
    endian      2 = LittleEndian
    ext_fname     = 
    ext_offset    = 0
darray->meta nvpairs struct, len = 1 :
    nvpair: 'TimeStep' = '2500.000000'

darray->coordsys giiCoordSystem struct
    dataspace  = NIFTI_XFORM_UNKNOWN
    xformspace = NIFTI_XFORM_UNKNOWN
    xform[0] :  1.000000  0.000000  0.000000  0.000000
    xform[1] :  0.000000  1.000000  0.000000  0.000000
    xform[2] :  0.000000  0.000000  1.000000  0.000000
    xform[3] :  0.000000  0.000000  0.000000  1.000000
    data       = <set>
    nvals      = 10242
    nbyper     = 4
    numCS      = 1
darray->ex_atrs nvpairs struct, len = 0 :
--------------------------------------------------
[...]
--------------------------------------------------
gim->darray[183] giiDataArray struct
    intent   2001 = NIFTI_INTENT_TIME_SERIES
    datatype   16 = NIFTI_TYPE_FLOAT32
    ind_ord     1 = RowMajorOrder
    num_dim       = 1
    dims          = 10242, 0, 0, 0, 0, 0
    encoding    3 = GZipBase64Binary
    endian      2 = LittleEndian
    ext_fname     = 
    ext_offset    = 0
darray->meta nvpairs struct, len = 1 :
    nvpair: 'TimeStep' = '2500.000000'

darray->coordsys giiCoordSystem struct
    dataspace  = NIFTI_XFORM_UNKNOWN
    xformspace = NIFTI_XFORM_UNKNOWN
    xform[0] :  1.000000  0.000000  0.000000  0.000000
    xform[1] :  0.000000  1.000000  0.000000  0.000000
    xform[2] :  0.000000  0.000000  1.000000  0.000000
    xform[3] :  0.000000  0.000000  0.000000  1.000000
    data       = <set>
    nvals      = 10242
    nbyper     = 4
    numCS      = 1
darray->ex_atrs nvpairs struct, len = 0 :
--------------------------------------------------
gifti_image struct
    swapped    = 0
    compressed = 1
 -- darray totals: 7 MB
gim->ex_atrs nvpairs struct, len = 0 :
==================================================

bold_file = File(mandatory=True, exists=True, desc="input BOLD file")
volume_target = traits.Enum("MNI152NLin2009cAsym", mandatory=True, usedefault=True,
desc="CIFTI volumetric output space")
surface_target = traits.Enum("fsaverage5", "fsaverage6", "fsaverage", mandatory=True,
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

probably remove fsaverage for now as well since it's going to be really large.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems fine to me to leave, since fsaverage is a valid (though off-by-default) target for FMRIPREP.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i'm just worried about the file size - also the reason i suggested leaving out cerebellar wm voxels.

@mgxd
Copy link
Collaborator Author

mgxd commented Feb 27, 2018

sits on top of nipreps/niworkflows#228

@mgxd
Copy link
Collaborator Author

mgxd commented Feb 28, 2018

okay - this might be ready to go, but will require releasing and using current niworkflows to allow users to pull the label atlas

@effigies
Copy link
Member

effigies commented Feb 28, 2018

Oh, sorry. You'll need to pin the niworkflows commit. I can do that on your branch, as it's probably the easiest way to show you how.

Edit: Looks like you got the style fixes.

@satra
Copy link

satra commented Feb 28, 2018

@oesteban - while this is fine for now, i did ask @mgxd to put interfaces into nipype itself, so this may move over into misc.

@oesteban
Copy link
Member

I'll cut a new niworkflows release. When migration to nipype is done we'll release niworkflows again.

@oesteban
Copy link
Member

@mgxd niworkflows-0.3.5 has been released just now.

@effigies
Copy link
Member

effigies commented Apr 6, 2018

Sorry about not formatting as a single review. I'm done now.

@mgxd
Copy link
Collaborator Author

mgxd commented Apr 10, 2018

@effigies re: #1001 (comment), it looks like inputnode is never added to the workflow within init_sdc_wf if only epi fieldmap types are available, so the line https://github.com/poldracklab/fmriprep/blob/1755986606806207fdc665a2a469aff2cbfd2dba/fmriprep/workflows/bold/base.py#L387 should fail. I'll add a patch within this PR

@effigies
Copy link
Member

Oh, wow. Good catch.

@oesteban
Copy link
Member

the inputnode should be added in here:

https://github.com/poldracklab/fmriprep/blob/1755986606806207fdc665a2a469aff2cbfd2dba/fmriprep/workflows/fieldmap/base.py#L158

for the case no fieldmap is found. Maybe it is only a case of real fieldmaps?

@oesteban
Copy link
Member

Sorry, I misunderstood. You are right!

oesteban added a commit to oesteban/fmriprep that referenced this pull request Apr 10, 2018
As a result of refactoring, this node got disconnected for that particular use-case.

This PR fixes that (nipreps#1001 (comment))

Thanks @mgxd for spotting this problem.
@mgxd
Copy link
Collaborator Author

mgxd commented Apr 17, 2018

i think we're good to go 👍 !

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, one last iteration. Feel free to argue it out, but we feel fairly strongly that the resampling issue should be resolved before merging.

@satra Do you want to have a look at the BOLD->label resampling? I think you might have a little more CIFTI context and can make sure we haven't fooled ourselves into approaching this wrongly.

label_img = nb.load(label_file)
if not bold_img.shape[:3] == label_img.shape:
# we need these images to have the same resolution
bold_img = resample_to_img(bold_img, label_img, interpolation="nearest")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, we had a discussion over this line, and we think it should be fine if you just drop the interpolation='nearest'. The issue is that in many cases this is going to be upsampling BOLD data to 2mm^3 voxels, which should really be interpolated.


bold_img = nb.load(bold_file)
label_img = nb.load(label_file)
if not bold_img.shape[:3] == label_img.shape:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can drop this check. nilearn.image.resample_to_img will just return a copy if the shape and affine match.

So I'd actually write the whole loading section in two lines:

label_img = nb.load(label_file)
bold_img = resample_to_img(bold_file, label_img)

@effigies
Copy link
Member

Awesome. In it goes.

@effigies effigies merged commit 677c556 into nipreps:master Apr 20, 2018
1.3.0 automation moved this from In progress to Done Apr 20, 2018
@chrisgorgo
Copy link
Contributor

This is great! Thank you so much @mgxd and everyone else for the feedback.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
No open projects
1.3.0
  
Done
Development

Successfully merging this pull request may close these issues.

None yet

5 participants