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