Function trackParabolaCore

plugin/seq/plotPDF.cpp:2081–2136

    const double a = ( m11 * rhs0 - m01 * rhs1) / det;
    const double b = (-m10 * rhs0 + m00 * rhs1) / det;

    return (0 < a) && (a < 1) && (0 < b) && (b < 1) && (0 < a+b) && (a+b < 1);
}

void trackParabolaCore( std::vector< std::vector<double> > &Cx, std::vector< std::vector<double> > &Cy,
			const double a, const double b, std::vector<double> &x,
			const double *const Vx, const double *const Vy )
{
    // y = a * x^2 + b
    std::sort( x.begin(), x.end() );

    // control points: Px[4], Py[4]
    // Bezier curve : P(t) = (1-t)^3 P0 + 3(1-t)^2 t P1 + 3(1-t)t^2 P2 + t^3 P3
    //                     = (-P0+3P1-3P2+P3) t^3 + (3P0-6P1+3P2) t^2 + (-3P0+3P1) t + P0
    //
    // parabola between x0 <= x <= x1 : 
    // P(t) = [ (x1-x0)t + x0 ; a*( (x1-x0)t + x0 )^2 + b*( (x1-x0)t + x0 ) + c ]
    //      = [ (x1-x0)t + x0 ; a(x1-x0)^2 t^2 + (2*a*x0+b)*(x1-x0) t + (a x_0^2 + bx_0 + c) ]
    //
    // Therefore
    //   P0 = [ x0 ; ax_0^2 + bx_0 + c ]
    //   -3P0+3P1 = [ x1-x0 ; (2*a*x0+b)*(x1-x0) ]
    //   3P0-6P1+3P2 = [ 0 ; a(x1-x0)^2 ]
    //   -P0+3P1-3P2+P3 = [ 0 ; 0 ]
    // <=>
    //   P0 = [ x0 ; ax_0^2 + bx_0 + c ]
    //   P1 = P0 + [ x1-x0 ; (2*a*x0+b)*(x1-x0) ]/3
    //   P2 = -P0 + 2*P1 + [ 0 ; a(x1-x0)^2 ]/3
    //   P3 = P0-3P1+3P2
    // From this, P0x = x0, P1x = x0+(x1-x0)/3, P2x = P0x + 2/3*(x1-x0), P3x = x1
	
    for(int i = 0; i+1 < x.size(); i++){

	const double &X0 = x[i];
	const double &X1 = x[i+1];
	const double h   = X1 - X0;

	const double EPS = 1e-10;
	if( h < EPS)
	    continue;

	// examine the arc X[i]--X[i+1] is inside triangle or not
	const double mX0 = X0 + h/100;
	const double mY0 = a*mX0*mX0 + b;

	const double mX1 = X1 - h/100;
	const double mY1 = a*mX1*mX1 + b;
	if( !(isInsideTriangle( mX0, mY0, Vx, Vy ) && isInsideTriangle( mX1, mY1, Vx, Vy )) )
	    continue;

	const double P0y = a*X0*X0 + b;
	const double P1y =  P0y + (2*a*X0)*h/3;
	const double P2y = -P0y + 2*P1y + a*h*h/3;
	const double P3y =  P0y - 3*P1y + 3*P2y;

	Cx.push_back( std::vector<double> { X0,  X0+h/3, X1-h/3, X1  } );
	Cy.push_back( std::vector<double> { P0y, P1y,    P2y,    P3y } );