Hacker’s Guide to SPM12

MR image file handling

In SPM, MR image files are handled in the form of the volume struct. If associated with a NIfTI image file, it has a field private that contains a nifti object, which in turn has a field dat that contains a file_array object.¹

Volume struct

A volume struct represents a 3D MR image. Essential fields are:

fname	filename of the image, from D.fname
dim	size of the volume (voxels × voxels × voxels), derived from D.dim
mat	4 × 4 affine transformation matrix mapping from voxel indices to real-world coordinates, from N.mat
dt	data format
(1)	data type (see spm_type), from property D.dtype
(2)	endianess (0: little, 1: big), from property D.be
pinfo	linear transformation from raw values, and file offset; if two-dimensional, values are plane-specific
(1, :)	slope, from D.scl_slope
(2, :)	intercept, from D.scl_inter
(3, :)	offset into file to (volume/plane) data start (bytes), derived from D.offset

These further fields are normally present:

n	volume indices
(1)	sequence index (for 4D files, else 1)
(2)	vector element index (for 5D files, else 1)
descrip	description of the image, from N.descrip
private	underlying nifti object (= N)

mat (and mat0 below) is an affine transformation matrix, i.e. it has the augmented form

[   A   b
  0 0 0 1 ]

where A is a 3 × 3 linear transformation matrix and b is a translation. Voxel indices (starting from 1) are transformed into real-world coordinates (in mm) by

[x y z 1] .' ← mat * [i j k 1] .'

A NIfTI image can contain up to 7-dimensional data. The first three dimensions are interpreted as spatial (spanning a volume), the 4th as temporal or sequential, the 5th to 7th as indexing vector, matrix, or tensor elements.

SPM supports referencing volumes from a multi-volume file in up to two dimensions using the notation filename,#,#, where the # stands for indices referring to the 4th and 5th dimension. filename,# is equivalent to filename,#,1. filename is equivalent to a sequence of volume references, filename,1, filename,2, …, filename,n, where n is the length of the 4th dimension. For a single-volume file, filename is therefore equivalent to filename,1.

nifti object

A nifti object represents a NIfTI image, either in dual- (img+hdr) or single-file (nii) format. Important fields are:

hdr	struct containing fields of the NIfTI header (hidden)
dat	image data represented by a file_array object (= D), initialized from dim, vox_offset, scl_slope, scl_inter, and detected endianess.
raw	raw image data (not linearly transformed), i.e. a copy of dat with empty D.scl_slope and D.scl_inter (hidden, read-only)
mat0	4 × 4 affine transformation matrix derived from qform
mat	4 × 4 affine transformation matrix derived from sform if present, otherwise falls back on qform (= mat0)
descrip	description of the image

This description references the NIfTI header fields stored in hdr.

Further fields that are often unset or unreliable:

extras	additional data read from companion Matlab mat-file
mat0_intent	intention of mat0 ('UNKNOWN' | 'Scanner' | 'Aligned' | 'Talairach' | 'MNI152')
mat_intent	intention of mat
intent	interpretation of image data
diminfo	MR encoding of different dimensions
timing	timing information
cal	display calibration
aux_file	name of an auxiliary file

The nifti object contained in a volume struct represents the whole NIfTI image, not just a specific volume, and its fields dat and raw the whole (up to 7D) image data.

file_array object

A file_array object D provides memory-mapped access to (a part of) a data file. It can be syntactically treated like a regular multidimensional array. Reference to array elements, D(i, j, k), is transparently translated into file read operations, and assignment to array elements, D(i, j, k) = y, translated into file write operations. Its fields are:

fname	name of the data file
dim	dimensions
dtype	data type including endianness (string), based on properties dtype and be
offset	offset into file to data start (bytes)
scl_slope	slope
scl_inter	intercept

The byte sequence in the file, starting from offset, is interpreted as a sequence of prod(dim) numbers in the format given by dtype. These raw numbers are linearly transformed by multiplying with scl_slope and adding scl_inter. The resulting sequence of numbers is treated like the linearly indexed contents of the array D. If both scl_slope and scl_inter are empty, the linear transform is omitted, i.e. D provides access to the raw numbers.

Supported data types and corresponding Matlab types are

'INT8'	int8
'UINT8'	uint8
'INT16'	int16
'UINT16'	uint16
'INT32'	int32
'UINT32'	uint32
'FLOAT32'	single
'FLOAT64'	double

The dtype field consists of the data type string, plus the suffix '-LE' or '-BE' to indicate endianness.

SPM functions for volume structs, nifti, and file_array objects

V = spm_vol(P)

creates an array of volume structs V from image filename(s) P and reads the header. Each filename may be appended with a volume reference of the form ‘,#’ or ‘,#,#’, see above. private.dat can now be read from (and written to).

Y = spm_read_vols(V)

reads image data Y from an array of volume structs V. Internally it calls spm_slice_vol, which uses the volume pinfo for linear transformation (and ignores D.scl_slope and D.scl_inter).

V = spm_create_vol(V)

writes the image header (beginning of nii or separate hdr file) according to the fields of the volume struct V. Fields fname, dim, and mat need to be set. If fields dt, pinfo, n, descrip, and private are not present, they are added with default values. pinfo(3) is patched with the data offset into the file (usually 352 for nii and 0 for img). Multiple volumes (in separate files, or a multi-volume NIfTI image if n is set correctly) can be created at once if V is an array (but this fails if dt and pinfo are not set in all volumes). private.dat can now be written to (and read from as far as it has been written), which creates the img file if applicable.

V = spm_write_vol(V, Y)

writes 3D data Y to an image represented by the single volume struct V. pinfo(1) and pinfo(2) are set heuristically if pinfo is not present or either one is Inf, spm_create_vol is called, and writing is done via spm_write_plane (which uses private.dat).

spm_check_orientations(V)

determines whether all volumes in the struct array V have the same values of dim and mat; if not, an error is thrown. Alternatively, the result of the check can be queried by giving an output argument (false in case of mismatch).

Within SPM12, these functions are mainly used through the wrappers spm_data_hdr_read, spm_data_read, spm_data_hdr_write, and spm_data_write. Note that spm_slice_vol reads data via resampling based on an affine transform, which is unnecessary and makes spm_read_vols slower than direct data access via private.dat.

N = nifti(P)

creates an array of nifti objects N from image filenames P and reads the header. N.dat can now be read from (and written to).

N = nifti

creates a nifti object that is minimally initialized and not associated with an image file. After setting at least the field dat to a file_array object (thereby opening or creating a file), create can be used to write the header information.

create(N)

writes NIfTI header information if N.dat references an image file, either to the beginning of the file (nii) or to a separate header file (img → hdr).

Reading and writing of NIfTI files (representing up to 7D arrays) can be completely achieved using nifti objects. The purpose of SPM’s volume struct arrays is to reference sequences of single volumes (3D arrays) within and across image files.

D = file_array(fname, dim, dtype, offset, scl_slope, scl_inter)

creates a

file_array object referencing an existing or not existing file. If it does not exist, writing to D creates the file as far as necessary. D can be read from as far as the file is written (within the array index bounds). If the file_array object is added to a nifti object, the offset is automatically set correctly.

SPM’s file_array objects are tailored for accessing image data in NIfTI files, but they can just as well be used for memory-mapped access to arbitary data files, interpreting their contents as linearly indexed elements of a multidimensional array.

First-level analysis

Module ‘fMRI model specification’

call hierarchy:

spm_run_fmri_spec: Setting up the general linear model for fMRI time-series
    spm_fmri_spm_ui: Setting up the general linear model for fMRI time-series
        spm_fMRI_design: Assembles a design for fMRI studies
            1st save of SPM.mat
        modification of xY.VY(:).pinfo for global scaling
        2nd save of SPM.mat
    3rd save of SPM.mat (only adds explicit mask)

Fields of SPM.mat

xY	data
P	volume names (filename,#) across sessions; two-dimensional char array
VY(:)	volume struct for each scan
RT	repeat time
nscan(:)	number of scans for each session
xBF	basis functions
T	microtime resolution (time points per scan)
T0	microtime onset (reference for onset times)
UNITS	units in which onsets are specified ('scans'|'secs')
Volterra	order of convolution, 1: standard, 2: plus interactions
dt	time bin length (seconds)
name	description of basis functions
length	window length (seconds)
order	order of basis set
bf	set of basis functions (orthogonal but not normalized)
Sess(:)	sessions
U(:)	conditions
dt	same as xBF.dt
ons	onsets (in xBF.UNITS)
dur	durations (in xBF.UNITS)
P(:)	parametric modulations
name	name of modulator ('time'|...)
h	order of polynomial expansion (number of powers)
i	related causes (column indices into Sess.U.u)
u	set of causes (condition and modulated), as stick or boxcar functions of microtime
name	names of causes
orth	are modulated regressors orthogonalized?
pst	peristimulus times (seconds)
C	user-specified regressors
C	set of regressors
name	names of regressors
row	scans of session (rows of xX.X)
col	regressors for session (columns of xX.X)
Fc(:)	F-contrasts across all regressors for each condition or interaction (see xBF.Volterra)
name	name of the contrast
i	involved regressors (indices into Sess.col)
p	for each regressor, involved cause (column indices into Sess.U.u)
xX	design
X	design matrix (set of regressors; raw)
name{:}	names of regressors
iH	indices of indicator variables (unused?)
iC	indices of covariates (conditions / modulated / interactions & user-specified)
iB	indices of block effects (constant regressors for each session)
iG	indices of nuisance variables (unused?)
K(:)	high-pass filter for each session
HParam	cut-off period (seconds)
row	same as corresponding Sess.row
RT	same as xY.RT
X0	set of low-frequency components to be removed (orthonormal cosines)
xGX	global variate
iGXcalc	global scaling option, 'none': session specific, 'scaling': scan specific
sGXcalc	calculation method, always 'mean voxel value'
rg	estimated global variate for each scan
sGMsca	global scaling method, erroneously(?) always 'session specific'
GM	target value for global scaling, always 100
gSF	global scaling factor for each scan, constant within session if session specific
xM	masking
gMT	masking threshold parameter, fraction of global variate
TH	masking threshold for each scan, scaled version of gMT
VM(:)	volume struct for each explicit mask image
xs.Masking	description, either 'analysis threshold' or 'analysis threshold+explicit mask'
T	masking index (?) for each scan, always 1s
I	use implicit mask?
xVi	non-sphericity
form	form of non-sphericity ('i.i.d.' | 'FAST' | 'AR(0.2)' | 'AR(#)' user-specified parameter)
Vi{:}	covariance components (correlation matrices), for each session one (i.i.d), two (AR), or many
xsDes	design description
Interscan_interval	from xY.RT
Number_of_sessions	from length of nscan
Basis_functions	same as xBF.name
Trials_per_session	from lengths of Sess(:).U
High_pass_Filter	from minimum of xX.K(:).HParam
Global_normalisation	same as xGX.iGXcalc
Global_calculation	same as xGX.sGXcalc
Grand_mean_scaling	same as xGX.sGMsca
factor(:)	design factors, slowest first
name	name of factor
levels	number of levels
SPMid	SPM version and last function to write SPM.mat

Module ‘Model estimation’, method ‘Classical’

depends on ‘fMRI model specification’

call hierarchy:

spm_run_fmri_est: Estimate parameters of a specified model
    spm_spm: [Re]ML Estimation of a General Linear Model
        spm_est_non_sphericity: Non-sphericity estimation using ReML
        write beta_*, mask, and ResMS images
        also ResI_* images for smoothness estimation
        spm_est_smoothness: Estimation of smoothness based on [residual] images
            write RPV image
        delete ResI_* images
        1st save of SPM.mat
    spm_conman: Contrast manager: GUI for defining valid contrasts
    for each contrast
        spm_contrasts: Compute and store contrast parameters and inference SPM{.}
            write ess & spmF or con & spmT image
            2nd and following save of SPM.mat
    write Res_* images

Fields of SPM.mat

xVi	non-sphericity
h	hyperparameters (weights of covariance components xVi.Vi{:})
V	estimated serial non-sphericity (linear combination of covariance components, scaled to mean(diag(V)) = 1)
Cy	spatial covariance between scans
UFp	from defaults.stats.fmri.ufp
xX	design
W	whitening matrix (approx. sqrtm(inv(xVi.V)))
xKXs	space structure (see spm_sp) of K W X
X	filtered and whitened design matrix
ds	its singular values
u	its left singular vectors
v	its right singular vectors
tol	numerical tolerance (depends on max(abs(ds)))
rk	its rank (based on tol)
oP	orthogonal projector on X
oPp	orthogonal projector on X'
ups	space in which this one is embeded
sus	subspace
pKX	pseudoinverse of xX.xKXs.X
V	filtered and whitened nonsphericity K W xVi.V W' K'
trRV	trace(R V) after filtering and whitening, based on R from xX.xKXs.X and xX.V
trRVRV	trace(R V R V) after filtering and whitening
erdf	effective residual degrees of freedom after filtering and whitening, trace(R V) ^ 2 / trace(R V R V)
Bcov	nonsphericity of parameter estimates, gives estimation covariance if multiplied with voxel-specific ResMS
nKX	xX.xKXs.X scaled for display
Vbeta(:)	volume struct of parameter estimates image 'beta_*', for each regressor; data generated by xX.pKX * spm_filter(xX.K, W * Y)
VResMS	volume struct of residual mean of squares image 'ResMS' (unbiased estimator of error variance)2
VM	volume struct of analysis mask image, i.e. combination of explicit mask(s) xM.VM, threshold mask xM.TH, and implicit (zero/NaN) mask
xVol	search volume
XYZ	coordinates of in-mask voxels (in mm)
S	number of in-mask voxels
M	from xY.VY(1).mat
iM	inv(M)
DIM	from xY.VY(1).dim
VRpv	volume struct of resels per voxel image 'RPV'
FWHM	full width at half maximum of average Gaussian point spread function, three-dimensional; based on resels per voxel
R	“resel counts”, see spm_resels_vol; based on mask and FWHM
xCon(:)	contrast structure, for each contrast; based on factor
name	name of contrast
STAT	'F' or 'T'
c	contrast column vector or matrix
X0	reduced design matrix
ukX0	in the basis of xX.X
iX0	how contrast was specified
X1o	remaining design space
ukX1o	in the basis of xX.X
eidf	effective interest (numerator) degrees of freedom
Vcon	volume struct of 'con' or 'ess' image
Vspm	volume struct of 'spmT' or 'spmF' image
SPMid	SPM version and last function to write SPM.mat
swd	directory of SPM.mat

Module ‘Model estimation’, method ‘Bayesian 1st-level’

information in this section is preliminary!

depends on ‘fMRI model specification’

Fields of SPM.mat

PPM	parameters and some results of Bayesian estimation
space_type
block_type
priors
AR_P	order of autoregressive model of nonsphericity
VB
window
Sess(:)	for each session
VHp	volume struct of standard deviation error image
VAR(:)	volume struct of AR parameter image, for each parameter
block
update_F
AN_slices
maxits
tol
compute_det_D
Gamma
xCon	? contrast structure
VCbeta(:)	volume struct of posterior mean of beta image, for each regressor
VM	volume struct of analysis mask image
VPsd	volume struct of posterior standard deviation of beta image, for each regressor
xCon	? contrast structure (again)
xVol	search volume
XYZ	coordinates of in-mask voxels (in mm)
S	number of in-mask voxels
M	from xY.VY(1).mat
iM	inv(M)
DIM	from xY.VY(1).dim
VRpv	? set to []
FWHM	? “full width at half maximum”, but scalar
R	? “resel counts”, but scalar
labels	?
xX	design
W	“whitening matrix”, but set to the identity matrix
erdf	“effective residual degrees of freedom”, but set to the number of scans
xKXs	? space structure (see spm_sp) of K W X, but different from ‘Classical’, possibly no whitening
nKX	xX.xKXs.X scaled for display
swd	directory of SPM.mat

Other resources

Nikki Sullivan’s SPM data structures

NiBabel Image use-cases in SPM

SPM Wikibook Programming intro

Footnotes

1. The nifti and file_array objects are implemented using Matlab’s pre-2008a object-oriented syntax.]
   In the following, we refer to the contained nifti object as N and the contained file_array object as D.

2. Residuals res are obtained by applying the residual-forming matrix of xX.xKXs to spm_filter(xX.K, W * Y), then ResMS = sum(res .^ 2) / xX.trRV; later regularized by ResMS = ResMS + 1e-3 * max(ResMS(isfinite(ResMS)))