nipy / heudiconv

Flexible DICOM conversion into structured directory layouts
https://heudiconv.readthedocs.io
Other
234 stars 125 forks source link

Converting physiological recordings #532

Open tjhendrickson opened 2 years ago

tjhendrickson commented 2 years ago

Summary

I cannot seem to figure out how to convert physiological recordings present within the DICOM folder from a neuroimaging dataset. The physiological recordings are organized as one DICOM file within its own folder and with the same naming scheme as the BOLD it relates to with "PhysioLog" suffixed. The heuristic script below failed to find this and even when I run convertall.py as the heuristic the physiological recordings are not found.

Any ideas?

def create_key(template, outtype=('nii.gz','dicom'), annotation_classes=None): #), annotation_classes=None): if template is None or not template: raise ValueError('Template must be a valid format string') return (template, outtype, annotation_classes)

def infotodict(seqinfo): """Heuristic evaluator for determining which runs belong where

allowed template fields - follow python string module:

item: index within category
subject: participant id
seqitem: run number during scanning
subindex: sub index within group
""" 
t1 = create_key('sub-{subject}/{session}/anat/sub-{subject}_{session}_rec-{rec}_run-{item:02d}_T1w')
t2 = create_key('sub-{subject}/{session}/anat/sub-{subject}_{session}_rec-{rec}_run-{item:02d}_T2w')

rest_ten_minute = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-10minute_run-{item:02d}_part-{part}_bold')
rest_ten_minute_sbref = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-10minute_run-{item:02d}_sbref')
rest_ten_minute_physio = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-10minute_run-{item:02d}_physio')

rest_sixteen_minute = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-16minute_run-{item:02d}_part-{part}_bold')
rest_sixteen_minute_sbref = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-16minute_run-{item:02d}_sbref')
rest_sixteen_minute_physio = create_key('sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-16minute_run-{item:02d}_physio')

spinecho_fieldmap_ME_bold = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-SpinEchoME_dir-{dir}_run-{item:02d}_epi')
spinecho_fieldmap_ME_physio = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-SpinEchoME_dir-{dir}_run-{item:02d}_physio')

gradientecho_fieldmap_ME_bold_mag = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-GradientEchoME_dir-{dir}_run-{item:02d}_magnitude')
gradientecho_fieldmap_ME_bold_phase = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-GradientEchoME_dir-{dir}_run-{item:02d}_phase')
gradientecho_fieldmap_ME_bold_sbref = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-GradientEchoME_dir-{dir}_run-{item:02d}_sbref')
gradientecho_fieldmap_ME_bold_physio = create_key('sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-GradientEchoME_dir-{dir}_run-{item:02d}_physio')

info = {t1: [], t2: [], 
        rest_ten_minute: [], rest_ten_minute_sbref: [], rest_ten_minute_physio: [],
        rest_sixteen_minute: [], rest_sixteen_minute_sbref: [], rest_sixteen_minute_physio: [],
        spinecho_fieldmap_ME_bold: [], spinecho_fieldmap_ME_physio: [],
        gradientecho_fieldmap_ME_bold_mag: [], gradientecho_fieldmap_ME_bold_phase: [], gradientecho_fieldmap_ME_bold_sbref: [], gradientecho_fieldmap_ME_bold_physio: []}

for idx, s in enumerate(seqinfo):
    # retreive previous element in seqinfo
    if idx > 0:
        s_previous = seqinfo[idx-1]
    if idx - 1 > 0:
        s_previous_two = seqinfo[idx-2]
    # retreive next element in seqinfo
    if idx + 1 < len(seqinfo):
        s_next = seqinfo[idx+1]
    # retreive next next element in seqinfo
    if idx + 2 < len(seqinfo):
        s_next_two = seqinfo[idx+2]

    # find pre scan normalized anatomicals
    if (s.dim3 == 176) and ('NORM' in s.image_type):
        if 'ABCD_T1' in s.dcm_dir_name:
            rec = 'normalized'
            info[t1].append({'item': s.series_id, 'rec': rec})
        elif 'ABCD_T2' in s.dcm_dir_name:
            rec = 'normalized'
            info[t2].append({'item': s.series_id, 'rec': rec})
    # find resting state scans. Differentiate by mag or phase
    elif (s.dim4 > 5) and ('rest' in s.series_description):
        if '10MIN' in s.protocol_name:
            if s.image_type[2] == 'M':
                if (s_next.dim4 > 5) and ('rest' in s_next.series_description) and (s_next.image_type[2] == 'P'):
                    part = 'mag'
                    info[rest_ten_minute].append({'item': s.series_id,'part': part})
            elif s.image_type[2] == 'P':
                if (s_previous.dim4 > 5) and ('rest' in s_previous.series_description) and (s_previous.image_type[2] == 'M'):
                    part = 'phase'
                    info[rest_ten_minute].append({'item': s.series_id,'part': part})
        elif '16MIN' in s.protocol_name:
            if s.image_type[2] == 'M':
                if (s_next.dim4 > 5) and ('rest' in s_next.series_description) and (s_next.image_type[2] == 'P'):
                    part = 'mag'
                    info[rest_sixteen_minute].append({'item': s.series_id,'part': part})
            elif s.image_type[2] == 'P':
                if (s_previous.dim4 > 5) and ('rest' in s_previous.series_description) and (s_previous.image_type[2] == 'M'):
                    part = 'phase'
                    info[rest_sixteen_minute].append({'item': s.series_id,'part': part})
    # retreive field maps
    elif 'FieldMap' in s.series_description:
        if 'GEFieldMap' in s.series_description:  
            if 'AP' in s.series_description:
                if s.image_type[2] == 'M':
                    if (s_next.dim4 >= 15) and ('GEFieldMap' in s_next.series_description) and (s_next.image_type[2] == 'P'):
                        info[gradientecho_fieldmap_ME_bold_mag].append({'item': s.series_id, 'dir': 'AP'})
                elif s.image_type[2] == 'P':
                    if (s_previous.dim4 >= 15) and ('GEFieldMap' in s_previous.series_description) and (s_previous.image_type[2] == 'P'):
                        info[gradientecho_fieldmap_ME_bold_phase].append({'item': s.series_id, 'dir': 'AP'})
            elif 'PA' in s.series_description:
                if s.image_type[2] == 'M':
                    if (s_next.dim4 >= 15) and ('GEFieldMap' in s_next.series_description) and (s_next.image_type[2] == 'P'):
                        info[gradientecho_fieldmap_ME_bold_mag].append({'item': s.series_id, 'dir': 'PA'})
                elif s.image_type[2] == 'P':
                    if (s_previous.dim4 >= 15) and ('GEFieldMap' in s_previous.series_description) and (s_previous.image_type[2] == 'P'):
                        info[gradientecho_fieldmap_ME_bold_phase].append({'item': s.series_id, 'dir': 'PA'})
        elif 'SpinEchoFieldMap' in s.series_description:
            if s.dim4 == 15: 
                if 'AP' in s.series_description:
                    info[spinecho_fieldmap_ME_bold].append({'item': s.series_id, 'dir': 'AP'})
                elif 'PA' in s.series_description:
                    info[spinecho_fieldmap_ME_bold].append({'item': s.series_id, 'dir': 'PA'})

    # retreive sbref images
    elif 'SBRef' in s.series_description:
        if 'rest' in s.series_description:
            if '10MIN' in s.protocol_name:
                if (s_next.dim4 > 5) and ('10MIN' in s_next.series_description) \
                and (s_next.image_type[2] == 'M') and (s_next_two.dim4 > 5) \
                and ('10MIN' in s_next_two.series_description) \
                and (s_next_two.image_type[2] == 'P'):
                    info[rest_ten_minute_sbref].append({'item': s.series_id})
            elif '16MIN' in s.protocol_name:
                if (s_next.dim4 > 5) and ('16MIN' in s_next.series_description) \
                and (s_next.image_type[2] == 'M') and (s_next_two.dim4 > 5) \
                and ('16MIN' in s_next_two.series_description) \
                and (s_next_two.image_type[2] == 'P'):
                    info[rest_sixteen_minute_sbref].append({'item': s.series_id})
        elif 'GEFieldMap' in s.series_description:
            pdb.set_trace()
            if 'AP' in s.series_description:
                if s_next.image_type[2] == 'M' and (s_next.dim4 >= 15) and ('GEFieldMap' in s_next.series_description) and (s_next_two.dim4 >= 15) and ('GEFieldMap' in s_next_two.series_description) and (s_next_two.image_type[2] == 'P'):
                    info[gradientecho_fieldmap_ME_bold_sbref].append({'item': s.series_id, 'dir': 'AP'})
            elif 'PA' in s.series_description:
                if s_next.image_type[2] == 'M' and (s_next.dim4 >= 15) and ('GEFieldMap' in s_next.series_description) and (s_next_two.dim4 >= 15) and ('GEFieldMap' in s_next_two.series_description) and (s_next_two.image_type[2] == 'P'):
                    info[gradientecho_fieldmap_ME_bold_sbref].append({'item': s.series_id, 'dir': 'PA'})
    # retreive physiological recordings
    elif 'PhysioLog' in s.dcm_dir_name:
        if 'rest' in s.series_description:
            if '10MIN' in s.protocol_name:
                if (s_previous_two.dim4 > 5) and ('10MIN' in s_previous_two.series_description) \
                and (s_previous_two.image_type[2] == 'M') and (s_previous.dim4 > 5) \
                and ('10MIN' in s_previous.series_description) \
                and (s_previous.image_type[2] == 'P'):
                    info[rest_ten_minute_physio].append({'item': s.series_id})
            elif '16MIN' in s.protocol_name:
                if (s_previous_two.dim4 > 5) and ('16MIN' in s_previous_two.series_description) \
                and (s_previous_two.image_type[2] == 'M') and (s_previous.dim4 > 5) \
                and ('16MIN' in s_previous.series_description) \
                and (s_previous.image_type[2] == 'P'):
                    info[rest_sixteen_minute_physio].append({'item': s.series_id})
        elif 'GEFieldMap' in s.series_description:
            if 'AP' in s.series_description:
                if s_previous_two.image_type[2] == 'M' and (s_previous_two.dim4 >= 15) and ('GEFieldMap' in s_previous_two.series_description) and (s_previous.dim4 >= 15) and ('GEFieldMap' in s_previous.series_description) and (s_previous.image_type[2] == 'P'):
                    info[gradientecho_fieldmap_ME_bold_physio].append({'item': s.series_id, 'dir': 'AP'})
            elif 'PA' in s.series_description:
                if s_previous_two.image_type[2] == 'M' and (s_previous_two.dim4 >= 15) and ('GEFieldMap' in s_previous_two.series_description) and (s_previous.dim4 >= 15) and ('GEFieldMap' in s_previous.series_description) and (s_previous.image_type[2] == 'P'):
                    info[gradientecho_fieldmap_ME_bold_physio].append({'item': s.series_id, 'dir': 'PA'})

return info


### Platform details:

Choose one:
- [X ] Local environment
<!-- If selected, please provide OS and python version -->
python = 3.9.5
OS Ubuntu Xenial-20210429
- [ ] Container
<!-- If selected, please provide container name and tag"-->

- Heudiconv version: 
<!-- To check: run heudiconv with just the --version flag -->
heudiconv = 0.9.0
pvelasco commented 2 years ago

@tjhendrickson , could you show an example of the DICOM folder organization? Is it something like this?:

my_DICOMs \
          |- run_1 \
          |        |- run_1_volume_1.dcm
          |        |- run_1_volume_2.dcm
          |        | ...
          |- run_2 \
          |        |- run_2_volume_1.dcm
          |        |- run_2_volume_2.dcm
          |        | ...
          | ...
          |- run_1_PhysioLog \
          |        |- run_1_PhysioLog.dcm
          |- run_2_PhysioLog \
          |        |- run_2_PhysioLog.dcm
          | ...

Also, can you post the contents of the .heudiconv/<subID>/info/dicominfo.tsv and .heudiconv/<subID>/info/<subID>.auto.txt files (at the root of the output folder)?

tjhendrickson commented 2 years ago

@pvelasco Thanks for getting back to me. Yes, the DICOMs are organized in the way you suspected. Here is an example: `20211001-ST001-SUB1000101ses1/MR-SE030-rest_ep2d_bold_ME_16MIN_PhysioLog: MR-ST001-SE030-0001.dcm

20211001-ST001-SUB1000101ses1/MR-SE031-rest_ep2d_bold_ME_16MIN2_SBRef: MR-ST001-SE031-0001.dcm MR-ST001-SE031-EC3-0001.dcm MR-ST001-SE031-EC5-0001.dcm MR-ST001-SE031-EC2-0001.dcm MR-ST001-SE031-EC4-0001.dcm ` Attached are the files you asked for.

1000101_ses-1.auto.txt dicominfo_ses-1.txt

pvelasco commented 2 years ago

Actually, the physio data extraction has not been added yet to the master branch. PR #446 addresses it, but it is still under review. You can test that PR, to see if it works with your data (still, your heuristic.py file will need some changes, like adding an extra argument (outtype = ('physio',)) to all the create_key calls for physio files.

pvelasco commented 2 years ago

Looking at the dicominfo file, it didn't "see" series number 30 (the one with the PhysioLog DICOM). I'm not sure why. Do you mind posting the full heudiconv call? I assume these DICOMs with the physio in them come from the CMRR SMS sequence. Is that correct?

sajjadtorabian commented 2 years ago

I'm checking in to see if there has been any update on this. I guess since #446 is still open, the best way to convert physio is directly through bidsphysio, right? I'm thinking to run heudiconv on everything except physio, and then on top of it run bidsphysio