/**
 * @classdesc A HobbyCurve/HobbyPath calculation class: compute a set of optimal
 *            cubic Bézier curves from a sequence of vertices.
 *
 * This Hobby curve (path) implementation was strongly inspired by
 * the one by Prof. Dr. Edmund Weitz:
 * Here's the website:
 *  http://weitz.de/hobby/
 *
 * @date     2020-04-07
 * @author   Transformed to a JS-class by Ikaros Kappler
 * @modified 2020-08-19 Ported from vanilla JS to TypeScript.
 * @version 1.0.1
 *
 * @file HobbyPath
 * @public
 **/

import { CubicBezierCurve } from "../../CubicBezierCurve";
import { Vertex } from "../../Vertex";

interface IControlPoints {
  startControlPoints: Array<Vertex>;
  endControlPoints: Array<Vertex>;
}

/**
 * @classdesc A HobbyCurve/HobbyPath calculation class: compute a set of optimal
 *            cubic Bézier curves from a sequence of vertices.
 *
 * @requires CubicBezierCurve
 * @requires Vertex
 */
export class HobbyPath {
  /**
   * @member {Array<Vertex>} vertices
   * @memberof HobbyPath
   * @type {Array<Vertex>}
   * @instance
   **/
  private vertices: Array<Vertex>;

  /**
   * @constructor
   * @name HobbyPath
   * @param {Array<Vertex>=} vertices? - An optional array of vertices to initialize the path with.
   **/
  constructor(vertices?: Array<Vertex>) {
    this.vertices = vertices ? vertices : [];
  }

  /**
   * Add a new point to the end of the vertex sequence.
   *
   * @name addPoint
   * @memberof HobbyPath
   * @instance
   * @param {Vertex} p - The vertex (point) to add.
   **/
  addPoint(p: Vertex): void {
    this.vertices.push(p);
  }

  /**
   * Generate a sequence of cubic Bézier curves from the point set.
   *
   * @name generateCurve
   * @memberof HobbyPath
   * @instance
   * @param {boolean=} circular - Specify if the path should be closed.
   * @param {number=0} omega - (default=0) An optional tension parameter.
   * @return Array<CubicBezierCurve>
   **/
  generateCurve(circular?: boolean, omega?: number): Array<CubicBezierCurve> {
    let n: number = this.vertices.length;
    if (n > 1) {
      if (n == 2) {
        // for two points, just draw a straight line
        return [new CubicBezierCurve(this.vertices[0], this.vertices[1], this.vertices[0], this.vertices[1])];
      } else {
        const curves: Array<CubicBezierCurve> = [];
        let controlPoints = this.hobbyControls(circular, omega);
        for (let i = 0; i < n - (circular ? 0 : 1); i++) {
          // if i is n-1, the "next" point is the first one
          let j: number = (i + 1) % n; // Use a succ function here?
          curves.push(
            new CubicBezierCurve(
              this.vertices[i],
              this.vertices[j],
              controlPoints.startControlPoints[i],
              controlPoints.endControlPoints[i]
            )
          );
        }
        return curves;
      }
    } else {
      return [];
    }
  }

  /**
   * Computes the control point coordinates for a Hobby curve through
   * the points given.
   *
   * @name hobbyControls
   * @memberof HobbyPath
   * @instance
   * @param {boolean}  circular - If true, then the path will be closed.
   * @param {number=0} omega    - The 'curl' or the path.
   * @return {IControlPoints} An object with two members: startControlPoints and endControlPoints (Array<Vertex>).
   **/
  private hobbyControls(circular?: boolean, omega?: number): IControlPoints {
    // This is a version that works for both, closed and non-closed paths.

    if (typeof omega === "undefined") omega = 0;

    let n: number = this.vertices.length - (circular ? 0 : 1);
    let D: Array<number> = new Array<number>(n);

    let ds: Array<Vertex> = new Array<Vertex>(n);
    var succ = (i: number): number => {
      return circular ? (i + 1) % n : i + 1;
    };
    var pred = (i: number): number => {
      return circular ? (i + n - 1) % n : i - 1;
    };

    for (let i = 0; i < n; i++) {
      // the "next" point in a modular way
      let j: number = succ(i);
      ds[i] = this.vertices[i].difference(this.vertices[j]);
      D[i] = Math.sqrt(ds[i].x * ds[i].x + ds[i].y * ds[i].y);
    }
    let gamma: Array<number> = new Array(n + (circular ? 0 : 1));
    for (let i = circular ? 0 : 1; i < n; i++) {
      // the "previous" point in a modular way
      let k: number = pred(i);
      let sin = ds[k].y / D[k];
      let cos = ds[k].x / D[k];
      let vec = HobbyPath.utils.rotate(ds[i], -sin, cos);
      gamma[i] = Math.atan2(vec.y, vec.x);
    }
    if (!circular) gamma[n] = 0;
    let a: Array<number> = new Array(n + (circular ? 0 : 1));
    let b: Array<number> = new Array(n + (circular ? 0 : 1));
    let c: Array<number> = new Array(n + (circular ? 0 : 1));
    let d: Array<number> = new Array(n + (circular ? 0 : 1));
    for (let i = circular ? 0 : 1; i < n; i++) {
      // j is the "next" point, k the "previous" one
      let j: number = succ(i);
      let k: number = pred(i);
      // see video for the equations
      a[i] = 1 / D[k];
      b[i] = (2 * D[k] + 2 * D[i]) / (D[k] * D[i]);
      c[i] = 1 / D[i];
      d[i] = -(2 * gamma[i] * D[i] + gamma[j] * D[k]) / (D[k] * D[i]);
    }
    // make matrix tridiagonal in preparation for the "sherman" function
    var alpha: Array<number>;
    var beta: Array<number>;
    if (circular) {
      let s: number = a[0] * omega; // Use omega here?
      a[0] = 0;
      let t: number = c[n - 1] * omega; // Use omega here?
      c[n - 1] = 0;
      alpha = HobbyPath.utils.sherman(a, b, c, d, s, t);
      beta = new Array<number>(n);
      for (let i = 0; i < n - (circular ? 0 : 1); i++) {
        // "next" point
        let j: number = succ(i);
        beta[i] = -gamma[j] - alpha[j];
      }
    } else {
      // see the Jackowski article for the following values; the result
      // will be that the curvature at the first point is identical to the
      // curvature at the second point (and likewise for the last and
      // second-to-last)
      b[0] = 2 + omega;
      c[0] = 2 * omega + 1;
      d[0] = -c[0] * gamma[1];
      a[n] = 2 * omega + 1;
      b[n] = 2 + omega;
      d[n] = 0;
      // solve system for the angles called "alpha" in the video
      alpha = HobbyPath.utils.thomas(a, b, c, d);
      // compute "beta" angles from "alpha" angles
      beta = new Array<number>(n);
      for (let i = 0; i < n - 1; i++) beta[i] = -gamma[i + 1] - alpha[i + 1];
      // again, see Jackowski article
      beta[n - 1] = -alpha[n];
    }
    let startControlPoints = new Array<Vertex>(n);
    let endControlPoints = new Array<Vertex>(n);
    for (let i = 0; i < n; i++) {
      let j: number = succ(i);
      let a: number = (HobbyPath.utils.rho(alpha[i], beta[i]) * D[i]) / 3;
      let b: number = (HobbyPath.utils.rho(beta[i], alpha[i]) * D[i]) / 3;
      let v: Vertex = HobbyPath.utils.normalize(HobbyPath.utils.rotateAngle(ds[i], alpha[i]));
      startControlPoints[i] = new Vertex(this.vertices[i].x + a * v.x, this.vertices[i].y + a * v.y);
      v = HobbyPath.utils.normalize(HobbyPath.utils.rotateAngle(ds[i], -beta[i]));
      endControlPoints[i] = new Vertex(this.vertices[j].x - b * v.x, this.vertices[j].y - b * v.y);
    }
    return { startControlPoints: startControlPoints, endControlPoints: endControlPoints };
  }

  static utils = {
    // rotates a vector [x, y] about an angle; the angle is implicitly
    // determined by its sine and cosine
    rotate: (vert: Vertex, sin: number, cos: number): Vertex => {
      return new Vertex(vert.x * cos - vert.y * sin, vert.x * sin + vert.y * cos);
    },
    // rotates a vector [x, y] about the angle alpha
    rotateAngle: (vert: Vertex, alpha: number): Vertex => {
      return HobbyPath.utils.rotate(vert, Math.sin(alpha), Math.cos(alpha));
    },
    // returns a normalized version of the vector
    normalize: (vec: Vertex): Vertex => {
      let n = Math.hypot(vec.x, vec.y);
      if (n == 0) return new Vertex(0, 0);
      else return new Vertex(vec.x / n, vec.y / n); // TODO: do in-place
    },
    // the "velocity function" (also called rho in the video); a and b are
    // the angles alpha and beta, the return value is the distance between
    // a control point and its neighboring point; to compute sigma(a,b)
    // we'll simply use rho(b,a)
    rho: (a: number, b: number): number => {
      // see video for formula
      let sa = Math.sin(a);
      let sb = Math.sin(b);
      let ca = Math.cos(a);
      let cb = Math.cos(b);
      let s5 = Math.sqrt(5);
      let num = 4 + Math.sqrt(8) * (sa - sb / 16) * (sb - sa / 16) * (ca - cb);
      let den = 2 + (s5 - 1) * ca + (3 - s5) * cb;
      return num / den;
    },
    // Implements the Thomas algorithm for a tridiagonal system with i-th
    // row a[i]x[i-1] + b[i]x[i] + c[i]x[i+1] = d[i] starting with row
    // i=0, ending with row i=n-1 and with a[0] = c[n-1] = 0.  Returns the
    // values x[i] as an array.
    thomas: (a: Array<number>, b: Array<number>, c: Array<number>, d: Array<number>): Array<number> => {
      let n: number = a.length;
      let cc: Array<number> = new Array<number>(n);
      let dd: Array<number> = new Array<number>(n);
      // forward sweep
      cc[0] = c[0] / b[0];
      dd[0] = d[0] / b[0];
      for (let i = 1; i < n; i++) {
        let den: number = b[i] - cc[i - 1] * a[i];
        cc[i] = c[i] / den;
        dd[i] = (d[i] - dd[i - 1] * a[i]) / den;
      }
      let x: Array<number> = new Array<number>(n);
      // back substitution
      x[n - 1] = dd[n - 1];
      for (let i = n - 2; i >= 0; i--) x[i] = dd[i] - cc[i] * x[i + 1];
      return x;
    },
    // Solves an "almost" tridiagonal linear system with i-th row
    // a[i]x[i-1] + b[i]x[i] + c[i]x[i+1] = d[i] starting with row i=0,
    // ending with row i=n-1 and with a[0] = c[n-1] = 0.  Returns the
    // values x[i] as an array.  The system is not really tridiagonal
    // because the 0-th row is b[0]x[0] + c[0]x[1] + sx[n-1] = d[0] and
    // row n-1 is tx[0] + a[n-1]x[n-2] + b[n-1]x[n-1] = d[n-1].  The
    // Sherman-Morrison-Woodbury formula is used so that the function
    // "thomas" can be called to solve the system.
    sherman: (a: Array<number>, b: Array<number>, c: Array<number>, d: Array<number>, s: number, t: number): Array<number> => {
      const n: number = a.length;
      const u: Array<number> = new Array<number>(n);
      u.fill(0, 1, n - 1);
      u[0] = 1;
      u[n - 1] = 1;
      let v: Array<number> = new Array<number>(n);
      v.fill(0, 1, n - 1);
      v[0] = t;
      v[n - 1] = s;
      b[0] -= t;
      b[n - 1] -= s;
      // this would be more efficient if computed in parallel, but hey...
      const Td: Array<number> = HobbyPath.utils.thomas(a, b, c, d);
      const Tu: Array<number> = HobbyPath.utils.thomas(a, b, c, u);
      const factor: number = (t * Td[0] + s * Td[n - 1]) / (1 + t * Tu[0] + s * Tu[n - 1]);
      const x = new Array(n);
      for (let i = 0; i < n; i++) x[i] = Td[i] - factor * Tu[i];
      return x;
    }
  };
} // END class