From: horst.kraemer@snafu.de (Horst Kraemer)
Subject: Re: Root finding problem
Date: Thu, 06 May 1999 18:45:48 GMT
Newsgroups: sci.math.num-analysis
Keywords: intersection points of two cubic parametric curves

On Fri, 7 May 1999 16:28:21 -0700, "Des Grenfell"
<des@grenfell.demon.co.uk> wrote:

> I really need to find all the intersection points of two cubic parametric
> lines in x and y.
> I tried various ways to reduce this to a univariate problem but they all
> prove very time-consuming so I'm sticking with a pair of bivariate cubic
> equations.
 
> I find that I can get the first root very quickly using Newton Raphson but
> the problem is trying to get the other roots - one root often tends to act
> like a magnet for the iteration scheme no matter what starting points I
> take.
 
> Does anyone know of any robust root bracketing or estimating schemes that I
> could use to improve my starting guesses and hopefully improve my algorithm.
> I did think about some form of subdivision but it looks very time consuming.

There is a subdivision algorithm which may be surprisingly efficient.

You may transform either polynomial curve into BEZIER format using the
BERNSTEIN polynomials of order 4 as base polynomials:

P(t) = A1*(1-t)^3 + A2*3t(1-t)^2 + A3*3t^2(1-t) + A4*t^3 , 0<=t<=1

Then P(0) = A1 ; P(1) = A4

The key property of this representation is: the whole curve lies in
the convex hull of the "control points" A1,A2,A3,A4.

It is very easy to split this curve into two curves at the parameter
value t_mid = 1/2 by using CASTELJAUs evaluation algorithm for t=1/2

A1
    (A1+A2)/2
A2             (A1+2A2+A3)/4
    (A2+A3)/2                 (A1+3A2+3A3+A4)/8 = P(1/2)
A3             (A2+2A3+A4)/4
    (A3+A4)/2 
A4


Every entry is the average of the left neighbours.

The upper side of the triangle

A1 , (A1+A2)/2 , (A1+2A2+A3)/4 ,  (A1+3A2+3A3+A4)/8

is the set of BEZIER control points of the first part, the lower side
of the triangle

(A1+3A2+3A3+A4)/8 , (A2+2A3+A4)/4 , (A3+A4)/2 , A4

is the set of BEZIER control points of the second half of the original
curve.

Now you may enclose the convex hull of any control set (C1,C2,C3,C4)
by a coordinate box using minx[i],maxx[i],miny[i],maxy[i].

This could be used for a nice recursive algorithm along the lines

procedure intersect(A1,A2,A3,A4,B1,B2,B3,B4);
begin        
  if intersection(Hull(A1,A3,A3,A4),Hull(B1,B2,B3,B4)) = empty then
    return
  else if Hulls are smaller than Threshold then
    return Intersection
  else begin
    split(A1,A2,A3,A4,AA1,AA2,AA3,AA4,AA5,AA5,AA7);
    split(B1,B2,B3,B4,BB1,BB2,BB3,BB4,BB5,BB6,BB7);
    intersect(AA1,AA2,AA3,AA4,BB1,BB2,BB3,BB4);
    intersect(AA1,AA2,AA3,AA4,BB4,BB5,BB6,BB7);
    intersect(AA4,AA5,AA6,AA7,BB1,BB2,BB3,BB4);
    intersect(AA4,AA5,AA6,AA7,BB4,BB5,BB6,BB7);
  end;


Regards
Horst