1 /* Copyright (C) 1992, 2000 artofcode LLC. All rights reserved.
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2 of the License, or (at your
6 option) any later version.
7
8 This program is distributed in the hope that it will be useful, but
9 WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 General Public License for more details.
12
13 You should have received a copy of the GNU General Public License along
14 with this program; if not, write to the Free Software Foundation, Inc.,
15 59 Temple Place, Suite 330, Boston, MA, 02111-1307.
16
17 */
18
19 /*$Id: gxpcopy.c,v 1.3.6.1.2.1 2003/01/17 00:49:04 giles Exp $ */
20 /* Path copying and flattening */
21 #include "math_.h"
22 #include "gx.h"
23 #include "gserrors.h"
24 #include "gconfigv.h" /* for USE_FPU */
25 #include "gxfixed.h"
26 #include "gxfarith.h"
27 #include "gxistate.h" /* for access to line params */
28 #include "gzpath.h"
29
30 /* Forward declarations */
31 private void adjust_point_to_tangent(P3(segment *, const segment *,
32 const gs_fixed_point *));
33 private int monotonize_internal(P2(gx_path *, const curve_segment *));
34
35 /* Copy a path, optionally flattening or monotonizing it. */
36 /* If the copy fails, free the new path. */
37 int
gx_path_copy_reducing(const gx_path * ppath_old,gx_path * ppath,fixed fixed_flatness,const gs_imager_state * pis,gx_path_copy_options options)38 gx_path_copy_reducing(const gx_path *ppath_old, gx_path *ppath,
39 fixed fixed_flatness, const gs_imager_state *pis,
40 gx_path_copy_options options)
41 {
42 const segment *pseg;
43 fixed flat = fixed_flatness;
44 gs_fixed_point expansion;
45 /*
46 * Since we're going to be adding to the path, unshare it
47 * before we start.
48 */
49 int code = gx_path_unshare(ppath);
50
51 if (code < 0)
52 return code;
53 #ifdef DEBUG
54 if (gs_debug_c('P'))
55 gx_dump_path(ppath_old, "before reducing");
56 #endif
57 if (options & pco_for_stroke) {
58 /* Precompute the maximum expansion of the bounding box. */
59 double width = pis->line_params.half_width;
60
61 expansion.x =
62 float2fixed((fabs(pis->ctm.xx) + fabs(pis->ctm.yx)) * width) * 2;
63 expansion.y =
64 float2fixed((fabs(pis->ctm.xy) + fabs(pis->ctm.yy)) * width) * 2;
65 }
66 pseg = (const segment *)(ppath_old->first_subpath);
67 while (pseg) {
68 switch (pseg->type) {
69 case s_start:
70 code = gx_path_add_point(ppath,
71 pseg->pt.x, pseg->pt.y);
72 break;
73 case s_curve:
74 {
75 const curve_segment *pc = (const curve_segment *)pseg;
76
77 if (fixed_flatness == max_fixed) { /* don't flatten */
78 if (options & pco_monotonize)
79 code = monotonize_internal(ppath, pc);
80 else
81 code = gx_path_add_curve_notes(ppath,
82 pc->p1.x, pc->p1.y, pc->p2.x, pc->p2.y,
83 pc->pt.x, pc->pt.y, pseg->notes);
84 } else {
85 fixed x0 = ppath->position.x;
86 fixed y0 = ppath->position.y;
87 segment_notes notes = pseg->notes;
88 curve_segment cseg;
89 int k;
90
91 if (options & pco_for_stroke) {
92 /*
93 * When flattening for stroking, the flatness
94 * must apply to the outside of the resulting
95 * stroked region. We approximate this by
96 * dividing the flatness by the ratio of the
97 * expanded bounding box to the original
98 * bounding box. This is crude, but pretty
99 * simple to calculate, and produces reasonably
100 * good results.
101 */
102 fixed min01, max01, min23, max23;
103 fixed ex, ey, flat_x, flat_y;
104
105 #define SET_EXTENT(r, c0, c1, c2, c3)\
106 BEGIN\
107 if (c0 < c1) min01 = c0, max01 = c1;\
108 else min01 = c1, max01 = c0;\
109 if (c2 < c3) min23 = c2, max23 = c3;\
110 else min23 = c3, max23 = c2;\
111 r = max(max01, max23) - min(min01, min23);\
112 END
113 SET_EXTENT(ex, x0, pc->p1.x, pc->p2.x, pc->pt.x);
114 SET_EXTENT(ey, y0, pc->p1.y, pc->p2.y, pc->pt.y);
115 #undef SET_EXTENT
116 /*
117 * We check for the degenerate case specially
118 * to avoid a division by zero.
119 */
120 if (ex == 0 || ey == 0)
121 flat = 0;
122 else {
123 flat_x =
124 fixed_mult_quo(fixed_flatness, ex,
125 ex + expansion.x);
126 flat_y =
127 fixed_mult_quo(fixed_flatness, ey,
128 ey + expansion.y);
129 flat = min(flat_x, flat_y);
130 }
131 }
132 k = gx_curve_log2_samples(x0, y0, pc, flat);
133 if (options & pco_accurate) {
134 segment *start;
135 segment *end;
136
137 /*
138 * Add an extra line, which will become
139 * the tangent segment.
140 */
141 code = gx_path_add_line_notes(ppath, x0, y0,
142 notes);
143 if (code < 0)
144 break;
145 start = ppath->current_subpath->last;
146 notes |= sn_not_first;
147 cseg = *pc;
148 code = gx_flatten_sample(ppath, k, &cseg, notes);
149 if (code < 0)
150 break;
151 /*
152 * Adjust the first and last segments so that
153 * they line up with the tangents.
154 */
155 end = ppath->current_subpath->last;
156 if ((code = gx_path_add_line_notes(ppath,
157 ppath->position.x,
158 ppath->position.y,
159 pseg->notes | sn_not_first)) < 0)
160 break;
161 adjust_point_to_tangent(start, start->next,
162 &pc->p1);
163 adjust_point_to_tangent(end, end->prev,
164 &pc->p2);
165 } else {
166 cseg = *pc;
167 code = gx_flatten_sample(ppath, k, &cseg, notes);
168 }
169 }
170 break;
171 }
172 case s_line:
173 code = gx_path_add_line_notes(ppath,
174 pseg->pt.x, pseg->pt.y, pseg->notes);
175 break;
176 case s_line_close:
177 code = gx_path_close_subpath(ppath);
178 break;
179 default: /* can't happen */
180 code = gs_note_error(gs_error_unregistered);
181 }
182 if (code < 0) {
183 gx_path_new(ppath);
184 return code;
185 }
186 pseg = pseg->next;
187 }
188 if (path_last_is_moveto(ppath_old))
189 gx_path_add_point(ppath, ppath_old->position.x,
190 ppath_old->position.y);
191 if (ppath_old->bbox_set) {
192 if (ppath->bbox_set) {
193 ppath->bbox.p.x = min(ppath_old->bbox.p.x, ppath->bbox.p.x);
194 ppath->bbox.p.y = min(ppath_old->bbox.p.y, ppath->bbox.p.y);
195 ppath->bbox.q.x = max(ppath_old->bbox.q.x, ppath->bbox.q.x);
196 ppath->bbox.q.y = max(ppath_old->bbox.q.y, ppath->bbox.q.y);
197 } else {
198 ppath->bbox_set = true;
199 ppath->bbox = ppath_old->bbox;
200 }
201 }
202 #ifdef DEBUG
203 if (gs_debug_c('P'))
204 gx_dump_path(ppath, "after reducing");
205 #endif
206 return 0;
207 }
208
209 /*
210 * Adjust one end of a line (the first or last line of a flattened curve)
211 * so it falls on the curve tangent. The closest point on the line from
212 * (0,0) to (C,D) to a point (U,V) -- i.e., the point on the line at which
213 * a perpendicular line from the point intersects it -- is given by
214 * T = (C*U + D*V) / (C^2 + D^2)
215 * (X,Y) = (C*T,D*T)
216 * However, any smaller value of T will also work: the one we actually
217 * use is 0.25 * the value we just derived. We must check that
218 * numerical instabilities don't lead to a negative value of T.
219 */
220 private void
adjust_point_to_tangent(segment * pseg,const segment * next,const gs_fixed_point * p1)221 adjust_point_to_tangent(segment * pseg, const segment * next,
222 const gs_fixed_point * p1)
223 {
224 const fixed x0 = pseg->pt.x, y0 = pseg->pt.y;
225 const fixed fC = p1->x - x0, fD = p1->y - y0;
226
227 /*
228 * By far the commonest case is that the end of the curve is
229 * horizontal or vertical. Check for this specially, because
230 * we can handle it with far less work (and no floating point).
231 */
232 if (fC == 0) {
233 /* Vertical tangent. */
234 const fixed DT = arith_rshift(next->pt.y - y0, 2);
235
236 if (fD == 0)
237 return; /* anomalous case */
238 if_debug1('2', "[2]adjusting vertical: DT = %g\n",
239 fixed2float(DT));
240 if ((DT ^ fD) > 0)
241 pseg->pt.y = DT + y0;
242 } else if (fD == 0) {
243 /* Horizontal tangent. */
244 const fixed CT = arith_rshift(next->pt.x - x0, 2);
245
246 if_debug1('2', "[2]adjusting horizontal: CT = %g\n",
247 fixed2float(CT));
248 if ((CT ^ fC) > 0)
249 pseg->pt.x = CT + x0;
250 } else {
251 /* General case. */
252 const double C = fC, D = fD;
253 const double T = (C * (next->pt.x - x0) + D * (next->pt.y - y0)) /
254 (C * C + D * D);
255
256 if_debug3('2', "[2]adjusting: C = %g, D = %g, T = %g\n",
257 C, D, T);
258 if (T > 0) {
259 pseg->pt.x = arith_rshift((fixed) (C * T), 2) + x0;
260 pseg->pt.y = arith_rshift((fixed) (D * T), 2) + y0;
261 }
262 }
263 }
264
265 /* ---------------- Curve flattening ---------------- */
266
267 #define x1 pc->p1.x
268 #define y1 pc->p1.y
269 #define x2 pc->p2.x
270 #define y2 pc->p2.y
271 #define x3 pc->pt.x
272 #define y3 pc->pt.y
273
274 #ifdef DEBUG
275 private void
dprint_curve(const char * str,fixed x0,fixed y0,const curve_segment * pc)276 dprint_curve(const char *str, fixed x0, fixed y0, const curve_segment * pc)
277 {
278 dlprintf9("%s p0=(%g,%g) p1=(%g,%g) p2=(%g,%g) p3=(%g,%g)\n",
279 str, fixed2float(x0), fixed2float(y0),
280 fixed2float(pc->p1.x), fixed2float(pc->p1.y),
281 fixed2float(pc->p2.x), fixed2float(pc->p2.y),
282 fixed2float(pc->pt.x), fixed2float(pc->pt.y));
283 }
284 #endif
285
286 /* Initialize a cursor for rasterizing a monotonic curve. */
287 void
gx_curve_cursor_init(curve_cursor * prc,fixed x0,fixed y0,const curve_segment * pc,int k)288 gx_curve_cursor_init(curve_cursor * prc, fixed x0, fixed y0,
289 const curve_segment * pc, int k)
290 {
291 fixed v01, v12;
292 int k2 = k + k, k3 = k2 + k;
293
294 #define bits_fit(v, n)\
295 (any_abs(v) <= max_fixed >> (n))
296 /* The +2s are because of t3d and t2d, see below. */
297 #define coeffs_fit(a, b, c)\
298 (k3 <= sizeof(fixed) * 8 - 3 &&\
299 bits_fit(a, k3 + 2) && bits_fit(b, k2 + 2) && bits_fit(c, k + 1))
300
301 prc->k = k;
302 prc->p0.x = x0, prc->p0.y = y0;
303 prc->pc = pc;
304 /* Compute prc->a..c taking into account reversal of xy0/3 */
305 /* in curve_x_at_y. */
306 {
307 fixed w0, w1, w2, w3;
308
309 if (y0 < y3)
310 w0 = x0, w1 = x1, w2 = x2, w3 = x3;
311 else
312 w0 = x3, w1 = x2, w2 = x1, w3 = x0;
313 curve_points_to_coefficients(w0, w1, w2, w3,
314 prc->a, prc->b, prc->c, v01, v12);
315 }
316 prc->double_set = false;
317 prc->fixed_limit =
318 (coeffs_fit(prc->a, prc->b, prc->c) ? (1 << k) - 1 : -1);
319 /* Initialize the cache. */
320 prc->cache.ky0 = prc->cache.ky3 = y0;
321 prc->cache.xl = x0;
322 prc->cache.xd = 0;
323 }
324
325 /*
326 * Determine the X value on a monotonic curve at a given Y value. It turns
327 * out that we use so few points on the curve that it's actually faster to
328 * locate the desired point by recursive subdivision each time than to try
329 * to keep a cursor that we move incrementally. What's even more surprising
330 * is that if floating point arithmetic is reasonably fast, it's faster to
331 * compute the X value at the desired point explicitly than to do the
332 * recursive subdivision on X as well as Y.
333 */
334 #define SUBDIVIDE_X USE_FPU_FIXED
335 fixed
gx_curve_x_at_y(curve_cursor * prc,fixed y)336 gx_curve_x_at_y(curve_cursor * prc, fixed y)
337 {
338 fixed xl, xd;
339 fixed yd, yrel;
340
341 /* Check the cache before doing anything else. */
342 if (y >= prc->cache.ky0 && y <= prc->cache.ky3) {
343 yd = prc->cache.ky3 - prc->cache.ky0;
344 yrel = y - prc->cache.ky0;
345 xl = prc->cache.xl;
346 xd = prc->cache.xd;
347 goto done;
348 } {
349 #define x0 prc->p0.x
350 #define y0 prc->p0.y
351 const curve_segment *pc = prc->pc;
352
353 /* Reduce case testing by ensuring y3 >= y0. */
354 fixed cy0 = y0, cy1, cy2, cy3 = y3;
355 fixed cx0;
356
357 #if SUBDIVIDE_X
358 fixed cx1, cx2, cx3;
359
360 #else
361 int t = 0;
362
363 #endif
364 int k, i;
365
366 if (cy0 > cy3)
367 cx0 = x3,
368 #if SUBDIVIDE_X
369 cx1 = x2, cx2 = x1, cx3 = x0,
370 #endif
371 cy0 = y3, cy1 = y2, cy2 = y1, cy3 = y0;
372 else
373 cx0 = x0,
374 #if SUBDIVIDE_X
375 cx1 = x1, cx2 = x2, cx3 = x3,
376 #endif
377 cy1 = y1, cy2 = y2;
378 #define midpoint_fast(a,b)\
379 arith_rshift_1((a) + (b) + 1)
380 for (i = k = prc->k; i > 0; --i) {
381 fixed ym = midpoint_fast(cy1, cy2);
382 fixed yn = ym + arith_rshift(cy0 - cy1 - cy2 + cy3 + 4, 3);
383
384 #if SUBDIVIDE_X
385 fixed xm = midpoint_fast(cx1, cx2);
386 fixed xn = xm + arith_rshift(cx0 - cx1 - cx2 + cx3 + 4, 3);
387
388 #else
389 t <<= 1;
390 #endif
391
392 if (y < yn)
393 #if SUBDIVIDE_X
394 cx1 = midpoint_fast(cx0, cx1),
395 cx2 = midpoint_fast(cx1, xm),
396 cx3 = xn,
397 #endif
398 cy1 = midpoint_fast(cy0, cy1),
399 cy2 = midpoint_fast(cy1, ym),
400 cy3 = yn;
401 else
402 #if SUBDIVIDE_X
403 cx2 = midpoint_fast(cx2, cx3),
404 cx1 = midpoint_fast(xm, cx2),
405 cx0 = xn,
406 #else
407 t++,
408 #endif
409 cy2 = midpoint_fast(cy2, cy3),
410 cy1 = midpoint_fast(ym, cy2),
411 cy0 = yn;
412 }
413 #if SUBDIVIDE_X
414 xl = cx0;
415 xd = cx3 - cx0;
416 #else
417 {
418 fixed a = prc->a, b = prc->b, c = prc->c;
419
420 /* We must use (1 << k) >> 1 instead of 1 << (k - 1) in case k == 0. */
421 #define compute_fixed(a, b, c)\
422 arith_rshift(arith_rshift(arith_rshift(a * t3, k) + b * t2, k)\
423 + c * t + ((1 << k) >> 1), k)
424 #define compute_diff_fixed(a, b, c)\
425 arith_rshift(arith_rshift(arith_rshift(a * t3d, k) + b * t2d, k)\
426 + c, k)
427
428 /* use multiply if possible */
429 #define np2(n) (1.0 / (1L << (n)))
430 static const double k_denom[11] =
431 {np2(0), np2(1), np2(2), np2(3), np2(4),
432 np2(5), np2(6), np2(7), np2(8), np2(9), np2(10)
433 };
434 static const double k2_denom[11] =
435 {np2(0), np2(2), np2(4), np2(6), np2(8),
436 np2(10), np2(12), np2(14), np2(16), np2(18), np2(20)
437 };
438 static const double k3_denom[11] =
439 {np2(0), np2(3), np2(6), np2(9), np2(12),
440 np2(15), np2(18), np2(21), np2(24), np2(27), np2(30)
441 };
442 double den1, den2;
443
444 #undef np2
445
446 #define setup_floating(da, db, dc, a, b, c)\
447 (k >= countof(k_denom) ?\
448 (den1 = ldexp(1.0, -k),\
449 den2 = den1 * den1,\
450 da = (den2 * den1) * a,\
451 db = den2 * b,\
452 dc = den1 * c) :\
453 (da = k3_denom[k] * a,\
454 db = k2_denom[k] * b,\
455 dc = k_denom[k] * c))
456 #define compute_floating(da, db, dc)\
457 ((fixed)(da * t3 + db * t2 + dc * t + 0.5))
458 #define compute_diff_floating(da, db, dc)\
459 ((fixed)(da * t3d + db * t2d + dc))
460
461 if (t <= prc->fixed_limit) { /* We can do everything in integer/fixed point. */
462 int t2 = t * t, t3 = t2 * t;
463 int t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
464
465 xl = compute_fixed(a, b, c) + cx0;
466 xd = compute_diff_fixed(a, b, c);
467 #ifdef DEBUG
468 {
469 double fa, fb, fc;
470 fixed xlf, xdf;
471
472 setup_floating(fa, fb, fc, a, b, c);
473 xlf = compute_floating(fa, fb, fc) + cx0;
474 xdf = compute_diff_floating(fa, fb, fc);
475 if (any_abs(xlf - xl) > fixed_epsilon ||
476 any_abs(xdf - xd) > fixed_epsilon
477 )
478 dlprintf9("Curve points differ: k=%d t=%d a,b,c=%g,%g,%g\n xl,xd fixed=%g,%g floating=%g,%g\n",
479 k, t,
480 fixed2float(a), fixed2float(b), fixed2float(c),
481 fixed2float(xl), fixed2float(xd),
482 fixed2float(xlf), fixed2float(xdf));
483 /*xl = xlf, xd = xdf; */
484 }
485 #endif
486 } else { /*
487 * Either t3 (and maybe t2) won't fit in an int, or more
488 * likely the result of the multiplies won't fit.
489 */
490 #define fa prc->da
491 #define fb prc->db
492 #define fc prc->dc
493 if (!prc->double_set) {
494 setup_floating(fa, fb, fc, a, b, c);
495 prc->double_set = true;
496 }
497 if (t < 1L << ((sizeof(long) * 8 - 1) / 3)) { /*
498 * t3 (and maybe t2) might not fit in an int, but they
499 * will fit in a long. If we have slow floating point,
500 * do the computation in double-precision fixed point,
501 * otherwise do it in fixed point.
502 */
503 long t2 = (long)t * t, t3 = t2 * t;
504 long t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
505
506 xl = compute_floating(fa, fb, fc) + cx0;
507 xd = compute_diff_floating(fa, fb, fc);
508 } else { /*
509 * t3 (and maybe t2) don't even fit in a long.
510 * Do the entire computation in floating point.
511 */
512 double t2 = (double)t * t, t3 = t2 * t;
513 double t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
514
515 xl = compute_floating(fa, fb, fc) + cx0;
516 xd = compute_diff_floating(fa, fb, fc);
517 }
518 #undef fa
519 #undef fb
520 #undef fc
521 }
522 }
523 #endif /* (!)SUBDIVIDE_X */
524
525 /* Update the cache. */
526 prc->cache.ky0 = cy0;
527 prc->cache.ky3 = cy3;
528 prc->cache.xl = xl;
529 prc->cache.xd = xd;
530 yd = cy3 - cy0;
531 yrel = y - cy0;
532 #undef x0
533 #undef y0
534 }
535 done:
536 /*
537 * Now interpolate linearly between current and next.
538 * We know that 0 <= yrel < yd.
539 * It's unlikely but possible that cy0 = y = cy3:
540 * handle this case specially.
541 */
542 if (yrel == 0)
543 return xl;
544 /*
545 * Compute in fixed point if possible.
546 */
547 #define HALF_FIXED_BITS ((fixed)1 << (sizeof(fixed) * 4))
548 if (yrel < HALF_FIXED_BITS) {
549 if (xd >= 0) {
550 if (xd < HALF_FIXED_BITS)
551 return (ufixed)xd * (ufixed)yrel / (ufixed)yd + xl;
552 } else {
553 if (xd > -HALF_FIXED_BITS) {
554 /* Be careful to take the floor of the result. */
555 ufixed num = (ufixed)(-xd) * (ufixed)yrel;
556 ufixed quo = num / (ufixed)yd;
557
558 if (quo * (ufixed)yd != num)
559 quo += fixed_epsilon;
560 return xl - (fixed)quo;
561 }
562 }
563 }
564 #undef HALF_FIXED_BITS
565 return fixed_mult_quo(xd, yrel, yd) + xl;
566 }
567
568 #undef x1
569 #undef y1
570 #undef x2
571 #undef y2
572 #undef x3
573 #undef y3
574
575 /* ---------------- Monotonic curves ---------------- */
576
577 /* Test whether a path is free of non-monotonic curves. */
578 bool
gx_path_is_monotonic(const gx_path * ppath)579 gx_path_is_monotonic(const gx_path * ppath)
580 {
581 const segment *pseg = (const segment *)(ppath->first_subpath);
582 gs_fixed_point pt0;
583
584 while (pseg) {
585 switch (pseg->type) {
586 case s_start:
587 {
588 const subpath *psub = (const subpath *)pseg;
589
590 /* Skip subpaths without curves. */
591 if (!psub->curve_count)
592 pseg = psub->last;
593 }
594 break;
595 case s_curve:
596 {
597 const curve_segment *pc = (const curve_segment *)pseg;
598 double t[2];
599 int nz = gx_curve_monotonic_points(pt0.y,
600 pc->p1.y, pc->p2.y, pc->pt.y, t);
601
602 if (nz != 0)
603 return false;
604 nz = gx_curve_monotonic_points(pt0.x,
605 pc->p1.x, pc->p2.x, pc->pt.x, t);
606 if (nz != 0)
607 return false;
608 }
609 break;
610 default:
611 ;
612 }
613 pt0 = pseg->pt;
614 pseg = pseg->next;
615 }
616 return true;
617 }
618
619 /* Monotonize a curve, by splitting it if necessary. */
620 /* In the worst case, this could split the curve into 9 pieces. */
621 private int
monotonize_internal(gx_path * ppath,const curve_segment * pc)622 monotonize_internal(gx_path * ppath, const curve_segment * pc)
623 {
624 fixed x0 = ppath->position.x, y0 = ppath->position.y;
625 segment_notes notes = pc->notes;
626 double t[2];
627
628 #define max_segs 9
629 curve_segment cs[max_segs];
630 const curve_segment *pcs;
631 curve_segment *pcd;
632 int i, j, nseg;
633 int nz;
634
635 /* Monotonize in Y. */
636 nz = gx_curve_monotonic_points(y0, pc->p1.y, pc->p2.y, pc->pt.y, t);
637 nseg = max_segs - 1 - nz;
638 pcd = cs + nseg;
639 if (nz == 0)
640 *pcd = *pc;
641 else {
642 gx_curve_split(x0, y0, pc, t[0], pcd, pcd + 1);
643 if (nz == 2)
644 gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
645 (t[1] - t[0]) / (1 - t[0]),
646 pcd + 1, pcd + 2);
647 }
648
649 /* Monotonize in X. */
650 for (pcs = pcd, pcd = cs, j = nseg; j < max_segs; ++pcs, ++j) {
651 nz = gx_curve_monotonic_points(x0, pcs->p1.x, pcs->p2.x,
652 pcs->pt.x, t);
653
654 if (nz == 0)
655 *pcd = *pcs;
656 else {
657 gx_curve_split(x0, y0, pcs, t[0], pcd, pcd + 1);
658 if (nz == 2)
659 gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
660 (t[1] - t[0]) / (1 - t[0]),
661 pcd + 1, pcd + 2);
662 }
663 pcd += nz + 1;
664 x0 = pcd[-1].pt.x;
665 y0 = pcd[-1].pt.y;
666 }
667 nseg = pcd - cs;
668
669 /* Add the segment(s) to the output. */
670 #ifdef DEBUG
671 if (gs_debug_c('2')) {
672 int pi;
673 gs_fixed_point pp0;
674
675 pp0 = ppath->position;
676 if (nseg == 1)
677 dprint_curve("[2]No split", pp0.x, pp0.y, pc);
678 else {
679 dlprintf1("[2]Split into %d segments:\n", nseg);
680 dprint_curve("[2]Original", pp0.x, pp0.y, pc);
681 for (pi = 0; pi < nseg; ++pi) {
682 dprint_curve("[2] =>", pp0.x, pp0.y, cs + pi);
683 pp0 = cs[pi].pt;
684 }
685 }
686 }
687 #endif
688 for (pcs = cs, i = 0; i < nseg; ++pcs, ++i) {
689 int code = gx_path_add_curve_notes(ppath, pcs->p1.x, pcs->p1.y,
690 pcs->p2.x, pcs->p2.y,
691 pcs->pt.x, pcs->pt.y,
692 notes |
693 (i > 0 ? sn_not_first :
694 sn_none));
695
696 if (code < 0)
697 return code;
698 }
699
700 return 0;
701 }
702
703 /*
704 * Split a curve if necessary into pieces that are monotonic in X or Y as a
705 * function of the curve parameter t. This allows us to rasterize curves
706 * directly without pre-flattening. This takes a fair amount of analysis....
707 * Store the values of t of the split points in pst[0] and pst[1]. Return
708 * the number of split points (0, 1, or 2).
709 */
710 int
gx_curve_monotonic_points(fixed v0,fixed v1,fixed v2,fixed v3,double pst[2])711 gx_curve_monotonic_points(fixed v0, fixed v1, fixed v2, fixed v3,
712 double pst[2])
713 {
714 /*
715 Let
716 v(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
717 Then
718 dv(t) = 3*a*t^2 + 2*b*t + c.
719 v(t) has a local minimum or maximum (or inflection point)
720 precisely where dv(t) = 0. Now the roots of dv(t) = 0 (i.e.,
721 the zeros of dv(t)) are at
722 t = ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
723 = ( -b +/- sqrt(b^2 - 3*a*c) ) / 3*a
724 (Note that real roots exist iff b^2 >= 3*a*c.)
725 We want to know if these lie in the range (0..1).
726 (The endpoints don't count.) Call such a root a "valid zero."
727 Since computing the roots is expensive, we would like to have
728 some cheap tests to filter out cases where they don't exist
729 (i.e., where the curve is already monotonic).
730 */
731 fixed v01, v12, a, b, c, b2, a3;
732 fixed dv_end, b2abs, a3abs;
733
734 curve_points_to_coefficients(v0, v1, v2, v3, a, b, c, v01, v12);
735 b2 = b << 1;
736 a3 = (a << 1) + a;
737 /*
738 If a = 0, the only possible zero is t = -c / 2*b.
739 This zero is valid iff sign(c) != sign(b) and 0 < |c| < 2*|b|.
740 */
741 if (a == 0) {
742 if ((b ^ c) < 0 && any_abs(c) < any_abs(b2) && c != 0) {
743 *pst = (double)(-c) / b2;
744 return 1;
745 } else
746 return 0;
747 }
748 /*
749 Iff a curve is horizontal at t = 0, c = 0. In this case,
750 there can be at most one other zero, at -2*b / 3*a.
751 This zero is valid iff sign(a) != sign(b) and 0 < 2*|b| < 3*|a|.
752 */
753 if (c == 0) {
754 if ((a ^ b) < 0 && any_abs(b2) < any_abs(a3) && b != 0) {
755 *pst = (double)(-b2) / a3;
756 return 1;
757 } else
758 return 0;
759 }
760 /*
761 Similarly, iff a curve is horizontal at t = 1, 3*a + 2*b + c = 0.
762 In this case, there can be at most one other zero,
763 at -1 - 2*b / 3*a, iff sign(a) != sign(b) and 1 < -2*b / 3*a < 2,
764 i.e., 3*|a| < 2*|b| < 6*|a|.
765 */
766 else if ((dv_end = a3 + b2 + c) == 0) {
767 if ((a ^ b) < 0 &&
768 (b2abs = any_abs(b2)) > (a3abs = any_abs(a3)) &&
769 b2abs < a3abs << 1
770 ) {
771 *pst = (double)(-b2 - a3) / a3;
772 return 1;
773 } else
774 return 0;
775 }
776 /*
777 If sign(dv_end) != sign(c), at least one valid zero exists,
778 since dv(0) and dv(1) have opposite signs and hence
779 dv(t) must be zero somewhere in the interval [0..1].
780 */
781 else if ((dv_end ^ c) < 0);
782 /*
783 If sign(a) = sign(b), no valid zero exists,
784 since dv is monotonic on [0..1] and has the same sign
785 at both endpoints.
786 */
787 else if ((a ^ b) >= 0)
788 return 0;
789 /*
790 Otherwise, dv(t) may be non-monotonic on [0..1]; it has valid zeros
791 iff its sign anywhere in this interval is different from its sign
792 at the endpoints, which occurs iff it has an extremum in this
793 interval and the extremum is of the opposite sign from c.
794 To find this out, we look for the local extremum of dv(t)
795 by observing
796 d2v(t) = 6*a*t + 2*b
797 which has a zero only at
798 t1 = -b / 3*a
799 Now if t1 <= 0 or t1 >= 1, no valid zero exists.
800 Note that we just determined that sign(a) != sign(b), so we know t1 > 0.
801 */
802 else if (any_abs(b) >= any_abs(a3))
803 return 0;
804 /*
805 Otherwise, we just go ahead with the computation of the roots,
806 and test them for being in the correct range. Note that a valid
807 zero is an inflection point of v(t) iff d2v(t) = 0; we don't
808 bother to check for this case, since it's rare.
809 */
810 {
811 double nbf = (double)(-b);
812 double a3f = (double)a3;
813 double radicand = nbf * nbf - a3f * c;
814
815 if (radicand < 0) {
816 if_debug1('2', "[2]negative radicand = %g\n", radicand);
817 return 0;
818 } {
819 double root = sqrt(radicand);
820 int nzeros = 0;
821 double z = (nbf - root) / a3f;
822
823 /*
824 * We need to return the zeros in the correct order.
825 * We know that root is non-negative, but a3f may be either
826 * positive or negative, so we need to check the ordering
827 * explicitly.
828 */
829 if_debug2('2', "[2]zeros at %g, %g\n", z, (nbf + root) / a3f);
830 if (z > 0 && z < 1)
831 *pst = z, nzeros = 1;
832 if (root != 0) {
833 z = (nbf + root) / a3f;
834 if (z > 0 && z < 1) {
835 if (nzeros && a3f < 0) /* order is reversed */
836 pst[1] = *pst, *pst = z;
837 else
838 pst[nzeros] = z;
839 nzeros++;
840 }
841 }
842 return nzeros;
843 }
844 }
845 }
846
847 /*
848 * Split a curve at an arbitrary point t. The above midpoint split is a
849 * special case of this with t = 0.5.
850 */
851 void
gx_curve_split(fixed x0,fixed y0,const curve_segment * pc,double t,curve_segment * pc1,curve_segment * pc2)852 gx_curve_split(fixed x0, fixed y0, const curve_segment * pc, double t,
853 curve_segment * pc1, curve_segment * pc2)
854 { /*
855 * If the original function was v(t), we want to compute the points
856 * for the functions v1(T) = v(t * T) and v2(T) = v(t + (1 - t) * T).
857 * Straightforwardly,
858 * v1(T) = a*t^3*T^3 + b*t^2*T^2 + c*t*T + d
859 * i.e.
860 * a1 = a*t^3, b1 = b*t^2, c1 = c*t, d1 = d.
861 * Similarly,
862 * v2(T) = a*[t + (1-t)*T]^3 + b*[t + (1-t)*T]^2 +
863 * c*[t + (1-t)*T] + d
864 * = a*[(1-t)^3*T^3 + 3*t*(1-t)^2*T^2 + 3*t^2*(1-t)*T +
865 * t^3] + b*[(1-t)^2*T^2 + 2*t*(1-t)*T + t^2] +
866 * c*[(1-t)*T + t] + d
867 * = a*(1-t)^3*T^3 + [a*3*t + b]*(1-t)^2*T^2 +
868 * [a*3*t^2 + b*2*t + c]*(1-t)*T +
869 * a*t^3 + b*t^2 + c*t + d
870 * We do this in the simplest way, namely, we convert the points to
871 * coefficients, do the arithmetic, and convert back. It would
872 * obviously be faster to do the arithmetic directly on the points,
873 * as the midpoint code does; this is just an implementation issue
874 * that we can revisit if necessary.
875 */
876 double t2 = t * t, t3 = t2 * t;
877 double omt = 1 - t, omt2 = omt * omt, omt3 = omt2 * omt;
878 fixed v01, v12, a, b, c, na, nb, nc;
879
880 if_debug1('2', "[2]splitting at t = %g\n", t);
881 #define compute_seg(v0, v)\
882 curve_points_to_coefficients(v0, pc->p1.v, pc->p2.v, pc->pt.v,\
883 a, b, c, v01, v12);\
884 na = (fixed)(a * t3), nb = (fixed)(b * t2), nc = (fixed)(c * t);\
885 curve_coefficients_to_points(na, nb, nc, v0,\
886 pc1->p1.v, pc1->p2.v, pc1->pt.v);\
887 na = (fixed)(a * omt3);\
888 nb = (fixed)((a * t * 3 + b) * omt2);\
889 nc = (fixed)((a * t2 * 3 + b * 2 * t + c) * omt);\
890 curve_coefficients_to_points(na, nb, nc, pc1->pt.v,\
891 pc2->p1.v, pc2->p2.v, pc2->pt.v)
892 compute_seg(x0, x);
893 compute_seg(y0, y);
894 #undef compute_seg
895 }
896