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
 
          

...to the top

Gem 2: Reslicing images

To reslice an image at 1mm cubic voxels, axial orientation, use this reorient.m script (from an email from John dated Mon, 5 Jun 2006 13:02:05 +0100; see also Gem3 and SPM99 Gem7).

...to the top

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.)

...to the top

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

...to the top

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).

...to the top

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.0499926
  

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.

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.050021

...to the top