Placing a ball in a CT volume

Beginner Questions
1

I’m trying to place a 25mm reference ball within a CT. I think I have a function that can create a ball in a volume that has the same shape and spacing as the CT. However, I tried adding, pasting and masking but the ball never appears in the CT correctly. Currently the ball gets placed in the corner and everything else is black. I have a feeling that setting the origin of the ball volume to the same origin as the CT did not really work.

Ideally I would like to do all that in ITK but I could also convert back from SimpleITK. I can also imagine that I don’t need to create a separate ball at all, since I could directly draw in the CT. But my attempts to do that failed completely.

Could someone give me a hint? Here is the code I’ve used to create the ball:

def create_ref_ball(pos, radius, spacing, size, fg_value=255, bg_value=0):
“”"Create a reference ball.

Parameters
----------
pos : list
    Position of the ball.
radius : float
    Radius of the ball.
spacing : list
    Spacing of the image.
size : list
    Size of the image.

Returns
-------
SimpleITK.Image
    Reference ball.
"""
# create a spherical gaussian blob
gaussian = sitk.GaussianSource(
    sitk.sitkUInt8,
    size=size,
    sigma=[radius, radius, radius],
    mean=pos,
    spacing=spacing,
)

# threshold to create a binary ball
ball = sitk.BinaryThreshold(gaussian, 150.0, 255.0, fg_value, bg_value)
# itk_ball = ConvertSimpleItkImageToItkImage(ball, itk.UC)
return ball

if name == “main”:
“”“25mm ball at [20, 20, 20] with spacing [0.4, 0.4, 1.0] in a empty
[100, 100, 100] volume”“”

ball = create_ref_ball([20, 30, 40], 12.5, [0.4, 0.4, 1.0], [100, 100, 100])
# itk.imwrite(ball, "ball.nii.gz")
sitk.WriteImage(ball, "ball.nii.gz")
2

You could try using AddImageFilte, or PasteImageFilter, or MaximumImageFilter to combine you ball.nii.gz with the original CT image.

To properly combine using add and max, your ball image should be first resampled to the grid of your original CT image (using identity transform). Otherwise results might be unexpected, or unintuitive.

1 Like
3

Hello,

For the GaussianSource, the CT’s meta-data should be copied to be used as a reference image for the output image. The origin, direction, spacing, and size should all be copied. This will ensure that the ball image and the CT are congruent in physical space.

In SimpleITK’s nightly builds, just merged was adding the ability to use a mask image in the index operator for assignments. So the following works in the SimpleITK 2.3.dev:

ct_image[ball_mask] = 0

This will enforce that the two images have the same meta-data and are congruent.

2 Likes
4

Hello @User0,

Naive code implementation below, will cause an exception if the sphere is partially outside the image. You can do better, but as long as the number of voxels occupied by the sphere isn’t too large and it is inside the image the naive implementation isn’t too horrible:

import SimpleITK as sitk
import numpy as np

def naive_draw_sphere(image, radius, location, value=-1000):
    rad2 = radius**2
    image_spacing = image.GetSpacing()
    min_coords, max_coords = zip(*[(loc-radius-spc, loc+radius+spc)for loc,spc in zip(location, image_spacing)])
    for x in np.arange(min_coords[0], max_coords[0], image_spacing[0]):
      for y in np.arange(min_coords[1], max_coords[1], image_spacing[1]):
          for z in np.arange(min_coords[2], max_coords[2], image_spacing[2]):
              if((x-location[0])**2 + (y-location[1])**2 + (z-location[2])**2)<=rad2:
                  image[image.TransformPhysicalPointToIndex([x,y,z])] = value

dcm_dir = "data/"

image = sitk.ReadImage(sitk.ImageSeriesReader.GetGDCMSeriesFileNames(dcm_dir))
image_center = image.TransformIndexToPhysicalPoint([sz//2 for sz in image.GetSize()])

#place sphere in center of image
sphere_center_location = image_center
sphere_radius = 25

sitk.Show(image, "before")
naive_draw_sphere(image, sphere_radius, sphere_center_location)
sitk.Show(image, "after")
3 Likes
5

Oh wow thank you all for those quick answers. This community is really great. I was able to solve the issue with zivy’s suggestion. I tried to install the dev 2.3dev version suggested by blowekamp too, but I couldn’t find a release on GitHub or a branch. Must have looked in the wrong place.

Anyway, thank you all very much.

2 Likes