3d体素数据斜向切片和曲面数据重建

1根据一串点的坐标生成贝塞尔曲线

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline


x = np.array([48.096, 77.488, 89.512, 100.2, 112.224, 128.256, 146.96, 164.328, 181.696,
              197.728, 216.432, 233.8, 245.824, 256.512, 264.528, 297.928])
y = np.array([203.072, 173.68, 136.272, 105.544, 81.496, 62.792, 52.104, 50.768, 52.104,
              53.44, 64.128, 82.832, 108.216, 137.608, 172.344, 208.416])

spl = make_interp_spline(x, y, k=2)

# 定义用于绘制曲线的参数
param = np.linspace(min(x), max(x), 500)
splined_x = param
splined_y = spl(param)

# 绘制曲线和等距点
plt.figure(figsize=(10, 6))
plt.scatter(x, y, label='Original Points')
plt.plot(splined_x, splined_y, label='bezier_')
plt.legend()
plt.show()

2斜向切片算法

计算斜向切片区域的索引，使用插值法索引

切片时超出范围的索引可以设置填充方式，现为使用0填充，计算切片索引的参数为：

水平面上的坐标 ：x， y 轴所在平面的坐标

角度： 切面与y轴的弧度值

宽和高：宽的范围为(-w/2+ x , w/2 + x) 高为（0， h）

import numpy as np
from scipy.ndimage import map_coordinates


def get_slice_position(point, angle, w, h):
    # 中心点坐标
    slice_position = np.repeat(point[:, np.newaxis], w, axis=1)
    # 生成中心位于原点的线段，并获得旋转后的坐标
    arr = np.zeros((2, w))
    arr[1, :] = arr[1, :] + np.arange(-w/2, w/2)
    arr[0, :], arr[1, :] = np.sin(angle) * arr[1, :], np.cos(angle) * arr[1, :]
    # 得到裁剪线段旋后的坐标
    slice_position = slice_position + arr
    # 线段坐标扩充成切面坐标
    slice_position = np.repeat(slice_position[:, :, np.newaxis], h, axis=2)
    h_ = np.repeat(np.arange(0, h, )[::-1].reshape((1, 1, h)), w, axis=1)
    slice_position = np.concatenate((slice_position, h_), axis=0)
    return slice_position



volume_data = np.zeros((400, 400, 400))

# 定义切片的角度和位置
angle = np.pi / 4  # 切片的角度（以弧度表示）
slice_position = np.array([89.512, 136.272])  # 切片位置的坐标 [x, y]

# 获取斜向的切片索引
grid = get_slice_position(slice_position, angle, 100, 300)

# 利用插值方法获取切片数据
slice_data = map_coordinates(volume_data, grid, order=1, cval=0)

3曲面数据重建

效果举例：画出一条长度为400的曲线，然后定下厚度100和高度200，得到一个400*100*200的新体素数据

实现过程：

1创建一条贝塞尔曲线

2获取曲线上的等距离点，距离为1

3获取每个点的与该点斜率垂直的斜面的切片索引（使用和斜率垂直的切面可能不够精确，可考虑其他的角度给出方式）

4使用得到的索引插值索引得到重建的数据

import numpy as np
from scipy.interpolate import make_interp_spline
from scipy.ndimage import map_coordinates


def get_derivative_at_point(spline, point):
    # 生成样条的一阶导数
    spline_derivative = spline.derivative()
    # 计算指定点的导数
    return spline_derivative(point)


# 获得曲线的某个位置的斜率，参数为x坐标
def get_perpendicular_angle_with_y_axis(spline, point):
    # 获取该点的斜率
    slope = get_derivative_at_point(spline, point)
    angle_with_x_axis = np.arctan(slope)
    return -angle_with_x_axis


# 获得曲线上间距为1的点的坐标和斜率
def get_equidistant_points(spl, x, distance=1):
    # 对曲线进行细分
    num_points = 200
    fine_x = np.linspace(min(x), max(x), num_points)
    fine_y = spl(fine_x)

    # 计算曲线上每一点的累积长度
    distances = np.sqrt(np.diff(fine_x)**2 + np.diff(fine_y)**2)
    cumulative_length = np.cumsum(distances)
    cumulative_length = np.insert(cumulative_length, 0, 0)  # 在开始处插入0

    # 插值等距点
    equidistant_points = []
    current_length = distance
    while current_length < cumulative_length[-1]:
        # 找到累积长度中最接近当前长度的点
        index = np.searchsorted(cumulative_length, current_length)
        if index >= len(fine_x):
            break

        # 插值找到准确位置
        lambda_ = (current_length - cumulative_length[index-1]) / (cumulative_length[index] - cumulative_length[index-1])
        new_x = fine_x[index-1] + lambda_ * (fine_x[index] - fine_x[index-1])
        new_y = fine_y[index-1] + lambda_ * (fine_y[index] - fine_y[index-1])
        angle = get_perpendicular_angle_with_y_axis(spl, new_x)
        equidistant_points.append((new_x, new_y, angle))
        current_length += distance
    return equidistant_points


# 获得某个点的位置的斜向切片索引
def get_slice_position(points, w, h):
    # 中心点坐标
    angle = points[:, 2]
    points = points[:, :2]
    slice_position = np.repeat(points[:, :, np.newaxis], w, axis=2)
    slice_position = slice_position.transpose((1, 0, 2))
    # 生成中心位于原点的线段，并获得旋转后的坐标
    arr = np.zeros((2, len(points), w))
    arr[1] = arr[1] + np.arange(-w/2, w/2)
    arr[0] = np.sin(angle)[:, None] * arr[1]
    arr[1] = np.cos(angle)[:, None] * arr[1]
    # 得到裁剪线段旋后的坐标
    slice_position = slice_position + arr
    # 线段坐标扩充成切面坐标
    slice_position = np.repeat(slice_position[:, :, :, np.newaxis], h, axis=3)
    h_ = np.repeat(np.arange(0, h)[::-1].reshape((1, 1, h)), w, axis=1)
    h_ = np.repeat(h_[:, np.newaxis, :, :], len(points), axis=1)
    slice_position = np.concatenate((slice_position, h_), axis=0)
    return slice_position


# 定义已知点
x = np.array([48.096, 77.488, 89.512, 100.2, 112.224, 128.256, 146.96, 164.328, 181.696,
              197.728, 216.432, 233.8, 245.824, 256.512, 264.528, 297.928])
y = np.array([203.072, 173.68, 136.272, 105.544, 81.496, 62.792, 52.104, 50.768, 52.104,
              53.44, 64.128, 82.832, 108.216, 137.608, 172.344, 208.416])

# 创建插值曲线
spl = make_interp_spline(x, y, k=2)
# 获得曲线上等距离的坐标点, 默认相邻距离为1
equidistant_points = np.array(get_equidistant_points(spl, x))
# 获得曲面区域的索引
grid = get_slice_position(equidistant_points, 100, 300)
# 创建体素数据
volume_data = np.zeros((400, 400, 400))
# 使用插值法索引得到该区域
slice_data = map_coordinates(volume_data, grid, order=1, cval=0)
print(slice_data)
print(slice_data.shape)