// Written by Jarek Rossignac, Georgia Tech, Sept 2005.
int count =5;                    // number of control vertices
pt[] P = new pt [count];         // control vertices 
int bi=-1;                       // index of selected vertex, -1 if none selected
pt Mouse = new pt(0,0);          // current mouse position

color red = color(200, 10, 10); color blue = color(10, 10, 200); color green = color(10, 200, 20); 
color black = color(10, 10, 10); color magenta = color(250, 150, 200); 

vec a,b,c;                      // cubic coeffs

float tt;                      // parameter of the point C on the cubic
boolean straight = false;      // toggled with "c" for making the fit independent of C
pt BF;

void setup() 
{ 
  size(600, 600);   
  for (int i=0; i<count; i++) {P[i]=new pt(-sin((i+0.5)*TWO_PI/count)*(250)+300,cos((i+0.5)*TWO_PI/count)*(250)+300);};
  strokeJoin(ROUND); strokeCap(ROUND);
  PFont font = loadFont("Courier-14.vlw"); textFont(font, 12);
} 
 
void draw() { 
  background(255); 
  if (bi!= -1) {P[bi].setToMouse();};                    // snap selected vertex to mouse position during dragging
  stroke(black); strokeWeight(1); drawControlCurve();
  stroke(magenta); strokeWeight(2);   tt=CubicFit(P[0],P[1],P[2],P[3],P[4]);  DrawCubic(); 
  stroke(green); strokeWeight(5.0);  DrawCorrection(tt);
  stroke(red); strokeWeight(2.0); BF = betterFit(P[0],P[1],P[2],P[3],P[4], sqrt(2.0/3.0)); BF.showLineTo(P[2]); 
  };
  
//**************************************
//**** BETTER FIT OF PROPORTIONAL BI-LAPLACIAN
//**************************************
//*** Proportional fit    
pt betterFit(pt A, pt B, pt C, pt D, pt E, float s) { 
   pt P, PPC;
     if (straight) {                                          // proportional bi laplacian
     P = midPt(A,C);   pt PB = B.make(); PB.moveTowards(P,s);
     P = midPt(B,D);   pt PC = C.make(); PC.moveTowards(P,s);
     P = midPt(C,E);   pt PD = D.make(); PD.moveTowards(P,s);
     P = midPt(PB,PD);   PPC = PC.make();   PPC.moveTowards(P,-s);
     }
    else {                                                  // standard bi laplacian
     P = B.make().propFit(A,C);   pt PB = B.make(); PB.moveTowards(P,s);
     P = C.make().propFit(B,D);   pt PC = C.make(); PC.moveTowards(P,s);
     P = D.make().propFit(C,E);   pt PD = D.make(); PD.moveTowards(P,s);
     P = PC.make().propFit(PB,PD);  PPC = PC.make();   PPC.moveTowards(P,-s);
   };
  return(PPC);
  };
  
//**************************************
//**** PROPORTIONAL CUBIC FIT
//**************************************
//*** Fits cubic through A, B, D, E using parameters that are proportional to edge lengths of AB, CB, CD, DE   
float CubicFit(pt A, pt B, pt C, pt D, pt E) {   // to compute using 4 cproportional onstraints ***
    float s, t, u, ct1, ct2, ct3, l;
    vec AB, AD, AE;
    if (straight) {                              // use |BD| for |BC| + |CD|
      l = A.vecTo(B).norm() + B.vecTo(D).norm() + D.vecTo(E).norm();
      s = A.vecTo(B).norm() / l;   u = s + B.vecTo(D).norm() / l; 
      t =  s + B.vecTo(C).norm() / (B.vecTo(C).norm()+C.vecTo(D).norm())  * B.vecTo(D).norm() / l;     
    }
    else {                                      // use |BC| + |CD|, so the cubic depends on C
      l = A.vecTo(B).norm() + B.vecTo(C).norm() + C.vecTo(D).norm() + D.vecTo(E).norm();
      s = A.vecTo(B).norm() / l;   t = s + B.vecTo(C).norm() / l;  u = t + C.vecTo(D).norm() / l;    
      };
    ct1 = (s*u + s - u - s*s) * s;     ct2 = (u*u + s - s*u - u) * u;     ct3 = s*u + 1 - s - u;
    AB = A.vecTo(B); AB.div(ct1);      AD = A.vecTo(D); AD.div(ct2);       AE = A.vecTo(E); AE.div(ct3);

    a = AB.make(); a.mul(-1.0); a.add(AD); a.add(AE);   
    b = AB.make(); b.mul(u+1); b.addScaled(-(s+1),AD); b.addScaled(-(u+s),AE);
    c = AB.make(); c.mul(-u);  c.addScaled(s,AD); c.addScaled(u*s,AE); 

      return(t);
   };
   
void DrawCubic() { for (float t=0.0; t<1.0; t+=0.01) { evaluateCubic(t).show(1); };   };

void DrawCorrection(float t) {pt CP = evaluateCubic(t); CP.showLineTo(P[2]); };   
   
pt evaluateCubic(float t) {
   vec V = a.make(); V.mul(t*t*t); V.addScaled(t*t,b); V.addScaled(t,c);
   pt R = P[0].make();  R.addVec(V);
   return (R);
   };

void drawControlCurve() {       
       fill(220); beginShape();  for (int i=0; i<count; i++) {P[i].vert(); P[i].show(2);}; P[0].vert(); endShape(); noFill();
       for (int i=0; i<count; i++) {vec V = P[ip(i)].vecTo(P[in(i)]); V.unit(); V.left(); V.mul(15); P[i].label(str(i),V );};
        };    
       

//**************************************
// **** GUI FOR EDITING THE CONTROL CURVE
//**************************************
void mousePressed() {                                                   // to do when mouse is pressed  
  float bd=900;                                                     // init square of smallest distance to selected point
  Mouse.setToMouse();                                                   // save current mouse location
   for (int i=0; i<count; i++) { if (d2(i)<bd) {bd=d2(i); bi=i;  };};       // select closeest vertex
  if (bd>10) {                                                         // if closest vertex is too far
      bd=600;                                                       // reinitilize distance squared
      for (int i=0; i<count; i++) {if (dm(i)<bd) {bd=dm(i); bi=i;};};   // closest mid-edge point
      if (bd<20) {  
          for (int i=count-1; i>bi; i--) {P[i+1].setTo(P[i]);};                // shift down the rest 
          bi++;  P[bi].setTo(Mouse); count++;                                  // insert new vertex at mouse position
          }
        else {bi=-1;};                                                 // nothing selected
      };
  }

float d2(int j) {return (Mouse.disTo(P[j]));};   //  squared distance from mouse to vertex P[i]
float dm(int j) {return (Mouse.vecToMid(P[j],P[in(j)]).norm());};

void mouseReleased() {                                      // do this when mouse released
  if ( (bi!=-1) &&  P[bi].isInWindow() ) {                  // if outside of port
      for (int i=bi; i<count; i++) {P[i].setTo(P[in(i)]);};        // shift up to delete selected vertex
      count--; };                                           // reduce vertex count
  bi= -1;   
  };
 

//**********************************
//***      CONTROL WITH KEYS
//**********************************
void keyPressed() {  
  if (key=='c') {straight=!straight;};
  if (key=='w') {writePts();};
  if (key=='s') {savePts();};
  if (key=='l') {loadPts();};
  if (key=='a') {                          // axis aligns the edges
      for (int i=0; i<count; i++) {
         if ( abs(P[in(i)].x-P[i].x)>abs(P[in(i)].y-P[i].y)) {float dy=(P[in(i)].y-P[i].y)/2; P[in(i)].y-=dy; P[i].y+=dy;}
          else {float dx=(P[in(i)].x-P[i].x)/2; P[in(i)].x-=dx; P[i].x+=dx;}
         };
     };
    };                                                                 

//**************************************
//**** utilities for polyloops
//**************************************
int in(int j) {  if (j==count-1) {return (0);}  else {return(j+1);}  };  // next vertex in control loop
int ip(int j) {  if (j==0) {return (count-1);}  else {return(j-1);}  };  // next vertex in control loop

void writePts() { for (int i=0; i<count; i++) {P[i].write();};};

void savePts() {
  String [] inppts = new String [count+1];
  int s=0;
  inppts[s++]=str(count);
  for (int i=0; i<count; i++) {inppts[s++]=str(P[i].x)+","+str(P[i].y);};
  saveStrings("loop.txt",inppts);
  };

void loadPts() {
  String [] ss = loadStrings("loop.txt");
  String subpts;
  int s=0;
  int comma;
  count = int(ss[s]);
  for(int i=0;i<count; i++) {
      comma=ss[++s].indexOf(',');
      P[i]=new pt (float(ss[s].substring(0, comma)), float(ss[s].substring(comma+1, ss[s].length())));
    };
  };