1/* $id: shapes.c,v 1.82 2007/12/24 04:50:36 ellson Exp $ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2012 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14/* This code is derived from the Java implementation by Luc Maisonobe */
15/* Copyright (c) 2003-2004, Luc Maisonobe
16 * All rights reserved.
17 *
18 * Redistribution and use in source and binary forms, with
19 * or without modification, are permitted provided that
20 * the following conditions are met:
21 *
22 * Redistributions of source code must retain the
23 * above copyright notice, this list of conditions and
24 * the following disclaimer.
25 * Redistributions in binary form must reproduce the
26 * above copyright notice, this list of conditions and
27 * the following disclaimer in the documentation
28 * and/or other materials provided with the
29 * distribution.
30 * Neither the names of spaceroots.org, spaceroots.com
31 * nor the names of their contributors may be used to
32 * endorse or promote products derived from this
33 * software without specific prior written permission.
34 *
35 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
36 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
37 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
38 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
39 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
40 * THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
41 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
43 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
44 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
45 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
46 * IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
47 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
48 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 * POSSIBILITY OF SUCH DAMAGE.
50 */
51
52#if STANDALONE
53#include <limits.h>
54#include <math.h>
55#include <stdlib.h>
56#include <stdio.h>
57
58#define MAX(a,b) ((a)>(b)?(a):(b))
59#define MIN(a,b) ((a)<(b)?(a):(b))
60
61#define NEW(t) ((t*)calloc(1,sizeof(t)))
62#define N_NEW(n,t) ((t*)calloc(n,sizeof(t)))
63
64#define PI 3.14159265358979323846
65
66#define TRUE 1
67#define FALSE 0
68typedef unsigned char boolean;
69
70typedef struct pointf_s {
71 double x, y;
72} pointf;
73typedef struct Ppoly_t {
74 pointf *ps;
75 int pn;
76} Ppoly_t;
77
78typedef Ppoly_t Ppolyline_t;
79#else
80#include "render.h"
81#include "pathplan.h"
82#endif
83
84#define TWOPI (2*M_PI)
85
86typedef struct {
87 double cx, cy; /* center */
88 double a, b; /* semi-major and -minor axes */
89
90 /* Orientation of the major axis with respect to the x axis. */
91 double theta, cosTheta, sinTheta;
92
93 /* Start and end angles of the arc. */
94 double eta1, eta2;
95
96 /* Position of the start and end points. */
97 double x1, y1, x2, y2;
98
99 /* Position of the foci. */
100 double xF1, yF1, xF2, yF2;
101
102 /* x of the leftmost point of the arc. */
103 double xLeft;
104
105 /* y of the highest point of the arc. */
106 double yUp;
107
108 /* Horizontal width and vertical height of the arc. */
109 double width, height;
110
111 double f, e2, g, g2;
112} ellipse_t;
113
114static void computeFoci(ellipse_t * ep)
115{
116 double d = sqrt(ep->a * ep->a - ep->b * ep->b);
117 double dx = d * ep->cosTheta;
118 double dy = d * ep->sinTheta;
119
120 ep->xF1 = ep->cx - dx;
121 ep->yF1 = ep->cy - dy;
122 ep->xF2 = ep->cx + dx;
123 ep->yF2 = ep->cy + dy;
124}
125
126 /* Compute the locations of the endpoints. */
127static void computeEndPoints(ellipse_t * ep)
128{
129 double aCosEta1 = ep->a * cos(ep->eta1);
130 double bSinEta1 = ep->b * sin(ep->eta1);
131 double aCosEta2 = ep->a * cos(ep->eta2);
132 double bSinEta2 = ep->b * sin(ep->eta2);
133
134 // start point
135 ep->x1 = ep->cx + aCosEta1 * ep->cosTheta - bSinEta1 * ep->sinTheta;
136 ep->y1 = ep->cy + aCosEta1 * ep->sinTheta + bSinEta1 * ep->cosTheta;
137
138 // end point
139 ep->x2 = ep->cx + aCosEta2 * ep->cosTheta - bSinEta2 * ep->sinTheta;
140 ep->y2 = ep->cy + aCosEta2 * ep->sinTheta + bSinEta2 * ep->cosTheta;
141}
142
143 /* Compute the bounding box. */
144static void computeBounds(ellipse_t * ep)
145{
146 double bOnA = ep->b / ep->a;
147 double etaXMin, etaXMax, etaYMin, etaYMax;
148
149 if (fabs(ep->sinTheta) < 0.1) {
150 double tanTheta = ep->sinTheta / ep->cosTheta;
151 if (ep->cosTheta < 0) {
152 etaXMin = -atan(tanTheta * bOnA);
153 etaXMax = etaXMin + M_PI;
154 etaYMin = 0.5 * M_PI - atan(tanTheta / bOnA);
155 etaYMax = etaYMin + M_PI;
156 } else {
157 etaXMax = -atan(tanTheta * bOnA);
158 etaXMin = etaXMax - M_PI;
159 etaYMax = 0.5 * M_PI - atan(tanTheta / bOnA);
160 etaYMin = etaYMax - M_PI;
161 }
162 } else {
163 double invTanTheta = ep->cosTheta / ep->sinTheta;
164 if (ep->sinTheta < 0) {
165 etaXMax = 0.5 * M_PI + atan(invTanTheta / bOnA);
166 etaXMin = etaXMax - M_PI;
167 etaYMin = atan(invTanTheta * bOnA);
168 etaYMax = etaYMin + M_PI;
169 } else {
170 etaXMin = 0.5 * M_PI + atan(invTanTheta / bOnA);
171 etaXMax = etaXMin + M_PI;
172 etaYMax = atan(invTanTheta * bOnA);
173 etaYMin = etaYMax - M_PI;
174 }
175 }
176
177 etaXMin -= (TWOPI * floor((etaXMin - ep->eta1) / TWOPI));
178 etaYMin -= (TWOPI * floor((etaYMin - ep->eta1) / TWOPI));
179 etaXMax -= (TWOPI * floor((etaXMax - ep->eta1) / TWOPI));
180 etaYMax -= (TWOPI * floor((etaYMax - ep->eta1) / TWOPI));
181
182 ep->xLeft = (etaXMin <= ep->eta2)
183 ? (ep->cx + ep->a * cos(etaXMin) * ep->cosTheta -
184 ep->b * sin(etaXMin) * ep->sinTheta)
185 : MIN(ep->x1, ep->x2);
186 ep->yUp = (etaYMin <= ep->eta2)
187 ? (ep->cy + ep->a * cos(etaYMin) * ep->sinTheta +
188 ep->b * sin(etaYMin) * ep->cosTheta)
189 : MIN(ep->y1, ep->y2);
190 ep->width = ((etaXMax <= ep->eta2)
191 ? (ep->cx + ep->a * cos(etaXMax) * ep->cosTheta -
192 ep->b * sin(etaXMax) * ep->sinTheta)
193 : MAX(ep->x1, ep->x2)) - ep->xLeft;
194 ep->height = ((etaYMax <= ep->eta2)
195 ? (ep->cy + ep->a * cos(etaYMax) * ep->sinTheta +
196 ep->b * sin(etaYMax) * ep->cosTheta)
197 : MAX(ep->y1, ep->y2)) - ep->yUp;
198
199}
200
201static void
202initEllipse(ellipse_t * ep, double cx, double cy, double a, double b,
203 double theta, double lambda1, double lambda2)
204{
205 ep->cx = cx;
206 ep->cy = cy;
207 ep->a = a;
208 ep->b = b;
209 ep->theta = theta;
210
211 ep->eta1 = atan2(sin(lambda1) / b, cos(lambda1) / a);
212 ep->eta2 = atan2(sin(lambda2) / b, cos(lambda2) / a);
213 ep->cosTheta = cos(theta);
214 ep->sinTheta = sin(theta);
215
216 // make sure we have eta1 <= eta2 <= eta1 + 2*PI
217 ep->eta2 -= TWOPI * floor((ep->eta2 - ep->eta1) / TWOPI);
218
219 // the preceding correction fails if we have exactly eta2 - eta1 = 2*PI
220 // it reduces the interval to zero length
221 if ((lambda2 - lambda1 > M_PI) && (ep->eta2 - ep->eta1 < M_PI)) {
222 ep->eta2 += TWOPI;
223 }
224
225 computeFoci(ep);
226 computeEndPoints(ep);
227 computeBounds(ep);
228
229 /* Flatness parameters */
230 ep->f = (ep->a - ep->b) / ep->a;
231 ep->e2 = ep->f * (2.0 - ep->f);
232 ep->g = 1.0 - ep->f;
233 ep->g2 = ep->g * ep->g;
234}
235
236typedef double erray_t[2][4][4];
237
238 // coefficients for error estimation
239 // while using quadratic Bezier curves for approximation
240 // 0 < b/a < 1/4
241static erray_t coeffs2Low = {
242 {
243 {3.92478, -13.5822, -0.233377, 0.0128206},
244 {-1.08814, 0.859987, 0.000362265, 0.000229036},
245 {-0.942512, 0.390456, 0.0080909, 0.00723895},
246 {-0.736228, 0.20998, 0.0129867, 0.0103456}
247 },
248 {
249 {-0.395018, 6.82464, 0.0995293, 0.0122198},
250 {-0.545608, 0.0774863, 0.0267327, 0.0132482},
251 {0.0534754, -0.0884167, 0.012595, 0.0343396},
252 {0.209052, -0.0599987, -0.00723897, 0.00789976}
253 }
254};
255
256 // coefficients for error estimation
257 // while using quadratic Bezier curves for approximation
258 // 1/4 <= b/a <= 1
259static erray_t coeffs2High = {
260 {
261 {0.0863805, -11.5595, -2.68765, 0.181224},
262 {0.242856, -1.81073, 1.56876, 1.68544},
263 {0.233337, -0.455621, 0.222856, 0.403469},
264 {0.0612978, -0.104879, 0.0446799, 0.00867312}
265 },
266 {
267 {0.028973, 6.68407, 0.171472, 0.0211706},
268 {0.0307674, -0.0517815, 0.0216803, -0.0749348},
269 {-0.0471179, 0.1288, -0.0781702, 2.0},
270 {-0.0309683, 0.0531557, -0.0227191, 0.0434511}
271 }
272};
273
274 // safety factor to convert the "best" error approximation
275 // into a "max bound" error
276static double safety2[] = {
277 0.02, 2.83, 0.125, 0.01
278};
279
280 // coefficients for error estimation
281 // while using cubic Bezier curves for approximation
282 // 0 < b/a < 1/4
283static erray_t coeffs3Low = {
284 {
285 {3.85268, -21.229, -0.330434, 0.0127842},
286 {-1.61486, 0.706564, 0.225945, 0.263682},
287 {-0.910164, 0.388383, 0.00551445, 0.00671814},
288 {-0.630184, 0.192402, 0.0098871, 0.0102527}
289 },
290 {
291 {-0.162211, 9.94329, 0.13723, 0.0124084},
292 {-0.253135, 0.00187735, 0.0230286, 0.01264},
293 {-0.0695069, -0.0437594, 0.0120636, 0.0163087},
294 {-0.0328856, -0.00926032, -0.00173573, 0.00527385}
295 }
296};
297
298 // coefficients for error estimation
299 // while using cubic Bezier curves for approximation
300 // 1/4 <= b/a <= 1
301static erray_t coeffs3High = {
302 {
303 {0.0899116, -19.2349, -4.11711, 0.183362},
304 {0.138148, -1.45804, 1.32044, 1.38474},
305 {0.230903, -0.450262, 0.219963, 0.414038},
306 {0.0590565, -0.101062, 0.0430592, 0.0204699}
307 },
308 {
309 {0.0164649, 9.89394, 0.0919496, 0.00760802},
310 {0.0191603, -0.0322058, 0.0134667, -0.0825018},
311 {0.0156192, -0.017535, 0.00326508, -0.228157},
312 {-0.0236752, 0.0405821, -0.0173086, 0.176187}
313 }
314};
315
316 // safety factor to convert the "best" error approximation
317 // into a "max bound" error
318static double safety3[] = {
319 0.001, 4.98, 0.207, 0.0067
320};
321
322/* Compute the value of a rational function.
323 * This method handles rational functions where the numerator is
324 * quadratic and the denominator is linear
325 */
326#define RationalFunction(x,c) ((x * (x * c[0] + c[1]) + c[2]) / (x + c[3]))
327
328/* Estimate the approximation error for a sub-arc of the instance.
329 * degree specifies degree of the Bezier curve to use (1, 2 or 3)
330 * tA and tB give the start and end angle of the subarc
331 * Returns upper bound of the approximation error between the Bezier
332 * curve and the real ellipse
333 */
334static double
335estimateError(ellipse_t * ep, int degree, double etaA, double etaB)
336{
337 double c0, c1, eta = 0.5 * (etaA + etaB);
338
339 if (degree < 2) {
340
341 // start point
342 double aCosEtaA = ep->a * cos(etaA);
343 double bSinEtaA = ep->b * sin(etaA);
344 double xA =
345 ep->cx + aCosEtaA * ep->cosTheta - bSinEtaA * ep->sinTheta;
346 double yA =
347 ep->cy + aCosEtaA * ep->sinTheta + bSinEtaA * ep->cosTheta;
348
349 // end point
350 double aCosEtaB = ep->a * cos(etaB);
351 double bSinEtaB = ep->b * sin(etaB);
352 double xB =
353 ep->cx + aCosEtaB * ep->cosTheta - bSinEtaB * ep->sinTheta;
354 double yB =
355 ep->cy + aCosEtaB * ep->sinTheta + bSinEtaB * ep->cosTheta;
356
357 // maximal error point
358 double aCosEta = ep->a * cos(eta);
359 double bSinEta = ep->b * sin(eta);
360 double x =
361 ep->cx + aCosEta * ep->cosTheta - bSinEta * ep->sinTheta;
362 double y =
363 ep->cy + aCosEta * ep->sinTheta + bSinEta * ep->cosTheta;
364
365 double dx = xB - xA;
366 double dy = yB - yA;
367
368 return fabs(x * dy - y * dx + xB * yA - xA * yB)
369 / sqrt(dx * dx + dy * dy);
370
371 } else {
372
373 double x = ep->b / ep->a;
374 double dEta = etaB - etaA;
375 double cos2 = cos(2 * eta);
376 double cos4 = cos(4 * eta);
377 double cos6 = cos(6 * eta);
378
379 // select the right coefficient's set according to degree and b/a
380 double (*coeffs)[4][4];
381 double *safety;
382 if (degree == 2) {
383 coeffs = (x < 0.25) ? coeffs2Low : coeffs2High;
384 safety = safety2;
385 } else {
386 coeffs = (x < 0.25) ? coeffs3Low : coeffs3High;
387 safety = safety3;
388 }
389
390 c0 = RationalFunction(x, coeffs[0][0])
391 + cos2 * RationalFunction(x, coeffs[0][1])
392 + cos4 * RationalFunction(x, coeffs[0][2])
393 + cos6 * RationalFunction(x, coeffs[0][3]);
394
395 c1 = RationalFunction(x, coeffs[1][0])
396 + cos2 * RationalFunction(x, coeffs[1][1])
397 + cos4 * RationalFunction(x, coeffs[1][2])
398 + cos6 * RationalFunction(x, coeffs[1][3]);
399
400 return RationalFunction(x, safety) * ep->a * exp(c0 + c1 * dEta);
401 }
402}
403
404/* Non-reentrant code to append points to a Bezier path
405 * Assume initial call to moveTo to initialize, followed by
406 * calls to curveTo and lineTo, and finished with endPath.
407 */
408static int bufsize;
409
410static void moveTo(Ppolyline_t * path, double x, double y)
411{
412 bufsize = 100;
413 path->ps = N_NEW(bufsize, pointf);
414 path->ps[0].x = x;
415 path->ps[0].y = y;
416 path->pn = 1;
417}
418
419static void
420curveTo(Ppolyline_t * path, double x1, double y1,
421 double x2, double y2, double x3, double y3)
422{
423 if (path->pn + 3 >= bufsize) {
424 bufsize *= 2;
425 path->ps = realloc(path->ps, bufsize * sizeof(pointf));
426 }
427 path->ps[path->pn].x = x1;
428 path->ps[path->pn++].y = y1;
429 path->ps[path->pn].x = x2;
430 path->ps[path->pn++].y = y2;
431 path->ps[path->pn].x = x3;
432 path->ps[path->pn++].y = y3;
433}
434
435static void lineTo(Ppolyline_t * path, double x, double y)
436{
437 pointf curp = path->ps[path->pn - 1];
438 curveTo(path, curp.x, curp.y, x, y, x, y);
439}
440
441static void endPath(Ppolyline_t * path, boolean close)
442{
443 if (close) {
444 pointf p0 = path->ps[0];
445 lineTo(path, p0.x, p0.y);
446 }
447
448 path->ps = realloc(path->ps, path->pn * sizeof(pointf));
449 bufsize = 0;
450}
451
452/* genEllipticPath:
453 * Approximate an elliptical arc via Beziers of given degree
454 * threshold indicates quality of approximation
455 * if isSlice is true, the path begins and ends with line segments
456 * to the center of the ellipse.
457 * Returned path must be freed by the caller.
458 */
459static Ppolyline_t *genEllipticPath(ellipse_t * ep, int degree,
460 double threshold, boolean isSlice)
461{
462 double dEta;
463 double etaB;
464 double cosEtaB;
465 double sinEtaB;
466 double aCosEtaB;
467 double bSinEtaB;
468 double aSinEtaB;
469 double bCosEtaB;
470 double xB;
471 double yB;
472 double xBDot;
473 double yBDot;
474 double t;
475 double alpha;
476 Ppolyline_t *path = NEW(Ppolyline_t);
477
478 // find the number of Bezier curves needed
479 boolean found = FALSE;
480 int i, n = 1;
481 while ((!found) && (n < 1024)) {
482 double dEta = (ep->eta2 - ep->eta1) / n;
483 if (dEta <= 0.5 * M_PI) {
484 double etaB = ep->eta1;
485 found = TRUE;
486 for (i = 0; found && (i < n); ++i) {
487 double etaA = etaB;
488 etaB += dEta;
489 found =
490 (estimateError(ep, degree, etaA, etaB) <= threshold);
491 }
492 }
493 n = n << 1;
494 }
495
496 dEta = (ep->eta2 - ep->eta1) / n;
497 etaB = ep->eta1;
498
499 cosEtaB = cos(etaB);
500 sinEtaB = sin(etaB);
501 aCosEtaB = ep->a * cosEtaB;
502 bSinEtaB = ep->b * sinEtaB;
503 aSinEtaB = ep->a * sinEtaB;
504 bCosEtaB = ep->b * cosEtaB;
505 xB = ep->cx + aCosEtaB * ep->cosTheta - bSinEtaB * ep->sinTheta;
506 yB = ep->cy + aCosEtaB * ep->sinTheta + bSinEtaB * ep->cosTheta;
507 xBDot = -aSinEtaB * ep->cosTheta - bCosEtaB * ep->sinTheta;
508 yBDot = -aSinEtaB * ep->sinTheta + bCosEtaB * ep->cosTheta;
509
510 if (isSlice) {
511 moveTo(path, ep->cx, ep->cy);
512 lineTo(path, xB, yB);
513 } else {
514 moveTo(path, xB, yB);
515 }
516
517 t = tan(0.5 * dEta);
518 alpha = sin(dEta) * (sqrt(4 + 3 * t * t) - 1) / 3;
519
520 for (i = 0; i < n; ++i) {
521
522 double xA = xB;
523 double yA = yB;
524 double xADot = xBDot;
525 double yADot = yBDot;
526
527 etaB += dEta;
528 cosEtaB = cos(etaB);
529 sinEtaB = sin(etaB);
530 aCosEtaB = ep->a * cosEtaB;
531 bSinEtaB = ep->b * sinEtaB;
532 aSinEtaB = ep->a * sinEtaB;
533 bCosEtaB = ep->b * cosEtaB;
534 xB = ep->cx + aCosEtaB * ep->cosTheta - bSinEtaB * ep->sinTheta;
535 yB = ep->cy + aCosEtaB * ep->sinTheta + bSinEtaB * ep->cosTheta;
536 xBDot = -aSinEtaB * ep->cosTheta - bCosEtaB * ep->sinTheta;
537 yBDot = -aSinEtaB * ep->sinTheta + bCosEtaB * ep->cosTheta;
538
539 if (degree == 1) {
540 lineTo(path, xB, yB);
541#if DO_QUAD
542 } else if (degree == 2) {
543 double k = (yBDot * (xB - xA) - xBDot * (yB - yA))
544 / (xADot * yBDot - yADot * xBDot);
545 quadTo(path, (xA + k * xADot), (yA + k * yADot), xB, yB);
546#endif
547 } else {
548 curveTo(path, (xA + alpha * xADot), (yA + alpha * yADot),
549 (xB - alpha * xBDot), (yB - alpha * yBDot), xB, yB);
550 }
551
552 }
553
554 endPath(path, isSlice);
555
556 return path;
557}
558
559/* ellipticWedge:
560 * Return a cubic Bezier for an elliptical wedge, with center ctr, x and y
561 * semi-axes xsemi and ysemi, start angle angle0 and end angle angle1.
562 * This includes beginning and ending line segments to the ellipse center.
563 * Calling function must free storage of returned path.
564 */
565Ppolyline_t *ellipticWedge(pointf ctr, double xsemi, double ysemi,
566 double angle0, double angle1)
567{
568 ellipse_t ell;
569 Ppolyline_t *pp;
570
571 initEllipse(&ell, ctr.x, ctr.y, xsemi, ysemi, 0, angle0, angle1);
572 pp = genEllipticPath(&ell, 3, 0.00001, 1);
573 return pp;
574}
575
576#ifdef STANDALONE
577main()
578{
579 ellipse_t ell;
580 Ppolyline_t *pp;
581 int i;
582
583 initEllipse(&ell, 200, 200, 100, 50, 0, M_PI / 4, 3 * M_PI / 2);
584 pp = genEllipticPath(&ell, 3, 0.00001, 1);
585
586 printf("newpath %.02lf %.02lf moveto\n", pp->ps[0].x, pp->ps[0].y);
587 for (i = 1; i < pp->pn; i += 3) {
588 printf("%.02lf %.02lf %.02lf %.02lf %.02lf %.02lf curveto\n",
589 pp->ps[i].x, pp->ps[i].y,
590 pp->ps[i + 1].x, pp->ps[i + 1].y,
591 pp->ps[i + 2].x, pp->ps[i + 2].y);
592 }
593 printf("stroke showpage\n");
594
595}
596#endif
597