/**
    @file     geo2d.pde
    @brief    This is the soul of this program, you might laugh, but it's true.\n
              This file contains all of the basis functions and classes for action.pde.\n
              We added one important class to this file. The <B>triangle</B> class.
*/

//!************************************************************************
//!**** 2D GEOMETRY CLASSES AND UTILITIES, Jarek Rossignac, REVISED MARCH 2007
//!**** ADD TRIANGLE CLASS, Sungbae Kim, Jeonggyu Lee, AUGUST 2007
//!************************************************************************

//!************************************************************************
//!**** POINTS
//!************************************************************************
/**
    @brief    \class pt

               This class represent a point on 2 dimensional space.
*/
public class pt
{
  float x=0,y=0; //!!< The default value of x and y is both 0.
  pt () {}


/**
    @brief     A constructer with two specific parameters.
    @see       pt ()
    @see       pt (pt P)
    @param px float value of x-coordinate.
    @param px float value of y-coordinate.
    @return pt
*/
  //!! Basic constructure without any parameter.
  pt (float px, float py)
  {
	x = px;
        y = py;
  };

/**
    @brief     A constructer with two specific parameters.
    @see       pt ()
    @see       pt (float px, float py)
    @param P  pt of any point.
    @return pt
*/
  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 setTo(pt P, vec v, float t) { x = P.x + v.x * t; y = P.y + v.y * t; }
  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)); };


   // 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)); };
  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); }; 
  void showDashTo(pt P, float len) 
  {
    vec v = new vec();
    v.setTo(this, P);
    float norm = v.norm();
    float lenSum = 0;
    v.normalize();
    pt S = new pt();
    pt E = new pt();
    S.setTo(this);
    while ( lenSum < norm )
    {
      E.setTo(S);
      E.translateBy(len, v);
      line(S.x, S.y, E.x, E.y);
      E.translateBy(len*2, v);
      S.setTo(E);
      lenSum = lenSum + 3*len;
    }
  }
  } // 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)); };

/**
    @brief     This function calculate the average of three points A, B and C by add them then divide by 3.
    @param A point A
    @param B point B
    @param C point C
    @return pt
*/
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));}; 

/**
    @brief     This is the Professor Jarek's function to calculate the center of the screen.
    @return		pt of center position.
*/
pt screenCenter() {return(new pt(width/2,height/2));}

/**
    @brief     This is the Professor Jarek's notorious function, isLeftTurn.
               First we make a vector AB. Then we turn it left and dot-product vector BC.
	       By inspecting result of it, we know it's left-turn or not(Remember there's cosine in dot-product).
    @param A point A
    @param B point B
    @param C point C
    @see   dot(vec U, vec V)
    @return true if it's left-turn, false otherwise.
*/
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));}
float trapezeArea(pt A, pt B) {return((B.x-A.x)*(B.y+A.y)/2.);}

//************************************************************************
//**** VECTORS
//************************************************************************
/**
    @brief    \class vec

               This class represent a vector.
*/
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)); };

/**
    @brief     This is the Professor Jarek's dot-product function..
               As we all know, in a given two vector U(x,y) and V(x,y) dot-product
			   can be calculate by following equation.\n
			   U dot V = Ux*Vx + Uy*Vy

    @param U vector U
    @param V vector V
    @return float
*/
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
//************************************************************************
/**
    @brief    \class frame

               This class represent a frame.
*/
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); 
     };
   };
/**
    @brief    computes the center of a circumscirbing circle to triangle (A,B,C)
*/
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);
    }

//************************************************************************
//**** TRIANGLES
//************************************************************************
//!constant variable
float nonlinear = -1000000;
//!constant variable
float tolerence = 0.1;

/**
    @brief    \class intersect

               This class for calculate intersection.
*/
//This class for intersection points
class intersect
{
  pt logicalPt = new pt();
  pt realStartPt = new pt();  // this means the start point to display
  pt realEndPt = new pt();    // this means the end point to display 

/**
    @brief    constructor
*/  
  intersect() { }
  
  void reset()
  {
    logicalPt.setTo(-1,-1);
    realStartPt.setTo(-1,-1);
    realEndPt.setTo(-1,-1);
  }
  void setTo(intersect inter)
  {
    logicalPt.setTo(inter.logicalPt);
    realStartPt.setTo(inter.realStartPt);
    realEndPt.setTo(inter.realEndPt);
  } 
}




/**
    @brief    \class triangle

               This class represent a triangle.
*/
class triangle
{
  pt[] Pt = new pt[3];                       // vertices
  String[] Label = new String[3];            // labels of vertices 
  intersect[] Intersect = new intersect[2];  // intersection points with a line
  int  weight = 1;                           // weight of edges 
  
/**
    @brief    constructor
*/
  triangle()
  {
	for ( int i = 0; i < 3; i++ ) { Pt[i] = new pt(); } for ( int i = 0; i < 2; i++ ) { Intersect[i] = new intersect(); }
  }
  
  // set functions
  void setPt (pt A, pt B, pt C) 
  { 
    Pt[0].setTo(A);
    Pt[1].setTo(B);
    Pt[2].setTo(C);
  }
  void setLabel (String A, String B, String C)
  {
    Label[0] = A; Label[1] = B; Label[2] = C;
  }
  
  void setWeight(int w) { weight = w; }
  

  /**
    @brief   check whether P is on a edge of this triangle
		@param P Point P
		@return TRUE when P is on a edge of this triangle
  */
  boolean isOnEdge(pt P)
  {
    for ( int i = 0; i < 3; i++ )
    {
      if (isOnLine(Pt[i%3], Pt[(i+1)%3], P))return true;
    }
    return false;
  }
  
  /**
    @brief   check whether P is in this triangle
		@param P Point P
		@return TRUE when P is in this triangle
  */
  boolean isIn(pt P)
  {
    // from pt class
    return (P.isInTriangle(Pt[0], Pt[1], Pt[2]));
  }

/**
    @brief   display two triangles accodring to the operation such as union, intersection, difference
		@param Triangle2 a triangle
		@param lineColor the color of edges
		@param fillColor the color to fill in triangles
		@param intersectionColor the color to fill in the intersection area
		@param op operation : 0 = nothing, 1 = uniion, 2 = intersection, 3 = difference(this-Triangle2)
		@return none
*/
  void show(triangle Triangle2, color lineColor, color fillColor, color intersectionColor, int op)
  {
    // fill trinagles
    if ( op == 0 )
    {
      show(lineColor, fillColor);
      Triangle2.show(lineColor, fillColor);
    }
    if ( op == 1 ) // union
    {
      show(lineColor, intersectionColor);
      Triangle2.show(lineColor, intersectionColor);
    }
    else if ( op == 2 ) // intersection
    {
      show(lineColor, fillColor);
      Triangle2.show(lineColor, fillColor);
      showIntersection(Triangle2, intersectionColor);
    }
    else if ( op == 3 ) // difference
    {
      // to diplay difference, we use just the intersection area
      show(lineColor, intersectionColor);
      Triangle2.show(lineColor, fillColor);
      showIntersection(Triangle2, fillColor);
    }
    
    // draw edges
    for (int i = 0; i < 3; i++)
    {  
      drawLine(Triangle2.Pt[i%3], Triangle2.Pt[(i+1)%3]);
      Triangle2.drawLine(Pt[i%3], Pt[(i+1)%3]);
    }
  }

/**
    @brief   diplay the intersection area of two triangles
		@param Triangle2 a triangle
		@param intersectionColor the color to fill in the interection area
		@return none
*/
  void showIntersection(triangle Triangle2, color intersectionColor)
  {
    pt[] vertexPoly = new pt[6];
    for ( int i = 0; i < 6; i++ )
    {
      vertexPoly[i] = new pt();
    }
    
    // Find out vertices of the intersection area
    // vertices of the intersection area are intersection points and each vertices are in the other triangle.
    int numPoly = 0;
    for ( int i = 0; i < 3; i++ )
    {
      pt S = Triangle2.Pt[i%3];
      pt E = Triangle2.Pt[(i+1)%3];
      if ( isIn(S) ) { vertexPoly[numPoly].setTo(S); numPoly++; }
      
      findIntersect( S, E );
      int num = getNumOfIntersect();
      for ( int k = 0; k < num; k++ )
      {
        vertexPoly[numPoly].setTo(Intersect[k].logicalPt);
        numPoly++;
      }
    }
    
    for ( int i = 0; i < 3; i++ )
    {
      if ( Triangle2.isIn(Pt[i]) )
      {
        vertexPoly[numPoly].setTo(Pt[i]);
        numPoly++;
      }  
    }
   
    // visualize
    if ( numPoly >= 3 )
    {
      // To keep the rotaion direction right or left, sort vertices according to cosine of angle
      // 
      float[] dotproduct = new float[numPoly];
      int[] ptIndex = new int[numPoly];
      vec vec1 = new vec();
      vec1.setTo(vertexPoly[0], vertexPoly[1]);
      vec1.normalize();
      
      dotproduct[0] = 1;
      dotproduct[1] = 1;
      ptIndex[0] = 0;
      ptIndex[1] = 1;
      for ( int i = 2; i < numPoly; i++ )
      {
        vec vec2 = new vec();
        vec2.setTo(vertexPoly[0], vertexPoly[i]);
        vec2.normalize();
        dotproduct[i] = dot(vec1, vec2);
        ptIndex[i] = i;
      }
      
      // just bubble sort
      for( int i = 0; i < numPoly-1; i++ )
      {
        for ( int j = 0; j < numPoly-1-i; j++ )
        {
          if ( dotproduct[j] < dotproduct[j+1] )
          {
            float tempProduct = dotproduct[j];
            dotproduct[j] = dotproduct[j+1];
            dotproduct[j+1] = tempProduct;
             
            int tempIndex = ptIndex[j];
            ptIndex[j] = ptIndex[j+1];
            ptIndex[j+1] = tempIndex; 
          }
        }
     }
      
     // fill area
     noStroke();
     fill(intersectionColor); 
     beginShape(); 
     for (int i=0; i<numPoly; i++)
     { 
       vertexPoly[ptIndex[i]].v(); 
       //println("(index, dot prodcut) = " + "(" + ptIndex[i] + "," + dotproduct[i] + ")");
     } 
     endShape(CLOSE);
   }
  }
  
  	/**
    @brief     This function show this triangle
		@param lineColor the color of edges
		@param fillColor the color to fill
		@return none
*/
  void show (color lineColor, color fillColor)
  {
    strokeWeight(weight);
    stroke(lineColor);   
    fill(fillColor); beginShape(); for (int i=0; i<3; i++){ Pt[i].v(); } endShape(CLOSE);
  }
  
	/**
    @brief     This function show labels of vertices
		@param labelColor the color of label
		@return none
*/ 
  void showLabel (color labelColor) 
  {
    fill(labelColor); for (int i = 0; i<3; i++) Pt[i].showLabel(Label[i]);
  }

/**
    @brief     This function draw a line PQ considering this triangle
		@param P start point of a line
		@param Q end point of a line
		@return none
*/  
  void drawLine(pt P, pt Q)
  {
    strokeWeight(weight);
    
    // if PQ is on an edge
    if ( findLineOnEdge( P, Q ) )
    {
      //draw line from A to first pt in green color
      if ( P.x != Intersect[0].logicalPt.x || P.y != Intersect[0].logicalPt.y)
      {
        stroke(green); P.showSegmentTo(Intersect[0].logicalPt);
      }      
      //darw line from 1st pt to 2nd pt in oragne color
      stroke(orange); Intersect[0].logicalPt.showSegmentTo(Intersect[1].logicalPt);
      
      //draw line from 2nd pt to B in green color
      if ( Q.x != Intersect[1].logicalPt.x || Q.y != Intersect[1].logicalPt.y)
      {
        stroke(green); Intersect[1].logicalPt.showSegmentTo(Q);
      }  
    }
    else
    {
      color[] linecolor = new color[3];
      if ( isIn(P) )
      {
        linecolor[0] = red;
        linecolor[1] = green;
        linecolor[2] = green; //dummy
      }
      else
      {
        linecolor[0] = green;
        linecolor[1] = red;
        linecolor[2] = green;
      }
  
      findIntersect(P, Q);
  
      int num;
      num = getNumOfIntersect();
 
      pt start = P;
      for ( int i = 0; i < num; i++ )
      {
        stroke(linecolor[i]); start.showSegmentTo(Intersect[i].realStartPt);
        stroke(orange); Intersect[i].realStartPt.showSegmentTo(Intersect[i].realEndPt);
        start = Intersect[i].realEndPt;
      } 
      stroke(linecolor[num]); start.showSegmentTo(Q);
    }
    
    // display vertices
    stroke(green);
    if (isOnEdge(P)) fill(orange);
    else if (isIn(P)) fill(red);
    else fill(green); 
    P.show(4);
  
    if (isOnEdge(Q)) fill(orange);
    else if (isIn(Q)) fill(red);
    else fill (green); 
    Q.show(4);
  }
  
/**
    @brief     This function check whether a line PQ is on an edge of this triangle
		@brief     and  get start and end points of overlapped part.
    @param P start point of a line
		@param Q end point of a line
		@return TRUE if a line PQ is on an edge of this triangle
*/
  boolean findLineOnEdge(pt P, pt Q)
  {
    vec v1 = new vec();
    v1.setTo( P, Q );
    
    resetIntersect();
    for ( int i = 0; i < 3; i++ )
    {
      pt S = Pt[i%3];        // start point of edge
      pt E = Pt[(i+1)%3];    // end point of edge
      
      // if PQ is on an edge, both of P and Q are on the extension of an edge
      float t1 = getRelativePosOnLine(P, Q, S);
      float t2 = getRelativePosOnLine(P, Q, E);
      if ( t1 != nonlinear && t2 != nonlinear )
      {
        if ( t1 < 0 ) t1 = 0;
        if ( t1 > 1 ) t1 = 1;
        if ( t2 < 0 ) t2 = 0;
        if ( t2 > 1 ) t2 = 1;
        
        // sort 
        float st, et;
        if ( t1 < t2 ) { st = t1; et = t2; }
        else { st = t2; et = t1; }
        Intersect[0].logicalPt.setTo(P, v1, st);
        Intersect[1].logicalPt.setTo(P, v1, et);
        return true;
      }
    }
    return false;  
  }
  void resetIntersect ()
  {
    for ( int i = 0; i < 2; i++ )
    {
      Intersect[i].reset(); 
    }
  }
  
	/**
    @brief     This function find intersection points of PQ and this triangle.
		@brief			
    @param P Start point of a line
    @param Q End point of a line
    @return none
*/  
  void findIntersect (pt P, pt Q)
  {
    resetIntersect();  // reset intersection points
    
    // make a vector for line (PQ)
    vec v1 = new vec();
    v1.setTo( P, Q );
    float normPQ = v1.norm();
    v1.normalize();                // normalize

		// find intersection points
    int i = 0;
    int j = 0;
    float oldt = -1;
    float t = -1;
    for ( i = 0, j = 0; i < 3 && j < 2; i++ )
    {
      pt S = Pt[i%3];        // start point of edge
      pt E = Pt[(i+1)%3];    // end point of edge
      if(isCross(P,Q, S, E))
      {
        t = findPoint(P,Q,S,E);
        if ( oldt != t ) // check whether this vertex has been already used.
        {
          int index = 0;

					// arrage order of intersection points according to distance from P
          if ( t < oldt && j > 0)
          { 
            Intersect[j].setTo(Intersect[j-1]);
            index = j-1;
          }
          else
          {
            index = j;
          }
          
          Intersect[index].logicalPt.setTo(P, v1, t);
          
          // calculate boundies of the intersetions
          vec v2 = new vec();
          v2.setTo( S, E );
          v2.normalize();                // normalize
          float cosAngle = abs(dot(v1, v2.makeTurnedLeft())); 
          if (cosAngle < 0.000001 )
          {
             cosAngle = 0.0000001;
          }   
          
          vec vecSI = new vec();
          vecSI.setTo( S, Intersect[index].logicalPt );
          
          vec vecEI = new vec();
          vecEI.setTo( E, Intersect[index].logicalPt );
          
          float delta = weight*0.5 / cosAngle;
          
          float realStart = t-delta;
          float realEnd = t+delta;       
 
          float edge1Norm;
          float edge2Norm;
          if ( dot(v1, v2) > 0 )
          {
            edge1Norm = vecSI.norm();
            edge2Norm = vecEI.norm();
          }
          else
          {
            edge1Norm = vecEI.norm();
            edge2Norm = vecSI.norm();
          }
          
          // if two lines becomes colinear, the boundary go to infinity
          // so we have to force the boundary to be in an edge or a line PQ.
          if ( t < edge1Norm && delta > t )
          {
            realStart = 0;
          }
          else if ( t >= edge1Norm && delta > edge1Norm )
          {
            realStart = t-edge1Norm;
          }
              
          if ( normPQ - t < edge2Norm && delta > normPQ - t)
          {
            realEnd = normPQ;
          }
          else if ( normPQ - t >= edge2Norm && delta > edge2Norm)
          {
            realEnd = t+edge2Norm;
          }
          Intersect[index].realStartPt.setTo(P, v1, realStart);
          Intersect[index].realEndPt.setTo(P, v1, realEnd);
          j++;
        }
        oldt = t;
      }
    }                   
  }
  /**
    @brief     This function returns the number of intersection points.
    @return		the number of intersection points.
*/ 
  int getNumOfIntersect()
  {
    int num = 0;
    for ( int i = 0; i < 2; i++ )
    {
      if ( Intersect[i].logicalPt.x != -1 || Intersect[i].logicalPt.y != -1 )
      {  
        num++;
      }
    }
    return num;      
  }
  
  intersect getIntersect(int i)
  { 
    return Intersect[i];
  }  
}  
// end of triangle class

/**
    @brief     This function decide whether two of lines are crossed or not.
    @see       <A href = "geometry.ppt">To clearly understand this function,
               please refer Slide 10 of this file.</A>
    @param A Start point of line 1
    @param B End point of line 1
    @param P Start point of line 2
    @param Q End point of line 2
    @return True if two of lines are crossed, false otherwise.
*/  
boolean isCross( pt A, pt B, pt P, pt Q )
{
  if( ( isLeftTurn( A, B, P ) != isLeftTurn( A, B, Q ) ) && ( isLeftTurn( P, Q, A ) != isLeftTurn( P, Q, B ) ) )
  {
    return true;
  }
  else
  {
    return false;  
  }
}

/**
    @brief     This function calculates the distance from the pt A to the crossed point.
    @see       <A href = "geometry.ppt">To clearly understand this function,
               please refer Slide 11 of this file.</A>
    @param A Start point of line 1
    @param B End point of line 1
    @param P Start point of line 2
    @param Q End point of line 2
    @return distance from S to t in float value.
*/
float findPoint( pt A, pt B, pt P, pt Q )
{
  //! First, we make a vector for line 1(AB)
  vec v1 = new vec();
  v1.setTo( A, B );
  v1.normalize();                //! normalize
  
  //! And then, we need another vector.. so make a vector for line 2(PQ)
  vec v2 = new vec();
  v2.setTo( P, Q );
  v2.turnLeft();                 //! turn left to find the N
  v2.normalize();                //! normalize
  
  //! We also need to make a temporary vector SQ
  vec SQ = new vec();
  SQ.setTo( A, Q );

  //! And calculate all together using Professor Jarek's lecture.
  float SQN = dot( SQ, v2 ); //! We calculate SQ.N
  float TN = dot( v1, v2 );  //! And T.N

  //! YAY! We got the "t" which is the distance from S.
  float t = SQN/TN;
  return t;
}


/**
    @brief     This function get the relative position of P on line AB.\n
               @see       <A href = "geometry.ppt">To clearly understand this function,
               please refer Slide 11 of this file.</A>
    @param A Point A
    @param B Point B
    @param P Point P
    @return -1000000 to indicate nonlinear when P is out of bound.
		@return the relative position of P when P in on AB
*/

float getRelativePosOnLine( pt A, pt B, pt P )
{
  // using AB line equation : A + t(B-A) ( 0 <= t <= 1 )
  if ( abs(B.x - A.x)  < tolerence && abs(B.y - A.y) < tolerence )
  {
    return nonlinear;
  }
  
  float t1, t2;  
	// AB is pararell to y-axis
  if ( abs(B.x - A.x) < tolerence)
  {
    if ( (P.x >= A.x && P.x <= B.x) || (P.x >= B.x && P.x <= A.x) ) 
    {
      t2 = (P.y - A.y ) / (B.y - A.y);
      t1 = t2;
    }
    else
    {
      return nonlinear;
    }
  }
	// AB is pararell to x-axis
  else if ( abs(B.y-A.y) < tolerence)
  {
    if ( (P.y >= A.y && P.y <= B.y) || (P.y >= B.y && P.y <= A.y) ) 
    {
      t1 = (P.x - A.x ) / (B.x - A.x);
      t2 = t1;
    }
    else
    {
      return nonlinear;
    }    
  }
  else 
  {
    t1 = (P.x - A.x ) / (B.x - A.x);
    t2 = (P.y - A.y ) / (B.y - A.y);
  }
  if ( abs(t1 - t2) < tolerence )
  {
    return t1;
  }
  else
  {
    // not linear
    return nonlinear;
  }
}

/**
    @brief   check whether P is on a line AB
*/
boolean isOnLine(pt A, pt B, pt P)
{
   float t = getRelativePosOnLine( A, B, P );
   //println("t="+t);
   
   if ( t >= 0 &&  t <= 1) return true;
   else return false;
}