#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,
														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(b.x, b.y));
		else
			return pair<TPointD, TPointD>(TPointD(b.x, b.y),
										  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(b.x, b.y));
			else
				return pair<TPointD, TPointD>(TPointD(b.x, b.y),
											  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));
			else
				return pair<TPointD, TPointD>(TPointD(x_1, y_1), TPointD(x_0, y_0));
		}
	}
}

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