Given two 3D Bezier surfaces, determine their intersections.
**algorithm** extended-newton-method **is
input:** Bezier surface F, G
**output:** initial point on F ∩ G
create arbitrary guess point p
S <- F - G
**do**
J <- Jacobian(S)
pnew <- p - pinv(J(p)) · S(p) // pinv is pseudo inverse
t <- norm(pnew, p)
**while** t >= e // e is threshold
**return** pnew or failure to converge
The two Bezier surfaces F and G are two parametric equations. To find their intersection, we need to solve for the roots to the system of equations yielded by F - G = 0 . Extended Newton’s method offers an iterative way to converge to the roots starting from some initial guess points.
**alorithm** trace **is**
**input:** start point p,
Bezier surface F(u, v), G(s, t)
**output:** vector of intersection points vec
N1 <- Fu × Fv / |Fu × Fv| // N1 is unit normal vector for F
N2 <- Gs × Gt / |Gs × Gt| // N2 is unit normal vector for G
T <- N1 × N2 / |N1 × N2| // tangent unit vector at point p
**do**
pnew = root of equation below
t <- norm(pnew, p)
vec <- pnew
**while** t < e // e is threshold
**return** vec
The equation mentioned in above pseudo code is:
$$ \begin{bmatrix}\frac{\partial F}{\partial u} & \frac{\partial F}{\partial v} & -\frac{\partial G}{\partial s} & -\frac{\partial G}{\partial t}\\ T \cdot \frac{\partial F}{\partial u}& T \cdot \frac{\partial F}{\partial v} & 0 & 0\end{bmatrix}\begin{bmatrix}u_1 - u_0\\ v_1 - v_0\\ s_1 - s_0\\ t_1 - t_0\end{bmatrix} = \begin{bmatrix}0\\ 0\\ 0\\ d\end{bmatrix} $$
We solve for the new point’s coordinate [u1, v1, s1, t1] with the marching distance d. Note here we can tune d according to local geometry in order to achieve different rate of marching update. For example, when there is high curvature, we can choose a smaller d to get a finer approximation.
