1 /* Copyright (C) 2001-2019 Artifex Software, Inc.
2 All Rights Reserved.
3
4 This software is provided AS-IS with no warranty, either express or
5 implied.
6
7 This software is distributed under license and may not be copied,
8 modified or distributed except as expressly authorized under the terms
9 of the license contained in the file LICENSE in this distribution.
10
11 Refer to licensing information at http://www.artifex.com or contact
12 Artifex Software, Inc., 1305 Grant Avenue - Suite 200, Novato,
13 CA 94945, U.S.A., +1(415)492-9861, for further information.
14 */
15
16
17 /* Rendering for non-mesh shadings */
18 #include "math_.h"
19 #include "memory_.h"
20 #include "gx.h"
21 #include "gserrors.h"
22 #include "gsmatrix.h" /* for gscoord.h */
23 #include "gscoord.h"
24 #include "gspath.h"
25 #include "gsptype2.h"
26 #include "gxcspace.h"
27 #include "gxdcolor.h"
28 #include "gxfarith.h"
29 #include "gxfixed.h"
30 #include "gxgstate.h"
31 #include "gxpath.h"
32 #include "gxshade.h"
33 #include "gxdevcli.h"
34 #include "gxshade4.h"
35 #include "gsicc_cache.h"
36
37 /* ---------------- Function-based shading ---------------- */
38
39 typedef struct Fb_frame_s { /* A rudiment of old code. */
40 gs_rect region;
41 gs_client_color cc[4]; /* colors at 4 corners */
42 int state;
43 } Fb_frame_t;
44
45 typedef struct Fb_fill_state_s {
46 shading_fill_state_common;
47 const gs_shading_Fb_t *psh;
48 gs_matrix_fixed ptm; /* parameter space -> device space */
49 Fb_frame_t frame;
50 } Fb_fill_state_t;
51 /****** NEED GC DESCRIPTOR ******/
52
53 static inline void
make_other_poles(patch_curve_t curve[4])54 make_other_poles(patch_curve_t curve[4])
55 {
56 int i, j;
57
58 for (i = 0; i < 4; i++) {
59 j = (i + 1) % 4;
60 curve[i].control[0].x = (curve[i].vertex.p.x * 2 + curve[j].vertex.p.x) / 3;
61 curve[i].control[0].y = (curve[i].vertex.p.y * 2 + curve[j].vertex.p.y) / 3;
62 curve[i].control[1].x = (curve[i].vertex.p.x + curve[j].vertex.p.x * 2) / 3;
63 curve[i].control[1].y = (curve[i].vertex.p.y + curve[j].vertex.p.y * 2) / 3;
64 curve[i].straight = true;
65 }
66 }
67
68 /* Transform a point with a fixed-point result. */
69 static void
gs_point_transform2fixed_clamped(const gs_matrix_fixed * pmat,double x,double y,gs_fixed_point * ppt)70 gs_point_transform2fixed_clamped(const gs_matrix_fixed * pmat,
71 double x, double y, gs_fixed_point * ppt)
72 {
73 gs_point fpt;
74
75 gs_point_transform(x, y, (const gs_matrix *)pmat, &fpt);
76 ppt->x = clamp_coord(fpt.x);
77 ppt->y = clamp_coord(fpt.y);
78 }
79
80 static int
Fb_fill_region(Fb_fill_state_t * pfs,const gs_fixed_rect * rect)81 Fb_fill_region(Fb_fill_state_t * pfs, const gs_fixed_rect *rect)
82 {
83 patch_fill_state_t pfs1;
84 patch_curve_t curve[4];
85 Fb_frame_t * fp = &pfs->frame;
86 int code;
87
88 memcpy(&pfs1, (shading_fill_state_t *)pfs, sizeof(shading_fill_state_t));
89 pfs1.Function = pfs->psh->params.Function;
90 code = init_patch_fill_state(&pfs1);
91 if (code < 0)
92 return code;
93 pfs1.maybe_self_intersecting = false;
94 pfs1.n_color_args = 2;
95 pfs1.rect = *rect;
96 gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.p.y, &curve[0].vertex.p);
97 gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.p.y, &curve[1].vertex.p);
98 gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.q.y, &curve[2].vertex.p);
99 gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.q.y, &curve[3].vertex.p);
100 make_other_poles(curve);
101 curve[0].vertex.cc[0] = fp->region.p.x; curve[0].vertex.cc[1] = fp->region.p.y;
102 curve[1].vertex.cc[0] = fp->region.q.x; curve[1].vertex.cc[1] = fp->region.p.y;
103 curve[2].vertex.cc[0] = fp->region.q.x; curve[2].vertex.cc[1] = fp->region.q.y;
104 curve[3].vertex.cc[0] = fp->region.p.x; curve[3].vertex.cc[1] = fp->region.q.y;
105 code = patch_fill(&pfs1, curve, NULL, NULL);
106 if (term_patch_fill_state(&pfs1))
107 return_error(gs_error_unregistered); /* Must not happen. */
108 return code;
109 }
110
111 int
gs_shading_Fb_fill_rectangle(const gs_shading_t * psh0,const gs_rect * rect,const gs_fixed_rect * rect_clip,gx_device * dev,gs_gstate * pgs)112 gs_shading_Fb_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,
113 const gs_fixed_rect * rect_clip,
114 gx_device * dev, gs_gstate * pgs)
115 {
116 const gs_shading_Fb_t * const psh = (const gs_shading_Fb_t *)psh0;
117 gs_matrix save_ctm;
118 int xi, yi, code;
119 float x[2], y[2];
120 Fb_fill_state_t state;
121
122 code = shade_init_fill_state((shading_fill_state_t *) & state, psh0, dev, pgs);
123 if (code < 0)
124 return code;
125 state.psh = psh;
126 /****** HACK FOR FIXED-POINT MATRIX MULTIPLY ******/
127 gs_currentmatrix((gs_gstate *) pgs, &save_ctm);
128 gs_concat((gs_gstate *) pgs, &psh->params.Matrix);
129 state.ptm = pgs->ctm;
130 gs_setmatrix((gs_gstate *) pgs, &save_ctm);
131 /* Compute the parameter X and Y ranges. */
132 {
133 gs_rect pbox;
134
135 code = gs_bbox_transform_inverse(rect, &psh->params.Matrix, &pbox);
136 if (code < 0)
137 return code;
138 x[0] = max(pbox.p.x, psh->params.Domain[0]);
139 x[1] = min(pbox.q.x, psh->params.Domain[1]);
140 y[0] = max(pbox.p.y, psh->params.Domain[2]);
141 y[1] = min(pbox.q.y, psh->params.Domain[3]);
142 }
143 if (x[0] > x[1] || y[0] > y[1]) {
144 /* The region is outside the shading area. */
145 if (state.icclink != NULL) gsicc_release_link(state.icclink);
146 return 0;
147 }
148 for (xi = 0; xi < 2; ++xi)
149 for (yi = 0; yi < 2; ++yi) {
150 float v[2];
151
152 v[0] = x[xi], v[1] = y[yi];
153 gs_function_evaluate(psh->params.Function, v,
154 state.frame.cc[yi * 2 + xi].paint.values);
155 }
156 state.frame.region.p.x = x[0];
157 state.frame.region.p.y = y[0];
158 state.frame.region.q.x = x[1];
159 state.frame.region.q.y = y[1];
160 code = Fb_fill_region(&state, rect_clip);
161 if (state.icclink != NULL) gsicc_release_link(state.icclink);
162 return code;
163 }
164
165 /* ---------------- Axial shading ---------------- */
166
167 typedef struct A_fill_state_s {
168 const gs_shading_A_t *psh;
169 gs_point delta;
170 double length;
171 double t0, t1;
172 double v0, v1, u0, u1;
173 } A_fill_state_t;
174 /****** NEED GC DESCRIPTOR ******/
175
176 /* Note t0 and t1 vary over [0..1], not the Domain. */
177
178 static int
A_fill_region(A_fill_state_t * pfs,patch_fill_state_t * pfs1)179 A_fill_region(A_fill_state_t * pfs, patch_fill_state_t *pfs1)
180 {
181 const gs_shading_A_t * const psh = pfs->psh;
182 double x0 = psh->params.Coords[0] + pfs->delta.x * pfs->v0;
183 double y0 = psh->params.Coords[1] + pfs->delta.y * pfs->v0;
184 double x1 = psh->params.Coords[0] + pfs->delta.x * pfs->v1;
185 double y1 = psh->params.Coords[1] + pfs->delta.y * pfs->v1;
186 double h0 = pfs->u0, h1 = pfs->u1;
187 patch_curve_t curve[4];
188 int code;
189
190 code = gs_point_transform2fixed(&pfs1->pgs->ctm, x0 + pfs->delta.y * h0, y0 - pfs->delta.x * h0, &curve[0].vertex.p);
191 if (code < 0)
192 return code;
193 code = gs_point_transform2fixed(&pfs1->pgs->ctm, x1 + pfs->delta.y * h0, y1 - pfs->delta.x * h0, &curve[1].vertex.p);
194 if (code < 0)
195 return code;
196 code = gs_point_transform2fixed(&pfs1->pgs->ctm, x1 + pfs->delta.y * h1, y1 - pfs->delta.x * h1, &curve[2].vertex.p);
197 if (code < 0)
198 return code;
199 code = gs_point_transform2fixed(&pfs1->pgs->ctm, x0 + pfs->delta.y * h1, y0 - pfs->delta.x * h1, &curve[3].vertex.p);
200 if (code < 0)
201 return code;
202 curve[0].vertex.cc[0] = pfs->t0; /* The element cc[1] is set to a dummy value against */
203 curve[1].vertex.cc[0] = pfs->t1; /* interrupts while an idle priocessing in gxshade.6.c . */
204 curve[2].vertex.cc[0] = pfs->t1;
205 curve[3].vertex.cc[0] = pfs->t0;
206 curve[0].vertex.cc[1] = 0; /* The element cc[1] is set to a dummy value against */
207 curve[1].vertex.cc[1] = 0; /* interrupts while an idle priocessing in gxshade.6.c . */
208 curve[2].vertex.cc[1] = 0;
209 curve[3].vertex.cc[1] = 0;
210 make_other_poles(curve);
211 return patch_fill(pfs1, curve, NULL, NULL);
212 }
213
214 static inline int
gs_shading_A_fill_rectangle_aux(const gs_shading_t * psh0,const gs_rect * rect,const gs_fixed_rect * clip_rect,gx_device * dev,gs_gstate * pgs)215 gs_shading_A_fill_rectangle_aux(const gs_shading_t * psh0, const gs_rect * rect,
216 const gs_fixed_rect *clip_rect,
217 gx_device * dev, gs_gstate * pgs)
218 {
219 const gs_shading_A_t *const psh = (const gs_shading_A_t *)psh0;
220 gs_function_t * const pfn = psh->params.Function;
221 gs_matrix cmat;
222 gs_rect t_rect;
223 A_fill_state_t state;
224 patch_fill_state_t pfs1;
225 float d0 = psh->params.Domain[0], d1 = psh->params.Domain[1];
226 float dd = d1 - d0;
227 double t0, t1;
228 gs_point dist;
229 int code;
230
231 state.psh = psh;
232 code = shade_init_fill_state((shading_fill_state_t *)&pfs1, psh0, dev, pgs);
233 if (code < 0)
234 return code;
235 pfs1.Function = pfn;
236 pfs1.rect = *clip_rect;
237 code = init_patch_fill_state(&pfs1);
238 if (code < 0)
239 goto fail;
240 pfs1.maybe_self_intersecting = false;
241 pfs1.function_arg_shift = 1;
242 /*
243 * Compute the parameter range. We construct a matrix in which
244 * (0,0) corresponds to t = 0 and (0,1) corresponds to t = 1,
245 * and use it to inverse-map the rectangle to be filled.
246 */
247 cmat.tx = psh->params.Coords[0];
248 cmat.ty = psh->params.Coords[1];
249 state.delta.x = psh->params.Coords[2] - psh->params.Coords[0];
250 state.delta.y = psh->params.Coords[3] - psh->params.Coords[1];
251 cmat.yx = state.delta.x;
252 cmat.yy = state.delta.y;
253 cmat.xx = cmat.yy;
254 cmat.xy = -cmat.yx;
255 code = gs_bbox_transform_inverse(rect, &cmat, &t_rect);
256 if (code < 0) {
257 code = 0; /* Swallow this silently */
258 goto fail;
259 }
260 t0 = min(max(t_rect.p.y, 0), 1);
261 t1 = max(min(t_rect.q.y, 1), 0);
262 state.v0 = t0;
263 state.v1 = t1;
264 state.u0 = t_rect.p.x;
265 state.u1 = t_rect.q.x;
266 state.t0 = t0 * dd + d0;
267 state.t1 = t1 * dd + d0;
268 code = gs_distance_transform(state.delta.x, state.delta.y, &ctm_only(pgs),
269 &dist);
270 if (code < 0)
271 goto fail;
272 state.length = hypot(dist.x, dist.y); /* device space line length */
273 code = A_fill_region(&state, &pfs1);
274 if (psh->params.Extend[0] && t0 > t_rect.p.y) {
275 if (code < 0)
276 goto fail;
277 /* Use the general algorithm, because we need the trapping. */
278 state.v0 = t_rect.p.y;
279 state.v1 = t0;
280 state.t0 = state.t1 = t0 * dd + d0;
281 code = A_fill_region(&state, &pfs1);
282 }
283 if (psh->params.Extend[1] && t1 < t_rect.q.y) {
284 if (code < 0)
285 goto fail;
286 /* Use the general algorithm, because we need the trapping. */
287 state.v0 = t1;
288 state.v1 = t_rect.q.y;
289 state.t0 = state.t1 = t1 * dd + d0;
290 code = A_fill_region(&state, &pfs1);
291 }
292 fail:
293 gsicc_release_link(pfs1.icclink);
294 if (term_patch_fill_state(&pfs1))
295 return_error(gs_error_unregistered); /* Must not happen. */
296 return code;
297 }
298
299 int
gs_shading_A_fill_rectangle(const gs_shading_t * psh0,const gs_rect * rect,const gs_fixed_rect * rect_clip,gx_device * dev,gs_gstate * pgs)300 gs_shading_A_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,
301 const gs_fixed_rect * rect_clip,
302 gx_device * dev, gs_gstate * pgs)
303 {
304 return gs_shading_A_fill_rectangle_aux(psh0, rect, rect_clip, dev, pgs);
305 }
306
307 /* ---------------- Radial shading ---------------- */
308
309 /* Some notes on what I have struggled to understand about the following
310 * function. This function renders the 'tube' given by interpolating one
311 * circle to another.
312 *
313 * The first circle is at (x0, y0) with radius r0, and has 'color' t0.
314 * The other circle is at (x1, y1) with radius r1, and has 'color' t1.
315 *
316 * We perform this rendering by approximating each quadrant of the 'tube'
317 * by a tensor patch. The tensor patch is formed by taking a curve along
318 * 1/4 of the circumference of the first circle, a straight line to the
319 * equivalent point on the circumference of the second circle, a curve
320 * back along the circumference of the second circle, and then a straight
321 * line back to where we started.
322 *
323 * There is additional logic in this function that forms the directions of
324 * the curves differently for different quadrants. This is done to ensure
325 * that we always paint 'around' the tube from the back towards the front,
326 * so we don't get unexpected regions showing though. This is explained more
327 * below.
328 *
329 * The original code here examined the position change between the two
330 * circles dx and dy. Based upon this vector it would pick which quadrant/
331 * tensor patch to draw first. It would draw the quadrants/tensor patches
332 * in anticlockwise order. Presumably this was intended to be done so that
333 * the 'top' quadrant would be drawn last.
334 *
335 * Unfortunately this did not always work; see bug 692513. If the quadrants
336 * were rendered in the order 0,1,2,3, the rendering of 1 was leaving traces
337 * on top of 0, which was unexpected.
338 *
339 * I have therefore altered the code slightly; rather than picking a start
340 * quadrant and moving anticlockwise, we now draw the 'undermost' quadrant,
341 * then the two adjacent quadrants, then the topmost quadrant.
342 *
343 * For the purposes of explanation, we shall label the octants as below:
344 *
345 * \2|1/ and Quadrants as: |
346 * 3\|/0 Q1 | Q0
347 * ---+--- ----+----
348 * 4/|\7 Q2 | Q3
349 * /5|6\ |
350 *
351 * We find (dx,dy), the difference between the centres of the circles.
352 * We look to see which octant this falls in. Firstly, this tells us which
353 * quadrant of the circle we need to draw first (Octant n, starts with
354 * Quadrant floor(n/2)). Secondly, it tells us which direction to form the
355 * tensor patch in; we always want to draw from the side 'closest' to
356 * dx/dy to the side further away. This ensures that we don't overwrite
357 * pixels in the incorrect order as the patch decomposes.
358 */
359 static int
R_tensor_annulus(patch_fill_state_t * pfs,double x0,double y0,double r0,double t0,double x1,double y1,double r1,double t1)360 R_tensor_annulus(patch_fill_state_t *pfs,
361 double x0, double y0, double r0, double t0,
362 double x1, double y1, double r1, double t1)
363 {
364 double dx = x1 - x0, dy = y1 - y0;
365 double d = hypot(dx, dy);
366 gs_point p0, p1, pc0, pc1;
367 int k, j, code, dirn;
368 bool inside = 0;
369
370 /* pc0 and pc1 are the centres of the respective circles. */
371 pc0.x = x0, pc0.y = y0;
372 pc1.x = x1, pc1.y = y1;
373 /* Set p0 up so it's a unit vector giving the direction of 90 degrees
374 * to the right of the major axis as we move from p0c to p1c. */
375 if (r0 + d <= r1 || r1 + d <= r0) {
376 /* One circle is inside another one.
377 Use any subdivision,
378 but don't depend on dx, dy, which may be too small. */
379 p0.x = 0, p0.y = -1, dirn = 0;
380 /* Align stripes along radii for faster triangulation : */
381 inside = 1;
382 pfs->function_arg_shift = 1;
383 } else {
384 /* Must generate canonic quadrangle arcs,
385 because we approximate them with curves. */
386 if(dx >= 0) {
387 if (dy >= 0)
388 p0.x = 1, p0.y = 0, dirn = (dx >= dy ? 1 : 0);
389 else
390 p0.x = 0, p0.y = -1, dirn = (dx >= -dy ? 0 : 1);
391 } else {
392 if (dy >= 0)
393 p0.x = 0, p0.y = 1, dirn = (-dx >= dy ? 1 : 0);
394 else
395 p0.x = -1, p0.y = 0, dirn = (-dx >= -dy ? 0 : 1);
396 }
397 pfs->function_arg_shift = 0;
398 }
399 /* fixme: wish: cut invisible parts off.
400 Note : when r0 != r1 the invisible part is not a half circle. */
401 for (k = 0; k < 4; k++) {
402 gs_point p[12];
403 patch_curve_t curve[4];
404
405 /* Set p1 to be 90 degrees anticlockwise from p0 */
406 p1.x = -p0.y; p1.y = p0.x;
407 if (dirn == 0) { /* Clockwise */
408 make_quadrant_arc(p + 0, &pc0, &p1, &p0, r0);
409 make_quadrant_arc(p + 6, &pc1, &p0, &p1, r1);
410 } else { /* Anticlockwise */
411 make_quadrant_arc(p + 0, &pc0, &p0, &p1, r0);
412 make_quadrant_arc(p + 6, &pc1, &p1, &p0, r1);
413 }
414 p[4].x = (p[3].x * 2 + p[6].x) / 3;
415 p[4].y = (p[3].y * 2 + p[6].y) / 3;
416 p[5].x = (p[3].x + p[6].x * 2) / 3;
417 p[5].y = (p[3].y + p[6].y * 2) / 3;
418 p[10].x = (p[9].x * 2 + p[0].x) / 3;
419 p[10].y = (p[9].y * 2 + p[0].y) / 3;
420 p[11].x = (p[9].x + p[0].x * 2) / 3;
421 p[11].y = (p[9].y + p[0].y * 2) / 3;
422 for (j = 0; j < 4; j++) {
423 int jj = (j + inside) % 4;
424
425 if (gs_point_transform2fixed(&pfs->pgs->ctm, p[j*3 + 0].x, p[j*3 + 0].y, &curve[jj].vertex.p) < 0)
426 gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 0].x, p[j*3 + 0].y, &curve[jj].vertex.p);
427
428 if (gs_point_transform2fixed(&pfs->pgs->ctm, p[j*3 + 1].x, p[j*3 + 1].y, &curve[jj].control[0]) < 0)
429 gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 1].x, p[j*3 + 1].y, &curve[jj].control[0]);
430
431 if (gs_point_transform2fixed(&pfs->pgs->ctm, p[j*3 + 2].x, p[j*3 + 2].y, &curve[jj].control[1]) < 0)
432 gs_point_transform2fixed_clamped(&pfs->pgs->ctm, p[j*3 + 2].x, p[j*3 + 2].y, &curve[jj].control[1]);
433 curve[j].straight = (((j + inside) & 1) != 0);
434 }
435 curve[(0 + inside) % 4].vertex.cc[0] = t0;
436 curve[(1 + inside) % 4].vertex.cc[0] = t0;
437 curve[(2 + inside) % 4].vertex.cc[0] = t1;
438 curve[(3 + inside) % 4].vertex.cc[0] = t1;
439 curve[0].vertex.cc[1] = curve[1].vertex.cc[1] = 0; /* Initialize against FPE. */
440 curve[2].vertex.cc[1] = curve[3].vertex.cc[1] = 0; /* Initialize against FPE. */
441 code = patch_fill(pfs, curve, NULL, NULL);
442 if (code < 0)
443 return code;
444 /* Move p0 to be ready for the next position */
445 if (k == 0) {
446 /* p0 moves clockwise */
447 p1 = p0;
448 p0.x = p1.y; p0.y = -p1.x;
449 dirn = 0;
450 } else if (k == 1) {
451 /* p0 flips sides */
452 p0.x = -p0.x; p0.y = -p0.y;
453 dirn = 1;
454 } else if (k == 2) {
455 /* p0 moves anti-clockwise */
456 p1 = p0;
457 p0.x = -p1.y; p0.y = p1.x;
458 dirn = 0;
459 }
460 }
461 return 0;
462 }
463
464 /* Find the control points for two points on the arc of a circle
465 * the points must be within the same quadrant.
466 */
find_arc_control_points(gs_point * from,gs_point * to,gs_point * from_control,gs_point * to_control,gs_point * centre)467 static int find_arc_control_points(gs_point *from, gs_point *to, gs_point *from_control, gs_point *to_control, gs_point *centre)
468 {
469 double from_tan_alpha, to_tan_alpha, from_alpha, to_alpha;
470 double half_inscribed_angle, intersect_x, intersect_y, intersect_dist;
471 double radius = sqrt(((from->x - centre->x) * (from->x - centre->x)) + ((from->y - centre->y) * (from->y - centre->y)));
472 double tangent_intersect_dist;
473 double F;
474 int quadrant;
475
476 /* Quadrant 0 is upper right, numbered anti-clockwise.
477 * If the direction of the from->to is atni-clockwise, add 4
478 */
479 if (from->x > to->x) {
480 if (from->y > to->y) {
481 if (to->y >= centre->y)
482 quadrant = 1 + 4;
483 else
484 quadrant = 3;
485 } else {
486 if (to->x >= centre->x)
487 quadrant = 0 + 4;
488 else
489 quadrant = 2;
490 }
491 } else {
492 if (from->y > to->y) {
493 if (from->x >= centre->x)
494 quadrant = 0;
495 else
496 quadrant = 2 + 4;
497 } else {
498 if (from->x >= centre->x)
499 quadrant = 3 + 4;
500 else
501 quadrant = 1;
502 }
503 }
504
505 switch(quadrant) {
506 /* quadrant 0, arc goes clockwise */
507 case 0:
508 if (from->x == centre->x) {
509 from_alpha = M_PI / 2;
510 } else {
511 from_tan_alpha = (from->y - centre->y) / (from->x - centre->x);
512 from_alpha = atan(from_tan_alpha);
513 }
514 to_tan_alpha = (to->y - centre->y) / (to->x - centre->x);
515 to_alpha = atan(to_tan_alpha);
516
517 half_inscribed_angle = (from_alpha - to_alpha) / 2;
518 intersect_dist = radius / cos(half_inscribed_angle);
519 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
520
521 intersect_x = centre->x + cos(to_alpha + half_inscribed_angle) * intersect_dist;
522 intersect_y = centre->y + sin(to_alpha + half_inscribed_angle) * intersect_dist;
523 break;
524 /* quadrant 1, arc goes clockwise */
525 case 1:
526 from_tan_alpha = (from->y - centre->y) / (centre->x - from->x);
527 from_alpha = atan(from_tan_alpha);
528
529 if (to->x == centre->x) {
530 to_alpha = M_PI / 2;
531 } else {
532 to_tan_alpha = (to->y - centre->y) / (centre->x - to->x);
533 to_alpha = atan(to_tan_alpha);
534 }
535
536 half_inscribed_angle = (to_alpha - from_alpha) / 2;
537 intersect_dist = radius / cos(half_inscribed_angle);
538 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
539
540 intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;
541 intersect_y = centre->y + sin(from_alpha + half_inscribed_angle) * intersect_dist;
542 break;
543 /* quadrant 2, arc goes clockwise */
544 case 2:
545 if (from->x == centre->x) {
546 from_alpha = M_PI / 2;
547 } else {
548 from_tan_alpha = (centre->y - from->y) / (centre->x - from->x);
549 from_alpha = atan(from_tan_alpha);
550 }
551
552 to_tan_alpha = (centre->y - to->y) / (centre->x - to->x);
553 to_alpha = atan(to_tan_alpha);
554
555 half_inscribed_angle = (to_alpha - from_alpha) / 2;
556 intersect_dist = radius / cos(half_inscribed_angle);
557 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
558
559 intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;
560 intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;
561 break;
562 /* quadrant 3, arc goes clockwise */
563 case 3:
564 from_tan_alpha = (centre->y - from->y) / (from->x - centre->x);
565 from_alpha = atan(from_tan_alpha);
566
567 if (to->x == centre->x) {
568 to_alpha = M_PI / 2;
569 } else {
570 to_tan_alpha = (centre->y - to->y) / (to->x - centre->x);
571 to_alpha = atan(to_tan_alpha);
572 }
573
574 half_inscribed_angle = (to_alpha - from_alpha) / 2;
575 intersect_dist = radius / cos(half_inscribed_angle);
576 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
577
578 intersect_x = centre->x + cos(from_alpha + half_inscribed_angle) * intersect_dist;
579 intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;
580 break;
581 /* quadrant 0, arc goes anti-clockwise */
582 case 4:
583 from_tan_alpha = (from->y - centre->y) / (from->x - centre->x);
584 from_alpha = atan(from_tan_alpha);
585
586 if (to->y == centre->y)
587 to_alpha = M_PI / 2;
588 else {
589 to_tan_alpha = (to->y - centre->y) / (to->x - centre->x);
590 to_alpha = atan(to_tan_alpha);
591 }
592
593 half_inscribed_angle = (to_alpha - from_alpha) / 2;
594 intersect_dist = radius / cos(half_inscribed_angle);
595 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
596
597 intersect_x = centre->x + cos(from_alpha + half_inscribed_angle) * intersect_dist;
598 intersect_y = centre->y + sin(from_alpha + half_inscribed_angle) * intersect_dist;
599 break;
600 /* quadrant 1, arc goes anti-clockwise */
601 case 5:
602 from_tan_alpha = (centre->x - from->x) / (from->y - centre->y);
603 from_alpha = atan(from_tan_alpha);
604
605 if (to->y == centre->y) {
606 to_alpha = M_PI / 2;
607 }
608 else {
609 to_tan_alpha = (centre->x - to->x) / (to->y - centre->y);
610 to_alpha = atan(to_tan_alpha);
611 }
612
613 half_inscribed_angle = (to_alpha - from_alpha) / 2;
614 intersect_dist = radius / cos(half_inscribed_angle);
615 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
616
617 intersect_x = centre->x - sin(from_alpha + half_inscribed_angle) * intersect_dist;
618 intersect_y = centre->y + cos(from_alpha + half_inscribed_angle) * intersect_dist;
619 break;
620 /* quadrant 2, arc goes anti-clockwise */
621 case 6:
622 from_tan_alpha = (from->y - centre->y) / (centre->x - from->x);
623 from_alpha = atan(from_tan_alpha);
624
625 if (to->x == centre->x) {
626 to_alpha = M_PI / 2;
627 } else {
628 to_tan_alpha = (centre->y - to->y) / (centre->x - to->x);
629 to_alpha = atan(to_tan_alpha);
630 }
631
632 half_inscribed_angle = (to_alpha - from_alpha) / 2;
633 intersect_dist = radius / cos(half_inscribed_angle);
634 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
635
636 intersect_x = centre->x - cos(from_alpha + half_inscribed_angle) * intersect_dist;
637 intersect_y = centre->y - sin(from_alpha + half_inscribed_angle) * intersect_dist;
638 break;
639 /* quadrant 3, arc goes anti-clockwise */
640 case 7:
641 if (from->x == centre->x) {
642 from_alpha = M_PI / 2;
643 } else {
644 from_tan_alpha = (centre->y - from->y) / (from->x - centre->x);
645 from_alpha = atan(from_tan_alpha);
646 }
647 to_tan_alpha = (centre->y - to->y) / (to->x - centre->x);
648 to_alpha = atan(to_tan_alpha);
649
650 half_inscribed_angle = (from_alpha - to_alpha) / 2;
651 intersect_dist = radius / cos(half_inscribed_angle);
652 tangent_intersect_dist = tan(half_inscribed_angle) * radius;
653
654 intersect_x = centre->x + cos(to_alpha + half_inscribed_angle) * intersect_dist;
655 intersect_y = centre->y - sin(to_alpha + half_inscribed_angle) * intersect_dist;
656 break;
657 }
658
659 F = (4.0 / 3.0) / (1 + sqrt(1 + ((tangent_intersect_dist / radius) * (tangent_intersect_dist / radius))));
660
661 from_control->x = from->x - ((from->x - intersect_x) * F);
662 from_control->y = from->y - ((from->y - intersect_y) * F);
663 to_control->x = to->x - ((to->x - intersect_x) * F);
664 to_control->y = to->y - ((to->y - intersect_y) * F);
665
666 return 0;
667 }
668
669 /* Create a 'patch_curve' element whch is a straight line between two points */
patch_lineto(gs_matrix_fixed * ctm,gs_point * from,gs_point * to,patch_curve_t * p,float t)670 static int patch_lineto(gs_matrix_fixed *ctm, gs_point *from, gs_point *to, patch_curve_t *p, float t)
671 {
672 double x_1third, x_2third, y_1third, y_2third;
673
674 x_1third = (to->x - from->x) / 3;
675 x_2third = x_1third * 2;
676 y_1third = (to->y - from->y) / 3;
677 y_2third = y_1third * 2;
678
679 gs_point_transform2fixed(ctm, from->x, from->y, &p->vertex.p);
680 gs_point_transform2fixed(ctm, from->x + x_1third, from->y + y_1third, &p->control[0]);
681 gs_point_transform2fixed(ctm, from->x + x_2third, from->y + y_2third, &p->control[1]);
682
683 p->vertex.cc[0] = t;
684 p->vertex.cc[1] = t;
685 p->straight = 1;
686
687 return 0;
688 }
689
patch_curveto(gs_matrix_fixed * ctm,gs_point * centre,gs_point * from,gs_point * to,patch_curve_t * p,float t)690 static int patch_curveto(gs_matrix_fixed *ctm, gs_point *centre, gs_point *from, gs_point *to, patch_curve_t *p, float t)
691 {
692 gs_point from_control, to_control;
693
694 find_arc_control_points(from, to, &from_control, &to_control, centre);
695
696 gs_point_transform2fixed(ctm, from->x, from->y, &p->vertex.p);
697 gs_point_transform2fixed(ctm, from_control.x, from_control.y, &p->control[0]);
698 gs_point_transform2fixed(ctm, to_control.x, to_control.y, &p->control[1]);
699 p->vertex.cc[0] = t;
700 p->vertex.cc[1] = t;
701 p->straight = 0;
702
703 return 0;
704 }
705
draw_quarter_annulus(patch_fill_state_t * pfs,gs_point * centre,double radius,gs_point * corner,float t)706 static int draw_quarter_annulus(patch_fill_state_t *pfs, gs_point *centre, double radius, gs_point *corner, float t)
707 {
708 gs_point p0, p1, initial;
709 patch_curve_t p[4];
710 int code;
711
712 if (corner->x > centre->x) {
713 initial.x = centre->x + radius;
714 }
715 else {
716 initial.x = centre->x - radius;
717 }
718 initial.y = centre->y;
719
720 p1.x = initial.x;
721 p1.y = corner->y;
722 patch_lineto(&pfs->pgs->ctm, &initial, &p1, &p[0], t);
723 p0.x = centre->x;
724 p0.y = p1.y;
725 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &p[1], t);
726 p1.x = centre->x;
727 if (centre->y > corner->y) {
728 p1.y = centre->y - radius;
729 } else {
730 p1.y = centre->y + radius;
731 }
732 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &p[2], t);
733 patch_curveto(&pfs->pgs->ctm, centre, &p1, &initial, &p[3], t);
734 code = patch_fill(pfs, (const patch_curve_t *)&p, NULL, NULL);
735 if (code < 0)
736 return code;
737
738 if (corner->x > centre->x)
739 initial.x = corner->x - (corner->x - (centre->x + radius));
740 else
741 initial.x = centre->x - radius;
742 initial.y = corner->y;
743 patch_lineto(&pfs->pgs->ctm, corner, &initial, &p[0], t);
744
745 p0.x = initial.x;
746 p0.y = centre->y;
747 patch_lineto(&pfs->pgs->ctm, &initial, &p0, &p[1], t);
748
749 p1.y = p0.y;
750 p1.x = corner->x;
751 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &p[2], t);
752 patch_lineto(&pfs->pgs->ctm, &p1, corner, &p[3], t);
753
754 return (patch_fill(pfs, (const patch_curve_t *)&p, NULL, NULL));
755 }
756
R_tensor_annulus_extend_tangent(patch_fill_state_t * pfs,double x0,double y0,double r0,double t0,double x1,double y1,double r1,double t1,double r2)757 static int R_tensor_annulus_extend_tangent(patch_fill_state_t *pfs,
758 double x0, double y0, double r0, double t0,
759 double x1, double y1, double r1, double t1, double r2)
760 {
761 patch_curve_t curve[4];
762 gs_point p0, p1;
763 int code = 0, q = 0;
764
765 /* special case axis aligned circles. Its quicker to handle these specially as it
766 * avoid lots of trigonometry in the general case code, and avoids us
767 * having to watch out for infinity as the result of tan() operations.
768 */
769 if (x0 == x1 || y0 == y1) {
770 if (x0 == x1 && y0 > y1) {
771 /* tangent at top of circles */
772 p0.x = x1, p0.y = y1;
773 p1.x = x1 + r2, p1.y = y1 - r2;
774 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
775 p1.x = x1 - r2, p1.y = y1 - r2;
776 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
777 p1.x = x1 + r2, p1.y = y1 + r1;
778 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
779 p1.x = x1 - r2, p1.y = y1 + r1;
780 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
781 }
782 if (x0 == x1 && y0 < y1) {
783 /* tangent at bottom of circles */
784 p0.x = x1, p0.y = y1;
785 p1.x = x1 + r2, p1.y = y1 + r2;
786 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
787 p1.x = x1 - r2, p1.y = y1 + r2;
788 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
789 p1.x = x1 + r2, p1.y = y1 - r1;
790 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
791 p1.x = x1 - r2, p1.y = y1 - r1;
792 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
793 }
794 if (y0 == y1 && x0 > x1) {
795 /* tangent at right of circles */
796 p0.x = x1, p0.y = y1;
797 p1.x = x1 - r2, p1.y = y1 - r2;
798 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
799 p1.x = x1 - r2, p1.y = y1 + r2;
800 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
801 p1.x = x1 + r1, p1.y = y1 + r2;
802 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
803 p1.x = x1 + r1, p1.y = y1 - r2;
804 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
805 }
806 if (y0 == y1 && x0 < x1) {
807 /* tangent at left of circles */
808 p0.x = x1, p0.y = y1;
809 p1.x = x1 + r2, p1.y = y1 - r2;
810 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
811 p1.x = x1 + r2, p1.y = y1 + r2;
812 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
813 p1.x = x1 - r1, p1.y = y1 + r2;
814 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
815 p1.x = x1 - r1, p1.y = y1 - r2;
816 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
817 }
818 }
819 else {
820 double tx, ty, endx, endy, intersectx, intersecty, alpha, sinalpha, cosalpha, tanalpha;
821 gs_point centre;
822
823 /* First lets figure out which quadrant the smaller circle is in (we always
824 * get called to fill from the larger circle), x0, y0, r0 is the smaller circle.
825 */
826 if (x0 < x1) {
827 if (y0 < y1)
828 q = 2;
829 else
830 q = 1;
831 } else {
832 if (y0 < y1)
833 q = 3;
834 else
835 q = 0;
836 }
837 switch(q) {
838 case 0:
839 /* We have two four-sided elements, from the tangent point
840 * each side, to the point where the tangent crosses an
841 * axis of the larger circle. A line back to the edge
842 * of the larger circle, a line to the point where an axis
843 * crosses the smaller circle, then an arc back to the starting point.
844 */
845 /* Figure out the tangent point */
846 /* sin (angle) = y1 - y0 / r1 - r0
847 * ty = ((y1 - y0) / (r1 - r0)) * r1
848 */
849 ty = y1 + ((y0 - y1) / (r1 - r0)) * r1;
850 tx = x1 + ((x0 - x1) / (r1 - r0)) * r1;
851 /* Now actually calculating the point where the tangent crosses the axis of the larger circle
852 * So we need to know the angle the tangent makes with the axis of the smaller circle
853 * as its the same angle where it crosses the axis of the larger circle.
854 * We know the centres and the tangent are co-linear, so sin(a) = y0 - y1 / r1 - r0
855 * We know the tangent is r1 from the centre of the larger circle, so the hypotenuse
856 * is r0 / cos(a). That gives us 'x' and we already know y as its the centre of the larger
857 * circle
858 */
859 sinalpha = (y0 - y1) / (r1 - r0);
860 alpha = asin(sinalpha);
861 cosalpha = cos(alpha);
862 intersectx = x1 + (r1 / cosalpha);
863 intersecty = y1;
864
865 p0.x = tx, p0.y = ty;
866 p1.x = tx + (intersectx - tx) / 2, p1.y = ty - (ty - intersecty) / 2;
867 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
868 p0.x = intersectx, p0.y = intersecty;
869 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
870 p1.x = x1 + r1, p1.y = y1;
871 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
872 p0.x = tx, p0.y = ty;
873 centre.x = x1, centre.y = y1;
874 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
875 code = patch_fill(pfs, curve, NULL, NULL);
876 if (code < 0)
877 return code;
878
879 if (intersectx < x1 + r2) {
880 /* didn't get all the way to the edge, quadrant 3 is composed of 2 quads :-(
881 * An 'annulus' where the right edge is less than the normal extent and a
882 * quad which is a rectangle with one corner chopped of at an angle.
883 */
884 p0.x = x1, p0.y = y1;
885 p1.x = intersectx, p1.y = y1 - r2;
886 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
887 endx = x1 + r2;
888 endy = y1 - (tan ((M_PI / 2) - alpha)) * (endx - intersectx);
889 p0.x = intersectx, p0.y = y1;
890 p1.x = x1 + r2, p1.y = endy;
891 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
892 p0.x = x1 + r2, p0.y = y0 - r2;
893 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
894 p1.x = intersectx, p1.y = p0.y;
895 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
896 p0.x = intersectx, p0.y = y1;
897 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
898 code = patch_fill(pfs, curve, NULL, NULL);
899 if (code < 0)
900 return code;
901
902 } else {
903 /* Quadrant 3 is a normal quarter annulua */
904 p0.x = x1, p0.y = y1;
905 p1.x = x1 + r2, p1.y = y1 - r2;
906 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
907 }
908
909 /* Q2 is always a full annulus... */
910 p0.x = x1, p0.y = y1;
911 p1.x = x1 - r2, p1.y = y1 - r2;
912 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
913
914 /* alpha is now the angle between the x axis and the tangent to the
915 * circles.
916 */
917 alpha = (M_PI / 2) - alpha;
918 cosalpha = cos(alpha);
919 endy = y1 + (r1 / cosalpha);
920 endx = x1;
921
922 p0.x = tx, p0.y = ty;
923 p1.x = endx - ((endx - tx) / 2), p1.y = endy - ((endy - ty) / 2);
924 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
925 p0.x = endx, p0.y = endy;
926 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
927 p1.x = x1, p1.y = y1 + r1;
928 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
929 p0.x = tx, p0.y = ty;
930 centre.x = x1, centre.y = y1;
931 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
932 code = patch_fill(pfs, curve, NULL, NULL);
933 if (code < 0)
934 return code;
935
936 /* Q1 is simimlar to Q3, either a full quarter annulus
937 * or a partial one, depending on where the tangent crosses
938 * the y axis
939 */
940 tanalpha = tan(alpha);
941 intersecty = y1 + tanalpha * (r2 + (intersectx - x1));
942 intersectx = x1 - r2;
943
944 if (endy < y1 + r2) {
945 /* didn't get all the way to the edge, quadrant 1 is composed of 2 quads :-(
946 * An 'annulus' where the right edge is less than the normal extent and a
947 * quad which is a rectangle with one corner chopped of at an angle.
948 */
949 p0.x = x1, p0.y = y1;
950 p1.x = x1 - r2, p1.y = endy;
951 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
952 p0.x = x1, p0.y = y1 + r1;
953 p1.x = x1, p1.y = endy;
954 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
955 p0.x = x1 - r2, p0.y = intersecty;
956 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
957 p1.x = p0.x, p1.y = y1 + r1;
958 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
959 p0.x = x1, p0.y = y1 + r1;
960 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
961 code = patch_fill(pfs, curve, NULL, NULL);
962 if (code < 0)
963 return code;
964 } else {
965 /* Quadrant 1 is a normal quarter annulua */
966 p0.x = x1, p0.y = y1;
967 p1.x = x1 - r2, p1.y = y1 + r2;
968 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
969 }
970 break;
971 case 1:
972 /* We have two four-sided elements, from the tangent point
973 * each side, to the point where the tangent crosses an
974 * axis of the larger circle. A line back to the edge
975 * of the larger circle, a line to the point where an axis
976 * crosses the smaller circle, then an arc back to the starting point.
977 */
978 /* Figure out the tangent point */
979 /* sin (angle) = y1 - y0 / r1 - r0
980 * ty = ((y1 - y0) / (r1 - r0)) * r1
981 */
982 ty = y1 + ((y0 - y1) / (r1 - r0)) * r1;
983 tx = x1 - ((x1 - x0) / (r1 - r0)) * r1;
984 /* Now actually calculating the point where the tangent crosses the axis of the larger circle
985 * So we need to know the angle the tangent makes with the axis of the smaller circle
986 * as its the same angle where it crosses the axis of the larger circle.
987 * We know the centres and the tangent are co-linear, so sin(a) = y0 - y1 / r1 - r0
988 * We know the tangent is r1 from the centre of the larger circle, so the hypotenuse
989 * is r0 / cos(a). That gives us 'x' and we already know y as its the centre of the larger
990 * circle
991 */
992 sinalpha = (y0 - y1) / (r1 - r0);
993 alpha = asin(sinalpha);
994 cosalpha = cos(alpha);
995 intersectx = x1 - (r1 / cosalpha);
996 intersecty = y1;
997
998 p0.x = tx, p0.y = ty;
999 p1.x = tx - (tx - intersectx) / 2, p1.y = ty - (ty - intersecty) / 2;
1000 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1001 p0.x = intersectx, p0.y = intersecty;
1002 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1003 p1.x = x1 - r1, p1.y = y1;
1004 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1005 p0.x = tx, p0.y = ty;
1006 centre.x = x1, centre.y = y1;
1007 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1008 code = patch_fill(pfs, curve, NULL, NULL);
1009 if (code < 0)
1010 return code;
1011
1012 if (intersectx > x1 - r2) {
1013 /* didn't get all the way to the edge, quadrant 2 is composed of 2 quads :-(
1014 * An 'annulus' where the right edge is less than the normal extent and a
1015 * quad which is a rectangle with one corner chopped of at an angle.
1016 */
1017 p0.x = x1, p0.y = y1;
1018 p1.x = intersectx, p1.y = y1 - r2;
1019 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1020 endx = x1 - r2;
1021 endy = y1 - (tan ((M_PI / 2) - alpha)) * (intersectx - endx);
1022 p0.x = intersectx, p0.y = y1;
1023 p1.x = x1 - r2, p1.y = endy;
1024 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1025 p0.x = x1 - r2, p0.y = y0 - r2;
1026 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1027 p1.x = intersectx, p1.y = p0.y;
1028 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1029 p0.x = intersectx, p0.y = y1;
1030 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1031 code = patch_fill(pfs, curve, NULL, NULL);
1032 if (code < 0)
1033 return code;
1034
1035 } else {
1036 /* Quadrant 2 is a normal quarter annulua */
1037 p0.x = x1, p0.y = y1;
1038 p1.x = x1 - r2, p1.y = y1 - r2;
1039 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1040 }
1041
1042 /* Q3 is always a full annulus... */
1043 p0.x = x1, p0.y = y1;
1044 p1.x = x1 + r2, p1.y = y1 - r2;
1045 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1046
1047 /* alpha is now the angle between the x axis and the tangent to the
1048 * circles.
1049 */
1050 alpha = (M_PI / 2) - alpha;
1051 cosalpha = cos(alpha);
1052 endy = y1 + (r1 / cosalpha);
1053 endx = x1;
1054
1055 p0.x = tx, p0.y = ty;
1056 p1.x = endx + ((tx - endx) / 2), p1.y = endy - ((endy - ty) / 2);
1057 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1058 p0.x = endx, p0.y = endy;
1059 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1060 p1.x = x1, p1.y = y1 + r1;
1061 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1062 p0.x = tx, p0.y = ty;
1063 centre.x = x1, centre.y = y1;
1064 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1065 code = patch_fill(pfs, curve, NULL, NULL);
1066 if (code < 0)
1067 return code;
1068
1069 /* Q0 is simimlar to Q2, either a full quarter annulus
1070 * or a partial one, depending on where the tangent crosses
1071 * the y axis
1072 */
1073 tanalpha = tan(alpha);
1074 intersecty = y1 + tanalpha * (r2 + (x1 - intersectx));
1075 intersectx = x1 + r2;
1076
1077 if (endy < y1 + r2) {
1078 /* didn't get all the way to the edge, quadrant 0 is composed of 2 quads :-(
1079 * An 'annulus' where the right edge is less than the normal extent and a
1080 * quad which is a rectangle with one corner chopped of at an angle.
1081 */
1082 p0.x = x1, p0.y = y1;
1083 p1.x = x1 + r2, p1.y = endy;
1084 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1085 p0.x = x1, p0.y = y1 + r1;
1086 p1.x = x1, p1.y = endy;
1087 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1088 p0.x = x1 + r2, p0.y = intersecty;
1089 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1090 p1.x = p0.x, p1.y = y1 + r1;
1091 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1092 p0.x = x1, p0.y = y1 + r1;
1093 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1094 code = patch_fill(pfs, curve, NULL, NULL);
1095 if (code < 0)
1096 return code;
1097 } else {
1098 /* Quadrant 0 is a normal quarter annulua */
1099 p0.x = x1, p0.y = y1;
1100 p1.x = x1 + r2, p1.y = y1 + r2;
1101 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1102 }
1103 break;
1104 case 2:
1105 /* We have two four-sided elements, from the tangent point
1106 * each side, to the point where the tangent crosses an
1107 * axis of the larger circle. A line back to the edge
1108 * of the larger circle, a line to the point where an axis
1109 * crosses the smaller circle, then an arc back to the starting point.
1110 */
1111 /* Figure out the tangent point */
1112 /* sin (angle) = y1 - y0 / r1 - r0
1113 * ty = ((y1 - y0) / (r1 - r0)) * r1
1114 */
1115 ty = y1 - ((y1 - y0) / (r1 - r0)) * r1;
1116 tx = x1 - ((x1 - x0) / (r1 - r0)) * r1;
1117 /* Now actually calculating the point where the tangent crosses the axis of the larger circle
1118 * So we need to know the angle the tangent makes with the axis of the smaller circle
1119 * as its the same angle where it crosses the axis of the larger circle.
1120 * We know the centres and the tangent are co-linear, so sin(a) = y0 - y1 / r1 - r0
1121 * We know the tangent is r1 from the centre of the larger circle, so the hypotenuse
1122 * is r0 / cos(a). That gives us 'x' and we already know y as its the centre of the larger
1123 * circle
1124 */
1125 sinalpha = (y1 - y0) / (r1 - r0);
1126 alpha = asin(sinalpha);
1127 cosalpha = cos(alpha);
1128 intersectx = x1 - (r1 / cosalpha);
1129 intersecty = y1;
1130
1131 p0.x = tx, p0.y = ty;
1132 p1.x = tx + (intersectx - tx) / 2, p1.y = ty - (ty - intersecty) / 2;
1133 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1134 p0.x = intersectx, p0.y = intersecty;
1135 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1136 p1.x = x1 - r1, p1.y = y1;
1137 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1138 p0.x = tx, p0.y = ty;
1139 centre.x = x1, centre.y = y1;
1140 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1141 code = patch_fill(pfs, curve, NULL, NULL);
1142 if (code < 0)
1143 return code;
1144 if (intersectx > x1 - r2) {
1145 /* didn't get all the way to the edge, quadrant 1 is composed of 2 quads :-(
1146 * An 'annulus' where the right edge is less than the normal extent and a
1147 * quad which is a rectangle with one corner chopped of at an angle.
1148 */
1149 p0.x = x1, p0.y = y1;
1150 p1.x = intersectx, p1.y = y1 + r2;
1151 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1152 endy = y1+r2;
1153 endx = intersectx - ((endy - intersecty) / (tan ((M_PI / 2) - alpha)));
1154 p0.x = intersectx, p0.y = y1;
1155 p1.x = endx, p1.y = endy;
1156 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1157 p0.x = x1 - r1, p0.y = endy;
1158 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1159 p1.x = x1 - r1, p1.y = y1;
1160 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1161 p0.x = intersectx, p0.y = y1;
1162 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1163 code = patch_fill(pfs, curve, NULL, NULL);
1164 if (code < 0)
1165 return code;
1166 } else {
1167 /* Quadrant 1 is a normal quarter annulua */
1168 p0.x = x1, p0.y = y1;
1169 p1.x = x1 - r2, p1.y = y1 + r2;
1170 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1171 }
1172
1173 /* Q0 is always a full annulus... */
1174 p0.x = x1, p0.y = y1;
1175 p1.x = x1 + r2, p1.y = y1 + r2;
1176 if (p1.y < 0)
1177 p1.y = 0;
1178 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1179
1180 /* alpha is now the angle between the x axis and the tangent to the
1181 * circles.
1182 */
1183 alpha = (M_PI / 2) - alpha;
1184 cosalpha = cos(alpha);
1185 endy = y1 - (r1 / cosalpha);
1186 endx = x1;
1187
1188 p0.x = tx, p0.y = ty;
1189 p1.x = endx + ((endx - tx) / 2), p1.y = endy - ((ty - endy) / 2);
1190 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1191 p0.x = endx, p0.y = endy;
1192 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1193 p1.x = x1, p1.y = y1 - r1;
1194 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1195 p0.x = tx, p0.y = ty;
1196 centre.x = x1, centre.y = y1;
1197 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1198 code = patch_fill(pfs, curve, NULL, NULL);
1199 if (code < 0)
1200 return code;
1201
1202 /* Q3 is simimlar to Q1, either a full quarter annulus
1203 * or a partial one, depending on where the tangent crosses
1204 * the y axis
1205 */
1206 tanalpha = tan(alpha);
1207 intersecty = y1 - tanalpha * (r2 + (x1 - intersectx));
1208 intersectx = x1 + r2;
1209
1210 if (endy > y1 - r2) {
1211 /* didn't get all the way to the edge, quadrant 3 is composed of 2 quads :-(
1212 * An 'annulus' where the right edge is less than the normal extent and a
1213 * quad which is a rectangle with one corner chopped of at an angle.
1214 */
1215 p0.x = x1, p0.y = y1;
1216 p1.x = x1 + r2, p1.y = endy;
1217 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1218 p0.x = x1, p0.y = y1 - r1;
1219 p1.x = x1, p1.y = endy;
1220 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1221 p0.x = x1 + r2, p0.y = intersecty;
1222 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1223 p1.x = p0.x, p1.y = y1 - r1;
1224 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1225 p0.x = x1, p0.y = y1 - r1;
1226 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1227 code = patch_fill(pfs, curve, NULL, NULL);
1228 if (code < 0)
1229 return code;
1230 } else {
1231 /* Quadrant 1 is a normal quarter annulua */
1232 p0.x = x1, p0.y = y1;
1233 p1.x = x1 + r2, p1.y = y1 - r2;
1234 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1235 }
1236 break;
1237 case 3:
1238 /* We have two four-sided elements, from the tangent point
1239 * each side, to the point where the tangent crosses an
1240 * axis of the larger circle. A line back to the edge
1241 * of the larger circle, a line to the point where an axis
1242 * crosses the smaller circle, then an arc back to the starting point.
1243 */
1244 /* Figure out the tangent point */
1245 /* sin (angle) = y1 - y0 / r1 - r0
1246 * ty = ((y1 - y0) / (r1 - r0)) * r1
1247 */
1248 ty = y1 - ((y1 - y0) / (r1 - r0)) * r1;
1249 tx = x1 + ((x0 - x1) / (r1 - r0)) * r1;
1250 /* Now actually calculating the point where the tangent crosses the axis of the larger circle
1251 * So we need to know the angle the tangent makes with the axis of the smaller circle
1252 * as its the same angle where it crosses the axis of the larger circle.
1253 * We know the centres and the tangent are co-linear, so sin(a) = y0 - y1 / r1 - r0
1254 * We know the tangent is r1 from the centre of the larger circle, so the hypotenuse
1255 * is r0 / cos(a). That gives us 'x' and we already know y as its the centre of the larger
1256 * circle
1257 */
1258 sinalpha = (y1 - y0) / (r1 - r0);
1259 alpha = asin(sinalpha);
1260 cosalpha = cos(alpha);
1261 intersectx = x1 + (r1 / cosalpha);
1262 intersecty = y1;
1263
1264 p0.x = tx, p0.y = ty;
1265 p1.x = tx + (intersectx - tx) / 2, p1.y = ty + (intersecty - ty) / 2;
1266 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1267 p0.x = intersectx, p0.y = intersecty;
1268 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1269 p1.x = x1 + r1, p1.y = y1;
1270 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1271 p0.x = tx, p0.y = ty;
1272 centre.x = x1, centre.y = y1;
1273 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1274 code = patch_fill(pfs, curve, NULL, NULL);
1275 if (code < 0)
1276 return code;
1277 if (intersectx < x1 + r2) {
1278 /* didn't get all the way to the edge, quadrant 0 is composed of 2 quads :-(
1279 * An 'annulus' where the right edge is less than the normal extent and a
1280 * quad which is a rectangle with one corner chopped of at an angle.
1281 */
1282 p0.x = x1, p0.y = y1;
1283 p1.x = intersectx, p1.y = y1 + r2;
1284 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1285 endy = y1 + r2;
1286 endx = intersectx + ((endy - intersecty) / (tan ((M_PI / 2) - alpha)));
1287 p0.x = intersectx, p0.y = y1;
1288 p1.x = endx, p1.y = endy;
1289 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1290 p0.x = x1 + r1, p0.y = endy;
1291 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1292 p1.x = x1 + r1, p1.y = y1;
1293 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1294 p0.x = intersectx, p0.y = y1;
1295 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1296 code = patch_fill(pfs, curve, NULL, NULL);
1297 if (code < 0)
1298 return code;
1299
1300 } else {
1301 /* Quadrant 0 is a normal quarter annulua */
1302 p0.x = x1, p0.y = y1;
1303 p1.x = x1 + r2, p1.y = y1 + r2;
1304 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1305 }
1306 /* Q1 is always a full annulus... */
1307 p0.x = x1, p0.y = y1;
1308 p1.x = x1 - r2, p1.y = y1 + r2;
1309 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1310
1311 /* alpha is now the angle between the x axis and the tangent to the
1312 * circles.
1313 */
1314 alpha = (M_PI / 2) - alpha;
1315 cosalpha = cos(alpha);
1316 endy = y1 - (r1 / cosalpha);
1317 endx = x1;
1318
1319 p0.x = tx, p0.y = ty;
1320 p1.x = endx + ((tx - endx) / 2), p1.y = endy + ((ty - endy) / 2);
1321 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1322 p0.x = endx, p0.y = endy;
1323 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1324 p1.x = x1, p1.y = y1 - r1;
1325 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1326 p0.x = tx, p0.y = ty;
1327 centre.x = x1, centre.y = y1;
1328 patch_curveto(&pfs->pgs->ctm, ¢re, &p1, &p0, &curve[3], t0);
1329 code = patch_fill(pfs, curve, NULL, NULL);
1330 if (code < 0)
1331 return code;
1332
1333 /* Q3 is simimlar to Q1, either a full quarter annulus
1334 * or a partial one, depending on where the tangent crosses
1335 * the y axis
1336 */
1337 tanalpha = tan(alpha);
1338 intersecty = y1 - tanalpha * (r2 + (intersectx - x1));
1339 intersectx = x1 - r2;
1340
1341 if (endy > y1 - r2) {
1342 /* didn't get all the way to the edge, quadrant 3 is composed of 2 quads :-(
1343 * An 'annulus' where the right edge is less than the normal extent and a
1344 * quad which is a rectangle with one corner chopped of at an angle.
1345 */
1346 p0.x = x1, p0.y = y1;
1347 p1.x = x1 - r2, p1.y = endy;
1348 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1349 p0.x = x1, p0.y = y1 - r1;
1350 p1.x = x1, p1.y = endy;
1351 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[0], t0);
1352 p0.x = x1 - r2, p0.y = intersecty;
1353 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[1], t0);
1354 p1.x = p0.x, p1.y = y1 - r1;
1355 patch_lineto(&pfs->pgs->ctm, &p0, &p1, &curve[2], t0);
1356 p0.x = x1, p0.y = y1 - r1;
1357 patch_lineto(&pfs->pgs->ctm, &p1, &p0, &curve[3], t0);
1358 code = patch_fill(pfs, curve, NULL, NULL);
1359 if (code < 0)
1360 return code;
1361 } else {
1362 /* Quadrant 1 is a normal quarter annulua */
1363 p0.x = x1, p0.y = y1;
1364 p1.x = x1 - r2, p1.y = y1 - r2;
1365 draw_quarter_annulus(pfs, &p0, r1, &p1, t0);
1366 }
1367 break;
1368 }
1369 }
1370 return 0;
1371 }
1372
1373 static int
R_outer_circle(patch_fill_state_t * pfs,const gs_rect * rect,double x0,double y0,double r0,double x1,double y1,double r1,double * x2,double * y2,double * r2)1374 R_outer_circle(patch_fill_state_t *pfs, const gs_rect *rect,
1375 double x0, double y0, double r0,
1376 double x1, double y1, double r1,
1377 double *x2, double *y2, double *r2)
1378 {
1379 double dx = x1 - x0, dy = y1 - y0;
1380 double sp, sq, s;
1381
1382 /* Compute a cone circle, which contacts the rect externally. */
1383 /* Don't bother with all 4 sides of the rect,
1384 just do with the X or Y span only,
1385 so it's not an exact contact, sorry. */
1386 if (any_abs(dx) > any_abs(dy)) {
1387 /* Solving :
1388 x0 + (x1 - x0) * sq + r0 + (r1 - r0) * sq == bbox_px
1389 (x1 - x0) * sp + (r1 - r0) * sp == bbox_px - x0 - r0
1390 sp = (bbox_px - x0 - r0) / (x1 - x0 + r1 - r0)
1391
1392 x0 + (x1 - x0) * sq - r0 - (r1 - r0) * sq == bbox_qx
1393 (x1 - x0) * sq - (r1 - r0) * sq == bbox_x - x0 + r0
1394 sq = (bbox_x - x0 + r0) / (x1 - x0 - r1 + r0)
1395 */
1396 if (x1 - x0 + r1 - r0 == 0) /* We checked for obtuse cone. */
1397 return_error(gs_error_unregistered); /* Must not happen. */
1398 if (x1 - x0 - r1 + r0 == 0) /* We checked for obtuse cone. */
1399 return_error(gs_error_unregistered); /* Must not happen. */
1400 sp = (rect->p.x - x0 - r0) / (x1 - x0 + r1 - r0);
1401 sq = (rect->q.x - x0 + r0) / (x1 - x0 - r1 + r0);
1402 } else {
1403 /* Same by Y. */
1404 if (y1 - y0 + r1 - r0 == 0) /* We checked for obtuse cone. */
1405 return_error(gs_error_unregistered); /* Must not happen. */
1406 if (y1 - y0 - r1 + r0 == 0) /* We checked for obtuse cone. */
1407 return_error(gs_error_unregistered); /* Must not happen. */
1408 sp = (rect->p.y - y0 - r0) / (y1 - y0 + r1 - r0);
1409 sq = (rect->q.y - y0 + r0) / (y1 - y0 - r1 + r0);
1410 }
1411 if (sp >= 1 && sq >= 1)
1412 s = max(sp, sq);
1413 else if(sp >= 1)
1414 s = sp;
1415 else if (sq >= 1)
1416 s = sq;
1417 else {
1418 /* The circle 1 is outside the rect, use it. */
1419 s = 1;
1420 }
1421 if (r0 + (r1 - r0) * s < 0) {
1422 /* Passed the cone apex, use the apex. */
1423 s = r0 / (r0 - r1);
1424 *r2 = 0;
1425 } else
1426 *r2 = r0 + (r1 - r0) * s;
1427 *x2 = x0 + (x1 - x0) * s;
1428 *y2 = y0 + (y1 - y0) * s;
1429 return 0;
1430 }
1431
1432 static double
R_rect_radius(const gs_rect * rect,double x0,double y0)1433 R_rect_radius(const gs_rect *rect, double x0, double y0)
1434 {
1435 double d, dd;
1436
1437 dd = hypot(rect->p.x - x0, rect->p.y - y0);
1438 d = hypot(rect->p.x - x0, rect->q.y - y0);
1439 dd = max(dd, d);
1440 d = hypot(rect->q.x - x0, rect->q.y - y0);
1441 dd = max(dd, d);
1442 d = hypot(rect->q.x - x0, rect->p.y - y0);
1443 dd = max(dd, d);
1444 return dd;
1445 }
1446
1447 static int
R_fill_triangle_new(patch_fill_state_t * pfs,const gs_rect * rect,double x0,double y0,double x1,double y1,double x2,double y2,double t)1448 R_fill_triangle_new(patch_fill_state_t *pfs, const gs_rect *rect,
1449 double x0, double y0, double x1, double y1, double x2, double y2, double t)
1450 {
1451 shading_vertex_t p0, p1, p2;
1452 patch_color_t *c;
1453 int code;
1454 reserve_colors(pfs, &c, 1); /* Can't fail */
1455
1456 p0.c = c;
1457 p1.c = c;
1458 p2.c = c;
1459 code = gs_point_transform2fixed(&pfs->pgs->ctm, x0, y0, &p0.p);
1460 if (code >= 0)
1461 code = gs_point_transform2fixed(&pfs->pgs->ctm, x1, y1, &p1.p);
1462 if (code >= 0)
1463 code = gs_point_transform2fixed(&pfs->pgs->ctm, x2, y2, &p2.p);
1464 if (code >= 0) {
1465 c->t[0] = c->t[1] = t;
1466 patch_resolve_color(c, pfs);
1467 code = mesh_triangle(pfs, &p0, &p1, &p2);
1468 }
1469 release_colors(pfs, pfs->color_stack, 1);
1470 return code;
1471 }
1472
1473 static int
R_obtuse_cone(patch_fill_state_t * pfs,const gs_rect * rect,double x0,double y0,double r0,double x1,double y1,double r1,double t0,double r_rect)1474 R_obtuse_cone(patch_fill_state_t *pfs, const gs_rect *rect,
1475 double x0, double y0, double r0,
1476 double x1, double y1, double r1, double t0, double r_rect)
1477 {
1478 double dx = x1 - x0, dy = y1 - y0, dr = any_abs(r1 - r0);
1479 double d = hypot(dx, dy);
1480 /* Assuming dr > d / 3 && d > dr + 1e-7 * (d + dr), see the caller. */
1481 double r = r_rect * 1.4143; /* A few bigger than sqrt(2). */
1482 double ax, ay, as; /* Cone apex. */
1483 double g0; /* The distance from apex to the tangent point of the 0th circle. */
1484 int code;
1485
1486 as = r0 / (r0 - r1);
1487 ax = x0 + (x1 - x0) * as;
1488 ay = y0 + (y1 - y0) * as;
1489 g0 = sqrt(dx * dx + dy * dy - dr * dr) * as;
1490 if (g0 < 1e-7 * r0) {
1491 /* Nearly degenerate, replace with half-plane. */
1492 /* Restrict the half plane with triangle, which covers the rect. */
1493 gs_point p0, p1, p2; /* Right tangent limit, apex limit, left tangent linit,
1494 (right, left == when looking from the apex). */
1495
1496 p0.x = ax - dy * r / d;
1497 p0.y = ay + dx * r / d;
1498 p1.x = ax - dx * r / d;
1499 p1.y = ay - dy * r / d;
1500 p2.x = ax + dy * r / d;
1501 p2.y = ay - dx * r / d;
1502 /* Split into 2 triangles at the apex,
1503 so that the apex is preciselly covered.
1504 Especially important when it is not exactly degenerate. */
1505 code = R_fill_triangle_new(pfs, rect, ax, ay, p0.x, p0.y, p1.x, p1.y, t0);
1506 if (code < 0)
1507 return code;
1508 return R_fill_triangle_new(pfs, rect, ax, ay, p1.x, p1.y, p2.x, p2.y, t0);
1509 } else {
1510 /* Compute the "limit" circle so that its
1511 tangent points are outside the rect. */
1512 /* Note: this branch is executed when the condition above is false :
1513 g0 >= 1e-7 * r0 .
1514 We believe that computing this branch with doubles
1515 provides enough precision after converting coordinates into 'fixed',
1516 and that the limit circle radius is not dramatically big.
1517 */
1518 double es, er; /* The limit circle parameter, radius. */
1519 double ex, ey; /* The limit circle centrum. */
1520
1521 es = as - as * r / g0; /* Always negative. */
1522 er = r * r0 / g0 ;
1523 ex = x0 + dx * es;
1524 ey = y0 + dy * es;
1525 /* Fill the annulus: */
1526 code = R_tensor_annulus(pfs, x0, y0, r0, t0, ex, ey, er, t0);
1527 if (code < 0)
1528 return code;
1529 /* Fill entire ending circle to ensure entire rect is covered. */
1530 return R_tensor_annulus(pfs, ex, ey, er, t0, ex, ey, 0, t0);
1531 }
1532 }
1533
1534 static int
R_tensor_cone_apex(patch_fill_state_t * pfs,const gs_rect * rect,double x0,double y0,double r0,double x1,double y1,double r1,double t)1535 R_tensor_cone_apex(patch_fill_state_t *pfs, const gs_rect *rect,
1536 double x0, double y0, double r0,
1537 double x1, double y1, double r1, double t)
1538 {
1539 double as = r0 / (r0 - r1);
1540 double ax = x0 + (x1 - x0) * as;
1541 double ay = y0 + (y1 - y0) * as;
1542
1543 return R_tensor_annulus(pfs, x1, y1, r1, t, ax, ay, 0, t);
1544 }
1545
1546 /*
1547 * A map of this code:
1548 *
1549 * R_extensions
1550 * |-> (R_rect_radius)
1551 * |-> (R_outer_circle)
1552 * |-> R_obtuse_cone
1553 * | |-> R_fill_triangle_new
1554 * | | '-> mesh_triangle
1555 * | | '-> mesh_triangle_rec <--.
1556 * | | |--------------------'
1557 * | | |-> small_mesh_triangle
1558 * | | | '-> fill_triangle
1559 * | | | '-> triangle_by_4 <--.
1560 * | | | |----------------'
1561 * | | | |-> constant_color_triangle
1562 * | | | |-> make_wedge_median (etc)
1563 * | | '-----------+--------------------.
1564 * | '-------------------. |
1565 * |-> R_tensor_cone_apex | |
1566 * | '-------------------+ |
1567 * '-> R_tensor_annulus <--' \|/
1568 * |-> (make_quadrant_arc) |
1569 * '-> patch_fill |
1570 * |-> fill_patch <--. |
1571 * | |-------------' |
1572 * | |------------------------------------+
1573 * | '-> fill_stripe |
1574 * | |-----------------------. |
1575 * | \|/ | |
1576 * |-> fill_wedges | |
1577 * '-> fill_wedges_aux <--. | |
1578 * |------------------' \|/ |
1579 * |----------------> mesh_padding '
1580 * | '----------------------------------.
1581 * '-> wedge_by_triangles <--. . |
1582 * |---------------------' | |
1583 * '-> fill_triangle_wedge <----' |
1584 * '-> fill_triangle_wedge_aux |
1585 * '-> fill_wedge_trap |
1586 * '-> wedge_trap_decompose |
1587 * '-> linear_color_trapezoid |
1588 * '-> decompose_linear_color <--|
1589 * |-------------------------'
1590 * '-> constant_color_trapezoid
1591 */
1592 static int
R_extensions(patch_fill_state_t * pfs,const gs_shading_R_t * psh,const gs_rect * rect,double t0,double t1,bool Extend0,bool Extend1)1593 R_extensions(patch_fill_state_t *pfs, const gs_shading_R_t *psh, const gs_rect *rect,
1594 double t0, double t1, bool Extend0, bool Extend1)
1595 {
1596 float x0 = psh->params.Coords[0], y0 = psh->params.Coords[1];
1597 double r0 = psh->params.Coords[2];
1598 float x1 = psh->params.Coords[3], y1 = psh->params.Coords[4];
1599 double r1 = psh->params.Coords[5];
1600 double dx = x1 - x0, dy = y1 - y0, dr = any_abs(r1 - r0);
1601 double d = hypot(dx, dy), r;
1602 int code;
1603
1604 /* In order for the circles to be nested, one end circle
1605 * needs to be sufficiently large to cover the entirety
1606 * of the other end circle. i.e.
1607 *
1608 * max(r0,r1) >= d + min(r0,r1)
1609 * === min(r0,r1) + dr >= d + min(r0,r1)
1610 * === dr >= d
1611 *
1612 * This, plus a fudge factor for FP operation is what we use below.
1613 *
1614 * An "Obtuse Cone" is defined to be one for which the "opening
1615 * angle" is obtuse.
1616 *
1617 * Consider two circles; one at (r0,r0) of radius r0, and one at
1618 * (r1,r1) of radius r1. These clearly lie on the acute/obtuse
1619 * boundary. The distance between the centres of these two circles
1620 * is d = sqr(2.(r0-r1)^2) by pythagoras. Thus d = sqr(2).dr.
1621 * By observation if d gets longer, we become acute, shorter, obtuse.
1622 * i.e. if sqr(2).dr > d we are obtuse, if d > sqr(2).dr we are acute.
1623 * (Thanks to Paul Gardiner for this reasoning).
1624 *
1625 * The code below tests (dr > d/3) (i.e. 3.dr > d). This
1626 * appears to be a factor of 2 and a bit out, so I am confused
1627 * by it.
1628 *
1629 * Either Igor meant something different to the standard meaning
1630 * of "Obtuse Cone", or he got his maths wrong. Or he was more
1631 * cunning than I can understand. Leave it as it until we find
1632 * an actual example that goes wrong.
1633 */
1634
1635 /* Tests with Acrobat seem to indicate that it uses a fudge factor
1636 * of around .0001. (i.e. [1.0001 0 0 0 0 1] is accepted as a
1637 * non nested circle, but [1.00009 0 0 0 0 1] is a nested one.
1638 * Approximate the same sort of value here to appease bug 690831.
1639 */
1640 if (any_abs (dr - d) < 0.001) {
1641 if ((r0 > r1 && Extend0) || (r1 > r0 && Extend1)) {
1642 r = R_rect_radius(rect, x0, y0);
1643 if (r0 < r1)
1644 code = R_tensor_annulus_extend_tangent(pfs, x0, y0, r0, t1, x1, y1, r1, t1, r);
1645 else
1646 code = R_tensor_annulus_extend_tangent(pfs, x1, y1, r1, t0, x0, y0, r0, t0, r);
1647 if (code < 0)
1648 return code;
1649 } else {
1650 if (r0 > r1) {
1651 if (Extend1 && r1 > 0)
1652 return R_tensor_annulus(pfs, x1, y1, r1, t1, x1, y1, 0, t1);
1653 }
1654 else {
1655 if (Extend0 && r0 > 0)
1656 return R_tensor_annulus(pfs, x0, y0, r0, t0, x0, y0, 0, t0);
1657 }
1658 }
1659 } else
1660 if (dr > d - 1e-4 * (d + dr)) {
1661 /* Nested circles, or degenerate. */
1662 if (r0 > r1) {
1663 if (Extend0) {
1664 r = R_rect_radius(rect, x0, y0);
1665 if (r > r0) {
1666 code = R_tensor_annulus(pfs, x0, y0, r, t0, x0, y0, r0, t0);
1667 if (code < 0)
1668 return code;
1669 }
1670 }
1671 if (Extend1 && r1 > 0)
1672 return R_tensor_annulus(pfs, x1, y1, r1, t1, x1, y1, 0, t1);
1673 } else {
1674 if (Extend1) {
1675 r = R_rect_radius(rect, x1, y1);
1676 if (r > r1) {
1677 code = R_tensor_annulus(pfs, x1, y1, r, t1, x1, y1, r1, t1);
1678 if (code < 0)
1679 return code;
1680 }
1681 }
1682 if (Extend0 && r0 > 0)
1683 return R_tensor_annulus(pfs, x0, y0, r0, t0, x0, y0, 0, t0);
1684 }
1685 } else if (dr > d / 3) {
1686 /* Obtuse cone. */
1687 if (r0 > r1) {
1688 if (Extend0) {
1689 r = R_rect_radius(rect, x0, y0);
1690 code = R_obtuse_cone(pfs, rect, x0, y0, r0, x1, y1, r1, t0, r);
1691 if (code < 0)
1692 return code;
1693 }
1694 if (Extend1 && r1 != 0)
1695 return R_tensor_cone_apex(pfs, rect, x0, y0, r0, x1, y1, r1, t1);
1696 return 0;
1697 } else {
1698 if (Extend1) {
1699 r = R_rect_radius(rect, x1, y1);
1700 code = R_obtuse_cone(pfs, rect, x1, y1, r1, x0, y0, r0, t1, r);
1701 if (code < 0)
1702 return code;
1703 }
1704 if (Extend0 && r0 != 0)
1705 return R_tensor_cone_apex(pfs, rect, x1, y1, r1, x0, y0, r0, t0);
1706 }
1707 } else {
1708 /* Acute cone or cylinder. */
1709 double x2, y2, r2, x3, y3, r3;
1710
1711 if (Extend0) {
1712 code = R_outer_circle(pfs, rect, x1, y1, r1, x0, y0, r0, &x3, &y3, &r3);
1713 if (code < 0)
1714 return code;
1715 if (x3 != x1 || y3 != y1) {
1716 code = R_tensor_annulus(pfs, x0, y0, r0, t0, x3, y3, r3, t0);
1717 if (code < 0)
1718 return code;
1719 }
1720 }
1721 if (Extend1) {
1722 code = R_outer_circle(pfs, rect, x0, y0, r0, x1, y1, r1, &x2, &y2, &r2);
1723 if (code < 0)
1724 return code;
1725 if (x2 != x0 || y2 != y0) {
1726 code = R_tensor_annulus(pfs, x1, y1, r1, t1, x2, y2, r2, t1);
1727 if (code < 0)
1728 return code;
1729 }
1730 }
1731 }
1732 return 0;
1733 }
1734
1735 static int
R_fill_rect_with_const_color(patch_fill_state_t * pfs,const gs_fixed_rect * clip_rect,float t)1736 R_fill_rect_with_const_color(patch_fill_state_t *pfs, const gs_fixed_rect *clip_rect, float t)
1737 {
1738 #if 0 /* Disabled because the clist writer device doesn't pass
1739 the clipping path with fill_recatangle. */
1740 patch_color_t pc;
1741 const gs_color_space *pcs = pfs->direct_space;
1742 gx_device_color dc;
1743 int code;
1744
1745 code = gs_function_evaluate(pfs->Function, &t, pc.cc.paint.values);
1746 if (code < 0)
1747 return code;
1748 pcs->type->restrict_color(&pc.cc, pcs);
1749 code = patch_color_to_device_color(pfs, &pc, &dc);
1750 if (code < 0)
1751 return code;
1752 return gx_fill_rectangle_device_rop(fixed2int_pixround(clip_rect->p.x), fixed2int_pixround(clip_rect->p.y),
1753 fixed2int_pixround(clip_rect->q.x) - fixed2int_pixround(clip_rect->p.x),
1754 fixed2int_pixround(clip_rect->q.y) - fixed2int_pixround(clip_rect->p.y),
1755 &dc, pfs->dev, pfs->pgs->log_op);
1756 #else
1757 /* Can't apply fill_rectangle, because the clist writer device doesn't pass
1758 the clipping path with fill_recatangle. Convert into trapezoids instead.
1759 */
1760 quadrangle_patch p;
1761 shading_vertex_t pp[2][2];
1762 const gs_color_space *pcs = pfs->direct_space;
1763 patch_color_t pc;
1764 int code;
1765
1766 code = gs_function_evaluate(pfs->Function, &t, pc.cc.paint.values);
1767 if (code < 0)
1768 return code;
1769 pcs->type->restrict_color(&pc.cc, pcs);
1770 pc.t[0] = pc.t[1] = t;
1771 pp[0][0].p = clip_rect->p;
1772 pp[0][1].p.x = clip_rect->q.x;
1773 pp[0][1].p.y = clip_rect->p.y;
1774 pp[1][0].p.x = clip_rect->p.x;
1775 pp[1][0].p.y = clip_rect->q.y;
1776 pp[1][1].p = clip_rect->q;
1777 pp[0][0].c = pp[0][1].c = pp[1][0].c = pp[1][1].c = &pc;
1778 p.p[0][0] = &pp[0][0];
1779 p.p[0][1] = &pp[0][1];
1780 p.p[1][0] = &pp[1][0];
1781 p.p[1][1] = &pp[1][1];
1782 return constant_color_quadrangle(pfs, &p, false);
1783 #endif
1784 }
1785
1786 typedef struct radial_shading_attrs_s {
1787 double x0, y0;
1788 double x1, y1;
1789 double span[2][2];
1790 double apex;
1791 bool have_apex;
1792 bool have_root[2]; /* ongoing contact, outgoing contact. */
1793 bool outer_contact[2];
1794 gs_point p[6]; /* 4 corners of the rectangle, p[4] = p[0], p[5] = p[1] */
1795 } radial_shading_attrs_t;
1796
1797 #define Pw2(a) ((a)*(a))
1798
1799 static void
radial_shading_external_contact(radial_shading_attrs_t * rsa,int point_index,double t,double r0,double r1,bool at_corner,int root_index)1800 radial_shading_external_contact(radial_shading_attrs_t *rsa, int point_index, double t, double r0, double r1, bool at_corner, int root_index)
1801 {
1802 double cx = rsa->x0 + (rsa->x1 - rsa->x0) * t;
1803 double cy = rsa->y0 + (rsa->y1 - rsa->y0) * t;
1804 double rx = rsa->p[point_index].x - cx;
1805 double ry = rsa->p[point_index].y - cy;
1806 double dx = rsa->p[point_index - 1].x - rsa->p[point_index].x;
1807 double dy = rsa->p[point_index - 1].y - rsa->p[point_index].y;
1808
1809 if (at_corner) {
1810 double Dx = rsa->p[point_index + 1].x - rsa->p[point_index].x;
1811 double Dy = rsa->p[point_index + 1].y - rsa->p[point_index].y;
1812 bool b1 = (dx * rx + dy * ry >= 0);
1813 bool b2 = (Dx * rx + Dy * ry >= 0);
1814
1815 if (b1 & b2)
1816 rsa->outer_contact[root_index] = true;
1817 } else {
1818 if (rx * dy - ry * dx < 0)
1819 rsa->outer_contact[root_index] = true;
1820 }
1821 }
1822
1823 static void
store_roots(radial_shading_attrs_t * rsa,const bool have_root[2],const double t[2],double r0,double r1,int point_index,bool at_corner)1824 store_roots(radial_shading_attrs_t *rsa, const bool have_root[2], const double t[2], double r0, double r1, int point_index, bool at_corner)
1825 {
1826 int i;
1827
1828 for (i = 0; i < 2; i++) {
1829 bool good_root;
1830
1831 if (!have_root[i])
1832 continue;
1833 good_root = (!rsa->have_apex || (rsa->apex <= 0 || r0 == 0 ? t[i] >= rsa->apex : t[i] <= rsa->apex));
1834 if (good_root) {
1835 radial_shading_external_contact(rsa, point_index, t[i], r0, r1, at_corner, i);
1836 if (!rsa->have_root[i]) {
1837 rsa->span[i][0] = rsa->span[i][1] = t[i];
1838 rsa->have_root[i] = true;
1839 } else {
1840 if (rsa->span[i][0] > t[i])
1841 rsa->span[i][0] = t[i];
1842 if (rsa->span[i][1] < t[i])
1843 rsa->span[i][1] = t[i];
1844 }
1845 }
1846 }
1847 }
1848
1849 static void
compute_radial_shading_span_extended_side(radial_shading_attrs_t * rsa,double r0,double r1,int point_index)1850 compute_radial_shading_span_extended_side(radial_shading_attrs_t *rsa, double r0, double r1, int point_index)
1851 {
1852 double cc, c;
1853 bool have_root[2] = {false, false};
1854 double t[2];
1855 bool by_x = (rsa->p[point_index].x != rsa->p[point_index + 1].x);
1856 int i;
1857
1858 /* As t moves from 0 to 1, the circles move from r0 to r1, and from
1859 * from position p0 to py. For simplicity, adjust so that p0 is at
1860 * the origin. Consider the projection of the circle drawn at any given
1861 * time onto the x axis. The range of points would be:
1862 * p1x*t +/- (r0+(r1-r0)*t). We are interested in the first (and last)
1863 * moments when the range includes a point c on the x axis. So solve for:
1864 * p1x*t +/- (r0+(r1-r0)*t) = c. Let cc = p1x.
1865 * So p1x*t0 + (r1-r0)*t0 = c - r0 => t0 = (c - r0)/(p1x + r1 - r0)
1866 * p1x*t1 - (r1-r0)*t1 = c + r0 => t1 = (c + r0)/(p1x - r1 + r0)
1867 */
1868 if (by_x) {
1869 c = rsa->p[point_index].x - rsa->x0;
1870 cc = rsa->x1 - rsa->x0;
1871 } else {
1872 c = rsa->p[point_index].y - rsa->y0;
1873 cc = rsa->y1 - rsa->y0;
1874 }
1875 t[0] = (c - r0) / (cc + r1 - r0);
1876 t[1] = (c + r0) / (cc - r1 + r0);
1877 if (t[0] > t[1]) {
1878 c = t[0];
1879 t[0] = t[1];
1880 t[1] = c;
1881 }
1882 for (i = 0; i < 2; i++) {
1883 double d, d0, d1;
1884
1885 if (by_x) {
1886 d = rsa->y1 - rsa->y0 + r0 + (r1 - r0) * t[i];
1887 d0 = rsa->p[point_index].y;
1888 d1 = rsa->p[point_index + 1].y;
1889 } else {
1890 d = rsa->x1 - rsa->x0 + r0 + (r1 - r0) * t[i];
1891 d0 = rsa->p[point_index].x;
1892 d1 = rsa->p[point_index + 1].x;
1893 }
1894 if (d1 > d0 ? d0 <= d && d <= d1 : d1 <= d && d <= d0)
1895 have_root[i] = true;
1896 }
1897 store_roots(rsa, have_root, t, r0, r1, point_index, false);
1898 }
1899
1900 static int
compute_radial_shading_span_extended_point(radial_shading_attrs_t * rsa,double r0,double r1,int point_index)1901 compute_radial_shading_span_extended_point(radial_shading_attrs_t *rsa, double r0, double r1, int point_index)
1902 {
1903 /* As t moves from 0 to 1, the circles move from r0 to r1, and from
1904 * from position p0 to py. At any given time t, therefore, we
1905 * paint the points that are distance r0+(r1-r0)*t from point
1906 * (p0x+(p1x-p0x)*t,p0y+(p1y-p0y)*t) = P(t).
1907 *
1908 * To simplify our algebra, adjust so that (p0x, p0y) is at the origin.
1909 * To find the time(s) t at which the a point q is painted, we therefore
1910 * solve for t in:
1911 *
1912 * |q-P(t)| = r0+(r1-r0)*t
1913 *
1914 * (qx-p1x*t)^2 + (qy-p1y*t)^2 - (r0+(r1-r0)*t)^2 = 0
1915 * = qx^2 - 2qx.p1x.t + p1x^2.t^2 + qy^2 - 2qy.p1y.t + p1y^2.t^2 -
1916 * (r0^2 + 2r0(r1-r0)t + (r1-r0)^2.t^2)
1917 * = qx^2 + qy^2 - r0^2
1918 * + -2(qx.p1x + qy.p1y + r0(r1-r0)).t
1919 * + (p1x^2 + p1y^2 - (r1-r0)^2).t^2
1920 *
1921 * So solve using the usual t = (-b +/- SQRT(b^2 - 4ac)) where
1922 * a = p1x^2 + p1y^2 - (r1-r0)^2
1923 * b = -2(qx.p1x + qy.p1y + r0(r1-r0))
1924 * c = qx^2 + qy^2 - r0^2
1925 */
1926 double p1x = rsa->x1 - rsa->x0;
1927 double p1y = rsa->y1 - rsa->y0;
1928 double qx = rsa->p[point_index].x - rsa->x0;
1929 double qy = rsa->p[point_index].y - rsa->y0;
1930 double a = (Pw2(p1x) + Pw2(p1y) - Pw2(r0 - r1));
1931 bool have_root[2] = {false, false};
1932 double t[2];
1933
1934 if (fabs(a) < 1e-8) {
1935 /* Linear equation. */
1936 /* This case is always the ongoing ellipse contact. */
1937 double cx = rsa->x0 - (rsa->x1 - rsa->x0) * r0 / (r1 - r0);
1938 double cy = rsa->y0 - (rsa->y1 - rsa->y0) * r0 / (r1 - r0);
1939
1940 t[0] = (Pw2(qx) + Pw2(qy))/(cx*qx + cy*qy) / 2;
1941 have_root[0] = true;
1942 } else {
1943 /* Square equation. No solution if b^2 - 4ac = 0. Equivalently if
1944 * (b^2)/4 -a.c = 0 === (b/2)^2 - a.c = 0 === (-b/2)^2 - a.c = 0 */
1945 double minushalfb = r0*(r1-r0) + p1x*qx + p1y*qy;
1946 double c = Pw2(qx) + Pw2(qy) - Pw2(r0);
1947 double desc2 = Pw2(minushalfb) - a*c; /* desc2 = 1/4 (b^2-4ac) */
1948
1949 if (desc2 < 0) {
1950 return -1; /* The point is outside the shading coverage.
1951 Do not shorten, because we didn't observe it in practice. */
1952 } else {
1953 double desc1 = sqrt(desc2); /* desc1 = 1/2 SQRT(b^2-4ac) */
1954
1955 if (a > 0) {
1956 t[0] = (minushalfb - desc1) / a;
1957 t[1] = (minushalfb + desc1) / a;
1958 } else {
1959 t[0] = (minushalfb + desc1) / a;
1960 t[1] = (minushalfb - desc1) / a;
1961 }
1962 have_root[0] = have_root[1] = true;
1963 }
1964 }
1965 store_roots(rsa, have_root, t, r0, r1, point_index, true);
1966 if (have_root[0] && have_root[1])
1967 return 15;
1968 if (have_root[0])
1969 return 15 - 4;
1970 if (have_root[1])
1971 return 15 - 2;
1972 return -1;
1973 }
1974
1975 #undef Pw2
1976
1977 static int
compute_radial_shading_span_extended(radial_shading_attrs_t * rsa,double r0,double r1)1978 compute_radial_shading_span_extended(radial_shading_attrs_t *rsa, double r0, double r1)
1979 {
1980 int span_type0, span_type1;
1981
1982 span_type0 = compute_radial_shading_span_extended_point(rsa, r0, r1, 1);
1983 if (span_type0 == -1)
1984 return -1;
1985 span_type1 = compute_radial_shading_span_extended_point(rsa, r0, r1, 2);
1986 if (span_type0 != span_type1)
1987 return -1;
1988 span_type1 = compute_radial_shading_span_extended_point(rsa, r0, r1, 3);
1989 if (span_type0 != span_type1)
1990 return -1;
1991 span_type1 = compute_radial_shading_span_extended_point(rsa, r0, r1, 4);
1992 if (span_type0 != span_type1)
1993 return -1;
1994 compute_radial_shading_span_extended_side(rsa, r0, r1, 1);
1995 compute_radial_shading_span_extended_side(rsa, r0, r1, 2);
1996 compute_radial_shading_span_extended_side(rsa, r0, r1, 3);
1997 compute_radial_shading_span_extended_side(rsa, r0, r1, 4);
1998 return span_type0;
1999 }
2000
2001 static int
compute_radial_shading_span(radial_shading_attrs_t * rsa,float x0,float y0,double r0,float x1,float y1,double r1,const gs_rect * rect)2002 compute_radial_shading_span(radial_shading_attrs_t *rsa, float x0, float y0, double r0, float x1, float y1, double r1, const gs_rect * rect)
2003 {
2004 /* If the shading area is much larger than the path bbox,
2005 we want to shorten the shading for a faster rendering.
2006 If any point of the path bbox falls outside the shading area,
2007 our math is not applicable, and we render entire shading.
2008 If the path bbox is inside the shading area,
2009 we compute 1 or 2 'spans' - the shading parameter intervals,
2010 which covers the bbox. For doing that we need to resolve
2011 a square eqation by the shading parameter
2012 for each corner of the bounding box,
2013 and for each side of the shading bbox.
2014 Note the equation to be solved in the user space.
2015 Since each equation gives 2 roots (because the points are
2016 strongly inside the shading area), we will get 2 parameter intervals -
2017 the 'lower' one corresponds to the first (ongoing) contact of
2018 the running circle, and the second one corresponds to the last (outgoing) contact
2019 (like in a sun eclipse; well our sun is rectangular).
2020
2021 Here are few exceptions.
2022
2023 First, the equation degenerates when the distance sqrt((x1-x0)^2 + (y1-y0)^2)
2024 appears equal to r0-r1. In this case the base circles do contact,
2025 and the running circle does contact at the same point.
2026 The equation degenerates to a linear one.
2027 Since we don't want float precision noize to affect the result,
2028 we compute this condition in 'fixed' coordinates.
2029
2030 Second, Postscript approximates any circle with 3d order beziers.
2031 This approximation may give a 2% error.
2032 Therefore using the precise roots may cause a dropout.
2033 To prevetn them, we slightly modify the base radii.
2034 However the sign of modification smartly depends
2035 on the relative sizes of the base circles,
2036 and on the contact number. Currently we don't want to
2037 define and debug the smart optimal logic for that,
2038 so we simply try all 4 variants for each source equation,
2039 and use the union of intervals.
2040
2041 Third, we could compute which quarter of the circle
2042 really covers the path bbox. Using it we could skip
2043 rendering of uncovering quarters. Currently we do not
2044 implement this optimization. The general tensor patch algorithm
2045 will skip uncovering parts.
2046
2047 Fourth, when one base circle is (almost) inside the other,
2048 the parameter interval must include the shading apex.
2049 To know that, we determine whether the contacting circle
2050 is outside the rectangle (the "outer" contact),
2051 or it is (partially) inside the rectangle.
2052
2053 At last, a small shortening of a shading won't give a
2054 sensible speedup, but it may replace a symmetric function domain
2055 with an assymmetric one, so that the rendering
2056 would be asymmetyric for a symmetric shading.
2057 Therefore we do not perform a small sortening.
2058 Instead we shorten only if the shading span
2059 is much smaller that the shading domain.
2060 */
2061 const double extent = 1.02;
2062 int span_type0, span_type1, span_type;
2063
2064 memset(rsa, 0, sizeof(*rsa));
2065 rsa->x0 = x0;
2066 rsa->y0 = y0;
2067 rsa->x1 = x1;
2068 rsa->y1 = y1;
2069 rsa->p[0] = rsa->p[4] = rect->p;
2070 rsa->p[1].x = rsa->p[5].x = rect->p.x;
2071 rsa->p[1].y = rsa->p[5].y = rect->q.y;
2072 rsa->p[2] = rect->q;
2073 rsa->p[3].x = rect->q.x;
2074 rsa->p[3].y = rect->p.y;
2075 rsa->have_apex = any_abs(r1 - r0) > 1e-7 * any_abs(r1 + r0);
2076 rsa->apex = (rsa->have_apex ? -r0 / (r1 - r0) : 0);
2077 span_type0 = compute_radial_shading_span_extended(rsa, r0 / extent, r1 * extent);
2078 if (span_type0 == -1)
2079 return -1;
2080 span_type1 = compute_radial_shading_span_extended(rsa, r0 / extent, r1 / extent);
2081 if (span_type0 != span_type1)
2082 return -1;
2083 span_type1 = compute_radial_shading_span_extended(rsa, r0 * extent, r1 * extent);
2084 if (span_type0 != span_type1)
2085 return -1;
2086 span_type1 = compute_radial_shading_span_extended(rsa, r0 * extent, r1 / extent);
2087 if (span_type1 == -1)
2088 return -1;
2089 if (r0 < r1) {
2090 if (rsa->have_root[0] && !rsa->outer_contact[0])
2091 rsa->span[0][0] = rsa->apex; /* Likely never happens. Remove ? */
2092 if (rsa->have_root[1] && !rsa->outer_contact[1])
2093 rsa->span[1][0] = rsa->apex;
2094 } else if (r0 > r1) {
2095 if (rsa->have_root[0] && !rsa->outer_contact[0])
2096 rsa->span[0][1] = rsa->apex;
2097 if (rsa->have_root[1] && !rsa->outer_contact[1])
2098 rsa->span[1][1] = rsa->apex; /* Likely never happens. Remove ? */
2099 }
2100 span_type = 0;
2101 if (rsa->have_root[0] && rsa->span[0][0] < 0)
2102 span_type |= 1;
2103 if (rsa->have_root[1] && rsa->span[1][0] < 0)
2104 span_type |= 1;
2105 if (rsa->have_root[0] && rsa->span[0][1] > 0 && rsa->span[0][0] < 1)
2106 span_type |= 2;
2107 if (rsa->have_root[1] && rsa->span[1][1] > 0 && rsa->span[1][0] < 1)
2108 span_type |= 4;
2109 if (rsa->have_root[0] && rsa->span[0][1] > 1)
2110 span_type |= 8;
2111 if (rsa->have_root[1] && rsa->span[1][1] > 1)
2112 span_type |= 8;
2113 return span_type;
2114 }
2115
2116 static bool
shorten_radial_shading(float * x0,float * y0,double * r0,float * d0,float * x1,float * y1,double * r1,float * d1,double span_[2])2117 shorten_radial_shading(float *x0, float *y0, double *r0, float *d0, float *x1, float *y1, double *r1, float *d1, double span_[2])
2118 {
2119 double s0 = span_[0], s1 = span_[1], w;
2120
2121 if (s0 < 0)
2122 s0 = 0;
2123 if (s1 < 0)
2124 s1 = 0;
2125 if (s0 > 1)
2126 s0 = 1;
2127 if (s1 > 1)
2128 s1 = 1;
2129 w = s1 - s0;
2130 if (w == 0)
2131 return false; /* Don't pass a degenerate shading. */
2132 if (w > 0.3)
2133 return false; /* The span is big, don't shorten it. */
2134 { /* Do shorten. */
2135 double R0 = *r0, X0 = *x0, Y0 = *y0, D0 = *d0;
2136 double R1 = *r1, X1 = *x1, Y1 = *y1, D1 = *d1;
2137
2138 *r0 = R0 + (R1 - R0) * s0;
2139 *x0 = X0 + (X1 - X0) * s0;
2140 *y0 = Y0 + (Y1 - Y0) * s0;
2141 *d0 = D0 + (D1 - D0) * s0;
2142 *r1 = R0 + (R1 - R0) * s1;
2143 *x1 = X0 + (X1 - X0) * s1;
2144 *y1 = Y0 + (Y1 - Y0) * s1;
2145 *d1 = D0 + (D1 - D0) * s1;
2146 }
2147 return true;
2148 }
2149
2150 static bool inline
is_radial_shading_large(double x0,double y0,double r0,double x1,double y1,double r1,const gs_rect * rect)2151 is_radial_shading_large(double x0, double y0, double r0, double x1, double y1, double r1, const gs_rect * rect)
2152 {
2153 const double d = hypot(x1 - x0, y1 - y0);
2154 const double area0 = M_PI * r0 * r0 / 2;
2155 const double area1 = M_PI * r1 * r1 / 2;
2156 const double area2 = (r0 + r1) / 2 * d;
2157 const double arbitrary = 8;
2158 double areaX, areaY;
2159
2160 /* The shading area is not equal to area0 + area1 + area2
2161 when one circle is (almost) inside the other.
2162 We believe that the 'arbitrary' coefficient recovers that
2163 when it is set greater than 2. */
2164 /* If one dimension is large enough, the shading parameter span is wide. */
2165 areaX = (rect->q.x - rect->p.x) * (rect->q.x - rect->p.x);
2166 if (areaX * arbitrary < area0 + area1 + area2)
2167 return true;
2168 areaY = (rect->q.y - rect->p.y) * (rect->q.y - rect->p.y);
2169 if (areaY * arbitrary < area0 + area1 + area2)
2170 return true;
2171 return false;
2172 }
2173
2174 static int
gs_shading_R_fill_rectangle_aux(const gs_shading_t * psh0,const gs_rect * rect,const gs_fixed_rect * clip_rect,gx_device * dev,gs_gstate * pgs)2175 gs_shading_R_fill_rectangle_aux(const gs_shading_t * psh0, const gs_rect * rect,
2176 const gs_fixed_rect *clip_rect,
2177 gx_device * dev, gs_gstate * pgs)
2178 {
2179 const gs_shading_R_t *const psh = (const gs_shading_R_t *)psh0;
2180 float d0 = psh->params.Domain[0], d1 = psh->params.Domain[1];
2181 float x0 = psh->params.Coords[0], y0 = psh->params.Coords[1];
2182 double r0 = psh->params.Coords[2];
2183 float x1 = psh->params.Coords[3], y1 = psh->params.Coords[4];
2184 double r1 = psh->params.Coords[5];
2185 radial_shading_attrs_t rsa;
2186 int span_type; /* <0 - don't shorten, 1 - extent0, 2 - first contact, 4 - last contact, 8 - extent1. */
2187 int code;
2188 patch_fill_state_t pfs1;
2189
2190 if (r0 == 0 && r1 == 0)
2191 return 0; /* PLRM requires to paint nothing. */
2192 code = shade_init_fill_state((shading_fill_state_t *)&pfs1, psh0, dev, pgs);
2193 if (code < 0)
2194 return code;
2195 pfs1.Function = psh->params.Function;
2196 code = init_patch_fill_state(&pfs1);
2197 if (code < 0) {
2198 if (pfs1.icclink != NULL) gsicc_release_link(pfs1.icclink);
2199 return code;
2200 }
2201 pfs1.function_arg_shift = 0;
2202 pfs1.rect = *clip_rect;
2203 pfs1.maybe_self_intersecting = false;
2204 if (is_radial_shading_large(x0, y0, r0, x1, y1, r1, rect))
2205 span_type = compute_radial_shading_span(&rsa, x0, y0, r0, x1, y1, r1, rect);
2206 else
2207 span_type = -1;
2208 if (span_type < 0) {
2209 code = R_extensions(&pfs1, psh, rect, d0, d1, psh->params.Extend[0], false);
2210 if (code >= 0)
2211 code = R_tensor_annulus(&pfs1, x0, y0, r0, d0, x1, y1, r1, d1);
2212 if (code >= 0)
2213 code = R_extensions(&pfs1, psh, rect, d0, d1, false, psh->params.Extend[1]);
2214 } else if (span_type == 1) {
2215 code = R_fill_rect_with_const_color(&pfs1, clip_rect, d0);
2216 } else if (span_type == 8) {
2217 code = R_fill_rect_with_const_color(&pfs1, clip_rect, d1);
2218 } else {
2219 bool second_interval = true;
2220
2221 code = 0;
2222 if (span_type & 1)
2223 code = R_extensions(&pfs1, psh, rect, d0, d1, psh->params.Extend[0], false);
2224 if ((code >= 0) && (span_type & 2)) {
2225 float X0 = x0, Y0 = y0, D0 = d0, X1 = x1, Y1 = y1, D1 = d1;
2226 double R0 = r0, R1 = r1;
2227
2228 if ((span_type & 4) && rsa.span[0][1] >= rsa.span[1][0]) {
2229 double united[2];
2230
2231 united[0] = rsa.span[0][0];
2232 united[1] = rsa.span[1][1];
2233 shorten_radial_shading(&X0, &Y0, &R0, &D0, &X1, &Y1, &R1, &D1, united);
2234 second_interval = false;
2235 } else {
2236 second_interval = shorten_radial_shading(&X0, &Y0, &R0, &D0, &X1, &Y1, &R1, &D1, rsa.span[0]);
2237 }
2238 code = R_tensor_annulus(&pfs1, X0, Y0, R0, D0, X1, Y1, R1, D1);
2239 }
2240 if (code >= 0 && second_interval) {
2241 if (span_type & 4) {
2242 float X0 = x0, Y0 = y0, D0 = d0, X1 = x1, Y1 = y1, D1 = d1;
2243 double R0 = r0, R1 = r1;
2244
2245 shorten_radial_shading(&X0, &Y0, &R0, &D0, &X1, &Y1, &R1, &D1, rsa.span[1]);
2246 code = R_tensor_annulus(&pfs1, X0, Y0, R0, D0, X1, Y1, R1, D1);
2247 }
2248 }
2249 if (code >= 0 && (span_type & 8))
2250 code = R_extensions(&pfs1, psh, rect, d0, d1, false, psh->params.Extend[1]);
2251 }
2252 if (pfs1.icclink != NULL) gsicc_release_link(pfs1.icclink);
2253 if (term_patch_fill_state(&pfs1))
2254 return_error(gs_error_unregistered); /* Must not happen. */
2255 return code;
2256 }
2257
2258 int
gs_shading_R_fill_rectangle(const gs_shading_t * psh0,const gs_rect * rect,const gs_fixed_rect * rect_clip,gx_device * dev,gs_gstate * pgs)2259 gs_shading_R_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,
2260 const gs_fixed_rect * rect_clip,
2261 gx_device * dev, gs_gstate * pgs)
2262 {
2263 return gs_shading_R_fill_rectangle_aux(psh0, rect, rect_clip, dev, pgs);
2264 }
2265