通过CBCT投影得到的MIP图像拟合牙弓曲线

1获得牙齿曲线

通过3d体素数据投影得到图像和高亮区域蒙版的代码：

import numpy as np
import cv2

# volume_data为np格式的三维数据，选取一个高度范围进行投影
img_bg = volume_data[:, :, 90:205]
img_bg[img_bg < 0] = 0
img_bg = img_bg.max(axis=2)
img_bg = (img_bg + abs(img_bg.min())) / (img_bg.max() - img_bg.min()) * 255
img_bg = img_bg.astype(np.uint8)
img_bg = img_bg.transpose((1, 0))

# 阈值为150，获得蒙版

_, teeth_mask_bg = cv2.threshold(image_bg, 150, 255, cv2.THRESH_BINARY)

得到的投影和牙齿区域蒙版：

拟合牙齿曲线的代码，实现思路为以图片中心为中心点，向切蛋糕一样将牙齿蒙版切成小块，计算每个小块的中心点，用这一串中心点拟合一条曲线：

import cv2
import numpy as np
from PIL import Image
from matplotlib import pyplot as plt

image_path_bg = r'D:\Tooth-cls_feature\img_bg.png'
image_bg = cv2.imread(image_path_bg, 0)


def find_middle_point(mask, center, angle_low, angle_high):
    """
    Find the middle point of the teeth in a specified angular section of the mask.
    """
    # Extract the region of interest
    coords = np.column_stack(np.where(mask > 0))
    angles = np.arctan2(coords[:, 0] - center[1], coords[:, 1] - center[0]) % (2 * np.pi)
    # Filter points that are within the specified angular range
    index = (angles >= angle_low) & (angles < angle_high)
    if not np.any(index):
        return None
    section_points = coords[index]
    # Find the middle point by averaging the coordinates
    middle_point = np.mean(section_points, axis=0)
    return middle_point


_, teeth_mask_bg = cv2.threshold(image_bg, 150, 255, cv2.THRESH_BINARY)
num_bins_bg = 100


image_center_bg = (int(teeth_mask_bg.shape[1]/2), int(teeth_mask_bg.shape[0]/2))
middle_points_bg = []

angular_width_bg = 2 * np.pi / num_bins_bg

for i in range(num_bins_bg):
    angle_low_bg = i * angular_width_bg
    angle_high_bg = (i + 1) * angular_width_bg
    middle_point_bg = find_middle_point(teeth_mask_bg, image_center_bg, angle_low_bg, angle_high_bg)
    if middle_point_bg is not None:
        middle_points_bg.append(middle_point_bg)

middle_points_bg = np.array(middle_points_bg)
p1 = np.polyfit(middle_points_bg[:, 1], middle_points_bg[:, 0], 4)
plt.figure(figsize=(10, 6))
x_range = np.linspace(50, 275, 500)
y_range = np.polyval(p1, x_range)
img_bg = Image.fromarray(image_bg)
plt.imshow(img_bg)
plt.scatter(middle_points_bg[:, 1], middle_points_bg[:, 0], color='red', label='Known Points')  # 绘制已知点
plt.plot(x_range, y_range, color='blue', label='Interpolation Curve')  # 绘制插值曲线
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Quadratic Polynomial Interpolation')
plt.legend()
plt.grid(True)
plt.show()

效果如图：

2获得颌骨曲线

实现思路：取得颌骨的左下角和右下角坐标和牙齿的左下角和右下角坐标，用这四个点拟合一条曲线；

这段代码紧跟拟合牙齿曲线的代码：

_, bone_mask_bg = cv2.threshold(image_bg, 20, 255, cv2.THRESH_BINARY)

contours, _ = cv2.findContours(bone_mask_bg, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

# 寻找最大轮廓
max_contour = max(contours, key=cv2.contourArea)

x, y, w, h = cv2.boundingRect(max_contour)

# 左下角坐标 (x, y + h)
left_bottom = (x, y + h)

# 右下角坐标 (x + w, y + h)
right_bottom = (x + w, y + h)
bone_points_x = np.array([left_bottom[0], middle_points_bg[:, 1].min(), middle_points_bg[:, 1].max(), right_bottom[0]])
bone_points_y = np.array([left_bottom[1], middle_points_bg[:, 0].max(), middle_points_bg[:, 0].max(), right_bottom[1]])
p2 = np.polyfit(bone_points_x, bone_points_y, 2)

x_range1 = np.linspace(0, 334, 500)
y_range1 = np.polyval(p2, x_range1)
img_bg = Image.fromarray(image_bg)
plt.imshow(img_bg)
plt.scatter(bone_points_x, bone_points_y, color='blue', label='Known Points')  # 绘制已知点
plt.plot(x_range1, y_range1, color='red', label='Interpolation Curve')  # 绘制插值曲线
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Quadratic Polynomial Interpolation')
plt.legend()
plt.grid(True)
plt.show()

颌骨曲线效果如图：

3两条曲线交叉区域平滑连接

先看看现在的效果：

和需要达到的效果：

实现原理：取交点左右固定一段距离的点，计算出终点x2处两条曲线的差值h，将x1和x2之间的这段曲线加上一个和h有关的权重，使其平滑连接至两一条曲线的x2处，这里采用的是正弦函数结合一个指数函数。

代码如下：

# 余弦加权函数
def cosine_weight(x, x_start, x_end):
    return 0.1 * np.sin(np.pi * (x - x_start) / (x_end - x_start))


def add_weight(x, x_start, x_end, y_d, L=True):
    if L:
        return ((x - x_start) / (x_end - x_start)) ** 3 * y_d
    else:
        return (1 - (x - x_start) / (x_end - x_start)) ** 3 * y_d



def curves_smooth(cross_x, x_range_bone, y_range_tooth, y_range_bone, p1, p2):
    """
    将两条抛物线交叉处平滑连接
    cross_x: 交叉点x坐标
    cross_y: 交叉点y坐标
    x_range_bone:抛物线的x坐标，两条抛物线的x坐标相同
    y_range_tooth：牙齿曲线的y坐标
    y_range_bone：颌骨曲线的y坐标
    p1: 牙齿曲线多项式的参数
    p2: 颌骨曲线多项式的参数
    """
    # 假设交点的x坐标大约在x_intersect范围内，我们需要确定这个范围
    x_intersect_start = cross_x[0] - 15  # 这个值需要您根据实际情况调整
    x_intersect_end = cross_x[0] + 15  # 这个值需要您根据实际情况调整

    # 计算在交点附近的权重
    weights = cosine_weight(x_range_bone, x_intersect_start, x_intersect_end)
    y_intersect_end_tooth = np.polyval(p1, x_intersect_end)
    y_intersect_end_bone = np.polyval(p2, x_intersect_end)
    y_d = -abs(y_intersect_end_bone - y_intersect_end_tooth)
    weights_add = add_weight(x_range_bone, x_intersect_start, x_intersect_end, y_d)
    y_range_weighted = y_d * weights + weights_add

    # 合并两条曲线
    y_combined = y_range_bone + y_range_weighted

    # 确保在交点范围外的曲线保持不变
    y_combined[x_range_bone < x_intersect_start] = y_range_bone[x_range_bone < x_intersect_start]
    y_combined[x_range_bone > x_intersect_end] = y_range_tooth[x_range_bone > x_intersect_end]

    # 处理第二个交点
    x_intersect_start = cross_x[1] - 15  
    x_intersect_end = cross_x[1] + 15  

    # 计算在交点附近的权重
    weights = cosine_weight(x_range_bone, x_intersect_start, x_intersect_end)
    y_intersect_start_tooth = np.polyval(p1, x_intersect_start)
    y_intersect_start_bone = np.polyval(p2, x_intersect_start)
    y_d = -abs(y_intersect_start_bone - y_intersect_start_tooth)
    weights_add = add_weight(x_range_bone, x_intersect_start, x_intersect_end, y_d, L=False)
    y_range_weighted = y_d * weights + weights_add

    # 合并两条曲线
    y_combined1 = y_range_bone + y_range_weighted

    # 确保在交点范围外的曲线保持不变
    y_combined1[x_range_bone < x_intersect_start] = y_combined[x_range_bone < x_intersect_start]
    y_combined1[x_range_bone > x_intersect_end] = y_range_bone[x_range_bone > x_intersect_end]
    return y_combined1

附上完整代码：

from PIL import Image
import cv2
import numpy as np
from matplotlib import pyplot as plt
from numpy.polynomial.polynomial import Polynomial


# 余弦加权函数
def cosine_weight(x, x_start, x_end):
    return 0.1 * np.sin(np.pi * (x - x_start) / (x_end - x_start))


def add_weight(x, x_start, x_end, y_d, L=True):
    if L:
        return ((x - x_start) / (x_end - x_start)) ** 3 * y_d
    else:
        return (1 - (x - x_start) / (x_end - x_start)) ** 3 * y_d


def find_middle_point(mask, center, angle_low, angle_high):
    """
    Find the middle point of the teeth in a specified angular section of the mask.
    """
    # Extract the region of interest
    coords = np.column_stack(np.where(mask > 0))
    angles = np.arctan2(coords[:, 0] - center[1], coords[:, 1] - center[0]) % (2 * np.pi)
    # Filter points that are within the specified angular range
    index = (angles >= angle_low) & (angles < angle_high)
    if not np.any(index):
        return None
    section_points = coords[index]
    # Find the middle point by averaging the coordinates
    middle_point = np.mean(section_points, axis=0)
    return middle_point


def find_cross_point(p1, p2):
    """
    参数为两条曲线的参数
    return:交点的x坐标和交点的y坐标
    """
    # 创建两个多项式对象
    poly1 = Polynomial(p1[::-1])  # np.polyfit的系数是从最高次开始的，Polynomial需要从常数项开始
    poly2 = Polynomial(p2[::-1])

    # 计算两个多项式的差
    poly_diff = poly1 - poly2

    # 寻找差多项式的根，即两曲线的交点的x坐标
    roots = poly_diff.roots()

    # 因为根可能包括复数，我们只想要实数根，因此需要过滤
    real_roots = roots[np.isreal(roots)].real

    # 计算交点的y坐标
    intersections_y = poly1(real_roots)
    return real_roots, intersections_y


def get_tooth_curves(img):
    """
    获得牙齿区域的曲线
    """
    _, teeth_mask_bg = cv2.threshold(img, 150, 255, cv2.THRESH_BINARY)
    _, bone_mask_bg = cv2.threshold(img, 20, 255, cv2.THRESH_BINARY)
    contours, _ = cv2.findContours(bone_mask_bg, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    # Find the center of the new image
    image_center_bg = (int(teeth_mask_bg.shape[1] / 2), int(teeth_mask_bg.shape[0] / 2))
    # Initialize list to hold the middle points for the new image
    middle_points_bg = []
    num_bins_bg = 100
    # Calculate the angular width of each bin for the new image
    angular_width_bg = 2 * np.pi / num_bins_bg

    # Find the middle points for each bin in the new image
    for i in range(num_bins_bg):
        angle_low_bg = i * angular_width_bg
        angle_high_bg = (i + 1) * angular_width_bg
        middle_point_bg = find_middle_point(teeth_mask_bg, image_center_bg, angle_low_bg, angle_high_bg)
        if middle_point_bg is not None:
            middle_points_bg.append(middle_point_bg)

    middle_points_bg = np.array(middle_points_bg)
    # 寻找最大轮廓
    max_contour = max(contours, key=cv2.contourArea)
    x, y, w, h = cv2.boundingRect(max_contour)
    # 左下角坐标 (x, y + h)
    left_bottom = (x, y + h)
    # 右下角坐标 (x + w, y + h)
    right_bottom = (x + w, y + h)
    bone_points_x = np.array(
        [left_bottom[0], middle_points_bg[:, 1].min(), middle_points_bg[:, 1].max(), right_bottom[0]])
    bone_points_y = np.array(
        [left_bottom[1], middle_points_bg[:, 0].max(), middle_points_bg[:, 0].max(), right_bottom[1]])
    p1 = np.polyfit(middle_points_bg[:, 1], middle_points_bg[:, 0], 4)
    p2 = np.polyfit(bone_points_x, bone_points_y, 2)
    return p1, p2


def curves_smooth(cross_x, x_range_bone, y_range_tooth, y_range_bone, p1, p2):
    """
    将两条抛物线交叉处平滑连接
    cross_x: 交叉点x坐标
    cross_y: 交叉点y坐标
    x_range_bone:抛物线的x坐标，两条抛物线的x坐标相同
    y_range_tooth：牙齿曲线的y坐标
    y_range_bone：颌骨曲线的y坐标
    p1: 牙齿曲线
    p2: 颌骨曲线
    """
    # 假设交点的x坐标大约在x_intersect范围内，我们需要确定这个范围
    x_intersect_start = cross_x[0] - 15  # 这个值需要您根据实际情况调整
    x_intersect_end = cross_x[0] + 15  # 这个值需要您根据实际情况调整

    # 计算在交点附近的权重
    weights = cosine_weight(x_range_bone, x_intersect_start, x_intersect_end)
    y_intersect_end_tooth = np.polyval(p1, x_intersect_end)
    y_intersect_end_bone = np.polyval(p2, x_intersect_end)
    y_d = -abs(y_intersect_end_bone - y_intersect_end_tooth)
    weights_add = add_weight(x_range_bone, x_intersect_start, x_intersect_end, y_d)
    y_range_weighted = y_d * weights + weights_add

    # 合并两条曲线
    y_combined = y_range_bone + y_range_weighted

    # 确保在交点范围外的曲线保持不变
    y_combined[x_range_bone < x_intersect_start] = y_range_bone[x_range_bone < x_intersect_start]
    y_combined[x_range_bone > x_intersect_end] = y_range_tooth[x_range_bone > x_intersect_end]

    # 处理第二个交点
    x_intersect_start = cross_x[1] - 15
    x_intersect_end = cross_x[1] + 15

    # 计算在交点附近的权重
    weights = cosine_weight(x_range_bone, x_intersect_start, x_intersect_end)
    y_intersect_start_tooth = np.polyval(p1, x_intersect_start)
    y_intersect_start_bone = np.polyval(p2, x_intersect_start)
    y_d = -abs(y_intersect_start_bone - y_intersect_start_tooth)
    weights_add = add_weight(x_range_bone, x_intersect_start, x_intersect_end, y_d, L=False)
    y_range_weighted = y_d * weights + weights_add

    # 合并两条曲线
    y_combined1 = y_range_bone + y_range_weighted

    # 确保在交点范围外的曲线保持不变
    y_combined1[x_range_bone < x_intersect_start] = y_combined[x_range_bone < x_intersect_start]
    y_combined1[x_range_bone > x_intersect_end] = y_range_bone[x_range_bone > x_intersect_end]
    return y_combined1

if __name__ == '__main__':
    image_path_bg = r'D:\Tooth-cls_feature\img_bg.png'
    image_bg = cv2.imread(image_path_bg, 0)
    # 得到牙齿和颌骨拟合出的抛物线
    p1, p2 = get_tooth_curves(image_bg)
    # 得到两条曲线的交叉点坐标
    cross_x, cross_y = find_cross_point(p1, p2)

    # 绘制已知点和插值曲线
    x_range_tooth = np.linspace(0, 334, 500)
    y_range_tooth = np.polyval(p1, x_range_tooth)
    x_range_bone = np.linspace(0, 334, 500)
    y_range_bone = np.polyval(p2, x_range_bone)

    y_combined1 = curves_smooth(cross_x, x_range_bone, y_range_tooth, y_range_bone, p1, p2)

    img_bg = Image.fromarray(image_bg)
    plt.imshow(img_bg)
    plt.plot(x_range_tooth, y_combined1, color='red', label='Combined Curve', linestyle='--')
    plt.legend()
    plt.show()

    plt.imshow(img_bg)
    plt.plot(x_range_tooth[80:400], y_range_tooth[80:400], color='blue', label='Tooth Curve')  # 绘制插值曲线
    plt.plot(x_range_bone, y_range_bone, color='red', label='Bone Curve')  # 绘制插值曲线
    plt.legend()
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('Quadratic Polynomial Interpolation')
    plt.legend()
    plt.grid(True)
    plt.show()