From 4f43a32e88599ca331c5d88b9433665af327113a Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Mon, 24 Aug 2026 13:02:11 +0200 Subject: [PATCH 1/8] Added `curves.hlsl` extension --- .../nbl/builtin/hlsl/ext/curves/curves.hlsl | 898 ++++++++++++++++++ 1 file changed, 898 insertions(+) create mode 100644 include/nbl/builtin/hlsl/ext/curves/curves.hlsl 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..a4456682dc --- /dev/null +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -0,0 +1,898 @@ +#ifndef _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ +#define _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ + +#ifdef USE_NABLA + +#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); + } + + //! used in special cases when parametric curves need to find inflection point by using solvnig for root of signed curvature + virtual float_t2 computeSecondOrderDifferential(float_t t) const + { + return float_t2(num_traits::quiet_NaN(), num_traits::quiet_NaN()); + } + + // gets you the t point of inflection with errorThreshold accuracy + // the curves we deal with have at most 1 inflection point + // if there is no inflection point this function will return NaN + virtual float_t computeInflectionPoint(float_t errorThreshold) const + { + return num_traits::quiet_NaN(); + } +}; + +// 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 (hlsl::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 final : public ExplicitCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b, c; + + Parabola(float_t a, float_t b, float_t c) + : a(a), b(b), c(c) + {} + + static Parabola fromThreePoints(const float_t2& P0, const float_t2& P1, const 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 override + { + return ((a * x) + b) * x + c; + } + + float_t derivative(float_t x) const override + { + return 2.0 * a * x + b; + } +}; + +template +struct CubicCurve final : ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + using float_t4 = nbl::hlsl::portable_vector_t4; + + float_t4 X; + float_t4 Y; + + CubicCurve(const float_t4& X, const float_t4& Y) + : X(X), Y(Y) + {} + + 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 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 num_traits::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; + + CircularArc(float_t r, float_t originY, float_t startAngle, float_t sweepAngle) + : r(r), originY(originY), startAngle(startAngle), sweepAngle(sweepAngle) + {} + + // from circle center (0, -v.y) to start pos (v.x, 0) + CircularArc(float_t2 v, float_t sweepAngle) + : originY(-v.y), sweepAngle(sweepAngle) + { + r = length(v); + startAngle = getSign(v.y) * acos(v.x / r); + } + + // from circle center (0, -v.y) to start pos (v.x, 0) + CircularArc(float_t2 v) + : originY(-v.y) + { + r = length(v); + startAngle = getSign(v.y) * acos(v.x / r); + sweepAngle = -2.0 * getSign(v.y) * acos(abs(originY) / r); + } + + 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) + ); + } + +private: + 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 final : public ExplicitCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b; + ExplicitEllipse(float_t a, float_t b) + : a(a), b(b) + {} + + float_t y(float_t x) const override + { + return a * sqrt(1.0 - pow((x / b), 2.0)); + } + + float_t derivative(float_t x) const override + { + return (-a * x) / ((b * b) * sqrt(1.0 - pow((x / b), 2.0))); + } +}; + +// Centered at (0,0), aligned with x +template +struct AxisAlignedEllipse final : public ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + float_t a, b; + float_t start, end; + AxisAlignedEllipse(float_t a, float_t b, float_t start, float_t end) + : a(a), b(b), start(start), end(end) + {} + + 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 : public ParametricCurve +{ + using float_t2 = nbl::hlsl::portable_vector_t2; + + nbl::hlsl::shapes::Quadratic quadratic; + float_t offset; + + OffsettedBezier(const nbl::hlsl::shapes::QuadraticBezier& quadBezier, float_t offset) + : offset(offset) + { + quadratic = nbl::hlsl::shapes::Quadratic::constructFromBezier(quadBezier.P0, quadBezier.P1, quadBezier.P2); + } + + 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(); + } +}; + +// TODO: should it stay in this header? +#include +// declare concept +#define NBL_CONCEPT_NAME BezierAdder +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename)(typename)(int32_t)(int32_t) +#define NBL_CONCEPT_TPLT_PRM_NAMES (U)(T)(Dims)(Components) +// not the greatest syntax but works +#define NBL_CONCEPT_PARAM_0 (a,U) +#define NBL_CONCEPT_PARAM_1 (uv,vector) +#define NBL_CONCEPT_PARAM_2 (layer,uint16_t) +#define NBL_CONCEPT_PARAM_3 (data,vector) +// start concept +NBL_CONCEPT_BEGIN(4) +// need to be defined AFTER the cocnept begins +#define a NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define uv NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +#define layer NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_2 +#define data NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_3 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_EXPR)(a.template set(uv,layer,data))) +); +#undef data +#undef layer +#undef uv +#undef a +#include + +template +class QuadraticBezierFitter +{ +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; + + // TODO [Przemek]: remove + 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. + // TODO [Przemek]: add constraints on the `AddBezierFunc` type + template + 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); + } + + // TODO [Przemek]: add constraints on the `AddBezierFunc` type + template + 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(lenghtMajor, lenghtMinor, ellipse.angleBounds.x, ellipse.angleBounds.y); + adaptive(aaEllipse, 0.0, 1.0, targetMaxError, addTransformedBezier, maxDepth); + } + } + + // TODO [Przemek]: add constraints on the `AddBezierFunc` type + template + 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)); + } + } +}; + +} + +#include "Shaders/globals.hlsl" + +namespace en::nabla2d::curves +{ +// Mixes/Interpolation of two Parametric Curves t from 0 to 1 +template +struct MixedParametricCurves final : public nbl::ext::curves::ParametricCurve +{ + 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; + + using ParametricCurve = nbl::ext::curves::ParametricCurve; + + const ParametricCurve* curve1; + const ParametricCurve* curve2; + + MixedParametricCurves(const ParametricCurve* curve1, const ParametricCurve* curve2) + : curve1(curve1), curve2(curve2) + {} + + float_t2 computePosition(float_t t) const override + { + const float_t2 curve1Pos = curve1->computePosition(t); + const float_t2 curve2Pos = curve2->computePosition(t); + return t * (curve2Pos - curve1Pos) + curve1Pos; + } + + //! compute unnormalized tangent vector at t + float_t2 computeTangent(float_t t) const override + { + const float_t2 curve1Pos = curve1->computePosition(t); + const float_t2 curve2Pos = curve2->computePosition(t); + const float_t2 curve1Tan = curve1->computeTangent(t); + const float_t2 curve2Tan = curve2->computeTangent(t); + return (1 - t) * curve1Tan - curve1Pos + (t)*curve2Tan + curve2Pos; + } + + //! compute second order differential at t + float_t2 computeSecondOrderDifferential(float_t t) const override + { + const float_t2 curve1Tan = curve1->computeTangent(t); + const float_t2 curve2Tan = curve2->computeTangent(t); + const float_t2 curve1SecondDiff = curve1->computeSecondOrderDifferential(t); + const float_t2 curve2SecondDiff = curve2->computeSecondOrderDifferential(t); + return (1 - t) * curve1SecondDiff + 2.0 * (curve2Tan - curve1Tan) + t * curve2SecondDiff; + } + + float_t computeInflectionPoint(float_t errorThreshold) const override + { + auto signedCurvatureUnnormalized = [&](float_t t) + { + const float_t2 first = computeTangent(t); + const float_t2 second = computeSecondOrderDifferential(t); + return float_t(first.x * second.y - second.x * first.y); + }; + + constexpr uint16_t MaxIterations = 32u; + float_t low = 0.0; + float_t high = 1.0; + float_t valLow = signedCurvatureUnnormalized(low); + float_t valHigh = signedCurvatureUnnormalized(high); + if (getSign(valLow) != getSign(valHigh)) + { + if (valLow > valHigh) + hlsl::swap(low, high); + + float_t guess = 0.0; + for (uint16_t i = 0u; i < MaxIterations; ++i) + { + guess = (low + high) / 2.0; + float_t valGuess = signedCurvatureUnnormalized(guess); + if (abs(valGuess) < errorThreshold) + return guess; + + if (valGuess < 0.0) + low = guess; + else + high = guess; + } + + return guess; + } + else + { + return num_traits::quiet_NaN; + } + } + +private: + static float_t getSign(float_t x) + { + return static_cast((x > 0.0)) - static_cast((x <= 0.0)); + } +}; + +// Mix between two parabolas from 0 to len +template +struct MixedParabola final : public nbl::ext::curves::ExplicitCurve +{ + 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; + + float_t a, b, c, d; + + MixedParabola(const nbl::ext::curves::Parabola& parabola1, const nbl::ext::curves::Parabola& parabola2, float_t chordLen) + { + a = (parabola2.a - parabola1.a) / chordLen; + b = (parabola2.b - parabola1.b) / chordLen + parabola1.a; + c = (parabola2.c - parabola1.c) / chordLen + parabola1.b; + d = parabola1.c; + } + + static MixedParabola fromFourPoints(const float_t2& P0, const float_t2& P1, const float_t2& P2, const float_t2& P3) + { + assert(P1.x == 0); + assert(P1.y == 0 && P2.y == 0); + auto parabola1 = nbl::ext::curves::Parabola::fromThreePoints(P0, P1, P2); + auto parabola2 = nbl::ext::curves::Parabola::fromThreePoints(P1, P2, P3); + return MixedParabola(parabola1, parabola2, abs(P2.x - P1.x)); + } + + float_t y(float_t x) const override + { + return (((a * x) + b) * x + c) * x + d; + } + + float_t derivative(float_t x) const override + { + return ((3.0 * a * x) + 2.0 * b) * x + c; + } + + float_t computeInflectionPoint(float_t errorThreshold) const override + { + return -b / (3.0 * a); + } +}; + +// Centered at (0, 0), P1 and P2 on x axis and P1.x = -P2.x +template +struct ExplicitMixedCircle final : public nbl::ext::curves::ExplicitCurve +{ + 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 ExplicitCircle + { + float_t2 origin; + float_t radius; + + ExplicitCircle(const float_t2& origin, float_t radius) + : origin(origin), radius(radius) + {} + + static ExplicitCircle fromThreePoints(float_t2 P0, float_t2 P1, float_t2 P2) + { + const float_t2 Mid0 = (P0 + P1) / 2.0; + const float_t2 Normal0 = float_t2(P1.y - P0.y, P0.x - P1.x); + + const float_t2 Mid1 = (P1 + P2) / 2.0; + const float_t2 Normal1 = float_t2(P2.y - P1.y, P1.x - P2.x); + + const float_t2 origin = shapes::util::LineLineIntersection(Mid0, Normal0, Mid1, Normal1); + const float_t radius = glm::length(P0 - origin); + return ExplicitCircle(origin, radius); + } + }; + + float_t origin1Y; + float_t origin2Y; + float_t radius1; + float_t radius2; + float_t chordLen; + + static ExplicitMixedCircle fromFourPoints(const float_t2& P0, const float_t2& P1, const float_t2& P2, const float_t2& P3); + + float_t y(float_t x) const override + { + // https://herbie.uwplse.org/demo/5966aab2781d6c55c07105f5b900f1506cad301e.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html + const float_t s1 = -1.0 * getSign(origin1Y); + const float_t s2 = -1.0 * getSign(origin2Y); + + const float_t t1 = s1 * sqrt(radius1 * radius1 - x * x) + origin1Y; + const float_t t2 = s2 * sqrt(radius2 * radius2 - x * x) + origin2Y; + const float_t ret = (x / chordLen + 0.5) * (t2 - t1) + t1; + return ret; + } + + float_t secondDerivative(float_t x) const + { + // https://www.wolframalpha.com/input?i=second+derivative+of+%28x%2Fl%2B1%2F2%29*%28s2*sqrt%28r_2%5E2-x%5E2%29%2Bo2-s1*sqrt%28r_1%5E2-x%5E2%29-o1%29%2Bs1*sqrt%28r_1%5E2-x%5E2%29%2Bo1 + const float_t s1 = -1.0 * getSign(origin1Y); + const float_t s2 = -1.0 * getSign(origin2Y); + + const float_t t1 = sqrt(radius1 * radius1 - x * x); + const float_t t2 = sqrt(radius2 * radius2 - x * x); + + const float_t u1 = (s1 * x) / t1; + const float_t u2 = (s2 * x) / t2; + + const float_t q1 = (-1.0 * (x * x) / pow(radius1 * radius1 - x * x, 1.5)) - (1.0 / t1); + const float_t q2 = (-1.0 * (x * x) / pow(radius2 * radius2 - x * x, 1.5)) - (1.0 / t2); + + const float_t ret = ((2.0 * (u1 - u2)) / chordLen) + (x / chordLen + 0.5) * (s2 * q2 - s1 * q1) + s1 * q1; + return ret; + } + + float_t derivative(float_t x) const override + { + // https://www.wolframalpha.com/input?i=derivative+%28x%2Fl%2B1%2F2%29*%28s2*sqrt%28r_2%5E2-x%5E2%29%2Borigin2Y-s1*sqrt%28r_1%5E2-x%5E2%29-origin1Y%29%2Bs1*sqrt%28r_1%5E2-x%5E2%29%2Borigin1Y + // https ://herbie.uwplse.org/demo/d48eba9a858160f5f8e4283f7e6211d41215d354.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html + const float_t s1 = -1.0 * getSign(origin1Y); + const float_t s2 = -1.0 * getSign(origin2Y); + + const float_t t0 = sqrt(radius1 * radius1 - x * x); + const float_t t1 = (s1 * x) / t0; + const float_t t2 = sqrt(radius2 * radius2 - x * x); + const float_t ret = ((x / chordLen + 0.5) * (t1 - ((x * s2) / t2)) + (((origin2Y - origin1Y) + (s2 * t2) - (s1 * t0)) / chordLen)) - t1; + return ret; + } + + float_t computeInflectionPoint(float_t errorThreshold) const override + { + // bisection search to find inflection point + // by seeing the graph of second derivative over wide range of values we have deduced that the inflection point exists iff the secondDerivative has opposite signs at begin and end + constexpr uint16_t MaxIterations = 64u; + float_t low = -chordLen / 2.0; + float_t high = chordLen / 2.0; + float_t valLow = secondDerivative(low + errorThreshold / 2.0); + float_t valHigh = secondDerivative(high - errorThreshold / 2.0); + if (getSign(valLow) != getSign(valHigh)) + { + if (valLow > valHigh) + hlsl::swap(low, high); + + float_t guess = 0.0; + for (uint16_t i = 0u; i < MaxIterations; ++i) + { + guess = (low + high) / 2.0; + float_t valGuess = secondDerivative(guess); + if (abs(valGuess) < errorThreshold) + return guess; + + if (valGuess < 0.0) + low = guess; + else + high = guess; + } + + return guess; + } + else + { + return num_traits::quiet_NaN; + } + } + +private: + static float_t getSign(float_t x) + { + return static_cast((x > 0.0)) - static_cast((x <= 0.0)); + } +}; + +} // namespace nbl +} // namespace hlsl +} // namespace ext +} // namespace curves +#endif + +#endif \ No newline at end of file From 6fae648ffde79c5a39f7a2725da73015bb2af468 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Tue, 25 Aug 2026 18:34:39 +0200 Subject: [PATCH 2/8] Initial `curves.hlsl` refactor --- .../nbl/builtin/hlsl/ext/curves/curves.hlsl | 406 ++++++++++-------- 1 file changed, 236 insertions(+), 170 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index a4456682dc..97fba0f44b 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -22,47 +22,38 @@ namespace ext namespace curves { +namespace impl +{ // 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) - {} + // TODO: maybe there is a better solution + ParametricCurve curve; - inline float_t operator()(const float_t t) const + inline float_t operator()(const float_t t) { - return m_curve->differentialArcLen(t); + return curve.differentialArcLen(t); } }; //! compute arc length by gauss legendere integration - float_t arcLen(float_t t0, float_t t1) const + float_t arcLen(float_t t0, float_t t1) { - constexpr uint16_t IntegrationOrder = 10u; - return nbl::hlsl::math::quadrature::GaussLegendreIntegration::calculateIntegral(ArcLenIntegrand(this), t0, t1); + const uint16_t IntegrationOrder = 10u; + + // TODO: maybe there is a better solution + ArcLenIntegrand arcLenIntegrand; + //arcLenIntegrand.curve = *this; // TODO: copy this + return nbl::hlsl::math::quadrature::GaussLegendreIntegration::calculateIntegral(arcLenIntegrand, 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 inverseArcLen_BisectionSearch(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4, const uint16_t iterationThreshold = 16u) { float_t xi = 0.0; float_t low = min; @@ -91,65 +82,70 @@ struct ParametricCurve } //! compute inverse arc len - float_t inverseArcLen(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4) const + float_t inverseArcLen(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4) { return inverseArcLen_BisectionSearch(targetLen, min, max, cdfAccuracyThreshold); } - //! used in special cases when parametric curves need to find inflection point by using solvnig for root of signed curvature - virtual float_t2 computeSecondOrderDifferential(float_t t) const + float_t differentialArcLen(float_t t) + { + return nbl::hlsl::length(computeTangent(t)); + } + float_t2 computeSecondOrderDifferential(float_t t) { return float_t2(num_traits::quiet_NaN(), num_traits::quiet_NaN()); } - - // gets you the t point of inflection with errorThreshold accuracy - // the curves we deal with have at most 1 inflection point - // if there is no inflection point this function will return NaN - virtual float_t computeInflectionPoint(float_t errorThreshold) const + float_t computeInflectionPoint(float_t errorThreshold) { return num_traits::quiet_NaN(); } }; -// 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; +// TODO: make an `ExplicitCurve` concept, it should require for a type to define `y` and `derivative` - 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 (hlsl::isinf(deriv)) - v = float_t2(0.0, 1.0); - return v; - } +// TODO: rationale this decision +// It's when t = x in a Parametric Curve +#define DEFINE_EXPLICIT_CURVE_FUNCTIONS \ +float_t differentialArcLen(float_t x) \ +{ \ + float_t deriv = derivative(x); \ + return sqrt(1.0 + deriv * deriv); \ +} \ + \ +nbl::hlsl::portable_vector_t2 computeTangent(float_t x) \ +{ \ + const float_t deriv = derivative(x); \ + nbl::hlsl::portable_vector_t2 v = nbl::hlsl::portable_vector_t2(1.0, deriv); \ + if (nbl::hlsl::isinf(deriv)) \ + v = nbl::hlsl::portable_vector_t2(0.0, 1.0); \ + return v; \ +} \ + \ +inline nbl::hlsl::portable_vector_t2 computePosition(float_t x) { return nbl::hlsl::portable_vector_t2(x, y(x)); } - inline float_t2 computePosition(float_t x) const override { return float_t2(x, y(x)); } -}; template -struct Parabola final : public ExplicitCurve +struct Parabola : ParametricCurve { + DEFINE_EXPLICIT_CURVE_FUNCTIONS + using float_t2 = nbl::hlsl::portable_vector_t2; float_t a, b, c; - Parabola(float_t a, float_t b, float_t c) - : a(a), b(b), c(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(const float_t2& P0, const float_t2& P1, const float_t2& P2) + 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), @@ -160,19 +156,19 @@ struct Parabola final : public ExplicitCurve return Parabola(M[0], M[1], M[2]); } - float_t y(float_t x) const override + float_t y(float_t x) { return ((a * x) + b) * x + c; } - float_t derivative(float_t x) const override + float_t derivative(float_t x) { return 2.0 * a * x + b; } }; template -struct CubicCurve final : ParametricCurve +struct CubicCurve : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; using float_t4 = nbl::hlsl::portable_vector_t4; @@ -180,11 +176,16 @@ struct CubicCurve final : ParametricCurve float_t4 X; float_t4 Y; - CubicCurve(const float_t4& X, const float_t4& Y) - : X(X), Y(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 + float_t2 computePosition(float_t t) { return float_t2( ((X[0] * t + X[1]) * t + X[2]) * t + X[3], @@ -193,7 +194,7 @@ struct CubicCurve final : ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) const override + float_t2 computeTangent(float_t t) { return float_t2( (3.0 * X[0] * t + 2.0 * X[1]) * t + X[2], @@ -202,7 +203,7 @@ struct CubicCurve final : ParametricCurve } //! compute second order differential at t - float_t2 computeSecondOrderDifferential(float_t t) const override + float_t2 computeSecondOrderDifferential(float_t t) { return float_t2( 6.0 * X[0] * t + 2.0 * X[1], @@ -210,7 +211,12 @@ struct CubicCurve final : ParametricCurve ); } - float_t computeInflectionPoint(float_t errorThreshold) const override + float_t differentialArcLen(float_t t) + { + return nbl::hlsl::length(computeTangent(t)); + } + + float_t computeInflectionPoint(float_t errorThreshold) { // solve for signed curvature root // when x'*y''-x''*y' = 0 @@ -241,28 +247,44 @@ struct CircularArc : ParametricCurve float_t startAngle; float_t sweepAngle; - CircularArc(float_t r, float_t originY, float_t startAngle, float_t sweepAngle) - : r(r), originY(originY), startAngle(startAngle), sweepAngle(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) - CircularArc(float_t2 v, float_t sweepAngle) - : originY(-v.y), sweepAngle(sweepAngle) + static CircularArc create(float_t2 v, float_t sweepAngle) { - r = length(v); - startAngle = getSign(v.y) * acos(v.x / r); + 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) - CircularArc(float_t2 v) - : originY(-v.y) + static CircularArc create(float_t2 v) { - r = length(v); - startAngle = getSign(v.y) * acos(v.x / r); - sweepAngle = -2.0 * getSign(v.y) * acos(abs(originY) / r); + 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 + float_t2 computePosition(float_t t) { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -272,7 +294,7 @@ struct CircularArc : ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) const override + float_t2 computeTangent(float_t t) { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -281,7 +303,7 @@ struct CircularArc : ParametricCurve ); } - float_t2 computeSecondOrderDifferential(float_t t) const override + float_t2 computeSecondOrderDifferential(float_t t) { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -290,7 +312,6 @@ struct CircularArc : ParametricCurve ); } -private: static float_t getSign(float_t x) { return static_cast((x > 0.0)) - static_cast((x <= 0.0)); @@ -299,21 +320,28 @@ private: // Centered at (0,0), aligned with x axis template -struct ExplicitEllipse final : public ExplicitCurve +struct ExplicitEllipse : ParametricCurve { + DEFINE_EXPLICIT_CURVE_FUNCTIONS + using float_t2 = nbl::hlsl::portable_vector_t2; float_t a, b; - ExplicitEllipse(float_t a, float_t b) - : a(a), b(b) - {} + static ExplicitEllipse create(float_t a, float_t b) + { + ExplicitEllipse output; + output.a = a; + output.b = b; - float_t y(float_t x) const override + return output; + } + + float_t y(float_t x) { return a * sqrt(1.0 - pow((x / b), 2.0)); } - float_t derivative(float_t x) const override + float_t derivative(float_t x) { return (-a * x) / ((b * b) * sqrt(1.0 - pow((x / b), 2.0))); } @@ -321,23 +349,30 @@ struct ExplicitEllipse final : public ExplicitCurve // Centered at (0,0), aligned with x template -struct AxisAlignedEllipse final : public ParametricCurve +struct AxisAlignedEllipse : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; float_t a, b; float_t start, end; - AxisAlignedEllipse(float_t a, float_t b, float_t start, float_t end) - : a(a), b(b), start(start), end(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; - float_t2 computePosition(float_t t) const override + return output; + } + + float_t2 computePosition(float_t t) { 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 + float_t2 computeTangent(float_t t) { const float_t theta = start + (end - start) * t; const float_t dThetaDt = end - start; @@ -355,7 +390,7 @@ struct EllipticalArcInfo float_t2 angleBounds; // [0, 2Pi) float_t eccentricity; // (0, 1] - inline bool isValid() const + inline bool isValid() { if (eccentricity > 1.0 || eccentricity <= 0.0) return false; @@ -368,20 +403,24 @@ struct EllipticalArcInfo }; template -struct OffsettedBezier : public ParametricCurve +struct OffsettedBezier : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; nbl::hlsl::shapes::Quadratic quadratic; float_t offset; - OffsettedBezier(const nbl::hlsl::shapes::QuadraticBezier& quadBezier, float_t offset) - : offset(offset) + static OffsettedBezier create(NBL_CONST_REF_ARG(nbl::hlsl::shapes::QuadraticBezier) quadBezier, float_t offset) { - quadratic = nbl::hlsl::shapes::Quadratic::constructFromBezier(quadBezier.P0, quadBezier.P1, quadBezier.P2); + 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 + float_t2 computePosition(float_t t) { const float_t2 deriv = quadratic.derivative(t); const float_t2 normal = normalize(float_t2(deriv.y, -deriv.x)); @@ -389,7 +428,7 @@ struct OffsettedBezier : public ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) const override + float_t2 computeTangent(float_t t) { const float_t2 ddt = quadratic.derivative(t); const float_t2 d2dt2 = quadratic.secondDerivative(t); @@ -398,7 +437,7 @@ struct OffsettedBezier : public ParametricCurve } //! if offset is more than minimum radius of curvature then we get an unwanted gouging/cusp - float_t2 findCusps() const + float_t2 findCusps() { // 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); @@ -438,24 +477,20 @@ NBL_CONCEPT_END( #include template -class QuadraticBezierFitter +struct QuadraticBezierFitter { -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; - // TODO [Przemek]: remove - 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. // TODO [Przemek]: add constraints on the `AddBezierFunc` type template - static void adaptive(const ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + static void adaptive(NBL_CONST_REF_ARG(ParametricCurve) curve, float_t min, float_t max, float_t targetMaxError, NBL_REF_ARG(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) @@ -470,7 +505,7 @@ public: // TODO [Przemek]: add constraints on the `AddBezierFunc` type template - static void adaptive(const EllipticalArcInfo& ellipse, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + static void adaptive(NBL_CONST_REF_ARG(EllipticalArcInfo) ellipse, float_t targetMaxError, NBL_REF_ARG(AddBezierFunc) addBezierFunc, uint32_t maxDepth = 12) { using namespace nbl::hlsl; @@ -533,7 +568,7 @@ public: } } -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 @@ -549,6 +584,7 @@ private: } } + template static void adaptive_impl(const ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t depth) { if (min == max) @@ -620,13 +656,22 @@ private: } +} // namespace curves +} // namespace ext +} // namespace hlsl +} // namespace nbl + #include "Shaders/globals.hlsl" -namespace en::nabla2d::curves +namespace en +{ +namespace nabla2d +{ +namespace curves { // Mixes/Interpolation of two Parametric Curves t from 0 to 1 template -struct MixedParametricCurves final : public nbl::ext::curves::ParametricCurve +struct MixedParametricCurves : nbl::ext::curves::ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; using float_t3 = nbl::hlsl::portable_vector_t3; @@ -635,50 +680,56 @@ struct MixedParametricCurves final : public nbl::ext::curves::ParametricCurve; - const ParametricCurve* curve1; - const ParametricCurve* curve2; + // TODO: rememeber these used to be pointers + ParametricCurve curve1; + ParametricCurve curve2; - MixedParametricCurves(const ParametricCurve* curve1, const ParametricCurve* curve2) - : curve1(curve1), curve2(curve2) - {} + static MixedParametricCurves create(NBL_CONST_REF_ARG(ParametricCurve) curve1, NBL_CONST_REF_ARG(ParametricCurve) curve2) + { + MixedParametricCurves output; + output.curve1 = curve1; + output.curve2 = curve2; + + return output; + } - float_t2 computePosition(float_t t) const override + float_t2 computePosition(float_t t) { - const float_t2 curve1Pos = curve1->computePosition(t); - const float_t2 curve2Pos = curve2->computePosition(t); + const float_t2 curve1Pos = curve1.computePosition(t); + const float_t2 curve2Pos = curve2.computePosition(t); return t * (curve2Pos - curve1Pos) + curve1Pos; } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) const override + float_t2 computeTangent(float_t t) { - const float_t2 curve1Pos = curve1->computePosition(t); - const float_t2 curve2Pos = curve2->computePosition(t); - const float_t2 curve1Tan = curve1->computeTangent(t); - const float_t2 curve2Tan = curve2->computeTangent(t); + const float_t2 curve1Pos = curve1.computePosition(t); + const float_t2 curve2Pos = curve2.computePosition(t); + const float_t2 curve1Tan = curve1.computeTangent(t); + const float_t2 curve2Tan = curve2.computeTangent(t); return (1 - t) * curve1Tan - curve1Pos + (t)*curve2Tan + curve2Pos; } //! compute second order differential at t - float_t2 computeSecondOrderDifferential(float_t t) const override + float_t2 computeSecondOrderDifferential(float_t t) { - const float_t2 curve1Tan = curve1->computeTangent(t); - const float_t2 curve2Tan = curve2->computeTangent(t); - const float_t2 curve1SecondDiff = curve1->computeSecondOrderDifferential(t); - const float_t2 curve2SecondDiff = curve2->computeSecondOrderDifferential(t); + const float_t2 curve1Tan = curve1.computeTangent(t); + const float_t2 curve2Tan = curve2.computeTangent(t); + const float_t2 curve1SecondDiff = curve1.computeSecondOrderDifferential(t); + const float_t2 curve2SecondDiff = curve2.computeSecondOrderDifferential(t); return (1 - t) * curve1SecondDiff + 2.0 * (curve2Tan - curve1Tan) + t * curve2SecondDiff; } - float_t computeInflectionPoint(float_t errorThreshold) const override + float_t signedCurvatureUnnormalized(float_t t) { - auto signedCurvatureUnnormalized = [&](float_t t) - { - const float_t2 first = computeTangent(t); - const float_t2 second = computeSecondOrderDifferential(t); - return float_t(first.x * second.y - second.x * first.y); - }; + const float_t2 first = computeTangent(t); + const float_t2 second = computeSecondOrderDifferential(t); + return float_t(first.x * second.y - second.x * first.y); + }; - constexpr uint16_t MaxIterations = 32u; + float_t computeInflectionPoint(float_t errorThreshold) + { + const uint16_t MaxIterations = 32u; float_t low = 0.0; float_t high = 1.0; float_t valLow = signedCurvatureUnnormalized(low); @@ -710,17 +761,18 @@ struct MixedParametricCurves final : public nbl::ext::curves::ParametricCurve((x > 0.0)) - static_cast((x <= 0.0)); + return _static_cast((x > 0.0)) - _static_cast((x <= 0.0)); } }; // Mix between two parabolas from 0 to len template -struct MixedParabola final : public nbl::ext::curves::ExplicitCurve +struct MixedParabola : nbl::ext::hlsl::curves::ExplicitCurve { + DEFINE_EXPLICIT_CURVE_FUNCTIONS + using float_t2 = nbl::hlsl::portable_vector_t2; using float_t3 = nbl::hlsl::portable_vector_t3; using float_t4 = nbl::hlsl::portable_vector_t4; @@ -728,34 +780,37 @@ struct MixedParabola final : public nbl::ext::curves::ExplicitCurve float_t a, b, c, d; - MixedParabola(const nbl::ext::curves::Parabola& parabola1, const nbl::ext::curves::Parabola& parabola2, float_t chordLen) + static MixedParabola create(NBL_CONST_REF_ARG(nbl::ext::hlsl::curves::Parabola) parabola1, NBL_CONST_REF_ARG(nbl::ext::hlsl::curves::Parabola) parabola2, float_t chordLen) { - a = (parabola2.a - parabola1.a) / chordLen; - b = (parabola2.b - parabola1.b) / chordLen + parabola1.a; - c = (parabola2.c - parabola1.c) / chordLen + parabola1.b; - d = parabola1.c; + MixedParabola output; + output.a = (parabola2.a - parabola1.a) / chordLen; + output.b = (parabola2.b - parabola1.b) / chordLen + parabola1.a; + output.c = (parabola2.c - parabola1.c) / chordLen + parabola1.b; + output.d = parabola1.c; + + return output; } - static MixedParabola fromFourPoints(const float_t2& P0, const float_t2& P1, const float_t2& P2, const float_t2& P3) + static MixedParabola fromFourPoints(NBL_CONST_REF_ARG(float_t2) P0, NBL_CONST_REF_ARG(float_t2) P1, NBL_CONST_REF_ARG(float_t2) P2, NBL_CONST_REF_ARG(float_t2) P3) { assert(P1.x == 0); assert(P1.y == 0 && P2.y == 0); - auto parabola1 = nbl::ext::curves::Parabola::fromThreePoints(P0, P1, P2); - auto parabola2 = nbl::ext::curves::Parabola::fromThreePoints(P1, P2, P3); + nbl::hlsl::ext::curves::Parabola parabola1 = nbl::hlsl::ext::curves::Parabola::fromThreePoints(P0, P1, P2); + nbl::hlsl::ext::curves::Parabola parabola2 = nbl::hlsl::ext::curves::Parabola::fromThreePoints(P1, P2, P3); return MixedParabola(parabola1, parabola2, abs(P2.x - P1.x)); } - float_t y(float_t x) const override + float_t y(float_t x) { return (((a * x) + b) * x + c) * x + d; } - float_t derivative(float_t x) const override + float_t derivative(float_t x) { return ((3.0 * a * x) + 2.0 * b) * x + c; } - float_t computeInflectionPoint(float_t errorThreshold) const override + float_t computeInflectionPoint(float_t errorThreshold) { return -b / (3.0 * a); } @@ -763,8 +818,10 @@ struct MixedParabola final : public nbl::ext::curves::ExplicitCurve // Centered at (0, 0), P1 and P2 on x axis and P1.x = -P2.x template -struct ExplicitMixedCircle final : public nbl::ext::curves::ExplicitCurve +struct ExplicitMixedCircle : nbl::hlsl::ext::curves::ParametricCurve { + DEFINE_EXPLICIT_CURVE_FUNCTIONS + using float_t2 = nbl::hlsl::portable_vector_t2; using float_t3 = nbl::hlsl::portable_vector_t3; using float_t4 = nbl::hlsl::portable_vector_t4; @@ -775,9 +832,14 @@ struct ExplicitMixedCircle final : public nbl::ext::curves::ExplicitCurve((x > 0.0)) - static_cast((x <= 0.0)); + return _static_cast((x > 0.0)) - _static_cast((x <= 0.0)); } }; -} // namespace nbl -} // namespace hlsl -} // namespace ext +#undef DEFINE_EXPLICIT_CURVE_FUNCTIONS + } // namespace curves +} // namespace nabla2d +} // namespace en + #endif #endif \ No newline at end of file From b4c82bad9caf60fbaab3179d661ec1f5ab04f973 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Tue, 25 Aug 2026 19:32:39 +0200 Subject: [PATCH 3/8] Reverted some changes --- .../nbl/builtin/hlsl/ext/curves/curves.hlsl | 182 +++++++++--------- 1 file changed, 90 insertions(+), 92 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index 97fba0f44b..819e9703ab 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -1,7 +1,8 @@ #ifndef _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ #define _NBL_BUILTIN_HLSL_EXT_CURVES_CURVES_INCLUDED_ -#ifdef USE_NABLA +// TODO: make this file HLSL compatible +#ifndef __HLSL_VERSION #include #include @@ -30,26 +31,37 @@ 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 { - // TODO: maybe there is a better solution - ParametricCurve curve; + const ParametricCurve* m_curve; + + ArcLenIntegrand(const ParametricCurve* curve) + : m_curve(curve) + {} - inline float_t operator()(const float_t t) + inline float_t operator()(const float_t t) const { - return curve.differentialArcLen(t); + return m_curve->differentialArcLen(t); } }; //! compute arc length by gauss legendere integration - float_t arcLen(float_t t0, float_t t1) + float_t arcLen(float_t t0, float_t t1) const { - const uint16_t IntegrationOrder = 10u; - - // TODO: maybe there is a better solution - ArcLenIntegrand arcLenIntegrand; - //arcLenIntegrand.curve = *this; // TODO: copy this - return nbl::hlsl::math::quadrature::GaussLegendreIntegration::calculateIntegral(arcLenIntegrand, t0, t1); + constexpr uint16_t IntegrationOrder = 10u; + return nbl::hlsl::math::quadrature::GaussLegendreIntegration::calculateIntegral(ArcLenIntegrand(this), t0, t1); } //! compute inverse arc len using bisection search @@ -87,15 +99,12 @@ struct ParametricCurve return inverseArcLen_BisectionSearch(targetLen, min, max, cdfAccuracyThreshold); } - float_t differentialArcLen(float_t t) - { - return nbl::hlsl::length(computeTangent(t)); - } - float_t2 computeSecondOrderDifferential(float_t t) + virtual float_t2 computeSecondOrderDifferential(float_t t) const { return float_t2(num_traits::quiet_NaN(), num_traits::quiet_NaN()); } - float_t computeInflectionPoint(float_t errorThreshold) + + virtual float_t computeInflectionPoint(float_t errorThreshold) const { return num_traits::quiet_NaN(); } @@ -105,16 +114,15 @@ struct ParametricCurve // TODO: make an `ExplicitCurve` concept, it should require for a type to define `y` and `derivative` -// TODO: rationale this decision // It's when t = x in a Parametric Curve #define DEFINE_EXPLICIT_CURVE_FUNCTIONS \ -float_t differentialArcLen(float_t x) \ +float_t differentialArcLen(float_t x) const override\ { \ float_t deriv = derivative(x); \ return sqrt(1.0 + deriv * deriv); \ } \ \ -nbl::hlsl::portable_vector_t2 computeTangent(float_t x) \ +nbl::hlsl::portable_vector_t2 computeTangent(float_t x) const override\ { \ const float_t deriv = derivative(x); \ nbl::hlsl::portable_vector_t2 v = nbl::hlsl::portable_vector_t2(1.0, deriv); \ @@ -123,11 +131,11 @@ nbl::hlsl::portable_vector_t2 computeTangent(float_t x) \ return v; \ } \ \ -inline nbl::hlsl::portable_vector_t2 computePosition(float_t x) { return nbl::hlsl::portable_vector_t2(x, y(x)); } +inline nbl::hlsl::portable_vector_t2 computePosition(float_t x) const override { return nbl::hlsl::portable_vector_t2(x, y(x)); } template -struct Parabola : ParametricCurve +struct Parabola : impl::ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -156,19 +164,19 @@ struct Parabola : ParametricCurve return Parabola(M[0], M[1], M[2]); } - float_t y(float_t x) + float_t y(float_t x) const { return ((a * x) + b) * x + c; } - float_t derivative(float_t x) + float_t derivative(float_t x) const { return 2.0 * a * x + b; } }; template -struct CubicCurve : ParametricCurve +struct CubicCurve : impl::ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; using float_t4 = nbl::hlsl::portable_vector_t4; @@ -185,7 +193,7 @@ struct CubicCurve : ParametricCurve return output; } - float_t2 computePosition(float_t t) + float_t2 computePosition(float_t t) const override { return float_t2( ((X[0] * t + X[1]) * t + X[2]) * t + X[3], @@ -194,7 +202,7 @@ struct CubicCurve : ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) + float_t2 computeTangent(float_t t) const override { return float_t2( (3.0 * X[0] * t + 2.0 * X[1]) * t + X[2], @@ -203,7 +211,7 @@ struct CubicCurve : ParametricCurve } //! compute second order differential at t - float_t2 computeSecondOrderDifferential(float_t t) + float_t2 computeSecondOrderDifferential(float_t t) const override { return float_t2( 6.0 * X[0] * t + 2.0 * X[1], @@ -211,12 +219,12 @@ struct CubicCurve : ParametricCurve ); } - float_t differentialArcLen(float_t t) + float_t differentialArcLen(float_t t) const override { return nbl::hlsl::length(computeTangent(t)); } - float_t computeInflectionPoint(float_t errorThreshold) + float_t computeInflectionPoint(float_t errorThreshold) const override { // solve for signed curvature root // when x'*y''-x''*y' = 0 @@ -238,7 +246,7 @@ struct CubicCurve : ParametricCurve // 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 +struct CircularArc : impl::ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -284,7 +292,7 @@ struct CircularArc : ParametricCurve return output; } - float_t2 computePosition(float_t t) + float_t2 computePosition(float_t t) const override { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -294,7 +302,7 @@ struct CircularArc : ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) + float_t2 computeTangent(float_t t) const override { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -303,7 +311,7 @@ struct CircularArc : ParametricCurve ); } - float_t2 computeSecondOrderDifferential(float_t t) + float_t2 computeSecondOrderDifferential(float_t t) const override { const float_t actualT = t * sweepAngle + startAngle; return float_t2( @@ -320,7 +328,7 @@ struct CircularArc : ParametricCurve // Centered at (0,0), aligned with x axis template -struct ExplicitEllipse : ParametricCurve +struct ExplicitEllipse : impl::ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -336,12 +344,12 @@ struct ExplicitEllipse : ParametricCurve return output; } - float_t y(float_t x) + float_t y(float_t x) const { return a * sqrt(1.0 - pow((x / b), 2.0)); } - float_t derivative(float_t x) + float_t derivative(float_t x) const { return (-a * x) / ((b * b) * sqrt(1.0 - pow((x / b), 2.0))); } @@ -349,7 +357,7 @@ struct ExplicitEllipse : ParametricCurve // Centered at (0,0), aligned with x template -struct AxisAlignedEllipse : ParametricCurve +struct AxisAlignedEllipse : impl::ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -366,13 +374,13 @@ struct AxisAlignedEllipse : ParametricCurve return output; } - float_t2 computePosition(float_t t) + 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) + float_t2 computeTangent(float_t t) const override { const float_t theta = start + (end - start) * t; const float_t dThetaDt = end - start; @@ -403,7 +411,7 @@ struct EllipticalArcInfo }; template -struct OffsettedBezier : ParametricCurve +struct OffsettedBezier : impl::ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -420,7 +428,7 @@ struct OffsettedBezier : ParametricCurve return output; } - float_t2 computePosition(float_t t) + 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)); @@ -428,7 +436,7 @@ struct OffsettedBezier : ParametricCurve } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) + float_t2 computeTangent(float_t t) const override { const float_t2 ddt = quadratic.derivative(t); const float_t2 d2dt2 = quadratic.secondDerivative(t); @@ -477,20 +485,21 @@ NBL_CONCEPT_END( #include template -struct QuadraticBezierFitter +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. - // TODO [Przemek]: add constraints on the `AddBezierFunc` type - template - static void adaptive(NBL_CONST_REF_ARG(ParametricCurve) curve, float_t min, float_t max, float_t targetMaxError, NBL_REF_ARG(AddBezierFunc) addBezierFunc, uint32_t maxDepth = 12) + static void adaptive(const impl::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) @@ -503,9 +512,7 @@ struct QuadraticBezierFitter adaptive_impl(curve, min, max, targetMaxError, addBezierFunc, maxDepth); } - // TODO [Przemek]: add constraints on the `AddBezierFunc` type - template - static void adaptive(NBL_CONST_REF_ARG(EllipticalArcInfo) ellipse, float_t targetMaxError, NBL_REF_ARG(AddBezierFunc) addBezierFunc, uint32_t maxDepth = 12) + static void adaptive(const EllipticalArcInfo& ellipse, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) { using namespace nbl::hlsl; @@ -542,8 +549,6 @@ struct QuadraticBezierFitter } } - // TODO [Przemek]: add constraints on the `AddBezierFunc` type - template static void adaptive(const OffsettedBezier& curve, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) { const float_t2 cusps = curve.findCusps(); @@ -568,7 +573,7 @@ struct QuadraticBezierFitter } } - +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 @@ -584,8 +589,7 @@ struct QuadraticBezierFitter } } - template - static void adaptive_impl(const ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t depth) + static void adaptive_impl(const impl::ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t depth) { if (min == max) return; @@ -671,52 +675,46 @@ namespace curves { // Mixes/Interpolation of two Parametric Curves t from 0 to 1 template -struct MixedParametricCurves : nbl::ext::curves::ParametricCurve +struct MixedParametricCurves : nbl::hlsl::ext::curves::impl::ParametricCurve { 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; - using ParametricCurve = nbl::ext::curves::ParametricCurve; + using ParametricCurve = nbl::hlsl::ext::curves::impl::ParametricCurve; - // TODO: rememeber these used to be pointers - ParametricCurve curve1; - ParametricCurve curve2; + const ParametricCurve* curve1; + const ParametricCurve* curve2; - static MixedParametricCurves create(NBL_CONST_REF_ARG(ParametricCurve) curve1, NBL_CONST_REF_ARG(ParametricCurve) curve2) - { - MixedParametricCurves output; - output.curve1 = curve1; - output.curve2 = curve2; - - return output; - } + MixedParametricCurves(const ParametricCurve* curve1, const ParametricCurve* curve2) + : curve1(curve1), curve2(curve2) + {} - float_t2 computePosition(float_t t) + float_t2 computePosition(float_t t) const override { - const float_t2 curve1Pos = curve1.computePosition(t); - const float_t2 curve2Pos = curve2.computePosition(t); + const float_t2 curve1Pos = curve1->computePosition(t); + const float_t2 curve2Pos = curve2->computePosition(t); return t * (curve2Pos - curve1Pos) + curve1Pos; } //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) + float_t2 computeTangent(float_t t) const override { - const float_t2 curve1Pos = curve1.computePosition(t); - const float_t2 curve2Pos = curve2.computePosition(t); - const float_t2 curve1Tan = curve1.computeTangent(t); - const float_t2 curve2Tan = curve2.computeTangent(t); + const float_t2 curve1Pos = curve1->computePosition(t); + const float_t2 curve2Pos = curve2->computePosition(t); + const float_t2 curve1Tan = curve1->computeTangent(t); + const float_t2 curve2Tan = curve2->computeTangent(t); return (1 - t) * curve1Tan - curve1Pos + (t)*curve2Tan + curve2Pos; } //! compute second order differential at t - float_t2 computeSecondOrderDifferential(float_t t) + float_t2 computeSecondOrderDifferential(float_t t) const override { - const float_t2 curve1Tan = curve1.computeTangent(t); - const float_t2 curve2Tan = curve2.computeTangent(t); - const float_t2 curve1SecondDiff = curve1.computeSecondOrderDifferential(t); - const float_t2 curve2SecondDiff = curve2.computeSecondOrderDifferential(t); + const float_t2 curve1Tan = curve1->computeTangent(t); + const float_t2 curve2Tan = curve2->computeTangent(t); + const float_t2 curve1SecondDiff = curve1->computeSecondOrderDifferential(t); + const float_t2 curve2SecondDiff = curve2->computeSecondOrderDifferential(t); return (1 - t) * curve1SecondDiff + 2.0 * (curve2Tan - curve1Tan) + t * curve2SecondDiff; } @@ -727,7 +725,7 @@ struct MixedParametricCurves : nbl::ext::curves::ParametricCurve return float_t(first.x * second.y - second.x * first.y); }; - float_t computeInflectionPoint(float_t errorThreshold) + float_t computeInflectionPoint(float_t errorThreshold) const override { const uint16_t MaxIterations = 32u; float_t low = 0.0; @@ -769,7 +767,7 @@ struct MixedParametricCurves : nbl::ext::curves::ParametricCurve // Mix between two parabolas from 0 to len template -struct MixedParabola : nbl::ext::hlsl::curves::ExplicitCurve +struct MixedParabola : nbl::hlsl::ext::curves::impl::ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -780,7 +778,7 @@ struct MixedParabola : nbl::ext::hlsl::curves::ExplicitCurve float_t a, b, c, d; - static MixedParabola create(NBL_CONST_REF_ARG(nbl::ext::hlsl::curves::Parabola) parabola1, NBL_CONST_REF_ARG(nbl::ext::hlsl::curves::Parabola) parabola2, float_t chordLen) + static MixedParabola create(NBL_CONST_REF_ARG(nbl::hlsl::ext::curves::Parabola) parabola1, NBL_CONST_REF_ARG(nbl::hlsl::ext::curves::Parabola) parabola2, float_t chordLen) { MixedParabola output; output.a = (parabola2.a - parabola1.a) / chordLen; @@ -800,17 +798,17 @@ struct MixedParabola : nbl::ext::hlsl::curves::ExplicitCurve return MixedParabola(parabola1, parabola2, abs(P2.x - P1.x)); } - float_t y(float_t x) + float_t y(float_t x) const { return (((a * x) + b) * x + c) * x + d; } - float_t derivative(float_t x) + float_t derivative(float_t x) const { return ((3.0 * a * x) + 2.0 * b) * x + c; } - float_t computeInflectionPoint(float_t errorThreshold) + float_t computeInflectionPoint(float_t errorThreshold) const override { return -b / (3.0 * a); } @@ -818,7 +816,7 @@ struct MixedParabola : nbl::ext::hlsl::curves::ExplicitCurve // Centered at (0, 0), P1 and P2 on x axis and P1.x = -P2.x template -struct ExplicitMixedCircle : nbl::hlsl::ext::curves::ParametricCurve +struct ExplicitMixedCircle : nbl::hlsl::ext::curves::impl::ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -866,7 +864,7 @@ struct ExplicitMixedCircle : nbl::hlsl::ext::curves::ParametricCurve //{ //} - float_t y(float_t x) + float_t y(float_t x) const { // https://herbie.uwplse.org/demo/5966aab2781d6c55c07105f5b900f1506cad301e.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html const float_t s1 = -1.0 * getSign(origin1Y); @@ -897,7 +895,7 @@ struct ExplicitMixedCircle : nbl::hlsl::ext::curves::ParametricCurve return ret; } - float_t derivative(float_t x) + float_t derivative(float_t x) const { // https://www.wolframalpha.com/input?i=derivative+%28x%2Fl%2B1%2F2%29*%28s2*sqrt%28r_2%5E2-x%5E2%29%2Borigin2Y-s1*sqrt%28r_1%5E2-x%5E2%29-origin1Y%29%2Bs1*sqrt%28r_1%5E2-x%5E2%29%2Borigin1Y // https ://herbie.uwplse.org/demo/d48eba9a858160f5f8e4283f7e6211d41215d354.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html @@ -911,7 +909,7 @@ struct ExplicitMixedCircle : nbl::hlsl::ext::curves::ParametricCurve return ret; } - float_t computeInflectionPoint(float_t errorThreshold) + float_t computeInflectionPoint(float_t errorThreshold) const override { // bisection search to find inflection point // by seeing the graph of second derivative over wide range of values we have deduced that the inflection point exists iff the secondDerivative has opposite signs at begin and end From cf93a6ef67a89dc43c1a17553a37d0e54d335330 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Wed, 26 Aug 2026 15:32:33 +0200 Subject: [PATCH 4/8] Added `Hatch.h` --- include/nbl/ext/Hatch/Hatch.h | 884 ++++++++++++++++++++++++++++++++++ 1 file changed, 884 insertions(+) create mode 100644 include/nbl/ext/Hatch/Hatch.h 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 From a8132a687f661ea5e0a66853f8a6a5c468bc98ad Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Thu, 27 Aug 2026 16:51:16 +0200 Subject: [PATCH 5/8] Refactored curves.hlsl --- .../nbl/builtin/hlsl/ext/curves/curves.hlsl | 349 +----------------- 1 file changed, 13 insertions(+), 336 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index 819e9703ab..8a550f8e4e 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -23,8 +23,6 @@ namespace ext namespace curves { -namespace impl -{ // Base class for all our curves template struct ParametricCurve @@ -101,17 +99,15 @@ struct ParametricCurve virtual float_t2 computeSecondOrderDifferential(float_t t) const { - return float_t2(num_traits::quiet_NaN(), num_traits::quiet_NaN()); + return float_t2(numeric_limits::quiet_NaN, numeric_limits::quiet_NaN); } virtual float_t computeInflectionPoint(float_t errorThreshold) const { - return num_traits::quiet_NaN(); + 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 @@ -135,7 +131,7 @@ inline nbl::hlsl::portable_vector_t2 computePosition(float_t x) const o template -struct Parabola : impl::ParametricCurve +struct Parabola : ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -176,7 +172,7 @@ struct Parabola : impl::ParametricCurve }; template -struct CubicCurve : impl::ParametricCurve +struct CubicCurve : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; using float_t4 = nbl::hlsl::portable_vector_t4; @@ -239,14 +235,14 @@ struct CubicCurve : impl::ParametricCurve return roots[0]; if (roots[1] <= 1.0 && roots[1] >= 0.0) return roots[1]; - return num_traits::quiet_NaN; + 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 : impl::ParametricCurve +struct CircularArc : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -328,7 +324,7 @@ struct CircularArc : impl::ParametricCurve // Centered at (0,0), aligned with x axis template -struct ExplicitEllipse : impl::ParametricCurve +struct ExplicitEllipse : ParametricCurve { DEFINE_EXPLICIT_CURVE_FUNCTIONS @@ -346,18 +342,18 @@ struct ExplicitEllipse : impl::ParametricCurve float_t y(float_t x) const { - return a * sqrt(1.0 - pow((x / b), 2.0)); + return a * sqrt(float_t(1.0) - pow(float_t(x / b), float_t(2.0))); } float_t derivative(float_t x) const { - return (-a * x) / ((b * b) * sqrt(1.0 - pow((x / b), 2.0))); + return (-a * x) / ((b * b) * sqrt(float_t(1.0) - pow(float_t(x / b), float_t(2.0)))); } }; // Centered at (0,0), aligned with x template -struct AxisAlignedEllipse : impl::ParametricCurve +struct AxisAlignedEllipse : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -411,7 +407,7 @@ struct EllipticalArcInfo }; template -struct OffsettedBezier : impl::ParametricCurve +struct OffsettedBezier : ParametricCurve { using float_t2 = nbl::hlsl::portable_vector_t2; @@ -457,33 +453,6 @@ struct OffsettedBezier : impl::ParametricCurve } }; -// TODO: should it stay in this header? -#include -// declare concept -#define NBL_CONCEPT_NAME BezierAdder -#define NBL_CONCEPT_TPLT_PRM_KINDS (typename)(typename)(int32_t)(int32_t) -#define NBL_CONCEPT_TPLT_PRM_NAMES (U)(T)(Dims)(Components) -// not the greatest syntax but works -#define NBL_CONCEPT_PARAM_0 (a,U) -#define NBL_CONCEPT_PARAM_1 (uv,vector) -#define NBL_CONCEPT_PARAM_2 (layer,uint16_t) -#define NBL_CONCEPT_PARAM_3 (data,vector) -// start concept -NBL_CONCEPT_BEGIN(4) -// need to be defined AFTER the cocnept begins -#define a NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 -#define uv NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 -#define layer NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_2 -#define data NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_3 -NBL_CONCEPT_END( - ((NBL_CONCEPT_REQ_EXPR)(a.template set(uv,layer,data))) -); -#undef data -#undef layer -#undef uv -#undef a -#include - template class QuadraticBezierFitter final { @@ -499,7 +468,7 @@ public: //! 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 impl::ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t maxDepth = 12) + 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) @@ -589,7 +558,7 @@ private: } } - static void adaptive_impl(const impl::ParametricCurve& curve, float_t min, float_t max, float_t targetMaxError, AddBezierFunc& addBezierFunc, uint32_t depth) + 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; @@ -658,305 +627,13 @@ private: } }; -} - } // namespace curves } // namespace ext } // namespace hlsl } // namespace nbl -#include "Shaders/globals.hlsl" - -namespace en -{ -namespace nabla2d -{ -namespace curves -{ -// Mixes/Interpolation of two Parametric Curves t from 0 to 1 -template -struct MixedParametricCurves : nbl::hlsl::ext::curves::impl::ParametricCurve -{ - 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; - - using ParametricCurve = nbl::hlsl::ext::curves::impl::ParametricCurve; - - const ParametricCurve* curve1; - const ParametricCurve* curve2; - - MixedParametricCurves(const ParametricCurve* curve1, const ParametricCurve* curve2) - : curve1(curve1), curve2(curve2) - {} - - float_t2 computePosition(float_t t) const override - { - const float_t2 curve1Pos = curve1->computePosition(t); - const float_t2 curve2Pos = curve2->computePosition(t); - return t * (curve2Pos - curve1Pos) + curve1Pos; - } - - //! compute unnormalized tangent vector at t - float_t2 computeTangent(float_t t) const override - { - const float_t2 curve1Pos = curve1->computePosition(t); - const float_t2 curve2Pos = curve2->computePosition(t); - const float_t2 curve1Tan = curve1->computeTangent(t); - const float_t2 curve2Tan = curve2->computeTangent(t); - return (1 - t) * curve1Tan - curve1Pos + (t)*curve2Tan + curve2Pos; - } - - //! compute second order differential at t - float_t2 computeSecondOrderDifferential(float_t t) const override - { - const float_t2 curve1Tan = curve1->computeTangent(t); - const float_t2 curve2Tan = curve2->computeTangent(t); - const float_t2 curve1SecondDiff = curve1->computeSecondOrderDifferential(t); - const float_t2 curve2SecondDiff = curve2->computeSecondOrderDifferential(t); - return (1 - t) * curve1SecondDiff + 2.0 * (curve2Tan - curve1Tan) + t * curve2SecondDiff; - } - - float_t signedCurvatureUnnormalized(float_t t) - { - const float_t2 first = computeTangent(t); - const float_t2 second = computeSecondOrderDifferential(t); - return float_t(first.x * second.y - second.x * first.y); - }; - - float_t computeInflectionPoint(float_t errorThreshold) const override - { - const uint16_t MaxIterations = 32u; - float_t low = 0.0; - float_t high = 1.0; - float_t valLow = signedCurvatureUnnormalized(low); - float_t valHigh = signedCurvatureUnnormalized(high); - if (getSign(valLow) != getSign(valHigh)) - { - if (valLow > valHigh) - hlsl::swap(low, high); - - float_t guess = 0.0; - for (uint16_t i = 0u; i < MaxIterations; ++i) - { - guess = (low + high) / 2.0; - float_t valGuess = signedCurvatureUnnormalized(guess); - if (abs(valGuess) < errorThreshold) - return guess; - - if (valGuess < 0.0) - low = guess; - else - high = guess; - } - - return guess; - } - else - { - return num_traits::quiet_NaN; - } - } - - static float_t getSign(float_t x) - { - return _static_cast((x > 0.0)) - _static_cast((x <= 0.0)); - } -}; - -// Mix between two parabolas from 0 to len -template -struct MixedParabola : nbl::hlsl::ext::curves::impl::ParametricCurve -{ - DEFINE_EXPLICIT_CURVE_FUNCTIONS - - 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; - - float_t a, b, c, d; - - static MixedParabola create(NBL_CONST_REF_ARG(nbl::hlsl::ext::curves::Parabola) parabola1, NBL_CONST_REF_ARG(nbl::hlsl::ext::curves::Parabola) parabola2, float_t chordLen) - { - MixedParabola output; - output.a = (parabola2.a - parabola1.a) / chordLen; - output.b = (parabola2.b - parabola1.b) / chordLen + parabola1.a; - output.c = (parabola2.c - parabola1.c) / chordLen + parabola1.b; - output.d = parabola1.c; - - return output; - } - - static MixedParabola fromFourPoints(NBL_CONST_REF_ARG(float_t2) P0, NBL_CONST_REF_ARG(float_t2) P1, NBL_CONST_REF_ARG(float_t2) P2, NBL_CONST_REF_ARG(float_t2) P3) - { - assert(P1.x == 0); - assert(P1.y == 0 && P2.y == 0); - nbl::hlsl::ext::curves::Parabola parabola1 = nbl::hlsl::ext::curves::Parabola::fromThreePoints(P0, P1, P2); - nbl::hlsl::ext::curves::Parabola parabola2 = nbl::hlsl::ext::curves::Parabola::fromThreePoints(P1, P2, P3); - return MixedParabola(parabola1, parabola2, abs(P2.x - P1.x)); - } - - float_t y(float_t x) const - { - return (((a * x) + b) * x + c) * x + d; - } - - float_t derivative(float_t x) const - { - return ((3.0 * a * x) + 2.0 * b) * x + c; - } - - float_t computeInflectionPoint(float_t errorThreshold) const override - { - return -b / (3.0 * a); - } -}; - -// Centered at (0, 0), P1 and P2 on x axis and P1.x = -P2.x -template -struct ExplicitMixedCircle : nbl::hlsl::ext::curves::impl::ParametricCurve -{ - DEFINE_EXPLICIT_CURVE_FUNCTIONS - - 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 ExplicitCircle - { - float_t2 origin; - float_t radius; - - static ExplicitCircle create(NBL_CONST_REF_ARG(float_t2) origin, float_t radius) - { - ExplicitCircle output; - output.origin = origin; - output.radius = radius; - - return output; - } - - static ExplicitCircle fromThreePoints(float_t2 P0, float_t2 P1, float_t2 P2) - { - const float_t2 Mid0 = (P0 + P1) / 2.0; - const float_t2 Normal0 = float_t2(P1.y - P0.y, P0.x - P1.x); - - const float_t2 Mid1 = (P1 + P2) / 2.0; - const float_t2 Normal1 = float_t2(P2.y - P1.y, P1.x - P2.x); - - const float_t2 origin = shapes::util::LineLineIntersection(Mid0, Normal0, Mid1, Normal1); - const float_t radius = glm::length(P0 - origin); - return ExplicitCircle(origin, radius); - } - }; - - float_t origin1Y; - float_t origin2Y; - float_t radius1; - float_t radius2; - float_t chordLen; - - // TODO: I haven't found an implementation - //static ExplicitMixedCircle fromFourPoints(NBL_CONST_REF_ARG(float_t2) P0, NBL_CONST_REF_ARG(float_t2) P1, NBL_CONST_REF_ARG(float_t2) P2, NBL_CONST_REF_ARG(float_t2) P3) - //{ - //} - - float_t y(float_t x) const - { - // https://herbie.uwplse.org/demo/5966aab2781d6c55c07105f5b900f1506cad301e.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html - const float_t s1 = -1.0 * getSign(origin1Y); - const float_t s2 = -1.0 * getSign(origin2Y); - - const float_t t1 = s1 * sqrt(radius1 * radius1 - x * x) + origin1Y; - const float_t t2 = s2 * sqrt(radius2 * radius2 - x * x) + origin2Y; - const float_t ret = (x / chordLen + 0.5) * (t2 - t1) + t1; - return ret; - } - - float_t secondDerivative(float_t x) - { - // https://www.wolframalpha.com/input?i=second+derivative+of+%28x%2Fl%2B1%2F2%29*%28s2*sqrt%28r_2%5E2-x%5E2%29%2Bo2-s1*sqrt%28r_1%5E2-x%5E2%29-o1%29%2Bs1*sqrt%28r_1%5E2-x%5E2%29%2Bo1 - const float_t s1 = -1.0 * getSign(origin1Y); - const float_t s2 = -1.0 * getSign(origin2Y); - - const float_t t1 = sqrt(radius1 * radius1 - x * x); - const float_t t2 = sqrt(radius2 * radius2 - x * x); - - const float_t u1 = (s1 * x) / t1; - const float_t u2 = (s2 * x) / t2; - - const float_t q1 = (-1.0 * (x * x) / pow(radius1 * radius1 - x * x, 1.5)) - (1.0 / t1); - const float_t q2 = (-1.0 * (x * x) / pow(radius2 * radius2 - x * x, 1.5)) - (1.0 / t2); - - const float_t ret = ((2.0 * (u1 - u2)) / chordLen) + (x / chordLen + 0.5) * (s2 * q2 - s1 * q1) + s1 * q1; - return ret; - } - - float_t derivative(float_t x) const - { - // https://www.wolframalpha.com/input?i=derivative+%28x%2Fl%2B1%2F2%29*%28s2*sqrt%28r_2%5E2-x%5E2%29%2Borigin2Y-s1*sqrt%28r_1%5E2-x%5E2%29-origin1Y%29%2Bs1*sqrt%28r_1%5E2-x%5E2%29%2Borigin1Y - // https ://herbie.uwplse.org/demo/d48eba9a858160f5f8e4283f7e6211d41215d354.10eb16397304ba9e003d59ffd34803200ce7cfa6/graph.html - const float_t s1 = -1.0 * getSign(origin1Y); - const float_t s2 = -1.0 * getSign(origin2Y); - - const float_t t0 = sqrt(radius1 * radius1 - x * x); - const float_t t1 = (s1 * x) / t0; - const float_t t2 = sqrt(radius2 * radius2 - x * x); - const float_t ret = ((x / chordLen + 0.5) * (t1 - ((x * s2) / t2)) + (((origin2Y - origin1Y) + (s2 * t2) - (s1 * t0)) / chordLen)) - t1; - return ret; - } - - float_t computeInflectionPoint(float_t errorThreshold) const override - { - // bisection search to find inflection point - // by seeing the graph of second derivative over wide range of values we have deduced that the inflection point exists iff the secondDerivative has opposite signs at begin and end - const uint16_t MaxIterations = 64u; - float_t low = -chordLen / 2.0; - float_t high = chordLen / 2.0; - float_t valLow = secondDerivative(low + errorThreshold / 2.0); - float_t valHigh = secondDerivative(high - errorThreshold / 2.0); - if (getSign(valLow) != getSign(valHigh)) - { - if (valLow > valHigh) - hlsl::swap(low, high); - - float_t guess = 0.0; - for (uint16_t i = 0u; i < MaxIterations; ++i) - { - guess = (low + high) / 2.0; - float_t valGuess = secondDerivative(guess); - if (abs(valGuess) < errorThreshold) - return guess; - - if (valGuess < 0.0) - low = guess; - else - high = guess; - } - - return guess; - } - else - { - return num_traits::quiet_NaN; - } - } - - static float_t getSign(float_t x) - { - return _static_cast((x > 0.0)) - _static_cast((x <= 0.0)); - } -}; - #undef DEFINE_EXPLICIT_CURVE_FUNCTIONS -} // namespace curves -} // namespace nabla2d -} // namespace en - #endif #endif \ No newline at end of file From da1099e60a5017e342f6f31cf9c4cfdf0ff7a544 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Thu, 27 Aug 2026 16:55:36 +0200 Subject: [PATCH 6/8] Additional refactor --- include/nbl/builtin/hlsl/ext/curves/curves.hlsl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index 8a550f8e4e..76d7e49724 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -342,12 +342,12 @@ struct ExplicitEllipse : ParametricCurve float_t y(float_t x) const { - return a * sqrt(float_t(1.0) - pow(float_t(x / b), float_t(2.0))); + 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(float_t(1.0) - pow(float_t(x / b), float_t(2.0)))); + return (-a * x) / ((b * b) * sqrt(static_cast(1.0) - pow((x / b), static_cast(2.0)))); } }; From 2cade9602da7b7be9e4ebcf8ca06fd50ac6263f5 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Thu, 27 Aug 2026 19:53:00 +0200 Subject: [PATCH 7/8] Restored the `ExplicitCurve` struct --- .../nbl/builtin/hlsl/ext/curves/curves.hlsl | 51 ++++++++++--------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index 76d7e49724..975d9d229b 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -111,30 +111,35 @@ struct ParametricCurve // 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 -#define DEFINE_EXPLICIT_CURVE_FUNCTIONS \ -float_t differentialArcLen(float_t x) const override\ -{ \ - float_t deriv = derivative(x); \ - return sqrt(1.0 + deriv * deriv); \ -} \ - \ -nbl::hlsl::portable_vector_t2 computeTangent(float_t x) const override\ -{ \ - const float_t deriv = derivative(x); \ - nbl::hlsl::portable_vector_t2 v = nbl::hlsl::portable_vector_t2(1.0, deriv); \ - if (nbl::hlsl::isinf(deriv)) \ - v = nbl::hlsl::portable_vector_t2(0.0, 1.0); \ - return v; \ -} \ - \ -inline nbl::hlsl::portable_vector_t2 computePosition(float_t x) const override { return nbl::hlsl::portable_vector_t2(x, y(x)); } +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 : ParametricCurve +struct Parabola : ExplicitCurve { - DEFINE_EXPLICIT_CURVE_FUNCTIONS - using float_t2 = nbl::hlsl::portable_vector_t2; float_t a, b, c; @@ -324,10 +329,8 @@ struct CircularArc : ParametricCurve // Centered at (0,0), aligned with x axis template -struct ExplicitEllipse : ParametricCurve +struct ExplicitEllipse : ExplicitCurve { - DEFINE_EXPLICIT_CURVE_FUNCTIONS - using float_t2 = nbl::hlsl::portable_vector_t2; float_t a, b; @@ -632,8 +635,6 @@ private: } // namespace hlsl } // namespace nbl -#undef DEFINE_EXPLICIT_CURVE_FUNCTIONS - #endif #endif \ No newline at end of file From f61a2b0ffb92cceeae7fa029dfe2730d604fc779 Mon Sep 17 00:00:00 2001 From: Przemog1 Date: Thu, 27 Aug 2026 20:32:07 +0200 Subject: [PATCH 8/8] Fixed bugs --- include/nbl/builtin/hlsl/ext/curves/curves.hlsl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl index 975d9d229b..06e8471af4 100644 --- a/include/nbl/builtin/hlsl/ext/curves/curves.hlsl +++ b/include/nbl/builtin/hlsl/ext/curves/curves.hlsl @@ -63,7 +63,7 @@ struct ParametricCurve } //! 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) + 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; @@ -92,7 +92,7 @@ struct ParametricCurve } //! compute inverse arc len - float_t inverseArcLen(float_t targetLen, float_t min, float_t max, const float_t cdfAccuracyThreshold = 1e-4) + 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); } @@ -397,7 +397,7 @@ struct EllipticalArcInfo float_t2 angleBounds; // [0, 2Pi) float_t eccentricity; // (0, 1] - inline bool isValid() + inline bool isValid() const { if (eccentricity > 1.0 || eccentricity <= 0.0) return false; @@ -444,7 +444,7 @@ struct OffsettedBezier : ParametricCurve } //! if offset is more than minimum radius of curvature then we get an unwanted gouging/cusp - float_t2 findCusps() + 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); @@ -516,7 +516,7 @@ public: if (ellipse.angleBounds.x != ellipse.angleBounds.y) { - AxisAlignedEllipse aaEllipse(lenghtMajor, lenghtMinor, 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); } }