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;