
// cubicInflectionPoint.hpp - templated routine for determining inflections points of Bezier cubic curve

#ifndef __cubicInflectionPoint_hpp__
#define __cubicInflectionPoint_hpp__

#pragma once

#include "cubic_solve.hpp"

template <typename Tvec>
static int cubicInflectionPoint(const Tvec cp[4], double key_point[])
{
    // Convert control points to cubic monomial polynomial coefficients
    const Tvec A = cp[3] - cp[0] + (cp[1] - cp[2]) * Tvec(3),
               B = (cp[0] - cp[1] * Tvec(2) + cp[2]) * Tvec(3),
               C = (cp[1] - cp[0]) * Tvec(3),
               D = cp[0];

    // "Inflection points of a cubic Bezier"
    // http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html
    // 0 = cross(P'(t),P"(t)) = Px*Pyy-Py*Pxx
    //
    // Developed in Maple:
    // 
    // with(linalg):
    // > unprotect(D);
    // > A:=vector([Ax,Ay,0]):
    // > B:=vector([Bx,By,0]):
    // > C:=vector([Cx,Cy,0]):
    // > D:=vector([Dx,Dy,0]):
    // > P:=A*t^3+B*t^2+C*t+D;
    // P := A*t^3+B*t^2+C*t+D
    // > Pt:=(diff(P,t));Ptt:=(diff(P,t,t));
    // Pt := 3*A*t^2+2*B*t+C
    // Ptt := 6*A*t+2*B
    // > v:=crossprod(Pt,Ptt);
    // v := vector([0, 0, (3*t^2*Ax+2*t*Bx+Cx)*(6*t*Ay+2*By)-(3*t^2*Ay+2*t*By+Cy)*(6*t*Ax+2*Bx)])
    // > collect(v[3]/2,t);  % dividing by 2 just scales coefficients without changing polynomial's roots
    // (3*Bx*Ay-3*Ax*By)*t^2+(3*Cx*Ay-3*Cy*Ax)*t+Cx*By-Cy*Bx
    //
    // So the inflection points are at the solutions (roots) of
    //
    //   0 = (3*Bx*Ay-3*Ax*By)*t^2+(3*Cx*Ay-3*Cy*Ax)*t+Cx*By-Cy*Bx
    //     = a*t^2+b*t+c
    //
    // where a, b, & c are:
    //   a = 3*(Bx*Ay-Ax*By)
    //   b = 3*(Cx*Ay-Cy*Ax)
    //   c = Cx*By-Cy*Bx
    double a = 3*(B.x*A.y-A.x*B.y),
           b = 3*(C.x*A.y-C.y*A.x),
           c = C.x*B.y-C.y*B.x;

    double roots[2];
    int solutions;
    // Is the quadratic really a degenerate line?
#if 0
    if (a == 0) {
#else
    if (fabs(a*1000000) < fabs(b)) {
#endif
        // Is the line really a degenerate point?
        if (b == 0) {
            solutions = 0;
        } else {
            solutions = 1;
            roots[0] = -c/b;
            roots[1] = sqrt(-1.0);
        }
    } else {
        solutions = quadratic(b/a, c/a, roots);
    }
    for (int i=0; i<solutions; i++) {
        double t = roots[i];

        key_point[i] = t;
    }
    return solutions;
}

#endif // __cubicInflectionPoint_hpp__
