Эффективное моделирование полой сферы в трехмерной воксельной сетке для моделирования компьютерного томографа ⇐ Python

[Anonymous]

Я работаю над проектом, который имитирует среду компьютерного томографа с пластиковыми полыми сферами в воде. Мне нужно смоделировать полую сферу в трехмерной сетке вокселей и вычислить интенсивность каждого воксела на основе его отношения к стенке сферы.

Подробности:
[*]Сетка: сетка вокселей имеет размер 10x21x21, каждый воксель имеет размеры [2, 1,5, 1,5]. [*]Сфера: расположена в нецелочисленных координатах (x0, y0, z0) с определенным внутренним радиусом и толщиной стенок. [*]Расчет интенсивности: если воксель полностью находится внутри стены сферы, его интенсивность равна wall_intensity. Если воксел частично находится внутри стены, его интенсивность равна доле wall_intensity, пропорциональной объему внутри стенки сферы. [*]Цель: расчет должен быть высокоэффективным и точным, поскольку функция вызывается часто.
Проблемы:
[*]Определение того, находится ли воксель полностью или частично внутри стенки сферы, принимая во внимание, что сфера не обязательно находится в центре воксела. [*]Точный расчет значения интенсивности для частично заполненных вокселей.
Вопросы:
[*]Какой алгоритм будет наиболее эффективным для определения отношения вокселя (полностью внутри, частично внутри или снаружи) к стенке сферы? [*]Как точно рассчитать интенсивность частично заполненных вокселей?
Дополнительный контекст:
[*]Wall_intensity представляет собой единицы Хаунсфилда пластика, используемого в сферах. [*]На границе сферы находятся воксели, частично состоящие из воды и частично из пластика, поэтому происходит частичное усреднение по объёму.
Пример процентных значений приведен на рисунке ниже. Обратите внимание, что это 2D-изображение с квадратными пикселями вместо настоящего прямоугольного кубоида.

введите здесь описание изображения

Я попробовал это, но это работает только в том случае, если центр вокселя находится внутри полой сферы.
импортировать numpy как np импортировать matplotlib.pyplot как plt def create_hollow_sphere (grid_shape, Pixel_size, x0, y0, z0, hounsfield_wall, внутренний_радиус, стена_толщина): """ Улучшенная версия функции create_hollow_sphere для правильного расчета эффекта частичного объема. """ сетка = np.zeros(grid_shape) x = np.arange(grid_shape[0]) * размер_пикселя[0] y = np.arange(grid_shape[1]) * размер_пикселя[1] z = np.arange(grid_shape[2]) * размер_пикселя[2] x_grid, y_grid, z_grid = np.meshgrid(x, y, z, indexing='ij') center_x = x0 * размер_пикселя[0] center_y = y0 * размер_пикселя[1] center_z = z0 * размер_пикселя[2] # Вычисляем расстояние каждого воксела от центра расстояния = np.sqrt((x_grid - center_x)**2 + (y_grid - center_y)**2 + (z_grid - center_z)**2) # Определим внешний радиус сферы внешний_радиус = внутренний_радиус + толщина_стены # Рассчитать воксельное включение в полую сферу для меня в диапазоне (grid_shape[0]): для j в диапазоне (grid_shape[1]): для k в диапазоне (grid_shape[2]): # Вычисляем минимальное и максимальное расстояния от центра до углов вокселей min_dist = np.sqrt(max((i * размер пикселя[0] - center_x)**2 - размер_пикселя[0]**2/4, 0) + max((j * размер_пикселя[1] - центр_y)**2 - размер_пикселя[1]**2/4, 0) + max((k * размер_пикселя[2] - центр_z)**2 - размер_пикселя[2]**2/4, 0)) max_dist = np.sqrt((i * размер пикселя[0] - center_x)**2 + размер_пикселя[0]**2/4 + (j * размер_пикселя[1] - центр_y)**2 + размер_пикселя[1]**2/4 + (k * размер_пикселя[2] - центр_z)**2 + размер_пикселя[2]**2/4) # Проверяем, находится ли воксель внутри сферы if (min_dist внутренний_радиус): # Вычисляем долю воксела, заполняемого стенкой сферы пропорция = мин(внешний_радиус, макс_дист) - макс(внутренний_радиус, мин_дист) сетка[i, j, k] = пропорция * hounsfield_wall / (max_dist - min_dist) обратная сетка # Установить параметры сетка_форма = (13, 21, 21) размер пикселя = [2,0, 1,5, 1,5] x0, y0, z0 = 6,5, 10,5, 10,5 # Центр сетки хаунсфилд_стена = 1000 внутренний_радиус = 10 # мм wall_thickness = 1 # мм # Создайте полую сферу с улучшенной функцией сетка = create_hollow_sphere(grid_shape, Pixel_size, x0, y0, z0, hounsfield_wall, внутренний_радиус, толщина_стены) # Постройте срезы с помощью цветных полос рис, топор = plt.subplots(1, 3, figsize=(15, 5)) срезы = [сетка[int(round(x0)), :, :], сетка[:, int(round(y0)), :], сетка[:, :, int(round(z0))]] titles = ['Разрез по x0', 'Разрез по y0', 'Разрез по z0'] для меня в диапазоне (3): im = ax.imshow(slices, cmap='gray') топор.set_title(титулы) fig.colorbar(im, ax=ax) plt.show()