[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Devel] Problems with bbox code and cubic bezier curves
From: |
Werner LEMBERG |
Subject: |
Re: [Devel] Problems with bbox code and cubic bezier curves |
Date: |
Mon, 23 Apr 2001 08:20:00 +0200 (CEST) |
> 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);
}
- Re: [Devel] Problems with bbox code and cubic bezier curves, (continued)
- Re: [Devel] Problems with bbox code and cubic bezier curves, Toby J Sargeant, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Toby J Sargeant, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Nathan Hurst, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/24
- Re: [Devel] Problems with bbox code and cubic bezier curves,
Werner LEMBERG <=
- Re: [Devel] Problems with bbox code and cubic bezier curves, Nathan Hurst, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Werner LEMBERG, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Toby J Sargeant, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, David Turner, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, David Turner, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Tom Kacvinsky, 2001/04/23
- Re: [Devel] Problems with bbox code and cubic bezier curves, Werner LEMBERG, 2001/04/23