/**
 * Cubic bezier curve class with all the fancies:
 *
 *   - Subdivision (de casteljau)
 *   - Parametric value computation
 *   - Bezier bounding box computation
 *       1) aligned to view coordinate system
 *       2) aligned to the curve's start-end line
 *
 */
class BezierCurve
{
  // when true, we show the result of applying de Casteljau's algorithm
  boolean showextras = true;
  void setShowExtras(boolean b) { showextras = b; }

  // when true, we only show the base curve, transparent
  boolean shadow = false;
  void setShadow(boolean b) {
    setShowExtras(!b);
    shadow = b; }

  // Main bezier parameters
  float start_x, start_y;
  float start_x_control, start_y_control;
  float end_x_control, end_y_control;
  float end_x, end_y;

  // getters for these values
  float getStartX() { return start_x - start_x_offset; }
  float getStartY() { return start_y - start_y_offset; }
  float getStartXControl() { return start_x_control - start_x_control_offset; }
  float getStartYControl() { return start_y_control - start_y_control_offset; }
  float getEndXControl() { return end_x_control - end_x_control_offset; }
  float getEndYControl() { return end_y_control - end_y_control_offset; }
  float getEndX() { return end_x - end_x_offset; }
  float getEndY() { return end_y - end_y_offset; }

  // Offsets that are used while repositioning points
  float start_x_offset = 0;
  float start_y_offset = 0;
  float start_x_control_offset = 0;
  float start_y_control_offset = 0;
  float end_x_control_offset = 0;
  float end_y_control_offset = 0;
  float end_x_offset = 0;
  float end_y_offset = 0;

  // Subdivision 't' value
  float sub_t = -1;
  float getSubT() { return sub_t; }

  // The "green" points
  float mid_x_1, mid_y_1;
  float mid_x_2, mid_y_2;
  float mid_x_3, mid_y_3;

  // The "yellow" points
  float mid_x_4, mid_y_4;
  float mid_x_5, mid_y_5;

  // The coordinates used in the split curves
  float prev_x_control, prev_y_control;
  float sub_x, sub_y;
  float next_x_control, next_y_control;

  // getters for these values
  float getPrevXControl() { return prev_x_control; }
  float getPrevYControl() { return prev_y_control; }
  float getSubX() { return sub_x; }
  float getSubY() { return sub_y; }
  float getNextXControl() { return next_x_control; }
  float getNextYControl() { return next_y_control; }

  // full curve bounding box
  float[] bounding_box;
  
  // bouding box modes
  int STANDARD = BoundingBoxRegulator.STANDARD;
  int STRAIGHTEN = BoundingBoxRegulator.STRAIGHTEN;
  int bounding_mode = BoundingBoxRegulator.boundingMode;
  void updateBoundingMode() { 
    bounding_mode = BoundingBoxRegulator.boundingMode;
    bounding_box = calculate_bbox(start_x,start_y,start_x_control,start_y_control,end_x_control,end_y_control,end_x,end_y); }

  // constants for using the bounding_box
  int MIN_X  = 0;
  int MIN_Y  = 1;
  int MAX_X = 2;
  int MAX_Y = 3;

  /**
   * Constructor - simply set up the initial curve and compute its bounding box
   */
  BezierCurve(float x1, float y1, float cx1, float cy1, float cx2, float cy2, float x2, float y2) {
    start_x = x1; start_y = y1;
    start_x_control = cx1; start_y_control = cy1;
    end_x_control = cx2; end_y_control = cy2;
    end_x = x2; end_y = y2;
    // compute the bounding box while we're at it.
    bounding_box = calculate_bbox(x1,y1,cx1,cy1,cx2,cy2,x2,y2); }

  /**
   * expand the x-bounds, if the value lies outside the bounding box
   */
  void exandXBounds(float[] bounds, float value) {
    if(bounds[MIN_X] > value) { bounds[MIN_X] = value; }
    else if(bounds[MAX_X] < value) { bounds[MAX_X] = value; }}

  /**
   * expand the x-bounds, if the value lies outside the bounding box
   */
  void exandYBounds(float[] bounds, float value) {
    if(bounds[MIN_Y] > value) { bounds[MIN_Y] = value; }
    else if(bounds[MAX_Y] < value) { bounds[MAX_Y] = value; }}

  /**
   * Calculate the bounding box. There are several "boxes"
   * that can be made, so this function's main purpose is to
   * check which one to make, and then delegate the job to
   * the function that knows how to perform the desired
   * calculations.
   */
  float[] calculate_bbox(float x1, float y1, float cx1, float cy1, float cx2, float cy2, float x2, float y2)
  {
    if(bounding_mode==STANDARD) {
      return calculate_standard_bbox(x1,y1,cx1,cy1,cx2,cy2,x2,y2); }
    else if(bounding_mode==STRAIGHTEN) {
      return calculate_straightened_bbox(x1,y1,cx1,cy1,cx2,cy2,x2,y2); }
    return null;
  }

  // specific angles used in straightening
  float d1 = 0.0;
  float d2 = 90.0;
  float d3 = 180.0;
  float d4 = 270.0;

  /**
   * rotate bezier so that {start->end is a horizontal} line,
   * then compute standard bbox, and counter-rotate it.
   */
  float[] calculate_straightened_bbox(float x1, float y1, float cx1, float cy1, float cx2, float cy2, float x2, float y2)
  {
    float angle = 0.0;
    float dx = x2-x1;
    float dy = y2-y1;
    //already horizontal!
    if(dy==0) { return calculate_standard_bbox(x1,y1,cx1,cy1,cx2,cy2,x2,y2); }
    
    // straightening is required!
    float adx = abs(dx);
    float ady = abs(dy);
    // special case: {start->end} is a vertical line
    if(dx==0) { angle = (dy>=0? d2 : d4); }
    // normal case: the end point is in one of four quadrants, relative to the starting point
    else if(dx>0 && dy>0) { angle = d1 + (atan(ady/adx) * (180.0/PI)); } // X+, Y+
    else if(dx<0 && dy<0) { angle = d3 + (atan(ady/adx) * (180.0/PI)); } // X-, Y-
    else if(dx<0 && dy>0) { angle = d2 + (atan(adx/ady) * (180.0/PI)); } // X-, Y+
    else if(dx>0 && dy<0) { angle = d4 + (atan(adx/ady) * (180.0/PI)); } // X+, Y-

    // Now that we have the angle, get the straightened bounding box.
    // We need to know the angle in radians for that, because that lets us
    // very quickly compute what the rotated coordinates should be:
    float phi = -(angle*PI/180.0);

    // translate by -x1/-y1, so the curve has 0/0 as origin
    cx1 -= x1; cy1 -= y1; cx2 -= x1; cy2 -= y1; x2 -= x1; y2 -= y1;
    
    // rotate over our angle
    float ncx1 = cx1*cos(phi) - cy1*sin(phi);
    float ncy1 = cx1*sin(phi) + cy1*cos(phi);
    float ncx2 = cx2*cos(phi) - cy2*sin(phi);
    float ncy2 = cx2*sin(phi) + cy2*cos(phi);
    float  nx2 =  x2*cos(phi) -  y2*sin(phi);
    float  ny2 =  x2*sin(phi) +  y2*cos(phi);

    // now we can compute the regular bounding box!
    float[] bounds = calculate_standard_bbox(0.0, 0.0, ncx1, ncy1, ncx2, ncy2, nx2, ny2);

    // if we rotate that straight bounding box back, and then move it to x1/y1 again,
    // we have our bounding box in the direction of {start->end}:
    phi = (angle*PI/180.0);
    float[] newbounds = { x1 + (bounds[0]*cos(phi) - bounds[1]*sin(phi)),
                          y1 + (bounds[0]*sin(phi) + bounds[1]*cos(phi)),
                          x1 + (bounds[2]*cos(phi) - bounds[3]*sin(phi)),
                          y1 + (bounds[2]*sin(phi) + bounds[3]*cos(phi)),
                          x1 + (bounds[4]*cos(phi) - bounds[5]*sin(phi)),
                          y1 + (bounds[4]*sin(phi) + bounds[5]*cos(phi)),
                          x1 + (bounds[6]*cos(phi) - bounds[7]*sin(phi)),
                          y1 + (bounds[6]*sin(phi) + bounds[7]*cos(phi))};
    return newbounds;
  }

  /**
   * Calculate the bezier value for one dimension at distance 't'
   */
  float calculate_bezier(float t, float p0, float p1, float p2, float p3) {
    float mt = (1-t);
    return (mt*mt*mt*p0) + (3*mt*mt*t*p1) + (3*mt*t*t*p2) + (t*t*t*p3); }
  
  /**
   * Calculate the bounding box for this bezier curve. The next bit is technical. See the comment on this topic on
   * http://newsgroups.derkeiler.com/Archive/Comp/comp.graphics.algorithms/2005-07/msg00334.html
   * and the worked out mathematics at http://pomax.nihongoresources.com/downloads/bezierbounds.html
   * for an explanation of why the following code is being used, and why it works.
   */
  float[] calculate_standard_bbox(float x1, float y1, float cx1, float cy1, float cx2, float cy2, float x2, float y2)
  {
    // compute linear bounds first
    float[] bounds = {min(x1,x2), min(y1,y2), max(x1,x2), max(y1,y2)};
  
    float dcx0 = cx1 - x1;
    float dcy0 = cy1  -y1;
    float dcx1 = cx2 - cx1;
    float dcy1 = cy2 - cy1;
    float dcx2 = x2 - cx2;
    float dcy2 = y2 - cy2;
  
    // Recompute bounds projected on the x-axis, if the control points lie outside the bounding box x-bounds
    if(cx1<bounds[MIN_X] || cx1>bounds[MAX_X] || cx2<bounds[MIN_X] || cx2>bounds[MAX_X]) {
      // we don't need to do this, but 'a', 'b' and 'c' are easier to read.
      float a = dcx0; float b = dcx1; float c = dcx2;
  
      // Do we have a problematic discriminator if we use these values?
      // If we do, because we're computing at sub-pixel level anyway, simply salt 'b' a tiny bit.
      if(a+c != 2*b) { b+=0.01; }
  
      float numerator = 2*(a - b);
      float denominator = 2*(a - 2*b + c);
      float quadroot = (2*b-2*a)*(2*b-2*a) - 2*a*denominator;
      float root = sqrt(quadroot);
  
      // there are two possible values for 't' that yield inflection points,
      // and each of these inflection points might be outside the linear bounds
      float t1 =  (numerator + root) / denominator;
      float t2 =  (numerator - root) / denominator;
  
      // so, which of these is the useful point? (remember that t must lie
      // in [0,1], and that t=0 and t=1 are already the linear bounds)
      if(0<t1 && t1<1) { exandXBounds(bounds, calculate_bezier(t1, x1, cx1, cx2, x2)); }
      if(0<t2 && t2<1) { exandXBounds(bounds, calculate_bezier(t2, x1, cx1, cx2, x2)); }
    }

    // Recompute bounds projected on the y-axis, if the control points lie outside the bounding box
    // y-bounds. No comments this time, because the code is the same, just with 'y' instead of 'x'
    if(cy1<bounds[MIN_Y] || cy1>bounds[MAX_Y] || cy2<bounds[MIN_Y] || cy2>bounds[MAX_Y]) {
      float a = dcy0; float b = dcy1; float c = dcy2;
      if(a+c != 2*b) { b+=0.01; }
      float numerator = 2*(a - b);
      float denominator = 2*(a - 2*b + c);
      float quadroot = (2*b-2*a)*(2*b-2*a) - 2*a*denominator;
      float root = sqrt(quadroot);
      float t1 =  (numerator + root) / denominator;
      float t2 =  (numerator - root) / denominator;
      if(0<t1 && t1<1) { exandYBounds(bounds, calculate_bezier(t1, y1, cy1, cy2, y2)); }
      if(0<t2 && t2<1) { exandYBounds(bounds, calculate_bezier(t2, y1, cy1, cy2, y2)); }
    }

    // and we're done, form this box's rectangle and return
    float[] bbox = {bounds[0],bounds[1], bounds[0],bounds[3], bounds[2],bounds[3], bounds[2],bounds[1] };
    return bbox;
  }

  /**
   * Split our curve at 't', using de Casteljau's algorithm.
   */
  BezierCurve[] subdivide(float t)
  {
    // --- 1st division ---

    sub_t = t;
    // start -> start_control
    float dx = (start_x_control-start_x_control_offset) - (start_x-start_x_offset);
    float dy = (start_y_control-start_y_control_offset) - (start_y-start_y_offset);
    // get first midpoint
    mid_x_1 = (start_x-start_x_offset) + t*dx;
    mid_y_1 = (start_y-start_y_offset) + t*dy;
    // start_control -> end_control
    dx = (end_x_control-end_x_control_offset) - (start_x_control-start_x_control_offset);
    dy = (end_y_control-end_y_control_offset) - (start_y_control-start_y_control_offset);
    // get second midpoint
    mid_x_2 = (start_x_control-start_x_control_offset) + t*dx;
    mid_y_2 = (start_y_control-start_y_control_offset) + t*dy;
    // end_control -> end
    dx = (end_x-end_x_offset) - (end_x_control-end_x_control_offset);
    dy = (end_y-end_y_offset) - (end_y_control-end_y_control_offset);
    // get third midpoint
    mid_x_3 = (end_x_control-end_x_control_offset) + t*dx;
    mid_y_3 = (end_y_control-end_y_control_offset) + t*dy;

    // --- 2nd division ---

    // first midpoint -> second midpoint
    dx = mid_x_2 - mid_x_1;
    dy = mid_y_2 - mid_y_1;
    // get first tangent line point
    mid_x_4 = mid_x_1 + t*dx;
    mid_y_4 = mid_y_1 + t*dy;
    // second midpoint -> third midpoint
    dx = mid_x_3 - mid_x_2;
    dy = mid_y_3 - mid_y_2;
    // get second tangent line point
    mid_x_5 = mid_x_2 + t*dx;
    mid_y_5 = mid_y_2 + t*dy;

    // --- final values ---

    // first tangent midpoint -> second tangent midpoint
    dx = mid_x_5 - mid_x_4;
    dy = mid_y_5 - mid_y_4;
    // sub-curve parameters
    prev_x_control = mid_x_4;
    prev_y_control = mid_y_4;
    sub_x = mid_x_4 + t*dx;
    sub_y = mid_y_4 + t*dy;
    next_x_control = mid_x_5;
    next_y_control = mid_y_5;
    
    // generate the two sub-curves
    BezierCurve[] curves = new BezierCurve[2];
    curves[0] = new BezierCurve(start_x-start_x_offset, start_y-start_y_offset, mid_x_1, mid_y_1, prev_x_control, prev_y_control, sub_x, sub_y);
    curves[1] = new BezierCurve(sub_x, sub_y, next_x_control, next_y_control, mid_x_3, mid_y_3, end_x-end_x_offset, end_y-end_y_offset);
    return curves;
  }

  /**
   * Draw this bezier curve
   */
  void draw() { drawOffset(0,0); }

  /**
   * Draw this bezier curve, offset by some x/y value. We also draw
   * the most recent splitting parameters, if 'showextras' is true.
   */
  void drawOffset(int xoff, int yoff)
  {
    noFill();
    // bounding box
    if(showextras && bounding_box!=null) {
      stroke(200,0,0,100);
      line(xoff+bounding_box[0], yoff+bounding_box[1], xoff+bounding_box[2], yoff+bounding_box[3]);
      line(xoff+bounding_box[2], yoff+bounding_box[3], xoff+bounding_box[4], yoff+bounding_box[5]);
      line(xoff+bounding_box[4], yoff+bounding_box[5], xoff+bounding_box[6], yoff+bounding_box[7]);
      line(xoff+bounding_box[6], yoff+bounding_box[7], xoff+bounding_box[0], yoff+bounding_box[1]); }

    // the bezier curve
    if(shadow) { stroke(0,0,127,60); } else { stroke(0); }
    bezier(	xoff+start_x-start_x_offset,
		yoff+start_y-start_y_offset,
		xoff+start_x_control-start_x_control_offset,
		yoff+start_y_control-start_y_control_offset,
		xoff+end_x_control-end_x_control_offset,
		yoff+end_y_control-end_y_control_offset,
		xoff+end_x-end_x_offset,
		yoff+end_y-end_y_offset);

    // curve parameters and control box
    if(showextras) {
      // control box
      stroke(0,0,255,100);
      line(	xoff+start_x-start_x_offset,
		yoff+start_y-start_y_offset,
		xoff+start_x_control-start_x_control_offset,
		yoff+start_y_control-start_y_control_offset);
      line(	xoff+end_x_control-end_x_control_offset,
		yoff+end_y_control-end_y_control_offset,
		xoff+end_x-end_x_offset,
		yoff+end_y-end_y_offset);
      line(	xoff+start_x_control-start_x_control_offset,
		yoff+start_y_control-start_y_control_offset,
		xoff+end_x_control-end_x_control_offset,
		yoff+end_y_control-end_y_control_offset);
      // end points
      stroke(0,0,255);
      fill(0,0,255);
      ellipse(	xoff+start_x-start_x_offset,
		yoff+start_y-start_y_offset,
		5, 5);
      ellipse(	xoff+end_x-end_x_offset,
		yoff+end_y-end_y_offset,
		5, 5);
      // control points
      stroke(255,0,0);
      fill(255,0,0);
      ellipse(	xoff+start_x_control-start_x_control_offset,
		yoff+start_y_control-start_y_control_offset,
		5, 5);
      ellipse(	xoff+end_x_control-end_x_control_offset,
		yoff+end_y_control-end_y_control_offset,
		5, 5);
      // de Casteljau's parameters
      if(sub_t>-1)
      {
        // "green" control box
        stroke(0,255,0,150);
        line(xoff+mid_x_1, yoff+mid_y_1, xoff+mid_x_2, yoff+mid_y_2);
        line(xoff+mid_x_2, yoff+mid_y_2, xoff+mid_x_3, yoff+mid_y_3);
        // "yellow" control box
        stroke(255,255,0);
        line(xoff+mid_x_4, yoff+mid_y_4, xoff+mid_x_5, yoff+mid_y_5);
        // "green" control points
        stroke(0,255,0);
        fill(0,255,0);
        ellipse(xoff+mid_x_1, yoff+mid_y_1, 5, 5);
        ellipse(xoff+mid_x_2, yoff+mid_y_2, 5, 5);
        ellipse(xoff+mid_x_3, yoff+mid_y_3, 5, 5);
        // "yellow" control points
        stroke(255,255,0);
        fill(255,255,0);
        ellipse(xoff+mid_x_4, yoff+mid_y_4, 5, 5);
        ellipse(xoff+mid_x_5, yoff+mid_y_5, 5, 5);
        // "black" point
        stroke(0);
        fill(0);
        ellipse(xoff+sub_x, yoff+sub_y, 5, 5);
      }
    }
  }

  // -----------------------
  // MOUSE EVENT HANLDING
  // -----------------------

  int activepoint = -1;
  int START = 1;
  int START_CONTROL = 2;
  int END_CONTROL = 3;
  int END = 4;
  float xmark = -1;
  float ymark = -1;

  /**
   * is this bezier curve being manipulated?
   */
  boolean isActive() { return activepoint>0; }

  /**
   * Signal that the mouse is over something interesting when
   * it is above one of the four points on the bezier curve.
   */
  boolean mouseOver(float mx, float my) {
    return  (abs(start_x-mx)<5 && abs(start_y-my)<5) ||
            (abs(start_x_control-mx)<5 && abs(start_y_control-my)<5) ||
            (abs(end_x_control-mx)<5 && abs(end_y_control-my)<5) ||
            (abs(end_x-mx)<5 && abs(end_y-my)<5); }

  /**
   * If the mouse is over one of the four points, and it is clicked, 
   * register on which point that occurred, and set that point to "active".
   */
  boolean mousePressed(float mx, float my) {
    if(abs(start_x-mx)<5 && abs(start_y-my)<5) { activepoint = START; }
    else if(abs(start_x_control-mx)<5 && abs(start_y_control-my)<5) { activepoint = START_CONTROL; }
    else if(abs(end_x_control-mx)<5 && abs(end_y_control-my)<5) { activepoint = END_CONTROL; }
    else if(abs(end_x-mx)<5 && abs(end_y-my)<5) { activepoint = END; }
    boolean active = isActive();
    if(active) {
      xmark=mx;
      ymark=my; }
    return active;
  }

  /**
   * When the mouse is dragged, we record the x/y offsets the mouse
   * is describing for the currently active point.
   */
  void mouseDragged(float mx, float my) {
    float xdiff = xmark - mx;
    float ydiff = ymark - my;
    if(isActive()) {
      if(activepoint==START) { start_x_offset = xdiff; start_y_offset = ydiff; }
      else if(activepoint==START_CONTROL) { start_x_control_offset = xdiff; start_y_control_offset = ydiff; }
      else if(activepoint==END_CONTROL) { end_x_control_offset = xdiff; end_y_control_offset = ydiff; }
      else if(activepoint==END) { end_x_offset = xdiff; end_y_offset = ydiff; }
      // recompute bounding box
      bounding_box = calculate_bbox((start_x - start_x_offset),
                                    (start_y - start_y_offset),
                                    (start_x_control - start_x_control_offset),
                                    (start_y_control - start_y_control_offset),
                                    (end_x_control - end_x_control_offset),
                                    (end_y_control - end_y_control_offset),
                                    (end_x - end_x_offset),
                                    (end_y - end_y_offset));
      // recompute subdivisions, if this curve was split at some 't' value
      if(sub_t>-1) { subdivide(sub_t); }}
  }

  /**
   * When the mouse is released, we commit the current active offset
   * to the current active point, then set everything back to inactive.
   */
  boolean mouseReleased(float mx, float my)
  {
    boolean handled = false;
    float xdiff = xmark - mx;
    float ydiff = ymark - my;
    if(activepoint==START) { start_x -= xdiff; start_y -= ydiff; }
    else if(activepoint==START_CONTROL) { start_x_control -= xdiff; start_y_control -= ydiff; }
    else if(activepoint==END_CONTROL) { end_x_control -= xdiff; end_y_control -= ydiff; }
    else if(activepoint==END) { end_x -= xdiff; end_y -= ydiff; }
    if(isActive()) {
      handled = true;
      bounding_box = calculate_bbox(start_x,start_y,start_x_control,start_y_control,end_x_control,end_y_control,end_x,end_y); }
    xmark = -1;
    ymark = -1;
    activepoint = -1;
    start_x_offset = 0;
    start_y_offset = 0;
    start_x_control_offset = 0;
    start_y_control_offset = 0;
    end_x_control_offset = 0;
    end_y_control_offset = 0;
    end_x_offset = 0;
    end_y_offset = 0;
    return handled;
  }
}