#ifndef lint static char sccsid[] = "@(#)arctochain.c 9.5 88/01/19 Copyright 1985 Sun Micro"; #endif /* * Copyright (c) 1984 by Sun Microsystems, Inc. */ /*- Convert a conic arc to a curve arctochain.c, Tue Jun 11 14:23:17 1985 "Elegance and truth are inversely related." -- Becker's Razor Vaughan Pratt Sun Microsystems This code is a version of the conix package that has been heavily massaged by James Gosling. */ #include #ifdef REF #include #include #endif #include "shape.h" struct shape_triangle { /* A triangle in rounded coordinates */ struct spoint A, B, C; }; /* * The conic equation for an arc gets transformed into the implicit * equation which is a homogeneous quadratic in two variables: * * alpha*x**2+2*beta*x*y+gamma*y**2+delta*x+epsilon*y+zeta == 0 * */ struct shape_implicit { int alpha; int beta; int gamma; int delta; int epsilon; int zeta; }; static short gshape_xpos, gshape_ypos; /* Some unpleasant global * variables */ #define MAXNUM fracti(15) #define MAXDEN fracti(25) #define EPSILON fraction(1,256) #define subxy(p, q) {(p).x -= (q).x; (p).y -= (q).y;} #define sprod(a,b) ((short)(a) * (short)(b)) /* macros for the conic machines */ #define out(b) { \ if (b) \ xpos += dir; \ else \ *ypos++ = xpos; \ } #define across(rel,q1) { \ out(1); \ Fx += Fxx, Fy += Fyx; \ if ((F += Fx) rel 0) \ goto q1; \ } #define down(rel,q1) { \ out(0); \ if (--rem0 < 0) \ goto final; \ Fx -= Fxy, Fy -= Fyy; \ if ((F -= Fy) rel 0) \ goto q1; \ } #define trans1(axis, q1) \ for (;;) \ axis(>=, q1); #define trans2(axis, q1, rel, q2) \ for (;;) { \ axis(<, q1); \ if (Fx rel Fx_err) \ goto q2; \ } int arc_subdivision_depth; /* * Converts an arc and adds it to a chain. Assumes that the arc goes down * (A.y<=B.y<=C.y) and is no more than a quadrant (A.x>=B.x>=C.x || * A.x<=B.x<=C.x). Also assumes that the chain has enough room. */ arc_to_chain(a, curve) struct shape_arc *a; struct shape_dcurve *curve; { struct shape_triangle M; struct shape_implicit imp; short u, v; struct spoint cent, size; M.A.x = roundfr(a->A.x), M.A.y = -roundfr(a->A.y); M.B.x = roundfr(a->B.x), M.B.y = -roundfr(a->B.y); M.C.x = roundfr(a->C.x), M.C.y = -roundfr(a->C.y); size.x = M.C.x - M.A.x, size.y = M.C.y - M.A.y; if (size.y >= 0) return; assert(M.A.x >= M.B.x && M.B.x >= M.C.x || M.A.x <= M.B.x && M.B.x <= M.C.x || a->sh == 0); /*- assert(M.A.y >= M.B.y && M.B.y >= M.C.y); This assertion may fail when there are numerical errors in the subdivision process */ if (!(M.A.y >= M.B.y && M.B.y >= M.C.y)) M.B.y = M.A.y; /* Assuming: size.y<=0 */ /* subdivide down to implicitizable size */ if (a->sh != 0 && (abs(size.x) > 150 || size.y < -150 || a->sh > fraction(138, 100) || a->sh < fraction(1, 5))) { struct shape_arc first, second; extern short fract_overflows; if (arc_subdivision_depth > 10) { this_curve_is_trash: curve->y0 = 0; curve->y1 = -1; curve->x0 = 0; curve->x1 = -1; return; } fract_overflows = 0; ++arc_subdivision_depth; halfdivide(a, &first, &second); if (fract_overflows) { --arc_subdivision_depth; goto this_curve_is_trash; } arc_to_chain(&first, curve); if (curve->y0 <= curve->y1) arc_to_chain(&second, curve); --arc_subdivision_depth; return; } if (a->sh != 0) rationalize(a->sh, &u, &v); else { u = 0; v = 1; } /* translate origin to center of triangle */ cent.x = frrsh(M.A.x + M.C.x, 1); /* find center of tri */ cent.y = frrsh(M.A.y + M.C.y, 1); gshape_xpos = M.A.x; gshape_ypos = -M.A.y; M.A.x -= cent.x, M.B.x -= cent.x, M.C.x -= cent.x; M.A.y -= cent.y, M.B.y -= cent.y, M.C.y -= cent.y; shape_implicitize(&M, u, v, &imp); /* implicitize */ shape_pitteway(&imp, M.A, size, imp.delta + imp.alpha + (imp.beta >> 1) + (M.A.y == M.B.y ? 2 * imp.alpha * M.C.x + imp.beta * M.C.y : 2 * imp.alpha * M.A.x + imp.beta * M.A.y) ,curve); return; } static rationalize(s, u, v) fract s; short *u, *v; { register fract ul, um, ur, sl, sm, sr; #ifdef mc68000 #ifdef SUN register *vlA, *vmA, *vrA; # define vl ((int)vlA) # define vm ((int)vmA) # define vr ((int)vrA) #else register vl, vm, vr; #endif #else register vl, vm, vr; #endif ul = fracti(0), ur = fracti(1); vl = 1, vr = 0; sr = fracti(0), sl = frmul(s, s); for (;;) { um = ul, um += ur, vm = vl, vm += vr, sm = sl, sm += sr; if (um > MAXNUM || vm > MAXDEN) { /* bounded rational */ if (sm < um) *u = floorfr(ul), *v = vl; else *u = floorfr(ur), *v = vr; return; } if (sm < um - EPSILON) ur = um, vr = vm, sr = sm; else if (sm > um + EPSILON) ul = um, vl = vm, sl = sm; else { *u = floorfr(um), *v = vm; return; } } } #define FUZZ fraction(1,50) #define close(a,b) (((unsigned)(a-b+FUZZ))<2*FUZZ) halfdivide(a, a1, a2) register struct shape_arc *a, *a1, *a2; { fract sBx, sBy, ss; struct fpoint D, E, F; a1->next = a2, a2->next = 0; a2->C = a->C, a1->A = a->A; sBx = vfrmul(a->sh, a->B.x); sBy = vfrmul(a->sh, a->B.y); D.x = frrsh(a->A.x + sBx, 1); D.y = frrsh(a->A.y + sBy, 1); E.x = frrsh(a->C.x + sBx, 1); E.y = frrsh(a->C.y + sBy, 1); F.x = frrsh(D.x + E.x, 1); F.y = frrsh(D.y + E.y, 1); ss = frrsh(fracti(1) + a->sh, 1); a1->sh = a2->sh = frsqrt(ss); a1->B.x = close(a->B.x, a->A.x) ? a->A.x : vfrdiv(D.x, ss); a1->B.y = close(a->B.y, a->A.y) ? a->A.y : vfrdiv(D.y, ss); a1->C.x = vfrdiv(F.x, ss); a1->C.y = vfrdiv(F.y, ss); if (!(a1->A.x >= a1->B.x && a1->B.x >= a1->C.x || a1->A.x <= a1->B.x && a1->B.x <= a1->C.x)) if (close(a1->A.x, a1->B.x)) a1->B.x = a1->A.x; else a1->B.x = a1->C.x; if (!(a1->A.y <= a1->B.y && a1->B.y <= a1->C.y)) if (close(a1->A.y, a1->B.y)) a1->B.y = a1->A.y; else a1->B.y = a1->C.y; a2->A = a1->C; a2->B.x = close(a->B.x, a->C.x) ? a->C.x : vfrdiv(E.x, ss); a2->B.y = close(a->B.y, a->C.y) ? a->C.y : vfrdiv(E.y, ss); if (!(a2->A.x >= a2->B.x && a2->B.x >= a2->C.x || a2->A.x <= a2->B.x && a2->B.x <= a2->C.x)) if (close(a2->A.x, a2->B.x)) a2->B.x = a2->A.x; else a2->B.x = a2->C.x; if (!(a2->A.y <= a2->B.y && a2->B.y <= a2->C.y)) if (close(a2->A.y, a2->B.y)) a2->B.y = a2->A.y; else a2->B.y = a2->C.y; } /* * assume arc fits in box [-150,150]x[-150,150], s <= 3 Result is exact, * essential to preserve exact semantics when modifying code. */ shape_implicitize(m, u, v, imp) register struct shape_triangle *m; short u, v; register struct shape_implicit *imp; { short ax, ay, az, bx, by, bz, cx, cy, cz; register short uax, uay, vbx, vby, ucx, ucy; int uaz, vbz, uaxcy, uaycx; /* NB: a>>1 = floor(a/2) even for negative a, depended on below */ if (u == 0) { /* straight line: special case */ straight:imp->alpha = 0, imp->beta = 0, imp->gamma = 0; imp->delta = m->A.y - m->C.y, imp->epsilon = m->C.x - m->A.x; /* should reverse sign of delta HERE for upward arcs */ imp->zeta = ((-imp->delta - imp->epsilon) >> 1) - (short) imp->delta * m->C.x - (short) imp->epsilon * m->C.y; return; } /* invert m treated as a 3x3 matrix with Az = Bz = Cz = 1 */ ax = m->B.y - m->C.y, ay = m->C.x - m->B.x; /* in [-150,150] */ bx = m->C.y - m->A.y, by = m->A.x - m->C.x; /* in [-150,150] */ cx = m->A.y - m->B.y, cy = m->B.x - m->A.x; /* in [-150,150] */ az = m->B.x * m->C.y - m->C.x * m->B.y; /* in [-11250,11250] */ bz = m->C.x * m->A.y - m->A.x * m->C.y; /* in [-11250,11250] */ cz = m->A.x * m->B.y - m->B.x * m->A.y; /* in [-11250,11250] */ uax = u * ax, uay = u * ay; ucx = u * cx, ucy = u * cy; /* in [-2250,2250] */ vbx = v * bx, vby = v * by; /* in [-3750,3750] */ uaz = u * az, vbz = v * bz; /* long */ uaxcy = uax * cy, uaycx = uay * cx; if (uaxcy == uaycx) /* s != 0 yet arc may still degenerate * into */ goto straight; /* straight line when discretizing A,B,C */ /* expand v*b**2-u*4*a*c using inverse */ imp->alpha = vbx * bx - 4 * (uax * cx); imp->gamma = vby * by - 4 * (uay * cy); imp->beta = 2 * (vbx * by) - 4 * (uaxcy + uaycx); imp->delta = vbx * (short) (2 * bz - bx - by) - sprod(2 * uax, 2 * cz - cx - cy) - sprod(2 * ucx, 2 * az - ax - ay); imp->epsilon = vby * (short) (2 * bz - bx - by) - sprod(2 * uay, 2 * cz - cx - cy) - sprod(2 * ucy, 2 * az - ax - ay); imp->zeta = (vbz * (bz - bx - by) + (((vbx + vby) * (bx + by)) >> 2)) - (2 * uaz - uax - uay) * (2 * cz - cx - cy); /* force +-+ sandwich. */ if (imp->alpha < 0) { /* negate implicit form */ imp->alpha = -imp->alpha, imp->beta = -imp->beta; imp->gamma = -imp->gamma, imp->delta = -imp->delta; imp->epsilon = -imp->epsilon, imp->zeta = -imp->zeta; } } shape_pitteway(impl, pos, size, side, curve) struct shape_implicit *impl; struct spoint pos; int side; struct spoint size; struct shape_dcurve *curve; { /* initialize F,Fx,Fy,Fxx,Fyy,Fxy to start of curve */ int alx, bex, gay, bey, rem0; register F, Fx, Fx_err; register xpos; register short *ypos; #ifndef fast register Fy, Fxx, Fyy, Fxy; register struct shape_implicit *imp; #else register *FyA, *FxxA, *FyyA, *FxyA; # define Fy ((int)FyA) # define Fxx ((int)FxxA) # define Fyy ((int)FyyA) # define Fxy ((int)FxyA) # define imp ((struct shape_implicit *)FxyA) #endif # define Fyx Fxy register dir; imp = impl; dir = -1; if (size.x > 0) { pos.x++; dir = 1; } size.y = -size.y; ypos = curve->x + (gshape_ypos - curve->y0); if (size.y <= 0) return; if (gshape_ypos < curve->y0 && gshape_ypos > curve->y1) return; xpos = gshape_xpos; pos.y++; alx = imp->alpha * pos.x, bex = imp->beta * pos.x; gay = imp->gamma * pos.y, bey = imp->beta * pos.y; F = (alx + imp->delta) * pos.x + (bex + gay + imp->epsilon) * pos.y + imp->zeta; Fx = (2 * alx) + bey - imp->alpha + imp->delta; Fy = (2 * gay) + bex + imp->gamma + imp->epsilon; Fxx = imp->alpha << 1; Fyy = imp->gamma << 1; Fxy = imp->beta; Fx_err = -Fxx >> 1; if (size.x <= 0) /* to make the across macro work for * leftbound */ Fx += Fxx, Fx = -Fx, Fxy = -Fxy; /* enter conic machine at appropriate state */ rem0 = size.y; if ((size.x > 0) == (side < 0)) goto EOD; else goto OED; /* OE state machine */ ODE:trans1(across, OED); DOE:trans2(across, ODE, >, OED); OED:trans2(down, ODE, <, DOE); /* EO state machine */ EDO:trans1(down, DEO); DEO:trans2(across, EDO, >, EOD); EOD:trans2(down, EDO, <, DEO); final:; /*- clearasil(buf0+1, bitlen); /* clean out 0011/1100 pimples */ assert(ypos <= curve->x + (curve->y1 - curve->y0 + 1)); }