SPM for programmers

Gerrit from the blog Sustainable Research posted a nice article about the screen-command for the linux console that allows you to run scripts in the background without an open console window. This can be also useful, when you want to run SPM-scripts without a user interface and remote on a server. You just start the script within a screen, exit the screen and close the remote connection. And with Putty for Android or iPhone you can supervise your scripts from everywhere :).

I suppose that this is not the most efficient way to do this, but it works. With the following commands you can easily write a script that displays any kind of nifti-image (raw files, mask, beta-maps) across all subjects and saves it as a JPG image. This can be handy for quality assurance:

spm_image('init','niftifile.nii') % display a nifti-file
spm_orthviews('Reposition', [22 0 -21]) % jump to a certain coordinate
fs = get(0,’Children’);
res = getframe(fs(1));
imwrite(res.cdata, ‘niftifile.jpg’); % save the window as JPG

This should work for SPM5 and SPM8. Any suggestions how to improve the “save as jpg” procedure?

I guess I am not the only one who uses matlab scripts to analyze large amounts of fMRI data. To get rid of the popups of SPM during the processing, you can start Matlab either with

matlab -nodesktop


to solely start matlab on the command line, or you start it with

matlab -noFigureWindows


which starts Matlab as normal but prevents the depiction of any additional windows.

This post was inspired by a question that I received by email. Thanks!

The high-pass filter can have a great effect on the subsequent parameter estimation. To explore, how a certain setting (default is 128) influences the time course of the regressor you can use the spm_filter function. First, we build a HRF-regressor:

reg = [zeros(10,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1)];


bf = spm_get_bf %TR is 2.2 in our case
U.u = reg;
U.name = {'reg'};
% create the regressor convolved with the HRF
convreg = spm_Volterra(U, bf.bf);


In your SPM.mat you can find this kind of regressors under SPM.xX.X.

spm_filter is a function that either creates the low frequencies to be removed or removes those frequencies from the time course. The input variable K defines all necessary parameters:

K.RT = 2.2;
K.row = 1:130; % 130 corresponds to the length of convreg
K.HParam = 128; % cut-off period in seconds


Here, I defined K.row as it also happens internally. But I guess, ones(1,130) should also work, as it seems that K.row just serves as a mask. Now, we compute the low frequencies that should be removed from convreg.

nK = spm_filter(K);


nK.X0 contains the frequencies that can be removed. A second call of spm_filter removes these frequencies.

Y = spm_filter(nK, convreg);


In SPM.mat these regressors are stored in SPM.xX.nKX and are depicted in the design matrix. As you can see, the high-pass filter cut-off period can greatly affect the regressors that are used for the parameter estimation:

K.HParam = 256;
figure; plot(spm_filter(spm_filter(K), convreg));


K.HParam = 64;
figure; plot(spm_filter(spm_filter(K), convreg));


If you experiment with DCM, it might be interesting to know that you can build synthetic input data for DCM with the spm_dcm_create-function. Just type spm_dcm_create on the console, press enter and follow the steps. Other interesting units for DCM are: spm_dcm_estimate (the estimation procedure) and spm_dcm_compare (SPM8) for comparing two or more models.

With the help of two SPM-functions you can manually convolve time series in Matlab. This might help to understand how the convolution with the HRF (or any other function) affects the time series. Here is an example:

% create a time series
reg = [zeros(10,1); ones(10,1); zeros(10,1); ones(10,1); zeros(10,1)];
plot(reg)


% create a basis function, in this case the hrf
bf = spm_get_bf
% with a tr of 2.2 and HRF as basis function you receive
dt: 2.2000
name: 'hrf'
length: 33
order: 1
bf: [15x1 double]
% these are necessary definitions for the spm_Volterra-function
U.u = reg;
U.name = {’reg’}; % U needs a name, but the string has no meaning
% convolve reg with the hrf
convreg = spm_Volterra(U, bf.bf);
plot(convreg)


To convert matrix coordinates for a given Nifti/Analyze-image into the image coordinates and vice versa, you need to load the transformation matrix of the image into the workspace:

vol = spm_vol('image.img')
vol.mat % displays the transformation matrix


% trasform from matrix/voxel coordinates to image coordinates (the fourth element has to be 1)
vol.mat * [x y z 1]‘
% transform image coordinates (e.g. mni space) to matrix coordinates
inv(vol.mat) * [x y z 1]‘


The “extent threshold” function deletes all clusters in a given contrast image that consist of less than k voxels. The corresponding matlab code can be found in spm_getSPM.m (line 674f). In case you would like to use this code manually for other imaging data, here is a matlab-function called extentThreshold that allows you to apply the extent threshold on a 3d volume of your choice.

The function can be used in the following way: if data is a 3d-matrix (of let’s say voxel values) and you want to display clusters with voxels that are greater than zero and with a cluster size of at least 20, write:

indices = find(data>0);
newdata = extentThreshold(data, indices, 20);


newdata includes only the elements of data that are greater than zero and that belong to a cluster with at least 20 other elements exceeding this threshold.

SPM has a nice feature to plot the frequency domain of a signal. You can use it, for example, when you click on “Review” and - after you selected an SPM.mat file - on “Design” - “Explore” - “Session 1″ - Regressor. The corresponding function can be found in spm_fMRI_design_show.m in line 90ff.

To call this plot directly from the Matlab console e.g. in order to visualize the frequency domain of extracted VOIs, I copied the code into a separate m-file function, which looks like this:

function fft_gui(signal, rt, HPF)
gX = abs(fft(signal)).^2;
gX = gX*diag(1./sum(gX));
q = size(gX,1);
Hz = [0:(q - 1)]/(q*rt);
q = 2:fix(q/2);
plot(Hz(q),gX(q,:))
patch([0 1 1 0]/HPF,[0 0 1 1]*max(max(gX)),[1 1 1]*.9);
xlabel(’Frequency (Hz)’)
ylabel(’relative spectral density’)
title(['Frequency domain',sprintf('\n'), ' {\bf',num2str(HPF),'}', ...
' second High-pass filter'],’Interpreter’,'Tex’);
grid on
axis tight

To display the frequency domain of a signal, just call this function with the appropriate RT and high-pass filter value, e.g.:

signal = (sin((1:100)*0.8)+sin((1:100)*1.4) + sin((1:100)*2))';
fft_gui(signal, 3, 128)

Here are some useful tips for displaying results in SPM:

Link