% Sets up an example for a 3D grid resampling operation

% Old grid dimensions
M = 512;
N = 512;
Q = 10;

dx = 0.1;
dy = 0.1;
dz = 0.2

% resolution after resampling (if empty, resolution is unchanged
dxNew = 2*dx;
dyNew = 2*dy;
dzNew = 0.5*dz;

% create the phantom
Pslice = phantom(M);
P = zeros(M,N,Q,'single');
for k=1:Q
    P(:,:,k) = Pslice;
end

% original data coordinates
x = [-(M-1)/2:(M-1)/2]*dx;
y = [-(N-1)/2:(N-1)/2]*dy;
z = [-(Q-1)/2:(Q-1)/2]*dz;

% new coordinate system
xNew = [x(1):dxNew:x(end)];
yNew = [y(1):dyNew:y(end)];
zNew = [z(1):dzNew:z(end)];

% do the resampling
Pp = gridResample3D(x,y,z,P,xNew,yNew,zNew,'linear');

% plot results
figure; imagesc(P(:,:,ceil(Q/2))); colormap gray;
title('Original image','FontSize',14);

figure; imagesc(Pp(:,:,ceil(numel(zNew)/2))); colormap gray;
title('Resampled image','FontSize',14);