John Ashburner is an author of SPM and frequent contributor to the SPM email list. Many of his answers consist of short snippets of code that are incredibly useful for SPM5 analyses. Below are the most some of the best (the 'gems') in an accessible form. (The inspiration for this page is "Graphics Gems", a series of books containing short bits of useful code for computer graphics.)
We don't yet know about any undiscovered bugs, so we aren't sure if they exist or not.
-JA, Mar 21, 2003.
I haven't tested them all but I'll make bug fixes as I find them. I hope you find this useful.
Gem 1: Introduction to SPM5 scripting
Date: Tue, 13 Dec 2005 17:03:13 +0100 From: John Ashburner <john@FIL.ION.UCL.AC.UK> To: SPM@JISCMAIL.AC.UK Subject: Re: [SPM] Where can we find some development materials for SPM ? > As you know,we usually need to modify the code of SPM to fit our > problem.but we can not find some relevant development tutorials. Would you > please tell me how to learn the framework of SPM step by step ? > In addition, I want to know where I can find the details of the SPM > structure. It may be easiest to learn by example. If you want to develop a new user-interface for SPM5, then you would create a file called spm_config_*.m, similar to the other spm_config.m files (if you strip out the documentation parts, you will see that these are actually quite small). Your spm_config* file can then be added to the toolbox subdirectory and accessed through the TOOLS pulldown. The help button allows you to navigate through the documentation of each of the Matlab functions, which you may find useful. For reading and writing images, you would use these functions. spm_vol spm_slice_vol spm_sample_vol spm_create_vol spm_write_plane spm_write_vol spm_get_space There is Matlab help on all these functions. Alternatively, you could use the NIFTI routines directly. There is no documentation on this, but here are a few examples of how you can use them: ====================================== % Example of creating a simulated .nii file. dat = file_array; dat.fname = 'junk.nii'; dat.dim = [64 64 32]; dat.dtype = 'FLOAT64-BE'; dat.offset = ceil(348/8)*8; % alternatively: % dat = file_array( 'junk.nii',dim,dtype,off,scale,inter) disp(dat) % Create an empty NIFTI structure N = nifti; fieldnames(N) % Dump fieldnames % Creating all the NIFTI header stuff N.dat = dat; N.mat = [2 0 0 -110 ; 0 2 0 -110; 0 0 -2 92; 0 0 0 1]; N.mat_intent = 'xxx' N.mat_intent = 'Scanner'; N.mat0 = N.mat; N.mat0_intent = 'Aligned'; N.diminfo.slice = 3; N.diminfo.phase = 2; N.diminfo.frequency = 2; N.diminfo.slice_time.code='xxx'; N.diminfo.slice_time.code = 'sequential_increasing'; N.diminfo.slice_time.start = 1; N.diminfo.slice_time.end = 32; N.diminfo.slice_time.duration = 3/32; N.intent.code='xxx' ; % dump possibilities N.intent.code='FTEST'; % or N.intent.code=4; N.intent.param = [4 8]; N.timing.toffset = 28800; N.timing.tspace=3; N.descrip = 'This is a NIFTI-1 file'; N.aux_file='aux-file-name.txt'; N.cal = [0 1]; create(N); % Writes hdr info dat(:,:,:)=0; [i,j,k] = ndgrid(1:64,1:64,1:32); dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 30^2))=1; dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 15^2))=2; % displaying a slice imagesc(dat(:,:,12));colorbar % get a handle to 'junk.nii'; M=nifti('junk.nii'); imagesc(M.dat(:,:,12)); ====================================== Best regards, -John
Gem 2: Reslicing images
Gem 3: Resizing images
A generalized version of John's reorient script (see Gem2) by Ged Ridgeway, which allows specification of arbitrary voxel dimensions: resize_img.m (this was current as of April 2007; see Ged's script directory for updates.)
Gem 4: Switching between SPM versions
This function will allow you to switch between different SPM versions. WARNING! As SPM depends on various global (and sometimes, local workspace) variables, this function clears all variables as part of the switch.
The function will need to be put in a directory in your Matlab path that does not contain SPM. SPMsel.m
Gem 5: Unnormalizing a point
This script of Johns will find the corresponding co-ordinate in the un-normalized image: get_orig_coord2.m (same script as for SPM2).
Gem 6: Corrected cluster size threshold
This is SPM5 version of SPM2 Gem 13.
This is a script that will tell you the corrected cluster size threshold for given cluster-defining threshold: CorrClusTh.m
The usage is pretty self explanatory:
Find the corrected cluster size threshold for a given alpha function [k,Pc] =CorrClusTh(SPM,u,alpha,guess) SPM - SPM data structure u - Cluster defining threshold If less than zero, u is taken to be uncorrected P-value alpha - FWE-corrected level (defaults to 0.05) guess - Set to NaN to use a Newton-Rhapson search (default) Or provide a explicit list (e.g. 1:1000) of cluster sizes to search over. If guess is a (non-NaN) scalar nothing happens, except the the corrected P-value of guess is printed. Finds the corrected cluster size (spatial extent) threshold for a given cluster defining threshold u and FWE-corrected level alpha.
To find the 0.05 (default alpha) corrected cluster size threshold for a 0.01 cluster-defining threshold:
>> load SPM >> CorrClusTh(SPM,0.01) For a cluster-defining threshold of 2.4671 the level 0.05 corrected cluster size threshold is 7860 and has size (corrected P-value) 0.0499926Notice that, due to the discreteness of cluster sizes you cannot get an exact 0.05 threshold.
The function uses an automated search which may sometimes fail. If you specify a 4th argument you can manually specify the cluster sizes to search over:
>> CorrClusTh(SPM,0.01,0.05,6000:7000) WARNING: Within the range of cluster sizes searched (6000...7000) a corrected P-value <= alpha was not found (smallest P: 0.0819589) Try increasing the range or an automatic search.Ooops... bad range.
Lastly, you can use it as a look up for a specific cluster size threshold. For example, how much over the 0.05 level would a cluster size of 7859 be?
>> CorrClusTh(SPM,0.01,0.05,7859) For a cluster-defining threshold of 2.4671 a cluster size threshold of 7859 has corrected P-value 0.050021Just a pinch!