#include "tcubicbezier.h"

//-----------------------------------------------------------------------------

DVAPI double invCubicBezierX(double x,

							 const TPointD &a,

							 const TPointD &aSpeed,

							 const TPointD &bSpeed,

							 const TPointD &b)

{

	double aSpeedX = aSpeed.x;

	double bSpeedX = bSpeed.x;

	if (aSpeedX == 0)

		aSpeedX = epsilon;

	if (bSpeedX == 0)

		bSpeedX = -epsilon;

	// a*u^3+b*u^2+c*u+d

	double x0 = a.x;

	double x1 = x0 + aSpeedX;

	double x3 = b.x;

	double x2 = x3 + bSpeedX;

	double aX = x3 + 3 * (x1 - x2) - x0;

	double bX = 3 * (x2 - 2 * x1 + x0);

	double cX = 3 * (x1 - x0);

	double dX = x0 - x;

	return cubicRoot(aX, bX, cX, dX);

}

//-----------------------------------------------------------------------------

DVAPI double getCubicBezierY(double x,

							 const TPointD &a,

							 const TPointD &aSpeed,

							 const TPointD &bSpeed,

							 const TPointD &b)

{

	double y0 = a.y;

	double y1 = y0 + aSpeed.y;

	double y3 = b.y;

	double y2 = y3 + bSpeed.y;

	double aY = y3 + 3 * (y1 - y2) - y0;

	double bY = 3 * (y2 - 2 * y1 + y0);

	double cY = 3 * (y1 - y0);

	double dY = y0;

	double u = invCubicBezierX(x, a, aSpeed, bSpeed, b);

	u = tcrop(u, 0.0, 1.0);

	return aY * u * u * u + bY * u * u + cY * u + dY;

}

//-----------------------------------------------------------------------------

DVAPI std::pair<tpointd, tpointd=""> getMinMaxCubicBezierY(const TPointD &a,</tpointd,>

														const TPointD &aSpeed,

														const TPointD &bSpeed,

														const TPointD &b)

{

	double y0 = a.y;

	double y1 = y0 + aSpeed.y;

	double y3 = b.y;

	double y2 = y3 + bSpeed.y;

	double aY = y3 + 3 * (y1 - y2) - y0;

	double bY = 3 * (y2 - 2 * y1 + y0);

	double cY = 3 * (y1 - y0);

	double dY = y0;

	double x0 = a.x;

	double x1 = x0 + aSpeed.x;

	double x3 = b.x;

	double x2 = x3 + bSpeed.x;

	double aX = x3 + 3 * (x1 - x2) - x0;

	double bX = 3 * (x2 - 2 * x1 + x0);

	double cX = 3 * (x1 - x0);

	double dX = x0;

	double aaY = 3 * (y1 - y2) - y0 + y3;

	double bbY = 2 * (y0 + y2 - 2 * y1);

	double ccY = y1 - y0;

	if (aaY == 0) {

		if (a.y < b.y)

			return std::pair<tpointd, tpointd="">(TPointD(a.x, a.y),</tpointd,>

											   TPointD(b.x, b.y));

		else

			return pair<tpointd, tpointd="">(TPointD(b.x, b.y),</tpointd,>

										  TPointD(a.x, a.y));

	} else {

		double discr = bbY * bbY - 4 * aaY * ccY;

		if (discr < 0) {

			if (a.y < b.y)

				return std::pair<tpointd, tpointd="">(TPointD(a.x, a.y),</tpointd,>

												   TPointD(b.x, b.y));

			else

				return pair<tpointd, tpointd="">(TPointD(b.x, b.y),</tpointd,>

											  TPointD(a.x, a.y));

		} else {

			double sqrt_discr = sqrt(discr);

			double inv_2aY = 1.0 / (2.0 * aaY);

			double u0 = (-bbY + sqrt_discr) * inv_2aY;

			double u1 = (-bbY - sqrt_discr) * inv_2aY;

			if (u0 > 1.0)

				u0 = 1.0;

			if (u0 < 0.0)

				u0 = 0.0;

			if (u1 > 1.0)

				u1 = 1.0;

			if (u1 < 0.0)

				u1 = 0.0;

			double y_0 = aY * u0 * u0 * u0 + bY * u0 * u0 + cY * u0 + dY;

			double y_1 = aY * u1 * u1 * u1 + bY * u1 * u1 + cY * u1 + dY;

			double x_0 = aX * u0 * u0 * u0 + bX * u0 * u0 + cX * u0 + dX;

			double x_1 = aX * u1 * u1 * u1 + bX * u1 * u1 + cX * u1 + dX;

			if (y_0 < y_1)

				return pair<tpointd, tpointd="">(TPointD(x_0, y_0), TPointD(x_1, y_1));</tpointd,>

			else

				return pair<tpointd, tpointd="">(TPointD(x_1, y_1), TPointD(x_0, y_0));</tpointd,>

		}

	}

}

//-----------------------------------------------------------------------------