volume regularization using ITK

Algorithms
, ,
1

Hello,

I’m working on a Python application that uses ITK, and I’m looking for guidance on how to regularize or resample a 3D volume when the slice spacing is not uniform—i.e., the distance between slices varies throughout the volume.

I know that 3D Slicer handles this kind of situation internally when loading and resampling volumes, but I need to implement a similar solution within my own application using ITK and Python.

Is there a recommended approach or existing method in ITK (or SimpleITK) for dealing with this? Ideally, I’d like to resample the volume to have isotropic or at least regular spacing across all dimensions.

Any suggestions, code examples, or references would be greatly appreciated.

Thanks in advance for your help!

Best regards,

2

A very generic approach would be to treat this problem as a scattered data approximation of a scalar field. For that, you could use this B-spline ITK filter. Unfortunately, I don’t see it implemented in SimpleITK. However, we do have it implemented in ANTsPy as the function ants.fit_bspline_object_to_scattered_data. Input is a numpy point set providing the data (in your case, this would be the intensity values), a numpy point set providing the physical location of those intensity values, and physical space definition of the bounding box (in terms of origin, spacing, direction). Output is an ITK image representing the sampled B-spline object (i.e., volume).

Here’s an example of this approach where I had 3 sets of ~10 orthogonal thick slices oriented in 3-D space sampling a particular . I wanted to use all three sets of data to create a single image volume.

NIck

2 Likes
3

Perhaps I didn’t explain my issue clearly.
I have a collection of DICOM images (a series, in DICOM terminology) with varying inter-slice spacing. I need to generate an isotropic ITK volume from this data. I’m wondering if there is a standard or recommended method to achieve this.

Thanks

4

Hello @Eva_Monclus_Lahoya,

Nick was pointing you to a generic scattered data implementation that can deal with your task and much more.

Below is code that does what you want but is much more limited due to assumptions on scan direction and usage of linear interpolation. There are probably more issues as I haven’t thought of corner cases. I suspect it will be easier for you to understand this code. The data referenced in the code is available here.

import SimpleITK as sitk
import random
import numpy as np

def get_uniformly_z_spaced_volume(file_names):
    """
    file_names - list which is sorted according to z scan direction.
    series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dir_name)
    file_names = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(dir_name, series_ids[0])
    """
    reader = sitk.ImageFileReader()
    reader.SetImageIO("GDCMImageIO")

    # Code assumes orthonormal axes and z scan direction
    original_z = []
    # Get the z coordinate for all images
    for file_name in file_names:
      reader.SetFileName(file_name)
      reader.ReadImageInformation()
      original_z.append(reader.GetOrigin()[2])

    # find the minimal z spacing in the possibly non-uniformly spaced dataset
    # and use that as the new uniform spacing
    new_z_spacing = np.sort((np.array(original_z[1:]) - np.array(original_z[0:-1])))[0]

    new_images = [sitk.ReadImage(file_names[0])]
    output_size = new_images[-1].GetSize()
    output_direction = new_images[-1].GetDirection()

    # converting from tuple to list because these are changed later on
    output_spacing = list(new_images[-1].GetSpacing())
    output_origin = list(new_images[-1].GetOrigin())
    origin_z = output_origin[2]
    # Extract as 2D slice
    new_images[-1] = new_images[-1][...,0]

    # read consecutive pairs of images and use this mini-volume to linearly
    # interpolate slices
    i = 0
    curr_vol = sitk.ReadImage(file_names[i:i+2])
    for j, curr_z in enumerate(np.arange(original_z[0]+new_z_spacing, original_z[-1]+0.5*new_z_spacing, new_z_spacing)):
      if curr_z > original_z[i]+curr_vol.GetSpacing()[2]:
          i+=1
          curr_vol = sitk.ReadImage(file_names[i:i+2])
      output_origin[2] = curr_z
      new_images.append(sitk.Resample(curr_vol,
                                      output_size,
                                      sitk.Transform(),
                                      sitk.sitkLinear,
                                      output_origin,
                                      output_spacing,
                                      output_direction)[...,0])
    image = sitk.JoinSeries(new_images)
    output_origin[2] = origin_z
    output_spacing[2] = new_z_spacing
    image.SetOrigin(output_origin)
    image.SetSpacing(output_spacing)
    return image

dir_name = "CIRS057A_MR_CT_DICOM"
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dir_name)
original_file_names = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(dir_name, series_ids[0])

# Create a non-uniformly spaced dataset from a uniformly spaced one.
num_images_to_remove = 20
image_indexes_to_remove = random.sample(range(0,len(original_file_names)), num_images_to_remove)



file_names = [original_file_names[i] for i in range(len(original_file_names)) if i not in image_indexes_to_remove]

sitk.Show(get_uniformly_z_spaced_volume(file_names), "res")
sitk.Show(sitk.ReadImage(original_file_names), "orig")
2 Likes