[metapost] Re: new is_clockwise routine

[Dan Luecking]

At 05:23 AM 11/27/2005, you wrote:
>Sunday, November 27, 2005 Paul Pichaureau wrote:
>
> > Le samedi 26 novembre 2005 à 19:51:46, vous écriviez :
>
>
>GB>> In general, to find the clockwiseness of a single Bézier
>GB>> cubics it suffices to check for the sign of
>
>GB>> (P1 - P0) × (P2 - P1)
>
>GB>> and
>
>GB>> (P2 - P1) × (P3 - P2)
>
>GB>> (P0, P1, P2, P3 are the points that define the curve: P0 ..
>GB>> controls P1 and P2 .. P3)
>
> > I've tried to implement an algorithm of this type.
>
> > Hélas, it's quite impossible in metapost, due to the limitation of the
> > magnitude of real number. For example, with the bezier curve
>
> > (-50,-50).. controls (50,150) and (150,50) ... (300,300)
>
> > metapost have to compute
>
> > (100,200) x (100,-100) =  -10 000 - 20 000 < 4 096   BANG !
>
> > and
>
> > (100,-100) x (150,250) = 25000 + 15 000 > 4 096 re-BANG !
>
> > The worst is the result can be really simple (e.g. 0 or -1) but the
> > metapost refuse even intermediate numbers greater than 4096...
>
>I suspected you would come across a problem like this.
>However, the solution is simple because we are interested in
>the *sign* of the cross product, so instead of trying to
>calculate a×b you just calculate dir(a) × dir(b): the sign
>is the same but the vectors have unitary length and thus the
>cross product will not bang.

You mean (unitvector a crossproduct unitvector b). This, of course,
gives (sind A) where A is the angle from a to b (A between -180 and 180).
That will have the same sign as A itself and you can also get that
angle directly:
   angle (b rotated - angle a);
(which is generally better than angle b - angle a).

If one really needs the cross product, one can find the cross product
of the unitvectors, then multiply that by (abs a) and then (perhaps
first checking for overflow) multiply by (abs b).

One can set
   warningcheck := 0;
and intermediate results can be as high as 32765. Try
   t := 1000; s := sind 10t;

Giuseppe: What you called a knot is what I would call a loop.
But what is the "knot factor"?  If that just the discriminant
of g'xg''?


Dan


Daniel H. Luecking
Department of Mathematical Sciences
University of Arkansas
"Be kind. Every person you meet is fighting a hard battle." - Anon.