Cleanup: whicespace
[blender.git] / extern / curve_fit_nd / intern / curve_fit_cubic.c
index 810cf92760d539198865142518e32e14cfb5d7a1..f07bb73429fdc0713f088210f3e580d0ddaa7808 100644 (file)
 
 #include "../curve_fit_nd.h"
 
+/* Take curvature into account when calculating the least square solution isn't usable. */
+#define USE_CIRCULAR_FALLBACK
+
 /* avoid re-calculating lengths multiple times */
 #define USE_LENGTH_CACHE
 
 /* store the indices in the cubic data so we can return the original indices,
- * useful when the caller has data assosiated with the curve. */
+ * useful when the caller has data associated with the curve. */
 #define USE_ORIG_INDEX_DATA
 
 typedef unsigned int uint;
@@ -109,9 +112,19 @@ typedef struct Cubic {
        *_p3 = _p2 + (dims); ((void)0)
 
 
+static size_t cubic_alloc_size(const uint dims)
+{
+       return sizeof(Cubic) + (sizeof(double) * 4 * dims);
+}
+
 static Cubic *cubic_alloc(const uint dims)
 {
-       return malloc(sizeof(Cubic) + (sizeof(double) * 4 * dims));
+       return malloc(cubic_alloc_size(dims));
+}
+
+static void cubic_copy(Cubic *cubic_dst, const Cubic *cubic_src, const uint dims)
+{
+       memcpy(cubic_dst, cubic_src, cubic_alloc_size(dims));
 }
 
 static void cubic_init(
@@ -278,7 +291,7 @@ static void cubic_calc_acceleration(
         double r_v[])
 {
        CUBIC_VARS_CONST(cubic, dims, p0, p1, p2, p3);
-    const double s = 1.0 - t;
+       const double s = 1.0 - t;
        for (uint j = 0; j < dims; j++) {
                r_v[j] = 6.0 * ((p2[j] - 2.0 * p1[j] + p0[j]) * s +
                                (p3[j] - 2.0 * p2[j] + p1[j]) * t);
@@ -286,20 +299,19 @@ static void cubic_calc_acceleration(
 }
 
 /**
- * Returns a 'measure' of the maximal discrepancy of the points specified
+ * Returns a 'measure' of the maximum distance (squared) of the points specified
  * by points_offset from the corresponding cubic(u[]) points.
  */
-static void cubic_calc_error(
+static double cubic_calc_error(
         const Cubic *cubic,
         const double *points_offset,
         const uint points_offset_len,
         const double *u,
         const uint dims,
 
-        double *r_error_sq_max,
         uint *r_error_index)
 {
-       double error_sq_max = 0.0;
+       double error_max_sq = 0.0;
        uint   error_index = 0;
 
        const double *pt_real = points_offset + dims;
@@ -313,14 +325,14 @@ static void cubic_calc_error(
                cubic_evaluate(cubic, u[i], dims, pt_eval);
 
                const double err_sq = len_squared_vnvn(pt_real, pt_eval, dims);
-               if (err_sq >= error_sq_max) {
-                       error_sq_max = err_sq;
+               if (err_sq >= error_max_sq) {
+                       error_max_sq = err_sq;
                        error_index = i;
                }
        }
 
-       *r_error_sq_max   = error_sq_max;
        *r_error_index = error_index;
+       return error_max_sq;
 }
 
 /**
@@ -388,12 +400,98 @@ static void points_calc_center_weighted(
        }
 }
 
+#ifdef USE_CIRCULAR_FALLBACK
+
+/**
+ * Return a scale value, used to calculate how much the curve handles should be increased,
+ *
+ * This works by placing each end-point on an imaginary circle,
+ * the placement on the circle is based on the tangent vectors,
+ * where larger differences in tangent angle cover a larger part of the circle.
+ *
+ * Return the scale representing how much larger the distance around the circle is.
+ */
+static double points_calc_circumference_factor(
+        const double  tan_l[],
+        const double  tan_r[],
+        const uint dims)
+{
+       const double dot = dot_vnvn(tan_l, tan_r, dims);
+       const double len_tangent = dot < 0.0 ? len_vnvn(tan_l, tan_r, dims) : len_negated_vnvn(tan_l, tan_r, dims);
+       if (len_tangent > DBL_EPSILON) {
+               double angle = acos(-fabs(dot));
+               /* Angle may be less than the length when the tangents define >180 degrees of the circle,
+                * (tangents that point away from each other).
+                * We could try support this but will likely cause extreme >1 scales which could cause other issues. */
+               // assert(angle >= len_tangent);
+               double factor = (angle / len_tangent) / (M_PI / 2);
+               factor = 1.0 - pow(1.0 - factor, 1.75);
+               assert(factor < 1.0 + DBL_EPSILON);
+               return factor;
+       }
+       else {
+               /* tangents are exactly aligned (think two opposite sides of a circle). */
+               return 1.0;
+       }
+}
+
+/**
+ * Calculate the scale the handles, which serves as a best-guess
+ * used as a fallback when the least-square solution fails.
+ */
+static double points_calc_cubic_scale(
+        const double v_l[], const double v_r[],
+        const double  tan_l[],
+        const double  tan_r[],
+        const double coords_length, uint dims)
+{
+       const double len_direct = len_vnvn(v_l, v_r, dims);
+       const double len_circle_factor = points_calc_circumference_factor(tan_l, tan_r, dims) * 1.75;
+       const double len_points = min(coords_length, len_circle_factor * len_direct);
+       return (len_direct + ((len_points - len_direct) * len_circle_factor)) / 3.0;
+}
+
+static void cubic_from_points_fallback(
+        const double *points_offset,
+        const uint    points_offset_len,
+        const double  tan_l[],
+        const double  tan_r[],
+        const uint dims,
+
+        Cubic *r_cubic)
+{
+       const double *p0 = &points_offset[0];
+       const double *p3 = &points_offset[(points_offset_len - 1) * dims];
+
+       double alpha = len_vnvn(p0, p3, dims) / 3.0;
+
+       double *p1 = CUBIC_PT(r_cubic, 1, dims);
+       double *p2 = CUBIC_PT(r_cubic, 2, dims);
+
+       copy_vnvn(CUBIC_PT(r_cubic, 0, dims), p0, dims);
+       copy_vnvn(CUBIC_PT(r_cubic, 3, dims), p3, dims);
+
+#ifdef USE_ORIG_INDEX_DATA
+       r_cubic->orig_span = (points_offset_len - 1);
+#endif
+
+       /* p1 = p0 - (tan_l * alpha_l);
+        * p2 = p3 + (tan_r * alpha_r);
+        */
+       msub_vn_vnvn_fl(p1, p0, tan_l, alpha, dims);
+       madd_vn_vnvn_fl(p2, p3, tan_r, alpha, dims);
+}
+#endif  /* USE_CIRCULAR_FALLBACK */
+
 /**
  * Use least-squares method to find Bezier control points for region.
  */
 static void cubic_from_points(
         const double *points_offset,
         const uint    points_offset_len,
+#ifdef USE_CIRCULAR_FALLBACK
+        const double  points_offset_coords_length,
+#endif
         const double *u_prime,
         const double  tan_l[],
         const double  tan_r[],
@@ -467,11 +565,20 @@ static void cubic_from_points(
         * so only problems absurd of approximation and not for bugs in the code.
         */
 
+       bool use_clamp = true;
+
        /* flip check to catch nan values */
        if (!(alpha_l >= 0.0) ||
            !(alpha_r >= 0.0))
        {
+#ifdef USE_CIRCULAR_FALLBACK
+               alpha_l = alpha_r = points_calc_cubic_scale(p0, p3, tan_l, tan_r, points_offset_coords_length, dims);
+#else
                alpha_l = alpha_r = len_vnvn(p0, p3, dims) / 3.0;
+#endif
+
+               /* skip clamping when we're using default handles */
+               use_clamp = false;
        }
 
        double *p1 = CUBIC_PT(r_cubic, 1, dims);
@@ -493,64 +600,69 @@ static void cubic_from_points(
        /* ------------------------------------
         * Clamping (we could make it optional)
         */
+       if (use_clamp) {
 #ifdef USE_VLA
-       double center[dims];
+               double center[dims];
 #else
-       double *center = alloca(sizeof(double) * dims);
+               double *center = alloca(sizeof(double) * dims);
 #endif
-       points_calc_center_weighted(points_offset, points_offset_len, dims, center);
+               points_calc_center_weighted(points_offset, points_offset_len, dims, center);
 
-       const double clamp_scale = 3.0;  /* clamp to 3x */
-       double dist_sq_max = 0.0;
+               const double clamp_scale = 3.0;  /* clamp to 3x */
+               double dist_sq_max = 0.0;
 
-       {
-               const double *pt = points_offset;
-               for (uint i = 0; i < points_offset_len; i++, pt += dims) {
+               {
+                       const double *pt = points_offset;
+                       for (uint i = 0; i < points_offset_len; i++, pt += dims) {
 #if 0
-                       double dist_sq_test = sq(len_vnvn(center, pt, dims) * clamp_scale);
+                               double dist_sq_test = sq(len_vnvn(center, pt, dims) * clamp_scale);
 #else
-                       /* do inline */
-                       double dist_sq_test = 0.0;
-                       for (uint j = 0; j < dims; j++) {
-                               dist_sq_test += sq((pt[j] - center[j]) * clamp_scale);
-                       }
+                               /* do inline */
+                               double dist_sq_test = 0.0;
+                               for (uint j = 0; j < dims; j++) {
+                                       dist_sq_test += sq((pt[j] - center[j]) * clamp_scale);
+                               }
 #endif
-                       dist_sq_max = max(dist_sq_max, dist_sq_test);
+                               dist_sq_max = max(dist_sq_max, dist_sq_test);
+                       }
                }
-       }
 
-       double p1_dist_sq = len_squared_vnvn(center, p1, dims);
-       double p2_dist_sq = len_squared_vnvn(center, p2, dims);
+               double p1_dist_sq = len_squared_vnvn(center, p1, dims);
+               double p2_dist_sq = len_squared_vnvn(center, p2, dims);
 
-       if (p1_dist_sq > dist_sq_max ||
-           p2_dist_sq > dist_sq_max)
-       {
+               if (p1_dist_sq > dist_sq_max ||
+                   p2_dist_sq > dist_sq_max)
+               {
+#ifdef USE_CIRCULAR_FALLBACK
+                       alpha_l = alpha_r = points_calc_cubic_scale(p0, p3, tan_l, tan_r, points_offset_coords_length, dims);
+#else
+                       alpha_l = alpha_r = len_vnvn(p0, p3, dims) / 3.0;
+#endif
 
-               alpha_l = alpha_r = len_vnvn(p0, p3, dims) / 3.0;
+                       /*
+                        * p1 = p0 - (tan_l * alpha_l);
+                        * p2 = p3 + (tan_r * alpha_r);
+                        */
+                       for (uint j = 0; j < dims; j++) {
+                               p1[j] = p0[j] - (tan_l[j] * alpha_l);
+                               p2[j] = p3[j] + (tan_r[j] * alpha_r);
+                       }
 
-               /*
-                * p1 = p0 - (tan_l * alpha_l);
-                * p2 = p3 + (tan_r * alpha_r);
-                */
-               for (uint j = 0; j < dims; j++) {
-                       p1[j] = p0[j] - (tan_l[j] * alpha_l);
-                       p2[j] = p3[j] + (tan_r[j] * alpha_r);
+                       p1_dist_sq = len_squared_vnvn(center, p1, dims);
+                       p2_dist_sq = len_squared_vnvn(center, p2, dims);
                }
 
-               p1_dist_sq = len_squared_vnvn(center, p1, dims);
-               p2_dist_sq = len_squared_vnvn(center, p2, dims);
-       }
-
-       /* clamp within the 3x radius */
-       if (p1_dist_sq > dist_sq_max) {
-               isub_vnvn(p1, center, dims);
-               imul_vn_fl(p1, sqrt(dist_sq_max) / sqrt(p1_dist_sq), dims);
-               iadd_vnvn(p1, center, dims);
-       }
-       if (p2_dist_sq > dist_sq_max) {
-               isub_vnvn(p2, center, dims);
-               imul_vn_fl(p2, sqrt(dist_sq_max) / sqrt(p2_dist_sq), dims);
-               iadd_vnvn(p2, center, dims);
+               /* clamp within the 3x radius */
+               if (p1_dist_sq > dist_sq_max) {
+                       isub_vnvn(p1, center, dims);
+                       imul_vn_fl(p1, sqrt(dist_sq_max) / sqrt(p1_dist_sq), dims);
+                       iadd_vnvn(p1, center, dims);
+               }
+               if (p2_dist_sq > dist_sq_max) {
+                       isub_vnvn(p2, center, dims);
+                       imul_vn_fl(p2, sqrt(dist_sq_max) / sqrt(p2_dist_sq), dims);
+                       iadd_vnvn(p2, center, dims);
+               }
        }
        /* end clamping */
 }
@@ -574,8 +686,10 @@ static void points_calc_coord_length_cache(
 }
 #endif  /* USE_LENGTH_CACHE */
 
-
-static void points_calc_coord_length(
+/**
+ * \return the accumulated length of \a points_offset.
+ */
+static double points_calc_coord_length(
         const double *points_offset,
         const uint    points_offset_len,
         const uint    dims,
@@ -608,6 +722,7 @@ static void points_calc_coord_length(
        for (uint i = 0; i < points_offset_len; i++) {
                r_u[i] /= w;
        }
+       return w;
 }
 
 /**
@@ -620,10 +735,10 @@ static void points_calc_coord_length(
  * \note Return value may be `nan` caller must check for this.
  */
 static double cubic_find_root(
-               const Cubic *cubic,
-               const double p[],
-               const double u,
-               const uint dims)
+        const Cubic *cubic,
+        const double p[],
+        const double u,
+        const uint dims)
 {
        /* Newton-Raphson Method. */
        /* all vectors */
@@ -695,7 +810,7 @@ static bool cubic_reparameterize(
 }
 
 
-static void fit_cubic_to_points(
+static bool fit_cubic_to_points(
         const double *points_offset,
         const uint    points_offset_len,
 #ifdef USE_LENGTH_CACHE
@@ -703,19 +818,15 @@ static void fit_cubic_to_points(
 #endif
         const double  tan_l[],
         const double  tan_r[],
-        const double  error_threshold,
+        const double  error_threshold_sq,
         const uint    dims,
-        /* fill in the list */
-        CubicList *clist)
+
+        Cubic *r_cubic, double *r_error_max_sq, uint *r_split_index)
 {
        const uint iteration_max = 4;
-       const double error_sq = sq(error_threshold);
-
-       Cubic *cubic;
 
        if (points_offset_len == 2) {
-               cubic = cubic_alloc(dims);
-               CUBIC_VARS(cubic, dims, p0, p1, p2, p3);
+               CUBIC_VARS(r_cubic, dims, p0, p1, p2, p3);
 
                copy_vnvn(p0, &points_offset[0 * dims], dims);
                copy_vnvn(p3, &points_offset[1 * dims], dims);
@@ -725,14 +836,16 @@ static void fit_cubic_to_points(
                madd_vn_vnvn_fl(p2, p3, tan_r, dist, dims);
 
 #ifdef USE_ORIG_INDEX_DATA
-               cubic->orig_span = 1;
+               r_cubic->orig_span = 1;
 #endif
-
-               cubic_list_prepend(clist, cubic);
-               return;
+               return true;
        }
 
        double *u = malloc(sizeof(double) * points_offset_len);
+
+#ifdef USE_CIRCULAR_FALLBACK
+       const double points_offset_coords_length  =
+#endif
        points_calc_coord_length(
                points_offset, points_offset_len, dims,
 #ifdef USE_LENGTH_CACHE
@@ -740,55 +853,127 @@ static void fit_cubic_to_points(
 #endif
                u);
 
-       cubic = cubic_alloc(dims);
-
-       double error_sq_max;
+       double error_max_sq;
        uint split_index;
 
        /* Parameterize points, and attempt to fit curve */
        cubic_from_points(
-               points_offset, points_offset_len, u, tan_l, tan_r, dims, cubic);
+               points_offset, points_offset_len,
+#ifdef USE_CIRCULAR_FALLBACK
+               points_offset_coords_length,
+#endif
+               u, tan_l, tan_r, dims, r_cubic);
 
        /* Find max deviation of points to fitted curve */
-       cubic_calc_error(
-               cubic, points_offset, points_offset_len, u, dims,
-               &error_sq_max, &split_index);
+       error_max_sq = cubic_calc_error(
+               r_cubic, points_offset, points_offset_len, u, dims,
+               &split_index);
+
+       Cubic *cubic_test = alloca(cubic_alloc_size(dims));
+
+#ifdef USE_CIRCULAR_FALLBACK
+       if (!(error_max_sq < error_threshold_sq)) {
+               /* Don't use the cubic calculated above, instead calculate a new fallback cubic,
+                * since this tends to give more balanced split_index along the curve.
+                * This is because the attempt to calcualte the cubic may contain spikes
+                * along the curve which may give a lop-sided maximum distance. */
+               cubic_from_points_fallback(
+                       points_offset, points_offset_len,
+                       tan_l, tan_r, dims, cubic_test);
+               const double error_max_sq_test = cubic_calc_error(
+                       cubic_test, points_offset, points_offset_len, u, dims,
+                       &split_index);
+
+               /* intentionally use the newly calculated 'split_index',
+                * even if the 'error_max_sq_test' is worse. */
+               if (error_max_sq > error_max_sq_test) {
+                       error_max_sq = error_max_sq_test;
+                       cubic_copy(r_cubic, cubic_test, dims);
+               }
+       }
+#endif
 
-       if (error_sq_max < error_sq) {
+       *r_error_max_sq = error_max_sq;
+       *r_split_index  = split_index;
+
+       if (error_max_sq < error_threshold_sq) {
                free(u);
-               cubic_list_prepend(clist, cubic);
-               return;
+               return true;
        }
        else {
+               cubic_copy(cubic_test, r_cubic, dims);
+
                /* If error not too large, try some reparameterization and iteration */
                double *u_prime = malloc(sizeof(double) * points_offset_len);
                for (uint iter = 0; iter < iteration_max; iter++) {
                        if (!cubic_reparameterize(
-                               cubic, points_offset, points_offset_len, u, dims, u_prime))
+                               cubic_test, points_offset, points_offset_len, u, dims, u_prime))
                        {
                                break;
                        }
 
                        cubic_from_points(
-                               points_offset, points_offset_len, u_prime,
-                               tan_l, tan_r, dims, cubic);
-                       cubic_calc_error(
-                               cubic, points_offset, points_offset_len, u_prime, dims,
-                               &error_sq_max, &split_index);
+                               points_offset, points_offset_len,
+#ifdef USE_CIRCULAR_FALLBACK
+                               points_offset_coords_length,
+#endif
+                               u_prime, tan_l, tan_r, dims, cubic_test);
+                       error_max_sq = cubic_calc_error(
+                               cubic_test, points_offset, points_offset_len, u_prime, dims,
+                               &split_index);
 
-                       if (error_sq_max < error_sq) {
+                       if (error_max_sq < error_threshold_sq) {
                                free(u_prime);
                                free(u);
-                               cubic_list_prepend(clist, cubic);
-                               return;
+
+                               cubic_copy(r_cubic, cubic_test, dims);
+                               *r_error_max_sq = error_max_sq;
+                               *r_split_index  = split_index;
+                               return true;
+                       }
+                       else if (error_max_sq < *r_error_max_sq) {
+                               cubic_copy(r_cubic, cubic_test, dims);
+                               *r_error_max_sq = error_max_sq;
+                               *r_split_index = split_index;
                        }
 
                        SWAP(double *, u, u_prime);
                }
                free(u_prime);
+               free(u);
+
+               return false;
        }
+}
+
+static void fit_cubic_to_points_recursive(
+        const double *points_offset,
+        const uint    points_offset_len,
+#ifdef USE_LENGTH_CACHE
+        const double *points_length_cache,
+#endif
+        const double  tan_l[],
+        const double  tan_r[],
+        const double  error_threshold_sq,
+        const uint    dims,
+        /* fill in the list */
+        CubicList *clist)
+{
+       Cubic *cubic = cubic_alloc(dims);
+       uint split_index;
+       double error_max_sq;
 
-       free(u);
+       if (fit_cubic_to_points(
+               points_offset, points_offset_len,
+#ifdef USE_LENGTH_CACHE
+               points_length_cache,
+#endif
+               tan_l, tan_r, error_threshold_sq, dims,
+               cubic, &error_max_sq, &split_index))
+       {
+               cubic_list_prepend(clist, cubic);
+               return;
+       }
        cubic_free(cubic);
 
 
@@ -814,21 +999,35 @@ static void fit_cubic_to_points(
                pt_a += dims;
        }
 
-       /* tan_center = (pt_a - pt_b).normalized() */
-       normalize_vn_vnvn(tan_center, pt_a, pt_b, dims);
+       {
+#ifdef USE_VLA
+               double tan_center_a[dims];
+               double tan_center_b[dims];
+#else
+               double *tan_center_a = alloca(sizeof(double) * dims);
+               double *tan_center_b = alloca(sizeof(double) * dims);
+#endif
+               const double *pt   = &points_offset[split_index * dims];
 
-       fit_cubic_to_points(
+               /* tan_center = ((pt_a - pt).normalized() + (pt - pt_b).normalized()).normalized() */
+               normalize_vn_vnvn(tan_center_a, pt_a, pt, dims);
+               normalize_vn_vnvn(tan_center_b, pt, pt_b, dims);
+               add_vn_vnvn(tan_center, tan_center_a, tan_center_b, dims);
+               normalize_vn(tan_center, dims);
+       }
+
+       fit_cubic_to_points_recursive(
                points_offset, split_index + 1,
 #ifdef USE_LENGTH_CACHE
                points_length_cache,
 #endif
-               tan_l, tan_center, error_threshold, dims, clist);
-       fit_cubic_to_points(
+               tan_l, tan_center, error_threshold_sq, dims, clist);
+       fit_cubic_to_points_recursive(
                &points_offset[split_index * dims], points_offset_len - split_index,
 #ifdef USE_LENGTH_CACHE
                points_length_cache + split_index,
 #endif
-               tan_center, tan_r, error_threshold, dims, clist);
+               tan_center, tan_r, error_threshold_sq, dims, clist);
 
 }
 
@@ -890,6 +1089,8 @@ int curve_fit_cubic_to_points_db(
                corner_index_array[corner_index++] = corners[0];
        }
 
+       const double error_threshold_sq = sq(error_threshold);
+
        for (uint i = 1; i < corners_len; i++) {
                const uint points_offset_len = corners[i] - corners[i - 1] + 1;
                const uint first_point = corners[i - 1];
@@ -919,12 +1120,12 @@ int curve_fit_cubic_to_points_db(
                                points_length_cache);
 #endif
 
-                       fit_cubic_to_points(
+                       fit_cubic_to_points_recursive(
                                &points[first_point * dims], points_offset_len,
 #ifdef USE_LENGTH_CACHE
                                points_length_cache,
 #endif
-                               tan_l, tan_r, error_threshold, dims, &clist);
+                               tan_l, tan_r, error_threshold_sq, dims, &clist);
                }
                else if (points_len == 1) {
                        assert(points_offset_len == 1);
@@ -1001,9 +1202,7 @@ int curve_fit_cubic_to_points_fl(
        const uint points_flat_len = points_len * dims;
        double *points_db = malloc(sizeof(double) * points_flat_len);
 
-       for (uint i = 0; i < points_flat_len; i++) {
-               points_db[i] = (double)points[i];
-       }
+       copy_vndb_vnfl(points_db, points, points_flat_len);
 
        double *cubic_array_db = NULL;
        float  *cubic_array_fl = NULL;
@@ -1031,4 +1230,100 @@ int curve_fit_cubic_to_points_fl(
        return result;
 }
 
+/**
+ * Fit a single cubic to points.
+ */
+int curve_fit_cubic_to_points_single_db(
+        const double *points,
+        const uint    points_len,
+        const uint    dims,
+        const double  error_threshold,
+        const double tan_l[],
+        const double tan_r[],
+
+        double  r_handle_l[],
+        double  r_handle_r[],
+        double  *r_error_max_sq)
+{
+       Cubic *cubic = alloca(cubic_alloc_size(dims));
+
+       uint split_index;
+
+       /* in this instance theres no advantage in using length cache,
+        * since we're not recursively calculating values. */
+#ifdef USE_LENGTH_CACHE
+       double *points_length_cache = malloc(sizeof(double) * points_len);
+       points_calc_coord_length_cache(
+               points, points_len, dims,
+               points_length_cache);
+#endif
+
+       fit_cubic_to_points(
+               points, points_len,
+#ifdef USE_LENGTH_CACHE
+               points_length_cache,
+#endif
+               tan_l, tan_r, error_threshold, dims,
+
+               cubic, r_error_max_sq, &split_index);
+
+#ifdef USE_LENGTH_CACHE
+       free(points_length_cache);
+#endif
+
+       copy_vnvn(r_handle_l, CUBIC_PT(cubic, 1, dims), dims);
+       copy_vnvn(r_handle_r, CUBIC_PT(cubic, 2, dims), dims);
+
+       return 0;
+}
+
+int curve_fit_cubic_to_points_single_fl(
+        const float  *points,
+        const uint    points_len,
+        const uint    dims,
+        const float   error_threshold,
+        const float   tan_l[],
+        const float   tan_r[],
+
+        float   r_handle_l[],
+        float   r_handle_r[],
+        float  *r_error_sq)
+{
+       const uint points_flat_len = points_len * dims;
+       double *points_db = malloc(sizeof(double) * points_flat_len);
+
+       copy_vndb_vnfl(points_db, points, points_flat_len);
+
+#ifdef USE_VLA
+       double tan_l_db[dims];
+       double tan_r_db[dims];
+       double r_handle_l_db[dims];
+       double r_handle_r_db[dims];
+#else
+       double *tan_l_db = alloca(sizeof(double) * dims);
+       double *tan_r_db = alloca(sizeof(double) * dims);
+       double *r_handle_l_db = alloca(sizeof(double) * dims);
+       double *r_handle_r_db = alloca(sizeof(double) * dims);
+#endif
+       double r_error_sq_db;
+
+       copy_vndb_vnfl(tan_l_db, tan_l, dims);
+       copy_vndb_vnfl(tan_r_db, tan_r, dims);
+
+       int result = curve_fit_cubic_to_points_single_db(
+               points_db, points_len, dims,
+               (double)error_threshold,
+               tan_l_db, tan_r_db,
+               r_handle_l_db, r_handle_r_db,
+               &r_error_sq_db);
+
+       free(points_db);
+
+       copy_vnfl_vndb(r_handle_l, r_handle_l_db, dims);
+       copy_vnfl_vndb(r_handle_r, r_handle_r_db, dims);
+       *r_error_sq = (float)r_error_sq_db;
+
+       return result;
+}
+
 /** \} */