Catmull-Rom bounding box 2D

gistfile1.txt

	//
	// CatmullRomBoundinBox2D()
	//
	//
	//*********************************************************************************************
	template <typename T_Vec2>
	void CatmullRomBoundingBox2D(const T_Vec2 ctrlPts[4], T_Vec2 &min, T_Vec2 &max)
	{
		// See this for description of concept:
		// https://stackoverflow.com/questions/24809978/calculating-the-bounding-box-of-cubic-bezier-curve
		//
		// Note however, that there is difference between how Bezier / Catmull-Rom polynomials are constructed
		// from control points !
			
		typedef typename T_Vec2::T_ValueType T;

		const T_Vec2 polyA = 1.5f * (ctrlPts[1] - ctrlPts[2]) + 0.5f * (ctrlPts[3] - ctrlPts[0]);
		const T_Vec2 polyB = ctrlPts[0] - 2.5f * ctrlPts[1] + 2.0f * ctrlPts[2] - ctrlPts[3] * 0.5f;
		const T_Vec2 polyC = (ctrlPts[2] - ctrlPts[0]) * 0.5f;
		const T_Vec2 polyD = ctrlPts[1];

		const T_Vec2 dpolyA = 4.5f * (ctrlPts[1] - ctrlPts[2]) + 1.5f * (ctrlPts[3] - ctrlPts[0]);
		const T_Vec2 dpolyB = 2.0f * ctrlPts[0] - 5.0f * ctrlPts[1] + 4.0f * ctrlPts[2] - ctrlPts[3];
		const T_Vec2 dpolyC = (ctrlPts[2] - ctrlPts[0]) * 0.5f;

		const T_Vec2 D = dpolyB * dpolyB - 4.0f * dpolyA * dpolyC;

		min[0] = FMin(ctrlPts[1][0], ctrlPts[2][0]);
		min[1] = FMin(ctrlPts[1][1], ctrlPts[2][1]);

		max[0] = FMax(ctrlPts[1][0], ctrlPts[2][0]);
		max[1] = FMax(ctrlPts[1][1], ctrlPts[2][1]);

		if (D[0] >= T(0.000001))
		{
			const T	Ds			= Sqrt(D[0]);
			const T oneOver2A	= T(1) /  (dpolyA[0] * 2);
			const T	root0		= (-dpolyB[0] + Ds) * oneOver2A;
			const T	root1		= (-dpolyB[0] - Ds) * oneOver2A;

			if (root0 > T(0) && root0 < T(1))
			{
				const T tmp = ((polyA[0] * root0 + polyB[0]) * root0 + polyC[0]) * root0 + polyD[0];

				min[0] = FMin(min[0], tmp);
				max[0] = FMax(max[0], tmp);
			}
			
			if (root1 > T(0) && root1 < T(1))
			{
				const T tmp = ((polyA[0] * root1 + polyB[0]) * root1 + polyC[0]) * root1 + polyD[0];

				min[0] = FMin(min[0], tmp);
				max[0] = FMax(max[0], tmp);
			}
		}

		if (D[1] >= T(0.000001))
		{
			const T	Ds			= Sqrt(D[1]);
			const T oneOver2A	= T(1) /  (dpolyA[1] * 2);
			const T	root0		= (-dpolyB[1] + Ds) * oneOver2A;
			const T	root1		= (-dpolyB[1] - Ds) * oneOver2A;

			if (root0 > T(0) && root0 < T(1))
			{
				const T tmp = ((polyA[1] * root0 + polyB[1]) * root0 + polyC[1]) * root0 + polyD[1];

				min[1] = FMin(min[1], tmp);
				max[1] = FMax(max[1], tmp);
			}

			if (root1 > T(0) && root1 < T(1))
			{
				const T tmp = ((polyA[1] * root1 + polyB[1]) * root1 + polyC[1]) * root1 + polyD[1];

				min[1] = FMin(min[1], tmp);
				max[1] = FMax(max[1], tmp);
			}
		}
	}
	//*********************************************************************************************