home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 2 / AACD 2.iso / AACD / Magazine / UsingPDF / GhostScript / source / gs5.10 / gxpflat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-11-10  |  14.2 KB  |  453 lines
  1. /* Copyright (C) 1997 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gxpflat.c */
  20. /* Path flattening algorithms */
  21. #include "gx.h"
  22. #include "gxarith.h"
  23. #include "gxfixed.h"
  24. #include "gzpath.h"
  25.  
  26. /* Define whether to merge nearly collinear line segments when flattening */
  27. /* curves.  This is very good for performance, but we feel a little */
  28. /* uneasy about its effects on character appearance. */
  29. #define MERGE_COLLINEAR_SEGMENTS 1
  30.  
  31. /* ---------------- Curve flattening ---------------- */
  32.  
  33. #define x1 pc->p1.x
  34. #define y1 pc->p1.y
  35. #define x2 pc->p2.x
  36. #define y2 pc->p2.y
  37. #define x3 pc->pt.x
  38. #define y3 pc->pt.y
  39.  
  40. /*
  41.  * To calculate how many points to sample along a path in order to
  42.  * approximate it to the desired degree of flatness, we define
  43.  *    dist((x,y)) = abs(x) + abs(y);
  44.  * then the number of points we need is
  45.  *    N = 1 + sqrt(3/4 * D / flatness),
  46.  * where
  47.  *    D = max(dist(p0 - 2*p1 + p2), dist(p1 - 2*p2 + p3)).
  48.  * Since we are going to use a power of 2 for the number of intervals,
  49.  * we can avoid the square root by letting
  50.  *    N = 1 + 2^(ceiling(log2(3/4 * D / flatness) / 2)).
  51.  * (Reference: DEC Paris Research Laboratory report #1, May 1989.)
  52.  *
  53.  * We treat two cases specially.  First, if the curve is very
  54.  * short, we halve the flatness, to avoid turning short shallow curves
  55.  * into short straight lines.  Second, if the curve forms part of a
  56.  * character (indicated by flatness = 0), we let
  57.  *    N = 1 + 2 * max(abs(x3-x0), abs(y3-y0)).
  58.  * This is probably too conservative, but it produces good results.
  59.  */
  60. int
  61. gx_curve_log2_samples(fixed x0, fixed y0, const curve_segment *pc,
  62.   fixed fixed_flat)
  63. {    fixed
  64.       x03 = x3 - x0,
  65.       y03 = y3 - y0;
  66.     int k;
  67.  
  68.     if ( x03 < 0 )
  69.       x03 = -x03;
  70.     if ( y03 < 0 )
  71.       y03 = -y03;
  72.     if ( (x03 | y03) < int2fixed(16) )
  73.       fixed_flat >>= 1;
  74.     if ( fixed_flat == 0 )
  75.     {    /* Use the conservative method. */
  76.         fixed m = max(x03, y03);
  77.         for ( k = 1; m > fixed_1; )
  78.           k++, m >>= 1;
  79.     }
  80.     else
  81.     {    const fixed
  82.           x12 = x1 - x2,
  83.           y12 = y1 - y2,
  84.           dx0 = x0 - x1 - x12,
  85.           dy0 = y0 - y1 - y12,
  86.           dx1 = x12 - x2 + x3,
  87.           dy1 = y12 - y2 + y3,
  88.           adx0 = any_abs(dx0),
  89.           ady0 = any_abs(dy0),
  90.           adx1 = any_abs(dx1),
  91.           ady1 = any_abs(dy1);
  92.         fixed
  93.           d = max(adx0, adx1) + max(ady0, ady1);
  94.         /*
  95.          * The following statement is split up to work around a
  96.          * bug in the gcc 2.7.2 optimizer on H-P RISC systems.
  97.          */
  98.         uint qtmp = d - (d >> 2) /* 3/4 * D */ + fixed_flat - 1;
  99.         uint q = qtmp / fixed_flat;
  100.  
  101.         if_debug6('2', "[2]d01=%g,%g d12=%g,%g d23=%g,%g\n",
  102.               fixed2float(x1 - x0), fixed2float(y1 - y0),
  103.               fixed2float(-x12), fixed2float(-y12),
  104.               fixed2float(x3 - x2), fixed2float(y3 - y2));
  105.         if_debug2('2', "     D=%f, flat=%f,",
  106.               fixed2float(d), fixed2float(fixed_flat));
  107.         /* Now we want to set k = ceiling(log2(q) / 2). */
  108.         for ( k = 0; q > 1; )
  109.           k++, q = (q + 3) >> 2;
  110.         if_debug1('2', " k=%d\n", k);
  111.     }
  112.     return k;
  113. }
  114.  
  115. /*
  116.  * Define the maximum number of points for sampling if we want accurate
  117.  * rasterizing.  2^(k_sample_max*3)-1 must fit into a uint with a bit
  118.  * to spare; also, we must be able to compute 1/2^(3*k) by table lookup.
  119.  */
  120. #define k_sample_max min((size_of(int) * 8 - 1) / 3, 10)
  121.  
  122. /*
  123.  * Split a curve segment into two pieces at the (parametric) midpoint.
  124.  * Algorithm is from "The Beta2-split: A special case of the Beta-spline
  125.  * Curve and Surface Representation," B. A. Barsky and A. D. DeRose, IEEE,
  126.  * 1985, courtesy of Crispin Goswell.
  127.  */
  128. private void
  129. split_curve_midpoint(fixed x0, fixed y0, const curve_segment *pc,
  130.   curve_segment *pc1, curve_segment *pc2)
  131. {    /*
  132.      * We have to define midpoint carefully to avoid overflow.
  133.      * (If it overflows, something really pathological is going
  134.      * on, but we could get infinite recursion that way....)
  135.      */
  136. #define midpoint(a,b)\
  137.   (arith_rshift_1(a) + arith_rshift_1(b) + ((a) & (b) & 1) + 1)
  138.     fixed x12 = midpoint(x1, x2);
  139.     fixed y12 = midpoint(y1, y2);
  140.  
  141.     /*
  142.      * pc1 or pc2 may be the same as pc, so we must be a little careful
  143.      * about the order in which we store the results.
  144.      */
  145.     pc1->p1.x = midpoint(x0, x1);
  146.     pc1->p1.y = midpoint(y0, y1);
  147.     pc2->p2.x = midpoint(x2, x3);
  148.     pc2->p2.y = midpoint(y2, y3);
  149.     pc1->p2.x = midpoint(pc1->p1.x, x12);
  150.     pc1->p2.y = midpoint(pc1->p1.y, y12);
  151.     pc2->p1.x = midpoint(x12, pc2->p2.x);
  152.     pc2->p1.y = midpoint(y12, pc2->p2.y);
  153.     if ( pc2 != pc )
  154.       pc2->pt.x = pc->pt.x,
  155.       pc2->pt.y = pc->pt.y;
  156.     pc1->pt.x = midpoint(pc1->p2.x, pc2->p1.x);
  157.     pc1->pt.y = midpoint(pc1->p2.y, pc2->p1.y);
  158. #undef midpoint
  159. }
  160.  
  161. /*
  162.  * Flatten a segment of the path by repeated sampling.
  163.  * 2^k is the number of lines to produce (i.e., the number of points - 1,
  164.  * including the endpoints); we require k >= 1.
  165.  * If k or any of the coefficient values are too large,
  166.  * use recursive subdivision to whittle them down.
  167.  */
  168. int
  169. gx_flatten_sample(gx_path *ppath, int k, curve_segment *pc,
  170.   segment_notes notes)
  171. {    fixed x0, y0;
  172.     /* x1 ... y3 were defined above */
  173.     fixed cx, bx, ax, cy, by, ay;
  174.     fixed ptx, pty;
  175.     fixed x, y;
  176.     /*
  177.      * We can compute successive values by finite differences,
  178.      * using the formulas:
  179.         x(t) =
  180.           a*t^3 + b*t^2 + c*t + d =>
  181.         dx(t) = x(t+e)-x(t) =
  182.           a*(3*t^2*e + 3*t*e^2 + e^3) + b*(2*t*e + e^2) + c*e =
  183.           (3*a*e)*t^2 + (3*a*e^2 + 2*b*e)*t + (a*e^3 + b*e^2 + c*e) =>
  184.         d2x(t) = dx(t+e)-dx(t) =
  185.           (3*a*e)*(2*t*e + e^2) + (3*a*e^2 + 2*b*e)*e =
  186.           (6*a*e^2)*t + (6*a*e^3 + 2*b*e^2) =>
  187.         d3x(t) = d2x(t+e)-d2x(t) =
  188.           6*a*e^3;
  189.         x(0) = d, dx(0) = (a*e^3 + b*e^2 + c*e),
  190.           d2x(0) = 6*a*e^3 + 2*b*e^2;
  191.      * In these formulas, e = 1/2^k; of course, there are separate
  192.      * computations for the x and y values.
  193.      *
  194.      * There is a tradeoff in doing the above computation in fixed
  195.      * point.  If we separate out the constant term (d) and require that
  196.      * all the other values fit in a long, then on a 32-bit machine with
  197.      * 12 bits of fraction in a fixed, k = 4 implies a maximum curve
  198.      * size of 128 pixels; anything larger requires subdividing the
  199.      * curve.  On the other hand, doing the computations in explicit
  200.      * double precision slows down the loop by a factor of 3 or so.  We
  201.      * found to our surprise that the latter is actually faster, because
  202.      * the additional subdivisions cost more than the slower loop.
  203.      *
  204.      * We represent each quantity as I+R/M, where I is an "integer" and
  205.      * the "remainder" R lies in the range 0 <= R < M=2^(3*k).  Note
  206.      * that R may temporarily exceed M; for this reason, we require that
  207.      * M have at least one free high-order bit.  To reduce the number of
  208.      * variables, we don't actually compute M, only M-1 (rmask).  */
  209.     uint i;
  210.     uint rmask;            /* M-1 */
  211.     fixed idx, idy, id2x, id2y, id3x, id3y;        /* I */
  212.     uint rx, ry, rdx, rdy, rd2x, rd2y, rd3x, rd3y;    /* R */
  213.     gs_fixed_point _ss *ppt;
  214. #define max_points 50            /* arbitrary */
  215.     gs_fixed_point points[max_points + 1];
  216.  
  217. top:    x0 = ppath->position.x;
  218.     y0 = ppath->position.y;
  219. #ifdef DEBUG
  220.     if ( gs_debug_c('3') )
  221.       dprintf4("[3]x0=%f y0=%f x1=%f y1=%f\n",
  222.            fixed2float(x0), fixed2float(y0),
  223.            fixed2float(x1), fixed2float(y1)),
  224.       dprintf5("   x2=%f y2=%f x3=%f y3=%f  k=%d\n",
  225.            fixed2float(x2), fixed2float(y2),
  226.            fixed2float(x3), fixed2float(y3), k);
  227. #endif
  228.     {    fixed x01, x12, y01, y12;
  229.         curve_points_to_coefficients(x0, x1, x2, x3, ax, bx, cx,
  230.                          x01, x12);
  231.         curve_points_to_coefficients(y0, y1, y2, y3, ay, by, cy,
  232.                          y01, y12);
  233.     }
  234.  
  235.     if_debug6('3', "[3]ax=%f bx=%f cx=%f\n   ay=%f by=%f cy=%f\n",
  236.           fixed2float(ax), fixed2float(bx), fixed2float(cx),
  237.           fixed2float(ay), fixed2float(by), fixed2float(cy));
  238. #define max_fast (max_fixed / 6)
  239. #define min_fast (-max_fast)
  240. #define in_range(v) (v < max_fast && v > min_fast)
  241.     if ( k == 0 )
  242.     {    /* The curve is very short, or anomalous in some way. */
  243.         /* Just add a line and exit. */
  244.         return gx_path_add_line_notes(ppath, x3, y3, notes);
  245.     }
  246.     if ( k <= k_sample_max &&
  247.          in_range(ax) && in_range(ay) &&
  248.          in_range(bx) && in_range(by) &&
  249.          in_range(cx) && in_range(cy)
  250.        )
  251.     {    x = x0, y = y0;
  252.         rx = ry = 0;
  253.         ppt = points;
  254.         /* Fast check for n == 3, a common special case */
  255.         /* for small characters. */
  256.         if ( k == 1 )
  257.         {
  258. #define poly2(a,b,c)\
  259.   arith_rshift_1(arith_rshift_1(arith_rshift_1(a) + b) + c)
  260.             x += poly2(ax, bx, cx);
  261.             y += poly2(ay, by, cy);
  262. #undef poly2
  263.             if_debug2('3', "[3]dx=%f, dy=%f\n",
  264.                   fixed2float(x - x0), fixed2float(y - y0));
  265.             if_debug3('3', "[3]%s x=%g, y=%g\n",
  266.                   (((x ^ x0) | (y ^ y0)) & float2fixed(-0.5) ?
  267.                    "add" : "skip"),
  268.                   fixed2float(x), fixed2float(y));
  269.             if ( ((x ^ x0) | (y ^ y0)) & float2fixed(-0.5) )
  270.               ppt->x = ptx = x,
  271.               ppt->y = pty = y,
  272.               ppt++;
  273.             goto last;
  274.         }
  275.         else
  276.         {    fixed bx2 = bx << 1, by2 = by << 1;
  277.             fixed ax6 = ((ax << 1) + ax) << 1,
  278.                   ay6 = ((ay << 1) + ay) << 1;
  279. #define adjust_rem(r, q)\
  280.   if ( r > rmask ) q ++, r &= rmask
  281.             const int k2 = k << 1;
  282.             const int k3 = k2 + k;
  283.             rmask = (1 << k3) - 1;
  284.             /* We can compute all the remainders as ints, */
  285.             /* because we know they don't exceed M. */
  286.             /* cx/y terms */
  287.             idx = arith_rshift(cx, k),
  288.               idy = arith_rshift(cy, k);
  289.             rdx = ((uint)cx << k2) & rmask,
  290.               rdy = ((uint)cy << k2) & rmask;
  291.             /* bx/y terms */
  292.             id2x = arith_rshift(bx2, k2),
  293.               id2y = arith_rshift(by2, k2);
  294.             rd2x = ((uint)bx2 << k) & rmask,
  295.               rd2y = ((uint)by2 << k) & rmask;
  296.             idx += arith_rshift_1(id2x),
  297.               idy += arith_rshift_1(id2y);
  298.             rdx += ((uint)bx << k) & rmask,
  299.               rdy += ((uint)by << k) & rmask;
  300.             adjust_rem(rdx, idx);
  301.             adjust_rem(rdy, idy);
  302.             /* ax/y terms */
  303.             idx += arith_rshift(ax, k3),
  304.               idy += arith_rshift(ay, k3);
  305.             rdx += (uint)ax & rmask,
  306.               rdy += (uint)ay & rmask;
  307.             adjust_rem(rdx, idx);
  308.             adjust_rem(rdy, idy);
  309.             id2x += id3x = arith_rshift(ax6, k3),
  310.               id2y += id3y = arith_rshift(ay6, k3);
  311.             rd2x += rd3x = (uint)ax6 & rmask,
  312.               rd2y += rd3y = (uint)ay6 & rmask;
  313.             adjust_rem(rd2x, id2x);
  314.             adjust_rem(rd2y, id2y);
  315. #undef adjust_rem
  316.         }
  317.     }
  318.     else
  319.     {    /*
  320.          * Curve is too long.  Break into two pieces and recur.
  321.          */
  322.         curve_segment cseg;
  323.         int code;
  324.  
  325.         k--;
  326.         split_curve_midpoint(x0, y0, pc, &cseg, pc);
  327.         code = gx_flatten_sample(ppath, k, &cseg, notes);
  328.         if ( code < 0 )
  329.           return code;
  330.         notes |= sn_not_first;
  331.         goto top;
  332.     }
  333.     if_debug1('2', "[2]sampling k=%d\n", k);
  334.     ptx = x0, pty = y0;
  335.     for ( i = (1 << k) - 1; ; )
  336.     {    int code;
  337. #ifdef DEBUG
  338.         if ( gs_debug_c('3') )
  339.           dprintf4("[3]dx=%f+%d, dy=%f+%d\n",
  340.                fixed2float(idx), rdx,
  341.                fixed2float(idy), rdy),
  342.           dprintf4("   d2x=%f+%d, d2y=%f+%d\n",
  343.                fixed2float(id2x), rd2x,
  344.                fixed2float(id2y), rd2y),
  345.           dprintf4("   d3x=%f+%d, d3y=%f+%d\n",
  346.                fixed2float(id3x), rd3x,
  347.                fixed2float(id3y), rd3y);
  348. #endif
  349. #define accum(i, r, di, dr)\
  350.   if ( (r += dr) > rmask ) r &= rmask, i += di + 1;\
  351.   else i += di
  352.         accum(x, rx, idx, rdx);
  353.         accum(y, ry, idy, rdy);
  354.         if_debug3('3', "[3]%s x=%g, y=%g\n",
  355.               (((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) ?
  356.                "add" : "skip"),
  357.               fixed2float(x), fixed2float(y));
  358.         /*
  359.          * Skip very short segments -- those that lie entirely within
  360.          * a square half-pixel.  Also merge nearly collinear
  361.          * segments -- those where one coordinate of all three points
  362.          * (the two endpoints and the midpoint) lie within the same
  363.          * half-pixel and both coordinates are monotonic.
  364.          * Note that ptx/y, the midpoint, is the same as ppt[-1].x/y;
  365.          * the previous point is ppt[-2].x/y.
  366.          */
  367. #define coord_near(v, ptv)\
  368.   (!( ((v) ^ (ptv)) & float2fixed(-0.5) ))
  369. #define coords_in_order(v0, v1, v2)\
  370.   ( (((v1) - (v0)) ^ ((v2) - (v1))) >= 0 )
  371.         if ( coord_near(x, ptx) )
  372.           {    /* X coordinates are within a half-pixel. */
  373.             if ( coord_near(y, pty) )
  374.               goto skip;        /* short segment */
  375. #if MERGE_COLLINEAR_SEGMENTS
  376.             /* Check for collinear segments. */
  377.             if ( ppt > points + 1 && coord_near(x, ppt[-2].x) &&
  378.                  coords_in_order(ppt[-2].x, ptx, x) &&
  379.                  coords_in_order(ppt[-2].y, pty, y)
  380.                )
  381.               --ppt;    /* remove middle point */
  382. #endif
  383.           }
  384.         else if ( coord_near(y, pty) )
  385.           {    /* Y coordinates are within a half-pixel. */
  386. #if MERGE_COLLINEAR_SEGMENTS
  387.             /* Check for collinear segments. */
  388.             if ( ppt > points + 1 && coord_near(y, ppt[-2].y) &&
  389.                  coords_in_order(ppt[-2].x, ptx, x) &&
  390.                  coords_in_order(ppt[-2].y, pty, y)
  391.                )
  392.               --ppt;    /* remove middle point */
  393. #endif
  394.           }
  395. #undef coord_near
  396. #undef coords_in_order
  397.         /* Add a line. */
  398.         if ( ppt == &points[max_points] )
  399.           { if ( notes & sn_not_first )
  400.               code = gx_path_add_lines_notes(ppath, points, max_points,
  401.                              notes);
  402.             else
  403.               { code = gx_path_add_line_notes(ppath, points[0].x,
  404.                               points[0].y, notes);
  405.                 if ( code < 0 )
  406.               return code;
  407.             code = gx_path_add_lines_notes(ppath, points,
  408.                     max_points - 1, notes | sn_not_first);
  409.               }
  410.             if ( code < 0 )
  411.               return code;
  412.             ppt = points;
  413.             notes |= sn_not_first;
  414.           }
  415.         ppt->x = ptx = x;
  416.         ppt->y = pty = y;
  417.         ppt++;
  418. skip:        if ( --i == 0 )
  419.           break;        /* don't bother with last accum */
  420.         accum(idx, rdx, id2x, rd2x);
  421.         accum(id2x, rd2x, id3x, rd3x);
  422.         accum(idy, rdy, id2y, rd2y);
  423.         accum(id2y, rd2y, id3y, rd3y);
  424. #undef accum
  425.     }
  426. last:    if_debug2('3', "[3]last x=%g, y=%g\n",
  427.           fixed2float(x3), fixed2float(y3));
  428.     if ( ppt > points )
  429.       { int count = ppt + 1 - points;
  430.         gs_fixed_point _ss *pts = points;
  431.  
  432.         if ( !(notes & sn_not_first) )
  433.           { int code = gx_path_add_line_notes(ppath,
  434.                           points[0].x, points[0].y,
  435.                           notes);
  436.             if ( code < 0 )
  437.           return code;
  438.         ++pts, --count;
  439.         notes |= sn_not_first;
  440.           }
  441.         ppt->x = x3, ppt->y = y3;
  442.         return gx_path_add_lines_notes(ppath, pts, count, notes);
  443.       }
  444.     return gx_path_add_line_notes(ppath, x3, y3, notes);
  445. }
  446.  
  447. #undef x1
  448. #undef y1
  449. #undef x2
  450. #undef y2
  451. #undef x3
  452. #undef y3
  453.