/usr/lib/python2.7/dist-packages/surfer/io.py is in python-surfer 0.7-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 | import os
from tempfile import mktemp
from subprocess import Popen, PIPE, check_output
import gzip
import numpy as np
import nibabel as nib
from nibabel.spatialimages import ImageFileError
from .utils import verbose
import logging
logger = logging.getLogger('surfer')
def read_scalar_data(filepath):
"""Load in scalar data from an image.
Parameters
----------
filepath : str
path to scalar data file
Returns
-------
scalar_data : numpy array
flat numpy array of scalar data
"""
try:
scalar_data = nib.load(filepath).get_data()
scalar_data = np.ravel(scalar_data, order="F")
return scalar_data
except ImageFileError:
ext = os.path.splitext(filepath)[1]
if ext == ".mgz":
openfile = gzip.open
elif ext == ".mgh":
openfile = open
else:
raise ValueError("Scalar file format must be readable "
"by Nibabel or .mg{hz} format")
fobj = openfile(filepath, "rb")
# We have to use np.fromstring here as gzip fileobjects don't work
# with np.fromfile; same goes for try/finally instead of with statement
try:
v = np.fromstring(fobj.read(4), ">i4")[0]
if v != 1:
# I don't actually know what versions this code will read, so to be
# on the safe side, let's only let version 1 in for now.
# Scalar data might also be in curv format (e.g. lh.thickness)
# in which case the first item in the file is a magic number.
raise NotImplementedError("Scalar data file version not supported")
ndim1 = np.fromstring(fobj.read(4), ">i4")[0]
ndim2 = np.fromstring(fobj.read(4), ">i4")[0]
ndim3 = np.fromstring(fobj.read(4), ">i4")[0]
nframes = np.fromstring(fobj.read(4), ">i4")[0]
datatype = np.fromstring(fobj.read(4), ">i4")[0]
# Set the number of bytes per voxel and numpy data type according to
# FS codes
databytes, typecode = {0: (1, ">i1"), 1: (4, ">i4"), 3: (4, ">f4"),
4: (2, ">h")}[datatype]
# Ignore the rest of the header here, just seek to the data
fobj.seek(284)
nbytes = ndim1 * ndim2 * ndim3 * nframes * databytes
# Read in all the data, keep it in flat representation
# (is this ever a problem?)
scalar_data = np.fromstring(fobj.read(nbytes), typecode)
finally:
fobj.close()
return scalar_data
def read_stc(filepath):
"""Read an STC file from the MNE package
STC files contain activations or source reconstructions
obtained from EEG and MEG data.
Parameters
----------
filepath: string
Path to STC file
Returns
-------
data: dict
The STC structure. It has the following keys:
tmin The first time point of the data in seconds
tstep Time between frames in seconds
vertices vertex indices (0 based)
data The data matrix (nvert * ntime)
"""
fid = open(filepath, 'rb')
stc = dict()
fid.seek(0, 2) # go to end of file
file_length = fid.tell()
fid.seek(0, 0) # go to beginning of file
# read tmin in ms
stc['tmin'] = float(np.fromfile(fid, dtype=">f4", count=1))
stc['tmin'] /= 1000.0
# read sampling rate in ms
stc['tstep'] = float(np.fromfile(fid, dtype=">f4", count=1))
stc['tstep'] /= 1000.0
# read number of vertices/sources
vertices_n = int(np.fromfile(fid, dtype=">u4", count=1))
# read the source vector
stc['vertices'] = np.fromfile(fid, dtype=">u4", count=vertices_n)
# read the number of timepts
data_n = int(np.fromfile(fid, dtype=">u4", count=1))
if ((file_length / 4 - 4 - vertices_n) % (data_n * vertices_n)) != 0:
raise ValueError('incorrect stc file size')
# read the data matrix
stc['data'] = np.fromfile(fid, dtype=">f4", count=vertices_n * data_n)
stc['data'] = stc['data'].reshape([data_n, vertices_n]).T
# close the file
fid.close()
return stc
@verbose
def project_volume_data(filepath, hemi, reg_file=None, subject_id=None,
projmeth="frac", projsum="avg", projarg=[0, 1, .1],
surf="white", smooth_fwhm=3, mask_label=None,
target_subject=None, verbose=None):
"""Sample MRI volume onto cortical manifold.
Note: this requires Freesurfer to be installed with correct
SUBJECTS_DIR definition (it uses mri_vol2surf internally).
Parameters
----------
filepath : string
Volume file to resample (equivalent to --mov)
hemi : [lh, rh]
Hemisphere target
reg_file : string
Path to TKreg style affine matrix file
subject_id : string
Use if file is in register with subject's orig.mgz
projmeth : [frac, dist]
Projection arg should be understood as fraction of cortical
thickness or as an absolute distance (in mm)
projsum : [avg, max, point]
Average over projection samples, take max, or take point sample
projarg : single float or sequence of three floats
Single float for point sample, sequence for avg/max specifying
start, stop, and step
surf : string
Target surface
smooth_fwhm : float
FWHM of surface-based smoothing to apply; 0 skips smoothing
mask_label : string
Path to label file to constrain projection; otherwise uses cortex
target_subject : string
Subject to warp data to in surface space after projection
verbose : bool, str, int, or None
If not None, override default verbose level (see surfer.verbose).
"""
env = os.environ
if 'FREESURFER_HOME' not in env:
raise RuntimeError('FreeSurfer environment not defined. Define the '
'FREESURFER_HOME environment variable.')
# Run FreeSurferEnv.sh if not most recent script to set PATH
if not env['PATH'].startswith(os.path.join(env['FREESURFER_HOME'], 'bin')):
cmd = ['bash', '-c', 'source {} && env'.format(
os.path.join(env['FREESURFER_HOME'], 'FreeSurferEnv.sh'))]
envout = check_output(cmd)
env = dict(line.split('=', 1)
for line in envout.decode('utf-8').split('\n')
if '=' in line)
# Set the basic commands
cmd_list = ["mri_vol2surf",
"--mov", filepath,
"--hemi", hemi,
"--surf", surf]
# Specify the affine registration
if reg_file is not None:
cmd_list.extend(["--reg", reg_file])
elif subject_id is not None:
cmd_list.extend(["--regheader", subject_id])
else:
raise ValueError("Must specify reg_file or subject_id")
# Specify the projection
proj_flag = "--proj" + projmeth
if projsum != "point":
proj_flag += "-"
proj_flag += projsum
if hasattr(projarg, "__iter__"):
proj_arg = list(map(str, projarg))
else:
proj_arg = [str(projarg)]
cmd_list.extend([proj_flag] + proj_arg)
# Set misc args
if smooth_fwhm:
cmd_list.extend(["--surf-fwhm", str(smooth_fwhm)])
if mask_label is not None:
cmd_list.extend(["--mask", mask_label])
if target_subject is not None:
cmd_list.extend(["--trgsubject", target_subject])
# Execute the command
out_file = mktemp(prefix="pysurfer-v2s", suffix='.mgz')
cmd_list.extend(["--o", out_file])
logger.info(" ".join(cmd_list))
p = Popen(cmd_list, stdout=PIPE, stderr=PIPE, env=env)
stdout, stderr = p.communicate()
out = p.returncode
if out:
raise RuntimeError(("mri_vol2surf command failed "
"with command-line: ") + " ".join(cmd_list))
# Read in the data
surf_data = read_scalar_data(out_file)
os.remove(out_file)
return surf_data
|