A 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 thereisa .mat file, then the origin is the first three values of

        M\[0 0 0 1]'
                

whereMis 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

OK, 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 blobs-on-brain 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.

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

Subject: Re: resample to lower resolution
From: John Ashburner <john@FIL.ION.UCL.AC.UK&gt
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-(d1-3)/2):d1*2;
else,         r1a = 1:d1/2;     r1b = d1/2+1; 
         r1c = d1*2-d1/2+1; r1d = (d1*2-d1/2+2):d1*2;
end;
if rem(d2,2), r2a = 1:(d2+1)/2; r2b = [];     
         r2c = [];          r2d = (d2*2-(d2-3)/2):d2*2;
else,         r2a = 1:d2/2;     r2b = d2/2+1; 
         r2c = d2*2-d2/2+1; r2d = (d2*2-d2/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*2-1)]),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;

Subject: Re: colour bar
From: John Ashburner &ltjohn@FIL.ION.UCL.AC.UK&gt
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

OK, 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&gt
Date: Tue, 27 Jun 2000 11:43:52 +0100
To: "Kevin J. Black" <kevin@npg.wustl.edu>, SPM <spm@mailbase.ac.uk&gt

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 non-positive 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 0-thresholding. 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

OK, this mail is actually from Matthew Brett, but it references code written by John.

Subject: Re: Scale factor
From: Matthew Brett 
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.mrc-cbu.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,

Matthew
  

dbedit tarball:dbedit.tar.gz

Subject: Re: rendering SPM96
From: John Ashburner <john@FIL.ION.UCL.AC.UK&gt
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,:) &lt (z1+20) &amp xyz(3,:) &gt (z1-5));
to something like:
        msk = find(xyz(3,:) &lt (z1+5) &amp xyz(3,:) &gt (z1-5));


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
-John
  

Attached file:fancy_rendering.m

Subject: 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.insa-lyon.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

This 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&gt
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,
-John

The 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;