xref: /plan9/sys/src/cmd/gs/src/gxpcopy.c (revision 593dc095aefb2a85c828727bbfa9da139a49bdf4)
1 /* Copyright (C) 1992, 2000 Aladdin Enterprises.  All rights reserved.
2 
3   This software is provided AS-IS with no warranty, either express or
4   implied.
5 
6   This software is distributed under license and may not be copied,
7   modified or distributed except as expressly authorized under the terms
8   of the license contained in the file LICENSE in this distribution.
9 
10   For more information about licensing, please refer to
11   http://www.ghostscript.com/licensing/. For information on
12   commercial licensing, go to http://www.artifex.com/licensing/ or
13   contact Artifex Software, Inc., 101 Lucas Valley Road #110,
14   San Rafael, CA  94903, U.S.A., +1(415)492-9861.
15 */
16 
17 /* $Id: gxpcopy.c,v 1.26 2005/08/30 06:38:44 igor Exp $ */
18 /* Path copying and flattening */
19 #include "math_.h"
20 #include "gx.h"
21 #include "gserrors.h"
22 #include "gconfigv.h"		/* for USE_FPU */
23 #include "gxfixed.h"
24 #include "gxfarith.h"
25 #include "gxistate.h"		/* for access to line params */
26 #include "gzpath.h"
27 #include "vdtrace.h"
28 
29 /* Forward declarations */
30 private void adjust_point_to_tangent(segment *, const segment *,
31 				     const gs_fixed_point *);
32 /* Copy a path, optionally flattening or monotonizing it. */
33 /* If the copy fails, free the new path. */
34 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)35 gx_path_copy_reducing(const gx_path *ppath_old, gx_path *ppath,
36 		      fixed fixed_flatness, const gs_imager_state *pis,
37 		      gx_path_copy_options options)
38 {
39     const segment *pseg;
40     fixed flat = fixed_flatness;
41     gs_fixed_point expansion;
42     /*
43      * Since we're going to be adding to the path, unshare it
44      * before we start.
45      */
46     int code = gx_path_unshare(ppath);
47 
48     if (code < 0)
49 	return code;
50 #ifdef DEBUG
51     if (gs_debug_c('P'))
52 	gx_dump_path(ppath_old, "before reducing");
53 #endif
54     if (options & pco_for_stroke) {
55 	/* Precompute the maximum expansion of the bounding box. */
56 	double width = pis->line_params.half_width;
57 
58 	expansion.x =
59 	    float2fixed((fabs(pis->ctm.xx) + fabs(pis->ctm.yx)) * width) * 2;
60 	expansion.y =
61 	    float2fixed((fabs(pis->ctm.xy) + fabs(pis->ctm.yy)) * width) * 2;
62     }
63     vd_setcolor(RGB(255,255,0));
64     pseg = (const segment *)(ppath_old->first_subpath);
65     while (pseg) {
66 	switch (pseg->type) {
67 	    case s_start:
68 		code = gx_path_add_point(ppath,
69 					 pseg->pt.x, pseg->pt.y);
70 		vd_moveto(pseg->pt.x, pseg->pt.y);
71 		break;
72 	    case s_curve:
73 		{
74 		    const curve_segment *pc = (const curve_segment *)pseg;
75 
76 		    if (fixed_flatness == max_fixed) {	/* don't flatten */
77 			if (options & pco_monotonize)
78 			    code = gx_curve_monotonize(ppath, pc);
79 			else
80 			    code = gx_path_add_curve_notes(ppath,
81 				     pc->p1.x, pc->p1.y, pc->p2.x, pc->p2.y,
82 					   pc->pt.x, pc->pt.y, pseg->notes);
83 		    } else {
84 			fixed x0 = ppath->position.x;
85 			fixed y0 = ppath->position.y;
86 			segment_notes notes = pseg->notes;
87 			curve_segment cseg;
88 			int k;
89 
90 			if (options & pco_for_stroke) {
91 			    /*
92 			     * When flattening for stroking, the flatness
93 			     * must apply to the outside of the resulting
94 			     * stroked region.  We approximate this by
95 			     * dividing the flatness by the ratio of the
96 			     * expanded bounding box to the original
97 			     * bounding box.  This is crude, but pretty
98 			     * simple to calculate, and produces reasonably
99 			     * good results.
100 			     */
101 			    fixed min01, max01, min23, max23;
102 			    fixed ex, ey, flat_x, flat_y;
103 
104 #define SET_EXTENT(r, c0, c1, c2, c3)\
105     BEGIN\
106 	if (c0 < c1) min01 = c0, max01 = c1;\
107 	else         min01 = c1, max01 = c0;\
108 	if (c2 < c3) min23 = c2, max23 = c3;\
109 	else         min23 = c3, max23 = c2;\
110 	r = max(max01, max23) - min(min01, min23);\
111     END
112 			    SET_EXTENT(ex, x0, pc->p1.x, pc->p2.x, pc->pt.x);
113 			    SET_EXTENT(ey, y0, pc->p1.y, pc->p2.y, pc->pt.y);
114 #undef SET_EXTENT
115 			    /*
116 			     * We check for the degenerate case specially
117 			     * to avoid a division by zero.
118 			     */
119 			    if (ex == 0 || ey == 0)
120 				k = 0;
121 			    else {
122 				flat_x =
123 				    fixed_mult_quo(fixed_flatness, ex,
124 						   ex + expansion.x);
125 				flat_y =
126 				    fixed_mult_quo(fixed_flatness, ey,
127 						   ey + expansion.y);
128 				flat = min(flat_x, flat_y);
129 				k = gx_curve_log2_samples(x0, y0, pc, flat);
130 			    }
131 			} else
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 			    vd_lineto(x0, y0);
146 			    start = ppath->current_subpath->last;
147 			    notes |= sn_not_first;
148 			    cseg = *pc;
149 			    code = gx_subdivide_curve(ppath, k, &cseg, notes);
150 			    if (code < 0)
151 				break;
152 			    /*
153 			     * Adjust the first and last segments so that
154 			     * they line up with the tangents.
155 			     */
156 			    end = ppath->current_subpath->last;
157 			    vd_lineto(ppath->position.x, ppath->position.y);
158 			    if ((code = gx_path_add_line_notes(ppath,
159 							  ppath->position.x,
160 							  ppath->position.y,
161 					    pseg->notes | sn_not_first)) < 0)
162 				break;
163 			    if (start->next->pt.x != pc->p1.x || start->next->pt.y != pc->p1.y)
164 				adjust_point_to_tangent(start, start->next, &pc->p1);
165 			    else if (start->next->pt.x != pc->p2.x || start->next->pt.y != pc->p2.y)
166 				adjust_point_to_tangent(start, start->next, &pc->p2);
167 			    else
168 				adjust_point_to_tangent(start, start->next, &end->prev->pt);
169 			    if (end->prev->pt.x != pc->p2.x || end->prev->pt.y != pc->p2.y)
170 				adjust_point_to_tangent(end, end->prev, &pc->p2);
171 			    else if (end->prev->pt.x != pc->p1.x || end->prev->pt.y != pc->p1.y)
172 				adjust_point_to_tangent(end, end->prev, &pc->p1);
173 			    else
174 				adjust_point_to_tangent(end, end->prev, &start->pt);
175 			} else {
176 			    cseg = *pc;
177 			    code = gx_subdivide_curve(ppath, k, &cseg, notes);
178 			}
179 		    }
180 		    break;
181 		}
182 	    case s_line:
183 		code = gx_path_add_line_notes(ppath,
184 				       pseg->pt.x, pseg->pt.y, pseg->notes);
185 		vd_lineto(pseg->pt.x, pseg->pt.y);
186 		break;
187 	    case s_line_close:
188 		code = gx_path_close_subpath(ppath);
189 		vd_closepath;
190 		break;
191 	    default:		/* can't happen */
192 		code = gs_note_error(gs_error_unregistered);
193 	}
194 	if (code < 0) {
195 	    gx_path_new(ppath);
196 	    return code;
197 	}
198 	pseg = pseg->next;
199     }
200     if (path_last_is_moveto(ppath_old))
201 	gx_path_add_point(ppath, ppath_old->position.x,
202 			  ppath_old->position.y);
203     if (ppath_old->bbox_set) {
204 	if (ppath->bbox_set) {
205 	    ppath->bbox.p.x = min(ppath_old->bbox.p.x, ppath->bbox.p.x);
206 	    ppath->bbox.p.y = min(ppath_old->bbox.p.y, ppath->bbox.p.y);
207 	    ppath->bbox.q.x = max(ppath_old->bbox.q.x, ppath->bbox.q.x);
208 	    ppath->bbox.q.y = max(ppath_old->bbox.q.y, ppath->bbox.q.y);
209 	} else {
210 	    ppath->bbox_set = true;
211 	    ppath->bbox = ppath_old->bbox;
212 	}
213     }
214 #ifdef DEBUG
215     if (gs_debug_c('P'))
216 	gx_dump_path(ppath, "after reducing");
217 #endif
218     return 0;
219 }
220 
221 /*
222  * Adjust one end of a line (the first or last line of a flattened curve)
223  * so it falls on the curve tangent.  The closest point on the line from
224  * (0,0) to (C,D) to a point (U,V) -- i.e., the point on the line at which
225  * a perpendicular line from the point intersects it -- is given by
226  *      T = (C*U + D*V) / (C^2 + D^2)
227  *      (X,Y) = (C*T,D*T)
228  * However, any smaller value of T will also work: the one we actually
229  * use is 0.25 * the value we just derived.  We must check that
230  * numerical instabilities don't lead to a negative value of T.
231  */
232 private void
adjust_point_to_tangent(segment * pseg,const segment * next,const gs_fixed_point * p1)233 adjust_point_to_tangent(segment * pseg, const segment * next,
234 			const gs_fixed_point * p1)
235 {
236     const fixed x0 = pseg->pt.x, y0 = pseg->pt.y;
237     const fixed fC = p1->x - x0, fD = p1->y - y0;
238 
239     /*
240      * By far the commonest case is that the end of the curve is
241      * horizontal or vertical.  Check for this specially, because
242      * we can handle it with far less work (and no floating point).
243      */
244     if (fC == 0) {
245 	/* Vertical tangent. */
246 	const fixed DT = arith_rshift(next->pt.y - y0, 2);
247 
248 	if (fD == 0)
249 	    return;		/* anomalous case */
250 	if_debug1('2', "[2]adjusting vertical: DT = %g\n",
251 		  fixed2float(DT));
252 	if ((DT ^ fD) > 0)
253 	    pseg->pt.y = DT + y0;
254     } else if (fD == 0) {
255 	/* Horizontal tangent. */
256 	const fixed CT = arith_rshift(next->pt.x - x0, 2);
257 
258 	if_debug1('2', "[2]adjusting horizontal: CT = %g\n",
259 		  fixed2float(CT));
260 	if ((CT ^ fC) > 0)
261 	    pseg->pt.x = CT + x0;
262     } else {
263 	/* General case. */
264 	const double C = fC, D = fD;
265 	double T = (C * (next->pt.x - x0) + D * (next->pt.y - y0)) /
266 	(C * C + D * D);
267 
268 	if_debug3('2', "[2]adjusting: C = %g, D = %g, T = %g\n",
269 		  C, D, T);
270 	if (T > 0) {
271 	    if (T > 1) {
272 		/* Don't go outside the curve bounding box. */
273 		T = 1;
274 	    }
275 	    pseg->pt.x = arith_rshift((fixed) (C * T), 2) + x0;
276 	    pseg->pt.y = arith_rshift((fixed) (D * T), 2) + y0;
277 	}
278     }
279 }
280 
281 /* ---------------- Monotonic curves ---------------- */
282 
283 /* Test whether a path is free of non-monotonic curves. */
284 bool
gx_path__check_curves(const gx_path * ppath,gx_path_copy_options options,fixed fixed_flat)285 gx_path__check_curves(const gx_path * ppath, gx_path_copy_options options, fixed fixed_flat)
286 {
287     const segment *pseg = (const segment *)(ppath->first_subpath);
288     gs_fixed_point pt0;
289 
290     while (pseg) {
291 	switch (pseg->type) {
292 	    case s_start:
293 		{
294 		    const subpath *psub = (const subpath *)pseg;
295 
296 		    /* Skip subpaths without curves. */
297 		    if (!psub->curve_count)
298 			pseg = psub->last;
299 		}
300 		break;
301 	    case s_curve:
302 		{
303 		    const curve_segment *pc = (const curve_segment *)pseg;
304 
305 		    if (options & pco_monotonize) {
306 			double t[2];
307 			int nz = gx_curve_monotonic_points(pt0.y,
308 					       pc->p1.y, pc->p2.y, pc->pt.y, t);
309 
310 			if (nz != 0)
311 			    return false;
312 			nz = gx_curve_monotonic_points(pt0.x,
313 					       pc->p1.x, pc->p2.x, pc->pt.x, t);
314 			if (nz != 0)
315 			    return false;
316 		    }
317 		    if (options & pco_small_curves) {
318 			fixed ax, bx, cx, ay, by, cy;
319 			int k = gx_curve_log2_samples(pt0.x, pt0.y, pc, fixed_flat);
320 
321 			if(!curve_coeffs_ranged(pt0.x, pc->p1.x, pc->p2.x, pc->pt.x,
322 				pt0.y, pc->p1.y, pc->p2.y, pc->pt.y,
323 				&ax, &bx, &cx, &ay, &by, &cy, k))
324 			    return false;
325 		    }
326 		}
327 		break;
328 	    default:
329 		;
330 	}
331 	pt0 = pseg->pt;
332 	pseg = pseg->next;
333     }
334     return true;
335 }
336 
337 /* Monotonize a curve, by splitting it if necessary. */
338 /* In the worst case, this could split the curve into 9 pieces. */
339 int
gx_curve_monotonize(gx_path * ppath,const curve_segment * pc)340 gx_curve_monotonize(gx_path * ppath, const curve_segment * pc)
341 {
342     fixed x0 = ppath->position.x, y0 = ppath->position.y;
343     segment_notes notes = pc->notes;
344     double t[4], tt = 1, tp;
345     int c[4];
346     int n0, n1, n, i, j, k = 0;
347     fixed ax, bx, cx, ay, by, cy, v01, v12;
348     fixed px, py, qx, qy, rx, ry, sx, sy;
349     const double delta = 0.0000001;
350 
351     /* Roots of the derivative : */
352     n0 = gx_curve_monotonic_points(x0, pc->p1.x, pc->p2.x, pc->pt.x, t);
353     n1 = gx_curve_monotonic_points(y0, pc->p1.y, pc->p2.y, pc->pt.y, t + n0);
354     n = n0 + n1;
355     if (n == 0)
356 	return gx_path_add_curve_notes(ppath, pc->p1.x, pc->p1.y,
357 		pc->p2.x, pc->p2.y, pc->pt.x, pc->pt.y, notes);
358     if (n0 > 0)
359 	c[0] = 1;
360     if (n0 > 1)
361 	c[1] = 1;
362     if (n1 > 0)
363 	c[n0] = 2;
364     if (n1 > 1)
365 	c[n0 + 1] = 2;
366     /* Order roots : */
367     for (i = 0; i < n; i++)
368 	for (j = i + 1; j < n; j++)
369 	    if (t[i] > t[j]) {
370 		int w;
371 		double v = t[i]; t[i] = t[j]; t[j] = v;
372 		w = c[i]; c[i] = c[j]; c[j] = w;
373 	    }
374     /* Drop roots near zero : */
375     for (k = 0; k < n; k++)
376 	if (t[k] >= delta)
377 	    break;
378     /* Merge close roots, and drop roots at 1 : */
379     if (t[n - 1] > 1 - delta)
380 	n--;
381     for (i = k + 1, j = k; i < n && t[k] < 1 - delta; i++)
382 	if (any_abs(t[i] - t[j]) < delta) {
383 	    t[j] = (t[j] + t[i]) / 2; /* Unlikely 3 roots are close. */
384 	    c[j] |= c[i];
385 	} else {
386 	    j++;
387 	    t[j] = t[i];
388 	    c[j] = c[i];
389 	}
390     n = j + 1;
391     /* Do split : */
392     curve_points_to_coefficients(x0, pc->p1.x, pc->p2.x, pc->pt.x, ax, bx, cx, v01, v12);
393     curve_points_to_coefficients(y0, pc->p1.y, pc->p2.y, pc->pt.y, ay, by, cy, v01, v12);
394     ax *= 3, bx *= 2; /* Coefficients of the derivative. */
395     ay *= 3, by *= 2;
396     px = x0;
397     py = y0;
398     qx = (fixed)((pc->p1.x - px) * t[0] + 0.5);
399     qy = (fixed)((pc->p1.y - py) * t[0] + 0.5);
400     tp = 0;
401     for (i = k; i < n; i++) {
402 	double ti = t[i];
403 	double t2 = ti * ti, t3 = t2 * ti;
404 	double omt = 1 - ti, omt2 = omt * omt, omt3 = omt2 * omt;
405 	double x = x0 * omt3 + 3 * pc->p1.x * omt2 * ti + 3 * pc->p2.x * omt * t2 + pc->pt.x * t3;
406 	double y = y0 * omt3 + 3 * pc->p1.y * omt2 * ti + 3 * pc->p2.y * omt * t2 + pc->pt.y * t3;
407 	double ddx = (c[i] & 1 ? 0 : ax * t2 + bx * ti + cx); /* Suppress noize. */
408 	double ddy = (c[i] & 2 ? 0 : ay * t2 + by * ti + cy);
409 	fixed dx = (fixed)(ddx + 0.5);
410 	fixed dy = (fixed)(ddy + 0.5);
411 	int code;
412 
413 	tt = (i + 1 < n ? t[i + 1] : 1) - ti;
414 	rx = (fixed)(dx * (t[i] - tp) / 3 + 0.5);
415 	ry = (fixed)(dy * (t[i] - tp) / 3 + 0.5);
416 	sx = (fixed)(x + 0.5);
417 	sy = (fixed)(y + 0.5);
418 	/* Suppress the derivative sign noize near a beak : */
419 	if ((double)(sx - px) * qx + (double)(sy - py) * qy < 0)
420 	    qx = -qx, qy = -qy;
421 	if ((double)(sx - px) * rx + (double)(sy - py) * ry < 0)
422 	    rx = -rx, ry = -qy;
423 	/* Do add : */
424 	code = gx_path_add_curve_notes(ppath, px + qx, py + qy, sx - rx, sy - ry, sx, sy, notes);
425 	if (code < 0)
426 	    return code;
427 	notes |= sn_not_first;
428 	px = sx;
429 	py = sy;
430 	qx = (fixed)(dx * tt / 3 + 0.5);
431 	qy = (fixed)(dy * tt / 3 + 0.5);
432 	tp = t[i];
433     }
434     sx = pc->pt.x;
435     sy = pc->pt.y;
436     rx = (fixed)((pc->pt.x - pc->p2.x) * tt + 0.5);
437     ry = (fixed)((pc->pt.y - pc->p2.y) * tt + 0.5);
438     /* Suppress the derivative sign noize near peaks : */
439     if ((double)(sx - px) * qx + (double)(sy - py) * qy < 0)
440 	qx = -qx, qy = -qy;
441     if ((double)(sx - px) * rx + (double)(sy - py) * ry < 0)
442 	rx = -rx, ry = -qy;
443     return gx_path_add_curve_notes(ppath, px + qx, py + qy, sx - rx, sy - ry, sx, sy, notes);
444 }
445 
446 /*
447  * Split a curve if necessary into pieces that are monotonic in X or Y as a
448  * function of the curve parameter t.  This allows us to rasterize curves
449  * directly without pre-flattening.  This takes a fair amount of analysis....
450  * Store the values of t of the split points in pst[0] and pst[1].  Return
451  * the number of split points (0, 1, or 2).
452  */
453 int
gx_curve_monotonic_points(fixed v0,fixed v1,fixed v2,fixed v3,double pst[2])454 gx_curve_monotonic_points(fixed v0, fixed v1, fixed v2, fixed v3,
455 			  double pst[2])
456 {
457     /*
458        Let
459        v(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
460        Then
461        dv(t) = 3*a*t^2 + 2*b*t + c.
462        v(t) has a local minimum or maximum (or inflection point)
463        precisely where dv(t) = 0.  Now the roots of dv(t) = 0 (i.e.,
464        the zeros of dv(t)) are at
465        t =  ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
466        =    ( -b +/- sqrt(b^2 - 3*a*c) ) / 3*a
467        (Note that real roots exist iff b^2 >= 3*a*c.)
468        We want to know if these lie in the range (0..1).
469        (The endpoints don't count.)  Call such a root a "valid zero."
470        Since computing the roots is expensive, we would like to have
471        some cheap tests to filter out cases where they don't exist
472        (i.e., where the curve is already monotonic).
473      */
474     fixed v01, v12, a, b, c, b2, a3;
475     fixed dv_end, b2abs, a3abs;
476 
477     curve_points_to_coefficients(v0, v1, v2, v3, a, b, c, v01, v12);
478     b2 = b << 1;
479     a3 = (a << 1) + a;
480     /*
481        If a = 0, the only possible zero is t = -c / 2*b.
482        This zero is valid iff sign(c) != sign(b) and 0 < |c| < 2*|b|.
483      */
484     if (a == 0) {
485 	if ((b ^ c) < 0 && any_abs(c) < any_abs(b2) && c != 0) {
486 	    *pst = (double)(-c) / b2;
487 	    return 1;
488 	} else
489 	    return 0;
490     }
491     /*
492        Iff a curve is horizontal at t = 0, c = 0.  In this case,
493        there can be at most one other zero, at -2*b / 3*a.
494        This zero is valid iff sign(a) != sign(b) and 0 < 2*|b| < 3*|a|.
495      */
496     if (c == 0) {
497 	if ((a ^ b) < 0 && any_abs(b2) < any_abs(a3) && b != 0) {
498 	    *pst = (double)(-b2) / a3;
499 	    return 1;
500 	} else
501 	    return 0;
502     }
503     /*
504        Similarly, iff a curve is horizontal at t = 1, 3*a + 2*b + c = 0.
505        In this case, there can be at most one other zero,
506        at -1 - 2*b / 3*a, iff sign(a) != sign(b) and 1 < -2*b / 3*a < 2,
507        i.e., 3*|a| < 2*|b| < 6*|a|.
508      */
509     else if ((dv_end = a3 + b2 + c) == 0) {
510 	if ((a ^ b) < 0 &&
511 	    (b2abs = any_abs(b2)) > (a3abs = any_abs(a3)) &&
512 	    b2abs < a3abs << 1
513 	    ) {
514 	    *pst = (double)(-b2 - a3) / a3;
515 	    return 1;
516 	} else
517 	    return 0;
518     }
519     /*
520        If sign(dv_end) != sign(c), at least one valid zero exists,
521        since dv(0) and dv(1) have opposite signs and hence
522        dv(t) must be zero somewhere in the interval [0..1].
523      */
524     else if ((dv_end ^ c) < 0);
525     /*
526        If sign(a) = sign(b), no valid zero exists,
527        since dv is monotonic on [0..1] and has the same sign
528        at both endpoints.
529      */
530     else if ((a ^ b) >= 0)
531 	return 0;
532     /*
533        Otherwise, dv(t) may be non-monotonic on [0..1]; it has valid zeros
534        iff its sign anywhere in this interval is different from its sign
535        at the endpoints, which occurs iff it has an extremum in this
536        interval and the extremum is of the opposite sign from c.
537        To find this out, we look for the local extremum of dv(t)
538        by observing
539        d2v(t) = 6*a*t + 2*b
540        which has a zero only at
541        t1 = -b / 3*a
542        Now if t1 <= 0 or t1 >= 1, no valid zero exists.
543        Note that we just determined that sign(a) != sign(b), so we know t1 > 0.
544      */
545     else if (any_abs(b) >= any_abs(a3))
546 	return 0;
547     /*
548        Otherwise, we just go ahead with the computation of the roots,
549        and test them for being in the correct range.  Note that a valid
550        zero is an inflection point of v(t) iff d2v(t) = 0; we don't
551        bother to check for this case, since it's rare.
552      */
553     {
554 	double nbf = (double)(-b);
555 	double a3f = (double)a3;
556 	double radicand = nbf * nbf - a3f * c;
557 
558 	if (radicand < 0) {
559 	    if_debug1('2', "[2]negative radicand = %g\n", radicand);
560 	    return 0;
561 	} {
562 	    double root = sqrt(radicand);
563 	    int nzeros = 0;
564 	    double z = (nbf - root) / a3f;
565 
566 	    /*
567 	     * We need to return the zeros in the correct order.
568 	     * We know that root is non-negative, but a3f may be either
569 	     * positive or negative, so we need to check the ordering
570 	     * explicitly.
571 	     */
572 	    if_debug2('2', "[2]zeros at %g, %g\n", z, (nbf + root) / a3f);
573 	    if (z > 0 && z < 1)
574 		*pst = z, nzeros = 1;
575 	    if (root != 0) {
576 		z = (nbf + root) / a3f;
577 		if (z > 0 && z < 1) {
578 		    if (nzeros && a3f < 0)	/* order is reversed */
579 			pst[1] = *pst, *pst = z;
580 		    else
581 			pst[nzeros] = z;
582 		    nzeros++;
583 		}
584 	    }
585 	    return nzeros;
586 	}
587     }
588 }
589 
590 /* ---------------- Path optimization for the filling algorithm. ---------------- */
591 
592 private bool
find_contacting_segments(const subpath * sp0,segment * sp0last,const subpath * sp1,segment * sp1last,segment ** sc0,segment ** sc1)593 find_contacting_segments(const subpath *sp0, segment *sp0last,
594 			 const subpath *sp1, segment *sp1last,
595 			 segment **sc0, segment **sc1)
596 {
597     segment *s0, *s1;
598     const segment *s0s, *s1s;
599     int count0, count1, search_limit = 50;
600     int min_length = fixed_1 * 1;
601 
602     /* This is a simplified algorithm, which only checks for quazi-colinear vertical lines.
603        "Quazi-vertical" means dx <= 1 && dy >= min_length . */
604     /* To avoid a big unuseful expence of the processor time,
605        we search the first subpath from the end
606        (assuming that it was recently merged near the end),
607        and restrict the search with search_limit segments
608        against a quadratic scanning of two long subpaths.
609        Thus algorithm is not necessary finds anything contacting.
610        Instead it either quickly finds something, or maybe not. */
611     for (s0 = sp0last, count0 = 0; count0 < search_limit && s0 != (segment *)sp0; s0 = s0->prev, count0++) {
612 	s0s = s0->prev;
613 	if (s0->type == s_line && (s0s->pt.x == s0->pt.x ||
614 	    (any_abs(s0s->pt.x - s0->pt.x) == 1 && any_abs(s0s->pt.y - s0->pt.y) > min_length))) {
615 	    for (s1 = sp1last, count1 = 0; count1 < search_limit && s1 != (segment *)sp1; s1 = s1->prev, count1++) {
616 		s1s = s1->prev;
617 		if (s1->type == s_line && (s1s->pt.x == s1->pt.x ||
618 		    (any_abs(s1s->pt.x - s1->pt.x) == 1 && any_abs(s1s->pt.y - s1->pt.y) > min_length))) {
619 		    if (s0s->pt.x == s1s->pt.x || s0->pt.x == s1->pt.x || s0->pt.x == s1s->pt.x || s0s->pt.x == s1->pt.x) {
620 			if (s0s->pt.y < s0->pt.y && s1s->pt.y > s1->pt.y) {
621 			    fixed y0 = max(s0s->pt.y, s1->pt.y);
622 			    fixed y1 = min(s0->pt.y, s1s->pt.y);
623 
624 			    if (y0 <= y1) {
625 				*sc0 = s0;
626 				*sc1 = s1;
627 				return true;
628 			    }
629 			}
630 			if (s0s->pt.y > s0->pt.y && s1s->pt.y < s1->pt.y) {
631 			    fixed y0 = max(s0->pt.y, s1s->pt.y);
632 			    fixed y1 = min(s0s->pt.y, s1->pt.y);
633 
634 			    if (y0 <= y1) {
635 				*sc0 = s0;
636 				*sc1 = s1;
637 				return true;
638 			    }
639 			}
640 		    }
641 		}
642 	    }
643 	}
644     }
645     return false;
646 }
647 
648 int
gx_path_merge_contacting_contours(gx_path * ppath)649 gx_path_merge_contacting_contours(gx_path *ppath)
650 {
651     /* Now this is a simplified algorithm,
652        which merge only contours by a common quazi-vertical line. */
653     int window = 5/* max spot holes */ * 6/* segments per subpath */;
654     subpath *sp0 = ppath->segments->contents.subpath_first;
655 
656     for (; sp0 != NULL; sp0 = (subpath *)sp0->last->next) {
657 	segment *sp0last = sp0->last;
658 	subpath *sp1 = (subpath *)sp0last->next, *spnext;
659 	subpath *sp1p = sp0;
660 	int count;
661 
662 	for (count = 0; sp1 != NULL && count < window; sp1 = spnext, count++) {
663 	    segment *sp1last = sp1->last;
664 	    segment *sc0, *sc1;
665 
666 	    spnext = (subpath *)sp1last->next;
667 	    if (find_contacting_segments(sp0, sp0last, sp1, sp1last, &sc0, &sc1)) {
668 		/* Detach the subpath 1 from the path: */
669 		sp1->prev->next = sp1last->next;
670 		if (sp1last->next != NULL)
671 		    sp1last->next->prev = sp1->prev;
672 		sp1->prev = 0;
673 		sp1last->next = 0;
674 		/* Change 'closepath' of the subpath 1 to a line (maybe degenerate) : */
675 		if (sp1last->type == s_line_close)
676 		    sp1last->type = s_line;
677 		/* Rotate the subpath 1 to sc1 : */
678 		{   segment *old_first = sp1->next;
679 
680 		    /* Detach s_start and make a loop : */
681 		    sp1last->next = old_first;
682 		    old_first->prev = sp1last;
683 		    /* Unlink before sc1 : */
684 		    sp1last = sc1->prev;
685 		    sc1->prev->next = 0;
686 		    sc1->prev = 0; /* Safety. */
687 		    /* sp1 is not longer in use. Free it : */
688 		    if (ppath->segments->contents.subpath_current == sp1) {
689 			ppath->segments->contents.subpath_current = sp1p;
690 		    }
691 		    gs_free_object(ppath->memory, sp1, "gx_path_merge_contacting_contours");
692 		    sp1 = 0; /* Safety. */
693 		}
694 		/* Insert the subpath 1 into the subpath 0 before sc0 :*/
695 		sc0->prev->next = sc1;
696 		sc1->prev = sc0->prev;
697 		sp1last->next = sc0;
698 		sc0->prev = sp1last;
699 		/* Remove degenearte "bridge" segments : (fixme: Not done due to low importance). */
700 		/* Edit the subpath count : */
701 		ppath->subpath_count--;
702 	    } else
703 		sp1p = sp1;
704 	}
705     }
706     return 0;
707 }
708 
709