John's Gems
Please note that since 2010 I no longer update John's Gems, and have instead integrated them into my blog Neuroimaging Tips & Tools. Please use the JohnGems tag in the blog to find all the Gems.
Essential code snippets by John Ashburner
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 SPM analyses or general data manipulation. 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.
See also the SPM2 gems page and the SPM5 gems page.
Tom Nichols
Gem 1: Adding blobs to an imageSay you have a statistic image that you created outside of SPM and you want to overlay it on a structural image. The bare bones solution is as follows. Create a thresholded version of the image, where all subthreshold voxels are set to NaN (See how to set NaNs below). Then... Subject: Re: Your Message Sent on Fri, 19 Jan 2001 23:28:45 0800 From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Mon, 22 Jan 2001 11:40:24 +0000 To: SPM@JISCMAIL.AC.UK I think it may take quite a bit of programming to get most of the routines working with your data. One thing to get you started would be to try an undocumented feature of spm_orthviews.m: spm_orthviews('image',spm_get(1,'*.img','Select background image')); spm_orthviews('addimage',1,spm_get(1,'*.img','select blobs image')); Best regards, JohnHere's a different approach, which I understand will show the blobs in a monotone color (instead of hot metal). Also, there are more bells and whistles (it's less of a hack). Subject: Re: overlaying resulting image From: John Ashburner <john@fil.ion.ucl.ac.uk> Date: Wed, 4 Oct 2000 10:37:01 +0100 (BST) To: spm@mailbase.ac.uk, duann@salk.edu  Would you please show me how to overlay a statistical map onto the brain  template SPM uses. Let's say I have an analyzeformated statistical  map obtained from other software. It was already normalized to the same  coordinates as the template images SPM has. How can I overlay this image  onto the SPM template just like the results shown in SPM convention. The following is supposed to work. It uses a hidden undocumented feature in SPM99, so it may contain bugs.... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P1 = spm_get(1,'*.img','Specify background image'); P2 = spm_get(1,'*.img','Specify blobs image'); % Clear graphics window.. spm_clf % Display background image.. h = spm_orthviews('Image', P1,[0.05 0.05 0.9 0.9]); % Display blobs in red. Use [0 1 0] for green, [0 0 1] for blue % [0.6 0 0.8] for purple etc.. spm_orthviews('AddColouredImage', h, P2,[1 0 0]); % Update the display.. spm_orthviews('Redraw'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Best regards, John Gem 2: Zeroing NaN valuesIn addition to setting NaN values to zero, this gem is also useful as a basic skeleton for reading each slice, doing something to it, and then writing it out. Subject: Re: NaN values in beta values From: John Ashburner <john@fil.ion.ucl.ac.uk> Date: Mon, 23 Oct 2000 17:27:24 +0100 (BST) To: spm@mailbase.ac.uk, steffejr@umdnj.edu I can't think of any nice way of doing this with ImCalc, so I'm afraid I'll have to show you a more labour intensive method that should be modified slightly (as the filenames need changing) before pasting into Matlab. VI = spm_vol('original_beta.img'); VO = VI; VO.fname = 'patched_beta.img'; VO = spm_create_image(VO); for i=1:VI.dim(3), img = spm_slice_vol(VI,spm_matrix([0 0 i]),VI.dim(1:2),0); tmp = find(isnan(img)); img(tmp) = 0; VO = spm_write_plane(VO,img,i); end; Gem 3: NaNing zero valuesSubject: Re: normalization of contrast images From: John Ashburner[...]  And, while I am at it, does anybody know how to convert a mask that  has zeros outside the brain to a mask that has NaN outside the brain?  I tried ImCalc with something like 'i1(find(i1==0))=NaN', but it  doesn't like it. By default, ImCalc outputs the data as 16 bit integer with scalefactors. There is no NaN representation for this, so the data would need to be written as floating point or double precision floating point. I think you can do this by typing something like: P = spm_get(1,'*.img'); Q = 'output.img'; f = 'i1.*(i1./i1)'; flags = {0,0,spm_type('float'),1}; spm_imcalc_ui(P,Q,f,flags); Best regards, JohnNote you can NaNout voxels below a threshold, say, 3:
Also note that this is a general way for making spm_imcalc write out images with float or double precision; concisely
Using [],[] instead of 0,0 ensures that default values will be used for those two flags. Gem 4: Histogram of an imageSubject: Re: tdistribution display (histogram) with spm99b From: john@fil.ion.ucl.ac.uk (John Ashburner) Date: Fri, 30 Jul 1999 14:06:20 +0100 To: spm@mailbase.ac.uk, jovicich@humc.edu The attached program should produce the histograms you are after. To call it, type: V = spm_vol(spm_get(1,'*.img','Select image...')); [n, x] = histvol(V, 100); figure; bar(x,n);[...] Regards, JohnThe attached function is here histvol.m and below
function [n, x]=histvol(V, nbins) % Create Histogram of an image volume % FORMAT [n, x]=histvol(V, nbins) % V  mapped image volume (see spm_vol) % nbins  number of bins to use. % n  number of counts in each bin % x  position of bin centres %_______________________________________________________________________ % @(#)JohnsGems.html 1.42 05/02/02 if nargin==1, nbins = 256; end; % determine range... mx = Inf; mn = Inf; for p=1:V.dim(3), img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1); msk = find(isfinite(img)); mx = max([max(img(msk)) mx]); mn = min([min(img(msk)) mn]); end; % compute histograms... x = [mn:(mxmn+1)/nbins:mx]; n = zeros(size(x)); for p=1:V.dim(3), img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1); msk = find(isfinite(img)); n = n+hist(img(msk),x); end; return; Gem 5: Mesh Plots of tmapsSubject: Re: Mesh Plots of tmaps From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Fri, 2 Feb 2001 11:14:20 +0000 To: SPM@JISCMAIL.AC.UK  I am wishing to create Mesh plots in Matlab of the tmap produced  from my SPM analysis. Being a somewhat Matlab virgin, I was wondering  whether  someone could possibly point me in the direction of exactly how to do this.   Which mat file contains all tvalues in the statistic image? The t values are stored in spmT_????.img files, where ???? refers to the contrast number. You can read the values from these images into Matlab something like: pl = 30; % plane 30 fname = spm_get(1,'*.img','Name of t image'); V = spm_vol(fname); M = spm_matrix([0 0 pl]); img = spm_slice_vol(V,M,V.dim(1:2),1); Displaying the values can be dome something like: surf(img); There are loads of commands for 3D visualisation in Matlab 5.x. You can find out what these are by typing: help graph3d All the best, John Gem 6: Reading raw dataSubject: Re: raw data From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Thu, 10 May 2001 14:33:19 +0100 To: SPM@JISCMAIL.AC.UK > I have a very simple question. Where does one find and how does one go > about extracting the raw data associated with a particular voxel? I would write a very short script to do this. Try modifying and copypasting the following.... * * * * * * * * * * * * * * * * * * * * * * * * * * * *  V=spm_vol(spm_get(Inf,'*.img')); x = 32; y = 32; z = 32; dat = zeros(length(V),1); for i=1:length(dat), dat(i) = spm_sample_vol(V(i),x,y,z,0); end; * * * * * * * * * * * * * * * * * * * * * * * * * * * *  The vector dat will contain the raw data from the voxel at 32,32,32. Best regards, John Gem 7: Reorienting imagesSubject: Re: ACPC positions From: John Ashburner <john@fil.ion.ucl.ac.uk> Date: Fri, 27 Oct 2000 15:00:05 +0100 (BST) To: spm@mailbase.ac.uk, spm@fil.ion.ucl.ac.uk[...] The best that I can suggest is you try manually reorienting your images via the <Display> button. Try different rotations and translations until the image is displayed how you want it. The attached Matlab function can then be used for reslicing the image(s) in the transverse orientation, with 1mm isotropic resolution. Best regards, JohnAttached file: reorient.m Error in c matrix fixed, June 6, 2001 TEN Subsequently, John offered modifications to use the native voxel size, which have been incorporated... Subject: Re: reorient.m From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Thu, 14 Dec 2000 17:04:38 +0000 To: SPM@JISCMAIL.AC.UK  I would like to reorient some images in two ways,  using the usefull function reorient() of John Ashburner.    first is to keep the voxel size intact,  but still reorienting in the transverse plane Without testing the code, and without too much thought, I think the modification to do this is involves something like changing from: mat = spm_matrix([mn1]); dim = (mat\[mx 1]')'; to something like: vox = spm_imatrix(V.mat); vox = vox(7:9); mat = spm_matrix([0 0 0 0 0 0 vox])*spm_matrix([mn1]); dim = (mat\[mx 1]')';    second is to reorient in the coronal plane,  with 1x1x1 mm resol. I think this is involves something like changing from: mat = spm_matrix([mn1]); to something like: mat = spm_matrix([0 0 0 pi/2])*spm_matrix([mn1]); or maybe: mat = spm_matrix([0 0 0 pi/2])*spm_matrix([mn1]);   In the two cases, the final image should have no .mat file  and be resliced using sinc ... To reslice using sinc interpolation, you change from: img = spm_slice_vol(V,M,dim(1:2),1); to something like: img = spm_slice_vol(V,M,dim(1:2),6); I hopehese suggestions work. Best regards, JohnIf you want to set arbitrary voxel size, just set vox as desired in the fix above, but then be sure to force dim to be an integer. For example, I wanted to increase the resolution of my images by a factor of three, so I did vox = spm_imatrix(V.mat); vox = vox(7:9)/3; mat = spm_matrix([0 0 0 0 0 0 vox])*spm_matrix([mn1]); dim = ceil(mat\[mx 1]')'); Gem 8: Reslicing imagesThis is for soley reslicing images, but it has a nice description of how to create a basic transformation matrix. Subject: Re: how to write out a resampled volume with new user specified voxel size From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Thu, 4 Jan 2001 10:57:52 +0000 To: SPM@JISCMAIL.AC.UK  I would like to batch resample some volume images to isotropic  voxels. Is there an SPM99 command line instruction (or menu option) that  would permit me to write out a resampled volume with a user specified voxel  size dimensions? I know this core function is implicit in several  modules. I considered using Coregister or Spatial Normalization modules  but these require a target or template volume of the desired resolution  which doesn't always exist in my application.  Many thanks for any suggestions. You could try the attached function that I've just scribbled together. It hasn't been tested, so I hope it works. Best regards, JohnThe attached function is here reslice.m and below
function reslice(PI,PO,dim,mat,hld) % FORMAT reslice(PI,PO,dim,mat,hld) % PI  input filename % PO  output filename % dim  1x3 matrix of image dimensions % mat  4x4 affine transformation matrix mapping % from vox to mm (for output image). % To define M from vox and origin, then % off = vox.*origin; % M = [vox(1) 0 0 off(1) % 0 vox(2) 0 off(2) % 0 0 vox(3) off(3) % 0 0 0 1]; % % hld  interpolation method. %___________________________________________________________________________ % @(#)JohnsGems.html 1.42 John Ashburner 05/02/02 VI = spm_vol(PI); VO = VI; VO.fname = deblank(PO); VO.mat = mat; VO.dim(1:3) = dim; VO = spm_create_image(VO); end; for x3 = 1:VO.dim(3), M = inv(spm_matrix([0 0 x3 0 0 0 1 1 1])*inv(VO.mat)*VI.mat); v = spm_slice_vol(VI,M,VO.dim(1:2),hld); VO = spm_write_plane(VO,v,x3); end; Gem 9: Rotations and translations from .mat filesSubject: Re: rot+trans From: John Ashburner <john@fil.ion.ucl.ac.uk> Date: Fri, 6 Oct 2000 16:15:23 +0100 (BST) To: spm@mailbase.ac.uk, pauna@creatis.insalyon.fr  1 how can I find the transformations  (rot+translations) from the transformation matrix? The relative translations and rotations between a pair of images is encoded in the .mat files of the images. So for images F.img and G.img, you would find the transformation by: % The mat file contents are loaded by: MF = spm_get_space('F.img'); MG = spm_get_space('G.img'); % A rigid body mapping is derived by: MR = MF/MG; % or MR = MG/MF; % From this, a set of rotations and translations % can be obtained via: params = spm_imatrix(MR) % See spm_matrix for an explanation of what these % parameters mean Gem 10: Surface Renderings: Rolling your own (w/ blobs!)Subject: Re: rendering SPM96 From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Tue, 20 Feb 2001 16:22:05 +0000 To: SPM@JISCMAIL.AC.UK SPM96 shows voxels that are considered to be between 5mm in front of and 20mm behind the surface that is displayed. This is the same as for the "old style" rendering of SPM99. To change this depth, try changing line 131 of spm_render.m from: msk = find(xyz(3,:) < (z1+20) & xyz(3,:) > (z15)); to something like: msk = find(xyz(3,:) < (z1+5) & xyz(3,:) > (z15)); I would rather see the rendering done on the brain of either the same subject that the data were acquired, or on a nice smooth average brain surface. Rendering can be done on a smooth average in SPM99, and there is also the capability of doing the rendering on to the individual subjects MR. The rendering in SPM96 is done onto a single subject brain, which can be quite misleading as it is not the same subject that the data comes from. Depths behind the surface are derived from the single subject brain, which can cause problems as spatial normalisation is not exact. An activation that could be on the brain surface may appear a few mm in front of or behind the single subject brain surface. Another option within SPM99 is to use the brain extraction feature to identify the approximate brain surface of a segmented structural image. This would be saved as surf_*.mat. A routine similar the one attached can then be used to render the blobs on to this surface. Best regards JohnAttached file: fancy_rendering.m Gem 11: Editing Analyze HeadersOK, this mail is actually from Matthew Brett, but it references code written by John. Subject: Re: Scale factor From: Matthew Brett <matthewbrett@YAHOO.COM> Date: Tue, 22 May 2001 18:47:15 0700 (21:47 EDT) To: SPM@JISCMAIL.AC.UK Dear Masahiro, > I would like to incorporate scaling factors to many header files. > It would be nice if I can incorporate multiple scaling factors saved > in a text file into multiple header files. Is there a way to do > this? For the same problem, I have used an utility written by John Ashburner many eons ago, called dbedit: http://www.mrccbu.cam.ac.uk/Imaging/dbedit.html You need to script it somehow, obviously, but the line in the file setting your scale factor might look like: dbedit myfile.hdr dime.funused1=$myval where your scale factor is in the variable $myval I've also been playing with some perl functions to do this kind of thing, written by Andrew Janke: http://www.cmr.uq.edu.au/~rotor/software/ Best, Matthewdbedit tarball: dbedit.tar.gz Gem 12: fMRI Analysis ThresholdOK, this mail isn't from John either, and it doesn't even reference John's code, but it's useful info that I've needed on more than one occation. The difficulty is that the fMRI interface doesn't querry about masking; this email spells out how to overcome this. Subject: Re: explicit masking From: Stefan Kiebel <skiebel@fil.ion.ucl.ac.uk> Date: Tue, 27 Jun 2000 11:43:52 +0100 To: "Kevin J. Black" <kevin@npg.wustl.edu>, SPM <spm@mailbase.ac.uk> Dear Kevin, > Is it possible to instruct spm99 to search all voxels within a given > mask image rather than all above a fixed or a %mean threshold? Yes, with SPM99 it's possible to use several masking options. To recap, there are 3 sorts of masks used in SPM99: 1. an analysis threshold 2. implicit masking 3. explicit masking 1: One can set this threshold for each image to Inf to switch off this threshold. 2: If the image allows this, NaN at a voxel position masks this voxel from the statistics, otherwise the mask value is zero (and the user can choose, whether implicit masking should be used at all). 3: Use mask image file(s), where NaN (when image format allows this) or a nonpositive value masks a voxel. On top of this, SPM automatically removes any voxels with constant values over time. So what you want is an analysis, where one only applies an explicit mask. In SPM99 for PET, you can do this by going for the Full Monty and choosing Inf for the implicit mask and no 0thresholding. Specify one or more mask images. (You could also define a new model structure, controlling the way SPM for PET asks questions). With fMRI data/models, SPM99 is fully capable of doing explicit masking, but the user interface for fMRI doesn't ask for it. One way to do this type of masking anyway is to specify your model, choose 'estimate later' and modify (in matlab) the resulting SPMcfg.mat file. (see spm_spm.m lines 27  39 and 688  713). 1. Load the SPMcfg.mat file, set the xM.TH values all to Inf, set xM.I to 0 (in case that you have an image format not allowing NaN). 2. Set xM.VM to a vector of structures, where each structure element is the output of spm_vol. For instance: xM.VM = spm_vol('Maskimage'); 3. Finally, save by save SPMcfg xM append > If so, does the program define a voxel to be used as one which has > nonzero value in the named mask image? Not nonzero, but any positive value and unequal NaN. Note that you can specify more than one mask image, where the resulting mask is then the intersection of all mask images. Stefan Gem 13: Setting Blob Color Bar WindowSubject: Re: colour bar From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Thu, 28 Mar 2002 11:23:16 +0000 (06:23 EST) To: SPM@JISCMAIL.AC.UK > I would like to use the same colour bar for two different VBM experiments.. > Please, does anyone know how can I do this? If this is for the display of orthogonal views, then you can tweek the values that the blobs are scaled to by using a little bit of extra Matlab code. First of all, display an image, with superimposed blobs. Then in Matlab, type: global st st.vols{1}.blobs{1}.mx This will give the maximum intensity of the set of blobs. Then do the same with the other set of results. Find the largest of both results (suppose it is 5.6) and scale the blobs so they are displayed with this maximum by: global st st.vols{1}.blobs{1}.mx = 5.6; Best regards, John Gem 14: Down SamplingSubject: Re: resample to lower resolution From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Wed, 3 Jul 2002 12:35:22 +0100 (07:35 EDT) To: SPM@JISCMAIL.AC.UK > Is there a way to resample images of 512x512 spatial resolution down to > 256x256 resolution? If you copy and paste the following into Matlab, then it should do the job for you. Note that I haven't fully tested the code, but it worked on the one image I tried it on. Note that it only reduces the data using Fourier transforms in two directions. Combining slices in the 3rd direction is just by averaging. Best regards, John %* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % Select image to reduce V = spm_vol(spm_get(1,'*.img','Select image')); % Create new image header information VO = V; VO.fname = 'reduced.img'; VO.dim(1:3) = floor(VO.dim(1:3)/2); VO.mat = V.mat*[2 0 0 0.5 ; 0 2 0 0.5 ; 0 0 2 0.5 ; 0 0 0 1]; % Write the header VO = spm_create_image(VO); % Work out which bits of the Fourier transform to retain d1 = VO.dim(1); d2 = VO.dim(2); if rem(d1,2), r1a = 1:(d1+1)/2; r1b = []; r1c = []; r1d = (d1*2(d13)/2):d1*2; else, r1a = 1:d1/2; r1b = d1/2+1; r1c = d1*2d1/2+1; r1d = (d1*2d1/2+2):d1*2; end; if rem(d2,2), r2a = 1:(d2+1)/2; r2b = []; r2c = []; r2d = (d2*2(d23)/2):d2*2; else, r2a = 1:d2/2; r2b = d2/2+1; r2c = d2*2d2/2+1; r2d = (d2*2d2/2+2):d2*2; end; for i=1:VO.dim(3), % Fourier transform of one slice f = fft2(spm_slice_vol(V,spm_matrix([0 0 (i*21)]),VO.dim(1:2)*2,0)); % Throw away the unwanted region f1 = [f(r1a,r2a) (f(r1a,r2b)+f(r1a,r2c))/2 f(r1a,r2d) ([f(r1b,r2a) (f(r1b,r2b)+f(r1b,r2c))/2 f(r1b,r2d)]+... [f(r1c,r2a) (f(r1c,r2b)+f(r1c,r2c))/2 f(r1c,r2d)])/2 f(r1d,r2a) (f(r1d,r2b)+f(r1d,r2c))/2 f(r1d,r2d)]/4; % Fourier transform of second slice f = fft2(spm_slice_vol(V,spm_matrix([0 0 (i*2 )]),VO.dim(1:2)*2,0)); % Throw away the unwanted region f2 = [f(r1a,r2a) (f(r1a,r2b)+f(r1a,r2c))/2 f(r1a,r2d) ([f(r1b,r2a) (f(r1b,r2b)+f(r1b,r2c))/2 f(r1b,r2d)]+... [f(r1c,r2a) (f(r1c,r2b)+f(r1c,r2c))/2 f(r1c,r2d)])/2 f(r1d,r2a) (f(r1d,r2b)+f(r1d,r2c))/2 f(r1d,r2d)]/4; % Create a simple average of two FTs and do an inverse FT f = real(ifft2((f1+f2)))/2; % Write the plane to disk VO = spm_write_plane(VO,f,i); end; Gem 15: Computing Cerebral Volume (VBM)Subject: Re: smoothed modulated image From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Fri, 12 Apr 2002 13:45:01 +0000 To: SPM@JISCMAIL.AC.UK > considering a smoothed modulated image, which is the right > interpretation of the matrix: "each value of the matrix denotes the > volume, measured in mm3, of gray matter within each voxel" or "each > value of the matrix is proportional to the volume, measured in mm3, > of gray matter within each voxel" or something else? The contents of a modulated image are a voxel compression map multiplied by tissue belonging probabilities (which range between zero and one). The units in the images are a bit tricky to explain easily (so I would suggest you say that intensities are proportional). To find the volume of tissue in a structure in one of the modulated images, you sum the voxels for that structure and multiply by the product of the voxel sizes of the modulated image. The total volume of grey matter in the original image can be determined by summing the voxels in the modulated, spatially normalised image and multiplying by the voxel volume (product of voxel size). For example, try the following code for an original image and the same image after spatial normalisation and modulation. Providing the bounding box of the normalised image is big enough, then both should give approximately the same answer. V = spm_vol(spm_get(1,'*.img')) tot = 0; for i=1:V(1).dim(3), img = spm_slice_vol(V(1),spm_matrix([0 0 i]),V(1).dim(1:2),0); tot = tot + sum(img(:)); end; voxvol = det(V(1).mat)/100^3; % volume of a voxel, in litres tot = tot; % integral of voxel intensities tot*voxvol I hope the above makes sense. All the best, John Gem 16: Scripting FiguresOK, this isn't a John email, but rather a tip of my own that uses one of John's functions. When preparing a manuscript you often want to display a "blobs on brain" image, where a reference image underlies a colored significance image. You can do this within the SPM Results facility, but since you never get a figure right the first time, I prefer to do it on the command line. The code snippet blow scriptizes the blobsonbrain figure. You'll get a large orthgonal viewer in the graphics window, so it's then easy to print (or grab a screen snapshot) to then create your figure.
% Make sure to first clear the graphics window % Select images Pbck = spm_get(1,'*.img','Select background image') Psta = spm_get(1,'*.img','Select statistic image') % Set the threshold value Th = 4; % Create a new image where all voxels below Th have value NaN PstaTh = [spm_str_manip(Psta,'s') 'Th']; spm_imcalc_ui(Psta,PstaTh,'i1+(0./(i1>=Th))',{[],[],spm_type('float')},Th); % Display! spm_orthviews('image',Pbck,[0.05 0.05 0.9 0.9]); spm_orthviews('addimage',1,PstaTh) % Possibly, set the crosshairs to your favorite location spm_orthviews('reposition',[0 10 10])This assumes that you just want to threshold your image based on a single intensity threshold. To make it totally scripted, replace the spm_get calls with hard assignments.
Gem 17: Origin maddnessA source of confusion is where the origin (the [0,0,0] location of an image) is stored. When there is no associated .mat file, the origin is read from the Analyze originator field. If this is zero it is assumed to match the center of the image field of view. If there is a .mat file, then the origin is the first three values of M\[0 0 0 1]'where M is the transformation matrix in the .mat file. One limitation is that the origin stored in the Analyze header is a (short) integer, and so cannot represent an origin with fractional values. To set the origin to specific, fractional value, use this code snippet: Orig = [ x y z ]; % Desired origin in units of voxels P = spm_get(Inf,'*.img'); % matrix of file names for i=1:size(P,1) M = spm_get_space(deblank(P(i,:))); R = M(1:3,1:3); % Set origin M(1:3,4) = R*Orig(:); spm_get_space(deblank(P(i,:)),M); end Gem 18: log10 Pvalues from T imagesPvalue images are difficult to visualize since "important" values are small and clumped near zero. A log10 transformation makes for much better visualization while still having interpretability (e.g. a value of 3 cooresponds to P=0.001). This function, T2nltP, will create log10 Pvalue image based on either a contrast number (which must be a T contrast) or a T statistic image and the degrees of freedom. (See also the equivalent SPM2 function.)
function T2nltP(a1,a2) % Write image of log10 Pvalues for a T image % % FORMAT T2nltP(c) % c Contrast number of a T constrast (assumes cwd is a SPM results dir) % % FORMAT T2nltP(Timg,df) % Timg Filename of T image % df Degrees of freedom % % % As per SPM convention, T images are zero masked, and so zeros will have % Pvalue NaN. % % @(#)T2nltP.m 1.2 T. Nichols 03/07/15 if nargin==1 c = a1; load xCon load SPM xX if xCon(c).STAT ~= 'T', error('Not a T contrast'); end Tnm = sprintf('spmT_%04d',c); df = xX.erdf; else Tnm = a1; df = a2; end Tvol = spm_vol(Tnm); Pvol = Tvol; Pvol.dim(4) = spm_type('float'); Pvol.fname = strrep(Tvol.fname,'spmT','spm_nltP'); if strcmp(Pvol.fname,Tvol.fname) Pvol.fname = fullfile(spm_str_manip(Tvol.fname,'H'), ... ['nltP' spm_str_manip(Tvol.fname,'t')]); end Pvol = spm_create_image(Pvol); for i=1:Pvol.dim(3), img = spm_slice_vol(Tvol,spm_matrix([0 0 i]),Tvol.dim(1:2),0); img(img==0) = NaN; tmp = find(isfinite(img)); if ~isempty(tmp) img(tmp) = log10(max(eps,1spm_Tcdf(img(tmp),df))); end Pvol = spm_write_plane(Pvol,img,i); end; Gem 19: VBM modulation scriptThis is the famed script to modulate spatially normalized probability images. For gray matter probability images, modulated images have units of gray matter volume per voxel, instead of gray matter concentration adjusted for differences in local brain size. In John's own words, here's a better explaination. Note that script below is SPM99 specific. In SPM2, it is done via spm_write_sn(V,prm,'modulate') (See spm_write_sn help fore more.) Date: Thu, 27 Jul 2000 14:28:39 +0100 (BST) Subject: Re: matlab question & VBM From: John Ashburner <john@fil.ion.ucl.ac.uk> To: spm@mailbase.ac.uk, x.chitnis@iop.kcl.ac.uk  Following on from John Ashburner's recent reply, is there a matlab function  that enables you to adjust spatially normalised images in order to preserve  original tissue volume for VBM? The function attached to this email will do this. Type the following bit of code into Matlab to run it: Mats = spm_get(Inf,'*_sn3d.mat','Select sn3d.mat files'); Images = spm_get(size(Mats,1),'*.img','Select images to modulate'); for i=1:size(Mats,1), spm_preserve_quantity(deblank(Mats(i,:)),deblank(Images(i,:))); end; [...] Best regards, JohnThe attached script is here spm_preserve_quantity.m
Gem 20: Creating customized templates and priorsFrom: John Ashburner <john@FIL.ION.UCL.AC.UK> Subject: Re: Help for constructing Template images Date: Wed, 18 Dec 2002 16:22:59 +0000 To: SPM@JISCMAIL.AC.UK > What are the advantages of customized template images > in VBM analysis? Customised templates are useful when: 1) The contrast in your MR images is not the same as the contrast used to generate the existing templates. If the contrast is different, then the mean squared cost function is not optimal. However, for "optimised" VBM this only really applies to the initial affine registration that is incorporated into the initial segmentation. Contrast differences are likely to have a relatively small effect on the final results. 2) The demographics of your subject population differ from those used to generate the existing templates and prior probability images. For example, serious problems can occur if your subjects have very large ventricles. In these data, there would be CSF in regions where the existing priors say CSF should not exist. This would force some of the CSF to be classified as white matter, seriously affecting the intensity distribution that is used to model white matter. This then has negative consequences for the whole of the segmentation. > Can any one please explain the detailed steps to > construct a customized template image (gray and white > matter images) for VBM analysis? The following script is one possible way of generating your own template image. Note that it takes a while to run, and does not save any intermediate images that could be useful for quality control. Also, if it crashes at any point then it is difficult to recover the work it has done so far. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *make_template.m
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * You may also wish to do some manual editing of the images afterwards  especially to remove extraskull CSF. When everything has finished, simply smooth the images by 8mm and call them templates and prior probability images. You can modify the default priors for the segmentation step in order that the customised ones are used. This can be done either by changing spm_defaults.m, or by typing the following in Matlab: spm_defaults global defaults defaults.segment.estimate.priors = ... spm_get(3,'*.IMAGE','Select GM,WM & CSF priors'); Note that this will be cleared if you reload the defaults. This could be done when you start spm, reset the defaults or if the optimised VBM script is run, as it calls spm_defaults.m. Alternatively the optimised VBM script could be modified to include the above. Note that I have only tried the script with three images, so I don't have a good feel for how robust it is likely to be. > > Please let me know the number of subjects required to > construct one? Its hard to say, but more is best. The 8mm smoothing means that you can get away with slightly fewer than otherwise. Best regards, JohnNote: The make_template is an updated version of what was originally posted. It is current as of Sep 9, 2003. Gem 21: ImCalc ScriptAs posted, this snippet is for SPM2; I've edited it to work with SPM99. Subject: Re: a script to use ImaCal From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Thu, 16 Oct 2003 11:24:17 +0000 (07:24 EDT) To: SPM@JISCMAIL.AC.UK > Brain image for each subject to mask out CSF signal was generated by > using MPR_seg1.img (i1) and MPR_seg2.img (i2) with (i1+i2)>0.5 in > ImaCal.(called brainmpr.img for each subject) > > Then I have more than twohundred maps, which need to mask out > CSF. I think I can use ImaCal again with selecting brainmpr.img (i1) > and FAmap.img (i1), and then calculating (i1.*i2) to generate a new > image named as bFAmap.img. Unfortunately, if I use the ImaCal, it > take so long time to finish all subjects. Could anyone have a script > to generate a multiplication imaging with choosing a brain > image(brainmpr.img, i1) and a map image (FAmap.img, i2) and writing > an output image (bFAmap.img) from i1.*i2? You can do this with a script in Matlab. Something along the lines of the following should do it: P1=spm_get(Inf,'*.img','Select i1'); P2=spm_get(size(P1,1),'*.img','Select i2'); for i=1:size(P1,1), P = strvcat(P1(i,:),P2(i,:))); Q = ['brainmpr_' num2str(i) '.img']; f = '(i1+i2)>0.5'; flags = {[],[],[],[]}; Q = spm_imcalc_ui(P,Q,f,flags); end; Note that I have not tested the above script. I'm sure you can fix it if it doesn't work. Best regards, John Gem 22: spm_orthviews tricksIMHO, after spatial normalization, John's key contribution to SPM is spm_orthviews, the function behind the 'Check Reg' button which let's you view many volumes simeltaneously. Some of the Gems use spm_orthviews (e.g. Gems 1 and 16) but listed here are some generally useful tricks for spm_orthviews. These tricks are useful anytime a threeview orthogonal slice view is shown, whether from useing 'Check Reg', 'Display' or when overlaying blobs from the 'Results' window.
Gem 23: Tabulating T StatisticsSubject: Fwd: Re: tabulating all statistics From: John Ashburner <john@FIL.ION.UCL.AC.UK> Date: Tue, 1 Jul 2003 11:47:07 +0000 (07:47 EDT) To: SPM@JISCMAIL.AC.UK > I was wondering if it would be possible to write t values for an > entire volume into a file? Try this: fid=fopen('tvalues.txt','w'); P=spm_get(1,'*.img','Select statistic image'); V=spm_vol(P); [x,y,z] = ndgrid(1:V.dim(1),1:V.dim(2),0); for i=1:V.dim(3), z = z + 1; tmp = spm_sample_vol(V,x,y,z,0); msk = find(tmp~=0 & finite(tmp)); if ~isempty(msk), tmp = tmp(msk); xyz1=[x(msk)'; y(msk)'; z(msk)'; ones(1,length(msk))]; xyzt=V.mat(1:3,:)*xyz1; for j=1:length(tmp), fprintf(fid,'%.4g %.4g %.4g\t%g\n',... xyzt(1,j),xyzt(2,j),xyzt(3,j),tmp(j)); end; end; end; fclose(fid); best regards, JohnAs noted in a 2 Feb 2004 email, to eliminate all voxels below a certain threshold, change msk = find(tmp~=0 & finite(tmp));to msk = find(tmp>0 & finite(tmp));
Gem 24: Unnormalizing a pointThis script of John's will find the corresponding coordinate in unnormalised image: get_orig_coord.m
