Re: [Devel] Problems with bbox code and cubic bezier curves

[Werner LEMBERG]

> Toby pointed out that the bounding box of a bezier can be computed by
> computing the max and min of the x and y coords of the parametric equation
> of the curve, which for a cubic would occur at either the end points or
> the solutions of:
> t^2[A-3B-C-3D] + t[2B+6D] + [C-3D] 
> where A,B,C,D are the x or y coords of the control points.
> 
> This requires 2 sqrts per curve, but I suspect that that would be
> cheaper than an unknown number of splits?

Exactly.  It is possible to construct third order Bézier curves which
need o(n) splits, where n is the length or height of the coordinate
bounding box.  Even in real life this can happen with rotated
outlines.  Attached is a code snipped from Richard Kinch (he was so
kind to send it to me) which solves the problem.  He uses real floats,
and I haven't found time to convert it to fixed floats (and,
admittedly, I'm a bit weak in doing that).


    Werner


======================================================================


/*
 *      Copyright 1998 by Richard J Kinch
 *      6994 Pebble Beach Ct
 *      Lake Worth FL 33467 USA
 *      Tel (561) 966-8400
 *      FAX (561) 966-0962
 *      address@hidden
 *      http://truetex.com
 */

/* Kinch, November 1998                                                 */

#include <stdio.h>
#include <math.h>
#include "trig.h"
#include "poly.h"

#define DEBUG 0

int
quadratic_roots(double c[3],double root[2]) {
    /* Solve, using closed form methods, the quadratic polynomial:      */
    /*          c[2]*x^2 + c[1]*x + c[0] = 0                            */
    /* for 2 real roots returned in root[0..1].  If error we return 0.  */
    /* We also return 0 or 1 real roots as appropriate, such as when    */
    /* the problem is actually linear.                                  */
    /* After _Numerical Recipes in C_, 2nd edition, Press et al.,       */
    /* page 183, although with some added case testing and forwarding.  */
    /* This is better than the _Graphics Gems_ technique, which admits  */
    /* the possibility of numerical errors cited in Press.              */
    double bb, q;
    int retval=0;
#if DEBUG
    fprintf(stderr,"q_r: quadratic is %lg*t^2 + %lg*t + %lg ; q^2=%lg\n",
            c[2],c[1],c[0],q);
#endif
    /* If problem is actually linear, return 0 or 1 easy roots          */
    if (fabs(c[2])<Epsilon) {
        if (fabs(c[1])>=Epsilon) root[retval++] = -c[0]/c[1];
        return(retval);
        }
    bb = c[1]*c[1];
    q = bb-4.0*c[2]*c[0];
    if (q<0.0) return(retval);
    q = sqrt(q);
    if (c[1]<0.0) q = -q;
    q = -0.5*(c[1]+q);
    if (fabs(q)>=Epsilon) root[retval++] = c[0] / q;
    if (fabs(c[2])>=Epsilon) root[retval++] = q / c[2];
    return(retval);
    }

int
cubic_roots(double c[4],double root[3]) {
    /* Solve, using closed form methods, the cubic polynomial:          */
    /*          c[3]*x^3 + c[2]*x^2 + c[1]*x + c[0] = 0                 */
    /* for 1 real root returned in root[0], or 3 real roots returned    */
    /* in root[0..2].  If error we return 0.  Note: we alter c[].       */
    /* If the polynomial is actually quadratic or linear (because       */
    /* coefficients c[3] or c[2] are zero), we forward the problem to   */
    /* the quadratic/linear solver and return the appropriate 1 or 2    */
    /* roots.                                                           */
    /* After _Numerical Recipes in C_, 2nd edition, Press et al.,       */
    /* page 184, although with some added case testing and forwarding.  */
    /* This is better than the _Graphics Gems_ technique, which admits  */
    /* the possibility of numerical errors cited in Press.              */
    double A, Q, R, QQQ, RR;
    double theta;
    /* Test for a quadratic or linear degeneracy                        */
    if (fabs(c[3])<Epsilon) {
#if DEBUG
        fprintf(stderr,"c_r: handing off to quadratic_roots()\n");
#endif
        return(quadratic_roots(c,root));
        }
    /* Normalize                                                        */
    c[2] /= c[3]; c[1] /= c[3]; c[0] /= c[3]; c[3] = 1.0;
    /* Compute discriminants                                            */
    Q = (c[2]*c[2]-3.0*c[1])/9.0; QQQ = Q*Q*Q;
    R = (2.0*c[2]*c[2]*c[2]-9.0*c[2]*c[1]+27.0*c[0])/54.0; RR = R*R;
    /* Three real roots                                                 */
    if (RR<QQQ) {
#if DEBUG
        fprintf(stderr,"c_r: three real roots\n");
#endif
        /* This sqrt and division is safe, since RR >= 0, so QQQ > RR,  */
        /* so QQQ > 0.  The acos is also safe, since RR/QQQ < 1, and    */
        /* thus R/sqrt(QQQ) < 1.                                        */
        theta = acos(R/sqrt(QQQ));
        /* This sqrt is safe, since QQQ >= 0, and thus Q >= 0           */
        root[0] = root[1] = root[2] = -2.0*sqrt(Q);
        root[0] *= cos(theta/3.0);
        root[1] *= cos((theta+2*PI)/3.0);
        root[2] *= cos((theta-2*PI)/3.0);
        root[0] -= c[2]/3.0;
        root[1] -= c[2]/3.0;
        root[2] -= c[2]/3.0;
        return(3);
        }
    /* One real root                                                    */
    else {
#if DEBUG
        fprintf(stderr,"c_r: one real root\n");
#endif
        A = -pow(fabs(R)+sqrt(RR-QQQ),1.0/3.0);
        if (A!=0.0) { if (R<0.0) A = -A ; root[0] = A + Q/A; }
        root[0] -= c[2]/3.0;
        return(1);
        }
    }


==============================================================================

/*
 *      Copyright 1994 Richard J Kinch
 *      6994 Pebble Beach Ct
 *      Lake Worth FL 33467 USA
 *      Tel (561) 966-8400
 *      FAX (561) 966-0962
 *      address@hidden
 *      http://truetex.com
 */

/* Kinch, December 1994                                                 */

#include <stdio.h>
#include <math.h>
#include "path.h"
#include "copy.h"
#include "trig.h"
#include "section.h"
#include "plot.h"                       /* For kid()                    */

int
find_extrema(double z0,double z0r,double z1l,double z1,
    double *t1,double *t2) {
    /* Finds extrema in the curve z0--z1, where the z's input are the   */
    /* x coordinates (to find vertical tangents) or y coordinates (to   */
    /* find horizontal tangents).  A curve can have zero, one, or two   */
    /* such tangents in each of the vertical or horizontal directions.  */
    /* We return the number of such tangents (0, 1, or 2) and the       */
    /* time of the first tangent (if any, retval >= 1) in t1 and of the */
    /* second (if any, retval == 2) in t2.  We maintain t1 < t2.        */
    /* We compute the tangent times using the closed form solution to   */
    /* quadratic polynomial which is the slope of the z(t) curve.       */
    /* We must take care of for singularities where: the quadratic      */
    /* coefficient(s) are zero; the solutions in t are outside the      */
    /* range (0,1), and to sort t1 and t2 when there are two solutions. */
    /* See the loose notes, "Inserting extrema's to Bezier's," 1/7/95.  */
    /* Scaffold: the quadratic solution should be more robustly done,   */
    /* as is implemented in inflect.c, after _Numerical Methods in C_.  */
    double A, B, C;             /* Intermediate coefficients            */
    double a, b, c;             /* Quadratic coefficients               */
    double q;                   /* Quadratic discriminant               */

    /* Compute quadratic coefficients after the derivation in notes.    */
    A = z1l - z1; B = z0r - z1l; C = z0 - z0r;
    a = A - 2.0*B + C; b = 2.0 * (B-C); c = C;
    q = b*b - 4.0*a*c;

    /* Avoid divide-by-zero when derivative was linear                  */
    if (fabs(a)<Epsilon) {
        /* Avoid divide-by-zero when derivative was constant. That is,  */
        /* the underlying curve is linear in z and thus has no tangents.*/
        if (fabs(b)<Epsilon) goto no_roots;
        /* Single linear slope solution                                 */
        *t1 = -(c/b);
        goto one_test;
        }
#if DEBUG
    fprintf(stderr,"f_e: A=%lg B=%lg C=%lg a=%lg b=%lg c=%lg q=%lg\n",
            A,B,C,a,b,c,q);
#endif
    /* Two real roots and two tangent times.                            */
    if (q>=Epsilon) {
        q = sqrt(q);
        /* Scaffold: divide by a=zero                                   */
        *t1 = -(b/(2.0*a)) + q/(2.0*a);
        *t2 = -(b/(2.0*a)) - q/(2.0*a);
        /* Sort tangent times                                           */
        if (*t1>*t2) { double t; t = *t2; *t2 = *t1; *t1 = t; }
#if DEBUG
        fprintf(stderr,"f_e: two roots: %lg %lg\n",*t1,*t2);
#endif
        /* Test that we don't return endpoints                          */
        /* Scaffold: numeric tolerances values here and below           */
        if (*t1 < 0.01 || *t1 > 0.99) { *t1 = *t2; goto one_test; }
        if (*t2 < 0.01 || *t2 > 0.99) goto one_test;
        return(2);
        }
    /* One real root and one tangent time                               */
    if (fabs(q)<Epsilon) {
        /* Scaffold: divide by a=zero                                   */
        *t1 = -(b/(2.0*a));
#if DEBUG
        fprintf(stderr,"f_e: one root: %lg\n",*t1);
#endif
        /* Test that single tangent is not at the endpoint              */
one_test:
        if (*t1 < 0.01 || *t1 > 0.99) goto no_roots;
        return(1);
        }
    /* No real roots and no tangents.                                   */
no_roots:
    return(0);
    }

int
extrema_times(struct knot *k,struct knot *kf,double *t) {
    /* Find and return sorted extrema times for curve k--kf in et[4],   */
    /* returning the count.                                             */
    double t1, t2;              /* Roots in dt/dx or dt/dy derivatives  */
    int n=0;                    /* Nr of H and/or V tangents in k--kf   */
    int count;                  /* Count of H or V tangents in k--kf    */
    int i,j;

    /* Must be curve                                                    */
    if (k->type&K_R_ONPOINT) return(0);
    /* Vertical tangents                                                */
    count = find_extrema(k->x,k->rx,kf->lx,kf->x,&t1,&t2);
    if (count) { t[n++] = t1; count--; }
    if (count) { t[n++] = t2; count--; }
    if (count) fprintf(stderr,"e_t: too many (%d) vert tangents!\n",count+2);
    /* Horizontal tangents                                              */
    count = find_extrema(k->y,k->ry,kf->ly,kf->y,&t1,&t2);
    if (count) { t[n++] = t1; count--; }
    if (count) { t[n++] = t2; count--; }
    if (count) fprintf(stderr,"e_t: too many (%d) horz tangents\n",count+2);
    /* Test early return if usual case of no extrema                    */
    if (n==0) return(0);
    /* Bubble sort the times.  Since rarely more than one entry, and    */
    /* never more than four, bubble sorting is appropriate.             */
    for (i=1; i<n; i++) for (j=0; j<i; j++) if (t[j]>t[i]) {
        double tt;
        tt = t[i]; t[i] = t[j]; t[j] = tt;
        }
    /* Remove duplicate or near-duplicate times; this can happen in the */
    /* special case where a curve has a cusp, resulting in a turn so    */
    /* fast that a simultaneous horizontal and vertical tangent exist.  */
    for (i=2; i<n; i++) if (t[i]-t[i-1]<Epsilon) t[i] = t[i-1];
                                        /* Equate any near-equal times  */
    for (i=j=0; i<n; ) {                /* Copy from index i to j       */
        t[j++] = t[i++];
        while (i<n && t[i-1]==t[i]) i++;
        }
    /* Return count of non-duplicate times                              */
    return(j);
    }

int
insert_extrema(struct path *c) {
    /* Insert horizontal/vertical tangent knots in contour c            */
    /* Return count of knots inserted.                                  */
    struct path *p;
    struct knot *k;
    double et[4];               /* Up to 4 extrema may be needed        */
    int n;                      /* Number of extrema in k--kf           */
    int nc;                     /* Sum of extrema inserted in c         */
    /* Do for each path in contour                                      */
    /* Do for each curve in path                                        */
    for (nc=0,p=c; p; p=p->next) for (k=p->knot_list; k && k->next; k=k->next) {
        /* Test if end of open path                                     */
        if (k->type&K_ENDPOINT) break;
        /* Gather the extrema times                                     */
        n = extrema_times(k,k->next,et);
#if DEBUG
        fprintf(stderr,"insert_extrema: (%s) t[%d] =",kid(k),n);
        for (i=0; i<n; i++) fprintf(stderr," %lg",t[i]);
        fprintf(stderr,"\n");
#endif
        /* N-section the curve at the tangency times.                   */
        /* Scaffold: should update k to skip inserted part.             */
        n_sect_curve(k,et,NULL,n);
        /* Sum the total tangents we've found                           */
        nc += n;
        if (k==p->knot_list->last) break;
        }
    return(nc);
    }