mripy package

Subpackages

• mripy.io package
  • Submodules
  • mripy.io.freesurfer module
    • read_fs_annotation()
    • read_fs_curv()
    • read_fs_int32()
    • read_fs_patch()
    • read_fs_str()
    • read_fs_surf()
    • read_fs_uint24()
    • write_fs_color_table()
    • write_fs_curv()
    • write_fs_int32()
    • write_fs_uint24()
  • mripy.io.gifti module
    • read_gii_dset()
  • mripy.io.niml module
    • parse_attr()
    • parse_data()
    • parse_data_format()
    • parse_ni_type()
    • parse_niml()
    • read_until()
  • Module contents
    • BallMask
    • CylinderMask
    • Mask
      • Mask.ball()
      • Mask.compatible()
      • Mask.concat()
      • Mask.constrain()
      • Mask.cylinder()
      • Mask.dump()
      • Mask.from_dict()
      • Mask.from_expr()
      • Mask.from_files()
      • Mask.ijk
      • Mask.infer_selector()
      • Mask.near()
      • Mask.pick()
      • Mask.slab()
      • Mask.to_dict()
      • Mask.to_file()
      • Mask.undump()
      • Mask.xyz
      • Mask.xyz_nifti
    • MaskDumper
      • MaskDumper.dump()
      • MaskDumper.undump()
    • SlabMask
    • change_dim_order()
    • change_space()
    • compress()
    • convert_dicom()
    • convert_dicoms()
    • decompress()
    • extract_physio()
    • filter_cluster()
    • filter_dicom_files()
    • generate_afni_idcode()
    • get_dim_order()
    • get_ni_type()
    • get_space()
    • hms2dt()
    • match_physio_with_series()
    • mmn2dt()
    • parse_dicom_header()
    • parse_physio_file()
    • parse_physio_files()
    • parse_series_info()
    • parse_slice_order()
    • read_affine()
    • read_afni()
    • read_asc()
    • read_gii()
    • read_label()
    • read_nii()
    • read_niml_bin_nodes()
    • read_niml_dset()
    • read_patch_asc()
    • read_register_dat()
    • read_stim()
    • read_surf_data()
    • read_surf_info()
    • read_surf_mesh()
    • read_txt()
    • read_vol()
    • read_warp()
    • sort_dicom_series()
    • unzipped
    • write_1D_nodes()
    • write_affine()
    • write_afni()
    • write_asc()
    • write_gii()
    • write_nii()
    • write_niml_bin_nodes()
    • write_surf_data()
    • write_surf_mesh()
    • write_vol()
• mripy.scripts package
  • Submodules
  • mripy.scripts.afni_viewer module
  • mripy.scripts.extract_physio module
  • mripy.scripts.mripy_1dplot module
  • mripy.scripts.mripy_run_iglesias18 module
  • mripy.scripts.mripy_run_wang15 module
  • mripy.scripts.report_parameters module
  • mripy.scripts.script_utils module
  • mripy.scripts.sort_dicom module
  • Module contents
• mripy.tests package
  • Submodules
  • mripy.tests.context module
  • mripy.tests.test_afni module
    • test_afni
      • test_afni.test_get_affine()
      • test_afni.test_get_prefix()
      • test_afni.test_get_suma_spec()
      • test_afni.test_substitute_hemi()
  • mripy.tests.test_io module
    • test_io
      • test_io.test_Mask()
  • mripy.tests.test_timecourse module
    • test_Attributes
      • test_Attributes.setUp()
      • test_Attributes.test_copy()
      • test_Attributes.test_pick()
    • test_Epochs_Attributes
      • test_Epochs_Attributes.setUp()
      • test_Epochs_Attributes.test_copy()
      • test_Epochs_Attributes.test_pick()
  • mripy.tests.test_utils module
    • test_utils
      • test_utils.test_fname_with_ext()
  • mripy.tests.test_utils_slow module
    • test_utils
      • test_utils.test_SharedMemoryArray_CoW()
      • test_utils.test_SharedMemoryArray_array()
      • test_utils.test_SharedMemoryArray_memory()
  • Module contents

Submodules

mripy.afni module

mripy.afni.add_colormap(cmap, name=None, cyclic=False, index=None, categorical=False)

cmap : list of RGB colors | matplotlib.colors.LinearSegmentedColormap

mripy.afni.call(cmd)

mripy.afni.check_output(cmd, tags=None, pattern=None, verbose=0, **kwargs)

The syntax of subprocess.check_output(shell=False) is tedious for long cmd. But for security reason, we don’t want to use shell=True for external cmd. This helper function allows you to execute a single cmd without shell=True.

Parameters:

• cmd (str) – A single command string packed with all options (but no wildcard)

• **kwargs – Go to subprocess.check_output(**kwargs)

Returns:

lines – Much easier to deal with compared with subprocess.check_output()

Return type:

list of lines

mripy.afni.filter_output(lines, tags=None, pattern=None, ex_tags=None, ex_pattern=None)

Filter output lines according to their initial tags (++, +, *, etc.) and/or a regex search pattern.

Parameters:

• tags (list of tags) – Default is [], which means all lines will pass the filter.

• pattern (str) –

• ex_tags (list of tags to exclude) –

• ex_pattern (str) –

mripy.afni.generate_spec(fname, surfs, ext=None, **kwargs)

mripy.afni.get_DELTA(fname)

mripy.afni.get_DIMENSION(fname)

[x, y, z, t, 0] Not work for bucket.

mripy.afni.get_ORIENT(fname, format='str')

Parameters:

format (str, {'code', 'str', 'mat', 'sorter'}) –

References

[1] https://afni.nimh.nih.gov/pub/dist/doc/program_help/README.attributes.html

#define ORI_R2L_TYPE 0 // Right to Left #define ORI_L2R_TYPE 1 // Left to Right #define ORI_P2A_TYPE 2 // Posterior to Anterior #define ORI_A2P_TYPE 3 // Anterior to Posterior #define ORI_I2S_TYPE 4 // Inferior to Superior #define ORI_S2I_TYPE 5 // Superior to Inferior

Thus “0 3 4” is standard DICOM Reference Coordinates System, i.e., RAI. The AFNI convention is also that R-L, A-P, and I-S are negative-to-positive, i.e., RAI.

[2] https://nipy.org/nibabel/nifti_images.html

On the other hand, NIFTI images have an affine relating the voxel coordinates to world coordinates in RAS+ space, or LPI in AFNI’s term.

mripy.afni.get_ORIGIN(fname)

mripy.afni.get_S2E_mat(fname, mat='S2E')

mripy.afni.get_TR(fname)

mripy.afni.get_affine(fname)

mripy.afni.get_affine_nifti(fname)

mripy.afni.get_attribute(fname, name, type=None)

mripy.afni.get_brick_labels(fname, label2index=False)

mripy.afni.get_censor_from_X(fname)

Return whether it is selected (“good” = not censored) for each volume.

Parameters:

fname (str) – X.prefix.1D (design matrix) file as generated by 3dDeconvolve.

mripy.afni.get_crop(fname, xyz, r)

For AFNI driver “crop=x1:x2,y1:y2” See https://afni.nimh.nih.gov/pub/dist/doc/program_help/README.driver.html

Parameters:

• fname (str) –

• xyz (array-like) – Cursor location [x, y, z] in mm (as shown in AFNI GUI).

• r (float) – Crop radius in mm.

mripy.afni.get_dims(fname)

Dimensions (number of voxels) of the data matrix. See also: get_head_dims

mripy.afni.get_head_delta(fname)

Resolution (voxel size) along R-L, A-P, I-S axes.

mripy.afni.get_head_dims(fname)

Dimensions (number of voxels) along R-L, A-P, I-S axes. See also: get_dims

mripy.afni.get_head_extents(fname)

Spatial extent along R, L, A, P, I and S.

mripy.afni.get_hemi(fname)

mripy.afni.get_nifti_field(fname, name, type=None)

mripy.afni.get_prefix(fname, with_path=False)

Return “dset” given “path/to/dset+orig.HEAD”, “dset+orig.”, “dset+tlrc”, “dsets”

mripy.afni.get_runs_from_X(fname, fmt='groups')

Return groups (run index started from 0) for each volume.

Parameters:

• fname (str) – X.prefix.1D (design matrix) file as generated by 3dDeconvolve.

• fmt (str, ‘starts’|’ends’|<'groups'>|’afni’) – The return format about runs information. - ‘starts’: The start volume index for each run, e.g., [0, 3] - ‘ends’: [start, end] volume index for each run, e.g., [[0,2],[3,5]] - ‘groups’: For use with sklearn, e.g., [0,0,0,1,1,1] - ‘afni’: For use with afni subbrick selection, e.g., [”’[0..2]’”,”’[3..5]’”]

mripy.afni.get_slice_timing(fname)

mripy.afni.get_subbrick_selector(good)

Convert boolean indexing to AFNI’s subbrick selection syntax. E.g., [True, False, True, True, True] -> “0,2..4”

mripy.afni.get_suma_info(suma_dir, suma_spec=None)

mripy.afni.get_suma_spec(suma_spec)

Infer other spec files from one spec file (either lh.spec, rh.spec, or both.spec).

Parameters:

suma_spec (str) – Either a .spec file or the suma_dir.

mripy.afni.get_suma_subj(suma_dir)

Infer SUMA subject given path to SUMA folder.

mripy.afni.get_surf_type(suma_dir)

Infer SUMA surface mesh file type (.gii vs .asc).

mripy.afni.get_surf_vol(suma_dir)

Infer SUMA SurfVol filename with full path (agnostic about file type: .nii vs +orig.HEAD/BRIK).

mripy.afni.infer_surf_dset_variants(fname, hemis=['lh', 'rh', 'both', 'mh', 'bh'])

>>> infer_surf_dset_variants('data.niml.dset')
{'lh': 'lh.data.niml.dset', 'rh': 'rh.data.niml.dset', 'both': 'both.data.niml.dset', mh': 'mh.data.niml.dset'}
>>> infer_surf_dset_variants('lh.data.niml.dset')
{'lh': 'lh.data.niml.dset'}

Parameters:

fname (str, list, or dict) –

mripy.afni.insert_suffix(fname, suffix)

mripy.afni.parse_patch(patch)

Notes

1. Each replacement is started with one or more comment lines. The last comment line is treated as replacement target, which may contain an optional replacement directive at the end: # This is an example <replace command=”1”/>

   Possible directives for replacing the original scripts includes:

   1. command=”n”: replace n commands

   2. line=”n”: replace n lines

   3. until=”regexp”: replace until a specific line (the regexp is the last line to be replaced)

2. Each replacement must end with two consecutive newlines.

mripy.afni.patch_afni_proc(original, patch, inplace=True)

mripy.afni.set_attribute(fname, name, value, type=None)

mripy.afni.set_brick_labels(fname, labels)

mripy.afni.set_nifti_field(fname, name, value, out_file=None)

mripy.afni.set_slice_timing(fname, times, TR)

We have to provide a TR because we don’t know whether the default value TR=1.0 is valid.

mripy.afni.split_out_file(out_file, split_path=False, trailing_slash=False)

Ensure that path.join(out_dir, prefix, ext) can be checked by path.exists().

>>> split_out_file('dset.nii')
('dset', '.nii')
>>> split_out_file('dset.1D')
('dset', '.1D')
>>> split_out_file('folder/dset')
('folder/dset', '+orig.HEAD')
>>> split_out_file('folder/dset+orig', split_path=True)
('folder', 'dset', '+orig.HEAD')
>>> split_out_file('dset+orig.', split_path=True)
('', 'dset', '+orig.HEAD')
>>> split_out_file('folder/dset+orig.HEAD', split_path=True, trailing_slash=True)
('folder/', 'dset', '+orig.HEAD')
>>> split_out_file('dset+tlrc.BRIK', split_path=True, trailing_slash=True)
('', 'dset', '+tlrc.HEAD')

mripy.afni.substitute_hemi(fname, hemi='{0}')

mripy.afni.update_afnirc(**kwargs)

mripy.afni.write_colorscale_file(fname, pal_name, colors, locations=None, interp=None)

Parameters:

• fname (*.pal file name) –

• pal_name (palette name (or title)) –

• colors (if you fill the colorscale file with a lot of) – first color (bottom) -> last color (top)

• locations (locations of the breakpoints where colors are defined) – 0 (bottom) -> 1 (top)

• interp ('linear'|'nearest') –

• colorscale." (AFNI document says "There are exactly 128 color locations on an AFNI) –

• details (For) –

• https (see) –

• fact (But in) –

• colors –

• used. (only the first 256 colors will be) –

mripy.afni.xyz2ijk(fname, xyz)

Convert RAI=dicom xyz to data matrix indices ijk

mripy.dcm module

class mripy.dcm.ExtendedCircle(width, pad=0.3)

Bases: Circle

An extended Circle BoxStyle that keeps a minimum predefined width parameter.

References

https://stackoverflow.com/questions/40796117/how-do-i-make-the-width-of-the-title-box-span-the-entire-plot https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/patches.py

transmute(x0, y0, width, height, mutation_size)

x0 and y0 are the lower left corner of original bbox. They are set automatically by matplotlib.

class mripy.dcm.ExtendedSimple(head_length=0.5, head_width=0.5, tail_width=0.2)

Bases: _Base

A simple arrow. Only works with a quadratic Bezier curve.

transmute(path, mutation_size, linewidth)

The transmute method is the very core of the ArrowStyle class and must be overridden in the subclasses. It receives the path object along which the arrow will be drawn, and the mutation_size, with which the arrow head etc. will be scaled. The linewidth may be used to adjust the path so that it does not pass beyond the given points. It returns a tuple of a Path instance and a boolean. The boolean value indicate whether the path can be filled or not. The return value can also be a list of paths and list of booleans of a same length.

mripy.dcm.draw_circle_arrow(xy, radius, patchA=None, ax=None)

References

https://stackoverflow.com/questions/37512502/how-to-make-arrow-that-loops-in-matplotlib/38224201

mripy.dcm.plot_dcm(rois, inputs, A, B, C, PA=None, PB=None, PC=None, roi_order=None, rads=None, ax=None)

Visualize a DCM (as well as estimated A, B, C parameters).

Parameters:

• A ([n_rois, n_rois]) – Baseline connectivity, dest <- src

• B ([n_rois, n_rois, n_inputs]) – Modulation on connectivity (non-zero for modulatory inputs only)

• C ([n_rois, n_inputs]) – Driving on nodes (non-zero for driving inputs only)

References

https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.annotate.html https://matplotlib.org/stable/gallery/userdemo/annotate_simple_coord01.html

https://matplotlib.org/3.3.4/gallery/lines_bars_and_markers/gradient_bar.html (gradient fill)

mripy.dcm.plot_peb_bma(GCM, **kwargs)

mripy.decoding module

class mripy.decoding.Demeaner

Bases: TransformerMixin, BaseEstimator

Remove the mean of each sample individually.

For fMRI, this removes mean activation within ROI for each sample.

fit(X, y=None)

Do nothing and return the estimator unchanged. This method is just there to implement the usual API and hence work in pipelines. :param X: The data to estimate the demean parameters. :type X: {array-like, sparse matrix} of shape (n_samples, n_features) :param y: Not used, present here for API consistency by convention. :type y: Ignored

Returns:

self – Fitted transformer.

Return type:

object

transform(X)

Subtract the mean of each row in X. :param X: The data to demean, row by row. :type X: {array-like, sparse matrix} of shape (n_samples, n_features)

Returns:

X_tr – Transformed array.

Return type:

{ndarray, sparse matrix} of shape (n_samples, n_features)

mripy.decoding.compute_critical_value(x, y, permute='permute', data=None, alpha=0.05, tail=2)

Get critical values based on permutation distribution, and account for multiple comparisons using extreme statistics.

Parameters:

• x (str, list of str) – Columns along which multiple comparisons occur (e.g., roi, time).

• y (str) – Column for performance measurement (e.g., test_accuracy, PC, RT).

• data (pd.DataFrame(x, y, permute)) – permute == 0 is originally observed data, >= 1 is permutation data.

mripy.decoding.cross_validate_ext(model, X, y, groups=None, cv=None, pred_kws=None, method=None)

mripy.decoding.cross_validate_with_permutation(model, X, y, groups, rois=None, n_permutations=1000, scoring=None, cv=None)

mripy.decoding.permute_within_group(y, groups)

mripy.decoding.standardize_within_group(X, groups, with_mean=True, with_std=True)

This is an extension of the “mean centering” method proposed in [1]. Can be used as a replacement for the training-set-wise standardization. Both with_mean and with_std may provide some extra performance.

References

[1] Lee, S., & Kable, J. W. (2018). Simple but robust improvement in

multivoxel pattern classification. PloS One, 13(11), e0207083.

mripy.dicom module

mripy.dicom.parse_SQ_data_element(fi)

We only support Data Element with Explicit VR at present (Table 7.1-2). We don’t support nested Item at present (2018-09-26).

References

[1] http://dicom.nema.org/Dicom/2013/output/chtml/part05/chapter_7.html

mripy.dicom.parse_Siemens_CSA(b)

mripy.dicom.parse_Siemens_CSA2(b)

mripy.dicom.parse_dicom_header(fname, search_for_tags=None, **kwargs)

Parameters:

• fname (str) –

• search_for_tags (set) – Search for specific dicom tags, and stop file scanning early if all tags of interest are seen. e.g., search_for_tags={‘0020,0011’, ‘0020,0013’} will search for SeriesNumber and InstanceNumber. This will save you some time, esp. when the remote file is accessed via slow data link.

• **kwargs – This is only for backward compatibility.

Notes

“Implicit and Explicit VR Data Elements shall not coexist in a Data Set and Data Sets nested within it (see Section 7.5). Whether a Data Set uses Explicit or Implicit VR, among other characteristics, is determined by the negotiated Transfer Syntax (see Section 10 and Annex A).” [1]

References

[1] http://dicom.nema.org/Dicom/2013/output/chtml/part05/chapter_7.html [2] https://stackoverflow.com/questions/119684/parse-dicom-files-in-native-python

mripy.dicom.parse_series_info(dicom_files, dicom_ext=None, parser=None, return_headers=False)

Parameters:

dicom_files (list or str) – A list of dicom files (e.g., as provided by sort_dicom_series), or a folder that contains a single series (e.g., “../raw_fmri/func01”), or a single dicom file.

mripy.dicom.sort_dicom_series(folder)

Parameters:

folder (string) – Path to the folder containing all the dicom files.

Returns:

studies – [{‘0001’: [file0, file1, …], ‘0002’: [files], …}, {study1}, …]

Return type:

list of dicts

mripy.dicom_report module

mripy.dicom_report.inspect_mp2rage(data_dir, subdir_pattern='T1??')

mripy.dicom_report.pares_patent_age(age)

mripy.dicom_report.print_subject_info(dir_pattern, exclude=None)

mripy.dicom_report.remove_duplicate_parts(folders)

mripy.dicom_report.report_parameters(dicom_folder, preset=None, return_preset=False, return_info=False)

mripy.encoding module

class mripy.encoding.BaseModel

Bases: Savable2

from_dict(d)

get_params(deep=True)

set_params(**parameters)

to_dict()

class mripy.encoding.BayesianChannelModel(n_channels='required', basis_func='required', stimulus_domain='required', circular=False, stimulus_prior=None, global_search=False, verbose=2)

Bases: BaseModel

property Omega_

property Omega_inv_

bayesian_inversion(X, stimulus_domain=None, stimulus_prior=None, density=True)

Parameters:

• X (2D array, n_trials * n_voxels) –

• stimulus_domain (1D array, n_domain) –

• stimulus_prior (2D array, n_trials * n_domain (or 1D array, n_domain)) – None for a flat stimulus prior, same for all trials.

Returns:

posterior

Return type:

2D array, n_trials * n_domain

fit(X, y)

Parameters:

• X (2D array) – n_trials * n_voxels BOLD response pattern (e.g., beta for each trial, or delayed and detrended time points within block plateau)

• y (1D array) – n_trials stimulus value (e.g., orientation, color)

get_params(deep=True)

loglikelihood(X, y)

predict(X, stimulus_domain=None, stimulus_prior=None, return_all=False)

class mripy.encoding.ChannelEncodingModel(n_channels, basis_func, stimulus_domain, circular=False, verbose=2)

Bases: BaseModel

channel_inversion(X, stimulus_domain=None)

correlation_inversion(X, stimulus_domain=None)

fit(X, y)

Parameters:

• X (2D array) – n_trials * n_voxels BOLD response pattern (e.g., beta for each trial, or delayed and detrended time points within block plateau)

• y (1D array) – n_trials stimulus value (e.g., orientation, color)

get_params(deep=True)

inverted_encoding(X)

pRF(stimulus_domain=None, method='ols', X=None, y=None)

predict(X, stimulus_domain=None, return_all=False)

voxel_inversion(X, stimulus_domain=None)

class mripy.encoding.EnsembleModel(n_ensemble=10, base_model='required', pred_method=None, pred_options=None)

Bases: BaseModel

fit(X, y)

from_dict(d)

get_params(deep=True)

predict(X, method=None, options=None, return_all=False, pred_kws=None)

to_dict()

mripy.encoding.basis_Sprague2013(s, n_channels=6, spacing=2, center=None, size=None, power=7, dim=1, intercept=False)

Parameters:

s (2D array, n_trials * dim | 1D array, n_trials) –

Returns:

fs

Return type:

2D array, n_channels * n_trials

References

{Sprague2013}

mripy.encoding.basis_vanBergen2015(s, n_channels=8, power=5)

Parameters:

s (1D array, n_trials) –

Returns:

fs

Return type:

2D array, n_channels * n_trials

References

{van Bergen2015}

mripy.encoding.circular_correct(y_true, y_pred, domain=None, n_targets=None, tolerance=None, return_dist=False)

mripy.encoding.discretize_prediction(y_pred, targets, circular_domain=None)

Discretize continous prediction to the nearest target. Can handle irregular target grid and also circular domain.

E.g., y_pred = encoding.discretize_prediction(y_hat, arange(8)/8*pi, circular_domain=[0, pi]) correct = encoding.circular_correct(y_true, y_hat, domain=[0, pi], n_targets=8) assert(allclose(mean(y_pred==y_true), mean(correct)))

mripy.encoding.shift_distribution(d, stimulus_domain, center_on=None, circular=True)

Parameters:

d (n_trials * n_domain) –

mripy.evaluation module

mripy.evaluation.afni_costs(base_file, in_file, mask=None)

In general, the alignment algorithm want to minimize the cost.

Some metrics in their canonical form are maximized when perfectly aligned (e.g., mi and hel). These are usually negated when output by -allcost, so that ALL cost (EXCEPT lss) values will decrease when alignment improves (given that the base/source contrasts are valid as required by the metric, e.g., lpc is ONLY meaningful if the base and source have opposite contrast). This behavior can be verified by manually applying a transform to a image.

BASED ON GLOBAL CORRELATION COEFFICIENT ls : 1 - abs(Pearson corrcoef), near zero sp : 1 - abs(Spearman corrcoef), near zero lss (maximized): Pearson corrcoef, near one (the ONLY exception)

BASED ON 2D JOINT HISTOGRAM mi (negated) : (negative) mutual information, large negative value nmi : 1/normalized MI, small positive value je : joint entropy, small positive value hel (negated) : (negative) Hellinger distance, large negative value crM : near zero crA : near zero crU : near zero

BASED ON LOCAL (A FEW MM) PEARSON CORRELATION COEFFICIENT (SIMILAR TO BOUNDARY-BASED METHODS) lpc : sum(w[i]*pc[i]*abs(pc[i])) / sum(w[i]), large negative value (FOR opposite contrast ONLY) lpa : 1 - abs(lpc), large negative value (FOR same contrast ONLY) lpc+ : lpc + hel*0.4 + crA*0.4 + nmi*0.2 + mi*0.2 + ov*0.4 lpa+ : lpa + hel*1.0 + crA*0.4 + nmi*0.2 + mi*0.0 + ov*0.0

References

[1] https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dAllineate.html

See “Cost functional descriptions (for use with -allcost output)” section

[2] https://www.youtube.com/watch?v=PaZinetFKGY&list=PL_CD549H9kgqJ1GDXAs1BWkgEimAHZeNX

mripy.evaluation.residual_motion(files, print_table=False)

mripy.math module

class mripy.math.DomainMapper(domain=None)

Bases: object

from2pi(y)

to2pi(x)

mripy.math.LPI2RAI_affine(mat)

mripy.math.RAI2LPI_affine(mat)

mripy.math.apply_affine(mat, xyz)

xyz : 3xN array

mripy.math.argsort_rows(rows, cols=None)

mripy.math.array2dataframe(a, factors, var_name='value')

a : array factors : dict (factor->levels)

mripy.math.circular_corrcoef(x1, x2, domain=None, n_perm=1000, ci=0.95, n_boot=None)

The complex corrcoef method used here is fundamentally different from the (Fisher & Lee, 1983) method which is implemented by pycircstat.corrcc(). For more details, read my note https://docs.google.com/document/d/1sl39YH3g3TFQu1zX-Ax477NULEyJE–MHMXg2gg0cyk/edit

mripy.math.circular_mean(x, domain=None, weight=None, axis=None)

Circular mean for values from arbitary circular domain (not necessarily angles).

mripy.math.circular_std(x, domain=None, weight=None, axis=None)

Following scipy.stats.circstd()’s definition of circular standard deviation that in the limit of small angles returns the ‘linear’ standard deviation. scipy.stats.circstd() doesn’t support weight. pycircstat.std() doesn’t support domain. astropy.stats.circvar() follows another definition.

mripy.math.concat_affine(mat2, mat1)

mat @ v = mat2 @ mat1 @ v

For afni’s “fetching” convention, the first transform goes first. mat1 = io.read_affine(‘SurfVoltoT1.aff12.1D’) # Fetching SurfVol from T1 grid mat2 = io.read_affine(‘T1toEPI.aff12.1D’) # Fetching T1 from EPI grid mat = math.concat_affine(mat1, mat2) # Note the order! The combined transform needs to fetch SurfVol from EPI io.write_affine(‘SurfVoltoEPI.aff12.1D’, mat) # Fetching SurfVol from EPI

mripy.math.corrcoef_along_axis(x, y, axis=None, norm=True)

mripy.math.gaussian_logpdf(x, mean, cov, cov_inv=None, axis=-1)

More efficient multivariate normal distribution log pdf than stats.multivariate_normal.logpdf()

mripy.math.invert_affine(mat)

mripy.math.median_argmax(x, axis=-1)

mripy.math.nearest(x, parity='odd', round=<function round_>)

mripy.math.normalize_logP(logP, axis=None)

mripy.math.pinv(x)

mripy.math.polyfit3d(x, y, z, f, deg, method=None)

mripy.math.tsarray2dataframe(tsarray, t=None, ts_name='value', t_name='time', trial_name='trial', trial_df=None)

Parameters:

tsarray (ndarray, n_trials * n_times) –

Returns:

df

Return type:

DataFrame

mripy.paraproc module

class mripy.paraproc.ArrayWrapper(name, bases, dct)

Bases: type

This is the metaclass for classes that wrap an np.ndarray and delegate non-reimplemented operators (among other magic functions) to the wrapped array.

class mripy.paraproc.PooledCaller(pool_size=None, verbose=1)

Bases: object

Execute multiple command line programs, as well as python callables, asynchronously and parallelly across a pool of processes.

all_successful(jobs=None, verbose=None)

dispatch()

idss(total, batch_size=None)

run(cmd, *args, _depends=None, _retry=None, _dispatch=False, _error_pattern=None, _error_whitelist=None, _suppress_warning=False, _block=False, **kwargs)

Asynchronously run command or callable (queued execution, return immediately).

See subprocess.Popen() for more information about the arguments.

Multiple commands can be separated with “;” and executed sequentially within a single subprocess in linux/mac, only if shell=True.

Python callable can also be executed in parallel via multiprocessing. Note that although return values of the callable are retrieved via PIPE, sometimes it could be advantageous to directly save the computation results into a shared file (e.g., an HDF5 file), esp. when they’re large. In the later case, a proper lock mechanism via multiprocessing.Lock() is required.

Parameters:

• cmd (list, str, or callable) – Computation in command line programs is handled with subprocess. Computation in python callable is handled with multiprocessing.

• shell (bool) – If provided, must be a keyword argument. If shell is True, the command will be executed through the shell.

• *args – If cmd is a callable, *args are passed to the callable as its arguments.

• **kwargs – If cmd is a callable, **kwargs are passed to the callable as its keyword arguments. If cmd is a list or str, **kwargs are passed to subprocess.Popen().

• _depends (list) – A list of jobs (identified by their uuid) that have to be done before this job can be scheduled.

• _retry (int) – Number of retry before accepting failure (if detecting non-zero return code).

• _dispatch (bool) – Dispatch the job immediately, which will run in the background without blocking.

• _error_pattern (str) –

• _suppress_warning (bool) –

• _block (bool) – if True, call wait() internally and block.

Returns:

_uuid – The uuid of current job (which can be used as future jobs’ dependency)

Return type:

str

run1(cmd, *args, _error_pattern=None, _error_whitelist=None, _suppress_warning=False, **kwargs)

wait(pool_size=None, return_codes=False, return_jobs=False, raise_when_failed=True)

Wait for all jobs in the queue to finish.

Returns:

• return_values (list) – Return values of executed python callable. Always None for command.

• codes (list (only when return_codes=True)) – The return code of the child process for each job.

• jobs (list (only when return_jobs=True)) – Detailed information about each child process, including captured stdout and stderr.

class mripy.paraproc.SharedMemoryArray(dtype, shape, initializer=None, lock=True)

Bases: object

This class can be used as a usual np.ndarray, but its data buffer is allocated in shared memory (under Cached Files in memory monitor), and can be passed across processes without any data copy/duplication, even when write access happens (which is lock-synchronized).

The idea is to allocate memory using multiprocessing.Array, and access it from current or another process via a numpy.ndarray view, without actually copying the data. So it is both convenient and efficient when used with multiprocessing.

This implementation also demonstrates the power of composition + metaclass, as opposed to the canonical multiple inheritance.

dtype2ctypes = {<class 'bool'>: <class 'ctypes.c_bool'>, <class 'int'>: <class 'ctypes.c_long'>, <class 'float'>: <class 'ctypes.c_double'>, dtype('bool'): <class 'ctypes.c_bool'>, dtype('int64'): <class 'ctypes.c_long'>, dtype('int32'): <class 'ctypes.c_int'>, dtype('int16'): <class 'ctypes.c_short'>, dtype('int8'): <class 'ctypes.c_byte'>, dtype('uint64'): <class 'ctypes.c_ulong'>, dtype('uint32'): <class 'ctypes.c_uint'>, dtype('uint16'): <class 'ctypes.c_ushort'>, dtype('uint8'): <class 'ctypes.c_ubyte'>, dtype('float64'): <class 'ctypes.c_double'>, dtype('float32'): <class 'ctypes.c_float'>}

classmethod from_array(arr, lock=True)

Initialize a new shared-memory array with an existing array.

classmethod zeros(shape, dtype=<class 'float'>, lock=True)

Return a new array of given shape and dtype, filled with zeros.

This is the preferred usage, which avoids holding two copies of the potentially very large data simultaneously in the memory.

class mripy.paraproc.TeeOut(err=False, tee=True)

Bases: StringIO

write(s)

Write string to file.

Returns the number of characters written, which is always equal to the length of the string.

mripy.paraproc.add_metaclass(metaclass)

Class decorator for creating a class with a metaclass.

mripy.paraproc.check_output_for_errors(output, error_pattern=None, error_whitelist=None, verbose=1, label='')

User can skip error checking by setting error_pattern=’’

mripy.paraproc.check_output_for_goal(output, goal_pattern=None)

mripy.paraproc.cmd_for_disp(cmd)

Format cmd for printing.

Parameters:

cmd (str, list, or callable) –

mripy.paraproc.cmd_for_exec(cmd, shell=False)

Format cmd appropriately for execution according to whether shell=True.

Split a cmd string into a list, if shell=False. Join a cmd list into a string, if shell=True. Do nothing to callable.

Parameters:

• cmd (str, list, or callable) –

• shell (bool) –

mripy.paraproc.format_duration(duration, format='standard')

Format duration (in seconds) in a more human friendly way.

mripy.paraproc.run(cmd, check=True, error_pattern=None, error_whitelist=None, goal_pattern=None, shell=False, verbose=2)

Run an external command line.

This function is similar to subprocess.run introduced in Python 3.5, but provides a slightly simpler and perhaps more convenient API.

Parameters:

cmd (str or list) –

mripy.preprocess module

class mripy.preprocess.ANTsTransform(transforms, base_file=None, source_file=None)

Bases: Transform

apply(in_file, out_file, base_file=None, interp=None, **kwargs)

For volumes, forward transform (from input/moving to base/fixed)

apply_inverse(in_file, out_file, base_file=None, interp=None, **kwargs)

For volumes, inverse transform (from base/fixed to input/moving)

apply_inverse_to_points(in_file, out_file)

For list of points, inverse transform (from base/fixed to input/moving)

Parameters:

• in_file (*.csv file with “x,y,z,t” header line.) –

• out_file (*.csv file with “x,y,z,t” header line.) –

apply_inverse_to_xyz(xyz, convention='DICOM')

Parameters:

• xyz (Nx3 array) –

• convention ('DICOM' | 'NIFTI') –

apply_to_points(in_file, out_file)

For list of points, forward transform (from input/moving to base/fixed)

Parameters:

• in_file (*.csv file with “x,y,z,t” header line.) –

• out_file (*.csv file with “x,y,z,t” header line.) –

apply_to_xyz(xyz, convention='DICOM')

Parameters:

• xyz (Nx3 array) –

• convention ('DICOM' | 'NIFTI') –

classmethod from_align_ants(outputs)

class mripy.preprocess.Transform(transforms, base_file=None, source_file=None)

Bases: object

apply(in_file, out_file, base_file=None, interp=None, **kwargs)

For volumes, forward transform (from input/moving to base/fixed)

apply_inverse(in_file, out_file, base_file=None, interp=None, **kwargs)

For volumes, inverse transform (from base/fixed to input/moving)

classmethod from_json(fname, replace_path=False)

Parameters:

replace_path (bool) – Replace path of the transform files according to json file path

inverse()

rebase(base_file)

replace_path(p)

to_json(fname)

mripy.preprocess.afni2ants_affine(afni_affine, ants_affine)

mripy.preprocess.afni2ants_warp(afni_warp, ants_warp)

mripy.preprocess.align_S2E(base_file, suma_dir, out_file=None, **kwargs)

mripy.preprocess.align_anat(base_file, in_file, out_file, strip=None, N4=None, init_shift=None, init_rotate=None, init_xform=None, method=None, cost=None, n_params=None, interp=None, max_rotate=None, max_shift=None, emask=None, save_weight=None)

emaskfname

Mask to exclude from analysis.

mripy.preprocess.align_anat2epi(anat_file, epi_file, out_file, base_file=None, init_oblique=None, init_epi_rotate=None, init_anat_rotate=None)

Different init methods are mutual exclusive (i.e., at most one init method could be used at a time).

mripy.preprocess.align_ants(base_file, in_file, out_file, strip=None, base_mask=None, in_mask=None, base_mask_SyN=None, in_mask_SyN=None, init_transform=None, preset=None, n_jobs=None)

Nonlinearly align in_file to base_file using ANTs’ SyN method via antsRegistration.

The default preset can be used for both within- and cross-modality nonlinear registration, i.e., for both T1-to-T1 and EPI-to-T1 (more precisely, T2*-to-T1) deformable alignment.

Examples

1. Align MNI template to T1_al using default preset. Skullstrip T1_al.

   Apply inversed mask (1-mask) to MNI template at the nonlinear (SyN) stage. >>> prep.align_ants(“T1_al.nii”, “MNI152_2009_template.nii.gz”, “MNI_al.nii”, strip=”base”, in_mask_SyN=”mask.nii -I”)

2. Align T1_ns_al.nii to MNI template using “test” preset.

   For a quick test of the parameters in a few minutes. The result will not be good, but should not be weird. >>> prep.align_ants(“MNI152_2009_template.nii.gz”, “T1_ns_al.nii”, “T1_ns_MNI.nii”, preset=”test”)

Parameters:

• strip (str) – ‘base’ | ‘input’ | ‘both’

• preset (str) – None | ‘default’ | ‘test’ | ‘path/to/my/preset.json’ For production, just leave it as None or use the ‘default’ presest (est. time: 3 hr). For quick test, use the ‘test’ presest (est. time: 3 min).

Returns:

outputs –

outputs[‘transform’]ANTsTransform object

You can apply the forward or inverse transform to other volumes. E.g., >>> outputs[‘transform’].apply(in_file, out_file) >>> outputs[‘transform’].apply_inverse(in_file, out_file)

Return type:

dict

mripy.preprocess.align_center(base_file, in_file, out_file)

mripy.preprocess.align_epi(in_files, out_files, best_reverse=None, blip_results=None, blip_kws=None, volreg_kws=None, template=None, template_pool=None, template_candidate_runs=None, final_resample=True, final_res=None)

Parameters:

blip_results (list or dict) – E.g., 1) {‘warp_file’: ‘.volreg.warp.nii’, ‘blip_file’: ‘.volreg.blip.nii’} 2) {‘warf_file’: ‘epi01.for2mid.warp.nii’}

mripy.preprocess.all_finished(outputs)

mripy.preprocess.ants2afni_affine(ants_affine, afni_affine)

mripy.preprocess.ants2afni_warp(ants_warp, afni_warp)

mripy.preprocess.apply_ants(transforms, base_file, in_file, out_file, interp=None, dim=None, image_type=None, n_jobs=None)

Parameters:

• transforms (list of file names) – Online matrix inversion is supported as *_0GenericAffine.mat -I.

• base_file (str) – If None, apply transforms to point list using antsApplyTransformsToPoints, and in_file is expected to be a *.csv file. Otherwise, apply transforms to image using antsApplyTransforms.

• interp (str) – ‘LanczosWindowedSinc’, ‘NearestNeighbor’, ‘Linear’, ‘BSpline[<order=3>]’, etc.

• volumes (Note that for) –

• grid) (last transform applies first (pulling from base) –

:param : :param as in AFNI: :param as well as in ANTs command line.: :param However for point lists: :param FIRST transform applies first (pushing input points): :param : :param and INVERSE transforms should be used compared with volume case: :param as in ANTs.:

mripy.preprocess.apply_transforms(transforms, base_file, in_file, out_file, interp=None, res=None, save_xform=None)

Note that last transform applys first, as in AFNI. Transforms can be modified as “xform.aff12.1D -I” for inversion.

mripy.preprocess.assign_mp2rage_labels(T1s, dicom_dirs, dicom_ext='.IMA')

mripy.preprocess.average_anat(T1s, out_file, template_idx=0, T1s_ns=None, weight=None)

mripy.preprocess.blip_unwarp(forward_file, reverse_file, reverse_loc, out_file, PE_axis='AP', min_patch=None, pre_blur=None)

abin/unWarpEPI.py do this through: unwarp estimate (before tshift) > tshift > unwarp apply > motion correction > concat apply

mripy.preprocess.calculate_min_patch(fname, warp_res=10)

warp_resfloat

Effective “warp resolution” in mm.

mripy.preprocess.clusterize(in_file, out_file, cluster_size, neighbor=2)

mripy.preprocess.collect_perblock_stats(mask_file, stats_file, n_conditions, n_repeats, n_runs=None, n_stats=2, stat_index=0, skip_first=1)

mripy.preprocess.combine_affine_transforms(transforms, out_file=None)

Note that:

1. Each transform can be modified by -I, -S, -P, e.g., transforms=[‘T1_al.aff12.1D -I’]

2. An oneline transform could contain multiple lines, one for each time point.

mripy.preprocess.combine_censors(censor_files, out_file, method='intersect')

Combine multiple censor files: 1=keep, 0=discard

Parameters:

method (str, {'intersect', 'concatinate'}) –

mripy.preprocess.copy_S2E_mat(src, dst)

mripy.preprocess.copy_dset(in_file, out_file)

mripy.preprocess.correct_bias_field(in_file, out_file, n_jobs=None)

mripy.preprocess.correct_motion(base_file, in_file, out_file, algorithm='3dvolreg', mode='rigid')

algorithm : {‘3dvolreg’, ‘3dAllineate’}

mripy.preprocess.count_outliers(in_file, out_file=None, censor_file=None, th=0.1)

Compute outlier voxel fraction for each volume.

Censor (=0, excluded from further analysis) a volume when too many voxels (default>10%) within the automask are outliers.

mripy.preprocess.create_brain_mask(in_files, out_file, clfrac=None)

clfracfloat, between 0.1 and 0.9, [default=0.5]

A small ‘clip level fraction’ will tend to make the mask larger.

mripy.preprocess.create_extent_mask(transforms, base_file, in_files, out_file)

mripy.preprocess.create_hcp_retinotopic_atlas(subj_dir, suma='SUMA', NIFTI=True)

Create HCP retinotopic atlas (benson14 template) using docker, and convert the volume and surface datasets into SUMA format.

mripy.preprocess.create_mp2rage_SNR_mask(T1s, out_file)

Need to call prep.assign_mp2rage_labels() first.

mripy.preprocess.create_suma_dir(subj_dir, suma_dir=None, NIFTI=True)

Notes

ABOUT NIFTI

If NIFTI=True, @SUMA_Make_Spec_FS will use mri_convert to convert $surf_dir/orig.mgz (which is averaged and hires) to SurfVol.nii. Volumes will be in .nii and surfaces will be in .gii format. This is the preferable way to go.

If NIFTI=False, @SUMA_Make_Spec_FS will use to3d to convert COR- files to SurfVol+orig.HEAD which is 1x1x1. Volumes will be in +orig.HEAD and surfaces will be in .asc format. This is only provided for backward compatibility.

mripy.preprocess.create_suma_script(spec_file, surf_vol, out_file, use_relpath=False)

mripy.preprocess.create_vessel_mask_BOLD(beta, out_file, th=10)

mripy.preprocess.create_vessel_mask_EPI(mean_epi, ribbon, out_file, corr_file=None, th=None)

Find vessel voxels based on bias-corrected EPI intensity (Kay’s method). The result may also highlight misalignment between EPI and T1, as well as errors in pial surface reconstruction.

Parameters:

• mean_epi (str) –

• ribbon (str) – Mask for gray matter which is aligned with mean_epi.

• th (float or 'auto', default 0.75 (as in [1])) –

References

[1] Kay, K., Jamison, K. W., Vizioli, L., Zhang, R., Margalit, E., & Ugurbil, K. (2019). A critical assessment of data quality and venous effects in sub-millimeter fMRI. NeuroImage, 189, 847–869.

mripy.preprocess.create_vessel_mask_PDGRE(PD, GRE, out_file, PD_div_GRE=None, ath=300, rth=5, nath=100, nrth=-5, base_file=None, strip_PD=True, strip_GRE=True)

Find vessel voxels based on PD/GRE contrast.

Parameters:

• PD (str, aligned with EXP) –

• GRE (str, aligned with PD) –

• ath (float) – Global absolute threshold for PD/GRE

• rth (float) – Local (within neighborhood) relative deviation for PD/GRE

mripy.preprocess.deoblique(in_file, out_file=None, template=None)

mripy.preprocess.detrend(motion_file, in_file, out_file, censor=True, motion_th=0.3, censor_file=None, regressor_file=None)

mripy.preprocess.find_best_reverse(seq_info, forward='func', reverse='reverse', for_out='epi{run}.nii', rev_out='reverse{run}.nii', return_single_best=False)

Find best forward-reverse pair (as a mapping dict) according to temporal proximity, assuming that head motion is usually smaller if the two images are closer in time.

seq_info : DataFrame

mripy.preprocess.find_epi_dark_voxel_threshold(v, method='Liu2020')

Find EPI intensity threshold for blood vessels on surface dset,

produced by Surface.vol2surf(func=’min’).

Parameters:

method (str, 'Liu2020' | 'Marquardt2018' | 'Kay2019') – Surprisingly, the three methods usually produce very similar results. The mixture of Gaussian method (‘Kay2019’) is robust with ‘N4’ or ‘Kay2019’ unifized EPI data (using preprocess.unifize_epi), but may fail with ‘3dUnifize’ or raw EPI data.

mripy.preprocess.fs_recon(T1s, out_dir, T2=None, FLAIR=None, NIFTI=True, hires=True, fs_ver=None, V1=True, HCP_atlas=True, n_jobs=None)

Parameters:

• T1s (list of str | 'brainmask_edit' | 'wm_edit') –

• fs_ver ({'v6', 'v6.hcp', 'skip'}) –

mripy.preprocess.glm(in_files, out_file, design, model='BLOCK', contrasts=None, TR=None, pick_runs=None, motion_files=None, censor=True, motion_th=0.3, censor_files=None, mask=None, regressor_file=None, poly=None, fitts=True, errts=True, REML=True, perblock=False, FDR=None, check_TR=True)

Parameters:

• design (OrderedDict(L=('../stimuli/L.txt', 24), R="../stimuli/R.txt 'BLOCK(24,1)'")) –

• model (str) –

• contrasts (OrderedDict([('L+R', '+0.5*L +0.5*R'), ...])) –

• pick_runs (array-like) – Indices of the runs you want to pick for analysis (start from zero)

Examples

1. Canonical GLM

   design = OrderedDict() design[‘L’] = (‘../stimuli/L_localizer.txt’, 24) design[‘R’] = (‘../stimuli/R_localizer.txt’, 24) contrasts = OrderedDict() contrasts[‘L+R’] = ‘+0.5*L +0.5*R’ contrasts[‘L-R’] = ‘+L -R’ model = ‘BLOCK’

2. FIR estimation

   design = OrderedDict() design[‘L’] = f”{stim_dir}/L_{condition}_deconv.txt ‘CSPLINzero(0,24,11)’” design[‘R’] = f”{stim_dir}/R_{condition}_deconv.txt ‘CSPLINzero(0,24,11)’” contrasts = None model = ‘BLOCK’ # This is not used for TENT/CSPLIN

mripy.preprocess.irregular_resample(transforms, xyz, in_file, order=3)

Parameters:

• transforms (list of str) – In order to project raw EPI in volume (in RAI) to surface coordinates (in LPI aka RAS+), i.e., surf_xyz = [dicom2nifti, exp2surf, volreg2template, unwarp2template, ijk2xyz] @ ijk, what we really do is mapping backward from xyz (assumed already in RAI): xyz -> E2A storing inv(exp2surf) -> … -> inv(MAT).

• xyz (Nx3 array, assumed in DICOM RAI as with AFNI volumes.) – Note that FreeSurfur surface vertices are in NIFTI LPI aka RAS+, whereas AFNI uses DICOM RAI internally.

mripy.preprocess.is_affine_transform(fname)

mripy.preprocess.manual_transform(in_file, out_file, shift=None, rotate=None, scale=None, shear=None, interp=None)

shift : [x, y, z] in mm rotate : [I, R, A] in degrees (i.e., -z, -x, -y axes), right-hand rule

For example, ‘shift=[1,0,0]’ shifts in_file to the left by 1 mm (not necessarily 1 voxel). The acutal applied param will be negated. Because if the source is shifted to the left compared with the base, the resulted param will be x=1 indicating that source is shifted to the left, and apply that param will actually cause the source to shift rightward by 1 to match the base. So to shift in_file leftward, the actual param will be x=-1.

The matrix is specified in DICOM-ordered (RAI) coordinates (x=-R+L,y=-A+P,z=-I+S). By default the shift is applied AFTER matrix transform, as in augmented 4x4 affine. By default the shear matrix is LOWER triangle.

For ‘-1Dmatrix_save’ and ‘-1Dmatrix_apply’, the matrix specifies coordinate transformation from base to source DICOM coordinates. In other words, with the estimated matrix at hand, you are ready to build the transformed image by filling the base grid with source data.

Refer to “DEFINITION OF AFFINE TRANSFORMATION PARAMETERS”@3dAllineate for details.

mripy.preprocess.nudge_cmd2mat(nudge_cmd, in_file, return_inverse=False)

Refer to “Example 4”@SUMA_AlignToExperiment for details.

The 1st return (MAT) is what you want to write as ‘src2base.aff12.1D’ and use apply_transforms to get the same effect as your manual nudge. Remember, the AFNI way is to store inverse matrix under forward name (to “pull” data from src). So the mathematical map from moveable to tempalte is in the 2nd return (INV).

mripy.preprocess.parse_seq_info(raw_dir, return_dataframe=True)

Assuming your rawdata folder hierarchy is like:

raw_fmri ├── func01 ├── func02 ├── reverse01 ├── reverse02 ├── T1

mripy.preprocess.prep_mp2rage(dicom_dirs, out_file='T1.nii', unwarp=False, dicom_ext='.IMA')

Convert dicom files and remove the noise pattern outside the brain.

dicom_dirslist or str

A list of dicom file folders, e.g., [‘T101’, ‘T102’, …], or a glob pattern like ‘raw_fmri/T1??’

mripy.preprocess.prep_mp2rages(data_dir, sessions=None, subdir_pattern='T1??', unwarp=True, **kwargs)

mripy.preprocess.resample(in_file, out_file, res=None, base_file=None, interp=None)

mripy.preprocess.resample_to_surface(transforms, surfaces, in_file, out_files=None, mask_file=None, n_jobs=1, **kwargs)

Parameters:

mask_file (str, dict) – Surface mask. Either a file name or dict(lh=’lh.mask.niml.dset’, rh=’rh.mask.niml.dset’). This will generate a partial surface dataset.

Examples

resample_to_surface(transforms=[f’SurfVol_Alnd_Exp.E2A.1D’, f’epi{run}.volreg.aff12.1D’, f’epi{run}.volreg.warp.nii’],

surfaces=[f’{suma_dir}/lh.pial.asc’, f’{suma_dir}/rh.smoothwm.asc’], in_file=f’epi{run}.tshift.nii’)

mripy.preprocess.retrieve_mp2rage_labels(dicom_dirs, dicom_ext='.IMA')

Retrieve mp2rage subvolume labels like UNI, ND, etc.

Parameters:

dicom_dirs (list or str) – A list of dicom file folders, e.g., [‘T101’, ‘T102’, …], or a glob pattern like ‘raw_fmri/T1??’

Returns:

label2dicom_dir – label -> (index, dicom_dir)

Return type:

OrderedDict

mripy.preprocess.scale(in_file, out_file, mask_file=None, dtype=None)

mripy.preprocess.skullstrip(in_file, out_file=None)

mripy.preprocess.unifize_epi(in_file, out_file, method='N4', ribbon=None)

Remove spatial trend/inhomogeneity in (mean) EPI data to better identify vessels using find_epi_dark_voxel_threshold().

The default ‘N4’ method requires you have ANTs installed.

Only ‘Kay2019’ method need an additional cortical ribbon mask that is aligned with the EPI volume.

Parameters:

method (str, 'N4' | 'Kay2019' | '3dUnifize') –

mripy.preprocess.zscore(in_file, out_file)

mripy.six module

Utilities for writing code that runs on Python 2 and 3

class mripy.six.Module_six_moves_urllib(name, doc=None)

Bases: module

Create a six.moves.urllib namespace that resembles the Python 3 namespace

error = <module 'mripy.six.moves.urllib.error'>

parse = <module 'mripy.six.moves.urllib_parse'>

request = <module 'mripy.six.moves.urllib.request'>

response = <module 'mripy.six.moves.urllib.response'>

robotparser = <module 'mripy.six.moves.urllib.robotparser'>

class mripy.six.Module_six_moves_urllib_error(name)

Bases: _LazyModule

Lazy loading of moved objects in six.moves.urllib_error

ContentTooShortError

HTTPError

URLError

class mripy.six.Module_six_moves_urllib_parse(name)

Bases: _LazyModule

Lazy loading of moved objects in six.moves.urllib_parse

ParseResult

SplitResult

parse_qs

parse_qsl

quote

quote_plus

splitquery

splittag

splituser

splitvalue

unquote

unquote_plus

unquote_to_bytes

urldefrag

urlencode

urljoin

urlparse

urlsplit

urlunparse

urlunsplit

uses_fragment

uses_netloc

uses_params

uses_query

uses_relative

class mripy.six.Module_six_moves_urllib_request(name)

Bases: _LazyModule

Lazy loading of moved objects in six.moves.urllib_request

AbstractBasicAuthHandler

AbstractDigestAuthHandler

BaseHandler

CacheFTPHandler

FTPHandler

FancyURLopener

FileHandler

HTTPBasicAuthHandler

HTTPCookieProcessor

HTTPDefaultErrorHandler

HTTPDigestAuthHandler

HTTPErrorProcessor

HTTPHandler

HTTPPasswordMgr

HTTPPasswordMgrWithDefaultRealm

HTTPRedirectHandler

HTTPSHandler

OpenerDirector

ProxyBasicAuthHandler

ProxyDigestAuthHandler

ProxyHandler

Request

URLopener

UnknownHandler

build_opener

getproxies

install_opener

pathname2url

proxy_bypass

url2pathname

urlcleanup

urlopen

urlretrieve

class mripy.six.Module_six_moves_urllib_response(name)

Bases: _LazyModule

Lazy loading of moved objects in six.moves.urllib_response

addbase

addclosehook

addinfo

addinfourl

class mripy.six.Module_six_moves_urllib_robotparser(name)

Bases: _LazyModule

Lazy loading of moved objects in six.moves.urllib_robotparser

RobotFileParser

class mripy.six.MovedAttribute(name, old_mod, new_mod, old_attr=None, new_attr=None)

Bases: _LazyDescr

class mripy.six.MovedModule(name, old, new=None)

Bases: _LazyDescr

mripy.six.add_metaclass(metaclass)

Class decorator for creating a class with a metaclass.

mripy.six.add_move(move)

Add an item to six.moves.

mripy.six.assertCountEqual(self, *args, **kwargs)

mripy.six.assertRaisesRegex(self, *args, **kwargs)

mripy.six.assertRegex(self, *args, **kwargs)

mripy.six.b(s)

Byte literal

mripy.six.create_unbound_method(func, cls)

mripy.six.get_unbound_function(unbound)

Get the function out of a possibly unbound function

mripy.six.int2byte()

S.pack(v1, v2, …) -> bytes

Return a bytes object containing values v1, v2, … packed according to the format string S.format. See help(struct) for more on format strings.

mripy.six.iteritems(d, **kw)

Return an iterator over the (key, value) pairs of a dictionary.

mripy.six.iterkeys(d, **kw)

Return an iterator over the keys of a dictionary.

mripy.six.iterlists(d, **kw)

Return an iterator over the (key, [values]) pairs of a dictionary.

mripy.six.itervalues(d, **kw)

Return an iterator over the values of a dictionary.

mripy.six.python_2_unicode_compatible(klass)

A decorator that defines __unicode__ and __str__ methods under Python 2. Under Python 3 it does nothing.

To support Python 2 and 3 with a single code base, define a __str__ method returning text and apply this decorator to the class.

mripy.six.raise_from(value, from_value)

mripy.six.remove_move(name)

Remove item from six.moves.

mripy.six.reraise(tp, value, tb=None)

Reraise an exception.

mripy.six.u(s)

Text literal

mripy.six.with_metaclass(meta, *bases)

Create a base class with a metaclass.

mripy.surface module

mripy.timecourse module

class mripy.timecourse.Attributes(shape)

Bases: object

add(name, value, axis)

classmethod concatinate(attributes_list, axis)

drop(name)

drop_all_with_axis(axis)

classmethod from_dict(d)

names_with_axis(axis)

pick(index, axis)

property shape

to_dict()

class mripy.timecourse.Epochs(raw, events, event_id=None, tmin=-5, tmax=15, baseline=(-2, 0), dt=0.1, interp='linear', hamm=None, conditions=None)

Bases: Savable, object

add_event_attr(name, value)

add_feature_attr(name, value)

aggregate(event=True, feature=False, time=False, method=<function nanmean>, keepdims=<no value>, return_index=False)

apply_baseline(baseline)

average(feature=True, time=False, method=<function nanmean>, error='bootstrap', ci=95, n_boot=1000, condition=None)

Average data over event (and optionally feature and/or time) dimensions, and return an Evoked object.

copy()

drop_events(ids)

classmethod from_array(data, TR=None, tmin=None, baseline=(-2, 0), events=None, event_id=None, conditions=None)

classmethod from_dict(d)

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish, or a proper initialization requires additional operations performed in __init__ but omitted here.

property n_events

property n_features

property n_times

pick(event=None, feature=None, time=None)

plot(hue=None, style=None, row=None, col=None, hue_order=None, style_order=None, row_order=None, col_order=None, palette=None, dashes=None, figsize=None, legend=True, bbox_to_anchor=None, subplots_kws=None, average_kws=None, axs=None, **kwargs)

property shape

summary(event=False, feature=True, time=False, method=<function nanmean>, attributes=None)

Summary data as a pandas DataFrame.

to_dict()

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish (i.e., other than int, float, str, list, tuple, dict, numpy array, pandas dataframe).

transform(feature_name, feature_values, transformer)

Transform the data into a new Epochs object. E.g., epochs.transform(‘depth’, depth, tc.wcutter(depth, linspace(0, 1, 20), win_size=0.2, win_func=’gaussian’, exclude_outer=True))

class mripy.timecourse.Evoked(info, data, nave, times, error=None, error_type=None, condition=None)

Bases: object

plot(color=None, error=True, info=True, error_kws=None, show_n='info', **kwargs)

property shape

class mripy.timecourse.Raw(fname, mask=None, TR=None)

Bases: Savable, object

property TR

copy()

classmethod from_array(data, TR)

data : 2D array, [n_features, n_times] TR : in sec

classmethod from_dict(d)

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish, or a proper initialization requires additional operations performed in __init__ but omitted here.

property n_features

property n_times

plot(events=None, event_id=None, color=None, palette=None, figsize=None, event_kws=None, **kwargs)

property shape

to_dict()

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish (i.e., other than int, float, str, list, tuple, dict, numpy array, pandas dataframe).

class mripy.timecourse.RawCache(fnames, mask, TR=None, cache_file=None, force_redo=False)

Bases: Savable, object

classmethod from_dict(d)

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish, or a proper initialization requires additional operations performed in __init__ but omitted here.

get_epochs(mask, events, event_id, ids=None, cache_file=None, **kwargs)

get_raws(mask, ids=None)

property n_runs

subset(mask, cache_file=None)

to_dict()

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish (i.e., other than int, float, str, list, tuple, dict, numpy array, pandas dataframe).

mripy.timecourse.concatinate_epochs(epochs_list, axis=0)

mripy.timecourse.convolve_HRF(starts, lens, TR=2, scan_time=None, HRF=None)

mripy.timecourse.create_ERP(t, x, events, tmin=-8, tmax=16, dt=0.1, baseline=[-2, 0], interp='linear')

t : time for each data point (can be non-contiguous) x : [feature, time] events : event onset time (can be on non-integer time point)

mripy.timecourse.create_base_corr_func(times, baseline=None, method=None)

Parameters:

• time (array-like) – Sampling time for your data.

• baseline ('none', 'all', or (tmin, tmax)) – Baseline time interval.

mripy.timecourse.create_ideal(stimuli, lens, **kwargs)

Parameters:

stimuli (list of fname) –

mripy.timecourse.create_times(tmin, tmax, dt)

mripy.timecourse.cut(data, val, bins, **kwargs)

mripy.timecourse.cutter(val, bins)

mripy.timecourse.events_from_dataframe(df, time, conditions, duration=None, run=None, n_runs=None, event_id=None)

Parameters:

• df (dataframe) –

• time (str) –

• conditions (list (names) or OrderedDict (name->levels)) –

• duration (str) –

• run (str) –

• event_id (dict (name->index)) –

Returns:

• events (array)

• event_id (dict (name->index))

mripy.timecourse.events_to_dataframe(events_list, event_id, conditions)

mripy.timecourse.first_peak(t, x, after=2, trough='before', sg_win=41, sg_order=3, visualize=False, ax=None)

Difference between first peak after after and first trough before it.

mripy.timecourse.get_HRF(name, TR=2, duration=32, normalize='area', **kwargs)

Get hemodynamic response function.

Parameters:

• name (str) –

  Name of the hemodynamic response function model. - ‘canonical’: difference of two gamma distribution functions [1]

  hrf = gamma.pdf(t, a1, scale=b1) - c*stats.gamma.pdf(t, a2, scale=b2) Acceptable kwargs are a1=6, b1=1, a2=16, b2=1, c=1/6. a1, a2 and b1, b2 are the shape and scale parameter of the two gammas, and c is the relative amplitude of the 2nd (negative) gamma. The default result matches the default spm_hrf result in SPM12.

  • ’SPM’: difference of two gamma, parameterized as in spm_hrf (SPM12) [2]

    Acceptable kwargs are P=[6,16,1,1,6,0,32], T=16. See spm_hrf for more. The conversion between ‘canonical’ and ‘SPM’ parameterizations: In scipy convention, larger scale parameter inplies wider distribution. So for gamma distribution parameterized as

    beta^alpha * x^(alpha-1) * np.exp(-beta*x) / special.gamma(alpha)

    we have scale = 1/beta, i.e., beta = 1/scale And since u (=0,1,2,…) is x * 1/dt, to keep the scale, beta should be multiplied by dt. So we have beta = dt/scale, and p3, p4 are actually the scale (disperse) parameter of gamma. p1, p2 are (standardized) shape parameters in units of scale parameter, p1 = alpha*p3, i.e., alpha = p1/p3. Note that the mean of the Gamma distribution is shape*scale and its mode is (shape-1)*scale. This means the positive and negative peaks of the hrf are approximately p1-1 and p2-1. Finally, p5 is the response/undershoot ratio, i.e., 1/c as in the ‘canonical’ parameterization. In summary, the following two expressions are equivalent:

    get_HRF(‘canonical’, a1=a1, a2=a2, b1=b1, b2=b2, c=c) get_HRF(‘SPM’, P=[a1*b1, a2*b2, b1, b2, 1/c, 0, 32])

  • ’GAM’: one gamma (no undershoot), from AFNI 3dDeconvolve [3]

    hrf = (t/(p*q))^p * exp(p-t/q)

    Acceptable kwargs are p=8.6, q=0.547. The peak of ‘GAM(p,q)’ is at time p*q after the stimulus. The FWHM is about 2.35*sqrt(p)*q.

  • ’TDM-early’ and ‘TDM-late’: two-gamma model for the early and late

    HRF components in temporal decomposition through manifold fitting [4] Note that the two HRF components should be fitted together in a GLM, although they have high correlation (~0.8).

• TR (float) – Sampling duration in seconds.

• duration (float) – Length of kernel in seconds.

• normalize (str) –

  • ‘area’ : normalize the (signed) area under the curve to one.

  • ’height’ : normalize the max respones to one.

• **kwargs – Additional parameters adjustable for the HRF model. See name for details.

Returns:

• t (array) – The time vector [0:TR:duration).

• hrf (array) – The HRF evaluated at each point of the time vector.

References

[1] Lindquist, M. A., Meng Loh, J., Atlas, L. Y., & Wager, T. D. (2009).

Modeling the hemodynamic response function in fMRI: Efficiency, bias and mis-modeling. NeuroImage, 45(1, Supplement 1), S187–S198.

[2] https://github.com/spm/spm/blob/main/spm_hrf.m [3] https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dDeconvolve.html [4] Kay, K., Jamison, K. W., Zhang, R.-Y., & Uğurbil, K. (2020). A temporal

decomposition method for identifying venous effects in task-based fMRI. Nature Methods.

Notes

Gamma distribution function and gamma function are different things.

mripy.timecourse.group_epochs(epochs_list, pool_feature=True)

mripy.timecourse.qcut(data, val, q, **kwargs)

mripy.timecourse.qcutter(val, q)

mripy.timecourse.read_events(event_files)

Read events from AFNI style (each row is a run, and each element is an occurance) stimulus timing files.

Parameters:

event_files (dict) – e.g., OrderedDict((‘Physical/Left’, ‘../stimuli/phy_left.txt’), (‘Physical/Right’, ‘../stimuli/phy_right.txt’), …)

Returns:

• events_list (list of n_events-by-3 arrays) – The three columns are [start_time(sec), reserved, event_id(int)]

• event_id (dict)

mripy.timecourse.wcut(data, val, v, win_size, win_func=None, exclude_outer=False, **kwargs)

mripy.timecourse.wcutter(val, v, win_size, win_func=None, exclude_outer=False)

mripy.utils module

class mripy.utils.CacheManager(persistent_file=None, ignore_init_run=True)

Bases: object

exists(fname, watch_files=None, force_redo=False, **kwargs)

kwargs_updated(fname, kwargs)

load_contexts()

save_contexts()

watch_files_updated(fname, watch_files)

class mripy.utils.FilenameManager(fmt, **kwargs)

Bases: object

classmethod fmt2kws(fmt, **kwargs)

format(fmt=None, keepdims=False, **kwargs)

classmethod from_glob(fmt, **kwargs)

glob(**kwargs)

parse(files, multi_value='list')

files : list of filenames multi_value : {‘wildcard’, ‘list’}

class mripy.utils.ParallelCaller

Bases: object

check_call(cmd, **kwargs)

Asynchronous check_call (launch and return immediately). Multiple commands can be stitched sequentially with “;” in linux/mac only if shell=True.

wait()

Wail for all parallel processes to finish (= wait for the slowest one).

class mripy.utils.Savable

Bases: object

classmethod from_dict(d)

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish, or a proper initialization requires additional operations performed in __init__ but omitted here.

classmethod load(fname)

save(fname)

to_dict()

Subclasses are expected to override this method if the data type of some attributes are not supported by deepdish (i.e., other than int, float, str, list, tuple, dict, numpy array, pandas dataframe).

class mripy.utils.Savable2

Bases: object

load(fname)

save(fname)

mripy.utils.cd(d)

mripy.utils.contain_wildcard(fname)

mripy.utils.exists(fname, force_redo=False)

Lazy evaluation the body only if fname not already exists.

mripy.utils.expand_index_list(index_list, format=None)

Expand a list of indices like: 1 3-5 7:10:2 8..10(2)

mripy.utils.factorize(n)

mripy.utils.fname_with_ext(fname, ext)

mripy.utils.has_N4()

mripy.utils.has_ants()

mripy.utils.has_hcp_retino_docker()

mripy.utils.iterable(x)

mripy.utils.parallel_1D(cmd, in_file, prefix, n_jobs=1, combine_output=True, **kwargs)

Parameters:

cmd (str) – The command for parallel execution, must contain two placeholders {prefix} and {in_file}, which will be substituded with ‘’.format(). As expected, other {} must be escaped as {{}}. For example, >>> parallel_1D(“3dTcat -prefix {prefix} -overwrite {in_file}’{{0:9}}’”, ‘xyz_list.1D’, ‘test’, n_jobs=4)

mripy.utils.parallel_3D(cmd, in_file, prefix, n_jobs=1, schema=None, fname_mapper=None, combine_output=True, **kwargs)

Parameters:

fname_mapper (dict or callable) –

mripy.utils.select_and_replace_affix(glob_pattern, old_affix, new_affix)

mripy.utils.temp_folder(parent=None)

mripy.utils.temp_prefix(prefix='tmp_', n=4, suffix='.')

mripy.vis module

mripy.vis.draw_color_circle(cmap, res=512)

Draw circular colorbar (color index from 0 to 1) clockwise from 0 to 12 o’clock

mripy.vis.get_color_list(cmap)

This function is still wrong!

mripy.vis.plot_volreg(dfiles, convention=None)

The six columns from left to right are:

• For 3dvolreg: roll(z), pitch(x), yaw(y), dS(z), dL(x), dP(y) Note that the interpretation for roll and yaw in afni is opposite from SPM (and common sense…) This program will follow the SPM convention.

• For 3dAllineate shift_rotate: x-shift y-shift z-shift$ z-angle x-angle$ y-angle$ Note that the dollar signs in the end indicate parameters that are fixed.

Module contents