Bezier交Bezier实现_思路

实现思路主要分为两个部分，一部分是求出两条曲线的局部最小距离点对，并以此作为后续算法的初始点；另一部分是使用牛顿迭代法精细化上一步提供的初始点，从而得到交点（若收敛则得到交点，不收敛则舍弃）。同时通过局部最小距离点对的距离与误差上限的关系来加速求交，排除明显无交的局部最小点对。

关于局部最小距离点对

使用遍历法查找局部最小点对的逻辑

  1. 输入两个样条线
  2. 将每个样条线分隔为 n 段
  3. 计算子样条的准控制多边形
  4. 遍历每个子样条求出原样条曲线的最近距离点对
  5. 检查子样条曲线的最近距离是否小于子样条和准控制多边形的误差上限之和：
    • 如果是，继续下一步
    • 如果否，排除该子样条对有交
  6. 检查最近距离是否小于初始点设置误差：
    • 如果是，继续下一步
    • 如果否，排除该最近距离点对
  7. 输出初始点

使用竞争流查找局部最小点对的逻辑

  1. 输入两条样条曲线
  2. 在曲线上选取离散点，每条曲线上均匀采样点
  3. 端点匹配检查，应该从相距最近的两个端点开始竞争流，否则竞争流会出现问题
  4. 调用竞争流算法，寻找局部最小距离的点对，并记录索引
    • 输入两组离散点，表示两条曲线的采样点，并初始化变量
    • 竞争流迭代过程
      1. 终止条件检查
      2. 计算相邻点的距离
        • 计算 (i+di, j) 和 (i, j+dj) 处的最短距离 dist_i 和 dist_j，用于判断局部最小值。
      3. 判断局部最小值
        • 如果 l 小于 dist_i 和 dist_j，且 is_find_min 标记为 true，则认为当前 (i, j) 是局部最小值。
      4. 更新索引：选择使距离最小的方向前进
      5. 处理索引越界和方向反转
    • 输出局部最小值点对的索引 (i, j) 及其对应的距离 l。
  5. 检查最近距离l是否小于误差上限之和：
    • 如果是，继续下一步
    • 如果否，排除该子样条对有交
  6. 检查最近距离是否小于初始点设置误差：
    • 如果是，继续下一步
    • 如果否，排除该最近距离点对
  7. 输出初始点

关于牛顿迭代

牛顿迭代法通常用于求解非线性方程组，其基本思想是通过迭代线性化方程组来逐步逼近方程组的解。对于两条贝塞尔曲线，我们希望找到参数 param1 和 param2，使得两条曲线在这些参数下的点尽可能接近。

下面介绍二维平面使得情况，三维类似。
$$
f(t,s)=x_1(t)-x_2(s)=0
$$

$$
g(t,s)=y_1(t)-y_2(s)=0
$$

在上式中$x_1(t)$是曲线的参数方程，
$$
x_1(t) = \sum_{i=0}^{n} \binom{n}{i} t^i (1-t)^{n-i} Px_i = \sum_{i=0}^{n} B_i(t) Px_i, \quad t \in [0,1]
$$
参数方程的函数值就是曲线$1$上的参数为$t$的一个点的$x$坐标，以此类推。

对于非线性方程组 $f(x,y)=0$ 和 $g(x,y)=0$，牛顿迭代法的更新公式可以表示为：
$$
\begin{bmatrix} x_{n+1} \ y_{n+1} \end{bmatrix} = \begin{bmatrix} x_n \ y_n \end{bmatrix} - J^{-1}(x_n) \begin{bmatrix} f(x_n, y_n) \ g(x_n, y_n) \end{bmatrix}
$$
其中，$J(x_n)$ 是雅可比矩阵，$J^{-1}(x_n)$ 是雅可比矩阵的逆矩阵，$x_n=\begin{bmatrix} x_n \ y_n \end{bmatrix}$ 是当前迭代点，$f(x_n, y_n)$ 和 $g(x_n, y_n)$ 是两个方程在 $(x_n,y_n)$ 处的函数值。

对于这个问题，雅可比矩阵 $J$ 是：
$$
J=\begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix}
$$
其中，$f$ 和 $g$ 分别是：
$$
f(x,y)=x_1-x_2 \quad g(x,y)=y_1-y_2
$$
这里，$(x_1,y_1)$ 和 $(x_2,y_2)$ 是两条贝塞尔曲线在参数 param1 和 param2 下的点。

代码中，迭代公式被实现为：

double det = jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0]
double delta_param1 = (f * jacobian[1][1] - g * jacobian[0][1]) / det;
double delta_param2 = (g * jacobian[0][0] - f * jacobian[1][0]) / det;

param1 -= delta_param1;
param2 -= delta_param2;