1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
|
path nullpath;
typedef guide interpolate(... guide[]);
// These numbers identify the side of a specifier in an operator spec or
// operator curl expression:
// a{out} .. {in}b
restricted int JOIN_OUT=0;
restricted int JOIN_IN=1;
// Define a.. tension t ..b to be equivalent to
// a.. tension t and t ..b
// and likewise with controls.
tensionSpecifier operator tension(real t, bool atLeast)
{
return operator tension(t,t,atLeast);
}
guide operator controls(pair z)
{
return operator controls(z,z);
}
guide[] operator cast(pair[] z)
{
return sequence(new guide(int i) {return z[i];},z.length);
}
path[] operator cast(pair[] z)
{
return sequence(new path(int i) {return z[i];},z.length);
}
path[] operator cast(guide[] g)
{
return sequence(new path(int i) {return g[i];},g.length);
}
guide[] operator cast(path[] g)
{
return sequence(new guide(int i) {return g[i];},g.length);
}
path[] operator cast(path p)
{
return new path[] {p};
}
path[] operator cast(guide g)
{
return new path[] {(path) g};
}
path[] operator ^^ (path p, path q)
{
return new path[] {p,q};
}
path[] operator ^^ (path p, explicit path[] q)
{
return concat(new path[] {p},q);
}
path[] operator ^^ (explicit path[] p, path q)
{
return concat(p,new path[] {q});
}
path[] operator ^^ (explicit path[] p, explicit path[] q)
{
return concat(p,q);
}
path[] operator * (transform t, explicit path[] p)
{
return sequence(new path(int i) {return t*p[i];},p.length);
}
pair[] operator * (transform t, pair[] z)
{
return sequence(new pair(int i) {return t*z[i];},z.length);
}
void write(file file, string s="", explicit path[] x, suffix suffix=none)
{
write(file,s);
if(x.length > 0) write(file,x[0]);
for(int i=1; i < x.length; ++i) {
write(file,endl);
write(file," ^^");
write(file,x[i]);
}
write(file,suffix);
}
void write(string s="", explicit path[] x, suffix suffix=endl)
{
write(stdout,s,x,suffix);
}
void write(file file, string s="", explicit guide[] x, suffix suffix=none)
{
write(file,s);
if(x.length > 0) write(file,x[0]);
for(int i=1; i < x.length; ++i) {
write(file,endl);
write(file," ^^");
write(file,x[i]);
}
write(file,suffix);
}
void write(string s="", explicit guide[] x, suffix suffix=endl)
{
write(stdout,s,x,suffix);
}
interpolate operator ..(tensionSpecifier t)
{
return new guide(... guide[] a) {
if(a.length == 0) return nullpath;
guide g=a[0];
for(int i=1; i < a.length; ++i)
g=g..t..a[i];
return g;
};
}
interpolate operator ::=operator ..(operator tension(1,true));
interpolate operator ---=operator ..(operator tension(infinity,true));
// return an arbitrary intersection point of paths p and q
pair intersectionpoint(path p, path q, real fuzz=-1)
{
real[] t=intersect(p,q,fuzz);
if(t.length == 0) abort("paths do not intersect");
return point(p,t[0]);
}
// return an array containing all intersection points of the paths p and q
pair[] intersectionpoints(path p, path q, real fuzz=-1)
{
real[][] t=intersections(p,q,fuzz);
return sequence(new pair(int i) {return point(p,t[i][0]);},t.length);
}
pair[] intersectionpoints(explicit path[] p, explicit path[] q, real fuzz=-1)
{
pair[] z;
for(int i=0; i < p.length; ++i)
for(int j=0; j < q.length; ++j)
z.append(intersectionpoints(p[i],q[j],fuzz));
return z;
}
struct slice {
path before,after;
}
slice cut(path p, path knife, int n)
{
slice s;
real[][] T=intersections(p,knife);
if(T.length == 0) {s.before=p; s.after=nullpath; return s;}
T.cyclic=true;
real t=T[n][0];
s.before=subpath(p,0,t);
s.after=subpath(p,t,length(p));
return s;
}
slice firstcut(path p, path knife)
{
return cut(p,knife,0);
}
slice lastcut(path p, path knife)
{
return cut(p,knife,-1);
}
pair dir(path p)
{
return dir(p,length(p));
}
pair dir(path p, path q)
{
return unit(dir(p)+dir(q));
}
// return the point on path p at arclength L
pair arcpoint(path p, real L)
{
return point(p,arctime(p,L));
}
// return the direction on path p at arclength L
pair arcdir(path p, real L)
{
return dir(p,arctime(p,L));
}
// return the time on path p at the relative fraction l of its arclength
real reltime(path p, real l)
{
return arctime(p,l*arclength(p));
}
// return the point on path p at the relative fraction l of its arclength
pair relpoint(path p, real l)
{
return point(p,reltime(p,l));
}
// return the direction of path p at the relative fraction l of its arclength
pair reldir(path p, real l)
{
return dir(p,reltime(p,l));
}
// return the initial point of path p
pair beginpoint(path p)
{
return point(p,0);
}
// return the point on path p at half of its arclength
pair midpoint(path p)
{
return relpoint(p,0.5);
}
// return the final point of path p
pair endpoint(path p)
{
return point(p,length(p));
}
path operator &(path p, cycleToken tok)
{
int n=length(p);
if(n < 0) return nullpath;
if(n == 0) return p--cycle;
if(cyclic(p)) return p;
return straight(p,n-1) ? subpath(p,0,n-1)--cycle :
subpath(p,0,n-1)..controls postcontrol(p,n-1) and precontrol(p,n)..cycle;
}
// return a cyclic path enclosing a region bounded by a list of two or more
// consecutively intersecting paths
path buildcycle(... path[] p)
{
int n=p.length;
if(n < 2) return nullpath;
real[] ta=new real[n];
real[] tb=new real[n];
if(n == 2) {
real[][] t=intersections(p[0],p[1]);
if(t.length < 2)
return nullpath;
int k=t.length-1;
ta[0]=t[0][0]; tb[0]=t[k][0];
ta[1]=t[k][1]; tb[1]=t[0][1];
} else {
int j=n-1;
for(int i=0; i < n; ++i) {
real[][] t=intersections(p[i],p[j]);
if(t.length == 0)
return nullpath;
ta[i]=t[0][0]; tb[j]=t[0][1];
j=i;
}
}
pair c;
for(int i=0; i < n ; ++i)
c += point(p[i],ta[i]);
c /= n;
path G;
for(int i=0; i < n ; ++i) {
real Ta=ta[i];
real Tb=tb[i];
if(cyclic(p[i])) {
int L=length(p[i]);
real t=Tb-L;
if(abs(c-point(p[i],0.5(Ta+t))) <
abs(c-point(p[i],0.5(Ta+Tb)))) Tb=t;
while(Tb < Ta) Tb += L;
}
G=G&subpath(p[i],Ta,Tb);
}
return G&cycle;
}
// return 1 if p strictly contains q,
// -1 if q strictly contains p,
// 0 otherwise.
int inside(path p, path q, pen fillrule=currentpen)
{
if(intersect(p,q).length > 0) return 0;
if(cyclic(p) && inside(p,point(q,0),fillrule)) return 1;
if(cyclic(q) && inside(q,point(p,0),fillrule)) return -1;
return 0;
}
// Return an arbitrary point strictly inside a cyclic path p according to
// the specified fill rule.
pair inside(path p, pen fillrule=currentpen)
{
if(!cyclic(p)) abort("path is not cyclic");
int n=length(p);
for(int i=0; i < n; ++i) {
pair z=point(p,i);
pair dir=dir(p,i);
if(dir == 0) continue;
real[] T=intersections(p,z,z+I*dir);
// Check midpoints of line segments formed between the
// corresponding intersection points and z.
for(int j=0; j < T.length; ++j) {
if(T[j] != i) {
pair w=point(p,T[j]);
pair m=0.5*(z+w);
if(interior(windingnumber(p,m),fillrule)) return m;
}
}
}
// cannot find an interior point: path is degenerate
return point(p,0);
}
// Return all intersection times of path g with the vertical line through (x,0).
real[] times(path p, real x, real fuzz=-1)
{
return intersections(p,(x,0),(x,1),fuzz);
}
// Return all intersection times of path g with the horizontal line through
// (0,z.y).
real[] times(path p, explicit pair z, real fuzz=-1)
{
return intersections(p,(0,z.y),(1,z.y),fuzz);
}
path randompath(int n, bool cumulate=true, interpolate join=operator ..)
{
guide g;
pair w;
for(int i=0; i <= n; ++i) {
pair z=(unitrand()-0.5,unitrand()-0.5);
if(cumulate) w += z;
else w=z;
g=join(g,w);
}
return g;
}
path[] strokepath(path g, pen p=currentpen)
{
path[] G=_strokepath(g,p);
if(G.length == 0) return G;
pair center(path g) {return 0.5*(min(g)+max(g));}
pair center(path[] g) {return 0.5*(min(g)+max(g));}
return shift(center(g)-center(G))*G;
}
real braceinnerangle=radians(60);
real braceouterangle=radians(70);
real bracemidangle=radians(0);
real bracedefaultratio=0.14;
path brace(pair a, pair b, real amplitude=bracedefaultratio*length(b-a))
{
real length=length(b-a);
real sign=sgn(amplitude);
real hamplitude=0.5*amplitude;
real hlength=0.5*length;
path brace;
if(abs(amplitude) < bracedefaultratio*length) {
real slope=2*bracedefaultratio;
real controldist=(abs(hamplitude))/slope;
brace=(0,0){expi(sign*braceouterangle)}::
{expi(sign*bracemidangle)}(controldist,hamplitude)::
{expi(sign*bracemidangle)}(hlength-controldist,hamplitude)::
{expi(sign*braceinnerangle)}(hlength,amplitude) {expi(-sign*braceinnerangle)}::
{expi(-sign*bracemidangle)}(hlength+controldist,hamplitude)::
{expi(-sign*bracemidangle)}(length-controldist,hamplitude)::
{expi(-sign*braceouterangle)}(length,0);
} else {
brace=(0,0){expi(sign*braceouterangle)}::
{expi(sign*bracemidangle)}(0.25*length,hamplitude)::
{expi(sign*braceinnerangle)}(hlength,amplitude){expi(-sign*braceinnerangle)}::
{expi(-sign*bracemidangle)}(0.75*length,hamplitude)::
{expi(-sign*braceouterangle)}(length,0);
}
return shift(a)*rotate(degrees(b-a,warn=false))*brace;
}
|