Subject: 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-(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;