diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl new file mode 100644 index 0000000000..06e8471af4 --- /dev/null +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -0,0 +1,640 @@ +#ifndef _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ +#define _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ + +// TODO: make this file HLSL compatible +#ifndef __HLSL_VERSION + +#include +#include +#include + +#include +#include +#include +#include + +namespace nbl +{ +namespace hlsl +{ +// TODO: do we want to keep it in the nbl::hlsl::ext namespace? this code is not generalized and used for specific scenario so maybe it should stay in the `ext` namespace? +namespace ext +{ +namespace curves +{ + +// Base class for all our curves +template +struct ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + //! compute position at t + virtual float_t2 computePosition(float_t t) const = 0; + + //! compute unnormalized tangent vector at t + virtual float_t2 computeTangent(float_t t) const = 0; + + //! compute differential arc length at t + virtual float_t differentialArcLen(float_t t) const + { + return nbl::hlsl::length(computeTangent(t)); + } + + struct ArcLenIntegrand + { + const ParametricCurve* m_curve; + + ArcLenIntegrand(const ParametricCurve* curve) + : m_curve(curve) + {} + + inline float_t operator()(const float_t t) const + { + return m_curve->differentialArcLen(t); + } + }; + + //! compute arc length by gauss legendere integration + float_t arcLen(float_t t0, float_t t1) const + { + constexpr uint16_t IntegrationOrder = 10u; + return nbl::hlsl::math::quadrature::GaussLegendreIntegration::calculateIntegral(ArcLenIntegrand(this), t0, t1); + } + + //! compute inverse arc len using bisection search + float_t inverseArcLen_BisectionSearch(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4, const uint16_t iterationThreshold = 16u) const + { + float_t xi = 0.0; + float_t low = min; + float_t high = max; + for (uint16_t i = 0; i < iterationThreshold; ++i) + { + xi = (low + high) / 2.0; + float_t sum = arcLen(min, xi); + float_t integral = sum + arcLen(xi, max); + + // we could've done sum/integral - targetLen, but this is more robust as it avoids a divsion + float_t valueAtParamGuess = sum - targetLen * integral; + + if (abs(valueAtParamGuess) < cdfAccuracyThreshold * integral) + return xi; // we found xi value that gives us a cdf of targetLen within cdfAccuracyThreshold + else + { + if (valueAtParamGuess > 0.0) + high = xi; + else + low = xi; + } + } + + return xi; + } + + //! compute inverse arc len + float_t inverseArcLen(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4) const + { + return inverseArcLen_BisectionSearch(targetLen, min, max, cdfAccuracyThreshold); + } + + virtual float_t2 computeSecondOrderDifferential(float_t t) const + { + return float_t2(numeric_limits::quiet_NaN, numeric_limits::quiet_NaN); + } + + virtual float_t computeInflectionPoint(float_t errorThreshold) const + { + return numeric_limits::quiet_NaN; + } +}; + +// TODO: make an `ExplicitCurve` concept, it should require for a type to define `y` and `derivative` + +// It's when t = x in a Parametric Curve +template +struct ExplicitCurve : public ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + virtual float_t y(float_t x) const = 0; + virtual float_t derivative(float_t x) const = 0; + + float_t differentialArcLen(float_t x) const override + { + float_t deriv = derivative(x); + return sqrt(1.0 + deriv * deriv); + } + + float_t2 computeTangent(float_t x) const override + { + const float_t deriv = derivative(x); + float_t2 v = float_t2(1.0, deriv); + if (std::isinf(deriv)) + v = float_t2(0.0, 1.0); + return v; + } + + inline float_t2 computePosition(float_t x) const override { return float_t2(x, y(x)); } +}; + +template +struct Parabola : ExplicitCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b, c; + + static Parabola create(float_t a, float_t b, float_t c) + { + Parabola output; + output.a = a; + output.b = b; + output.c = c; + + return output; + } + + static Parabola fromThreePoints(NBL_CONST_REF_ARG(float_t2) P0, NBL_CONST_REF_ARG(float_t2) P1, NBL_CONST_REF_ARG(float_t2) P2) + { + glm::dmat3 X = glm::dmat3( + glm::dvec3(P0.x * P0.x, P0.x, 1.0), + glm::dvec3(P1.x * P1.x, P1.x, 1.0), + glm::dvec3(P2.x * P2.x, P2.x, 1.0) + ); + glm::dvec3 M = inverse(transpose(X)) * glm::dvec3(P0.y, P1.y, P2.y); + return Parabola(M[0], M[1], M[2]); + } + + float_t y(float_t x) const + { + return ((a * x) + b) * x + c; + } + + float_t derivative(float_t x) const + { + return 2.0 * a * x + b; + } +}; + +template +struct CubicCurve : ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + using float_t4 = nbl::hlsl::portable_vector_t4; + + float_t4 X; + float_t4 Y; + + static CubicCurve create(NBL_CONST_REF_ARG(float_t4) X, NBL_CONST_REF_ARG(float_t4) Y) + { + CubicCurve output; + output.X = X; + output.Y = Y; + + return output; + } + + float_t2 computePosition(float_t t) const override + { + return float_t2( + ((X[0] * t + X[1]) * t + X[2]) * t + X[3], + ((Y[0] * t + Y[1]) * t + Y[2]) * t + Y[3] + ); + } + + //! compute unnormalized tangent vector at t + float_t2 computeTangent(float_t t) const override + { + return float_t2( + (3.0 * X[0] * t + 2.0 * X[1]) * t + X[2], + (3.0 * Y[0] * t + 2.0 * Y[1]) * t + Y[2] + ); + } + + //! compute second order differential at t + float_t2 computeSecondOrderDifferential(float_t t) const override + { + return float_t2( + 6.0 * X[0] * t + 2.0 * X[1], + 6.0 * Y[0] * t + 2.0 * Y[1] + ); + } + + float_t differentialArcLen(float_t t) const override + { + return nbl::hlsl::length(computeTangent(t)); + } + + float_t computeInflectionPoint(float_t errorThreshold) const override + { + // solve for signed curvature root + // when x'*y''-x''*y' = 0 + // https://www.wolframalpha.com/input?i=cross+product+%283*x0*t%5E2%2B2*x1%2Bx2%2C3*y0*t%5E2%2B2*y1%2By2%29+and+%286*x0*t%2B2*x1%2C6*y0*t%2B2*y1%29 + const float_t a = 6.0 * (X[0] * Y[1] - X[1] * Y[0]); + const float_t b = 6.0 * (2.0 * X[1] * Y[0] - 2.0 * X[0] * Y[1] + X[2] * Y[0] - X[0] * Y[2]); + const float_t c = 2.0 * (X[2] * Y[1] - X[1] * Y[2]); + + nbl::hlsl::math::equations::Quadratic quadratic = nbl::hlsl::math::equations::Quadratic::construct(a, b, c); + const float_t2 roots = quadratic.computeRoots(); + if (roots[0] <= 1.0 && roots[0] >= 0.0) + return roots[0]; + if (roots[1] <= 1.0 && roots[1] >= 0.0) + return roots[1]; + return numeric_limits::quiet_NaN; + } +}; + +// specialized circular arc for the purpose of mixing it with another curve of the same type later +// (r*cos(t*sweep+start), r*sin(t*sweep+start) + originY) +template +struct CircularArc : ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t r; + float_t originY; // originX is 0 + float_t startAngle; + float_t sweepAngle; + + static CircularArc create(float_t r, float_t originY, float_t startAngle, float_t sweepAngle) + { + CircularArc output; + output.r = r; + output.originY = originY; + output.startAngle = startAngle; + output.sweepAngle = sweepAngle; + + return output; + } + + // from circle center (0, -v.y) to start pos (v.x, 0) + static CircularArc create(float_t2 v, float_t sweepAngle) + { + CircularArc output; + output.originY = -v.y; + output.sweepAngle = sweepAngle; + + output.r = length(v); + output.startAngle = getSign(v.y) * acos(v.x / output.r); + + return output; + } + + // from circle center (0, -v.y) to start pos (v.x, 0) + static CircularArc create(float_t2 v) + { + CircularArc output; + output.originY = -v.y; + + output.r = length(v); + output.startAngle = getSign(v.y) * acos(v.x / output.r); + output.sweepAngle = -2.0 * getSign(v.y) * acos(abs(output.originY) / output.r); + + return output; + } + + float_t2 computePosition(float_t t) const override + { + const float_t actualT = t * sweepAngle + startAngle; + return float_t2( + r * cos(actualT), + r * sin(actualT) + originY + ); + } + + //! compute unnormalized tangent vector at t + float_t2 computeTangent(float_t t) const override + { + const float_t actualT = t * sweepAngle + startAngle; + return float_t2( + -1.0 * r * sweepAngle * sin(actualT), + +1.0 * r * sweepAngle * cos(actualT) + ); + } + + float_t2 computeSecondOrderDifferential(float_t t) const override + { + const float_t actualT = t * sweepAngle + startAngle; + return float_t2( + -1.0 * r * sweepAngle * sweepAngle * cos(actualT), + -1.0 * r * sweepAngle * sweepAngle * sin(actualT) + ); + } + + static float_t getSign(float_t x) + { + return static_cast((x > 0.0)) - static_cast((x <= 0.0)); + } +}; + +// Centered at (0,0), aligned with x axis +template +struct ExplicitEllipse : ExplicitCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b; + static ExplicitEllipse create(float_t a, float_t b) + { + ExplicitEllipse output; + output.a = a; + output.b = b; + + return output; + } + + float_t y(float_t x) const + { + return a * sqrt(static_cast(1.0) - pow((x / b), static_cast(2.0))); + } + + float_t derivative(float_t x) const + { + return (-a * x) / ((b * b) * sqrt(static_cast(1.0) - pow((x / b), static_cast(2.0)))); + } +}; + +// Centered at (0,0), aligned with x +template +struct AxisAlignedEllipse : ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b; + float_t start, end; + static AxisAlignedEllipse create(float_t a, float_t b, float_t start, float_t end) + { + AxisAlignedEllipse output; + output.a = a; + output.b = b; + output.start = start; + output.end = end; + + return output; + } + + float_t2 computePosition(float_t t) const override + { + const float_t theta = start + (end - start) * t; + return float_t2(a * cos(theta), b * sin(theta)); + } + + float_t2 computeTangent(float_t t) const override + { + const float_t theta = start + (end - start) * t; + const float_t dThetaDt = end - start; + return float_t2(-a * dThetaDt * sin(theta), b * dThetaDt * cos(theta)); + } +}; + +template +struct EllipticalArcInfo +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t2 majorAxis; + float_t2 center; + float_t2 angleBounds; // [0, 2Pi) + float_t eccentricity; // (0, 1] + + inline bool isValid() const + { + if (eccentricity > 1.0 || eccentricity <= 0.0) + return false; + if (angleBounds.y == angleBounds.x) + return false; + if (abs(angleBounds.y - angleBounds.x) > 2.0 * nbl::core::PI()) + return false; + return true; + } +}; + +template +struct OffsettedBezier : ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + nbl::hlsl::shapes::Quadratic quadratic; + float_t offset; + + static OffsettedBezier create(NBL_CONST_REF_ARG(nbl::hlsl::shapes::QuadraticBezier) quadBezier, float_t offset) + { + OffsettedBezier output; + output.offset = offset; + + output.quadratic = nbl::hlsl::shapes::Quadratic::constructFromBezier(quadBezier.P0, quadBezier.P1, quadBezier.P2); + + return output; + } + + float_t2 computePosition(float_t t) const override + { + const float_t2 deriv = quadratic.derivative(t); + const float_t2 normal = normalize(float_t2(deriv.y, -deriv.x)); + return quadratic.evaluate(t) + offset * normal; + } + + //! compute unnormalized tangent vector at t + float_t2 computeTangent(float_t t) const override + { + const float_t2 ddt = quadratic.derivative(t); + const float_t2 d2dt2 = quadratic.secondDerivative(t); + const float_t g = offset * (ddt.x * d2dt2.y - ddt.y * d2dt2.x); + return ddt + (ddt * g) / glm::length(ddt); + } + + //! if offset is more than minimum radius of curvature then we get an unwanted gouging/cusp + float_t2 findCusps() const + { + // we're basically solving for t in "offset = radiusOfCurvature(t)" + const float_t lhs = pow(offset * 2.0 * abs(quadratic.B.x * quadratic.A.y - quadratic.B.y * quadratic.A.x), 2.0 / 3.0); + const float_t a = 4.0 * (quadratic.A.x * quadratic.A.x + quadratic.A.y * quadratic.A.y); + const float_t b = 4.0 * (quadratic.A.x * quadratic.B.x + quadratic.A.y * quadratic.B.y); + const float_t c = quadratic.B.x * quadratic.B.x + quadratic.B.y * quadratic.B.y - lhs; + nbl::hlsl::math::equations::Quadratic findCuspsQuadratic = nbl::hlsl::math::equations::Quadratic::construct(a, b, c); + return findCuspsQuadratic.computeRoots(); + } +}; + +template +class QuadraticBezierFitter final +{ +public: + using float_t2 = nbl::hlsl::portable_vector_t2; + using float_t3 = nbl::hlsl::portable_vector_t3; + using float_t4 = nbl::hlsl::portable_vector_t4; + using float_t2x2 = nbl::hlsl::portable_matrix_t2x2; + + typedef std::function&&)> AddBezierFunc; + + //! this subdivision algorithm works/converges for any x-monotonic curve (only 1 y for each x) over the [min, max] range and will continue until hits the `maxDepth` or `targetMaxError` threshold + //! this function will call the AddBezierFunc when the bezier is finalized, whether to render it directly, write it to file, add it to a vector, etc.. is up to the user. + //! the subdivision samples the points based on arc length and the error is computed by distance in y direction, so pre and post transform may be needed for your curve and the outputted beziers + //! it will first split at inflection point of the curve; curves are assumed to have at most 1 inflection point, and will get the best convergence rates. but it will work for curves with more inflection points as well. + static void adaptive(const ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + { + // The curves we're working with will have at most 1 inflection point. + const float_t inflectX = curve.computeInflectionPoint(targetMaxError); // if no inflection point then this will return NaN and the adaptive subdivision will continue as normal (from min to max) + if (inflectX > min && inflectX < max) + { + adaptive_impl(curve, min, inflectX, targetMaxError, addBezierFunc, maxDepth); + adaptive_impl(curve, inflectX, max, targetMaxError, addBezierFunc, maxDepth); + } + else + adaptive_impl(curve, min, max, targetMaxError, addBezierFunc, maxDepth); + } + + static void adaptive(const EllipticalArcInfo& ellipse, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + { + using namespace nbl::hlsl; + + if (!ellipse.isValid()) + { + _NBL_DEBUG_BREAK_IF(true); + return; + } + + float_t lenghtMajor = length(ellipse.majorAxis); + float_t lenghtMinor = lenghtMajor * ellipse.eccentricity; + float_t2 normalizedMajor = ellipse.majorAxis / lenghtMajor; + + float_t2x2 rotate = float_t2x2({ + float_t2(normalizedMajor.x, -normalizedMajor.y), + float_t2(normalizedMajor.y, normalizedMajor.x) + }); + + AddBezierFunc addTransformedBezier = [&](nbl::hlsl::shapes::QuadraticBezier&& quadBezier) + { + quadBezier.P0 = mul(rotate, quadBezier.P0); + quadBezier.P1 = mul(rotate, quadBezier.P1); + quadBezier.P2 = mul(rotate, quadBezier.P2); + quadBezier.P0 += ellipse.center; + quadBezier.P1 += ellipse.center; + quadBezier.P2 += ellipse.center; + addBezierFunc(std::move(quadBezier)); + }; + + if (ellipse.angleBounds.x != ellipse.angleBounds.y) + { + AxisAlignedEllipse aaEllipse = AxisAlignedEllipse::create(lenghtMajor, lenghtMinor, ellipse.angleBounds.x, ellipse.angleBounds.y); + adaptive(aaEllipse, 0.0, 1.0, targetMaxError, addTransformedBezier, maxDepth); + } + } + + static void adaptive(const OffsettedBezier& curve, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + { + const float_t2 cusps = curve.findCusps(); + + const float_t t0 = nbl::core::min(cusps[0], cusps[1]); + const float_t t1 = nbl::core::max(cusps[0], cusps[1]); + + const bool firstCusp = t0 > 0.0 && t0 < 1.0; + const bool secondCusp = t1 > 0.0 && t1 < 1.0; + + // if there are two cusps (offset = radius of curvature) then we have that unwanted gouging and we prefer to seperately subdivide those three sections + if (firstCusp && secondCusp) + { + adaptive_impl(curve, 0.0, t0, targetMaxError, addBezierFunc, maxDepth); + adaptive_impl(curve, t0, t1, targetMaxError, addBezierFunc, maxDepth); + adaptive_impl(curve, t1, 1.0, targetMaxError, addBezierFunc, maxDepth); + } + // otherwise just subdivide from start/0.0 to end/1.0 + else + { + adaptive_impl(curve, 0.0, 1.0, targetMaxError, addBezierFunc, maxDepth); + } + } + +private: + // Fix Bezier Hack for when P1 is "outside" P0 -> P2 + // We project P1 into P0->P2 line and see whether it lies inside. + // Because our curves shouldn't go back on themselves in the direction of the chord + static void fixBezierMidPoint(nbl::hlsl::shapes::QuadraticBezier& bezier) + { + const float_t2 localChord = bezier.P2 - bezier.P0; + const float_t localX = dot(normalize(localChord), bezier.P1 - bezier.P0); + const bool outside = localX<0 || localX>length(localChord); + if (outside || nbl::core::isnan(bezier.P1.x) || nbl::core::isnan(bezier.P1.y)) + { + // _NBL_DEBUG_BREAK_IF(true); // this shouldn't happen but we fix it just in case anyways + bezier.P1 = bezier.P0 * 0.4 + bezier.P2 * 0.6; + } + } + + static void adaptive_impl(const ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t depth) + { + if (min == max) + return; + assert(min < max); + + float_t split = curve.inverseArcLen_BisectionSearch(0.5, min, max); + + // Shouldn't happen but may happen if we use NewtonRaphson for non convergent inverse CDF + if (split <= min || split >= max) + { + _NBL_DEBUG_BREAK_IF(split < min || split > max); + split = (min + max) / 2.0; + } + + const float_t2 P0 = curve.computePosition(min); + const float_t2 V0 = curve.computeTangent(min); + const float_t2 P2 = curve.computePosition(max); + const float_t2 V2 = curve.computeTangent(max); + nbl::hlsl::shapes::QuadraticBezier bezier = nbl::hlsl::shapes::QuadraticBezier::constructBezierWithTwoPointsAndTangents(P0, V0, P2, V2); + + bool shouldSubdivide = false; + + // TODO: compare with certain threshold + if (depth > 0u && normalize(V0) == normalize(V2)) + { + shouldSubdivide = true; + } + else + { + fixBezierMidPoint(bezier); + if (depth > 0u) + { + if (glm::distance(P0, P2) < targetMaxError) + { + const float_t2 posAtSplit = curve.computePosition(split); + // If it came down to a bezier small that causes P0 P2 and the position at split smaller than targetMaxError then we stop + if (glm::distance(posAtSplit, P0) < targetMaxError) + shouldSubdivide = false; + // But sometimes when P0 and P2 are close together a split will fix them, like a full circle and needs further subdivision + else + shouldSubdivide = true; + } + else + { + const float_t2 curvePositionAtSplit = curve.computePosition(split); + float_t bezierYAtSplit = bezier.calcYatX(curvePositionAtSplit.x); + //_NBL_DEBUG_BREAK_IF(nbl::core::isnan(bezierYAtSplit)); + // TODO: maybe a better error comaprison is find the normal at split and intersect with the bezier + if (nbl::core::isnan(bezierYAtSplit) || abs(curvePositionAtSplit.y - bezierYAtSplit) > targetMaxError) + shouldSubdivide = true; + } + } + } + + if (shouldSubdivide) + { + adaptive_impl(curve, min, split, targetMaxError, addBezierFunc, depth - 1u); + adaptive_impl(curve, split, max, targetMaxError, addBezierFunc, depth - 1u); + } + else + { + const bool degenerate = (bezier.P0 == bezier.P2); + if (!degenerate) + addBezierFunc(std::move(bezier)); + } + } +}; + +} // namespace curves +} // namespace ext +} // namespace hlsl +} // namespace nbl + +#endif + +#endif \ No newline at end of file diff --git a/include/nbl/ext/Hatch/Hatch.h b/include/nbl/ext/Hatch/Hatch.h new file mode 100644 index 0000000000..0e15d0e8da --- /dev/null +++ b/include/nbl/ext/Hatch/Hatch.h @@ -0,0 +1,884 @@ +#ifndef _NBL_EXT_HATCH_H_ +#define _NBL_EXT_HATCH_H_ + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace nbl::ext::csg2d +{ + +using namespace nbl::core; +using nbl::hlsl::shapes::QuadraticBezier; +using nbl::hlsl::shapes::Quadratic; + +enum class SweepMajorAxis +{ + MAJOR_X = 0u, + MAJOR_Y = 1u, +}; + +/** + * @brief A geometric processing engine that decomposes closed Bezier polygons into fillable slabs. + * + * The HatchBuilder utilizes a algorithm (similar to sweep-line but for beziers) to process a collection of quadratic Bezier + * curves that form closed polygons. As the sweep-line moves along the defined MajorAxis, it + * splits and sorts the curves to resolve complex overlaps and self-intersections. + * * The algorithm outputs these resolved areas as `CurveBox` structures (often called slabs or + * monotone regions). Each CurveBox represents a discrete, continuous region bounded by two + * minor-axis curves (the left and right boundaries if SweepMajorAxis is MAJOR_Y). + * * Crucially, this decomposition is governed by the **XOR (Even-Odd) fill rule**. By identifying + * where curves begin, end, and intersect, the builder pairs the active boundaries so that the + * resulting CurveBoxes represent exactly the interior regions that should be filled to achieve a + * topologically correct XOR fill of the original polygons. + * + * @tparam float_t The floating-point precision type used for all internal geometric calculations (e.g., float, double). + * @tparam MajorAxis The primary axis direction for the sweep-line to travel across. Defaults to SweepMajorAxis::MAJOR_Y. + */ +template +class HatchBuilder +{ +public: + // Constants calculated at compile time + static constexpr int MajorIdx = static_cast(MajorAxis); + static constexpr int MinorIdx = 1 - MajorIdx; + + using float_t2 = nbl::hlsl::portable_vector_t2; + using float_t3 = nbl::hlsl::portable_vector_t3; + using float_t4 = nbl::hlsl::portable_vector_t4; + using float_t2x2 = nbl::hlsl::portable_matrix_t2x2; + + struct CurveBox + { + QuadraticBezier minCurve; // left curve if Major is Y + QuadraticBezier maxCurve; // right curve if Major is Y + float_t2 aabbMin; + float_t2 aabbMax; + }; + + struct Segment + { + const QuadraticBezier* originalBezier = nullptr; + // because beziers are broken down, depending on the type this is t_start or t_end + float_t t_start; + float_t t_end; // beziers get broken down + }; + + /** + * @brief Constructs a new HatchBuilder. + * @param logger An optional smart pointer to the Nabla system logger. + * @param majorAxis The primary axis for the sweep-line algorithm (0 for X, 1 for Y). + */ + HatchBuilder(nbl::system::logger_opt_smart_ptr logger) + : m_logger(std::move(logger)) + { + } + + + /** + * @brief Ensures a bezier is strictly monotonic along the major axis before adding it. + * Reorders control points if the curve is running backwards along the sweep. + * * @param bezier The pre-split, guaranteed monotonic bezier curve. + */ + void addMonotonicBezier(const QuadraticBezier& bezier) + { + QuadraticBezier outputBezier = bezier; + + // Ensure strictly increasing along major axis + if (outputBezier.P0[MajorIdx] > outputBezier.P2[MajorIdx]) + { + outputBezier.P2 = bezier.P0; + outputBezier.P0 = bezier.P2; + } + + // Minor precision fix + if (outputBezier.P1[MajorIdx] < outputBezier.P0[MajorIdx]) + outputBezier.P1[MajorIdx] = outputBezier.P0[MajorIdx]; + +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput) + { + uint32_t bezierIdx = beziers.size(); + float32_t4 colors[5] = { + float32_t4(33,150,243, 255) / float32_t4(255.0), + float32_t4(29,233,182, 255) / float32_t4(255.0), + float32_t4(238,255,65, 255) / float32_t4(255.0), + float32_t4(244,81,30, 255) / float32_t4(255.0), + float32_t4(211,47,47, 255) / float32_t4(255.0) + }; + //drawDebugBezier(bezier, colors[bezierIdx % 5]); + } +#endif + + m_beziers.push_back(outputBezier); + } + + + /** + * @brief Adds a single bezier curve to the builder's internal queue. + * * @param bezier The quadratic bezier curve to add. + */ + void addBezier(const QuadraticBezier& bezier) + { + std::array, 2> monotonicSegments; + + // Assuming Impl exists and provides this functionality + bool isMonotonic = splitIntoMajorMonotonicSegments(bezier, monotonicSegments); + + if (isMonotonic) + { + addMonotonicBezier(bezier); +#ifdef DEBUG_HATCH_VISUALLY + //if (debugOutput) + //drawDebugBezier(unsplitBezier, float32_t4(0.8, 0.8, 0.8, 1.0)); +#endif + } + else + { + addMonotonicBezier(monotonicSegments[0]); + addMonotonicBezier(monotonicSegments[1]); +#ifdef DEBUG_HATCH_VISUALLY + //if (debugOutput) + //{ drawDebugBezier(monotonicSegments.data()[0], float32_t4(1.0, 0.0, 0.0, 1.0)); drawDebugBezier(monotonicSegments.data()[1], float32_t4(0.0, 1.0, 0.0, 1.0)); } +#endif + } + } + + /** + * @brief Clears all queued segments. + * Useful for reusing the builder instance without reallocating memory. + */ + void clear() + { + m_beziers.clear(); + } + + // ======================================================================== + // REQUIRED DERIVED INTERFACE + // ======================================================================== + /** + * @brief Called by process() whenever a completed monotonic region (Curve Box/Slab) is generated. + * * This pure virtual function must be implemented by the derived class. It acts as the + * primary callback during the sweep-line algorithm, allowing the user to process, render, or store the resulting geometric slabs. + * @param curveBox The generated CurveBox + */ + virtual void consumeCurveBox(const CurveBox& curveBox) = 0; + + /** + * @brief Executes the CSG/Sweep-line algorithm on all added beziers. + * As CurveBoxes are discovered, this function will automatically trigger `consumeCurveBox()`. + */ + void build() + { + // this threshsold is used to decide when to consider MinorIdx position to be + // the same and check tangents because intersection algorithms has rounding errors + constexpr float_t MinorPositionComparisonThreshhold = 1e-3; + constexpr float_t TangentComparisonThreshhold = 1e-7; + + std::stack starts; // Next segments sorted by start points + std::stack ends; // Next end points + std::priority_queue, std::greater > intersections; // Next intersection points as MajorIdx coordinate + float_t maxMajor; + +#ifdef DEBUG_HATCH_VISUALLY + + /* old parameters int32_t* debugStepPtr, const std::function& debugOutput */ + + int32_t debugStepDefault = 0u; + int32_t& debugStep = (debugStepPtr) ? *debugStepPtr : debugStepDefault; + + auto drawDebugBezier = [&](QuadraticBezier bezier, float32_t4 color) + { + CPolyline outputPolyline; + std::vector> beziers; + shapes::QuadraticBezier tmpBezier; + tmpBezier.P0 = bezier.P0; + tmpBezier.P1 = bezier.P1; + tmpBezier.P2 = bezier.P2; + beziers.push_back(tmpBezier); + outputPolyline.addQuadBeziers(beziers); + + LineStyleInfo lineStyleInfo; + lineStyleInfo.screenSpaceLineWidth = 4.0f; + lineStyleInfo.worldSpaceLineWidth = 0.0f; + lineStyleInfo.color = color; + + debugOutput(outputPolyline, lineStyleInfo); + }; + + auto drawDebugLine = [&](float64_t2 start, float64_t2 end, float32_t4 color) + { + CPolyline outputPolyline; + std::vector points; + points.push_back(start); + points.push_back(end); + outputPolyline.addLinePoints(points); + + LineStyleInfo lineStyleInfo; + lineStyleInfo.screenSpaceLineWidth = 2.0f; + lineStyleInfo.worldSpaceLineWidth = 0.0f; + lineStyleInfo.color = color; + + debugOutput(outputPolyline, lineStyleInfo); + }; +#endif + + // Inititalize start segments in a sorted manner + { + std::vector segments; + + for (uint32_t bezierIdx = 0; bezierIdx < m_beziers.size(); bezierIdx++) + { + auto hatchBezier = &m_beziers[bezierIdx]; + Segment segment; + segment.originalBezier = hatchBezier; + segment.t_start = 0.0; + segment.t_end = 1.0; + segments.push_back(segment); + } + + if (segments.empty()) + { + m_logger.log("Empty Polylines with no segments were fed into the Hatch construction.", nbl::system::ILogger::ELL_WARNING); + return; + } + + std::sort(segments.begin(), segments.end(), [&](const Segment& a, const Segment& b) { return a.originalBezier->P0[MajorIdx] > b.originalBezier->P0[MajorIdx]; }); + for (Segment& segment : segments) + starts.push(segment); + + std::sort(segments.begin(), segments.end(), [&](const Segment& a, const Segment& b) { return a.originalBezier->P2[MajorIdx] > b.originalBezier->P2[MajorIdx]; }); + for (Segment& segment : segments) + ends.push(segment.originalBezier->P2[MajorIdx]); + maxMajor = segments.front().originalBezier->P2[MajorIdx]; + } + + +#ifdef DEBUG_HATCH_VISUALLY + int32_t step = 0; +#endif + + // Sweep line algorithm + std::vector activeCandidates; // Set of active candidates for neighbor search in sweep line + + // if we weren't spawning quads, we could just have unsorted `vector` + auto candidateComparator = [&](const Segment& lhs, const Segment& rhs) + { + // btw you probably want the beziers in Quadratic At^2+B+C form, not control points + double _lhs = lhs.originalBezier->evaluate(lhs.t_start)[MinorIdx]; + double _rhs = rhs.originalBezier->evaluate(rhs.t_start)[MinorIdx]; + + double lenLhs = glm::distance(lhs.originalBezier->P0, lhs.originalBezier->P2); + double lenRhs = glm::distance(rhs.originalBezier->P0, rhs.originalBezier->P2); + auto minLen = std::min(lenLhs, lenRhs); +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && step == debugStep) + { + //printf(std::format("comparison: lhs = ({}, {}), ({}, {}), ({}, {}) rhs = ({}, {}), ({}, {}), ({}, {})", + // lhs.originalBezier->P0.x, lhs.originalBezier->P0.y, + // lhs.originalBezier->P1.x, lhs.originalBezier->P1.y, + // lhs.originalBezier->P2.x, lhs.originalBezier->P2.y, + + // rhs.originalBezier->P0.x, rhs.originalBezier->P0.y, + // rhs.originalBezier->P1.x, rhs.originalBezier->P1.y, + // rhs.originalBezier->P2.x, rhs.originalBezier->P2.y + // ).c_str()); + //drawDebugLine(float64_t2(_lhs, -1000.0), float64_t2(_lhs, 1000.0), float64_t4(0.1, 0.1, 1.0, 1.0)); + //drawDebugLine(float64_t2(_rhs, -1000.0), float64_t2(_rhs, 1000.0), float64_t4(0.1, 1.0, 1.0, 1.0)); + //printf(std::format("(comparing MinorIdx) _lhs: {} (len: {}) _rhs: {} (len: {}) minLen: {} diff: {} ", + // _lhs, lenLhs, _rhs, lenRhs, minLen, abs(_lhs - _rhs)).c_str()); + } +#endif + + // Threshhold here for intersection points, where the MinorIdx values for the curves are + // very close but could be smaller, causing the curves to be in the wrong order + if (abs(_lhs - _rhs) < MinorPositionComparisonThreshhold * minLen) + { + // this is how you want to order the derivatives dmin/dmaj=-INF dmin/dmaj = 0 dmin/dmaj=INF + // also leverage the guarantee that `dmaj>=0` to ger numerically stable compare + // also leverage the guarantee that `dmaj>=0` to ger numerically stable compare + auto lhsQuadratic = Quadratic::constructFromBezier(*lhs.originalBezier); + auto rhsQuadratic = Quadratic::constructFromBezier(*rhs.originalBezier); + + float64_t2 lTan = lhs.originalBezier->derivative(lhs.t_start); + float64_t2 rTan = rhs.originalBezier->derivative(rhs.t_start); + _lhs = lTan[MinorIdx] * rTan[MajorIdx]; + _rhs = rTan[MinorIdx] * lTan[MajorIdx]; +#ifdef DEBUG_HATCH_VISUALLY + if (false) //(debugOutput && step == debugStep) + { + printf(std::format("(comparing tangent) lTan: {}, {} rTan: {}, {} _lhs: {} _rhs: {} abs(_lhs - _rhs): {} abs(_lhs - 0.0): {} ", + lTan.x, lTan.y, rTan.x, rTan.y, _lhs, _rhs, + abs(_lhs - _rhs), abs(_lhs - 0.0)).c_str()); + } +#endif + // negative values mess with the comparison operator when using multiplication + // they should be positive because of MajorIdx monotonicity + assert(lTan[MajorIdx] >= 0.0); + assert(rTan[MajorIdx] >= 0.0); + + if (abs(_lhs - _rhs) < TangentComparisonThreshhold) + { + float64_t2 lAcc = 2.0 * lhsQuadratic.A; + float64_t2 rAcc = 2.0 * rhsQuadratic.A; + + // In this branch, _lhs == _rhs == 0 (tangents are both 0) + if (abs(_lhs - 0.0) < TangentComparisonThreshhold) + { + // TODO https://discord.com/channels/593902898015109131/723305695046533151/1169377896658383008 + bool lTanSign = lTan[MinorIdx] >= 0.0; + bool rTanSign = rTan[MinorIdx] >= 0.0; + // CASE A + // + // If the signs of the horizontal tangents differ, we know thje negative + // one belongs to the left curve + if (lTanSign != rTanSign) + { +#ifdef DEBUG_HATCH_VISUALLY + if (false) //(debugOutput && step == debugStep) + { + printf(std::format("(comparing sign) lTanSign: {} rTanSign: {} ", + lTanSign ? "positive" : "negative", rTanSign ? "positive" : "negative").c_str()); + printf("\n"); + } +#endif + // We want to return true if lhs < rhs (lhs is to the left of rhs) + // For this to be false, rhs would need to be the left one (and therefore negative) + return rTanSign; + } + + // Otherwise (CASE B) + // + // In this case tangents are both on the same side + // so only the magnitude / abs of the d2Major / dMinor2 is important + _lhs = -abs(lAcc[MinorIdx] * lTan[MajorIdx] - lAcc[MajorIdx] * lTan[MinorIdx]) * pow(rTan[MinorIdx], 3.0); + _rhs = -abs(rAcc[MinorIdx] * rTan[MajorIdx] - rAcc[MajorIdx] * rTan[MinorIdx]) * pow(lTan[MinorIdx], 3.0); +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && step == debugStep) + { + printf(std::format("(comparing second derivatives) A1={},{} B1={},{} A2={},{} B2={},{} _lhs={} _rhs={}", + lhsQuadratic.A.x, lhsQuadratic.A.y, lhsQuadratic.B.x, lhsQuadratic.B.y, + rhsQuadratic.A.x, rhsQuadratic.A.y, rhsQuadratic.B.x, rhsQuadratic.B.y, + _lhs, _rhs + ).c_str()); + } +#endif + } + else + { + _lhs = (lAcc[MinorIdx] * lTan[MajorIdx] - lAcc[MajorIdx] * lTan[MinorIdx]) * pow(rTan[MajorIdx], 3.0); + _rhs = (rAcc[MinorIdx] * rTan[MajorIdx] - rAcc[MajorIdx] * rTan[MinorIdx]) * pow(lTan[MajorIdx], 3.0); + } + } + } +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && step == debugStep) + { + printf("\n"); + } +#endif + return _lhs < _rhs; + }; + auto addToCandidateSet = [&](const Segment& entry) + { + if (isSegmentStraightLineConstantMajor(entry)) + return; + // Look for intersections among active candidates + // this is a little O(n^2) but only in the `n=candidates.size()` + for (const auto& segment : activeCandidates) + { + // find intersections entry vs segment + auto intersectionPoints = intersectSegments(entry, segment); +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && step == debugStep) + { + for (uint32_t i = 0; i < intersectionPoints.size(); i++) + { + if (nbl::core::isnan(intersectionPoints[i])) + continue; + auto point = segment.originalBezier->evaluate(intersectionPoints[i]); + auto min = point - 0.3; + auto max = point + 0.3; + drawDebugLine(float64_t2(min.x, min.y), float64_t2(max.x, min.y), float32_t4(0.0, 0.3, 0.0, 0.8)); + drawDebugLine(float64_t2(max.x, min.y), float64_t2(max.x, max.y), float32_t4(0.0, 0.3, 0.0, 0.8)); + drawDebugLine(float64_t2(min.x, max.y), float64_t2(max.x, max.y), float32_t4(0.0, 0.3, 0.0, 0.8)); + drawDebugLine(float64_t2(min.x, min.y), float64_t2(min.x, max.y), float32_t4(0.0, 0.3, 0.0, 0.8)); + } + } +#endif + + for (uint32_t i = 0; i < intersectionPoints.size(); i++) + { + if (nbl::core::isnan(intersectionPoints[i])) + continue; + intersections.push(segment.originalBezier->evaluate(intersectionPoints[i])[MajorIdx]); + } + } + activeCandidates.push_back(entry); + }; + + double lastMajor = starts.top().originalBezier->evaluate(starts.top().t_start)[MajorIdx]; + while (lastMajor!=maxMajor) + { +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && step > debugStep) + break; + bool isCurrentDebugStep = step == debugStep; +#endif + + double newMajor; + bool addStartSegmentToCandidates = false; + + if (ends.empty()) + { + m_logger.log("Hatch Creation Failure: `ends` stack is empty in the main loop", nbl::system::ILogger::ELL_ERROR); + _NBL_DEBUG_BREAK_IF(true); // This shouldn't happen, TODO: LOG + break; + } + const double maxMajorEnds = ends.top(); + + const Segment nextStartEvent = starts.empty() ? Segment() : starts.top(); + const double minMajorStart = nextStartEvent.originalBezier ? nextStartEvent.originalBezier->evaluate(nextStartEvent.t_start)[MajorIdx] : 0.0; + + // We check which event, within start, end and intersection events have the smallest + // MajorIdx coordinate at this point + auto intersectionVisit = [&]() + { + const double newMajor = intersections.top(); +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && isCurrentDebugStep) + drawDebugLine(float64_t2(-1000.0, newMajor), float64_t2(1000.0, newMajor), float32_t4(0.0, 0.0, 0.8, 1.0)); +#endif + intersections.pop(); // O(n) + return newMajor; + }; + + // next start event is before next end event + if (nextStartEvent.originalBezier && minMajorStart < maxMajorEnds) + { + // next start event is before next intersection event + // (start event) + if (intersections.empty() || minMajorStart < intersections.top()) // priority queue top() is O(1) + { + starts.pop(); + newMajor = minMajorStart; + addStartSegmentToCandidates = true; +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && isCurrentDebugStep) + drawDebugLine(float64_t2(-1000.0, newMajor), float64_t2(1000.0, newMajor), float32_t4(0.0, 0.8, 0.0, 1.0)); +#endif + } + // (intersection event) + else newMajor = intersectionVisit(); + } + // next intersection event is before next end event + // (intersection event) + else if (!intersections.empty() && intersections.top() < maxMajorEnds) + newMajor = intersectionVisit(); + else + { + // (end event) + newMajor = maxMajorEnds; + ends.pop(); +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && isCurrentDebugStep) + drawDebugLine(float64_t2(-1000.0, newMajor), float64_t2(1000.0, newMajor), float32_t4(0.0, 0.0, 0.8, 1.0)); +#endif + //std::cout << "End event at " << newMajor << "\n"; + } + // spawn quads for the previous iterations if we advanced + + if (newMajor > lastMajor) + { + const auto candidatesSize = std::distance(activeCandidates.begin(),activeCandidates.end()); + // Because n4ce works on loops, this must be `true` in almost every case, but can fail at times, because we skip adding beziers (lines) almost constant in MajorIdx direction + if (candidatesSize % 2u == 0u) + { +#ifdef DEBUG_HATCH_VISUALLY + if (debugOutput && isCurrentDebugStep) + drawDebugLine(float64_t2(-1000.0, lastMajor), float64_t2(1000.0, lastMajor), float32_t4(0.1, 0.1, 0.0, 0.5)); +#endif + // trim + if ((candidatesSize % 2u) != 0u) + { + m_logger.log("Hatch Creation Failure: candidatesSize is odd", nbl::system::ILogger::ELL_ERROR); + _NBL_DEBUG_BREAK_IF(true); // input polyline/polygon + } +#ifdef DEBUG_HATCH_VISUALLY + if (candidatesSize % 2u == 1u) + { + for (uint32_t i = 0u; i < candidatesSize; i++) + { + const Segment& item = activeCandidates[i]; + auto curveMinEnd = intersectOrtho(*item.originalBezier, newMajor, MajorIdx); + auto splitCurveMin = *item.originalBezier; + splitCurveMin.splitCurveFromMinToMax(item.t_start, nbl::core::isnan(curveMinEnd) ? 1.0 : curveMinEnd); + + drawDebugBezier(splitCurveMin, (i == candidatesSize - 1) ? float32_t4(0.0, 0.0, 1.0, 1.0) : float32_t4(1.0, 0.0, 0.0, 1.0)); + if (i == candidatesSize - 1) + { + printf(std::format("problematic guy: ({}, {}), ({}, {}), ({}, {})", + splitCurveMin.P0.x, splitCurveMin.P0.y, + splitCurveMin.P1.x, splitCurveMin.P1.y, + splitCurveMin.P2.x, splitCurveMin.P2.y + ).c_str()); + } + } + } +#endif + for (auto i = 0u; i < (candidatesSize / 2) * 2;) + { + const Segment& left = activeCandidates[i++]; + const Segment& right = activeCandidates[i++]; + + CurveBox curveBox = {}; + + // Due to precision, if the curve is right at the end, intersectOrtho may return nan + auto curveMinEnd = intersectOrtho(*left.originalBezier, newMajor, MajorIdx); + auto curveMaxEnd = intersectOrtho(*right.originalBezier, newMajor, MajorIdx); + + auto splitCurveMin = *left.originalBezier; + splitCurveMin.splitFromMinToMax(left.t_start, nbl::core::isnan(curveMinEnd) ? 1.0 : curveMinEnd); + auto splitCurveMax = *right.originalBezier; + splitCurveMax.splitFromMinToMax(right.t_start, nbl::core::isnan(curveMaxEnd) ? 1.0 : curveMaxEnd); + + assert(splitCurveMin.evaluate(0.0)[MajorIdx] <= splitCurveMin.evaluate(1.0)[MajorIdx]); + assert(splitCurveMax.evaluate(0.0)[MajorIdx] <= splitCurveMax.evaluate(1.0)[MajorIdx]); + + auto curveMinAabb = getBezierBoundingBoxMinor(splitCurveMin); + auto curveMaxAabb = getBezierBoundingBoxMinor(splitCurveMax); + curveBox.aabbMin = float64_t2(std::min(curveMinAabb.first.x, curveMaxAabb.first.x), lastMajor); + curveBox.aabbMax = float64_t2(std::max(curveMinAabb.second.x, curveMaxAabb.second.x), newMajor); + +#ifdef DEBUG_HATCH_VISUALLY + if (isCurrentDebugStep) + { + drawDebugBezier(splitCurveMin, float64_t4(1.0, 0.0, 0.0, 1.0)); + drawDebugBezier(splitCurveMax, float64_t4(0.0, 1.0, 0.0, 1.0)); + + printf(std::format("AABB min: {}, {} max: {}, {} curve min: ({}, {}), ({}, {}), ({}, {}) curve max ({}, {}), ({}, {}), ({}, {})\n", + curveBox.aabbMin.x, curveBox.aabbMin.y, curveBox.aabbMax.x, curveBox.aabbMax.y, + + splitCurveMin.P0.x, splitCurveMin.P0.y, + splitCurveMin.P1.x, splitCurveMin.P1.y, + splitCurveMin.P2.x, splitCurveMin.P2.y, + splitCurveMax.P0.x, splitCurveMax.P0.y, + splitCurveMax.P1.x, splitCurveMax.P1.y, + splitCurveMax.P2.x, splitCurveMax.P2.y + ).c_str()); + } +#endif + + curveBox.minCurve = splitCurveMin; + curveBox.maxCurve = splitCurveMax; + consumeCurveBox(curveBox); + } + } + + // advance and trim all of the beziers in the candidate set + auto oit = activeCandidates.begin(); + for (auto iit = activeCandidates.begin(); iit != activeCandidates.end(); iit++) + { + const double evalAtMajor = iit->originalBezier->evaluate(iit->t_end)[MajorIdx]; + + auto origBez = iit->originalBezier; + // if we scrolled past the end of the segment, remove it + // (basically, we memcpy everything after something is different + // and we skip on the memcpy for any items that are also different) + // (this is supposedly a pattern with input/output operators) + if (newMajor < evalAtMajor) + { + const double new_t_start = intersectOrtho(*iit->originalBezier, newMajor, MajorIdx); + + // little optimization (don't memcpy anything before something was removed) + if (oit != iit) + *oit = *iit; + oit->t_start = new_t_start; + oit++; + } + } + // trim + const auto newSize = std::distance(activeCandidates.begin(), oit); + activeCandidates.resize(newSize); + } + + // If we had a start event, we need to add the candidate + if (addStartSegmentToCandidates) + { + addToCandidateSet(nextStartEvent); + } + + // We'll need to sort if we had a start event and added to the candidate set + // or if we have advanced our candidate set + if (addStartSegmentToCandidates || newMajor > lastMajor) + { + std::sort(activeCandidates.begin(), activeCandidates.end(), candidateComparator); + } + + if (newMajor > lastMajor) + lastMajor = newMajor; + +#ifdef DEBUG_HATCH_VISUALLY + step++; +#endif + } +#ifdef DEBUG_HATCH_VISUALLY + debugStep = debugStep - step; +#endif + + } + + static bool isLineSegment(const QuadraticBezier& bezier) + { + auto quadratic = Quadratic::constructFromBezier(bezier); + float_t lenSqA = dot(quadratic.A, quadratic.A); + return lenSqA < exp(-23.0f) * dot(quadratic.B, quadratic.B); + } + +private: + + float_t intersectOrtho(const QuadraticBezier& bezier, float_t lineConstant, int component) + { + // https://pomax.github.io/bezierinfo/#intersections + float_t points[3]; + points[0] = bezier.P0[component]; + points[1] = bezier.P1[component]; + points[2] = bezier.P2[component]; + + for (uint32_t i = 0; i < 3; i++) + points[i] -= lineConstant; + + float_t A = points[0] - 2.0 * points[1] + points[2]; + float_t B = 2.0 * (points[1] - points[0]); + float_t C = points[0]; + + float_t2 roots = nbl::hlsl::math::equations::Quadratic::construct(A, B, C).computeRoots(); + if (roots.x >= 0.0 && roots.x <= 1.0) return roots.x; + if (roots.y >= 0.0 && roots.y <= 1.0) return roots.y; + return nbl::core::nan(); + } + + // returns two possible values of t in the lhs curve where the curves intersect + std::array bezierBezierIntersections(const QuadraticBezier& lhs, const QuadraticBezier& rhs) + { + const auto quarticEquation = nbl::hlsl::shapes::getBezierBezierIntersectionEquation(lhs, rhs); + + using nbl::hlsl::math::equations::Quartic; + using nbl::hlsl::math::equations::Cubic; + using nbl::hlsl::math::equations::Quadratic; + constexpr float_t QUARTIC_THRESHHOLD = 1e-10; + + std::array t = { nbl::core::nan(), nbl::core::nan(), nbl::core::nan(), nbl::core::nan() }; // only two candidates in range, ever + + const float_t quadCoeffMag = std::max(std::abs(quarticEquation.d), std::abs(quarticEquation.e)); + const float_t cubCoeffMag = std::max(std::abs(quarticEquation.c), quadCoeffMag); + const float_t quartCoeffMag = std::max(std::abs(quarticEquation.b), cubCoeffMag); + + if (std::abs(quarticEquation.a) > quartCoeffMag * QUARTIC_THRESHHOLD) + { + auto res = quarticEquation.computeRoots(); + memcpy(&t[0], &res.x, sizeof(float_t) * 4); + } + else if (abs(quarticEquation.b) > quadCoeffMag * QUARTIC_THRESHHOLD) + { + auto res = Cubic::construct(quarticEquation.b, quarticEquation.c, quarticEquation.d, quarticEquation.e).computeRoots(); + memcpy(&t[0], &res.x, sizeof(float_t) * 3); + } + else + { + auto res = Quadratic::construct(quarticEquation.c, quarticEquation.d, quarticEquation.e).computeRoots(); + memcpy(&t[0], &res.x, sizeof(float_t) * 2); + } + + // TODO: why did we do this? + // if (t[0] == t[1] || nbl::core::isnan(t[0]) || nbl::core::isnan(t[1])) + // t[0] = (t[0] != 0.0) ? 0.0 : 1.0; + + return t; + } + + bool splitIntoMajorMonotonicSegments(const QuadraticBezier& bezier, std::array, 2>& out) + { + auto quadratic = Quadratic::constructFromBezier(bezier); + + // Getting derivatives for our quadratic bezier + float_t a = quadratic.A[MajorIdx]; + float_t b = quadratic.B[MajorIdx]; + + if (a == 0.0) // would cause inf and return false for a simple linear bezier + return true; + + // Calculus 101: Finding roots for the quadratic bezier derivatives, this give us the t where the bezier changes direction, it almost always does, but if it's outside [0,1] range, then we consider it monotonic + auto t = -b / (2.0 * a); + if (t <= 0.0 || t >= 1.0) return true; + QuadraticBezier lower = bezier; lower.splitFromStart(t); + QuadraticBezier upper = bezier; upper.splitToEnd(t); + out = {lower, upper}; + return false; + } + + // https://pomax.github.io/bezierinfo/#boundingbox + std::pair getBezierBoundingBoxMinor(const QuadraticBezier& bezier) + { + float_t A = bezier.P0[MinorIdx] - 2.0 * bezier.P1[MinorIdx] + bezier.P2[MinorIdx]; + float_t B = 2.0 * (bezier.P1[MinorIdx] - bezier.P0[MinorIdx]); + + const int searchTSize = 3; + float_t searchT[searchTSize]; + searchT[0] = 0.0; + searchT[1] = 1.0; + searchT[2] = -B / (2 * A); + + float64_t2 min = float64_t2(std::numeric_limits::infinity()); + float64_t2 max = float64_t2(-std::numeric_limits::infinity()); + + for (uint32_t i = 0; i < searchTSize; i++) + { + float_t t = searchT[i]; + if (t < 0.0 || t > 1.0 || nbl::core::isnan(t)) + continue; + float64_t2 value = bezier.evaluate(t); + min = float64_t2(std::min(min.x, value.x), std::min(min.y, value.y)); + max = float64_t2(std::max(max.x, value.x), std::max(max.y, value.y)); + } + + return std::pair(min, max); + } + + // checks if it's a straight line e.g. if you're sweeping along y axis the it's a line parallel to x + static bool isSegmentStraightLineConstantMajor(const Segment& segment) + { + const float_t p0 = segment.originalBezier->P0[MajorIdx], + p1 = segment.originalBezier->P1[MajorIdx], + p2 = segment.originalBezier->P2[MajorIdx]; + //assert(p0 <= p1 && p1 <= p2); (PRECISION ISSUES ARISE ONCE MORE) + return abs(p1 - p0) <= nbl::core::exp2(-24.0) && abs(p2 - p0) <= nbl::hlsl::exp(-24.0f); + } + + std::array intersectSegments(const Segment& lhs, const Segment& rhs) + { + std::array result = { nbl::core::nan(), nbl::core::nan() }; + int resultIdx = 0; + + // Use line intersections if one or both of the beziers are linear (a = 0) + const bool selfLinear = isLineSegment(*lhs.originalBezier); + const bool otherLinear = isLineSegment(*rhs.originalBezier); + if (selfLinear && otherLinear) + { + // Line/line intersection + //TODO: use cpp-compat hlsl builtin + auto intersectionPoint = nbl::hlsl::shapes::util::LineLineIntersection( + lhs.originalBezier->P0, lhs.originalBezier->P2 - lhs.originalBezier->P0, + rhs.originalBezier->P0, rhs.originalBezier->P2 - rhs.originalBezier->P0 + ); + const float_t x1 = lhs.originalBezier->P0.x, y1 = lhs.originalBezier->P0.y, + x2 = lhs.originalBezier->P2.x, y2 = lhs.originalBezier->P2.y, + x3 = rhs.originalBezier->P0.x, y3 = rhs.originalBezier->P0.y, + x4 = rhs.originalBezier->P2.x, y4 = rhs.originalBezier->P2.y; + + // Return if point is on the lines + if (std::min(x1, x2) <= intersectionPoint.x && y1 <= intersectionPoint.y && std::max(x1, x2) >= intersectionPoint.x && y2 >= intersectionPoint.y && + std::min(x3, x4) <= intersectionPoint.x && y3 <= intersectionPoint.y && std::max(x3, x4) >= intersectionPoint.x && y4 >= intersectionPoint.y) + { + // Gets t for "other" by using intersectOrtho + auto otherT = intersectOrtho(*rhs.originalBezier, intersectionPoint.y, MajorIdx); + auto intersectionMajor = rhs.originalBezier->evaluate(otherT)[MajorIdx]; + auto thisT = intersectOrtho(*lhs.originalBezier, intersectionMajor, MajorIdx); + + if (otherT >= rhs.t_start && otherT <= rhs.t_end && thisT >= lhs.t_start && thisT <= lhs.t_end) + { + result[0] = otherT; + } + } + } + else if (selfLinear || otherLinear) + { + // Line/curve intersection + const auto& line = selfLinear ? *lhs.originalBezier : *rhs.originalBezier; + const auto& curve = selfLinear ? *rhs.originalBezier : *lhs.originalBezier; + + float_t2 D = normalize(line.P2 - line.P0); + float_t2x2 rotation = { + {D.x, D.y}, + {-D.y, D.x} + }; + QuadraticBezier rotatedCurve = { + mul(rotation, curve.P0 - line.P0), + mul(rotation, curve.P1 - line.P0), + mul(rotation, curve.P2 - line.P0) + }; + + auto intersectionCurveT = intersectOrtho(rotatedCurve, 0, (int)SweepMajorAxis::MAJOR_Y /* Always in rotation to align with X Axis */); + auto intersectionMajor = curve.evaluate(intersectionCurveT)[MajorIdx]; + auto intersectionLineT = intersectOrtho(line, intersectionMajor, (int)MajorIdx); + + auto thisT = selfLinear ? intersectionLineT : intersectionCurveT; + auto otherT = selfLinear ? intersectionCurveT : intersectionLineT; + + if (otherT >= rhs.t_start && otherT <= rhs.t_end && thisT >= lhs.t_start && thisT <= lhs.t_end) + { + result[0] = otherT; + } + } + else + { + auto thisBezier = *lhs.originalBezier; + thisBezier.splitFromMinToMax(lhs.t_start, lhs.t_end); // to get correct P0, P1, P2 for intersection testing + + const auto p0 = thisBezier.P0; + const auto p1 = thisBezier.P1; + const auto p2 = thisBezier.P2; + const bool sideP1 = nbl::hlsl::cross2D(p2 - p0, p1 - p0) >= 0.0; + + const auto& otherBezier = *rhs.originalBezier; + const std::array intersections = bezierBezierIntersections(otherBezier, thisBezier); + + for (uint32_t i = 0; i < intersections.size(); i++) + { + auto t = intersections[i]; + if (nbl::core::isnan(t) || rhs.t_start >= t || t >= rhs.t_end) + continue; + + auto intersection = otherBezier.evaluate(t); + + // Optimization istead of doing SDF to find other T and check against bounds: + // If both P1 and the intersection point are on the same side of the P0 -> P2 line of thisBezier, it's a a valid intersection + const bool sideIntersection = nbl::hlsl::cross2D(p2 - p0, intersection - p0) >= 0.0; + if (sideP1 != sideIntersection) + continue; + + const bool duplicateT = (resultIdx > 0 && t == result[0]) || (resultIdx > 1 && t == result[1]); + if (!duplicateT) + { + if (resultIdx < 2) + { + result[resultIdx] = t; + resultIdx++; + } + else + { + _NBL_DEBUG_BREAK_IF(true); // more intersections that expected + } + } + } + } + + return result; + } + + nbl::system::logger_opt_smart_ptr m_logger; + std::vector> m_beziers; // Referenced into by the segments +}; + +} // namespace nbl::ext::csg2d + +#endif