How to interpolate with different voxel sizes?

Hi all,

I want to load in a PET image which has the dimensions 171x171x127 voxels, with each voxel being 4.17x4.17x2.03 mm. Within this image, there is for example a sphere with a radius of 30 voxels. I want to interpolate this to a grid in which each voxel has the dimensions of 3x3x3 mm, in which this sphere adjusts to the new grid, where the radius of the sphere then has a different number of voxels as a radius (in the new grid that would be 10 voxels then). I have the following code already:

x=171; y=171; z=121;
PETimage=zeros(x,y,z);
rad=30;
PETimage(round(.5*x)-rad:round(.5*x)+rad,round(.5*y)-rad:round(.5*y)+rad,round(.5*z)-rad:round(.5*z)+rad)=1;

size_x=4.17;
size_y=4.17;
size_z=2.03;

x_grid=x*size_x; x_newgrid=floor(x_grid/3);
y_grid=y*size_y; y_newgrid=floor(y_grid/3);
z_grid=z*size_z; z_newgrid=floor(z_grid/3);
newgrid=zeros(x_newgrid,y_newgrid,z_newgrid);

I think I now have to use either interpol3 or griddedInterpolant, but I don't know how. Can anybody help me out?

Greetings, Lieke

Answer 1

[William Rose]

@Lieke Pullen,

I would use interp3(). First you must create the grids X,Y,Z (the old coordinates of the voxels) and Xq,Yq,Zq (the coordinates of the new voxels) with meshgrid. Then you call interp3(). See code below.

nx=171; ny=171; nz=121;
dx=4.17; dy=4.17; dz=2.03;
x=(0:nx-1)*dx; %vector of x coords
y=(0:ny-1)*dy; %vector of y coords
z=(0:nz-1)*dz; %vector of z coords
[X,Y,Z]=meshgrid(x,y,z); %make three 3-D meshes X,Y,Z
dx2=3.0; dy2=3.0; dz2=3.0;
nx2=floor(nx*dx/dx2); 
ny2=floor(ny*dy/dy2);
nz2=floor(nz*dz/dz2);
x2=(0:nx2-1)*dx2; %vector new of x coords
y2=(0:ny2-1)*dy2; %vector of new y coords
z2=(0:nz2-1)*dz2; %vector of new z coords
[X2,Y2,Z2]=meshgrid(x2,y2,z2); %make three 3-D meshes X2,Y2,Z2
PETimage=zeros(nx,ny,nz);
rad=30;
PETimage(round(.5*nx)-rad:round(.5*nx)+rad,round(.5*ny)-rad:round(.5*ny)+rad,round(.5*nz)-rad:round(.5*nz)+rad)=1;
PETimage2=zeros(nx2,ny2,nz2);
PETimage2=interp3(X,Y,Z,PETimage,X2,Y2,Z2);
fprintf('Size of image 1=%dx%dx%d. Size of image 2 = %dx%dx%d.\n',size(PETimage),size(PETimage2));

Size of image 1=171x171x121. Size of image 2 = 237x237x81.

Try it. Good luck.

[William Rose]

@Lieke Pullen,

The line

PETimage2=zeros(nx2,ny2,nz2);

in my code is unnecessary, and can be deleted.

[William Rose]

@Lieke Pullen,

I used "floor()" when calculating nx2, ny2, nz2, to assure that the new coordinates would not extend byond the old coordinates. That will assure that we are interpolating, and not extrapolating.

[Lieke Pullen]

Thank you William, this was indeed what I was looking for! I didn't know how to incorporate the voxelsizes, but this works. Thanks again!

[Lieke Pullen]

Hi @William Rose, I have another question regarding the simulation of the sphere in the grid. For example, I have a sphere in the X,Y,Z meshgrid with a radius of 30 mm. I want to interpolate this on the 'new' grid, so that the radius is still 30 mm, but the sphere now consists of another amount of voxels with different dimensions with respect to the original grid. How do I do this. I now have the following code

nx=171; ny=171; nz=121;
dx=4.17; dy=4.17; dz=2.03;
x=(0:nx-1)*dx; %vector of x coords
y=(0:ny-1)*dy; %vector of y coords
z=(0:nz-1)*dz; %vector of z coords
[X,Y,Z]=meshgrid(x,y,z); %make three 3-D meshes X,Y,Z

dx2=3.0; dy2=3.0; dz2=3.0;
nx2=floor(nx*dx/dx2); 
ny2=floor(ny*dy/dy2);
nz2=floor(nz*dz/dz2);
x2=(0:nx2-1)*dx2; %vector new of x coords
y2=(0:ny2-1)*dy2; %vector of new y coords
z2=(0:nz2-1)*dz2; %vector of new z coords

radius=30 %in mm
sphere_tumor=sqrt(X.^2 + Y.^2 + Z.^2) <= radius; sphere_tumor=double(sphere_tumor); %creating the sphere
sphere_PET=interp3(X,Y,Z,sphere_tumor,X2,Y2,Z2);

I am mainly struggling with the different voxel sizes. Thank you in advance!

[William Rose]

@Lieke Pullen,

Your original code used the line

%PETimage(round(.5*nx)-rad:round(.5*nx)+rad,round(.5*ny)-rad:round(.5*ny)+rad,round(.5*nz)-rad:round(.5*nz)+rad)=1;

The line above creates a rectangular tumor, 2*rad+1 voxels in each direction, in the original coordinates. To create a spherical tumor of radius r, centered at (x0,y0,z0), in the original coordinates, create a function that returns True=1 if (x,y,z) is inside or on the sphere, and 0 otherwise. Then make a loop (actually 3 nested loops for x,y,z) to call the function for each voxel in the cube that bounds the sphere. Set PETimage to the tumor denisty in the voxels where the function returns True. The code below shows examples of calls to the function, folowed by the function itself.

After you create the tumor in the orignal image array, it will be correctly represented in the new array.

%examples of calls to isInside are below, followed by the function
isInside([1,0,0],[0,0,0],1)

ans = logical
   1

isInside([1,-1,0],[0,0,0],1)

ans = logical
   0

isInside([27,14,26],[20,20,20],11)

ans = logical
   1

function x=isInside(p,p0,r)
%ISINSIDE determines if a point is (in or on) a specified sphere, or not
%Inputs: p=[x,y,z]=position, p0=[x0,y0,z0]=sphere center, r=sphere radius
%Output: x=True if ||p-p0||<=r, else x=False
x=sum((p-p0).^2)<=r^2;
end

Try it.