From 12972ee5e2052a4d1dedfe3201899d9e34a8a9ec Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Mon, 11 Feb 2013 17:44:34 -0800 Subject: [PATCH 01/12] NF: make SPM processing example for Open FMRI ds107 Example of how to batch script SPM5 (!) preprocessing for the Open FMRI ds107 dataset. --- examples/interfaces/process_ds107.py | 234 +++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100755 examples/interfaces/process_ds107.py diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py new file mode 100755 index 0000000000..76ccf39851 --- /dev/null +++ b/examples/interfaces/process_ds107.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +''' Single subject analysis script for SPM / Open FMRI ds107 ''' +import sys +from os.path import join as pjoin, abspath, splitext, isfile +from glob import glob +from warnings import warn +import gzip + +import numpy as np + +import nipy.interfaces.matlab as nimat +from nipy.interfaces.spm import (spm_info, make_job, scans_for_fnames, + run_jobdef, fnames_presuffix, fname_presuffix, + fltcols) + + +# The batch scripts currently need SPM5 +nimat.matlab_cmd = 'matlab-2007a-spm5 -nodesktop -nosplash' + +TR = 3.0 +NSLICES = 37 +SLICE_ORDER = range(1, NSLICES, 2) + range(2, NSLICES, 2) + +def get_data(data_path, subj_id): + data_path = abspath(data_path) + data_def = {} + subject_path = pjoin(data_path, 'sub%03d' % subj_id) + functionals = sorted( + glob(pjoin(subject_path, 'BOLD', 'task*', 'bold*.nii.gz'))) + anatomicals = sorted( + glob(pjoin(subject_path, 'anatomy', 'highres001.nii.gz'))) + assert len(functionals) == 2 + assert len(anatomicals) == 1 + for flist in (anatomicals, functionals): + for i, fname in enumerate(flist): + nogz, _ = splitext(fname) + if not isfile(nogz): + contents = gzip.open(fname, 'rb').read() + with open(nogz, 'wb') as fobj: + fobj.write(contents) + if not isfile(nogz): + raise RuntimeError('No nii for ' + fname) + flist[i] = nogz + data_def['anatomical'] = anatomicals[0] + data_def['functionals'] = functionals + return data_def + + +def default_ta(tr, nslices): + slice_time = tr / float(nslices) + return slice_time * (nslices - 1) + + +def slicetime(data_def): + sess_scans = scans_for_fnames(data_def['functionals']) + stinfo = make_job('temporal', 'st', { + 'scans': sess_scans, + 'so':SLICE_ORDER, + 'tr':TR, + 'ta':default_ta(TR, NSLICES), + 'nslices':float(NSLICES), + 'refslice':1 + }) + run_jobdef(stinfo) + + +def realign(data_def): + sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) + rinfo = make_job('spatial', 'realign', [{ + 'estimate':{ + 'data':sess_scans, + 'eoptions':{ + 'quality':0.9, + 'sep':4.0, + 'fwhm':5.0, + 'rtm':True, + 'interp':2.0, + 'wrap':[0.0,0.0,0.0], + 'weight':[] + } + } + }]) + run_jobdef(rinfo) + + +def reslice(data_def): + sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) + rsinfo = make_job('spatial', 'realign', [{ + 'write':{ + 'data': np.vstack(sess_scans.flat), + 'roptions':{ + 'which':[2, 1], + 'interp':4.0, + 'wrap':[0.0,0.0,0.0], + 'mask':True, + } + } + }]) + run_jobdef(rsinfo) + + +def coregister(data_def): + func1 = data_def['functionals'][0] + mean_fname = fname_presuffix(func1, 'meana') + crinfo = make_job('spatial', 'coreg', [{ + 'estimate':{ + 'ref': [mean_fname], + 'source': [data_def['anatomical']], + 'other': [[]], + 'eoptions':{ + 'cost_fun':'nmi', + 'sep':[4.0, 2.0], + 'tol':np.array( + [0.02,0.02,0.02, + 0.001,0.001,0.001, + 0.01,0.01,0.01, + 0.001,0.001,0.001]).reshape(1,12), + 'fwhm':[7.0, 7.0] + } + } + }]) + run_jobdef(crinfo) + + +def segnorm(data_def): + def_tpms = np.zeros((3,1), dtype=np.object) + spm_path = spm_info.spm_path + def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), + def_tpms[1] = pjoin(spm_path, 'tpm', 'white.nii'), + def_tpms[2] = pjoin(spm_path, 'tpm', 'csf.nii') + data = np.zeros((1,), dtype=object) + data[0] = data_def['anatomical'] + sninfo = make_job('spatial', 'preproc', { + 'data': data, + 'output':{ + 'GM':fltcols([0,0,1]), + 'WM':fltcols([0,0,1]), + 'CSF':fltcols([0,0,0]), + 'biascor':1.0, + 'cleanup':False, + }, + 'opts':{ + 'tpm':def_tpms, + 'ngaus':fltcols([2,2,2,4]), + 'regtype':'mni', + 'warpreg':1.0, + 'warpco':25.0, + 'biasreg':0.0001, + 'biasfwhm':60.0, + 'samp':3.0, + 'msk':np.array([], dtype=object), + } + }) + run_jobdef(sninfo) + + +def norm_write(data_def): + sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) + matname = fname_presuffix(data_def['anatomical'], + suffix='_seg_sn.mat', + use_ext=False) + subj = { + 'matname': np.zeros((1,), dtype=object), + 'resample': np.vstack(sess_scans.flat), + } + subj['matname'][0] = matname + roptions = { + 'preserve':False, + 'bb':np.array([[-78,-112, -50],[78,76,85.0]]), + 'vox':fltcols([2.0,2.0,2.0]), + 'interp':1.0, + 'wrap':[0.0,0.0,0.0], + } + nwinfo = make_job('spatial', 'normalise', [{ + 'write':{ + 'subj': subj, + 'roptions': roptions, + } + }]) + run_jobdef(nwinfo) + # knock out the list of images, replacing with only one + subj['resample'] = np.zeros((1,), dtype=object) + subj['resample'][0] = data_def['anatomical'] + roptions['interp'] = 4.0 + run_jobdef(nwinfo) + + +def smooth(data_def, fwhm=8.0): + try: + len(fwhm) + except TypeError: + fwhm = [fwhm] * 3 + fwhm = np.asarray(fwhm, dtype=np.float).reshape(1,3) + sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'wa')) + sinfo = make_job('spatial', 'smooth', + {'data':np.vstack(sess_scans.flat), + 'fwhm':fwhm, + 'dtype':0}) + run_jobdef(sinfo) + + +def process_subject(ddef): + """ Process subject from subject data dict `ddef` + """ + if not ddef['anatomical']: + warn("No anatomical, aborting processing") + return + slicetime(ddef) + realign(ddef) + reslice(ddef) + coregister(ddef) + segnorm(ddef) + norm_write(ddef) + smooth(ddef) + + +def process_subjects(data_path, subj_ids): + for subj_id in subj_ids: + ddef = get_data(data_path, subj_id) + process_subject(ddef) + + +if __name__ == '__main__': + try: + data_path = sys.argv[1] + except IndexError: + raise OSError('Need ds107 data path as input') + if len(sys.argv) > 2: + subj_ids = [int(id) for id in sys.argv[2:]] + else: + subj_ids = range(1, 16) + process_subjects(data_path, subj_ids) From 31afdbb43712ff72d4edfee114811d060f73049b Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 15:59:35 -0700 Subject: [PATCH 02/12] RF: refactor SPM batch script for study params Allow passing of study and analysis parameter dictionaries. --- examples/interfaces/process_ds107.py | 346 +++++++++++++++------------ 1 file changed, 187 insertions(+), 159 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index 76ccf39851..8a86206f62 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -3,6 +3,7 @@ # vi: set ft=python sts=4 ts=4 sw=4 et: ''' Single subject analysis script for SPM / Open FMRI ds107 ''' import sys +from copy import deepcopy from os.path import join as pjoin, abspath, splitext, isfile from glob import glob from warnings import warn @@ -19,9 +20,13 @@ # The batch scripts currently need SPM5 nimat.matlab_cmd = 'matlab-2007a-spm5 -nodesktop -nosplash' -TR = 3.0 -NSLICES = 37 -SLICE_ORDER = range(1, NSLICES, 2) + range(2, NSLICES, 2) +N_SLICES = 37 +STUDY_DEF = dict( + TR = 3.0, + n_slices = N_SLICES, + time_to_space = range(1, N_SLICES, 2) + range(2, N_SLICES, 2) +) + def get_data(data_path, subj_id): data_path = abspath(data_path) @@ -53,173 +58,196 @@ def default_ta(tr, nslices): return slice_time * (nslices - 1) -def slicetime(data_def): - sess_scans = scans_for_fnames(data_def['functionals']) - stinfo = make_job('temporal', 'st', { - 'scans': sess_scans, - 'so':SLICE_ORDER, - 'tr':TR, - 'ta':default_ta(TR, NSLICES), - 'nslices':float(NSLICES), - 'refslice':1 - }) - run_jobdef(stinfo) - - -def realign(data_def): - sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) - rinfo = make_job('spatial', 'realign', [{ - 'estimate':{ - 'data':sess_scans, - 'eoptions':{ - 'quality':0.9, - 'sep':4.0, - 'fwhm':5.0, - 'rtm':True, - 'interp':2.0, - 'wrap':[0.0,0.0,0.0], - 'weight':[] +class SPMSubjectAnalysis(object): + """ Class to preprocess single subject in SPM + """ + def __init__(self, data_def, study_def, ana_def): + self.data_def = deepcopy(data_def) + self.study_def = self.add_study_defaults(study_def) + self.ana_def = self.add_ana_defaults(deepcopy(ana_def)) + + def add_study_defaults(self, study_def): + full_study_def = deepcopy(study_def) + if 'TA' not in full_study_def: + full_study_def['TA'] = default_ta( + full_study_def['TR'], full_study_def['n_slices']) + return full_study_def + + def add_ana_defaults(self, ana_def): + full_ana_def = deepcopy(ana_def) + if 'fwhm' not in full_ana_def: + full_ana_def['fwhm'] = 8.0 + return full_ana_def + + def slicetime(self): + sess_scans = scans_for_fnames(self.data_def['functionals']) + sdef = self.study_def + stinfo = make_job('temporal', 'st', { + 'scans': sess_scans, + 'so': sdef['time_to_space'], + 'tr': sdef['TR'], + 'ta': sdef['TA'], + 'nslices': float(sdef['n_slices']), + 'refslice':1 + }) + run_jobdef(stinfo) + + + def realign(self): + sess_scans = scans_for_fnames( + fnames_presuffix(self.data_def['functionals'], 'a')) + rinfo = make_job('spatial', 'realign', [{ + 'estimate':{ + 'data':sess_scans, + 'eoptions':{ + 'quality': 0.9, + 'sep': 4.0, + 'fwhm': 5.0, + 'rtm': True, + 'interp': 2.0, + 'wrap': [0.0,0.0,0.0], + 'weight': [] + } } - } - }]) - run_jobdef(rinfo) - - -def reslice(data_def): - sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) - rsinfo = make_job('spatial', 'realign', [{ - 'write':{ - 'data': np.vstack(sess_scans.flat), - 'roptions':{ - 'which':[2, 1], - 'interp':4.0, - 'wrap':[0.0,0.0,0.0], - 'mask':True, + }]) + run_jobdef(rinfo) + + def reslice(self): + sess_scans = scans_for_fnames( + fnames_presuffix(self.data_def['functionals'], 'a')) + rsinfo = make_job('spatial', 'realign', [{ + 'write':{ + 'data': np.vstack(sess_scans.flat), + 'roptions':{ + 'which':[2, 1], + 'interp':4.0, + 'wrap':[0.0,0.0,0.0], + 'mask':True, + } } - } - }]) - run_jobdef(rsinfo) - - -def coregister(data_def): - func1 = data_def['functionals'][0] - mean_fname = fname_presuffix(func1, 'meana') - crinfo = make_job('spatial', 'coreg', [{ - 'estimate':{ - 'ref': [mean_fname], - 'source': [data_def['anatomical']], - 'other': [[]], - 'eoptions':{ - 'cost_fun':'nmi', - 'sep':[4.0, 2.0], - 'tol':np.array( - [0.02,0.02,0.02, - 0.001,0.001,0.001, - 0.01,0.01,0.01, - 0.001,0.001,0.001]).reshape(1,12), - 'fwhm':[7.0, 7.0] + }]) + run_jobdef(rsinfo) + + def coregister(self): + func1 = self.data_def['functionals'][0] + mean_fname = fname_presuffix(func1, 'meana') + crinfo = make_job('spatial', 'coreg', [{ + 'estimate':{ + 'ref': [mean_fname], + 'source': [self.data_def['anatomical']], + 'other': [[]], + 'eoptions':{ + 'cost_fun':'nmi', + 'sep':[4.0, 2.0], + 'tol':np.array( + [0.02,0.02,0.02, + 0.001,0.001,0.001, + 0.01,0.01,0.01, + 0.001,0.001,0.001]).reshape(1,12), + 'fwhm':[7.0, 7.0] + } } - } - }]) - run_jobdef(crinfo) - - -def segnorm(data_def): - def_tpms = np.zeros((3,1), dtype=np.object) - spm_path = spm_info.spm_path - def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), - def_tpms[1] = pjoin(spm_path, 'tpm', 'white.nii'), - def_tpms[2] = pjoin(spm_path, 'tpm', 'csf.nii') - data = np.zeros((1,), dtype=object) - data[0] = data_def['anatomical'] - sninfo = make_job('spatial', 'preproc', { - 'data': data, - 'output':{ - 'GM':fltcols([0,0,1]), - 'WM':fltcols([0,0,1]), - 'CSF':fltcols([0,0,0]), - 'biascor':1.0, - 'cleanup':False, - }, - 'opts':{ - 'tpm':def_tpms, - 'ngaus':fltcols([2,2,2,4]), - 'regtype':'mni', - 'warpreg':1.0, - 'warpco':25.0, - 'biasreg':0.0001, - 'biasfwhm':60.0, - 'samp':3.0, - 'msk':np.array([], dtype=object), - } - }) - run_jobdef(sninfo) - - -def norm_write(data_def): - sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'a')) - matname = fname_presuffix(data_def['anatomical'], - suffix='_seg_sn.mat', - use_ext=False) - subj = { - 'matname': np.zeros((1,), dtype=object), - 'resample': np.vstack(sess_scans.flat), - } - subj['matname'][0] = matname - roptions = { - 'preserve':False, - 'bb':np.array([[-78,-112, -50],[78,76,85.0]]), - 'vox':fltcols([2.0,2.0,2.0]), - 'interp':1.0, - 'wrap':[0.0,0.0,0.0], - } - nwinfo = make_job('spatial', 'normalise', [{ - 'write':{ - 'subj': subj, - 'roptions': roptions, - } - }]) - run_jobdef(nwinfo) - # knock out the list of images, replacing with only one - subj['resample'] = np.zeros((1,), dtype=object) - subj['resample'][0] = data_def['anatomical'] - roptions['interp'] = 4.0 - run_jobdef(nwinfo) - - -def smooth(data_def, fwhm=8.0): - try: - len(fwhm) - except TypeError: - fwhm = [fwhm] * 3 - fwhm = np.asarray(fwhm, dtype=np.float).reshape(1,3) - sess_scans = scans_for_fnames(fnames_presuffix(data_def['functionals'], 'wa')) - sinfo = make_job('spatial', 'smooth', - {'data':np.vstack(sess_scans.flat), - 'fwhm':fwhm, - 'dtype':0}) - run_jobdef(sinfo) - - -def process_subject(ddef): + }]) + run_jobdef(crinfo) + + def segnorm(self): + def_tpms = np.zeros((3,1), dtype=np.object) + spm_path = spm_info.spm_path + def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), + def_tpms[1] = pjoin(spm_path, 'tpm', 'white.nii'), + def_tpms[2] = pjoin(spm_path, 'tpm', 'csf.nii') + data = np.zeros((1,), dtype=object) + data[0] = self.data_def['anatomical'] + sninfo = make_job('spatial', 'preproc', { + 'data': data, + 'output':{ + 'GM':fltcols([0,0,1]), + 'WM':fltcols([0,0,1]), + 'CSF':fltcols([0,0,0]), + 'biascor':1.0, + 'cleanup':False, + }, + 'opts':{ + 'tpm':def_tpms, + 'ngaus':fltcols([2,2,2,4]), + 'regtype':'mni', + 'warpreg':1.0, + 'warpco':25.0, + 'biasreg':0.0001, + 'biasfwhm':60.0, + 'samp':3.0, + 'msk':np.array([], dtype=object), + } + }) + run_jobdef(sninfo) + + def norm_write(self): + sess_scans = scans_for_fnames( + fnames_presuffix(self.data_def['functionals'], 'a')) + matname = fname_presuffix(self.data_def['anatomical'], + suffix='_seg_sn.mat', + use_ext=False) + subj = { + 'matname': np.zeros((1,), dtype=object), + 'resample': np.vstack(sess_scans.flat), + } + subj['matname'][0] = matname + roptions = { + 'preserve':False, + 'bb':np.array([[-78,-112, -50],[78,76,85.0]]), + 'vox':fltcols([2.0,2.0,2.0]), + 'interp':1.0, + 'wrap':[0.0,0.0,0.0], + } + nwinfo = make_job('spatial', 'normalise', [{ + 'write':{ + 'subj': subj, + 'roptions': roptions, + } + }]) + run_jobdef(nwinfo) + # knock out the list of images, replacing with only one + subj['resample'] = np.zeros((1,), dtype=object) + subj['resample'][0] = self.data_def['anatomical'] + roptions['interp'] = 4.0 + run_jobdef(nwinfo) + + def smooth(self): + fwhm = self.ana_def['fwhm'] + try: + len(fwhm) + except TypeError: + fwhm = [fwhm] * 3 + fwhm = np.asarray(fwhm, dtype=np.float).reshape(1,3) + sess_scans = scans_for_fnames( + fnames_presuffix(self.data_def['functionals'], 'wa')) + sinfo = make_job('spatial', 'smooth', + {'data':np.vstack(sess_scans.flat), + 'fwhm':fwhm, + 'dtype':0}) + run_jobdef(sinfo) + + +def process_subject(ddef, study_def, ana_def): """ Process subject from subject data dict `ddef` """ if not ddef['anatomical']: warn("No anatomical, aborting processing") return - slicetime(ddef) - realign(ddef) - reslice(ddef) - coregister(ddef) - segnorm(ddef) - norm_write(ddef) - smooth(ddef) + ana = SPMSubjectAnalysis(ddef, study_def, ana_def) + ana.slicetime() + ana.realign() + ana.reslice() + ana.coregister() + ana.segnorm() + ana.norm_def() + ana.smooth() -def process_subjects(data_path, subj_ids): +def process_subjects(data_path, subj_ids, study_def, ana_def): for subj_id in subj_ids: ddef = get_data(data_path, subj_id) - process_subject(ddef) + process_subject(ddef, study_def, ana_def) if __name__ == '__main__': @@ -231,4 +259,4 @@ def process_subjects(data_path, subj_ids): subj_ids = [int(id) for id in sys.argv[2:]] else: subj_ids = range(1, 16) - process_subjects(data_path, subj_ids) + process_subjects(data_path, subj_ids, STUDY_DEF, {}) From e9715930e63db2521898a64e36a001fc44dd8d0f Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 16:10:13 -0700 Subject: [PATCH 03/12] RF: pass through if is non-gzipped nifti already --- examples/interfaces/process_ds107.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index 8a86206f62..336fcdbb6e 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -40,14 +40,13 @@ def get_data(data_path, subj_id): assert len(anatomicals) == 1 for flist in (anatomicals, functionals): for i, fname in enumerate(flist): - nogz, _ = splitext(fname) - if not isfile(nogz): - contents = gzip.open(fname, 'rb').read() - with open(nogz, 'wb') as fobj: - fobj.write(contents) + nogz, gz_ext = splitext(fname) + if gz_ext == '.gz': if not isfile(nogz): - raise RuntimeError('No nii for ' + fname) - flist[i] = nogz + contents = gzip.open(fname, 'rb').read() + with open(nogz, 'wb') as fobj: + fobj.write(contents) + flist[i] = nogz data_def['anatomical'] = anatomicals[0] data_def['functionals'] = functionals return data_def From 48aef2880d124a2e2c981a666bf0a3650ee66057 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 16:13:43 -0700 Subject: [PATCH 04/12] Allow bare .nii files instead of .nii.gz --- examples/interfaces/process_ds107.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index 336fcdbb6e..1e41d817bc 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -33,9 +33,9 @@ def get_data(data_path, subj_id): data_def = {} subject_path = pjoin(data_path, 'sub%03d' % subj_id) functionals = sorted( - glob(pjoin(subject_path, 'BOLD', 'task*', 'bold*.nii.gz'))) + glob(pjoin(subject_path, 'BOLD', 'task*', 'bold*.nii*'))) anatomicals = sorted( - glob(pjoin(subject_path, 'anatomy', 'highres001.nii.gz'))) + glob(pjoin(subject_path, 'anatomy', 'highres001.nii*'))) assert len(functionals) == 2 assert len(anatomicals) == 1 for flist in (anatomicals, functionals): From ace4765d7ce2cfaef8dfe9b3bdcdf0e59e3b30e8 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 17:10:43 -0700 Subject: [PATCH 05/12] Add prefix input output --- examples/interfaces/process_ds107.py | 41 ++++++++++++++++------------ 1 file changed, 24 insertions(+), 17 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index 1e41d817bc..fd7c5af921 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -78,7 +78,7 @@ def add_ana_defaults(self, ana_def): full_ana_def['fwhm'] = 8.0 return full_ana_def - def slicetime(self): + def slicetime(self, prefix=''): sess_scans = scans_for_fnames(self.data_def['functionals']) sdef = self.study_def stinfo = make_job('temporal', 'st', { @@ -90,9 +90,10 @@ def slicetime(self): 'refslice':1 }) run_jobdef(stinfo) + return 'a' + prefix - def realign(self): + def realign(self, prefix=''): sess_scans = scans_for_fnames( fnames_presuffix(self.data_def['functionals'], 'a')) rinfo = make_job('spatial', 'realign', [{ @@ -110,10 +111,11 @@ def realign(self): } }]) run_jobdef(rinfo) + return prefix - def reslice(self): + def reslice(self, prefix=''): sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], 'a')) + fnames_presuffix(self.data_def['functionals'], prefix)) rsinfo = make_job('spatial', 'realign', [{ 'write':{ 'data': np.vstack(sess_scans.flat), @@ -126,10 +128,11 @@ def reslice(self): } }]) run_jobdef(rsinfo) + return 'r' + prefix - def coregister(self): + def coregister(self, prefix=''): func1 = self.data_def['functionals'][0] - mean_fname = fname_presuffix(func1, 'meana') + mean_fname = fname_presuffix(func1, 'mean' + prefix) crinfo = make_job('spatial', 'coreg', [{ 'estimate':{ 'ref': [mean_fname], @@ -148,8 +151,9 @@ def coregister(self): } }]) run_jobdef(crinfo) + return prefix - def segnorm(self): + def segnorm(self, prefix=''): def_tpms = np.zeros((3,1), dtype=np.object) spm_path = spm_info.spm_path def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), @@ -179,10 +183,11 @@ def segnorm(self): } }) run_jobdef(sninfo) + return prefix - def norm_write(self): + def norm_write(self, prefix=''): sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], 'a')) + fnames_presuffix(self.data_def['functionals'], prefix)) matname = fname_presuffix(self.data_def['anatomical'], suffix='_seg_sn.mat', use_ext=False) @@ -210,8 +215,9 @@ def norm_write(self): subj['resample'][0] = self.data_def['anatomical'] roptions['interp'] = 4.0 run_jobdef(nwinfo) + return 'w' + prefix - def smooth(self): + def smooth(self, prefix=''): fwhm = self.ana_def['fwhm'] try: len(fwhm) @@ -219,12 +225,13 @@ def smooth(self): fwhm = [fwhm] * 3 fwhm = np.asarray(fwhm, dtype=np.float).reshape(1,3) sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], 'wa')) + fnames_presuffix(self.data_def['functionals'], prefix)) sinfo = make_job('spatial', 'smooth', {'data':np.vstack(sess_scans.flat), 'fwhm':fwhm, 'dtype':0}) run_jobdef(sinfo) + return 's' + prefix def process_subject(ddef, study_def, ana_def): @@ -234,13 +241,13 @@ def process_subject(ddef, study_def, ana_def): warn("No anatomical, aborting processing") return ana = SPMSubjectAnalysis(ddef, study_def, ana_def) - ana.slicetime() - ana.realign() - ana.reslice() - ana.coregister() + st_prefix = ana.slicetime('') + ana.realign(st_prefix) + # ana.reslice(prefix) + ana.coregister(st_prefix) ana.segnorm() - ana.norm_def() - ana.smooth() + n_st_prefix = ana.norm_write(st_prefix) + ana.smooth(n_st_prefix) def process_subjects(data_path, subj_ids, study_def, ana_def): From 8feda39ad58506815e5d3a0bcc10c040e4c5aed7 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 18:21:27 -0700 Subject: [PATCH 06/12] RF: throw out duplicate filenames, nii; nii.gz --- examples/interfaces/process_ds107.py | 34 ++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index fd7c5af921..e22efde147 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -28,16 +28,28 @@ ) +def _sorted_prefer_nii(file_list): + """ Strip any filanames ending nii.gz if matching .nii filename in list + """ + preferred = [] + for fname in file_list: + if not fname.endswith('.gz'): + preferred.append(fname) + else: + nogz, ext = splitext(fname) + if not nogz in file_list: + preferred.append(fname) + return sorted(preferred) + + def get_data(data_path, subj_id): data_path = abspath(data_path) data_def = {} subject_path = pjoin(data_path, 'sub%03d' % subj_id) - functionals = sorted( + functionals = _sorted_prefer_nii( glob(pjoin(subject_path, 'BOLD', 'task*', 'bold*.nii*'))) - anatomicals = sorted( + anatomicals = _sorted_prefer_nii( glob(pjoin(subject_path, 'anatomy', 'highres001.nii*'))) - assert len(functionals) == 2 - assert len(anatomicals) == 1 for flist in (anatomicals, functionals): for i, fname in enumerate(flist): nogz, gz_ext = splitext(fname) @@ -243,17 +255,17 @@ def process_subject(ddef, study_def, ana_def): ana = SPMSubjectAnalysis(ddef, study_def, ana_def) st_prefix = ana.slicetime('') ana.realign(st_prefix) - # ana.reslice(prefix) ana.coregister(st_prefix) ana.segnorm() n_st_prefix = ana.norm_write(st_prefix) ana.smooth(n_st_prefix) -def process_subjects(data_path, subj_ids, study_def, ana_def): +def get_subjects(data_path, subj_ids, study_def, ana_def): + ddefs = [] for subj_id in subj_ids: - ddef = get_data(data_path, subj_id) - process_subject(ddef, study_def, ana_def) + ddefs.append(get_data(data_path, subj_id)) + return ddefs if __name__ == '__main__': @@ -265,4 +277,8 @@ def process_subjects(data_path, subj_ids, study_def, ana_def): subj_ids = [int(id) for id in sys.argv[2:]] else: subj_ids = range(1, 16) - process_subjects(data_path, subj_ids, STUDY_DEF, {}) + for subj_id in subj_ids: + ddef = get_data(data_path, subj_id) + assert len(ddef['functionals']) == 2 + assert len(ddef['anatomicals']) == 1 + process_subject(ddef, STUDY_DEF, {}) From c99e3968c4b2d830be1eca3951a6ff32f2de1fd9 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Wed, 29 May 2013 18:28:12 -0700 Subject: [PATCH 07/12] RF: allow zero anatomicals --- examples/interfaces/process_ds107.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds107.py index e22efde147..591e5aea5a 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds107.py @@ -59,7 +59,10 @@ def get_data(data_path, subj_id): with open(nogz, 'wb') as fobj: fobj.write(contents) flist[i] = nogz - data_def['anatomical'] = anatomicals[0] + if len(anatomicals) == 0: + data_def['anatomical'] = None + else: + data_def['anatomical'] = anatomicals[0] data_def['functionals'] = functionals return data_def @@ -280,5 +283,4 @@ def get_subjects(data_path, subj_ids, study_def, ana_def): for subj_id in subj_ids: ddef = get_data(data_path, subj_id) assert len(ddef['functionals']) == 2 - assert len(ddef['anatomicals']) == 1 process_subject(ddef, STUDY_DEF, {}) From 54053a75f8edd5b76f0130d0237616c8c344858a Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Thu, 24 Oct 2013 00:01:32 -0700 Subject: [PATCH 08/12] RF: adapt script to ds105 Fix some bugs in script too. --- .../{process_ds107.py => process_ds105.py} | 53 ++++++++++++++----- 1 file changed, 40 insertions(+), 13 deletions(-) rename examples/interfaces/{process_ds107.py => process_ds105.py} (88%) diff --git a/examples/interfaces/process_ds107.py b/examples/interfaces/process_ds105.py similarity index 88% rename from examples/interfaces/process_ds107.py rename to examples/interfaces/process_ds105.py index 591e5aea5a..d20a6430b0 100755 --- a/examples/interfaces/process_ds107.py +++ b/examples/interfaces/process_ds105.py @@ -1,7 +1,18 @@ #!/usr/bin/env python # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -''' Single subject analysis script for SPM / Open FMRI ds107 ''' +''' Single subject analysis script for SPM / Open FMRI ds105 + +https://openfmri.org/dataset/ds000105 + +Download and extract the ds105 archive to some directory. + +Run this script with:: + + process_ds105.py ~/data/ds105 + +where ``~/data/ds105`` is the directory containing the ds105 data. +''' import sys from copy import deepcopy from os.path import join as pjoin, abspath, splitext, isfile @@ -20,9 +31,12 @@ # The batch scripts currently need SPM5 nimat.matlab_cmd = 'matlab-2007a-spm5 -nodesktop -nosplash' -N_SLICES = 37 +# This definition is partly for slice timing. We can't do slice timing for this +# dataset because the slice dimension is the first, and SPM assumes it is the +# last. +N_SLICES = 40 # X slices STUDY_DEF = dict( - TR = 3.0, + TR = 2.5, n_slices = N_SLICES, time_to_space = range(1, N_SLICES, 2) + range(2, N_SLICES, 2) ) @@ -110,7 +124,7 @@ def slicetime(self, prefix=''): def realign(self, prefix=''): sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], 'a')) + fnames_presuffix(self.data_def['functionals'], prefix)) rinfo = make_job('spatial', 'realign', [{ 'estimate':{ 'data':sess_scans, @@ -128,14 +142,21 @@ def realign(self, prefix=''): run_jobdef(rinfo) return prefix - def reslice(self, prefix=''): + def reslice(self, prefix='', out=('1..n', 'mean')): + which = [0, 0] + if 'mean' in out: + which[1] = 1 + if '1..n' in out or 'all' in out: + which[0] = 2 + elif '2..n' in out: + which[0] = 1 sess_scans = scans_for_fnames( fnames_presuffix(self.data_def['functionals'], prefix)) rsinfo = make_job('spatial', 'realign', [{ 'write':{ 'data': np.vstack(sess_scans.flat), 'roptions':{ - 'which':[2, 1], + 'which': which, 'interp':4.0, 'wrap':[0.0,0.0,0.0], 'mask':True, @@ -168,7 +189,7 @@ def coregister(self, prefix=''): run_jobdef(crinfo) return prefix - def segnorm(self, prefix=''): + def seg_norm(self, prefix=''): def_tpms = np.zeros((3,1), dtype=np.object) spm_path = spm_info.spm_path def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), @@ -256,10 +277,12 @@ def process_subject(ddef, study_def, ana_def): warn("No anatomical, aborting processing") return ana = SPMSubjectAnalysis(ddef, study_def, ana_def) - st_prefix = ana.slicetime('') + # st_prefix = ana.slicetime('') # We can't run slice timing + st_prefix = '' ana.realign(st_prefix) + ana.reslice(st_prefix, out=('mean',)) ana.coregister(st_prefix) - ana.segnorm() + ana.seg_norm() n_st_prefix = ana.norm_write(st_prefix) ana.smooth(n_st_prefix) @@ -271,16 +294,20 @@ def get_subjects(data_path, subj_ids, study_def, ana_def): return ddefs -if __name__ == '__main__': +def main(): try: data_path = sys.argv[1] except IndexError: - raise OSError('Need ds107 data path as input') + raise OSError('Need ds105 data path as input') if len(sys.argv) > 2: subj_ids = [int(id) for id in sys.argv[2:]] else: - subj_ids = range(1, 16) + subj_ids = range(1, 6) for subj_id in subj_ids: ddef = get_data(data_path, subj_id) - assert len(ddef['functionals']) == 2 + assert len(ddef['functionals']) in (11, 12) process_subject(ddef, STUDY_DEF, {}) + + +if __name__ == '__main__': + main() From 0608a41b64102cae9d459b2e5332d3c1109096b1 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Fri, 25 Oct 2013 08:08:54 -0700 Subject: [PATCH 09/12] RF: update spm interface for SPM8 batch jobs SPM8 batch jobs need extra 'initcfg' step to spm_jobman Add SPM version check to decided whether to add this initcfg step. --- nipy/interfaces/spm.py | 51 ++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 14 deletions(-) diff --git a/nipy/interfaces/spm.py b/nipy/interfaces/spm.py index a6ffe0cf26..4e424083e0 100644 --- a/nipy/interfaces/spm.py +++ b/nipy/interfaces/spm.py @@ -10,28 +10,47 @@ from scipy.io import savemat from nibabel import load -from nibabel.onetime import setattr_on_read from nibabel.tmpdirs import InTemporaryDirectory from .matlab import run_matlab_script class SpmInfo(object): - @setattr_on_read - def spm_path(self): - with InTemporaryDirectory() as tmpdir: - run_matlab_script(""" + def __init__(self): + self._spm_path = None + self._spm_ver = None + + def _set_properties(self): + with InTemporaryDirectory(): + run_matlab_script(r""" spm_path = spm('dir'); -fid = fopen('spm_path.txt', 'wt'); -fprintf(fid, '%s', spm_path); +spm_ver = spm('ver'); +fid = fopen('spm_stuff.txt', 'wt'); +fprintf(fid, '%s\n', spm_path); +fprintf(fid, '%s\n', spm_ver); fclose(fid); """) - spm_path = open('spm_path.txt', 'rt').read() - return spm_path + with open('spm_stuff.txt', 'rt') as fobj: + lines = fobj.readlines() + self._spm_path = lines[0].strip() + self._spm_ver = lines[1].strip() + + @property + def spm_path(self): + if self._spm_path is None: + self._set_properties() + return self._spm_path + + @property + def spm_ver(self): + if self._spm_ver is None: + self._set_properties() + return self._spm_ver + spm_info = SpmInfo() - + def make_job(jobtype, jobname, contents): return {'jobs':[{jobtype:[{jobname:contents}]}]} @@ -43,12 +62,16 @@ def fltcols(vals): def run_jobdef(jobdef): - with InTemporaryDirectory(): - savemat('pyjobs.mat', jobdef) - run_matlab_script(""" + script = """ load pyjobs; spm_jobman('run', jobs); -""") +""" + # Need initcfg for SPM8 + if spm_info.spm_ver != 'SPM5': + script = "spm_jobman('initcfg');\n" + script + with InTemporaryDirectory(): + savemat('pyjobs.mat', jobdef) + run_matlab_script(script) def scans_for_fname(fname): From 5b6fc82e13fe85a59b15461725e60956ef824fa5 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Fri, 25 Oct 2013 08:10:02 -0700 Subject: [PATCH 10/12] RF: update ds105 example for SPM8 batching --- examples/interfaces/process_ds105.py | 63 +++++++++++++++------------- 1 file changed, 34 insertions(+), 29 deletions(-) diff --git a/examples/interfaces/process_ds105.py b/examples/interfaces/process_ds105.py index d20a6430b0..a2147a848a 100755 --- a/examples/interfaces/process_ds105.py +++ b/examples/interfaces/process_ds105.py @@ -29,7 +29,7 @@ # The batch scripts currently need SPM5 -nimat.matlab_cmd = 'matlab-2007a-spm5 -nodesktop -nosplash' +nimat.matlab_cmd = 'matlab-spm8 -nodesktop -nosplash' # This definition is partly for slice timing. We can't do slice timing for this # dataset because the slice dimension is the first, and SPM assumes it is the @@ -107,8 +107,9 @@ def add_ana_defaults(self, ana_def): full_ana_def['fwhm'] = 8.0 return full_ana_def - def slicetime(self, prefix=''): - sess_scans = scans_for_fnames(self.data_def['functionals']) + def slicetime(self, in_prefix='', out_prefix='a'): + sess_scans = scans_for_fnames( + fnames_presuffix(self.data_def['functionals'], in_prefix)) sdef = self.study_def stinfo = make_job('temporal', 'st', { 'scans': sess_scans, @@ -116,15 +117,16 @@ def slicetime(self, prefix=''): 'tr': sdef['TR'], 'ta': sdef['TA'], 'nslices': float(sdef['n_slices']), - 'refslice':1 + 'refslice':1, + 'prefix': out_prefix, }) run_jobdef(stinfo) - return 'a' + prefix + return out_prefix + in_prefix - def realign(self, prefix=''): + def realign(self, in_prefix=''): sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], prefix)) + fnames_presuffix(self.data_def['functionals'], in_prefix)) rinfo = make_job('spatial', 'realign', [{ 'estimate':{ 'data':sess_scans, @@ -140,9 +142,9 @@ def realign(self, prefix=''): } }]) run_jobdef(rinfo) - return prefix + return in_prefix - def reslice(self, prefix='', out=('1..n', 'mean')): + def reslice(self, in_prefix='', out_prefix='r', out=('1..n', 'mean')): which = [0, 0] if 'mean' in out: which[1] = 1 @@ -151,7 +153,7 @@ def reslice(self, prefix='', out=('1..n', 'mean')): elif '2..n' in out: which[0] = 1 sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], prefix)) + fnames_presuffix(self.data_def['functionals'], in_prefix)) rsinfo = make_job('spatial', 'realign', [{ 'write':{ 'data': np.vstack(sess_scans.flat), @@ -160,20 +162,22 @@ def reslice(self, prefix='', out=('1..n', 'mean')): 'interp':4.0, 'wrap':[0.0,0.0,0.0], 'mask':True, + 'prefix': out_prefix } } }]) run_jobdef(rsinfo) - return 'r' + prefix + return out_prefix + in_prefix - def coregister(self, prefix=''): + def coregister(self, in_prefix=''): func1 = self.data_def['functionals'][0] - mean_fname = fname_presuffix(func1, 'mean' + prefix) + mean_fname = fname_presuffix(func1, 'mean' + in_prefix) crinfo = make_job('spatial', 'coreg', [{ 'estimate':{ - 'ref': [mean_fname], - 'source': [self.data_def['anatomical']], - 'other': [[]], + 'ref': np.asarray(mean_fname, dtype=object), + 'source': np.asarray(self.data_def['anatomical'], + dtype=object), + 'other': [''], 'eoptions':{ 'cost_fun':'nmi', 'sep':[4.0, 2.0], @@ -187,9 +191,9 @@ def coregister(self, prefix=''): } }]) run_jobdef(crinfo) - return prefix + return in_prefix - def seg_norm(self, prefix=''): + def seg_norm(self, in_prefix=''): def_tpms = np.zeros((3,1), dtype=np.object) spm_path = spm_info.spm_path def_tpms[0] = pjoin(spm_path, 'tpm', 'grey.nii'), @@ -219,11 +223,11 @@ def seg_norm(self, prefix=''): } }) run_jobdef(sninfo) - return prefix + return in_prefix - def norm_write(self, prefix=''): + def norm_write(self, in_prefix='', out_prefix='w'): sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], prefix)) + fnames_presuffix(self.data_def['functionals'], in_prefix)) matname = fname_presuffix(self.data_def['anatomical'], suffix='_seg_sn.mat', use_ext=False) @@ -238,6 +242,7 @@ def norm_write(self, prefix=''): 'vox':fltcols([2.0,2.0,2.0]), 'interp':1.0, 'wrap':[0.0,0.0,0.0], + 'prefix': out_prefix, } nwinfo = make_job('spatial', 'normalise', [{ 'write':{ @@ -251,9 +256,9 @@ def norm_write(self, prefix=''): subj['resample'][0] = self.data_def['anatomical'] roptions['interp'] = 4.0 run_jobdef(nwinfo) - return 'w' + prefix + return out_prefix + in_prefix - def smooth(self, prefix=''): + def smooth(self, in_prefix='', out_prefix='s'): fwhm = self.ana_def['fwhm'] try: len(fwhm) @@ -261,13 +266,13 @@ def smooth(self, prefix=''): fwhm = [fwhm] * 3 fwhm = np.asarray(fwhm, dtype=np.float).reshape(1,3) sess_scans = scans_for_fnames( - fnames_presuffix(self.data_def['functionals'], prefix)) + fnames_presuffix(self.data_def['functionals'], in_prefix)) sinfo = make_job('spatial', 'smooth', {'data':np.vstack(sess_scans.flat), 'fwhm':fwhm, 'dtype':0}) run_jobdef(sinfo) - return 's' + prefix + return out_prefix + in_prefix def process_subject(ddef, study_def, ana_def): @@ -279,9 +284,9 @@ def process_subject(ddef, study_def, ana_def): ana = SPMSubjectAnalysis(ddef, study_def, ana_def) # st_prefix = ana.slicetime('') # We can't run slice timing st_prefix = '' - ana.realign(st_prefix) - ana.reslice(st_prefix, out=('mean',)) - ana.coregister(st_prefix) + ana.realign(in_prefix=st_prefix) + ana.reslice(in_prefix=st_prefix, out=('mean',)) + ana.coregister(in_prefix=st_prefix) ana.seg_norm() n_st_prefix = ana.norm_write(st_prefix) ana.smooth(n_st_prefix) @@ -302,7 +307,7 @@ def main(): if len(sys.argv) > 2: subj_ids = [int(id) for id in sys.argv[2:]] else: - subj_ids = range(1, 6) + subj_ids = range(1, 7) for subj_id in subj_ids: ddef = get_data(data_path, subj_id) assert len(ddef['functionals']) in (11, 12) From 26665b21c7096277b60d3209413c2d2e11d584d2 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Tue, 29 Oct 2013 12:33:10 -0700 Subject: [PATCH 11/12] DOC: point to nipype for extended use --- examples/interfaces/process_ds105.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/examples/interfaces/process_ds105.py b/examples/interfaces/process_ds105.py index a2147a848a..028928f9bd 100755 --- a/examples/interfaces/process_ds105.py +++ b/examples/interfaces/process_ds105.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -''' Single subject analysis script for SPM / Open FMRI ds105 +''' Single subject analysis script for SPM / Open FMRI ds105 https://openfmri.org/dataset/ds000105 @@ -12,6 +12,13 @@ process_ds105.py ~/data/ds105 where ``~/data/ds105`` is the directory containing the ds105 data. + +The example uses the very basic MATLAB / SPM interface routines in NIPY. + +If you need more than very basic use, please consider using nipype. nipype has +extended capabilities to interface with external tools and for dataflow +management. nipype can handle vanilla SPM in MATLAB or SPM run through the +MATLAB common runtime (free from MATLAB Licensing). ''' import sys from copy import deepcopy From da010df9a5319d9451d97932d6ac1385442c5601 Mon Sep 17 00:00:00 2001 From: Matthew Brett Date: Thu, 31 Oct 2013 12:34:32 -0700 Subject: [PATCH 12/12] BF: update call to savemat for compatibility Latest scipy assumes oned_as='row' default. Enforce in savemat call --- nipy/interfaces/spm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipy/interfaces/spm.py b/nipy/interfaces/spm.py index 4e424083e0..c27ae05a16 100644 --- a/nipy/interfaces/spm.py +++ b/nipy/interfaces/spm.py @@ -70,7 +70,7 @@ def run_jobdef(jobdef): if spm_info.spm_ver != 'SPM5': script = "spm_jobman('initcfg');\n" + script with InTemporaryDirectory(): - savemat('pyjobs.mat', jobdef) + savemat('pyjobs.mat', jobdef, oned_as='row') run_matlab_script(script)