Efficiently Modeling a Hollow Sphere in a 3D Voxel Grid for CT Scanner Simulation

Question

I am working on a project that simulates a CT scanner environment with plastic hollow spheres in water. I need to model a hollow sphere in a 3D grid of voxels and calculate the intensity of each voxel based on its relation to the sphere's wall.

Details:

• Grid: The voxel grid is 10x21x21, with each voxel having dimensions [2, 1.5, 1.5].
• Sphere: Placed at non-integer coordinates (x0, y0, z0), with a specific inner radius and wall thickness.
• Intensity Calculation: If a voxel is entirely within the sphere's wall, its intensity equals wall_intensity. If a voxel is partially within the wall, its intensity is a fraction of wall_intensity, proportional to the volume inside the sphere's wall.
• Objective: The calculation needs to be highly efficient and accurate, as the function is called frequently.

Challenges:

• Determining whether a voxel is completely or partially within the sphere's wall, considering the sphere is not necessarily centered within a voxel.
• Accurately calculating the intensity value for partially filled voxels.

Questions:

• What would be the most efficient algorithm to determine the voxel's relation (fully inside, partially inside, or outside) to the sphere's wall?
• How can I accurately calculate the intensity for partially filled voxels?

Additional Context:

• The wall_intensity represents the Hounsfield units of the plastic used in the spheres.
• At the border of the sphere there are voxels that are part water and part plastic, so a partial volume averaging takes places.

See the figure below for an example of the percentages. Note that this is in 2d and with square pixels instead of the actual rectangular cuboid.

enter image description here

This is what I tried, but it only works if the center of the voxel is inside of the hollow sphere.

import numpy as np
import matplotlib.pyplot as plt

def create_hollow_sphere(grid_shape, pixel_size, x0, y0, z0, hounsfield_wall, inner_radius, wall_thickness):
    """
    Improved version of the create_hollow_sphere function to correctly calculate the partial volume effect.
    """
    grid = np.zeros(grid_shape)
    x = np.arange(grid_shape[0]) * pixel_size[0]
    y = np.arange(grid_shape[1]) * pixel_size[1]
    z = np.arange(grid_shape[2]) * pixel_size[2]
    x_grid, y_grid, z_grid = np.meshgrid(x, y, z, indexing='ij')

    center_x = x0 * pixel_size[0]
    center_y = y0 * pixel_size[1]
    center_z = z0 * pixel_size[2]

    # Calculate the distance of each voxel from the center
    distances = np.sqrt((x_grid - center_x)**2 + (y_grid - center_y)**2 + (z_grid - center_z)**2)

    # Define the outer radius of the sphere
    outer_radius = inner_radius + wall_thickness

    # Calculate voxel inclusion in the hollow sphere
    for i in range(grid_shape[0]):
        for j in range(grid_shape[1]):
            for k in range(grid_shape[2]):
                # Calculate the minimum and maximum distances from the center to the voxel corners
                min_dist = np.sqrt(max((i * pixel_size[0] - center_x)**2 - pixel_size[0]**2/4, 0) +
                                   max((j * pixel_size[1] - center_y)**2 - pixel_size[1]**2/4, 0) +
                                   max((k * pixel_size[2] - center_z)**2 - pixel_size[2]**2/4, 0))
                max_dist = np.sqrt((i * pixel_size[0] - center_x)**2 + pixel_size[0]**2/4 +
                                   (j * pixel_size[1] - center_y)**2 + pixel_size[1]**2/4 +
                                   (k * pixel_size[2] - center_z)**2 + pixel_size[2]**2/4)

                # Check if the voxel is within the wall of the sphere
                if (min_dist < outer_radius and max_dist > inner_radius):
                    # Calculate the proportion of the voxel that is filled by the sphere's wall
                    proportion = min(outer_radius, max_dist) - max(inner_radius, min_dist)
                    grid[i, j, k] = proportion * hounsfield_wall / (max_dist - min_dist)

    return grid

# Set parameters
grid_shape = (13, 21, 21)
pixel_size = [2.0, 1.5, 1.5]
x0, y0, z0 = 6.5, 10.5, 10.5  # Center of the grid
hounsfield_wall = 1000
inner_radius = 10  # mm
wall_thickness = 1  # mm

# Create the hollow sphere with the improved function
grid = create_hollow_sphere(grid_shape, pixel_size, x0, y0, z0, hounsfield_wall, inner_radius, wall_thickness)

# Plot the slices with colorbars
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
slices = [grid[int(round(x0)), :, :], grid[:, int(round(y0)), :], grid[:, :, int(round(z0))]]
titles = ['Slice at x0', 'Slice at y0', 'Slice at z0']
for i in range(3):
    im = ax[i].imshow(slices[i], cmap='gray')
    ax[i].set_title(titles[i])
    fig.colorbar(im, ax=ax[i])

plt.show()

Answer 1

Ok, here is just simple idea how to do it. Use sampling a la monte Carlo to compute avg density per voxel. Code is completely untested, I'm from computer, and it's a New Year Eve

import numpy as np

rng = np.random.default_rng(135797537)

dx = dy = 1.5
dz = 2.0

nx = ny = 21
nz = 10

N = nx*ny*nz*1000 # nof samples

grid = np.zeros((nx, ny, nz), dtype=np.float32);

x = rng.uniform(low=0.0, high=nx*dx, size=N, dtype=np.float32)
y = rng.uniform(low=0.0, high=ny*dy, size=N, dtype=np.float32)
z = rng.uniform(low=0.0, high=nz*dz, size=N, dtype=np.float32)

r2 = np.square(x-x0) + np.square(y-y0) + np.square(z-z0)

mask = r2 <= r0*r0
grid[mask] += density0
grid[~mask & (r2 < R0*R0)] += density1

grid /= np.float32(N)

UPDATE

Actually points to be sampled not in whole phantom but in (x0-R0,x0+R0,y0-R0,yo+R0,z0-R0,z0+R0) subcube