Cohen-Sutherland 3D线裁剪算法陷入无限循环问题求助

问题背景 #

基于维基百科的Cohen–Sutherland算法实现3D齐次坐标下的线裁剪功能，多数场景运行正常，但部分线段会陷入无限循环。例如：

• 线段端点：
  P1=(1.8855625490365544, 2.159434556152527, 1.8855426244586575, 1.8855625490365544)
  P2=(0.29644864320281183, -1.8363609313964844, -1.5815350850164043, -1.5815150217554024)
• 循环现象：P1因y>w（TOP区域外）被裁剪，裁剪后又因x>w（RIGHT区域外）再次裁剪，之后又回到TOP区域外状态，反复循环无法终止。

前10次迭代输出 #

（此处补充实际迭代输出内容，示例格式：）
迭代1：裁剪P1（TOP）后坐标=(xxx, xxx, xxx, xxx)，区域码=TOP
迭代2：裁剪P1（RIGHT）后坐标=(xxx, xxx, xxx, xxx)，区域码=TOP
迭代3：裁剪P1（TOP）后坐标=(xxx, xxx, xxx, xxx)，区域码=RIGHT
...

实现代码 #

# 替换为你的实际实现代码
def cohen_sutherland_clip_3d(p1, p2):
    def compute_code(p):
        x, y, z, w = p
        code = 0
        # 区域码定义：RIGHT(8), LEFT(4), TOP(2), BOTTOM(1), FRONT(16), BACK(32)
        if x > w:
            code |= 0b1000
        if x < -w:
            code |= 0b0100
        if y > w:
            code |= 0b0010
        if y < -w:
            code |= 0b0001
        if z > w:
            code |= 0b10000
        if z < -w:
            code |= 0b01000
        return code

    code1 = compute_code(p1)
    code2 = compute_code(p2)

    while True:
        if code1 == 0 and code2 == 0:
            return (p1, p2)
        if (code1 & code2) != 0:
            return None

        code_out = code1 if code1 != 0 else code2
        p_out = p1 if code_out == code1 else p2
        x, y, z, w = p_out

        # 计算与裁剪平面的交点参数t
        t = 0.0
        if code_out & 0b0010:  # TOP: y = w
            t = (w - y) / ((p2[1] - p2[3]) - (p1[1] - p1[3]))
        elif code_out & 0b0001:  # BOTTOM: y = -w
            t = (-w - y) / ((p2[1] + p2[3]) - (p1[1] + p1[3]))
        elif code_out & 0b1000:  # RIGHT: x = w
            t = (w - x) / ((p2[0] - p2[3]) - (p1[0] - p1[3]))
        elif code_out & 0b0100:  # LEFT: x = -w
            t = (-w - x) / ((p2[0] + p2[3]) - (p1[0] + p1[3]))
        elif code_out & 0b10000:  # FRONT: z = w
            t = (w - z) / ((p2[2] - p2[3]) - (p1[2] - p1[3]))
        elif code_out & 0b01000:  # BACK: z = -w
            t = (-w - z) / ((p2[2] + p2[3]) - (p1[2] + p1[3]))

        # 生成新端点
        new_x = p1[0] + t * (p2[0] - p1[0])
        new_y = p1[1] + t * (p2[1] - p1[1])
        new_z = p1[2] + t * (p2[2] - p1[2])
        new_w = p1[3] + t * (p2[3] - p1[3])
        new_p = (new_x, new_y, new_z, new_w)

        # 更新端点和区域码
        if code_out == code1:
            p1 = new_p
            code1 = compute_code(p1)
        else:
            p2 = new_p
            code2 = compute_code(p2)

问题根源分析 #

1. 浮点数精度误差：齐次坐标计算交点时，浮点数舍入误差导致裁剪后的点仍被判定为区域外（比如y略大于w），触发重复裁剪。
2. 交点计算逻辑错误：3D齐次坐标下的参数t推导可能存在符号或公式错误，导致交点未真正落在裁剪平面上。
3. 区域码判断无容错：严格的>/<判断未考虑精度误差，微小数值偏差会导致误判区域。

解决方案 #

1. 添加精度阈值 #

修改区域码计算函数，引入极小epsilon（如1e-8）避免浮点误差误判：

def compute_code(p):
    x, y, z, w = p
    eps = 1e-8
    code = 0
    if x > w + eps:
        code |= 0b1000
    if x < -w - eps:
        code |= 0b0100
    if y > w + eps:
        code |= 0b0010
    if y < -w - eps:
        code |= 0b0001
    if z > w + eps:
        code |= 0b10000
    if z < -w - eps:
        code |= 0b01000
    return code

2. 修正交点计算 #

确保t的推导符合齐次坐标裁剪平面定义，以TOP平面（y = w）为例，正确公式应为：

# TOP平面：y(t) = w(t) → P1.y + t*(P2.y-P1.y) = P1.w + t*(P2.w-P1.w)
t = (p1[3] - p1[1]) / ((p2[1] - p2[3]) - (p1[1] - p1[3]))

检查所有裁剪平面的t计算，确保公式符号和分子分母正确。

3. 限制迭代次数 #

在循环中添加最大迭代次数限制，防止无限循环：

max_iterations = 20
iter_count = 0
while iter_count < max_iterations:
    iter_count += 1
    # 原有循环逻辑...
else:
    print("裁剪陷入循环，已终止")
    return None

4. 归一化齐次坐标 #

裁剪后对新点进行归一化（除以w），减少精度累积：

# 生成新端点后添加归一化步骤
if new_w != 0:
    norm_factor = 1.0 / new_w
    new_p = (new_x * norm_factor, new_y * norm_factor, new_z * norm_factor, 1.0)

内容的提问来源于stack exchange，提问作者Shail