Resample volume to specific voxel spacing - SimpleITK

Beginner Questions
1

Hi everyone!

I have a series of 3D, MR volumes acquired from different scanners. I would like to resample all volumes to their median voxel size [0.39, 0.39, 0.55] mm. Moreover, for each volume, I also have a binary mask (same size or MR) that delineates a ROI.

By looking at similar topics, I came up with a solution that seems to work. However, when I load the newly created volume with itksnap, this is not centered and I can only see part of it.

Can someone please tell me what is wrong? :slightly_smiling_face: Thanks a lot!

Here’s my code:

import numpy as np
import SimpleITK as sitk

def resample_volume(volume_path, mask=False):
    volume = sitk.ReadImage(volume_path)
    volume = sitk.Cast(volume, sitk.sitkFloat32)  # cast to float32
    resample_vol = sitk.ResampleImageFilter()  # create resampler object

    if mask:  # if mask is True
        resample_vol.SetInterpolator = sitk.sitkNearestNeighbor  # set interpolator for mask (it does not create new masks)
    else:
        resample_vol.SetInterpolator = sitk.sitkLinear  # set interpolator for volume

    resample_vol.SetOutputDirection = volume.GetDirection()  # set output volume direction equal to input volume direction
    resample_vol.SetOutputOrigin = volume.GetOrigin()  # set output volume origin equal to input volume origin

    new_spacing = [0.39, 0.39, 0.55]
    resample_vol.SetOutputSpacing(new_spacing)

    orig_size = np.array(volume.GetSize(), dtype=np.int)
    orig_spacing = volume.GetSpacing()
    new_size = orig_size * (np.divide(orig_spacing, new_spacing))
    new_size = np.ceil(new_size).astype(np.int)  # volume dimensions are in integers
    new_size = [int(s) for s in new_size]  # convert from np.array to list
    resample_vol.SetSize(new_size)

    resampled_volume = resample_vol.Execute(volume)
    return resampled_volume

path_to_MR_volume = "/home/newuser/.../mr.nii.gz"
path_to_mask = "/home/newuser/.../mask.nii.gz"
resampled_MR = resample_volume(path_to_MR_volume)  # resample volume
resampled_mask = resample_volume(path_to_mask, mask=True)  # resample mask
2

Hello @Tommaso_Di_Noto,

Didn’t look too deeply into the code as it is rather long for the task. Please find below a shorter version that will resample volumes. As the only difference between resampling a mask and a scalar volume is in the interpolator that is a required parameter and not a boolean flag. Also, not sure why the volumes are converted to sitkFloat32 after reading, but I assume you have your reasons, so combined it with the reading:

def resample_volume(volume_path, interpolator = sitk.sitkLinear, new_spacing = [0.39, 0.39, 0.55]):
    volume = sitk.ReadImage(volume_path, sitk.sitkFloat32) # read and cast to float32
    original_spacing = volume.GetSpacing()
    original_size = volume.GetSize()
    new_size = [int(round(osz*ospc/nspc)) for osz,ospc,nspc in zip(original_size, original_spacing, new_spacing)]
    return sitk.Resample(volume, new_size, sitk.Transform(), interpolator,
                         volume.GetOrigin(), new_spacing, volume.GetDirection(), 0,
                         volume.GetPixelID())

If this code doesn’t do what you wanted, please provide additional details and we’ll get you there.

5 Likes
3

Dear @zivy,

It works! Thank you very much! :smiley:

Just a question: what are the last two arguments for in sitk.Resample (namely, 0 and volume.GetPixelID())?

4

Hi @Tommaso_Di_Noto,

Happy to hear.

Please see the documentation for Resample and ResampleImageFilter. I always have the documentation pages open, there’s just too much to remember.

2 Likes
5

Dear @zivy,

What if I now have to transform (i.e. get the coordinates of) a set of points from original to resampled space?

E.g. [i,j,k] = [120, 150, 90] is an important landmark and I’d like to know its [i,j,k] after resampling.

Thanks a lot again! :slight_smile:

6

It should be something like this:

oldInd=[120, 150, 90]
p=oldImage.TransformIndexToPhysicalPoint(oldInd)
newInd=newImage.TransformPhysicalPointToIndex(p)
print(newInd)
2 Likes
7

Hello @Tommaso_Di_Noto,

To expand upon @dzenanz’s answer, the philosophy of ITK/SimpleITK is that we work in the physical world. That is, there is a physical point associated with the indexes into the image and you move between indexes and physical points using methods from the Image class. The relevant methods are TransformContinuousIndexToPhysicalPoint, TransformIndexToPhysicalPoint, TransformPhysicalPointToContinuousIndex, and TransformPhysicalPointToIndex.

2 Likes
8

Dear @dzenanz, @zivy,

Thank you both for your answers! Of course, it makes sense to leverage the physical space to link points :smiley: