//************************************************************************
//**** 2D GEOMETRY CLASSES AND UTILITIES, Jarek Rossignac, REVISED MARCH 2007
//************************************************************************

//************************************************************************
//**** POINTS
//************************************************************************
class pt { float x=0,y=0; 
  pt () {}
  pt (float px, float py) {x = px; y = py;};
  pt (pt P) {x = P.x; y = P.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 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 translateBy(float s, vec V) {x += s*V.x; y += s*V.y;};   
  void addScaledPt(float s, pt V) {x += s*V.x; y += s*V.y;};   
  void translateBy(float u, float v) {x += u; y += v;};
  void translateTowards(float s, pt P) {x=s*(x-P.x)+P.x;  y=s*(y-P.y)+P.y; };
  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; };
  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; };
  
  // OUTPUT POINT
  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 projectOnEdge(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 makeRotatedBy(float s) {         };
//  pt makeRotatedBy(float s, pt P) {         };

   // OUTPUT VEC
  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 makeVecToCenter (pt P, pt Q) {return(new vec((P.x+Q.x)/2.0-x,(P.y+Q.y)/2.0-y)); };
  vec makeVecToCenter (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(interpolate(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(interpolate(P,a/b,Q))); };
 
  // OUTPUT TEST 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)&&(x<height)&&(y>0)&&(y<height)));};
  boolean projectsBetween(pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q)), b=dot(P.makeVecTo(Q),P.makeVecTo(Q)); return((0<a)&&(a<b)); };
  float ratioOfProjectionBetween(pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q)), b=dot(P.makeVecTo(Q),P.makeVecTo(Q)); return(a/b); };
  float disToEdge(pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q).makeUnit().makeTurnedLeft()); return(abs(a)); };
  float triArea(pt P, pt Q) {float a=dot(P.makeVecTo(this),P.makeVecTo(Q).makeTurnedLeft()); return(a/2.); };
  boolean isLeftOf(pt P, pt Q) {boolean l=dot(P.makeVecTo(this),P.makeVecTo(Q).makeTurnedLeft())>0; 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) );};

  // DRAW , PRINT
  void show() {ellipse(x, y, height/200, height/200); }; 
  void show(float r) {ellipse(x, y, 2*r, 2*r); }; 
  void v() {vertex(x,y);};
  void write() {println("("+x+","+y+")");};
  void showLabel(String s, vec D) {text(s, x+D.x,y+D.y);  };
  void showLabel(String s) {text(s, x+5,y+4);  };
  void showLabel(int i) {text(str(i), x+5,y+4);  };
  void showLabel(String s, float u, float v) {text(s, x+u, y+v);  };
  void showSegmentTo (pt P) {line(x,y,P.x,P.y); }; 

  } // end of pt class
  
pt average(pt A, pt B) {return(new pt((A.x+B.x)/2.0,(A.y+B.y)/2.0)); };
pt average(pt A, pt B, pt C) {return(new pt((A.x+B.x+C.x)/3.0,(A.y+B.y+C.y)/3.0)); };
pt scaledAverage(pt A, pt B, pt C, float s) {return(new pt((A.x+B.x+C.x)/s,(A.y+B.y+C.y)/s)); };
pt interpolate(pt A, float s, pt B) {return(new pt(A.x+s*(B.x-A.x),A.y+s*(B.y-A.y))); };
pt mouse() {return(new pt(mouseX,mouseY));};
pt pmouse() {return(new pt(pmouseX,pmouseY));};
pt mouseInWindow() {float x=mouseX, y=mouseY; x=max(x,0); y=max(y,0); x=min(x,height); y=min(y,height);  return(new pt(x,y));}; 
pt screenCenter() {return(new pt(height/2,height/2));}
boolean isLeftTurn(pt A, pt B, pt C) {return(dot(A.makeVecTo(B).makeTurnedLeft() , B.makeVecTo(C) )>0); };
pt s(pt A, float s, pt B) {return(new pt(A.x+s*(B.x-A.x),A.y+s*(B.y-A.y))); };
pt b(pt A, pt B, pt C, float s) {return( s(s(B,s/4.,A),0.5,s(B,s/4.,C))); };                          // tucks in a vertex towards its neighbors
pt f(pt A, pt B, pt C, pt D, float s) {return( s(s(A,1.+(1.-s)/8.,B) ,0.5, s(D,1.+(1.-s)/8.,C))); };    // bulges out a mid-edge point 
pt weightedSum(float a, pt A, float b, pt B) {pt P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B); return(P);}
pt weightedSum(float a, pt A, float b, pt B, float c, pt C) {pt P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B);  P.addScaledPt(c,C); return(P);}
pt weightedSum(float a, pt A, float b, pt B, float c, pt C, float d, pt D) 
        {pt P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B); P.addScaledPt(c,C); P.addScaledPt(d,D); return(P);}
boolean leftTurn(pt A, pt B, pt C) {return(C.isLeftOf(A,B));}
//************************************************************************
//**** VECTORS
//************************************************************************
class vec { float x=0,y=0; 
  vec () {};
  vec (vec V) {x = V.x; y = V.y;};
  vec (float px, float py) {x = px; y = py;};
 
 // MODIFY
  void setTo(float px, float py) {x = px; y = py;}; 
  void setTo(pt P, pt Q) {x = Q.x-P.x; y = Q.y-P.y;}; 
  void setTo(vec V) {x = V.x; y = V.y;}; 
  void scaleBy(float f) {x*=f; y*=f;};
  void scaleBy(float u, float v) {x*=u; y*=v;};
  void normalize() {float n=sqrt(sq(x)+sq(y)); if (n>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 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 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 makeScaledBy(float s) {return(new vec(x*s, y*s));};
  vec makeTurnedLeft() {return(new vec(-y,x));};
  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) {return(new vec(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a))); };
 
  // 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 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); }; 
  } // end vec class
 
vec average(vec A, vec B) {return(new vec((A.x+B.x)/2.0,(A.y+B.y)/2.0)); };
float dot(vec U, vec V) {return(U.x*V.x+U.y*V.y); };
vec interpolate(vec A, float s, vec B) {return(new vec(A.x+s*(B.x-A.x),A.y+s*(B.y-A.y))); };

float angle(vec U, vec V) {return(atan2(dot(U.makeTurnedLeft(),V),dot(U,V))); };
float angle(vec V) {return(atan2(V.y,V.x)); };
float mPItoPIangle(float a) { if(a>PI) return(mPItoPIangle(a-2*PI)); if(a<-PI) return(mPItoPIangle(a+2*PI)); return(a);};
float toDeg(float a) {return(a*180/PI);}
float toRad(float a) {return(a*PI/180);}

 //************************************************************************
//**** FRAMES
//************************************************************************
class frame {
  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(average(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());}
  } // end frame class
 
frame makeMidEdgeFrame(pt A, pt B) {return(new frame(average(A,B),A.makeVecTo(B)));}
frame interpolate(frame A, float s, frame B) {
  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) {
  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 = average(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);
  }
  

//************************************************************************
//**** CURVES
//************************************************************************
pt cubicBezier(pt A, pt B, pt C, pt D, float t) {return( s( s( s(A,t,B) ,t, s(B,t,C) ) ,t, s( s(B,t,C) ,t, s(C,t,D) ) ) ); }
void drawCubicBezier(pt A, pt B, pt C, pt D) { beginShape();  for (float t=0; t<=1; t+=0.02) {cubicBezier(A,B,C,D,t).v(); };  endShape(); }
float cubicBezierAngle (pt A, pt B, pt C, pt D, float t) {pt P = s(s(A,t,B),t,s(B,t,C)); pt Q = s(s(B,t,C),t,s(C,t,D)); vec V=P.makeVecTo(Q); float a=atan2(V.y,V.x); return(a);}  

void drawParabolaInHat(pt A, pt B, pt C, int rec) {
   if (rec==0) { B.showSegmentTo(A); B.showSegmentTo(C); } 
   else { 
     float w = (A.makeVecTo(B).norm()+C.makeVecTo(B).norm())/2;
     float l = A.makeVecTo(C).norm()/2;
     float t = l/(w+l);
     pt L = new pt(A); 
     L.translateBy(t,A.makeVecTo(B)); 
     pt R = new pt(C); R.translateBy(t,C.makeVecTo(B)); 
     pt M = average(L,R);
     drawParabolaInHat(A,L, M,rec-1); drawParabolaInHat(M,R, C,rec-1); 
     };
   };

pt circumCenter (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);
  };

void showArcThrough (pt A, pt B, pt C) {
   pt O = circumCenter ( A,  B,  C);
   float r=2.*(O.disTo(A) + O.disTo(B)+ O.disTo(C))/3.;
   float a = O.makeVecTo(A).angle(); //average(O,A).showLabel(str(toDeg(a)));
   float c = O.makeVecTo(C).angle(); //average(O,C).showLabel(str(toDeg(c)));
   if(isLeftTurn(A,B,C)) arc(O.x,O.y,r,r,a,c); else arc(O.x,O.y,r,r,c,a);
    }