Graphviz 13.0.0~dev.20241220.2304
Loading...
Searching...
No Matches
ellipse.c
Go to the documentation of this file.
1
3/*************************************************************************
4 * Copyright (c) 2012 AT&T Intellectual Property
5 * All rights reserved. This program and the accompanying materials
6 * are made available under the terms of the Eclipse Public License v1.0
7 * which accompanies this distribution, and is available at
8 * https://www.eclipse.org/legal/epl-v10.html
9 *
10 * Contributors: Details at https://graphviz.org
11 *************************************************************************/
12
13/* This code is derived from the Java implementation by Luc Maisonobe */
14/* Copyright (c) 2003-2004, Luc Maisonobe
15 * All rights reserved.
16 *
17 * Redistribution and use in source and binary forms, with
18 * or without modification, are permitted provided that
19 * the following conditions are met:
20 *
21 * Redistributions of source code must retain the
22 * above copyright notice, this list of conditions and
23 * the following disclaimer.
24 * Redistributions in binary form must reproduce the
25 * above copyright notice, this list of conditions and
26 * the following disclaimer in the documentation
27 * and/or other materials provided with the
28 * distribution.
29 * Neither the names of spaceroots.org, spaceroots.com
30 * nor the names of their contributors may be used to
31 * endorse or promote products derived from this
32 * software without specific prior written permission.
33 *
34 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
35 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
36 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
37 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
38 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
39 * THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
40 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
41 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
42 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
43 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
44 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
45 * IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
46 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
47 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
48 * POSSIBILITY OF SUCH DAMAGE.
49 */
50
51#include <assert.h>
52#include <cgraph/list.h>
53#include <limits.h>
54#include <math.h>
55#include <stdbool.h>
56#include <common/render.h>
57#include <pathplan/pathplan.h>
58#include <util/alloc.h>
59
60#define TWOPI (2*M_PI)
61
62typedef struct {
63 double cx, cy; /* center */
64 double a, b; /* semi-major and -minor axes */
65
66 /* Start and end angles of the arc. */
67 double eta1, eta2;
68} ellipse_t;
69
70static void initEllipse(ellipse_t * ep, double cx, double cy, double a,
71 double b, double lambda1, double lambda2) {
72 ep->cx = cx;
73 ep->cy = cy;
74 ep->a = a;
75 ep->b = b;
76
77 ep->eta1 = atan2(sin(lambda1) / b, cos(lambda1) / a);
78 ep->eta2 = atan2(sin(lambda2) / b, cos(lambda2) / a);
79
80 // make sure we have eta1 <= eta2 <= eta1 + 2*PI
81 ep->eta2 -= TWOPI * floor((ep->eta2 - ep->eta1) / TWOPI);
82
83 // the preceding correction fails if we have exactly eta2 - eta1 = 2*PI
84 // it reduces the interval to zero length
85 if (lambda2 - lambda1 > M_PI && ep->eta2 - ep->eta1 < M_PI) {
86 ep->eta2 += TWOPI;
87 }
88}
89
90typedef double erray_t[2][4][4];
91
92 // coefficients for error estimation
93 // while using cubic Bezier curves for approximation
94 // 0 < b/a < 1/4
96 {
97 {3.85268, -21.229, -0.330434, 0.0127842},
98 {-1.61486, 0.706564, 0.225945, 0.263682},
99 {-0.910164, 0.388383, 0.00551445, 0.00671814},
100 {-0.630184, 0.192402, 0.0098871, 0.0102527}
101 },
102 {
103 {-0.162211, 9.94329, 0.13723, 0.0124084},
104 {-0.253135, 0.00187735, 0.0230286, 0.01264},
105 {-0.0695069, -0.0437594, 0.0120636, 0.0163087},
106 {-0.0328856, -0.00926032, -0.00173573, 0.00527385}
107 }
108};
109
110 // coefficients for error estimation
111 // while using cubic Bezier curves for approximation
112 // 1/4 <= b/a <= 1
114 {
115 {0.0899116, -19.2349, -4.11711, 0.183362},
116 {0.138148, -1.45804, 1.32044, 1.38474},
117 {0.230903, -0.450262, 0.219963, 0.414038},
118 {0.0590565, -0.101062, 0.0430592, 0.0204699}
119 },
120 {
121 {0.0164649, 9.89394, 0.0919496, 0.00760802},
122 {0.0191603, -0.0322058, 0.0134667, -0.0825018},
123 {0.0156192, -0.017535, 0.00326508, -0.228157},
124 {-0.0236752, 0.0405821, -0.0173086, 0.176187}
125 }
126};
127
128 // safety factor to convert the "best" error approximation
129 // into a "max bound" error
130static double safety3[] = {
131 0.001, 4.98, 0.207, 0.0067
132};
133
134/* Compute the value of a rational function.
135 * This method handles rational functions where the numerator is
136 * quadratic and the denominator is linear
137 */
138#define RationalFunction(x,c) ((x * (x * c[0] + c[1]) + c[2]) / (x + c[3]))
139
140/* Estimate the approximation error for a sub-arc of the instance.
141 * tA and tB give the start and end angle of the subarc
142 * Returns upper bound of the approximation error between the Bezier
143 * curve and the real ellipse
144 */
145static double estimateError(ellipse_t *ep, double etaA, double etaB) {
146 double c0, c1, eta = 0.5 * (etaA + etaB);
147
148 double x = ep->b / ep->a;
149 double dEta = etaB - etaA;
150 double cos2 = cos(2 * eta);
151 double cos4 = cos(4 * eta);
152 double cos6 = cos(6 * eta);
153
154 // select the right coefficient's set according to b/a
155 double (*coeffs)[4][4];
156 coeffs = x < 0.25 ? coeffs3Low : coeffs3High;
157
158 c0 = RationalFunction(x, coeffs[0][0])
159 + cos2 * RationalFunction(x, coeffs[0][1])
160 + cos4 * RationalFunction(x, coeffs[0][2])
161 + cos6 * RationalFunction(x, coeffs[0][3]);
162
163 c1 = RationalFunction(x, coeffs[1][0])
164 + cos2 * RationalFunction(x, coeffs[1][1])
165 + cos4 * RationalFunction(x, coeffs[1][2])
166 + cos6 * RationalFunction(x, coeffs[1][3]);
167
168 return RationalFunction(x, safety3) * ep->a * exp(c0 + c1 * dEta);
169}
170
171DEFINE_LIST(bezier_path, pointf)
172
173/* append points to a Bezier path
174 * Assume initial call to moveTo to initialize, followed by
175 * calls to curveTo and lineTo, and finished with endPath.
176 */
177
178static void moveTo(bezier_path_t *polypath, double x, double y) {
179 bezier_path_append(polypath, (pointf){.x = x, .y = y});
180}
181
182static void curveTo(bezier_path_t *polypath, double x1, double y1, double x2,
183 double y2, double x3, double y3) {
184 bezier_path_append(polypath, (pointf){.x = x1, .y = y1});
185 bezier_path_append(polypath, (pointf){.x = x2, .y = y2});
186 bezier_path_append(polypath, (pointf){.x = x3, .y = y3});
187}
188
189static void lineTo(bezier_path_t *polypath, double x, double y) {
190 const pointf curp = bezier_path_get(polypath, bezier_path_size(polypath) - 1);
191 curveTo(polypath, curp.x, curp.y, x, y, x, y);
192}
193
194static void endPath(bezier_path_t *polypath) {
195 const pointf p0 = bezier_path_get(polypath, 0);
196 lineTo(polypath, p0.x, p0.y);
197}
198
199/* genEllipticPath:
200 * Approximate an elliptical arc via Beziers of degree 3
201 * The path begins and ends with line segments to the center of the ellipse.
202 * Returned path must be freed by the caller.
203 */
205 double dEta;
206 double etaB;
207 double cosEtaB;
208 double sinEtaB;
209 double aCosEtaB;
210 double bSinEtaB;
211 double aSinEtaB;
212 double bCosEtaB;
213 double xB;
214 double yB;
215 double xBDot;
216 double yBDot;
217 double t;
218 double alpha;
219 Ppolyline_t *polypath = gv_alloc(sizeof(Ppolyline_t));
220
221 static const double THRESHOLD = 0.00001; // quality of approximation
222
223 // find the number of Bezier curves needed
224 bool found = false;
225 int i, n = 1;
226 while (!found && n < 1024) {
227 double diffEta = (ep->eta2 - ep->eta1) / n;
228 if (diffEta <= 0.5 * M_PI) {
229 double etaOne = ep->eta1;
230 found = true;
231 for (i = 0; found && i < n; ++i) {
232 double etaA = etaOne;
233 etaOne += diffEta;
234 found = estimateError(ep, etaA, etaOne) <= THRESHOLD;
235 }
236 }
237 n = n << 1;
238 }
239
240 dEta = (ep->eta2 - ep->eta1) / n;
241 etaB = ep->eta1;
242
243 cosEtaB = cos(etaB);
244 sinEtaB = sin(etaB);
245 aCosEtaB = ep->a * cosEtaB;
246 bSinEtaB = ep->b * sinEtaB;
247 aSinEtaB = ep->a * sinEtaB;
248 bCosEtaB = ep->b * cosEtaB;
249 xB = ep->cx + aCosEtaB;
250 yB = ep->cy + bSinEtaB;
251 xBDot = -aSinEtaB;
252 yBDot = bCosEtaB;
253
254 bezier_path_t bezier_path = {0};
255 moveTo(&bezier_path, ep->cx, ep->cy);
256 lineTo(&bezier_path, xB, yB);
257
258 t = tan(0.5 * dEta);
259 alpha = sin(dEta) * (sqrt(4 + 3 * t * t) - 1) / 3;
260
261 for (i = 0; i < n; ++i) {
262
263 double xA = xB;
264 double yA = yB;
265 double xADot = xBDot;
266 double yADot = yBDot;
267
268 etaB += dEta;
269 cosEtaB = cos(etaB);
270 sinEtaB = sin(etaB);
271 aCosEtaB = ep->a * cosEtaB;
272 bSinEtaB = ep->b * sinEtaB;
273 aSinEtaB = ep->a * sinEtaB;
274 bCosEtaB = ep->b * cosEtaB;
275 xB = ep->cx + aCosEtaB;
276 yB = ep->cy + bSinEtaB;
277 xBDot = -aSinEtaB;
278 yBDot = bCosEtaB;
279
280 curveTo(&bezier_path, xA + alpha * xADot, yA + alpha * yADot,
281 xB - alpha * xBDot, yB - alpha * yBDot, xB, yB);
282 }
283
284 endPath(&bezier_path);
285
286 polypath->pn = bezier_path_size(&bezier_path);
287 polypath->ps = bezier_path_detach(&bezier_path);
288
289 return polypath;
290}
291
292/* ellipticWedge:
293 * Return a cubic Bezier for an elliptical wedge, with center ctr, x and y
294 * semi-axes xsemi and ysemi, start angle angle0 and end angle angle1.
295 * This includes beginning and ending line segments to the ellipse center.
296 * Calling function must free storage of returned path.
297 */
298Ppolyline_t *ellipticWedge(pointf ctr, double xsemi, double ysemi,
299 double angle0, double angle1)
300{
301 ellipse_t ell;
302 Ppolyline_t *pp;
303
304 initEllipse(&ell, ctr.x, ctr.y, xsemi, ysemi, angle0, angle1);
305 pp = genEllipticPath(&ell);
306 return pp;
307}
Memory allocation wrappers that exit on failure.
static void * gv_alloc(size_t size)
Definition alloc.h:47
#define M_PI
Definition arith.h:41
static void lineTo(bezier_path_t *polypath, double x, double y)
Definition ellipse.c:189
static Ppolyline_t * genEllipticPath(ellipse_t *ep)
Definition ellipse.c:204
static double safety3[]
Definition ellipse.c:130
static erray_t coeffs3Low
Definition ellipse.c:95
#define TWOPI
Definition ellipse.c:60
static void initEllipse(ellipse_t *ep, double cx, double cy, double a, double b, double lambda1, double lambda2)
Definition ellipse.c:70
static erray_t coeffs3High
Definition ellipse.c:113
Ppolyline_t * ellipticWedge(pointf ctr, double xsemi, double ysemi, double angle0, double angle1)
Definition ellipse.c:298
static double estimateError(ellipse_t *ep, double etaA, double etaB)
Definition ellipse.c:145
static void endPath(bezier_path_t *polypath)
Definition ellipse.c:194
static void curveTo(bezier_path_t *polypath, double x1, double y1, double x2, double y2, double x3, double y3)
Definition ellipse.c:182
#define RationalFunction(x, c)
Definition ellipse.c:138
static void moveTo(bezier_path_t *polypath, double x, double y)
Definition ellipse.c:178
double erray_t[2][4][4]
Definition ellipse.c:90
#define DEFINE_LIST(name, type)
Definition list.h:26
finds and smooths shortest paths
#define alpha
Definition shapes.c:4058
size_t pn
Definition pathgeom.h:47
Ppoint_t * ps
Definition pathgeom.h:46
double a
Definition ellipse.c:64
double cx
Definition ellipse.c:63
double b
Definition ellipse.c:64
double eta1
Definition ellipse.c:67
double eta2
Definition ellipse.c:67
double cy
Definition ellipse.c:63
double x
Definition geom.h:29
double y
Definition geom.h:29