Dear everyone,
just a little addendum to Johns mail below.
The problem appears for "single scan per subject" designs such as e.g. when
comparing glucose consumption (FDG), blood flow or gray matter denstity between
different populations. In particular it appears when in these designs there is only
partial coverage of the brain in the axial direction.
When there are multiple conditions per subject the problem does NOT appear in the
"condition effects" thanks to the masking that is implicit in the realignment. Note
that this is the "normal" type of design for e.g. fMRI, so this problem is really
confined to a small subpopulation of studies.
In addition, for VBM type studies (where this could IN PRINCIPLE be a problem) it
would be rare indeed to find scans with only partial coverage of the brain.
So, the problem is mainly to be found for PET (old scanners with 10cm FOV, or crap
positioning) and SPECT using FDG or some rCBF marker.
PET, and in particular SPECT is associated with a poor signal to noise ratio
(remember that the masking operation should be performed prior to spatial
smoothing). Together with some of the particulars about FBP this means that there
will be "valid" zeros in the image volumes. These valid zeros originate in areas
where the expectation is close to zero (CSF) or where the expectation is at least
close to zero relative to the noise (white matter). Then we of course have the
"invalid" zeros, i.e. those originating from not having measured at a certain
location in the normalised space (MNI). The masking in spm_mask will remove from
consideration any voxels where the intensity value is zero in at least one of the
scans. What we would really like to do is to remove any voxels where there is an
INVALID zero in at least one of the scans. Depending on the specifics the
distinction may be of quite considerable importance. For a reasonably large set of
scans with poor SNR (even worse if it has at some stage been stored in an 8bit
format) the resulting mask might look more like a swiss cheese than a brain.
To avoid that one needs some means of distinguishing between valid and invalid
zeros. One distinguishing feature is that any zero that originates from FOV problems
will be found together with other zeros (the volume outside the scanned volume). So,
if we find a zero and find that it is connected to other zeros that is a strong
indication it is an invalid zero. That is the reasoning behind the modified
spm_mask.m that I have attached to this mail.
Hence, if you try spm_mask.m and find that your brains look like they all have BSE,
give the modified spm_mask.m a try. Usage is same as described by John below.
If you are not on a Linux machine you need to compile the get_connected_zeros_dtj.c
file. In Matlab just type
mex get_connected_zeros_dtj.c
at the command prompt. It should work.
All the usual disclaimers (I hear voices in my head and only do what they tell me,
etc etc) apply. I have tested it "a bit" though.
Good luck Jesper
> > Meanwhile, are you sure of the below?
> >
> > > > I am analyzing a group of SPECT scans consisting of about 50 scans. The
> > > > problem is that a portion of these scans (14 of them) are slightly to
> > > > mildly cropped (i.e. part of the cerebellum is cut off). I am
> > > > wondering how this will effect my results and/or what I can do about
> > > > this problem. Based on previous emails, it looks like one suggestion
> > > > might be to create an object mask, if that is a possible solution how
> > > > do I go about creating one?
> > >
> > > It shouldn't have too much of an effect, other than the region of brain
> > > included in the final analysis will be smaller. I would only suggest any
> > > kind of "object masking" if there is a region within the image volume
> > > that is not how it should be. I am assuming that the part of the
> > > cerebellum that is cropped is outside the image volume.
> >
> > Lets say you have n subjects in group A and m subjects in group B. Also say
> > that there has been a systematic difference in axial positioning such that
> > all subjects in group A has been scanned 1cm (2 slices) more cranially than
> > group B. After spatial normalisation this means there will be a few planes
> > containing non-zero values in group A, but with zero-values in group B (or
> > is there some way to apply masking in the spatial normalisation (similar to
> > the realignment masking, but across subjects rather than scans) ?).
> > If we did the statistics there and then it would not be a problem since the
> > implicit intensity based mask would remove all those voxels. But we don't.
> > Instead we spatially smooth the data. If we consider a voxel at the very
> > edge of the area where we have data for both subjects then we will have a
> > different mixture of brain and non-brain values within the filter-kernel
> > for the two groups. Also, especially for the "inner part of the edge" where
> > we have all brain values for group A and only a few non-brain values for
> > group B within the kernel it is likely that voxels will pass the conditions
> > for the implicit mask. Statistics on this will now compare apples and
> > pears.
> >
> > Knowing you I know you will have thought of this, and I am sure there is
> > some clever way in SPM to avoid this problem. However, I haven't found it
> > and I suspect perhaps some others haven't either. I have actually helped
> > some people (in Uppsala and Huddinge who did SPECT and FDG PET respectively
> > on different dementia groups) with this very problem. When they just
> > performed a "straightforward" analysis there were quite massive differences
> > at the edge of the "actual axial FOV" of the kind one would expect from the
> > reasoning above. By AND-masking all the images prior to spatial smoothing
> > these differences dissapeared.
> > Have I been doing something stupid (again) to "fix" something that SPM
> > could have handled?
>
> Good point. I forgot about this (got mixed up by object masking in spatial
> normalisation). The data probably should be masked between the spatially
> normalised data being written, and it being smoothed. The spm_mask.m
> routine can be used to do this:
>
> P1 = spm_get(Inf,'*.img','Images to define mask from');
> P2 = spm_get(Inf,'*.img','Images to mask');
> spm_mask(P1,P2);
>
> The 'Images to define mask from' would be an example spatially normalised image
> from each subject, whereas 'Images to mask' would be all the images to enter into
> the analysis. The output masked images are prefixed by 'm'.
>
> Best regards,
> -John
>
> --
> Dr John Ashburner.
> Functional Imaging Lab., 12 Queen Square, London WC1N 3BG, UK.
> tel: +44 (0)20 78337491 or +44 (0)20 78373611 x4381
> fax: +44 (0)20 78131420 http://www.fil.ion.ucl.ac.uk/~john
function indx = get_connected_zeros(f,ncon)
% FORMAT indx = get_connected_zeros(f,ncon)
% f - A 2D or 3D image(-volume)
% ncon - Requested no. of connected voxels
%
% indx - Array of indicies to zero-elements
% in f that are connected (six-con)
% to at least ncon other zero-elements.
%
% This routine takes a 2D or 3D image volume
% and returns an array of indicies to zero-
% elements. It can be used for example to create
% masks from realigned images without getting
% spurious zero-values.
% Note that if one plane has been shifted out
% completely any spurious zeros in the next
% plane will be connected to that plane of
% zeros, and will therefore be masked away.
% To avoid that one might want to use the
% 2D facility.
%
%______________________________________________
% Jesper Andersson 22/1-02
indx = get_connected_zeros_dtj(f,ncon);
return
/*
C-code for small utility that finds, in
an image volume, all zeros that are
connected to at least ncon other zeros.
Jesper Andersson
*/
#include "mex.h"
#include <math.h>
#define MAX_N 100
#ifndef MAX
#define MAX(a,b) (((a)>(b)) ? (a) : (b))
#endif
#ifndef MIN
#define MIN(a,b) (((a)<(b)) ? (a) : (b))
#endif
#ifndef INDX
#define INDX(i,j,k,dim) (((k)-1)*dim[0]*dim[1]+((j)-1)*dim[0]+((i)-1))
#endif
/* Utility function. */
int is_n_connected(double *f,
unsigned int dim[3],
int ii,
int ij,
int ik,
int in,
int recursive)
{
static int ijk[MAX_N]; /* Map of visited voxels. */
static int n;
int i = 0;
if (!recursive) /* Initiate recursive search. */
{
n = 0;
memset(ijk,0,MAX_N*sizeof(int));
}
if (n == in) return(n); /* Return if criteron fullfilled. */
for (i=0; i<n; i++) /* See if new place. */
{
if (ijk[i] == INDX(ii,ij,ik,dim)) break;
}
if (i<n)
{
return(n); /* Old place, return. */
}
else
{
ijk[n] = INDX(ii,ij,ik,dim);
n++;
if ((ii+1)<=dim[0] && f[INDX(ii+1,ij,ik,dim)] == 0)
{
is_n_connected(f,dim,ii+1,ij,ik,in,1);
}
if ((ii-1)>0 && f[INDX(ii-1,ij,ik,dim)] == 0)
{
is_n_connected(f,dim,ii-1,ij,ik,in,1);
}
if ((ij+1)<=dim[1] && f[INDX(ii,ij+1,ik,dim)] == 0)
{
is_n_connected(f,dim,ii,ij+1,ik,in,1);
}
if ((ij-1)>0 && f[INDX(ii,ij-1,ik,dim)] == 0)
{
is_n_connected(f,dim,ii,ij-1,ik,in,1);
}
if ((ik+1)<=dim[2] && f[INDX(ii,ij,ik+1,dim)] == 0)
{
is_n_connected(f,dim,ii,ij,ik+1,in,1);
}
if ((ik-1)>0 && f[INDX(ii,ij,ik-1,dim)] == 0)
{
is_n_connected(f,dim,ii,ij,ik-1,in,1);
}
}
return(n);
}
/* Function doing the job. */
int get_connected_zeros_dtj(unsigned int dim[3],
double *f,
unsigned int ncon,
double *zeros)
{
int nz = 0;
int i, j, k;
for (k=1; k<=dim[2]; k++)
{
for (j=1; j<=dim[1]; j++)
{
for (i=1; i<=dim[0]; i++)
{
if (f[INDX(i,j,k,dim)] == 0)
{
if (is_n_connected(f,dim,i,j,k,ncon,0)==ncon)
{
zeros[nz++] = INDX(i,j,k,dim) + 1;
}
}
}
}
}
return(nz);
}
/* Gateway function with error check. */
void mexFunction(int nlhs, /* No. of output arguments */
mxArray *plhs[], /* Output arguments. */
int nrhs, /* No. of input arguments. */
const mxArray *prhs[]) /* Input arguments. */
{
unsigned int ndim;
unsigned int dim[3];
int nzeros = 0;
int ncon = 0;
char str[256];
double *zeros = NULL;
if (nrhs == 0) mexErrMsgTxt("usage: indx = get_connected_zeros_dtj(f,ncon)");
if (nrhs != 2) mexErrMsgTxt("get_connected_zeros_dtj: 2 input arguments required");
if (nlhs != 1) mexErrMsgTxt("get_connected_zeros_dtj: 1 output argument required");
/* Get data. */
if (!mxIsNumeric(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]) || !mxIsDouble(prhs[0]))
{
mexErrMsgTxt("get_connected_zeros_dtj: f must be numeric, real, full and double");
}
if ((ndim = mxGetNumberOfDimensions(prhs[0])) > 3)
{
mexErrMsgTxt("get_connected_zeros_dtj: f must be 2 or 3 dimensional");
}
if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) || mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1]))
{
mexErrMsgTxt("get_connected_zeros_dtj: ncon must be numeric, real, full and double");
}
if (mxGetN(prhs[1])!=1 || mxGetM(prhs[1])!=1)
{
mexErrMsgTxt("get_connected_zeros_dtj: ncon must be a scalar");
}
ncon = ((int) mxGetScalar(prhs[1]));
if (ncon > MAX_N)
{
sprintf(str,"get_connected_zeros_dtj: ncon must be less than %d",MAX_N);
mexErrMsgTxt(str);
}
memcpy(dim,mxGetDimensions(prhs[0]),ndim*sizeof(int));
if (ndim == 2) dim[2] = 1;
zeros = (double *) malloc(dim[0]*dim[1]*dim[2]*sizeof(double));
if ((nzeros = get_connected_zeros_dtj(dim,mxGetPr(prhs[0]),ncon,zeros)) < 0)
{
free(zeros);
mexErrMsgTxt("get_connected_zeros_dtj: something went wrong");
}
/* Allocate memory for output and copy data. */
if (nzeros == 0)
{
plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
}
else
{
plhs[0] = mxCreateDoubleMatrix(nzeros,1,mxREAL);
memcpy(mxGetPr(plhs[0]),zeros,nzeros*sizeof(double));
free(zeros);
}
return;
}
function spm_mask(P1,P2,ncon)
% Masks Images.
% FORMAT spm_mask(P1,P2)
% P1 - matrix of input image filenames from which
% to compute the mask.
% P2 - matrix of input image filenames on which
% to apply the mask.
% ncon - Number of zeros a zero must be connected
% to to qualify as invalid.
% The masked images are prepended with the prefix `m'.
% In addition there will be a file named mmask.img
% (to separate it from the mask.img generated by the
% stats module in SPM) which shows the "largest common
% volume" for the study.
%
% If any voxel in the series of images is zero (for data types without
% a floating point representation), and is connected to a set of other
% zeros, or does not have a finite value (for floating point and double
% precision images), then that voxel is set to NaN or zero in all the images.
%
%_______________________________________________________________________
% @(#)spm_mask.m 2.7 John Ashburner 99/04/01
%
% Modfied 02/10/01 Jesper Andersson
%
if nargin==0,
P1=spm_get(Inf,'.img','Images to compute mask from');
P2=spm_get(Inf,'.img','Images to apply mask to');
ncon = 10;
end;
if nargin==1,
P2 = P1;
ncon = 10;
end;
if nargin==2
ncon = 10;
end
P1=spm_vol(P1);
P2=spm_vol(P2);
m1=prod(size(P1));
m2=prod(size(P2));
spm_progress_bar('Init',P1(1).dim(3),'Finding mask','planes completed')
if P1(1).dim(4)==2 | P1(1).dim(4)==4 | P1(1).dim(4)==8 |...
P1(1).dim(4)==512 | P1(1).dim(4)==1024 | P1(1).dim(4)==2048
msk = zeros(P1(1).dim(1:3));
[x,y] = ndgrid(1:P1(1).dim(1),1:P1(1).dim(2));
xyz = [x(:) y(:) ones(prod(P1(1).dim(1:2)),1)];
for sl=1:P1(1).dim(3)
tmp = ones(prod(P1(1).dim(1:2)),1);
for i=1:length(P1)
f = reshape(spm_sample_vol(P1(i),xyz(:,1),xyz(:,2),...
sl*xyz(:,3),0),P1(1).dim(1:2));
indx = get_connected_zeros(f,ncon);
tmp(indx) = 0;
end
msk(:,:,sl) = reshape(tmp,P1(1).dim(1:2));
spm_progress_bar('Set',sl);
end
else % Use John's original code
msk = zeros(P1(1).dim(1:3));
for j=1:dim(3),
tmp = zeros(dim(1:2));
M1 = spm_matrix([0 0 j]);
% Load slice j from all images
for i=1:m1
img = spm_slice_vol(P1(i),M1,dim(1:2),[1 0]);
tmp = tmp + finite(img);
end;
msk(:,:,j) = tmp;
spm_progress_bar('Set',j);
end
end
spm_progress_bar('Clear');
%
% Write mask
%
[mypath,myfile,myext] = fileparts(P1(1).fname);
mP = P1(1);
mP.fname = fullfile(mypath,'mmask.img');
mP.descrip = 'Mask generated by spm_mask';
mP = spm_write_vol(mP,msk);
%
% Write masked images.
%
spm_progress_bar('Init',length(P2),'Writing masked images','volumes completed')
%
% Write masked data
%
for i=1:length(P2)
oP = P2(i);
[mypath,myfile,myext] = fileparts(P2(i).fname);
oP.fname = fullfile(mypath,['m' myfile myext]);
oP.descrip = [oP.descrip ' masked'];
f = spm_read_vols(P2(i));
oP = spm_write_vol(oP,f.*msk);
spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');
return;
|
