1 /*
2 ppedit - A pattern plate editor for Spiro splines.
3 Copyright (C) 2007 Raph Levien
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 02110-1301, USA.
20 */
21 /* C implementation of third-order polynomial spirals. */
23 #include <math.h>
24 #include <stdlib.h>
25 #include <string.h>
27 #include "bezctx_intf.h"
28 #include "spiro.h"
30 struct spiro_seg_s {
31 double x;
32 double y;
33 char ty;
34 double bend_th;
35 double ks[4];
36 double seg_ch;
37 double seg_th;
38 double l;
39 };
41 typedef struct {
42 double a[11]; /* band-diagonal matrix */
43 double al[5]; /* lower part of band-diagonal decomposition */
44 } bandmat;
46 #ifndef M_PI
47 #define M_PI 3.14159265358979323846 /* pi */
48 #endif
50 int n = 4;
52 #ifndef ORDER
53 #define ORDER 12
54 #endif
56 /* Integrate polynomial spiral curve over range -.5 .. .5. */
57 void
58 integrate_spiro(const double ks[4], double xy[2])
59 {
60 #if 0
61 int n = 1024;
62 #endif
63 double th1 = ks[0];
64 double th2 = .5 * ks[1];
65 double th3 = (1./6) * ks[2];
66 double th4 = (1./24) * ks[3];
67 double x, y;
68 double ds = 1. / n;
69 double ds2 = ds * ds;
70 double ds3 = ds2 * ds;
71 double k0 = ks[0] * ds;
72 double k1 = ks[1] * ds;
73 double k2 = ks[2] * ds;
74 double k3 = ks[3] * ds;
75 int i;
76 double s = .5 * ds - .5;
78 x = 0;
79 y = 0;
81 for (i = 0; i < n; i++) {
83 #if ORDER > 2
84 double u, v;
85 double km0, km1, km2, km3;
87 if (n == 1) {
88 km0 = k0;
89 km1 = k1 * ds;
90 km2 = k2 * ds2;
91 } else {
92 km0 = (((1./6) * k3 * s + .5 * k2) * s + k1) * s + k0;
93 km1 = ((.5 * k3 * s + k2) * s + k1) * ds;
94 km2 = (k3 * s + k2) * ds2;
95 }
96 km3 = k3 * ds3;
97 #endif
99 {
101 #if ORDER == 4
102 double km0_2 = km0 * km0;
103 u = 24 - km0_2;
104 v = km1;
105 #endif
107 #if ORDER == 6
108 double km0_2 = km0 * km0;
109 double km0_4 = km0_2 * km0_2;
110 u = 24 - km0_2 + (km0_4 - 4 * km0 * km2 - 3 * km1 * km1) * (1./80);
111 v = km1 + (km3 - 6 * km0_2 * km1) * (1./80);
112 #endif
114 #if ORDER == 8
115 double t1_1 = km0;
116 double t1_2 = .5 * km1;
117 double t1_3 = (1./6) * km2;
118 double t1_4 = (1./24) * km3;
119 double t2_2 = t1_1 * t1_1;
120 double t2_3 = 2 * (t1_1 * t1_2);
121 double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
122 double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
123 double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
124 double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
125 double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
126 double t4_4 = t2_2 * t2_2;
127 double t4_5 = 2 * (t2_2 * t2_3);
128 double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
129 double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
130 double t6_6 = t4_4 * t2_2;
131 u = 1;
132 v = 0;
133 v += (1./12) * t1_2 + (1./80) * t1_4;
134 u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6;
135 v -= (1./480) * t3_4 + (1./2688) * t3_6;
136 u += (1./1920) * t4_4 + (1./10752) * t4_6;
137 v += (1./53760) * t5_6;
138 u -= (1./322560) * t6_6;
139 #endif
141 #if ORDER == 10
142 double t1_1 = km0;
143 double t1_2 = .5 * km1;
144 double t1_3 = (1./6) * km2;
145 double t1_4 = (1./24) * km3;
146 double t2_2 = t1_1 * t1_1;
147 double t2_3 = 2 * (t1_1 * t1_2);
148 double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
149 double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
150 double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
151 double t2_7 = 2 * (t1_3 * t1_4);
152 double t2_8 = t1_4 * t1_4;
153 double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
154 double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
155 double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
156 double t4_4 = t2_2 * t2_2;
157 double t4_5 = 2 * (t2_2 * t2_3);
158 double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
159 double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
160 double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
161 double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
162 double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
163 double t6_6 = t4_4 * t2_2;
164 double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
165 double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
166 double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
167 double t8_8 = t6_6 * t2_2;
168 u = 1;
169 v = 0;
170 v += (1./12) * t1_2 + (1./80) * t1_4;
171 u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
172 v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8;
173 u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8;
174 v += (1./53760) * t5_6 + (1./276480) * t5_8;
175 u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8;
176 v -= (1./1.16122e+07) * t7_8;
177 u += (1./9.28973e+07) * t8_8;
178 #endif
180 #if ORDER == 12
181 double t1_1 = km0;
182 double t1_2 = .5 * km1;
183 double t1_3 = (1./6) * km2;
184 double t1_4 = (1./24) * km3;
185 double t2_2 = t1_1 * t1_1;
186 double t2_3 = 2 * (t1_1 * t1_2);
187 double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
188 double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
189 double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
190 double t2_7 = 2 * (t1_3 * t1_4);
191 double t2_8 = t1_4 * t1_4;
192 double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
193 double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
194 double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
195 double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
196 double t4_4 = t2_2 * t2_2;
197 double t4_5 = 2 * (t2_2 * t2_3);
198 double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
199 double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
200 double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
201 double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
202 double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
203 double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
204 double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
205 double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
206 double t6_6 = t4_4 * t2_2;
207 double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
208 double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
209 double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
210 double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
211 double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
212 double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
213 double t8_8 = t6_6 * t2_2;
214 double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
215 double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
216 double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
217 double t10_10 = t8_8 * t2_2;
218 u = 1;
219 v = 0;
220 v += (1./12) * t1_2 + (1./80) * t1_4;
221 u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
222 v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10;
223 u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10;
224 v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10;
225 u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10;
226 v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10;
227 u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10;
228 v += (1./4.08748e+09) * t9_10;
229 u -= (1./4.08748e+10) * t10_10;
230 #endif
232 #if ORDER == 14
233 double t1_1 = km0;
234 double t1_2 = .5 * km1;
235 double t1_3 = (1./6) * km2;
236 double t1_4 = (1./24) * km3;
237 double t2_2 = t1_1 * t1_1;
238 double t2_3 = 2 * (t1_1 * t1_2);
239 double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
240 double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
241 double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
242 double t2_7 = 2 * (t1_3 * t1_4);
243 double t2_8 = t1_4 * t1_4;
244 double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
245 double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
246 double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
247 double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
248 double t3_12 = t2_8 * t1_4;
249 double t4_4 = t2_2 * t2_2;
250 double t4_5 = 2 * (t2_2 * t2_3);
251 double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
252 double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
253 double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
254 double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
255 double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
256 double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
257 double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
258 double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
259 double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
260 double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
261 double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
262 double t6_6 = t4_4 * t2_2;
263 double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
264 double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
265 double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
266 double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
267 double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
268 double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
269 double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
270 double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
271 double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
272 double t8_8 = t6_6 * t2_2;
273 double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
274 double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
275 double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
276 double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
277 double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
278 double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
279 double t10_10 = t8_8 * t2_2;
280 double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
281 double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
282 double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
283 double t12_12 = t10_10 * t2_2;
284 u = 1;
285 v = 0;
286 v += (1./12) * t1_2 + (1./80) * t1_4;
287 u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
288 v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10 + (1./319488) * t3_12;
289 u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10 + (1./1.27795e+06) * t4_12;
290 v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10 + (1./6.38976e+06) * t5_12;
291 u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10 + (1./3.83386e+07) * t6_12;
292 v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10 + (1./2.6837e+08) * t7_12;
293 u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10 + (1./2.14696e+09) * t8_12;
294 v += (1./4.08748e+09) * t9_10 + (1./1.93226e+10) * t9_12;
295 u -= (1./4.08748e+10) * t10_10 + (1./1.93226e+11) * t10_12;
296 v -= (1./2.12549e+12) * t11_12;
297 u += (1./2.55059e+13) * t12_12;
298 #endif
300 #if ORDER == 16
301 double t1_1 = km0;
302 double t1_2 = .5 * km1;
303 double t1_3 = (1./6) * km2;
304 double t1_4 = (1./24) * km3;
305 double t2_2 = t1_1 * t1_1;
306 double t2_3 = 2 * (t1_1 * t1_2);
307 double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
308 double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
309 double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
310 double t2_7 = 2 * (t1_3 * t1_4);
311 double t2_8 = t1_4 * t1_4;
312 double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
313 double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
314 double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
315 double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
316 double t3_12 = t2_8 * t1_4;
317 double t4_4 = t2_2 * t2_2;
318 double t4_5 = 2 * (t2_2 * t2_3);
319 double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
320 double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
321 double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
322 double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
323 double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
324 double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
325 double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
326 double t4_13 = 2 * (t2_5 * t2_8 + t2_6 * t2_7);
327 double t4_14 = 2 * (t2_6 * t2_8) + t2_7 * t2_7;
328 double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
329 double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
330 double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
331 double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
332 double t5_14 = t4_10 * t1_4 + t4_11 * t1_3 + t4_12 * t1_2 + t4_13 * t1_1;
333 double t6_6 = t4_4 * t2_2;
334 double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
335 double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
336 double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
337 double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
338 double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
339 double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
340 double t6_13 = t4_5 * t2_8 + t4_6 * t2_7 + t4_7 * t2_6 + t4_8 * t2_5 + t4_9 * t2_4 + t4_10 * t2_3 + t4_11 * t2_2;
341 double t6_14 = t4_6 * t2_8 + t4_7 * t2_7 + t4_8 * t2_6 + t4_9 * t2_5 + t4_10 * t2_4 + t4_11 * t2_3 + t4_12 * t2_2;
342 double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
343 double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
344 double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
345 double t7_14 = t6_10 * t1_4 + t6_11 * t1_3 + t6_12 * t1_2 + t6_13 * t1_1;
346 double t8_8 = t6_6 * t2_2;
347 double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
348 double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
349 double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
350 double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
351 double t8_13 = t6_6 * t2_7 + t6_7 * t2_6 + t6_8 * t2_5 + t6_9 * t2_4 + t6_10 * t2_3 + t6_11 * t2_2;
352 double t8_14 = t6_6 * t2_8 + t6_7 * t2_7 + t6_8 * t2_6 + t6_9 * t2_5 + t6_10 * t2_4 + t6_11 * t2_3 + t6_12 * t2_2;
353 double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
354 double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
355 double t9_14 = t8_10 * t1_4 + t8_11 * t1_3 + t8_12 * t1_2 + t8_13 * t1_1;
356 double t10_10 = t8_8 * t2_2;
357 double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
358 double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
359 double t10_13 = t8_8 * t2_5 + t8_9 * t2_4 + t8_10 * t2_3 + t8_11 * t2_2;
360 double t10_14 = t8_8 * t2_6 + t8_9 * t2_5 + t8_10 * t2_4 + t8_11 * t2_3 + t8_12 * t2_2;
361 double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
362 double t11_14 = t10_10 * t1_4 + t10_11 * t1_3 + t10_12 * t1_2 + t10_13 * t1_1;
363 double t12_12 = t10_10 * t2_2;
364 double t12_13 = t10_10 * t2_3 + t10_11 * t2_2;
365 double t12_14 = t10_10 * t2_4 + t10_11 * t2_3 + t10_12 * t2_2;
366 double t13_14 = t12_12 * t1_2 + t12_13 * t1_1;
367 double t14_14 = t12_12 * t2_2;
368 u = 1;
369 u -= 1./24 * t2_2 + 1./160 * t2_4 + 1./896 * t2_6 + 1./4608 * t2_8;
370 u += 1./1920 * t4_4 + 1./10752 * t4_6 + 1./55296 * t4_8 + 1./270336 * t4_10 + 1./1277952 * t4_12 + 1./5898240 * t4_14;
371 u -= 1./322560 * t6_6 + 1./1658880 * t6_8 + 1./8110080 * t6_10 + 1./38338560 * t6_12 + 1./176947200 * t6_14;
372 u += 1./92897280 * t8_8 + 1./454164480 * t8_10 + 4.6577500191e-10 * t8_12 + 1.0091791708e-10 * t8_14;
373 u -= 2.4464949595e-11 * t10_10 + 5.1752777990e-12 * t10_12 + 1.1213101898e-12 * t10_14;
374 u += 3.9206649992e-14 * t12_12 + 8.4947741650e-15 * t12_14;
375 u -= 4.6674583324e-17 * t14_14;
376 v = 0;
377 v += 1./12 * t1_2 + 1./80 * t1_4;
378 v -= 1./480 * t3_4 + 1./2688 * t3_6 + 1./13824 * t3_8 + 1./67584 * t3_10 + 1./319488 * t3_12;
379 v += 1./53760 * t5_6 + 1./276480 * t5_8 + 1./1351680 * t5_10 + 1./6389760 * t5_12 + 1./29491200 * t5_14;
380 v -= 1./11612160 * t7_8 + 1./56770560 * t7_10 + 1./268369920 * t7_12 + 8.0734333664e-10 * t7_14;
381 v += 2.4464949595e-10 * t9_10 + 5.1752777990e-11 * t9_12 + 1.1213101898e-11 * t9_14;
382 v -= 4.7047979991e-13 * t11_12 + 1.0193728998e-13 * t11_14;
383 v += 6.5344416654e-16 * t13_14;
384 #endif
386 }
388 if (n == 1) {
389 #if ORDER == 2
390 x = 1;
391 y = 0;
392 #else
393 x = u;
394 y = v;
395 #endif
396 } else {
397 double th = (((th4 * s + th3) * s + th2) * s + th1) * s;
398 double cth = cos(th);
399 double sth = sin(th);
401 #if ORDER == 2
402 x += cth;
403 y += sth;
404 #else
405 x += cth * u - sth * v;
406 y += cth * v + sth * u;
407 #endif
408 s += ds;
409 }
410 }
412 #if ORDER == 4 || ORDER == 6
413 xy[0] = x * (1./24 * ds);
414 xy[1] = y * (1./24 * ds);
415 #else
416 xy[0] = x * ds;
417 xy[1] = y * ds;
418 #endif
419 }
421 static double
422 compute_ends(const double ks[4], double ends[2][4], double seg_ch)
423 {
424 double xy[2];
425 double ch, th;
426 double l, l2, l3;
427 double th_even, th_odd;
428 double k0_even, k0_odd;
429 double k1_even, k1_odd;
430 double k2_even, k2_odd;
432 integrate_spiro(ks, xy);
433 ch = hypot(xy[0], xy[1]);
434 th = atan2(xy[1], xy[0]);
435 l = ch / seg_ch;
437 th_even = .5 * ks[0] + (1./48) * ks[2];
438 th_odd = .125 * ks[1] + (1./384) * ks[3] - th;
439 ends[0][0] = th_even - th_odd;
440 ends[1][0] = th_even + th_odd;
441 k0_even = l * (ks[0] + .125 * ks[2]);
442 k0_odd = l * (.5 * ks[1] + (1./48) * ks[3]);
443 ends[0][1] = k0_even - k0_odd;
444 ends[1][1] = k0_even + k0_odd;
445 l2 = l * l;
446 k1_even = l2 * (ks[1] + .125 * ks[3]);
447 k1_odd = l2 * .5 * ks[2];
448 ends[0][2] = k1_even - k1_odd;
449 ends[1][2] = k1_even + k1_odd;
450 l3 = l2 * l;
451 k2_even = l3 * ks[2];
452 k2_odd = l3 * .5 * ks[3];
453 ends[0][3] = k2_even - k2_odd;
454 ends[1][3] = k2_even + k2_odd;
456 return l;
457 }
459 static void
460 compute_pderivs(const spiro_seg *s, double ends[2][4], double derivs[4][2][4],
461 int jinc)
462 {
463 double recip_d = 2e6;
464 double delta = 1./ recip_d;
465 double try_ks[4];
466 double try_ends[2][4];
467 int i, j, k;
469 compute_ends(s->ks, ends, s->seg_ch);
470 for (i = 0; i < jinc; i++) {
471 for (j = 0; j < 4; j++)
472 try_ks[j] = s->ks[j];
473 try_ks[i] += delta;
474 compute_ends(try_ks, try_ends, s->seg_ch);
475 for (k = 0; k < 2; k++)
476 for (j = 0; j < 4; j++)
477 derivs[j][k][i] = recip_d * (try_ends[k][j] - ends[k][j]);
478 }
479 }
481 static double
482 mod_2pi(double th)
483 {
484 double u = th / (2 * M_PI);
485 return 2 * M_PI * (u - floor(u + 0.5));
486 }
488 static spiro_seg *
489 setup_path(const spiro_cp *src, int n)
490 {
491 int n_seg = src[0].ty == '{' ? n - 1 : n;
492 spiro_seg *r = (spiro_seg *)malloc((n_seg + 1) * sizeof(spiro_seg));
493 int i;
494 int ilast;
496 for (i = 0; i < n_seg; i++) {
497 r[i].x = src[i].x;
498 r[i].y = src[i].y;
499 r[i].ty = src[i].ty;
500 r[i].ks[0] = 0.;
501 r[i].ks[1] = 0.;
502 r[i].ks[2] = 0.;
503 r[i].ks[3] = 0.;
504 }
505 r[n_seg].x = src[n_seg % n].x;
506 r[n_seg].y = src[n_seg % n].y;
507 r[n_seg].ty = src[n_seg % n].ty;
509 for (i = 0; i < n_seg; i++) {
510 double dx = r[i + 1].x - r[i].x;
511 double dy = r[i + 1].y - r[i].y;
512 r[i].seg_ch = hypot(dx, dy);
513 r[i].seg_th = atan2(dy, dx);
514 }
516 ilast = n_seg - 1;
517 for (i = 0; i < n_seg; i++) {
518 if (r[i].ty == '{' || r[i].ty == '}' || r[i].ty == 'v')
519 r[i].bend_th = 0.;
520 else
521 r[i].bend_th = mod_2pi(r[i].seg_th - r[ilast].seg_th);
522 ilast = i;
523 }
524 return r;
525 }
527 static void
528 bandec11(bandmat *m, int *perm, int n)
529 {
530 int i, j, k;
531 int l;
533 /* pack top triangle to the left. */
534 for (i = 0; i < 5; i++) {
535 for (j = 0; j < i + 6; j++)
536 m[i].a[j] = m[i].a[j + 5 - i];
537 for (; j < 11; j++)
538 m[i].a[j] = 0.;
539 }
540 l = 5;
541 for (k = 0; k < n; k++) {
542 int pivot = k;
543 double pivot_val = m[k].a[0];
544 double pivot_scale;
546 l = l < n ? l + 1 : n;
548 for (j = k + 1; j < l; j++)
549 if (fabs(m[j].a[0]) > fabs(pivot_val)) {
550 pivot_val = m[j].a[0];
551 pivot = j;
552 }
554 perm[k] = pivot;
555 if (pivot != k) {
556 for (j = 0; j < 11; j++) {
557 double tmp = m[k].a[j];
558 m[k].a[j] = m[pivot].a[j];
559 m[pivot].a[j] = tmp;
560 }
561 }
563 if (fabs(pivot_val) < 1e-12) pivot_val = 1e-12;
564 pivot_scale = 1. / pivot_val;
565 for (i = k + 1; i < l; i++) {
566 double x = m[i].a[0] * pivot_scale;
567 m[k].al[i - k - 1] = x;
568 for (j = 1; j < 11; j++)
569 m[i].a[j - 1] = m[i].a[j] - x * m[k].a[j];
570 m[i].a[10] = 0.;
571 }
572 }
573 }
575 static void
576 banbks11(const bandmat *m, const int *perm, double *v, int n)
577 {
578 int i, k, l;
580 /* forward substitution */
581 l = 5;
582 for (k = 0; k < n; k++) {
583 i = perm[k];
584 if (i != k) {
585 double tmp = v[k];
586 v[k] = v[i];
587 v[i] = tmp;
588 }
589 if (l < n) l++;
590 for (i = k + 1; i < l; i++)
591 v[i] -= m[k].al[i - k - 1] * v[k];
592 }
594 /* back substitution */
595 l = 1;
596 for (i = n - 1; i >= 0; i--) {
597 double x = v[i];
598 for (k = 1; k < l; k++)
599 x -= m[i].a[k] * v[k + i];
600 v[i] = x / m[i].a[0];
601 if (l < 11) l++;
602 }
603 }
605 int compute_jinc(char ty0, char ty1)
606 {
607 if (ty0 == 'o' || ty1 == 'o' ||
608 ty0 == ']' || ty1 == '[')
609 return 4;
610 else if (ty0 == 'c' && ty1 == 'c')
611 return 2;
612 else if (((ty0 == '{' || ty0 == 'v' || ty0 == '[') && ty1 == 'c') ||
613 (ty0 == 'c' && (ty1 == '}' || ty1 == 'v' || ty1 == ']')))
614 return 1;
615 else
616 return 0;
617 }
619 int count_vec(const spiro_seg *s, int nseg)
620 {
621 int i;
622 int n = 0;
624 for (i = 0; i < nseg; i++)
625 n += compute_jinc(s[i].ty, s[i + 1].ty);
626 return n;
627 }
629 static void
630 add_mat_line(bandmat *m, double *v,
631 double derivs[4], double x, double y, int j, int jj, int jinc,
632 int nmat)
633 {
634 int k;
636 if (jj >= 0) {
637 int joff = (j + 5 - jj + nmat) % nmat;
638 if (nmat < 6) {
639 joff = j + 5 - jj;
640 } else if (nmat == 6) {
641 joff = 2 + (j + 3 - jj + nmat) % nmat;
642 }
643 #ifdef VERBOSE
644 printf("add_mat_line j=%d jj=%d jinc=%d nmat=%d joff=%d\n", j, jj, jinc, nmat, joff);
645 #endif
646 v[jj] += x;
647 for (k = 0; k < jinc; k++)
648 m[jj].a[joff + k] += y * derivs[k];
649 }
650 }
652 static double
653 spiro_iter(spiro_seg *s, bandmat *m, int *perm, double *v, int n)
654 {
655 int cyclic = s[0].ty != '{' && s[0].ty != 'v';
656 int i, j, jj;
657 int nmat = count_vec(s, n);
658 double norm;
659 int n_invert;
661 for (i = 0; i < nmat; i++) {
662 v[i] = 0.;
663 for (j = 0; j < 11; j++)
664 m[i].a[j] = 0.;
665 for (j = 0; j < 5; j++)
666 m[i].al[j] = 0.;
667 }
669 j = 0;
670 if (s[0].ty == 'o')
671 jj = nmat - 2;
672 else if (s[0].ty == 'c')
673 jj = nmat - 1;
674 else
675 jj = 0;
676 for (i = 0; i < n; i++) {
677 char ty0 = s[i].ty;
678 char ty1 = s[i + 1].ty;
679 int jinc = compute_jinc(ty0, ty1);
680 double th = s[i].bend_th;
681 double ends[2][4];
682 double derivs[4][2][4];
683 int jthl = -1, jk0l = -1, jk1l = -1, jk2l = -1;
684 int jthr = -1, jk0r = -1, jk1r = -1, jk2r = -1;
686 compute_pderivs(&s[i], ends, derivs, jinc);
688 /* constraints crossing left */
689 if (ty0 == 'o' || ty0 == 'c' || ty0 == '[' || ty0 == ']') {
690 jthl = jj++;
691 jj %= nmat;
692 jk0l = jj++;
693 }
694 if (ty0 == 'o') {
695 jj %= nmat;
696 jk1l = jj++;
697 jk2l = jj++;
698 }
700 /* constraints on left */
701 if ((ty0 == '[' || ty0 == 'v' || ty0 == '{' || ty0 == 'c') &&
702 jinc == 4) {
703 if (ty0 != 'c')
704 jk1l = jj++;
705 jk2l = jj++;
706 }
708 /* constraints on right */
709 if ((ty1 == ']' || ty1 == 'v' || ty1 == '}' || ty1 == 'c') &&
710 jinc == 4) {
711 if (ty1 != 'c')
712 jk1r = jj++;
713 jk2r = jj++;
714 }
716 /* constraints crossing right */
717 if (ty1 == 'o' || ty1 == 'c' || ty1 == '[' || ty1 == ']') {
718 jthr = jj;
719 jk0r = (jj + 1) % nmat;
720 }
721 if (ty1 == 'o') {
722 jk1r = (jj + 2) % nmat;
723 jk2r = (jj + 3) % nmat;
724 }
726 add_mat_line(m, v, derivs[0][0], th - ends[0][0], 1, j, jthl, jinc, nmat);
727 add_mat_line(m, v, derivs[1][0], ends[0][1], -1, j, jk0l, jinc, nmat);
728 add_mat_line(m, v, derivs[2][0], ends[0][2], -1, j, jk1l, jinc, nmat);
729 add_mat_line(m, v, derivs[3][0], ends[0][3], -1, j, jk2l, jinc, nmat);
730 add_mat_line(m, v, derivs[0][1], -ends[1][0], 1, j, jthr, jinc, nmat);
731 add_mat_line(m, v, derivs[1][1], -ends[1][1], 1, j, jk0r, jinc, nmat);
732 add_mat_line(m, v, derivs[2][1], -ends[1][2], 1, j, jk1r, jinc, nmat);
733 add_mat_line(m, v, derivs[3][1], -ends[1][3], 1, j, jk2r, jinc, nmat);
734 if (jthl >= 0)
735 v[jthl] = mod_2pi(v[jthl]);
736 if (jthr >= 0)
737 v[jthr] = mod_2pi(v[jthr]);
738 j += jinc;
739 }
740 if (cyclic) {
741 memcpy(m + nmat, m, sizeof(bandmat) * nmat);
742 memcpy(m + 2 * nmat, m, sizeof(bandmat) * nmat);
743 memcpy(v + nmat, v, sizeof(double) * nmat);
744 memcpy(v + 2 * nmat, v, sizeof(double) * nmat);
745 n_invert = 3 * nmat;
746 j = nmat;
747 } else {
748 n_invert = nmat;
749 j = 0;
750 }
751 #ifdef VERBOSE
752 for (i = 0; i < n; i++) {
753 int k;
754 for (k = 0; k < 11; k++)
755 printf(" %2.4f", m[i].a[k]);
756 printf(": %2.4f\n", v[i]);
757 }
758 printf("---\n");
759 #endif
760 bandec11(m, perm, n_invert);
761 banbks11(m, perm, v, n_invert);
762 norm = 0.;
763 for (i = 0; i < n; i++) {
764 char ty0 = s[i].ty;
765 char ty1 = s[i + 1].ty;
766 int jinc = compute_jinc(ty0, ty1);
767 int k;
769 for (k = 0; k < jinc; k++) {
770 double dk = v[j++];
772 #ifdef VERBOSE
773 printf("s[%d].ks[%d] += %f\n", i, k, dk);
774 #endif
775 s[i].ks[k] += dk;
776 norm += dk * dk;
777 }
778 s[i].ks[0] = 2.0*mod_2pi(s[i].ks[0]/2.0);
779 }
780 return norm;
781 }
783 int
784 solve_spiro(spiro_seg *s, int nseg)
785 {
786 bandmat *m;
787 double *v;
788 int *perm;
789 int nmat = count_vec(s, nseg);
790 int n_alloc = nmat;
791 double norm;
792 int i;
794 if (nmat == 0)
795 return 0;
796 if (s[0].ty != '{' && s[0].ty != 'v')
797 n_alloc *= 3;
798 if (n_alloc < 5)
799 n_alloc = 5;
800 m = (bandmat *)malloc(sizeof(bandmat) * n_alloc);
801 v = (double *)malloc(sizeof(double) * n_alloc);
802 perm = (int *)malloc(sizeof(int) * n_alloc);
804 for (i = 0; i < 10; i++) {
805 norm = spiro_iter(s, m, perm, v, nseg);
806 #ifdef VERBOSE
807 printf("%% norm = %g\n", norm);
808 #endif
809 if (norm < 1e-12) break;
810 }
812 free(m);
813 free(v);
814 free(perm);
815 return 0;
816 }
818 static void
819 spiro_seg_to_bpath(const double ks[4],
820 double x0, double y0, double x1, double y1,
821 bezctx *bc, int depth)
822 {
823 double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
824 fabs((1./48) * ks[3]);
826 if (!bend > 1e-8) {
827 bezctx_lineto(bc, x1, y1);
828 } else {
829 double seg_ch = hypot(x1 - x0, y1 - y0);
830 double seg_th = atan2(y1 - y0, x1 - x0);
831 double xy[2];
832 double ch, th;
833 double scale, rot;
834 double th_even, th_odd;
835 double ul, vl;
836 double ur, vr;
838 integrate_spiro(ks, xy);
839 ch = hypot(xy[0], xy[1]);
840 th = atan2(xy[1], xy[0]);
841 scale = seg_ch / ch;
842 rot = seg_th - th;
843 if (depth > 5 || bend < 1.) {
844 th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
845 th_odd = (1./48) * ks[2] + .5 * ks[0];
846 ul = (scale * (1./3)) * cos(th_even - th_odd);
847 vl = (scale * (1./3)) * sin(th_even - th_odd);
848 ur = (scale * (1./3)) * cos(th_even + th_odd);
849 vr = (scale * (1./3)) * sin(th_even + th_odd);
850 bezctx_curveto(bc, x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
851 } else {
852 /* subdivide */
853 double ksub[4];
854 double thsub;
855 double xysub[2];
856 double xmid, ymid;
857 double cth, sth;
859 ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
860 ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
861 ksub[2] = .125 * ks[2] - (1./32) * ks[3];
862 ksub[3] = (1./16) * ks[3];
863 thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
864 cth = .5 * scale * cos(thsub);
865 sth = .5 * scale * sin(thsub);
866 integrate_spiro(ksub, xysub);
867 xmid = x0 + cth * xysub[0] - sth * xysub[1];
868 ymid = y0 + cth * xysub[1] + sth * xysub[0];
869 spiro_seg_to_bpath(ksub, x0, y0, xmid, ymid, bc, depth + 1);
870 ksub[0] += .25 * ks[1] + (1./384) * ks[3];
871 ksub[1] += .125 * ks[2];
872 ksub[2] += (1./16) * ks[3];
873 spiro_seg_to_bpath(ksub, xmid, ymid, x1, y1, bc, depth + 1);
874 }
875 }
876 }
878 spiro_seg *
879 run_spiro(const spiro_cp *src, int n)
880 {
881 int nseg = src[0].ty == '{' ? n - 1 : n;
882 spiro_seg *s = setup_path(src, n);
883 if (nseg > 1)
884 solve_spiro(s, nseg);
885 return s;
886 }
888 void
889 free_spiro(spiro_seg *s)
890 {
891 free(s);
892 }
894 void
895 spiro_to_bpath(const spiro_seg *s, int n, bezctx *bc)
896 {
897 int i;
898 int nsegs = s[n - 1].ty == '}' ? n - 1 : n;
900 for (i = 0; i < nsegs; i++) {
901 double x0 = s[i].x;
902 double y0 = s[i].y;
903 double x1 = s[i + 1].x;
904 double y1 = s[i + 1].y;
906 if (i == 0)
907 bezctx_moveto(bc, x0, y0, s[0].ty == '{');
908 bezctx_mark_knot(bc, i);
909 spiro_seg_to_bpath(s[i].ks, x0, y0, x1, y1, bc, 0);
910 }
911 }
913 double
914 get_knot_th(const spiro_seg *s, int i)
915 {
916 double ends[2][4];
918 if (i == 0) {
919 compute_ends(s[i].ks, ends, s[i].seg_ch);
920 return s[i].seg_th - ends[0][0];
921 } else {
922 compute_ends(s[i - 1].ks, ends, s[i - 1].seg_ch);
923 return s[i - 1].seg_th + ends[1][0];
924 }
925 }
927 #ifdef UNIT_TEST
928 #include <stdio.h>
929 #include <sys/time.h> /* for gettimeofday */
931 static double
932 get_time (void)
933 {
934 struct timeval tv;
935 struct timezone tz;
937 gettimeofday (&tv, &tz);
939 return tv.tv_sec + 1e-6 * tv.tv_usec;
940 }
942 int
943 test_integ(void) {
944 double ks[] = {1, 2, 3, 4};
945 double xy[2];
946 double xynom[2];
947 double ch, th;
948 int i, j;
949 int nsubdiv;
951 n = ORDER < 6 ? 4096 : 1024;
952 integrate_spiro(ks, xynom);
953 nsubdiv = ORDER < 12 ? 8 : 7;
954 for (i = 0; i < nsubdiv; i++) {
955 double st, en;
956 double err;
957 int n_iter = (1 << (20 - i));
959 n = 1 << i;
960 st = get_time();
961 for (j = 0; j < n_iter; j++)
962 integrate_spiro(ks, xy);
963 en = get_time();
964 err = hypot(xy[0] - xynom[0], xy[1] - xynom[1]);
965 printf("%d %d %g %g\n", ORDER, n, (en - st) / n_iter, err);
966 ch = hypot(xy[0], xy[1]);
967 th = atan2(xy[1], xy[0]);
968 #if 0
969 printf("n = %d: integ(%g %g %g %g) = %g %g, ch = %g, th = %g\n", n,
970 ks[0], ks[1], ks[2], ks[3], xy[0], xy[1], ch, th);
971 printf("%d: %g %g\n", n, xy[0] - xynom[0], xy[1] - xynom[1]);
972 #endif
973 }
974 return 0;
975 }
977 void
978 print_seg(const double ks[4], double x0, double y0, double x1, double y1)
979 {
980 double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
981 fabs((1./48) * ks[3]);
983 if (bend < 1e-8) {
984 printf("%g %g lineto\n", x1, y1);
985 } else {
986 double seg_ch = hypot(x1 - x0, y1 - y0);
987 double seg_th = atan2(y1 - y0, x1 - x0);
988 double xy[2];
989 double ch, th;
990 double scale, rot;
991 double th_even, th_odd;
992 double ul, vl;
993 double ur, vr;
995 integrate_spiro(ks, xy);
996 ch = hypot(xy[0], xy[1]);
997 th = atan2(xy[1], xy[0]);
998 scale = seg_ch / ch;
999 rot = seg_th - th;
1000 if (bend < 1.) {
1001 th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
1002 th_odd = (1./48) * ks[2] + .5 * ks[0];
1003 ul = (scale * (1./3)) * cos(th_even - th_odd);
1004 vl = (scale * (1./3)) * sin(th_even - th_odd);
1005 ur = (scale * (1./3)) * cos(th_even + th_odd);
1006 vr = (scale * (1./3)) * sin(th_even + th_odd);
1007 printf("%g %g %g %g %g %g curveto\n",
1008 x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
1010 } else {
1011 /* subdivide */
1012 double ksub[4];
1013 double thsub;
1014 double xysub[2];
1015 double xmid, ymid;
1016 double cth, sth;
1018 ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
1019 ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
1020 ksub[2] = .125 * ks[2] - (1./32) * ks[3];
1021 ksub[3] = (1./16) * ks[3];
1022 thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
1023 cth = .5 * scale * cos(thsub);
1024 sth = .5 * scale * sin(thsub);
1025 integrate_spiro(ksub, xysub);
1026 xmid = x0 + cth * xysub[0] - sth * xysub[1];
1027 ymid = y0 + cth * xysub[1] + sth * xysub[0];
1028 print_seg(ksub, x0, y0, xmid, ymid);
1029 ksub[0] += .25 * ks[1] + (1./384) * ks[3];
1030 ksub[1] += .125 * ks[2];
1031 ksub[2] += (1./16) * ks[3];
1032 print_seg(ksub, xmid, ymid, x1, y1);
1033 }
1034 }
1035 }
1037 void
1038 print_segs(const spiro_seg *segs, int nsegs)
1039 {
1040 int i;
1042 for (i = 0; i < nsegs; i++) {
1043 double x0 = segs[i].x;
1044 double y0 = segs[i].y;
1045 double x1 = segs[i + 1].x;
1046 double y1 = segs[i + 1].y;
1048 if (i == 0)
1049 printf("%g %g moveto\n", x0, y0);
1050 printf("%% ks = [ %g %g %g %g ]\n",
1051 segs[i].ks[0], segs[i].ks[1], segs[i].ks[2], segs[i].ks[3]);
1052 print_seg(segs[i].ks, x0, y0, x1, y1);
1053 }
1054 printf("stroke\n");
1055 }
1057 int
1058 test_curve(void)
1059 {
1060 spiro_cp path[] = {
1061 {334, 117, 'v'},
1062 {305, 176, 'v'},
1063 {212, 142, 'c'},
1064 {159, 171, 'c'},
1065 {224, 237, 'c'},
1066 {347, 335, 'c'},
1067 {202, 467, 'c'},
1068 {81, 429, 'v'},
1069 {114, 368, 'v'},
1070 {201, 402, 'c'},
1071 {276, 369, 'c'},
1072 {218, 308, 'c'},
1073 {91, 211, 'c'},
1074 {124, 111, 'c'},
1075 {229, 82, 'c'}
1076 };
1077 spiro_seg *segs;
1078 int i;
1080 n = 1;
1081 for (i = 0; i < 1000; i++) {
1082 segs = setup_path(path, 15);
1083 solve_spiro(segs, 15);
1084 }
1085 printf("100 800 translate 1 -1 scale 1 setlinewidth\n");
1086 print_segs(segs, 15);
1087 printf("showpage\n");
1088 return 0;
1089 }
1091 int main(int argc, char **argv)
1092 {
1093 return test_curve();
1094 }
1095 #endif