Isotropic resampling of a CT scan

Beginner Questions
1

Hi,

I am working on resampling and interpolating some CT scans to have isotropic voxel spacing ([1, 1, 1]). I have tried two approaches that both execute without error, but I am not sure why the first approach looks odd when viewed in Slicer.

The first approach uses the original CT as a reference image but overwrites to isotropic spacing. This approach keeps the same size (e.g. 512x512x134).

While the second approach derives a new size based on the desired spacing, such that the resampled image has the same volume in physical space (mm).

import SimpleITK as sitk

file = "example.nrrd"
img = sitk.ReadImage(file)

## First Approach (keep the same size)
rs = sitk.ResampleImageFilter()
rs.SetReferenceImage(img)
rs.SetOutputSpacing([1, 1, 1])
img2 = rs.Execute(img)
sitk.WriteImage(img2, "imgSameSize.nrrd")

## Second Approach (keep the same volume)
rs2 = sitk.ResampleImageFilter()

new_x = round(img.GetSize()[0] * img.GetSpacing()[0] / 1)
new_y = round(img.GetSize()[1] * img.GetSpacing()[1] / 1)
new_z = round(img.GetSize()[2] * img.GetSpacing()[2] / 1)

rs2.SetOutputSpacing([1, 1, 1])
rs2.SetSize([new_x, new_y, new_z])
rs2.SetOutputDirection(img.GetDirection())
rs2.SetOutputOrigin(img.GetOrigin())
rs2.SetTransform(sitk.Transform())
rs2.SetInterpolator(sitk.sitkLinear)

img3 = rs.Execute(img)

sitk.WriteImage(img3, "imgDiffSize.nrrd")

I’m hoping that someone here could offer insight into whether I have made any mistakes in either approach? I come from the R world so I apologize if my code is not adhering to python best practices.

2

The first approach creates a bigger image, creating gray bars (default value 0) beyond the physical extent of the original image. To make it look less weird, you could specify default value of -1000 in order to get black bars around. The second approach seems correct.

3

Thanks for the response, @dzenanz.

Adding rs.SetDefaultPixelValue(-1000) to the first approach didn’t seem to fix the “gray bar” issue; I’m not sure why. Also, the coronal and sagittal reconstructions are cropped, which makes me think there are other issues with this first approach. Any idea why this might occur?

Screen Shot 2021-05-26 at 11.24.49 AM
Screen Shot 2021-05-26 at 11.24.49 AM1325Ă—1078 121 KB

4

Actually, setting the default pixel value to -3024 did work to “remove” the gray bars. The cropped reconstructions issue remains.

5

Your original spacing along k/Z axis is larger than 1.0, so you get cropping along that axis instead of extra empty space.

1 Like
6

Hello @mattwarkentin,

The first approach is telling the resmapler to use the original image’s size and the new spacing. As a result the resampling grid does not cover the same physical extent as the original image. You need to compute a new output size when you have a new output spacing in order to cover the same physical region.

Take a look at the make_isotropic function in this Jupyter notebook.

3 Likes
7

Thanks both of you for the responses. This is very helpful.

You need to compute a new output size when you have a new output spacing in order to cover the same physical region.

So it seems like I have done this correctly in the second approach. It looks like my code (albeit less concise) is doing the same thing as the make_isotropic() function.

2 Likes
8

When we resample images, we usually determine the fill color (that will be passed to the image resampler to fill the empty regions) as the median value of the 8 corner values of the volume.

I guess you need this resampling to normalize your input data for deep learning. If this is the case then I would recommend to have a look at torchio and monai, as they may already have robust, well-tested, and optimized implementation of all the data normalization and augmentation steps you plan to perform. Torchio uses ITK internally, maybe monai, too.

3 Likes
9

Thanks for bringing these projects to my attention, @lassoan. At a quick glance they both look like they could be extremely helpful.

1 Like
10

My two cents, in five lines:

import torchio as tio

image = tio.ScalarImage("example.nrrd")
resample = tio.Resample()  # default is 1 mm isotropic
resampled = resample(image)
resampled.save("imgDiffSize.nrrd")

As @lassoan said, TorchIO uses SimpleITK under the hood for this type of operations.

2 Likes
11

Thanks for the input, @fepegar. I’ve actually spent most of today learning and playing around with torchio and I have found it extremely user-friendly and powerful. I am very glad @lassoan brought it to my attention.

3 Likes