1/*****
2 * runpath3.in
3 *
4 * Runtime functions for path3 operations.
5 *
6 *****/
7
8pair     => primPair()
9triple   => primTriple()
10path3     => primPath3()
11boolarray* => booleanArray()
12realarray* => realArray()
13realarray2* => realArray2()
14triplearray* => tripleArray()
15triplearray2* => tripleArray2()
16
17#include "path3.h"
18#include "array.h"
19#include "drawsurface.h"
20#include "predicates.h"
21
22using namespace camp;
23using namespace vm;
24
25typedef array boolarray;
26typedef array realarray;
27typedef array realarray2;
28typedef array triplearray;
29typedef array triplearray2;
30
31using types::booleanArray;
32using types::realArray;
33using types::realArray2;
34using types::tripleArray;
35using types::tripleArray2;
36
37// Autogenerated routines:
38
39
40path3 path3(triplearray *pre, triplearray *point, triplearray *post,
41            boolarray *straight, bool cyclic)
42{
43  size_t n=checkArrays(pre,point);
44  checkEqual(n,checkArray(post));
45  checkEqual(n,checkArray(straight));
46  mem::vector<solvedKnot3> nodes(n);
47  for(size_t i=0; i < n; ++i) {
48    nodes[i].pre=read<triple>(pre,i);
49    nodes[i].point=read<triple>(point,i);
50    nodes[i].post=read<triple>(post,i);
51    nodes[i].straight=read<bool>(straight,i);
52  }
53
54  return path3(nodes,(Int) n,cyclic);
55}
56
57path3 :nullPath3()
58{
59  return nullpath3;
60}
61
62bool ==(path3 a, path3 b)
63{
64  return a == b;
65}
66
67bool !=(path3 a, path3 b)
68{
69  return !(a == b);
70}
71
72triple point(path3 p, Int t)
73{
74  return p.point((Int) t);
75}
76
77triple point(path3 p, real t)
78{
79  return p.point(t);
80}
81
82triple precontrol(path3 p, Int t)
83{
84  return p.precontrol((Int) t);
85}
86
87triple precontrol(path3 p, real t)
88{
89  return p.precontrol(t);
90}
91
92triple postcontrol(path3 p, Int t)
93{
94  return p.postcontrol((Int) t);
95}
96
97triple postcontrol(path3 p, real t)
98{
99  return p.postcontrol(t);
100}
101
102triple dir(path3 p, Int t, Int sign=0, bool normalize=true)
103{
104  return p.dir(t,sign,normalize);
105}
106
107triple dir(path3 p, real t, bool normalize=true)
108{
109  return p.dir(t,normalize);
110}
111
112triple accel(path3 p, Int t, Int sign=0)
113{
114  return p.accel(t,sign);
115}
116
117triple accel(path3 p, real t)
118{
119  return p.accel(t);
120}
121
122real radius(path3 p, real t)
123{
124  triple v=p.dir(t,false);
125  triple a=p.accel(t);
126  real d=dot(a,v);
127  real v2=v.abs2();
128  real a2=a.abs2();
129  real denom=v2*a2-d*d;
130  real r=v2*sqrt(v2);
131  return denom > 0 ? r/sqrt(denom) : 0.0;
132}
133
134real radius(triple z0, triple c0, triple c1, triple z1, real t)
135{
136  triple v=(3.0*(z1-z0)+9.0*(c0-c1))*t*t+(6.0*(z0+c1)-12.0*c0)*t+3.0*(c0-z0);
137  triple a=6.0*(z1-z0+3.0*(c0-c1))*t+6.0*(z0+c1)-12.0*c0;
138  real d=dot(a,v);
139  real v2=v.abs2();
140  real a2=a.abs2();
141  real denom=v2*a2-d*d;
142  real r=v2*sqrt(v2);
143  return denom > 0 ? r/sqrt(denom) : 0.0;
144}
145
146path3 reverse(path3 p)
147{
148  return p.reverse();
149}
150
151path3 subpath(path3 p, Int a, Int b)
152{
153  return p.subpath((Int) a, (Int) b);
154}
155
156path3 subpath(path3 p, real a, real b)
157{
158  return p.subpath(a,b);
159}
160
161Int length(path3 p)
162{
163  return p.length();
164}
165
166bool cyclic(path3 p)
167{
168  return p.cyclic();
169}
170
171bool straight(path3 p, Int t)
172{
173  return p.straight(t);
174}
175
176path3 unstraighten(path3 p)
177{
178  return p.unstraighten();
179}
180
181// Return the maximum perpendicular deviation of segment i of path3 g
182// from a straight line.
183real straightness(path3 p, Int t)
184{
185  if(p.straight(t)) return 0;
186  triple z0=p.point(t);
187  triple u=unit(p.point(t+1)-z0);
188  return ::max(length(perp(p.postcontrol(t)-z0,u)),
189               length(perp(p.precontrol(t+1)-z0,u)));
190}
191
192// Return the maximum perpendicular deviation of z0..controls c0 and c1..z1
193// from a straight line.
194real straightness(triple z0, triple c0, triple c1, triple z1)
195{
196  triple u=unit(z1-z0);
197  return ::max(length(perp(c0-z0,u)),length(perp(c1-z0,u)));
198}
199
200bool piecewisestraight(path3 p)
201{
202  return p.piecewisestraight();
203}
204
205real arclength(path3 p)
206{
207  return p.arclength();
208}
209
210real arctime(path3 p, real dval)
211{
212  return p.arctime(dval);
213}
214
215realarray* intersect(path3 p, path3 q, real fuzz=-1)
216{
217  bool exact=fuzz <= 0.0;
218  if(fuzz < 0)
219    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
220                       ::max(length(q.max()),length(q.min())));
221
222  std::vector<real> S,T;
223  real s,t;
224  if(intersections(s,t,S,T,p,q,fuzz,true,exact)) {
225    array *V=new array(2);
226    (*V)[0]=s;
227    (*V)[1]=t;
228    return V;
229  } else
230    return new array(0);
231}
232
233realarray2* intersections(path3 p, path3 q, real fuzz=-1)
234{
235  bool exact=fuzz <= 0.0;
236  if(fuzz < 0)
237    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
238                       ::max(length(q.max()),length(q.min())));
239  bool single=!exact;
240
241  real s,t;
242  std::vector<real> S,T;
243  bool found=intersections(s,t,S,T,p,q,fuzz,single,exact);
244  if(!found) return new array(0);
245  array *V;
246  if(single) {
247    V=new array(1);
248    array *Vi=new array(2);
249    (*V)[0]=Vi;
250    (*Vi)[0]=s;
251    (*Vi)[1]=t;
252  } else {
253    size_t n=S.size();
254    V=new array(n);
255    for(size_t i=0; i < n; ++i) {
256      array *Vi=new array(2);
257      (*V)[i]=Vi;
258      (*Vi)[0]=S[i];
259      (*Vi)[1]=T[i];
260    }
261  }
262  stable_sort(V->begin(),V->end(),run::compare2<real>());
263  return V;
264}
265
266realarray* intersect(path3 p, triplearray2 *P, real fuzz=-1)
267{
268  triple *A;
269  copyArray2C(A,P,true,4);
270  if(fuzz <= 0) fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
271                                   norm(A,16));
272  std::vector<real> T,U,V;
273  bool found=intersections(T,U,V,p,A,fuzz,true);
274  delete[] A;
275  if(found) {
276    array *W=new array(3);
277    (*W)[0]=T[0];
278    (*W)[1]=U[0];
279    (*W)[2]=V[0];
280    return W;
281  } else
282    return new array(0);
283}
284
285realarray2* intersections(path3 p, triplearray2 *P, real fuzz=-1)
286{
287  triple *A;
288  copyArray2C(A,P,true,4);
289  if(fuzz <= 0) fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
290                                   norm(A,16));
291  std::vector<real> T,U,V;
292  intersections(T,U,V,p,A,fuzz,false);
293  delete[] A;
294  size_t n=T.size();
295  array *W=new array(n);
296  for(size_t i=0; i < n; ++i) {
297    array *Wi=new array(3);
298    (*W)[i]=Wi;
299    (*Wi)[0]=T[i];
300    (*Wi)[1]=U[i];
301    (*Wi)[2]=V[i];
302  }
303  return W; // Sorting will done in asy.
304}
305
306Int size(path3 p)
307{
308  return p.size();
309}
310
311path3 &(path3 p, path3 q)
312{
313  return camp::concat(p,q);
314}
315
316triple min(path3 p)
317{
318  return p.min();
319}
320
321triple max(path3 p)
322{
323  return p.max();
324}
325
326realarray *mintimes(path3 p)
327{
328  array *V=new array(3);
329  triple v=p.mintimes();
330  (*V)[0]=v.getx();
331  (*V)[1]=v.gety();
332  (*V)[2]=v.getz();
333  return V;
334}
335
336realarray *maxtimes(path3 p)
337{
338  array *V=new array(3);
339  triple v=p.maxtimes();
340  (*V)[0]=v.getx();
341  (*V)[1]=v.gety();
342  (*V)[2]=v.getz();
343  return V;
344}
345
346path3 Operator *(realarray2 *t, path3 g)
347{
348  return transformed(*t,g);
349}
350
351pair minratio(path3 g)
352{
353  return g.ratio(::min);
354}
355
356pair maxratio(path3 g)
357{
358  return g.ratio(::max);
359}
360
361// Return a negative (positive) value if a--b--c--cycle is oriented
362// counterclockwise (clockwise) when viewed from d or zero if all four
363// points are coplanar.
364// The value returned is the determinant
365// |a.x a.y a.z 1|
366// |b.x b.y b.z 1|
367// |c.x c.y c.z 1|
368// |d.x d.y d.z 1|
369real orient(triple a, triple b, triple c, triple d)
370{
371  real A[]={a.getx(),a.gety(),a.getz()};
372  real B[]={b.getx(),b.gety(),b.getz()};
373  real C[]={c.getx(),c.gety(),c.getz()};
374  real D[]={d.getx(),d.gety(),d.getz()};
375  return orient3d(A,B,C,D);
376}
377
378// Return a positive (negative) value if e lies inside (outside)
379// the sphere passing through the points a,b,c,d oriented so that
380// a--b--c--cycle appears in clockwise order when viewed from d
381// or zero if all five points are cospherical.
382// The value returned is the determinant
383// |a.x a.y a.z a.x^2+a.y^2+a.z^2 1|
384// |b.x b.y b.z b.x^2+b.y^2+b.z^2 1|
385// |c.x c.y c.z c.x^2+c.y^2+c.z^2 1|
386// |d.x d.y d.z d.x^2+d.y^2+d.z^2 1|
387// |e.x e.y e.z e.x^2+e.y^2+e.z^2 1|
388real insphere(triple a, triple b, triple c, triple d, triple e)
389{
390  real A[]={a.getx(),a.gety(),a.getz()};
391  real B[]={b.getx(),b.gety(),b.getz()};
392  real C[]={c.getx(),c.gety(),c.getz()};
393  real D[]={d.getx(),d.gety(),d.getz()};
394  real E[]={e.getx(),e.gety(),e.getz()};
395  return insphere(A,B,C,D,E);
396}
397