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

Conversation

Projects
5 participants
@mgxd
Copy link
Collaborator

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 (#1025)
  • support fsnative surface + subject space aparc+aseg
    • number of indices in fsaverage5 (~10k) vs fsnative (+100k) will make these files fairly large

mgxd added some commits Feb 21, 2018

@effigies
Copy link
Collaborator

effigies left a comment

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

fmriprep-docker -i poldracklab/fmriprep:latest \
--config $PWD/nipype.cfg -w /tmp/ds005/work \
/tmp/data/ds005 /tmp/ds005/derivatives participant \
--debug --write-graph --use-syn-sdc --use-aroma \
--ignore-aroma-denoising-errors --mem_mb 4096 \
--output-space T1w template fsaverage5 \
--nthreads 2 -vv

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")

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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

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

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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)?

This comment has been minimized.

@mgxd

mgxd Feb 22, 2018

Author Collaborator

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?

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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)

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

self.inputs.*


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

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

self.inputs.surface_target

raise NotImplementedError

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

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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')]),

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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]

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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]))

This comment has been minimized.

@effigies

effigies Feb 22, 2018

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

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

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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]))

This comment has been minimized.

@effigies

effigies Feb 22, 2018

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

This comment has been minimized.

Copy link
Contributor

oesteban commented Feb 22, 2018

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

This comment has been minimized.

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],

This comment has been minimized.

@satra

satra Feb 22, 2018

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],

This comment has been minimized.

@satra

satra Feb 22, 2018

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

This comment has been minimized.

@mgxd

mgxd Feb 23, 2018

Author Collaborator

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

This comment has been minimized.

@satra

satra Feb 23, 2018

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])")

This comment has been minimized.

@satra

satra Feb 22, 2018

@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.

This comment has been minimized.

@effigies

effigies Feb 22, 2018

Collaborator

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,

This comment has been minimized.

@satra

satra Feb 22, 2018

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

This comment has been minimized.

@effigies

effigies Feb 23, 2018

Collaborator

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

This comment has been minimized.

@satra

satra Feb 23, 2018

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

mgxd added some commits Feb 23, 2018

enh+fix: save output as valid CIFTI dtseries
remove unused CIFTI structures, add link for OSF label storage
@mgxd

This comment has been minimized.

Copy link
Collaborator Author

mgxd commented Feb 27, 2018

@mgxd

This comment has been minimized.

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

This comment has been minimized.

Copy link
Collaborator

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

This comment has been minimized.

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.

effigies and others added some commits Feb 28, 2018

@oesteban

This comment has been minimized.

Copy link
Contributor

oesteban commented Feb 28, 2018

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

@oesteban

This comment has been minimized.

Copy link
Contributor

oesteban commented Feb 28, 2018

@mgxd niworkflows-0.3.5 has been released just now.

@mgxd

This comment has been minimized.

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

bold_sdc_wf.inputs.inputnode.template = template
should fail. I'll add a patch within this PR

@effigies

This comment has been minimized.

Copy link
Collaborator

effigies commented Apr 10, 2018

Oh, wow. Good catch.

@oesteban

This comment has been minimized.

Copy link
Contributor

oesteban commented Apr 10, 2018

the inputnode should be added in here:

(inputnode, outputnode, [('bold_ref', 'bold_ref'),

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

@oesteban

This comment has been minimized.

Copy link
Contributor

oesteban commented Apr 10, 2018

Sorry, I misunderstood. You are right!

oesteban added a commit to oesteban/fmriprep that referenced this pull request Apr 10, 2018

[FIX] Connect inputnode to SDC for pepolar images
As a result of refactoring, this node got disconnected for that particular use-case.

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

Thanks @mgxd for spotting this problem.
@mgxd

This comment has been minimized.

Copy link
Collaborator Author

mgxd commented Apr 17, 2018

i think we're good to go 👍 !

mgxd added some commits Apr 19, 2018

@effigies
Copy link
Collaborator

effigies left a comment

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")

This comment has been minimized.

@effigies

effigies Apr 20, 2018

Collaborator

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:

This comment has been minimized.

@effigies

effigies Apr 20, 2018

Collaborator

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

This comment has been minimized.

Copy link
Collaborator

effigies commented Apr 20, 2018

Awesome. In it goes.

@effigies effigies merged commit 677c556 into poldracklab:master Apr 20, 2018

9 checks passed

ci/circleci: build Your tests passed on CircleCI!
Details
ci/circleci: build_docs Your tests passed on CircleCI!
Details
ci/circleci: ds005 Your tests passed on CircleCI!
Details
ci/circleci: ds054 Your tests passed on CircleCI!
Details
ci/circleci: ds210 Your tests passed on CircleCI!
Details
ci/circleci: get_data Your tests passed on CircleCI!
Details
ci/circleci: test_pytest Your tests passed on CircleCI!
Details
ci/circleci: update_cache Your tests passed on CircleCI!
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

1.3.0 automation moved this from In progress to Done Apr 20, 2018

@chrisfilo

This comment has been minimized.

Copy link
Contributor

chrisfilo commented Apr 21, 2018

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