| /* |
| * Copyright 2020 Google Inc. |
| * |
| * Use of this source code is governed by a BSD-style license that can be |
| * found in the LICENSE file. |
| */ |
| |
| #include "rive/math/math_types.hpp" |
| #include "rive/math/bezier_utils.hpp" |
| #include "rive/math/simd.hpp" |
| #include "common/rand.hpp" |
| #include <catch.hpp> |
| |
| // glsl cross-compiling requires the clang/gcc vector extension. |
| #ifndef _MSC_VER |
| #include "../renderer/cpp.glsl" |
| namespace glsl_cross |
| { |
| #include "generated/shaders/bezier_utils.minified.glsl" |
| } |
| #endif |
| |
| namespace rive |
| { |
| namespace math |
| { |
| constexpr static float kEpsilon = 1.f / (1 << 12); |
| |
| static bool fuzzy_equal(float a, float b, float tolerance = kEpsilon) |
| { |
| assert(tolerance >= 0); |
| return fabsf(a - b) <= tolerance; |
| } |
| |
| static float frand() { return rand() / static_cast<float>(RAND_MAX); } |
| |
| TEST_CASE("chop_cubic_at", "[bezier_utils]") |
| { |
| // Make sure src and dst can be the same pointer. |
| { |
| Vec2D pts[7]; |
| for (int i = 0; i < 7; ++i) |
| { |
| pts[i] = {(float)i, (float)i}; |
| } |
| math::chop_cubic_at(pts, pts, .5f); |
| for (int i = 0; i < 7; ++i) |
| { |
| CHECK(pts[i].x == pts[i].y); |
| CHECK(pts[i].x == i * .5f); |
| } |
| } |
| |
| static const float chopTs[] = { |
| 0, 3 / 83.f, 3 / 79.f, 3 / 73.f, 3 / 71.f, 3 / 67.f, |
| 3 / 61.f, 3 / 59.f, 3 / 53.f, 3 / 47.f, 3 / 43.f, 3 / 41.f, |
| 3 / 37.f, 3 / 31.f, 3 / 29.f, 3 / 23.f, 3 / 19.f, 3 / 17.f, |
| 3 / 13.f, 3 / 11.f, 3 / 7.f, 3 / 5.f, 1, |
| }; |
| float ones[] = {1, 1, 1, 1, 1}; |
| |
| // Ensure an odd number of T values so we exercise the single chop code at |
| // the end of SkChopCubicAt form multiple T. |
| static_assert(std::size(chopTs) % 2 == 1); |
| static_assert(std::size(ones) % 2 == 1); |
| |
| for (int iterIdx = 0; iterIdx < 5; ++iterIdx) |
| { |
| Vec2D pts[4] = {{frand(), frand()}, |
| {frand(), frand()}, |
| {frand(), frand()}, |
| {frand(), frand()}}; |
| |
| Vec2D allChops[4 + std::size(chopTs) * 3]; |
| math::chop_cubic_at(pts, allChops, chopTs, std::size(chopTs)); |
| int i = 3; |
| for (float chopT : chopTs) |
| { |
| // Ensure we chop at approximately the correct points when we chop |
| // an entire list. |
| Vec2D expectedPt = math::eval_cubic_at(pts, chopT); |
| CHECK(fuzzy_equal(allChops[i].x, expectedPt.x)); |
| CHECK(fuzzy_equal(allChops[i].y, expectedPt.y)); |
| if (chopT == 0) |
| { |
| CHECK(allChops[i] == pts[0]); |
| } |
| if (chopT == 1) |
| { |
| CHECK(allChops[i] == pts[3]); |
| } |
| i += 3; |
| |
| // Ensure the middle is exactly degenerate when we chop at two equal |
| // points. |
| Vec2D localChops[10]; |
| math::chop_cubic_at(pts, localChops, chopT, chopT); |
| CHECK(localChops[3] == localChops[4]); |
| CHECK(localChops[3] == localChops[5]); |
| CHECK(localChops[3] == localChops[6]); |
| if (chopT == 0) |
| { |
| // Also ensure the first curve is exactly p0 when we chop at |
| // T=0. |
| CHECK(localChops[0] == pts[0]); |
| CHECK(localChops[1] == pts[0]); |
| CHECK(localChops[2] == pts[0]); |
| CHECK(localChops[3] == pts[0]); |
| } |
| if (chopT == 1) |
| { |
| // Also ensure the last curve is exactly p3 when we chop at T=1. |
| CHECK(localChops[6] == pts[3]); |
| CHECK(localChops[7] == pts[3]); |
| CHECK(localChops[8] == pts[3]); |
| CHECK(localChops[9] == pts[3]); |
| } |
| } |
| |
| // Now test what happens when SkChopCubicAt does 0/0 and gets NaN |
| // values. |
| Vec2D oneChops[4 + std::size(ones) * 3]; |
| math::chop_cubic_at(pts, oneChops, ones, std::size(ones)); |
| CHECK(oneChops[0] == pts[0]); |
| CHECK(oneChops[1] == pts[1]); |
| CHECK(oneChops[2] == pts[2]); |
| for (size_t index = 3; index < std::size(oneChops); ++index) |
| { |
| CHECK(oneChops[index] == pts[3]); |
| } |
| } |
| } |
| |
| TEST_CASE("chop_cubic_at(tValues=null)", "[bezier_utils]") |
| { |
| constexpr int kMaxNumChops = 20; |
| float chopT[kMaxNumChops]; |
| Vec2D chopTChops[(kMaxNumChops + 1) * 3 + 1]; |
| Vec2D nullTChops[(kMaxNumChops + 1) * 3 + 1]; |
| for (int numChops = 1; numChops <= kMaxNumChops; ++numChops) |
| { |
| Vec2D pts[4] = {{frand(), frand()}, |
| {frand(), frand()}, |
| {frand(), frand()}, |
| {frand(), frand()}}; |
| float stepT = 1.f / (numChops + 1); |
| float t = 0; |
| for (size_t i = 0; i < numChops; ++i) |
| { |
| chopT[i] = t += stepT; |
| } |
| math::chop_cubic_at(pts, chopTChops, chopT, numChops); |
| math::chop_cubic_at(pts, nullTChops, nullptr, numChops); |
| size_t numChoptPts = numChops * 3 + 1; |
| for (size_t i = 0; i <= numChoptPts; ++i) |
| { |
| CHECK(nullTChops[i].x == Approx(chopTChops[i].x)); |
| CHECK(nullTChops[i].y == Approx(chopTChops[i].y)); |
| } |
| } |
| } |
| |
| static float measure_non_inflect_cubic_rotation(const Vec2D pts[4]) |
| { |
| Vec2D a = pts[1] - pts[0]; |
| Vec2D b = pts[2] - pts[1]; |
| Vec2D c = pts[3] - pts[2]; |
| float lengthA = a.length(); |
| float lengthB = b.length(); |
| float lengthC = c.length(); |
| float lengthMax = std::max(std::max(lengthA, lengthB), lengthC); |
| float zeroThreshold = std::max(lengthMax, 1.f) * 1e-4f; |
| if (lengthA <= zeroThreshold) |
| { |
| return math::measure_angle_between_vectors(b, c); |
| } |
| if (lengthB <= zeroThreshold) |
| { |
| return math::measure_angle_between_vectors(a, c); |
| } |
| if (lengthC <= zeroThreshold) |
| { |
| return math::measure_angle_between_vectors(a, b); |
| } |
| // Postulate: When no points are colocated and there are no inflection |
| // points in T=0..1, the rotation is: 360 degrees, minus the angle |
| // [p0,p1,p2], minus the angle [p1,p2,p3]. |
| return 2 * math::PI - math::measure_angle_between_vectors(a, -b) - |
| math::measure_angle_between_vectors(b, -c); |
| } |
| |
| TEST_CASE("measure_non_inflect_cubic_rotation", "[bezier_utils]") |
| { |
| static Vec2D kFlatCubic[4] = {{0, 0}, {0, 1}, {0, 2}, {0, 3}}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(kFlatCubic), 0)); |
| |
| static Vec2D kFlatCubic180_1[4] = {{0, 0}, {1, 0}, {3, 0}, {2, 0}}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(kFlatCubic180_1), |
| math::PI)); |
| |
| static Vec2D kFlatCubic180_2[4] = {{0, 1}, {0, 0}, {0, 2}, {0, 3}}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(kFlatCubic180_2), |
| math::PI)); |
| |
| static Vec2D kFlatCubic360[4] = {{0, 1}, {0, 0}, {0, 3}, {0, 2}}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(kFlatCubic360), |
| 2 * math::PI)); |
| |
| static Vec2D kSquare180[4] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; |
| CHECK( |
| fuzzy_equal(measure_non_inflect_cubic_rotation(kSquare180), math::PI)); |
| |
| auto checkQuadRotation = [](const Vec2D pts[3], float expectedRotation) { |
| #if 0 |
| float r = SkMeasureQuadRotation(pts); |
| CHECK(fuzzy_equal(r, expectedRotation)); |
| #endif |
| Vec2D cubic1[4] = {pts[0], pts[0], pts[1], pts[2]}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(cubic1), |
| expectedRotation)); |
| |
| Vec2D cubic2[4] = {pts[0], pts[1], pts[1], pts[2]}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(cubic2), |
| expectedRotation)); |
| |
| Vec2D cubic3[4] = {pts[0], pts[1], pts[2], pts[2]}; |
| CHECK(fuzzy_equal(measure_non_inflect_cubic_rotation(cubic3), |
| expectedRotation)); |
| }; |
| |
| static Vec2D kFlatQuad[4] = {{0, 0}, {0, 1}, {0, 2}}; |
| checkQuadRotation(kFlatQuad, 0); |
| |
| static Vec2D kFlatQuad180_1[4] = {{1, 0}, {0, 0}, {2, 0}}; |
| checkQuadRotation(kFlatQuad180_1, math::PI); |
| |
| static Vec2D kFlatQuad180_2[4] = {{0, 0}, {0, 2}, {0, 1}}; |
| checkQuadRotation(kFlatQuad180_2, math::PI); |
| |
| static Vec2D kTri120[3] = {{0, 0}, {.5f, std::sqrt(3.f) / 2}, {1, 0}}; |
| checkQuadRotation(kTri120, 2 * math::PI / 3); |
| } |
| |
| static bool is_linear(Vec2D p0, Vec2D p1, Vec2D p2) |
| { |
| return fuzzy_equal(Vec2D::cross(p0 - p1, p2 - p1), 0); |
| } |
| |
| static bool is_linear(const Vec2D p[4]) |
| { |
| return is_linear(p[0], p[1], p[2]) && is_linear(p[0], p[2], p[3]) && |
| is_linear(p[1], p[2], p[3]); |
| } |
| |
| static int valid_unit_divide(float numer, float denom, float* ratio) |
| { |
| assert(ratio); |
| |
| if (numer < 0) |
| { |
| numer = -numer; |
| denom = -denom; |
| } |
| |
| if (denom == 0 || numer == 0 || numer >= denom) |
| { |
| return 0; |
| } |
| |
| float r = numer / denom; |
| if (std::isnan(r)) |
| { |
| return 0; |
| } |
| assert(r >= 0 && r < 1.f); |
| if (r == 0) |
| { // catch underflow if numer <<<< denom |
| return 0; |
| } |
| *ratio = r; |
| return 1; |
| } |
| |
| // Just returns its argument, but makes it easy to set a break-point to know |
| // when SkFindUnitQuadRoots is going to return 0 (an error). |
| static int return_check_zero(int value) |
| { |
| if (value == 0) |
| { |
| return 0; |
| } |
| return value; |
| } |
| |
| static int SkFindUnitQuadRoots(float A, float B, float C, float roots[2]) |
| { |
| assert(roots); |
| |
| if (A == 0) |
| { |
| return return_check_zero(valid_unit_divide(-C, B, roots)); |
| } |
| |
| float* r = roots; |
| |
| // use doubles so we don't overflow temporarily trying to compute R |
| double dr = (double)B * B - 4 * (double)A * C; |
| if (dr < 0) |
| { |
| return return_check_zero(0); |
| } |
| dr = sqrt(dr); |
| float R = static_cast<float>(dr); |
| if (std::isnan(R) || std::isinf(R)) |
| { |
| return return_check_zero(0); |
| } |
| |
| float Q = (B < 0) ? -(B - R) / 2 : -(B + R) / 2; |
| r += valid_unit_divide(Q, A, r); |
| r += valid_unit_divide(C, Q, r); |
| if (r - roots == 2) |
| { |
| if (roots[0] > roots[1]) |
| { |
| using std::swap; |
| swap(roots[0], roots[1]); |
| } |
| else if (roots[0] == roots[1]) |
| { // nearly-equal? |
| r -= 1; // skip the double root |
| } |
| } |
| return return_check_zero((int)(r - roots)); |
| } |
| |
| /** http://www.faculty.idc.ac.il/arik/quality/appendixA.html |
| |
| Inflection means that curvature is zero. |
| Curvature is [F' x F''] / [F'^3] |
| So we solve F'x X F''y - F'y X F''y == 0 |
| After some canceling of the cubic term, we get |
| A = b - a |
| B = c - 2b + a |
| C = d - 3c + 3b - a |
| (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 |
| */ |
| static int SkFindCubicInflections(const Vec2D src[4], float tValues[2]) |
| { |
| float Ax = src[1].x - src[0].x; |
| float Ay = src[1].y - src[0].y; |
| float Bx = src[2].x - 2 * src[1].x + src[0].x; |
| float By = src[2].y - 2 * src[1].y + src[0].y; |
| float Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x; |
| float Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y; |
| |
| return SkFindUnitQuadRoots(Bx * Cy - By * Cx, |
| Ax * Cy - Ay * Cx, |
| Ax * By - Ay * Bx, |
| tValues); |
| } |
| |
| static void check_cubic_convex_180(const Vec2D p[4]) |
| { |
| bool areCusps = false; |
| float inflectT[2], convex180T[2]; |
| if (int inflectN = SkFindCubicInflections(p, inflectT)) |
| { |
| // The curve has inflections. FindCubicConvex180Chops should return the |
| // inflection points. |
| int convex180N = |
| math::find_cubic_convex_180_chops(p, convex180T, &areCusps); |
| REQUIRE(inflectN == convex180N); |
| if (!areCusps) |
| { |
| REQUIRE((inflectN == 1 || |
| fabsf(inflectT[0] - inflectT[1]) >= kEpsilon)); |
| } |
| for (int i = 0; i < convex180N; ++i) |
| { |
| REQUIRE(fuzzy_equal(inflectT[i], convex180T[i])); |
| } |
| } |
| else |
| { |
| float totalRotation = measure_non_inflect_cubic_rotation(p); |
| int convex180N = |
| math::find_cubic_convex_180_chops(p, convex180T, &areCusps); |
| Vec2D chops[10]; |
| math::chop_cubic_at(p, chops, convex180T, convex180N); |
| float radsSum = 0; |
| if (areCusps) |
| { |
| if (convex180N == 1) |
| { |
| // The cusp itself accounts for 180 degrees of rotation. |
| radsSum = math::PI; |
| // Rechop, this time straddling the cusp, to avoid fp32 |
| // precision issues from crossover. |
| Vec2D straddleChops[10]; |
| math::chop_cubic_at(p, |
| straddleChops, |
| std::max(convex180T[0] - math::PI, 0.f), |
| std::min(convex180T[0] + math::PI, 1.f)); |
| memcpy(chops + 1, straddleChops + 1, 2 * sizeof(Vec2D)); |
| memcpy(chops + 4, straddleChops + 7, 2 * sizeof(Vec2D)); |
| } |
| else if (convex180N == 2) |
| { |
| // The only way a cubic can have two cusps is if it's flat with |
| // two 180 degree turnarounds. |
| radsSum = math::PI * 2; |
| chops[1] = chops[0]; |
| chops[2] = chops[4] = chops[3]; |
| chops[5] = chops[7] = chops[6]; |
| chops[8] = chops[9]; |
| } |
| } |
| for (int i = 0; i <= convex180N; ++i) |
| { |
| float rads = measure_non_inflect_cubic_rotation(chops + i * 3); |
| assert(rads < math::PI + kEpsilon); |
| radsSum += rads; |
| } |
| if (totalRotation < math::PI - kEpsilon) |
| { |
| // The curve should never chop if rotation is <180 degrees. |
| REQUIRE(convex180N == 0); |
| } |
| else if (!is_linear(p)) |
| { |
| REQUIRE(fuzzy_equal(radsSum, totalRotation)); |
| if (totalRotation > math::PI + kEpsilon) |
| { |
| REQUIRE(convex180N == 1); |
| // This works because cusps take the "inflection" path above, so |
| // we don't get non-lilnear curves that lose rotation when |
| // chopped. |
| REQUIRE(fuzzy_equal(measure_non_inflect_cubic_rotation(chops), |
| math::PI)); |
| REQUIRE( |
| fuzzy_equal(measure_non_inflect_cubic_rotation(chops + 3), |
| totalRotation - math::PI)); |
| } |
| REQUIRE(!areCusps); |
| } |
| else |
| { |
| REQUIRE(areCusps); |
| } |
| } |
| } |
| |
| TEST_CASE("find_cubic_convex_180_chops", "[bezier_utils]") |
| { |
| // Test all combinations of corners from the square [0,0,1,1]. This covers |
| // every cubic type as well as a wide variety of special cases for cusps, |
| // lines, loops, and inflections. |
| for (int i = 0; i < (1 << 8); ++i) |
| { |
| Vec2D p[4] = {Vec2D((i >> 0) & 1, (i >> 1) & 1), |
| Vec2D((i >> 2) & 1, (i >> 3) & 1), |
| Vec2D((i >> 4) & 1, (i >> 5) & 1), |
| Vec2D((i >> 6) & 1, (i >> 7) & 1)}; |
| check_cubic_convex_180(p); |
| } |
| |
| { |
| // This cubic has a convex-180 chop at T=1-"epsilon" |
| static const uint32_t hexPts[] = {0x3ee0ac74, |
| 0x3f1e061a, |
| 0x3e0fc408, |
| 0x3f457230, |
| 0x3f42ac7c, |
| 0x3f70d76c, |
| 0x3f4e6520, |
| 0x3f6acafa}; |
| Vec2D p[4]; |
| memcpy(p, hexPts, sizeof(p)); |
| check_cubic_convex_180(p); |
| } |
| |
| // Now test an exact quadratic. |
| Vec2D quad[4] = {{0, 0}, {2, 2}, {4, 2}, {6, 0}}; |
| float T[2]; |
| bool areCusps; |
| CHECK(math::find_cubic_convex_180_chops(quad, T, &areCusps) == 0); |
| |
| // Now test that cusps and near-cusps get flagged as cusps. |
| Vec2D cusp[4] = {{0, 0}, {1, 1}, {1, 0}, {0, 1}}; |
| CHECK(math::find_cubic_convex_180_chops(cusp, T, &areCusps) == 1); |
| CHECK(areCusps == true); |
| |
| // Find the height of the right side of "cusp" at which the distance between |
| // its inflection points is epsilon (in parametric space). |
| constexpr double epsilon = 1.0 / (1 << 11); |
| constexpr double epsilonPow2 = epsilon * epsilon; |
| double h = (1 - epsilonPow2) / (3 * epsilonPow2 + 1); |
| double dy = (1 - h) / 2; |
| cusp[1].y = (float)(1 - dy); |
| cusp[2].y = (float)(0 + dy); |
| CHECK(SkFindCubicInflections(cusp, T) == 2); |
| CHECK(fuzzy_equal(T[1] - T[0], (float)epsilon, (float)epsilonPow2)); |
| |
| // Ensure two inflection points barely more than epsilon apart do not get |
| // flagged as cusps. |
| cusp[1].y = (float)(1 - 4 * dy); |
| cusp[2].y = (float)(0 + 4 * dy); |
| CHECK(math::find_cubic_convex_180_chops(cusp, T, &areCusps) == 2); |
| CHECK(areCusps == false); |
| |
| // Ensure two inflection points barely less than epsilon apart do get |
| // flagged as cusps. |
| cusp[1].y = (float)(1 - .9 * dy); |
| cusp[2].y = (float)(0 + .9 * dy); |
| CHECK(math::find_cubic_convex_180_chops(cusp, T, &areCusps) == 1); |
| CHECK(areCusps == true); |
| |
| // Ensure fast floatingpoint is not enbled. |
| float2 p[] = {{460, 1060}, |
| {774, 526}, |
| {60, 660}, |
| {460, 460}, |
| {667, 460}, |
| {1060, 460}, |
| {686, 460}, |
| {686, 660}, |
| {1042, 1020}}; |
| float2 c0 = simd::mix(p[3], p[4], float2(2 / 3.f)); |
| float2 c1 = simd::mix(p[5], p[4], float2(2 / 3.f)); |
| Vec2D cubic[] = {math::bit_cast<Vec2D>(p[3]), |
| math::bit_cast<Vec2D>(c0), |
| math::bit_cast<Vec2D>(c1), |
| math::bit_cast<Vec2D>(p[5])}; |
| CHECK(math::find_cubic_convex_180_chops(cubic, T, &areCusps) == 0); |
| } |
| |
| // Check that flat lines with ordered points never get detected as cusps. |
| TEST_CASE("find_cubic_convex_180_chops_lines", "[bezier_utils]") |
| { |
| Vec2D p0 = {123, 200}, p3 = {223, 432}; |
| for (float t0 = 1e-3f; t0 < 1; t0 += .12f) |
| { |
| for (float t1 = t0 + .097f; t1 < 1; t1 += .097f) |
| { |
| Vec2D line[4] = {p0, lerp(p0, p3, t0), lerp(p0, p3, t1), p3}; |
| float _[2]; |
| bool areCusps = true; |
| math::find_cubic_convex_180_chops(line, _, &areCusps); |
| CHECK_FALSE(areCusps); |
| } |
| } |
| } |
| |
| Vec2D eval_cubic(const Vec2D* p, float t) |
| { |
| return math::pow3(1 - t) * p[0] + 3 * math::pow2(1 - t) * t * p[1] + |
| 3 * (1 - t) * math::pow2(t) * p[2] + math::pow3(t) * p[3]; |
| } |
| |
| constexpr static Vec2D TEST_CUBICS[][4] = { |
| {{199, 1225}, {197, 943}, {349, 607}, {549, 427}}, |
| {{549, 427}, {349, 607}, {197, 943}, {199, 1225}}, |
| {{460, 1060}, {403, -320}, {60, 660}, {1181, 634}}, |
| {{1181, 634}, {60, 660}, {403, -320}, {460, 1060}}, |
| {{0, 0}, {0, 0}, {0, 0}, {0, 0}}, |
| {{0, 0}, {0, 0}, {0, 0}, {100, 100}}, |
| {{0, 0}, {0, 0}, {100, 100}, {100, 100}}, |
| {{0, 0}, {100, 100}, {100, 100}, {0, 0}}, |
| {{-100, -100}, {100, 100}, {100, -100}, {-100, 100}}, // Cusp |
| {{0, 0}, {0, 0}, {100, 100}, {100, 100}}, // Line |
| {{0, 0}, {-100, -100}, {200, 200}, {100, 100}}, // Line w/ 2 cusps |
| {{0, 0}, |
| {50 * 2.f / 3.f, 100 * 2.f / 3.f}, |
| {100 - 50 * 2.f / 3.f, 100 * 2.f / 3.f}, |
| {100, 0}}, // Quadratic |
| // The remaining cubics had some sort of issue during development. They're |
| // in the list to make sure they don't regress. |
| {{0, 0}, |
| {50 * 2.f / 3.f, 100 * 2.f / 3.f}, |
| {100 - 50 * 2.f / 3.f, 100 * 2.f / 3.f}, |
| {100, 100}}, |
| {{100, 0}, {0, 0}, {0, 0}, {0, 0}}, |
| }; |
| |
| static void check_eval_cubic(const EvalCubic& evalCubic, |
| const Vec2D* cubic, |
| float t0, |
| float t1) |
| { |
| float4 pp = evalCubic(float4{t0, t0, t1, t1}); |
| Vec2D p0ref = eval_cubic(cubic, t0); |
| Vec2D p1ref = eval_cubic(cubic, t1); |
| CHECK(pp.x == Approx(p0ref.x).margin(1e-3f)); |
| CHECK(pp.y == Approx(p0ref.y).margin(1e-3f)); |
| CHECK(pp.z == Approx(p1ref.x).margin(1e-3f)); |
| CHECK(pp.w == Approx(p1ref.y).margin(1e-3f)); |
| float2 p0 = evalCubic(t0); |
| float2 p1 = evalCubic(t1); |
| CHECK(p0.x == Approx(p0ref.x).margin(1e-3f)); |
| CHECK(p0.y == Approx(p0ref.y).margin(1e-3f)); |
| CHECK(p1.x == Approx(p1ref.x).margin(1e-3f)); |
| CHECK(p1.y == Approx(p1ref.y).margin(1e-3f)); |
| } |
| |
| TEST_CASE("EvalCubic", "[bezier_utils]") |
| { |
| for (auto cubic : TEST_CUBICS) |
| { |
| math::EvalCubic evalCubic(cubic); |
| check_eval_cubic(evalCubic, cubic, 0, 1); |
| for (float t = 0; t <= 1; t += 1e-2f) |
| { |
| check_eval_cubic(evalCubic, cubic, t, t + 3e-3f); |
| } |
| } |
| } |
| |
| // glsl cross-compiling requires the clang/gcc vector extension. |
| #ifndef _MSC_VER |
| static void check_cubic_coeffs_tangents_glsl(const Vec2D* pts) |
| { |
| // Check cubic coefficients. |
| math::CubicCoeffs coeffs(pts); |
| float2 p0 = simd::load2f(pts + 0); |
| float2 p1 = simd::load2f(pts + 1); |
| float2 p2 = simd::load2f(pts + 2); |
| float2 p3 = simd::load2f(pts + 3); |
| float2 A, B, C; |
| glsl_cross::find_cubic_coeffs(p0, p1, p2, p3, A, B, C); |
| CHECK(simd::all(A == coeffs.A)); |
| CHECK(simd::all(B == coeffs.B)); |
| CHECK(simd::all(C == coeffs.C)); |
| |
| // Check cubic tangents. |
| Vec2D tangents[2]; |
| math::find_cubic_tangents(pts, tangents); |
| std::array<float2, 2> glslTangents = |
| glsl_cross::find_cubic_tangents(p0, p1, p2, p3); |
| CHECK(simd::all(simd::load2f(&tangents[0]) == glslTangents[0])); |
| CHECK(simd::all(simd::load2f(&tangents[1]) == glslTangents[1])); |
| } |
| |
| TEST_CASE("find_cubic_coeffs_tangents_glsl", "[bezier_utils]") |
| { |
| for (auto pts : TEST_CUBICS) |
| { |
| check_cubic_coeffs_tangents_glsl(pts); |
| } |
| // Test all combinations of corners from the square [0,0,1,1]. This covers |
| // every cubic type as well as a wide variety of special cases for cusps, |
| // lines, loops, and inflections. |
| for (int i = 0; i < (1 << 8); ++i) |
| { |
| Vec2D pts[4] = {Vec2D((i >> 0) & 1, (i >> 1) & 1) * 100, |
| Vec2D((i >> 2) & 1, (i >> 3) & 1) * 100, |
| Vec2D((i >> 4) & 1, (i >> 5) & 1) * 100, |
| Vec2D((i >> 6) & 1, (i >> 7) & 1) * 100}; |
| check_cubic_coeffs_tangents_glsl(pts); |
| } |
| Rand rando; |
| rando.seed(0); |
| for (int i = 0; i < 100; ++i) |
| { |
| Vec2D randos[] = { |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| }; |
| check_cubic_coeffs_tangents_glsl(randos); |
| } |
| } |
| |
| TEST_CASE("clamped_divide_glsl", "[bezier_utils]") |
| { |
| constexpr static float INF = std::numeric_limits<float>::infinity(); |
| CHECK(glsl_cross::clamped_divide(1, 2) == .5f); |
| CHECK(glsl_cross::clamped_divide(-2, 0) == 0); |
| CHECK(glsl_cross::clamped_divide(1, 0) == 1); |
| CHECK(glsl_cross::clamped_divide(-INF, 1) == 0); |
| CHECK(glsl_cross::clamped_divide(INF, 1) == 1); |
| CHECK(glsl_cross::clamped_divide(INF, INF) == 1); |
| CHECK(glsl_cross::clamped_divide(-INF, -INF) == 1); |
| CHECK(glsl_cross::clamped_divide(-INF, INF) == 0); |
| CHECK(glsl_cross::clamped_divide(INF, -INF) == 0); |
| CHECK(glsl_cross::clamped_divide(0, 0) == 0); |
| CHECK(glsl_cross::clamped_divide(1, NAN) == 1); |
| CHECK(glsl_cross::clamped_divide(-1, NAN) == 0); |
| CHECK(glsl_cross::clamped_divide(NAN, 1) == 0); |
| CHECK(glsl_cross::clamped_divide(NAN, -1) == 0); |
| CHECK(glsl_cross::clamped_divide(NAN, NAN) == 0); |
| } |
| |
| static void check_cubic_max_height_glsl(const Vec2D* pts) |
| { |
| float t, h = glsl_cross::find_cubic_max_height(simd::load2f(pts + 0), |
| simd::load2f(pts + 1), |
| simd::load2f(pts + 2), |
| simd::load2f(pts + 3), |
| t); |
| CHECK(h >= 0); |
| CHECK(t >= 0); |
| CHECK(t <= 1); |
| Vec2D norm = (pts[3] - pts[0]).normalized(); |
| norm = {-norm.y, norm.x}; |
| float k = -Vec2D::dot(norm, pts[0]); |
| auto heightAt = [=](float t) { |
| return fabsf(Vec2D::dot(norm, math::eval_cubic_at(pts, t)) + k); |
| }; |
| CHECK(heightAt(t) == Approx(h).margin(1e-4f)); |
| float epsilon = std::max(h, 1.f) * 1e-3f; |
| CHECK(h + epsilon > heightAt(0)); |
| CHECK(h + epsilon > heightAt(1)); |
| for (float t2 = 0; t2 <= 1; t2 += .005137f) |
| { |
| CHECK(h + epsilon > heightAt(t2)); |
| } |
| } |
| |
| TEST_CASE("find_cubic_max_height_glsl", "[bezier_utils]") |
| { |
| for (auto pts : TEST_CUBICS) |
| { |
| check_cubic_max_height_glsl(pts); |
| } |
| // Test all combinations of corners from the square [0,0,1,1]. This covers |
| // every cubic type as well as a wide variety of special cases for cusps, |
| // lines, loops, and inflections. |
| for (int i = 0; i < (1 << 8); ++i) |
| { |
| Vec2D pts[4] = {Vec2D((i >> 0) & 1, (i >> 1) & 1) * 100, |
| Vec2D((i >> 2) & 1, (i >> 3) & 1) * 100, |
| Vec2D((i >> 4) & 1, (i >> 5) & 1) * 100, |
| Vec2D((i >> 6) & 1, (i >> 7) & 1) * 100}; |
| check_cubic_max_height_glsl(pts); |
| } |
| Rand rando; |
| rando.seed(0); |
| for (int i = 0; i < 100; ++i) |
| { |
| Vec2D randos[] = { |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| }; |
| check_cubic_max_height_glsl(randos); |
| } |
| } |
| |
| static float guess_cubic_local_curvature(const Vec2D* p, |
| const CubicCoeffs& coeffs, |
| float t, |
| float desiredSpread) |
| { |
| // Iteratively guess a spread and compare with the computed result. |
| float maxDT = fminf(t, 1 - t) * .9999f; |
| float2 tan = 3.f * ((coeffs.A * t + 2.f * coeffs.B) * t + coeffs.C); |
| float lengthTan = sqrtf(simd::dot(tan, tan)); |
| tan /= lengthTan; |
| float dt = 0; |
| math::EvalCubic evalCubic(coeffs, p[0]); |
| // This would converge faster if we used the derivative of the Spread(dt) |
| // function from measure_cubic_local_curvature, but since the whole point of |
| // this test is to validate that function, use a binary search instead. |
| for (int i = 0; i < 24; ++i) |
| { |
| float guessDT = (dt + maxDT) * .5f; |
| float guessT0 = t - guessDT; |
| float guessT1 = t + guessDT; |
| float4 endpts = evalCubic(float4{guessT0, guessT0, guessT1, guessT1}); |
| float spread = fabsf(simd::dot(endpts.zw - endpts.xy, tan)); |
| if (spread <= desiredSpread) |
| dt = guessDT; |
| else |
| maxDT = guessDT; |
| } |
| Vec2D chops[10]; |
| math::chop_cubic_at(p, chops, t - dt, t + dt); |
| return measure_non_inflect_cubic_rotation(chops + 3); |
| } |
| |
| constexpr static float FEATHERING_CUSP_PADDING = 1e-3f; |
| |
| static int find_cubic_convex_90_chops(const Vec2D pts[], |
| float outT[4], |
| float cuspPadding, |
| bool* areCusps) |
| { |
| assert(pts); |
| assert(outT); |
| assert(areCusps); |
| |
| // Now find the cubic's inflection function. |
| // There are inflections where F' x F'' == 0. |
| // |
| // We formulate this as a quadratic equation: |
| // |
| // F' x F'' == a * T^2 + b * T + c == 0. |
| // |
| // See: |
| // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf |
| // NOTE: We only need the roots, so a uniform scale factor does not affect |
| // the solution. |
| CubicCoeffs coeffs(pts); |
| float a = simd::cross(coeffs.A, coeffs.B); |
| float b_over_2 = simd::cross(coeffs.A, coeffs.C) * .5f; |
| float c = simd::cross(coeffs.B, coeffs.C); |
| float discr_over_4 = b_over_2 * b_over_2 - a * c; |
| |
| // If -cuspThreshold <= discr_over_4 <= cuspThreshold, it means the two |
| // roots are within TESS_EPSILON of one another (in parametric space). This |
| // is close enough for our purposes to consider them a single cusp. |
| constexpr static float TESS_EPSILON = 1.f / (1 << 10); |
| float cuspThreshold = a * (TESS_EPSILON / 2); |
| cuspThreshold *= cuspThreshold; |
| |
| // Find the first two chops, based on curve classification. Also fill in |
| // "tan90", which will define the second pair of chops as the two points |
| // perpendicular to "tan90". |
| float4 T; |
| float2 tan90; |
| if (discr_over_4 < -cuspThreshold || |
| // Check if it's quadratic. |
| std::max(fabs(a), fabs(b_over_2)) < fabs(c) * TESS_EPSILON) |
| { |
| // The curve is a loop or quadratic. |
| // One chop is where rotation == 180 deg (which happens at infinity if |
| // the curve is quadratic). |
| // (This is the 2nd root where the tangent is parallel to tan0.) |
| // |
| // Tangent_Direction(T) x tan0 == 0 |
| // (AT^2 x tan0) + (2BT x tan0) + (C x tan0) == 0 |
| // (A x C)T^2 + (2B x C)T + (C x C) == 0 |
| // [[because tan0 == P1 - P0 == C]] |
| // bT^2 + 2cT + 0 == 0 [[because A x C == b, B x C == c]] |
| // T = [0, -2c/b] |
| // |
| // NOTE: if C == 0, then C != tan0. But this is fine because the curve |
| // can only rotate 180 degrees if the endpoints are colocated, and this |
| // gets handled next. |
| T.xy = {-c / b_over_2, 1}; |
| |
| // Next chop 90 degrees from the starting tangent of the curve. |
| tan90 = simd::any(coeffs.C != 0.f) |
| ? coeffs.C |
| : math::bit_cast<float2>(pts[2] - pts[0]); |
| *areCusps = false; |
| } |
| else if (discr_over_4 > cuspThreshold) |
| { |
| // The curve is serpentine. Solve for the two inflection points. |
| float q = sqrtf(discr_over_4); |
| q = -b_over_2 - copysignf(q, b_over_2); |
| T.xy = float2{q, c} / float2{a, q}; |
| |
| // Next chop 90 degrees from the whichever inflection point is closest |
| // to the middle. |
| float t = fabsf(T.x - .5f) < fabsf(T.y - .5f) ? T.x : T.y; |
| tan90 = (coeffs.A * t + 2.f * coeffs.B) * t + coeffs.C; |
| *areCusps = false; |
| } |
| else |
| { |
| // The curve is a cusp. A proper cusp is at T=-b/2a, but just solving |
| // for 90 degrees from the starting tangent will also find it, in |
| // addition to finding cusps from degenerate flat lines reversing |
| // direction. Since 180 degrees of rotation is lost to the cusp, we only |
| // need to find 2 roots max. |
| T.xy = 1; |
| tan90 = simd::any(coeffs.C != 0.f) |
| ? coeffs.C |
| : math::bit_cast<float2>(pts[2] - pts[0]); |
| *areCusps = true; |
| } |
| |
| // Find a second set of chops where the curve is perpendicular to tan90. |
| // |
| // Tangent_Direction(T) dot tan90 == 0 |
| // (A dot tan90) * T^2 + (2B dot tan90) * T + (C dot tan90) == 0 |
| // |
| a = simd::dot(coeffs.A, tan90); |
| b_over_2 = simd::dot(coeffs.B, tan90); |
| c = simd::dot(coeffs.C, tan90); |
| discr_over_4 = b_over_2 * b_over_2 - a * c; |
| float q = sqrtf(discr_over_4); |
| q = -b_over_2 - copysignf(q, b_over_2); |
| T.zw = float2{q, c} / float2{a, q}; |
| |
| // Throw out T <= epsilon and T >= epsilon by converting them to 1. |
| // (Use logic such that NaN also converts to 1.) |
| T = simd::if_then_else((T > 0) & (T < 1), T, float4(1)); |
| assert(simd::all(T > 0)); |
| assert(simd::all(T <= 1)); |
| |
| // Sort the roots. |
| T = simd::if_then_else((float2{T.x, T.z} < float2{T.y, T.w}).xxyy, |
| T, |
| T.yxwz); |
| T = simd::if_then_else((float2{T.x, T.y} < float2{T.z, T.w}).xyxy, |
| T, |
| T.zwxy); |
| T = T.y < T.z ? T : T.xzyw; |
| |
| // Count the number of roots that != 1 and store T. |
| int4 n4 = (T != 1) & 1; |
| n4.xy += n4.zw; |
| int n = n4.x + n4.y; |
| simd::store(outT, T); |
| |
| if (*areCusps && n > 0) |
| { |
| // Generate padding around cusp points. Odd numbered chops are always |
| // padding sections that pass through a cusp. |
| assert(n <= 2); |
| for (int i = n - 1; i >= 0; --i) |
| { |
| float maxT = i == n - 1 ? 1 : outT[i * 2 + 1]; |
| float minT = i == 0 ? 0 : (outT[i - 1] + outT[i]) * .5f; |
| outT[i * 2 + 1] = std::min(outT[i] + cuspPadding, maxT); |
| outT[i * 2 + 0] = std::max(outT[i] - cuspPadding, minT); |
| } |
| n *= 2; |
| // Re-clip and re-sort n after adding cusp padding. This is a hack, but |
| // we leave it here for now becase we're about to remove this entire |
| // method in favor of something more robust. |
| if (outT[n - 1] == 1) |
| { |
| --n; |
| } |
| std::sort(outT, outT + n); |
| } |
| |
| return n; |
| } |
| |
| static void check_cubic_convex_90_chops(const Vec2D* pts) |
| { |
| float T[4]; |
| bool areCusps; |
| int n = |
| find_cubic_convex_90_chops(pts, T, FEATHERING_CUSP_PADDING, &areCusps); |
| CHECK(n <= 4); |
| |
| Vec2D chops[16]; |
| assert(n * 3 + 1 <= std::size(chops)); |
| math::chop_cubic_at(pts, chops, T, n); |
| Vec2D* p = chops; |
| for (int i = 0; i <= n; ++i, p += 3) |
| { |
| if (areCusps && (i & 1)) |
| { |
| // If the chops are around a cusp, odd-numbered chops are a padding |
| // section that passes through a cusp. |
| continue; |
| } |
| CHECK(measure_non_inflect_cubic_rotation(p) <= |
| Approx(math::PI / 2).margin(1e-2f)); |
| } |
| } |
| |
| TEST_CASE("find_cubic_convex_90_chops", "[bezier_utils]") |
| { |
| for (auto pts : TEST_CUBICS) |
| { |
| check_cubic_convex_90_chops(pts); |
| } |
| // Test all combinations of corners from the square [0,0,1,1]. This covers |
| // every cubic type as well as a wide variety of special cases for cusps, |
| // lines, loops, and inflections. |
| for (int i = 0; i < (1 << 8); ++i) |
| { |
| Vec2D pts[4] = {Vec2D((i >> 0) & 1, (i >> 1) & 1) * 100, |
| Vec2D((i >> 2) & 1, (i >> 3) & 1) * 100, |
| Vec2D((i >> 4) & 1, (i >> 5) & 1) * 100, |
| Vec2D((i >> 6) & 1, (i >> 7) & 1) * 100}; |
| check_cubic_convex_90_chops(pts); |
| } |
| Rand rando; |
| rando.seed(0); |
| for (int i = 0; i < 100; ++i) |
| { |
| Vec2D randos[] = { |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| }; |
| check_cubic_convex_90_chops(randos); |
| } |
| } |
| |
| static void check_cubic_local_curvature_glsl(const Vec2D* pts) |
| { |
| float T[4]; |
| bool areCusps; |
| int n = |
| find_cubic_convex_90_chops(pts, T, FEATHERING_CUSP_PADDING, &areCusps); |
| CHECK(n <= 4); |
| |
| Vec2D chops[16]; |
| assert(n * 3 + 1 <= std::size(chops)); |
| math::chop_cubic_at(pts, chops, T, n); |
| Vec2D* p = chops; |
| for (int i = 0; i <= n; ++i, p += 3) |
| { |
| if (areCusps && (i & 1)) |
| { |
| // If the chops are around a cusp, odd-numbered chops are a padding |
| // section that passes through a cusp. |
| continue; |
| } |
| |
| float2 p0 = simd::load2f(p + 0); |
| float2 p1 = simd::load2f(p + 1); |
| float2 p2 = simd::load2f(p + 2); |
| float2 p3 = simd::load2f(p + 3); |
| |
| float maxHeightT; |
| CHECK(measure_non_inflect_cubic_rotation(p) <= |
| Approx(math::PI / 2).margin(2.6e-3f)); |
| glsl_cross::find_cubic_max_height(p0, p1, p2, p3, maxHeightT); |
| CubicCoeffs coeffs(p); |
| for (float desiredSpread : {1.f, 10.f, 100.f}) |
| { |
| for (float t : {maxHeightT, .5f}) |
| { |
| float theta = |
| glsl_cross::measure_cubic_local_curvature(p0, |
| p1, |
| p2, |
| p3, |
| t, |
| desiredSpread); |
| CHECK(theta >= 0); |
| CHECK(theta <= math::PI); |
| |
| float guess = |
| guess_cubic_local_curvature(p, coeffs, t, desiredSpread); |
| CHECK(theta == Approx(guess).margin(1e-2f)); |
| } |
| } |
| } |
| } |
| |
| TEST_CASE("measure_cubic_local_curvature_glsl", "[bezier_utils]") |
| { |
| for (auto pts : TEST_CUBICS) |
| { |
| check_cubic_local_curvature_glsl(pts); |
| } |
| // Test all combinations of corners from the square [0,0,1,1]. This covers |
| // every cubic type as well as a wide variety of special cases for cusps, |
| // lines, loops, and inflections. |
| for (int i = 0; i < (1 << 8); ++i) |
| { |
| Vec2D pts[4] = {Vec2D((i >> 0) & 1, (i >> 1) & 1) * 100, |
| Vec2D((i >> 2) & 1, (i >> 3) & 1) * 100, |
| Vec2D((i >> 4) & 1, (i >> 5) & 1) * 100, |
| Vec2D((i >> 6) & 1, (i >> 7) & 1) * 100}; |
| check_cubic_local_curvature_glsl(pts); |
| } |
| Rand rando; |
| rando.seed(0); |
| for (int i = 0; i < 100; ++i) |
| { |
| Vec2D randos[] = { |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| {rando.f32(-100, 100), rando.f32(-100, 100)}, |
| }; |
| check_cubic_local_curvature_glsl(randos); |
| } |
| } |
| #endif |
| } // namespace math |
| } // namespace rive |
