How can I get an interpolated value "between" voxels?

Beginner Questions
, ,
1

Hello,
I have a 3D image in raw format, the .mhd file is similar to this one:

ObjectType = Image
NDims = 3
ElementSpacing = 0.031473 0.031473 0.031473
DimSize = 100 100 100
ElementType = MET_SHORT
ElementDataFile = phantom.raw

I would like to get an interpolated value at non-integer voxel coordinate.
With this program I can get the value for example at coordinate 25,50,50:

import SimpleITK as sitk
import numpy as np

# Reads the image using SimpleITK
itkimage = sitk.ReadImage('phantom.mhd')

# Convert the image to a  numpy array first and then shuffle the dimensions to get axis in the order z,y,x
ct_scan = sitk.GetArrayFromImage(itkimage)

print(ct_scan[25,50,50])

But what if I need the value say at 25.1,49.9,50.0?

I think I need to use a ResampleImageFilter and I would like to use the interpolator of type WindowedSincInterpolateImageFunction but I am not capable to adapt the examples I have found (i.e.
https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/Filtering/ResampleImageFilter8.cxx) because they seem to be always related to a transformation like a rotation or a translation of the whole image.

My final goal is to get not just one interpolated value but N values for example along a curve.

1 Like
2

I have a couple thoughts:

  • A new method could be added to the sitk::Image class. Perhaps something like GetInterpolatedPixel(std::vector &, InterpolatorEnum). It may need to have different return types for many of the pixel types, which adds a little complication.
  • For you case of interpolating along a curve, the ResampleImageFilter may be a good approach. If you consider your image as input, and the output a 3D image if Nx1x1 , where N is the number of values you want. The transform maps the outputs from the output to the input. If you can come up with this transform you can use the resample filter.
3

Thank you for your hints!

For the first thoughts, adding a method to Image class, I am afraid I can not help a lot as a newbie user :slight_smile:.

For the second thought: let’s say I am interested in the line L described here with N points on it:

N=10
L=[]
begin=np.array([50,50,50])
end = np.array([90,90,90])
ts=np.linspace(0,1,N)
for t in ts:
	L.append(begin+t*(end-begin))

now, how can I define a 3D image with shape Nx1x1?

4

For the case of sampling along a line, a 3d rotational matrix with translation can be used to perform the transformation. With a little research you should be able to find a formula, per construct such a transform. Keep in mind the transformation maps the point from physical space to physical space ( not between index spaces! ), so the output image’s geometry is important.

1 Like
5

Thank you.
This is what I understood, even if this is not strictly related to ITK maybe you can help me anyway:
given a line L in 3D (for example given two points in mm); first find the translation T that bring L passing through the origin, call L1 this new line passing through the origin; then find the rotation R that overlap L1 to, say, axis x; then apply T and R to my 3D image using a ResampleImageFilter with a WindowedSincInterpolateImageFunction interpolator (like in the example https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/Filtering/ResampleImageFilter8.cxx) to get a new image img2 and finally get from img2 all the voxel which have y=0 and z=0 (that’s because I overlapped L1 to axis x which is the locus of points where y=0 and z=0).

6

The ITK Python wrapping can be used to directly interpolate pixel values.

For example,

import itk

image = itk.imread('ImageFile.mhd')
interpolator = itk.LinearInterpolateImageFunction.New(image)
# Interpolate at the physical space location [1.2, 3.4, 4.4]
interpolated_pixel = interpolator.Evaluate([1.2, 3.4, 4.4])
# There is also .EvaluateAtIndex() and .EvaluateAtContinuousIndex() to interpolate with pixel coordinates

Also available in the ITK Python wrapping are Path objects to help define and use paths in space.

3 Likes
7

Per Brad’s suggestions, I will add this as a feature request to SimpleITK.

In the meanwhile you can use a broadcasting approach to obtain the intensity value at an arbitrary set of points:

import SimpleITK as sitk
import numpy as np

img = sitk.ReadImage('training_001_ct.mha')

# generate random samples in image
num_samples = 10

physical_points = []
for pnt in zip(*[list(np.random.random(num_samples)*sz) for sz in img.GetSize()]):
    physical_points.append(img.TransformContinuousIndexToPhysicalPoint(pnt))


interp_grid_img = sitk.Image((num_samples, *([1]*(img.GetDimension()-1))), sitk.sitkUInt8)
displacement_img = sitk.Image((num_samples, *([1]*(img.GetDimension()-1))), sitk.sitkVectorFloat64, img.GetDimension())

for i, pnt in enumerate(physical_points):
    displacement_img[(i, *([0]*(img.GetDimension()-1)))] = np.array(pnt) - np.array(interp_grid_img.TransformIndexToPhysicalPoint((i, *([0]*(img.GetDimension()-1)))))

resampled_points = sitk.Resample(img, interp_grid_img, sitk.DisplacementFieldTransform(displacement_img), sitk.sitkLinear, 0.0, sitk.sitkFloat32)

for i in range(resampled_points.GetWidth()):
    print(resampled_points[i,0,0])
3 Likes