在 3D 体素网格中高效建模空心球以进行 CT 扫描仪模拟

Question

我正在开展一个项目，该项目用水中的塑料空心球来模拟 CT 扫描仪环境。我需要在 3D 体素网格中对空心球体进行建模，并根据每个体素与球体壁的关系计算每个体素的强度。

详情：

• 网格：体素网格为 10x21x21，每个体素的尺寸为 [2, 1.5, 1.5]。
• 球体：放置在非整数坐标（x0，y0，z0）处，具有特定的内半径和壁厚。
• 强度计算：如果体素完全位于球体壁内，则其强度等于 wall_intensity。如果体素部分位于球体壁内，则其强度是 wall_intensity 的一部分，与球体壁内的体积成正比。
• 目的：由于函数调用频繁，计算需要高效、准确。

挑战：

• 考虑到球体不一定位于体素中心，确定体素是否完全或部分位于球体壁内。
• 准确计算部分填充体素的强度值。

问题：

• 确定体素与球体壁的关系（完全内部、部分内部或外部）的最有效算法是什么？
• 如何准确计算部分填充体素的强度？

其他背景：

• wall_intensity 代表球体中使用的塑料的 Hounsfield 单位。
• 在球体的边界处，存在部分是水、部分是塑料的体素，因此会发生部分体积平均。

请参阅下图了解百分比示例。请注意，这是二维的，并且具有正方形像素，而不是实际的长方体。

这是我尝试过的，但只有当体素的中心位于空心球体内部时它才有效。

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

好的，这只是如何做到这一点的简单想法。使用蒙特卡洛采样来计算每个体素的平均密度。代码完全未经测试，我来自计算机，而且现在是除夕夜

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)