Method FitWithConic

src/LNLib/Geometry/Curve/NurbsCurve.cpp:2806–2882

}

bool LNLib::NurbsCurve::FitWithConic(const std::vector<XYZ>& throughPoints, int startPointIndex, int endPointIndex, const XYZ& startTangent, const XYZ& endTangent, double maxError, std::vector<XYZW>& middleControlPoints)
{
	VALIDATE_ARGUMENT(throughPoints.size() > 0, "throughPoints", "ThroughPoints size must be greater than zero.");
	VALIDATE_ARGUMENT_RANGE(startPointIndex, 0, throughPoints.size() - 1);
	VALIDATE_ARGUMENT_RANGE(endPointIndex, startPointIndex + 1, throughPoints.size() - 1);
	VALIDATE_ARGUMENT(!startTangent.IsZero(), "startTangent", "StartTangent must not be zero vector.");
	VALIDATE_ARGUMENT(!endTangent.IsZero(), "endTangent", "EndTangent must not be zero vector.");

	XYZ startPoint = throughPoints[startPointIndex];
	XYZ endPoint = throughPoints[endPointIndex];
	if (endPointIndex - startPointIndex == 1)
	{
		return BezierCurve::ComputerMiddleControlPointsOnQuadraticCurve(startPoint, startTangent, endPoint, endTangent, middleControlPoints);
	}

	double alf1, alf2 = 0.0;
	XYZ R(0, 0, 0);
	CurveCurveIntersectionType type = Intersection::ComputeRays(startPoint, startTangent, endPoint, endTangent, alf1, alf2, R);
	if (type == CurveCurveIntersectionType::Coincident)
	{
		middleControlPoints.emplace_back(XYZW((startPoint + endPoint) / 2, 1));
		return true;
	}
	else if (type == CurveCurveIntersectionType::Skew || type == CurveCurveIntersectionType::Parallel)
	{
		return false;
	}
	if (MathUtils::IsLessThanOrEqual(alf1, 0.0) ||
		MathUtils::IsGreaterThanOrEqual(alf2, 0.0))
	{
		return false;
	}
	double s = 0.0;
	XYZ V = endPoint - startPoint;
	XYZ dummy;
	for (int i = startPointIndex + 1; i <= endPointIndex - 1; i++)
	{
		XYZ V1 = throughPoints[i] - R;
		type = Intersection::ComputeRays(startPoint, V, R, V1, alf1, alf2, dummy);
		if (type != CurveCurveIntersectionType::Intersecting || 
			MathUtils::IsLessThanOrEqual(alf1, 0.0) ||
			MathUtils::IsGreaterThanOrEqual(alf1, 1.0) ||
			MathUtils::IsLessThanOrEqual(alf2, 0.0))
		{
			return false;	
		}
		double wi = 0.0;
		if (CreateOneConicArc(startPoint, V, R, V1, throughPoints[i], dummy, wi))
		{
			s = s + wi / (1 + wi);
		}
	}
	
	s = s / (endPointIndex - startPointIndex - 1);
	double w = s / (1.0 - s);

	std::vector<XYZW> controlPoints = { XYZW(startPoint,1), XYZW(R,w), XYZW(endPoint,1) };
	std::vector<double> knotVectors = { 0,0,0,1.0,1.0,1.0 };