//************************************************************* //**** 2D GEOMETRY CLASSES AND UTILITIES, Jarek Rossignac ***** //**** LAST EDITED ON: March 06 2008 ***** //************************************************************* //************************************************************************ //**** POINTS //************************************************************************ // create pt P(float x, float y) {return new pt(x,y); }; // make point (x,y) pt P() {return P(0,0); }; // make point (0,0) pt P(pt P) {return P(P.x,P.y); }; // make copy of point P // combine pt S(pt A, pt B) {return new pt(A.x+B.x,A.y+B.y); }; // Weighted sum: A+B pt S(pt A, pt B, pt C) {return S(A,S(B,C)); }; // Weighted sum: A+B+C pt S(pt A, pt B, pt C, pt D) {return S(S(A,B),S(C,D)); }; // Weighted sum: A+B+C+D pt S(float s, pt A) {return new pt(s*A.x,s*A.y); }; // Weighted sum: sA pt S(pt A, float s, pt B) {return S(A,S(s,B)); }; // Weighted sum: A+sB pt S(float a, pt A, float b, pt B) {return S(S(a,A),S(b,B));} // Weighted sum: aA+bB pt S(float a, pt A, float b, pt B, float c, pt C) {return S(S(a,A),S(b,B),S(c,C));} // Weighted sum: aA+bB+cC pt S(float a, pt A, float b, pt B, float c, pt C, float d, pt D){return A(S(a,A,b,B),S(c,C,d,D));} // Weighted sum: aA+bB+cC+dD pt A(pt A, pt B) {return P((A.x+B.x)/2.0,(A.y+B.y)/2.0); }; // Average: (A+B)/2 pt A(pt A, pt B, pt C) {return P((A.x+B.x+C.x)/3.0,(A.y+B.y+C.y)/3.0); }; // Average: (A+B+C)/3 pt L(pt A, float s, pt B) {return P(A.x+s*(B.x-A.x),A.y+s*(B.y-A.y)); }; // Linear interpolation: A+sAB // transform pt R(pt Q, float a) {float dx=Q.x, dy=Q.y, c=cos(a), s=sin(a); return new pt(c*dx+s*dy, -s*dx+c*dy); }; // Rotated Q by a around origin pt R(pt Q, float a, pt P) {float dx=Q.x-P.x, dy=Q.y-P.y, c=cos(a), s=sin(a); return P(P.x+c*dx-s*dy, P.y+s*dx+c*dy); }; // Rotated Q by a around P pt T(pt P, vec V) {return P(P.x + V.x, P.y + V.y); } // Translate: P+V pt T(pt P, float s, vec V) {return T(P,S(s,V)); } // Translate: P+sV pt T(pt A, float s, pt B) { return T(A,s,U(V(A,B))); }; // Translate: P by s towards Q: P+sU(PQ) // measure boolean isSame(pt A, pt B) {return (A.x==B.x)&&(A.y==B.y) ;} // Equality: A==B boolean isSame(pt A, pt B, float e) {return ((abs(A.x-B.x)0)&&(mouseX0)&&(mouseYPI) return mPItoPIangle(a-2*PI); if(a<-PI) return mPItoPIangle(a+2*PI); return a; }; float angle(vec V) {return(atan2(V.y,V.x)); }; float angle(pt A, pt B, pt C) {return angle(V(B,A),V(B,C)); } // angle (A,B,C) float toDeg(float a) {return(a*180/PI);} float toRad(float a) {return(a*PI/180);} float mPItoPIangle(float a) { if(a>PI) return(mPItoPIangle(a-2*PI)); if(a<-PI) return(mPItoPIangle(a+2*PI)); return(a);}; float Oto2PIangle(float a) { if(a>TWO_PI) return(Oto2PIangle(a-2*PI)); if(a<0) return(Oto2PIangle(a+2*PI)); return(a);}; //************************************************************************ //**** EDGES //************************************************************************ //************************************************************************ //**** RAYS //************************************************************************ boolean isRightOf(pt A, pt Q, vec T) {return dot(R(T),V(Q,A)) > 0 ; }; // A is on right of ray(Q,T) (as seen on screen) //************************************************************************ //**** TRIANGLES //************************************************************************ // centers pt CenterMass (pt A, pt B, pt C) {return A(A,B,C);} pt CenterCircum (pt A, pt B, pt C) { vec AB = V(A,B); vec AC = R(V(A,C)); return T(A,1./2/dot(AB,AC),S(-n2(AC),R(AB),n2(AB),AC)); }; // circumcenter of tri(A,B,C) pt CenterIn (pt A, pt B, pt C) { float Z=area(A,B,C), a=B.disTo(C), b=C.disTo(A), c=A.disTo(B), s=a+b+c, r=2*Z/s, R=a*b*c/(2*r*s); return S(a/s,A,b/s,B,c/s,C); } pt CenterD (pt A, pt B, pt C) { float Z=area(A,B,C), a=B.disTo(C), b=C.disTo(A), c=A.disTo(B), ss=a*a+b*b+c*c; return S(a*a/ss,A,b*b/ss,B,c*c/ss,C); } pt CenterN (pt A, pt B, pt C) { float Z=area(A,B,C), a=B.disTo(C), b=C.disTo(A), c=A.disTo(B), s=a+b+c; return S((s-a)/2/s,A,(s-b)/2/s,B,(s-c)/2/s,C); } // measures float radiusCircum (pt A, pt B, pt C) { float a=B.disTo(C); float b=C.disTo(A); float c=A.disTo(B); float s=(a+b+c)/2; float d=sqrt(s*(s-a)*(s-b)*(s-c)); float r=a*b*c/4/d; return (r); }; float radiusMonotonic (pt A, pt B, pt C) { // size of bubble pushed through (A,C) and touching B, >0 when ABC ic clockwise float a=d(B,C), b=d(C,A), c=d(A,B); float s=(a+b+c)/2; float d=sqrt(s*(s-a)*(s-b)*(s-c)); float r=a*b*c/4/d; if (abs(angle(A,B,C))>PI/2) r=sq(d(A,C)/2)/r; if (abs(angle(C,A,B))>PI/2) r=sq(d(C,B)/2)/r; if (abs(angle(B,C,A))>PI/2) r=sq(d(B,A)/2)/r; if (!isRightTurn(A,B,C)) r=-r; return r; }; float radiusN (pt A, pt B, pt C) { float Z=area(A,B,C), a=B.disTo(C), b=C.disTo(A), c=A.disTo(B), s=a+b+c; return 2*Z/s; } float area(pt A, pt B, pt C) {return 0.5*dot(V(A,C),R(V(A,B))); } // signed area of triangle float thickness(pt A, pt B, pt C) {float a = abs(dot(V(B,A),U(R(V(B,C))))), b = abs(dot(V(C,B),U(R(V(C,A))))), c = abs(dot(V(A,C),U(R(V(A,B))))); return min (a,b,c); } boolean ccw(pt A, pt B, pt C) {return C.isLeftOf(A,B) ;} // render void show(pt A, pt B, pt C) {beginShape(); A.v(); B.v(); C.v(); endShape(CLOSE);} void show(pt A, pt B, pt C, float r) {if (thickness(A,B,C)<2*r) return; float s=r; if (!ccw(A,B,C)) s=-s; pt AA = A.makeOffset(B,C,s); pt BB = B.makeOffset(C,A,s); pt CC = C.makeOffset(A,B,s); beginShape(); AA.v(); BB.v(); CC.v(); endShape(CLOSE); } // constructions pt makeProjectionOnLine(pt P, pt B, pt Q) {float a=dot(V(P,B),V(P,Q)), b=d2(P,Q); return T(P,a/b,Q); }; pt makeOffset(pt A, pt B, pt C, float r) { float a = angle(V(B,A),V(B,C))/2; float d = r/sin(a); vec N = U(A(U(V(B,A)),U(V(B,C)))); return T(B,d,N); }; boolean isRightTurn(pt A, pt B, pt C) {return dot( R(V(A,B)),V(B,C)) > 0 ; }; // A-B-C makes right turn (as seen on screen) //************************************************************************ //**** CIRCLES //************************************************************************ // ????????????? pt circumCenterOLD (pt A, pt B, pt C) { // computes the center of a circumscirbing circle to triangle (A,B,C) vec AB = A.makeVecTo(B); float ab2 = dot(AB,AB); vec AC = A.makeVecTo(C); AC.turnLeft(); float ac2 = dot(AC,AC); float d = 2*dot(AB,AC); AB.turnLeft(); AB.scaleBy(-ac2); AC.scaleBy(ab2); AB.add(AC); AB.scaleBy(1./d); pt X = A.makeClone(); X.translateBy(AB); return(X); }; //************************************************************************ //**** ARCS //************************************************************************ //************************************************************************ //**** CURVES //************************************************************************ void retrofitBezier(pt[] PP, pt[] QQ) { // sets control polygon QQ so that tits Bezier curve interpolates PP QQ[0]=P(PP[0]); QQ[1]=S(-15./18.,PP[0],54./18.,PP[1],-27./18.,PP[2],6./18.,PP[3]); QQ[2]=S(-15./18.,PP[3],54./18.,PP[2],-27./18.,PP[1],6./18.,PP[0]); QQ[3]=P(PP[3]); } //************************************************************************ //**** POINTS //************************************************************************ class pt { float x=0,y=0; // CREATE pt () {} pt (float px, float py) {x = px; y = py;}; pt (pt P) {x = P.x; y = P.y;}; pt (pt P, vec V) {x = P.x+V.x; y = P.y+V.y;}; pt (pt P, float s, vec V) {x = P.x+s*V.x; y = P.y+s*V.y;}; pt (pt A, float s, pt B) {x = A.x+s*(B.x-A.x); y = A.y+s*(B.y-A.y);}; // MODIFY void setTo(float px, float py) {x = px; y = py;}; void setTo(pt P) {x = P.x; y = P.y;}; void setToMouse() { x = mouseX; y = mouseY; }; void moveWithMouse() { x += mouseX-pmouseX; y += mouseY-pmouseY; }; void translateToTrack(float s, pt P) {setTo(T(P,s,V(P,this)));}; void track(float s, pt P) {setTo(T(P,s,V(P,this)));}; void scaleBy(float f) {x*=f; y*=f;}; void scaleBy(float u, float v) {x*=u; y*=v;}; void translateBy(vec V) {x += V.x; y += V.y;}; void addVec(vec V) {x += V.x; y += V.y;}; void translateBy(float s, vec V) {x += s*V.x; y += s*V.y;}; void addScaled(float s, vec V) {x += s*V.x; y += s*V.y;}; void addScaled(float s, pt P) {x += s*P.x; y += s*P.y;}; void addScaledPt(float s, pt P) {x += s*P.x; y += s*P.y;}; // incorrect notation, but useful for computing weighted averages void translateBy(float u, float v) {x += u; y += v;}; void translateTowards(float s, pt P) {x+=s*(P.x-x); y+=s*(P.y-y); }; void translateTowardsBy(float s, pt P) {vec V = this.makeVecTo(P); V.normalize(); this.translateBy(s,V); }; void addPt(pt P) {x += P.x; y += P.y;}; // incorrect notation, but useful for computing weighted averages void rotateBy(float a) {float dx=x, dy=y, c=cos(a), s=sin(a); x=c*dx+s*dy; y=-s*dx+c*dy; }; // around origin void rotateBy(float a, pt P) {float dx=x-P.x, dy=y-P.y, c=cos(a), s=sin(a); x=P.x+c*dx+s*dy; y=P.y-s*dx+c*dy; }; // around point P void rotateBy(float s, float t, pt P) {float dx=x-P.x, dy=y-P.y; dx-=dy*t; dy+=dx*s; dx-=dy*t; x=P.x+dx; y=P.y+dy; }; // s=sin(a); t=tan(a/2); void clipToWindow() {x=max(x,0); y=max(y,0); x=min(x,height); y=min(y,height); } // OUTPUT POINT pt clone() {return new pt(x,y); }; pt make() {return new pt(x,y); }; pt makeClone() {return new pt(x,y); }; pt makeTranslatedBy(vec V) {return(new pt(x + V.x, y + V.y));}; pt makeTranslatedBy(float s, vec V) {return(new pt(x + s*V.x, y + s*V.y));}; pt makeTransaltedTowards(float s, pt P) {return(new pt(x + s*(P.x-x), y + s*(P.y-y)));}; pt makeTranslatedBy(float u, float v) {return(new pt(x + u, y + v));}; pt makeRotatedBy(float a, pt P) {float dx=x-P.x, dy=y-P.y, c=cos(a), s=sin(a); return(new pt(P.x+c*dx+s*dy, P.y-s*dx+c*dy)); }; pt makeRotatedBy(float a) {float dx=x, dy=y, c=cos(a), s=sin(a); return(new pt(c*dx+s*dy, -s*dx+c*dy)); }; pt makeProjectionOnLine(pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q)), b=dot(P.makeVecTo(Q),P.makeVecTo(Q)); return(P.makeTransaltedTowards(a/b,Q)); }; pt makeOffset(pt P, pt Q, float r) { float a = angle(vecTo(P),vecTo(Q))/2; float h = r/tan(a); vec T = vecTo(P); T.normalize(); vec N = T.left(); pt R = new pt(x,y); R.translateBy(h,T); R.translateBy(r,N); return R; }; // OUTPUT VEC vec vecTo(pt P) {return(new vec(P.x-x,P.y-y)); }; vec makeVecTo(pt P) {return(new vec(P.x-x,P.y-y)); }; vec makeVecToCenter () {return(new vec(x-height/2.,y-height/2.)); }; vec makeVecToAverage (pt P, pt Q) {return(new vec((P.x+Q.x)/2.0-x,(P.y+Q.y)/2.0-y)); }; vec makeVecToAverage (pt P, pt Q, pt R) {return(new vec((P.x+Q.x+R.x)/3.0-x,(P.y+Q.y+R.x)/3.0-y)); }; vec makeVecToMouse () {return(new vec(mouseX-x,mouseY-y)); }; vec makeVecToBisectProjection (pt P, pt Q) {float a=this.disTo(P), b=this.disTo(Q); return(this.makeVecTo(L(P,a/(a+b),Q))); }; vec makeVecToNormalProjection (pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q)), b=dot(P.makeVecTo(Q),P.makeVecTo(Q)); return(this.makeVecTo(L(P,a/b,Q))); }; vec makeVecTowards(pt P, float d) {vec V = makeVecTo(P); float n = V.norm(); V.normalize(); V.scaleBy(d-n); return V; }; // OUTPUT TEST OR MEASURE float disTo(pt P) {return(sqrt(sq(P.x-x)+sq(P.y-y))); }; float disToMouse() {return(sqrt(sq(x-mouseX)+sq(y-mouseY))); }; boolean isInWindow() {return(((x>0)&&(x0)&&(y0; return(l); }; boolean isLeftOf(pt P, pt Q, float e) {boolean l=dot(P.makeVecTo(this),P.makeVecTo(Q).makeTurnedLeft())>e; return(l); }; boolean isInTriangle(pt A, pt B, pt C) { boolean a = this.isLeftOf(B,C); boolean b = this.isLeftOf(C,A); boolean c = this.isLeftOf(A,B); return((a&&b&&c)||(!a&&!b&&!c));}; boolean isInCircle(pt C, float r) {return d(this,C)0.000001) {x/=n; y/=n;};}; void add(vec V) {x += V.x; y += V.y;}; void add(float s, vec V) {x += s*V.x; y += s*V.y;}; void addScaled(float s, vec V) {x += s*V.x; y += s*V.y;}; void add(float u, float v) {x += u; y += v;}; void turnLeft() {float w=x; x=-y; y=w;}; void rotateBy (float a) {float xx=x, yy=y; x=xx*cos(a)-yy*sin(a); y=xx*sin(a)+yy*cos(a); }; // OUTPUT VEC vec make() {return(new vec(x,y));}; vec clone() {return(new vec(x,y));}; vec makeClone() {return(new vec(x,y));}; vec makeUnit() {float n=sqrt(sq(x)+sq(y)); if (n<0.000001) n=1; return(new vec(x/n,y/n));}; vec unit() {float n=sqrt(sq(x)+sq(y)); if (n<0.000001) n=1; return(new vec(x/n,y/n));}; vec makeScaledBy(float s) {return(new vec(x*s, y*s));}; vec makeTurnedLeft() {return(new vec(-y,x));}; vec left() {return(new vec(-y,x));}; vec makeOffsetVec(vec V) {return(new vec(x + V.x, y + V.y));}; vec makeOffsetVec(float s, vec V) {return(new vec(x + s*V.x, y + s*V.y));}; vec makeOffsetVec(float u, float v) {return(new vec(x + u, y + v));}; vec makeRotatedBy(float a) {float c=cos(a), s=sin(a); return(new vec(x*c-y*s,x*s+y*c)); }; vec makeReflectedVec(vec N) { return makeOffsetVec(-2.*dot(this,N),N);}; // OUTPUT TEST MEASURE float norm() {return(sqrt(sq(x)+sq(y)));} boolean isNull() {return((abs(x)+abs(y)<0.000001));} float angle() {return(atan2(y,x)); } // DRAW, PRINT void write() {println("("+x+","+y+")");}; void show (pt P) {line(P.x,P.y,P.x+x,P.y+y); }; void showAt (pt P) {line(P.x,P.y,P.x+x,P.y+y); }; void showArrowAt (pt P) {line(P.x,P.y,P.x+x,P.y+y); float n=min(this.norm()/10.,height/50.); pt Q=P.makeTranslatedBy(this); vec U = this.makeUnit().makeScaledBy(-n); vec W = U.makeTurnedLeft().makeScaledBy(0.3); beginShape(); Q.makeTranslatedBy(U).makeTranslatedBy(W).v(); Q.v(); W.scaleBy(-1); Q.makeTranslatedBy(U).makeTranslatedBy(W).v(); endShape(CLOSE); }; void showLabel(String s, pt P) {pt Q = P.makeTranslatedBy(0.5,this); vec N = makeUnit().makeTurnedLeft(); Q.makeTranslatedBy(3,N).showLabel(s); }; } // end vec class //************************************************************************ //**** FRAMES //************************************************************************ class frame { // frame [O I J] pt O = new pt(); vec I = new vec(1,0); vec J = new vec(0,1); frame() {} frame(pt pO, vec pI, vec pJ) {O.setTo(pO); I.setTo(pI); J.setTo(pJ); } frame(pt A, pt B, pt C) {O.setTo(B); I=A.makeVecTo(C); I.normalize(); J=I.makeTurnedLeft();} frame(pt A, pt B) {O.setTo(A); I=A.makeVecTo(B).makeUnit(); J=I.makeTurnedLeft();} frame(pt A, vec V) {O.setTo(A); I=V.makeUnit(); J=I.makeTurnedLeft();} frame(float x, float y) {O.setTo(x,y);} frame(float x, float y, float a) {O.setTo(x,y); this.rotateBy(a);} frame(float a) {this.rotateBy(a);} frame makeClone() {return(new frame(O,I,J));} void reset() {O.setTo(0,0); I.setTo(1,0); J.setTo(0,1); } void setTo(frame F) {O.setTo(F.O); I.setTo(F.I); J.setTo(F.J); } void setTo(pt pO, vec pI, vec pJ) {O.setTo(pO); I.setTo(pI); J.setTo(pJ); } void show() {float d=height/20; O.show(); I.makeScaledBy(d).showArrowAt(O); J.makeScaledBy(d).showArrowAt(O); } void showLabels() {float d=height/20; O.makeTranslatedBy(A(I,J).makeScaledBy(-d/4)).showLabel("O",-3,5); O.makeTranslatedBy(d,I).makeTranslatedBy(-d/5.,J).showLabel("I",-3,5); O.makeTranslatedBy(d,J).makeTranslatedBy(-d/5.,I).showLabel("J",-3,5); } void translateBy(vec V) {O.translateBy(V);} void translateBy(float x, float y) {O.translateBy(x,y);} void rotateBy(float a) {I.rotateBy(a); J.rotateBy(a); } frame makeTranslatedBy(vec V) {frame F = this.makeClone(); F.translateBy(V); return(F);} frame makeTranslatedBy(float x, float y) {frame F = this.makeClone(); F.translateBy(x,y); return(F); } frame makeRotatedBy(float a) {frame F = this.makeClone(); F.rotateBy(a); return(F); } float angle() {return(I.angle());} void apply() {translate(O.x,O.y); rotate(angle());} // rigid body tansform, use between pushMatrix(); and popMatrix(); void moveTowards(frame B, float s) {O.translateTowards(s,B.O); rotateBy(s*(B.angle()-angle()));} // for chasing or interpolating frames } // end frame class frame makeMidEdgeFrame(pt A, pt B) {return(new frame(A(A,B),A.makeVecTo(B)));} // creates frame for edge frame interpolate(frame A, float s, frame B) { // creates a frame that is a linear interpolation between two other frames frame F = A.makeClone(); F.O.translateTowards(s,B.O); F.rotateBy(s*(B.angle()-A.angle())); return(F); } frame twist(frame A, float s, frame B) { // a circular interpolation float d=A.O.disTo(B.O); float b=mPItoPIangle(angle(A.I,B.I)); frame F = A.makeClone(); F.rotateBy(s*b); pt M = A(A.O,B.O); if ((abs(b)<0.000001) || (abs(b-PI)<0.000001)) F.O.translateTowards(s,B.O); else { float h=d/2/tan(b/2); //else print("/b"); vec W = A.O.makeVecTo(B.O); W.normalize(); vec L = W.makeTurnedLeft(); L.scaleBy(h); M.translateBy(L); // fill(0); M.show(6); L.scaleBy(-1); L.normalize(); if (abs(h)>=0.000001) L.scaleBy(abs(h+sq(d)/4/h)); //else print("/h"); pt N = M.makeClone(); N.translateBy(L); F.O.rotateBy(-s*b,M); }; return(F); } //************************************************************************ //**** EDGES //************************************************************************ EDGE E() {return new EDGE(); } EDGE E(pt A, pt B) {return new EDGE(A,B); } EDGE E(EDGE E) {return new EDGE(E.A(),E.B()); } class EDGE { pt P[]={P(),P()}; int p=0,q=1; // pt A=P[0],B=P[1]; vec Vec() {pt A=P[p]; pt B=P[q];return V(A,B);} float n2() {pt A=P[p]; pt B=P[q];return dot(V(A,B),V(A,B));} float n() {pt A=P[p]; pt B=P[q];return d(A,B);} EDGE() {}; EDGE(pt R, pt Q) {P[0].setTo(R); P[1].setTo(Q);}; pt A() {return P[p];} pt B() {return P[q];} void setTo(EDGE E) {P[0].setTo(E.A()); P[1].setTo(E.B());}; void setTo(pt A, pt B) {P[0].setTo(A); P[1].setTo(B);}; void show() {pt A=P[p]; pt B=P[q]; A.show(2); A.to(B);} void show(pt c) {pt A=P[p]; pt B=P[q];pt C=put(c); A.show(2); A.to(B); B.to(C); C.to(A);} pt grab(pt C) {pt A=P[p]; pt B=P[q]; return P(dot(V(A,C),V(A,B))/n2(),dot(V(A,C),R(V(A,B)))/n2() );} pt put(pt C) {pt A=P[p]; pt B=P[q];return T(T(A,C.x,V(A,B)),C.y,R(V(A,B))); } void pickClosestVertex(pt M) {if (M.disTo(P[0])<=M.disTo(P[1])) {p=0;q=1;} else {p=1;q=0;}; } void scaleBy(float s) {P[q].setTo(L(P[q],s,P[p]));}; void translateBy(vec V) {P[p].translateBy(V); P[q].translateBy(V); }; void rotateBy(float aa) {P[q].rotateBy(aa,P[p]);}; void drag(vec V) {pt A=P[p]; pt B=P[q]; float d=d(A,B); P[p].translateBy(V); P[q].setTo(T(B,d-d(A,B),A));}; int p() {return p;} } // end class EDGE EDGE S(EDGE E0, float t, EDGE E1) {return E(L(E0.A(),t,E1.A()),L(E0.B(),t,E1.B()));} EDGE R(EDGE E0, float t, EDGE E1) { float a = angle(E0.Vec(),E1.Vec()); float s = E1.n() / E0.n(); if ((abs(a)<0.01)&&(abs(s-1)<0.01)) return E(L(E0.A(),t,E1.A()),L(E0.B(),t,E1.B())); // pure translation, no rotation, no scaling ********* if (isSame(E0.A(),E1.A(),0.01)) return E( E0.A() , spiral(E0.B(),E0.A(),s,a,t) ); if (isSame(E0.B(),E1.B(),0.01)) return E( spiral(E0.A(),E0.B(),s,a,t) , E0.B() ); if (badSpiral(E0,E1)) {strokeWeight(3); stroke(green); E0.show(); stroke(blue); E1.show(); strokeWeight(1); }; pt G=spiralCenter(E0,E1); return E( spiral(E0.A(),G,s,a,t) , spiral(E0.B(),G,s,a,t) ); } EDGE R(EDGE E0, float t, EDGE E1, EDGE E) { float a = angle(E0.Vec(),E1.Vec()); float s = E1.n() / E0.n(); pt G=spiralCenter(E0,E1); return E( spiral(E.A(),G,s,a,t) , spiral(E.B(),G,s,a,t) ); } EDGE bezier(EDGE A, EDGE B, EDGE C, EDGE D, float t) {return R( R( R(A,t,B) ,t, R(B,t,C) ) ,t, R( R(B,t,C) ,t, R(C,t,D) ) ); } EDGE catmull(EDGE A, EDGE B, EDGE C, EDGE D, float t) {return bezier(B,R(A,1./6,C,B),R(D,1./6,B,C),C,t); } void catmullConstruction(EDGE A, EDGE B, EDGE C, EDGE D) {R(A,1./6,C,B).show(); R(D,1./6,B,C).show(); }; boolean edgesIntersect(pt A, pt B, pt C, pt D) {boolean hit=true; if (isRightTurn(A,B,C)==isRightTurn(A,B,D)) hit=false; if (isRightTurn(C,D,A)==isRightTurn(C,D,B)) hit=false; return hit; } boolean edgesIntersect(pt A, pt B, pt C, pt D,float e) { return ((A.isLeftOf(C,D,e) && B.isLeftOf(D,C,e))||(B.isLeftOf(C,D,e) && A.isLeftOf(D,C,e)))&& ((C.isLeftOf(A,B,e) && D.isLeftOf(B,A,e))||(D.isLeftOf(A,B,e) && C.isLeftOf(B,A,e))) ; } pt linesIntersection(pt A, pt B, pt C, pt D) {vec AB = A.makeVecTo(B); vec CD = C.makeVecTo(D); vec N=CD.makeTurnedLeft(); vec AC = A.makeVecTo(C); float s = dot(AC,N)/dot(AB,N); return A.makeTranslatedBy(s,AB); } //************************************************************************ //**** SPIRALS //************************************************************************ pt spiral(pt A, pt G, float s, float a) {return L(G,s,R(A,a,G));} pt spiral(pt A, pt G, float s, float a, float t) {return L(G,1+t*(s-1),R(A,t*a,G));} pt spiralCenter(EDGE A, EDGE B) {return spiralCenter(A.A(),A.B(),B.A(),B.B()); } pt spiralCenterM(pt A, pt B, pt C, pt D) { // derived from a formula provided by Michael Hale CS3451-S2008, but does not handle cases where AC.x==0 vec AB = V(A,B); vec CD = V(C,D); vec AC = V(A,C); vec M = S(AB,-1,CD); float c=dot(M,AC), s=dot(R(M),AC), h=sqrt(sq(c)+sq(s)); c/=h; s/=h; // to avoid atan2 call float t = 0; if(abs(M.x*c-M.y*s)>0.0000001) t=AC.x/( M.x*c-M.y*s); return T(A,t,V(AB.x*c-AB.y*s,AB.x*s+AB.y*c)); } pt spiralCenter(pt A, pt B, pt C, pt D) { // Jarek's version: G has the same relative coordinates in (A,B) as in (C,D) vec AB = V(A,B); vec CD = V(C,D); vec CA=V(C,A); float AB2=dot(AB,AB); float CD2=dot(CD,CD); vec ab=S(1./AB2,AB,-1./CD2,CD); vec de=S(1./AB2,R(AB),-1./CD2,R(CD)); float a=ab.x; float b=ab.y; float c=dot(CD,CA)/CD2; float d=de.x; float e=de.y; float f=dot(R(CD),CA)/CD2; float det=a*e-b*d; // if (abs(det)<0.00000000001) {println("spiralCenter: det==0"); return P(0,0); } float x=(c*e-b*f)/det; float y=(a*f-c*d)/det; // println(a+"*"+x+"+"+b+"*"+y+" = "+(a*x+b*y)+" = "+c); // println(d+"*"+x+"+"+e+"*"+y+" = "+(d*x+e*y)+" = "+f); // println(); // pt P=P(x,y); return put(A,B,P); return P(A.x+x,A.y+y); } boolean badSpiral(EDGE A, EDGE B) {return badSpiral(A.A(),A.B(),B.A(),B.B()); } boolean badSpiral(pt A, pt B, pt C, pt D) { vec AB = V(A,B); vec CD = V(C,D); vec CA=V(C,A); float AB2=dot(AB,AB); float CD2=dot(CD,CD); vec ab=S(1./AB2,AB,-1./CD2,CD); vec de=S(1./AB2,R(AB),-1./CD2,R(CD)); float a=ab.x; float b=ab.y; float c=dot(CD,CA)/CD2; float d=de.x; float e=de.y; float f=dot(R(CD),CA)/CD2; float det=a*e-b*d; return abs(det)<0.00000000001; } //************************************************************************ //**** RAYS //************************************************************************ RAY ray(pt A, pt B) {return new RAY(A,B); } RAY ray(pt Q, vec T) {return new RAY(Q,T); } RAY ray(pt Q, vec T, float d) {return new RAY(Q,T,d); } RAY ray(RAY R) {return new RAY(R.Q,R.T,R.r); } RAY leftTangentToCircle(pt P, pt C, float r) {return tangentToCircle(P,C,r,-1); } RAY rightTangentToCircle(pt P, pt C, float r) {return tangentToCircle(P,C,r,1); } RAY tangentToCircle(pt P, pt C, float r, float s) { float n=d(P,C); float w=sqrt(sq(n)-sq(r)); float h=r*w/n; float d=h*w/r; vec T = S(d,U(V(P,C)),s*h,R(U(V(P,C)))); return ray(P,T,w);} class RAY { pt Q = new pt(300,300); // start vec T = new vec(1,0); // direction float r = 50; // length for display arrow and dragging float d = 300; // distance to hit RAY () {}; RAY (pt pQ, vec pT) {Q.setTo(pQ);T.setTo(U(pT)); }; RAY (pt pQ, vec pT, float pd) {Q.setTo(pQ);T.setTo(U(pT)); d=max(0,pd);}; RAY(pt A, pt B) {Q.setTo(A); T.setTo(U(V(A,B))); d=d(A,B);} RAY(RAY B) {Q.setTo(B.Q); T.setTo(B.T); d=B.d; T.normalize();} void drag() {pt P=T(Q,r,T); if(P.disToMouse()0.000001) t = -dot(N,V(A,Q))/n; return t;} boolean hitsEdge(pt A, pt B) {boolean hit=isRightOf(A,Q,T)!=isRightOf(B,Q,T); if (isRightTurn(A,B,Q)==(dot(T,R(V(A,B)))>0)) hit=false; return hit;} // if hits float disToEdge(pt A, pt B) {vec N = U(R(V(A,B))); float t=0; float n=dot(N,T); if(abs(n)>0.000001) t=-dot(N,V(A,Q))/n; return t;} // distance to edge along ray if hits pt intersectionWithEdge(pt A, pt B) {return at(disToEdge(A,B));} // hit point if hits RAY reflectedOfEdge(pt A, pt B) {pt X=intersectionWithEdge(A,B); vec V =T.makeReflectedVec(R(U(V(A,B)))); float rd=d-disToEdge(A,B); return ray (X,V,rd); } // bounced ray RAY surfelOfEdge(pt A, pt B) {pt X=intersectionWithEdge(A,B); vec V = R(U(V(A,B))); float rd=d-disToEdge(A,B); return ray (X,V,rd); } // bounced ray float disToCircle(pt C, float r) { return rayCircleIntesectionParameter(Q,T,C,r);} // distance to circle along ray pt intersectionWithCircle(pt C, float r) {return at(disToCircle(C,r));} // intersection point if hits boolean hitsCircle(pt C, float r) {return disToCircle(C,r)!=-1;} // hit test RAY reflectedOfCircle(pt C, float r) {pt X=intersectionWithCircle(C,r); vec V =T.makeReflectedVec(U(V(C,X))); float rd=d-disToCircle(C,r); return ray (X,V,rd); } } //************************************************************************ //**** TRIANGLES //************************************************************************ //************************************************************************ //**** INTEGRAL PRPERTIES, AREA, CENTER //************************************************************************ float trapezeArea(pt A, pt B) {return((B.x-A.x)*(B.y+A.y)/2.);} pt trapezeCenter(pt A, pt B) { return(new pt(A.x+(B.x-A.x)*(A.y+2*B.y)/(A.y+B.y)/3., (A.y*A.y+A.y*B.y+B.y*B.y)/(A.y+B.y)/3.) ); } //************************************************************************ //**** CIRCLES //************************************************************************ void showArcThrough (pt A, pt B, pt C) { pt O = CenterCircum ( A, B, C); //O.show(3); float r=2.*(O.disTo(A) + O.disTo(B)+ O.disTo(C))/3.; float a = V(O,A).angle(); //average(O,A).showLabel(str(toDeg(a))); float c = V(O,C).angle(); //average(O,C).showLabel(str(toDeg(c))); if(!isRightTurn(A,B,C)) arc(O.x,O.y,r,r,a,c); else arc(O.x,O.y,r,r,c,a); } void showArcThroughXXX (pt A, pt B, pt C) { pt O = CenterCircum ( A, B, C); float r=2.*(O.disTo(A) + O.disTo(B)+ O.disTo(C))/3.; float a = O.vecTo(A).angle(); //average(O,A).showLabel(str(toDeg(a))); float c = O.vecTo(C).angle(); //average(O,C).showLabel(str(toDeg(c))); float w = angle(O.vecTo(C),O.vecTo(A)); if(!isRightTurn(A,B,C)) w=TWO_PI-w; if (w>TWO_PI) w-=TWO_PI; if (w<0) w+=TWO_PI; //println(w*180/PI); float d=w/2.; if(!isRightTurn(A,B,C)) arc(O.x,O.y,r,r,a-d,c+d); else arc(O.x,O.y,r,r,c-d,a+d); } float interpolateAngles(float a, float t, float b) {if(bPI) m-=TWO_PI; return m;} pt pointOnArcThrough (pt A, float t, pt B, pt C) { pt X=new pt(); pt O = CenterCircum ( A, B, C); float r=(O.disTo(A) + O.disTo(B)+ O.disTo(C))/3.; float a = V(O,A).angle(); float ab = Oto2PIangle(angle(V(O,A),V(O,B))); if(isRightTurn(A,B,C)) ab-=TWO_PI; return ptOnCircle(a+t*ab,O,r); } pt pointOnArcThrough (pt A, pt B, float t, pt C) { pt X=new pt(); pt O = CenterCircum ( A, B, C); float r=(O.disTo(A) + O.disTo(B)+ O.disTo(C))/3.; float b = V(O,B).angle(); float bc = Oto2PIangle(angle(V(O,B),V(O,C))); if(isRightTurn(A,B,C)) bc-=TWO_PI; return ptOnCircle(b+t*bc,O,r); } pt ptOnCircle(float a, pt O, float r) {return P(r*cos(a)+O.x,r*sin(a)+O.y);} float edgeCircleIntesectionParameter (pt A, pt B, pt C, float r) { // computes parameter t such that A+tAB is on circle(C,r) vec T = V(A,B); float n = n(T); T.normalize(); float t=rayCircleIntesectionParameter(A,T,C,r); return t*n; } float rayCircleIntesectionParameter (pt A, vec T, pt C, float r) { // computes parameter t such that A+tT is on circle(C,r) or -1 vec AC = V(A,C); float d=dot(AC,T); float h = dot(AC,R(T)); float t=-1; if (abs(h)