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

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

# 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

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.

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

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?

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

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).

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

For example,

``````import itk

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

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

# 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, *(*(img.GetDimension()-1))), sitk.sitkUInt8)
displacement_img = sitk.Image((num_samples, *(*(img.GetDimension()-1))), sitk.sitkVectorFloat64, img.GetDimension())

for i, pnt in enumerate(physical_points):
displacement_img[(i, *(*(img.GetDimension()-1)))] = np.array(pnt) - np.array(interp_grid_img.TransformIndexToPhysicalPoint((i, *(*(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