473e4ca5b8c793d4e60e5b11d25d96963c4802c8
[blender-staging.git] / extern / curve_fit_nd / intern / curve_fit_cubic.c
1 /*
2  * Copyright (c) 2016, DWANGO Co., Ltd.
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *     * Redistributions of source code must retain the above copyright
8  *       notice, this list of conditions and the following disclaimer.
9  *     * Redistributions in binary form must reproduce the above copyright
10  *       notice, this list of conditions and the following disclaimer in the
11  *       documentation and/or other materials provided with the distribution.
12  *     * Neither the name of the <organization> nor the
13  *       names of its contributors may be used to endorse or promote products
14  *       derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27
28 /** \file curve_fit_cubic.c
29  *  \ingroup curve_fit
30  */
31
32 #ifdef _MSC_VER
33 #  define _USE_MATH_DEFINES
34 #endif
35
36 #include <math.h>
37 #include <float.h>
38 #include <stdbool.h>
39 #include <assert.h>
40
41 #include <string.h>
42 #include <stdlib.h>
43
44 #include "../curve_fit_nd.h"
45
46 /* Take curvature into account when calculating the least square solution isn't usable. */
47 #define USE_CIRCULAR_FALLBACK
48
49 /* avoid re-calculating lengths multiple times */
50 #define USE_LENGTH_CACHE
51
52 /* store the indices in the cubic data so we can return the original indices,
53  * useful when the caller has data associated with the curve. */
54 #define USE_ORIG_INDEX_DATA
55
56 typedef unsigned int uint;
57
58 #include "curve_fit_inline.h"
59
60 #ifdef _MSC_VER
61 #  define alloca(size) _alloca(size)
62 #endif
63
64 #if !defined(_MSC_VER)
65 #  define USE_VLA
66 #endif
67
68 #ifdef USE_VLA
69 #  ifdef __GNUC__
70 #    pragma GCC diagnostic ignored "-Wvla"
71 #  endif
72 #else
73 #  ifdef __GNUC__
74 #    pragma GCC diagnostic error "-Wvla"
75 #  endif
76 #endif
77
78 #define SWAP(type, a, b)  {    \
79         type sw_ap;                \
80         sw_ap = (a);               \
81         (a) = (b);                 \
82         (b) = sw_ap;               \
83 } (void)0
84
85
86 /* -------------------------------------------------------------------- */
87
88 /** \name Cubic Type & Functions
89  * \{ */
90
91 typedef struct Cubic {
92         /* single linked lists */
93         struct Cubic *next;
94 #ifdef USE_ORIG_INDEX_DATA
95         uint orig_span;
96 #endif
97         /* 0: point_0, 1: handle_0, 2: handle_1, 3: point_1,
98          * each one is offset by 'dims' */
99         double pt_data[0];
100 } Cubic;
101
102 #define CUBIC_PT(cubic, index, dims) \
103         (&(cubic)->pt_data[(index) * (dims)])
104
105 #define CUBIC_VARS(c, dims, _p0, _p1, _p2, _p3) \
106         double \
107         *_p0 = (c)->pt_data, \
108         *_p1 = _p0 + (dims), \
109         *_p2 = _p1 + (dims), \
110         *_p3 = _p2 + (dims); ((void)0)
111 #define CUBIC_VARS_CONST(c, dims, _p0, _p1, _p2, _p3) \
112         const double \
113         *_p0 = (c)->pt_data, \
114         *_p1 = _p0 + (dims), \
115         *_p2 = _p1 + (dims), \
116         *_p3 = _p2 + (dims); ((void)0)
117
118
119 static size_t cubic_alloc_size(const uint dims)
120 {
121         return sizeof(Cubic) + (sizeof(double) * 4 * dims);
122 }
123
124 static Cubic *cubic_alloc(const uint dims)
125 {
126         return malloc(cubic_alloc_size(dims));
127 }
128
129 static void cubic_copy(Cubic *cubic_dst, const Cubic *cubic_src, const uint dims)
130 {
131         memcpy(cubic_dst, cubic_src, cubic_alloc_size(dims));
132 }
133
134 static void cubic_init(
135         Cubic *cubic,
136         const double p0[], const double p1[], const double p2[], const double p3[],
137         const uint dims)
138 {
139         copy_vnvn(CUBIC_PT(cubic, 0, dims), p0, dims);
140         copy_vnvn(CUBIC_PT(cubic, 1, dims), p1, dims);
141         copy_vnvn(CUBIC_PT(cubic, 2, dims), p2, dims);
142         copy_vnvn(CUBIC_PT(cubic, 3, dims), p3, dims);
143 }
144
145 static void cubic_free(Cubic *cubic)
146 {
147         free(cubic);
148 }
149
150 /** \} */
151
152
153 /* -------------------------------------------------------------------- */
154
155 /** \name CubicList Type & Functions
156  * \{ */
157
158 typedef struct CubicList {
159         struct Cubic *items;
160         uint          len;
161         uint          dims;
162 } CubicList;
163
164 static void cubic_list_prepend(CubicList *clist, Cubic *cubic)
165 {
166         cubic->next = clist->items;
167         clist->items = cubic;
168         clist->len++;
169 }
170
171 static double *cubic_list_as_array(
172         const CubicList *clist
173 #ifdef USE_ORIG_INDEX_DATA
174         ,
175         const uint index_last,
176         uint *r_orig_index
177 #endif
178         )
179 {
180         const uint dims = clist->dims;
181         const uint array_flat_len = (clist->len + 1) * 3 * dims;
182
183         double *array = malloc(sizeof(double) * array_flat_len);
184         const double *handle_prev = &((Cubic *)clist->items)->pt_data[dims];
185
186 #ifdef USE_ORIG_INDEX_DATA
187         uint orig_index_value = index_last;
188         uint orig_index_index = clist->len;
189         bool use_orig_index = (r_orig_index != NULL);
190 #endif
191
192         /* fill the array backwards */
193         const size_t array_chunk = 3 * dims;
194         double *array_iter = array + array_flat_len;
195         for (Cubic *citer = clist->items; citer; citer = citer->next) {
196                 array_iter -= array_chunk;
197                 memcpy(array_iter, &citer->pt_data[2 * dims], sizeof(double) * 2 * dims);
198                 memcpy(&array_iter[2 * dims], &handle_prev[dims], sizeof(double) * dims);
199                 handle_prev = citer->pt_data;
200
201 #ifdef USE_ORIG_INDEX_DATA
202                 if (use_orig_index) {
203                         r_orig_index[orig_index_index--] = orig_index_value;
204                         orig_index_value -= citer->orig_span;
205                 }
206 #endif
207         }
208
209 #ifdef USE_ORIG_INDEX_DATA
210         if (use_orig_index) {
211                 assert(orig_index_index == 0);
212                 assert(orig_index_value == 0 || index_last == 0);
213                 r_orig_index[orig_index_index] = index_last ? orig_index_value : 0;
214
215         }
216 #endif
217
218         /* flip tangent for first and last (we could leave at zero, but set to something useful) */
219
220         /* first */
221         array_iter -= array_chunk;
222         memcpy(&array_iter[dims], handle_prev, sizeof(double) * 2 * dims);
223         flip_vn_vnvn(&array_iter[0 * dims], &array_iter[1 * dims], &array_iter[2 * dims], dims);
224         assert(array == array_iter);
225
226         /* last */
227         array_iter += array_flat_len - (3 * dims);
228         flip_vn_vnvn(&array_iter[2 * dims], &array_iter[1 * dims], &array_iter[0 * dims], dims);
229
230         return array;
231 }
232
233 static void cubic_list_clear(CubicList *clist)
234 {
235         Cubic *cubic_next;
236         for (Cubic *citer = clist->items; citer; citer = cubic_next) {
237                 cubic_next = citer->next;
238                 cubic_free(citer);
239         }
240         clist->items = NULL;
241         clist->len  = 0;
242 }
243
244 /** \} */
245
246
247 /* -------------------------------------------------------------------- */
248
249 /** \name Cubic Evaluation
250  * \{ */
251
252 static void cubic_evaluate(
253         const Cubic *cubic, const double t, const uint dims,
254         double r_v[])
255 {
256         CUBIC_VARS_CONST(cubic, dims, p0, p1, p2, p3);
257         const double s = 1.0 - t;
258
259         for (uint j = 0; j < dims; j++) {
260                 const double p01 = (p0[j] * s) + (p1[j] * t);
261                 const double p12 = (p1[j] * s) + (p2[j] * t);
262                 const double p23 = (p2[j] * s) + (p3[j] * t);
263                 r_v[j] = ((((p01 * s) + (p12 * t))) * s) +
264                          ((((p12 * s) + (p23 * t))) * t);
265         }
266 }
267
268 static void cubic_calc_point(
269         const Cubic *cubic, const double t, const uint dims,
270         double r_v[])
271 {
272         CUBIC_VARS_CONST(cubic, dims, p0, p1, p2, p3);
273         const double s = 1.0 - t;
274         for (uint j = 0; j < dims; j++) {
275                 r_v[j] = p0[j] * s * s * s +
276                          3.0 * t * s * (s * p1[j] + t * p2[j]) + t * t * t * p3[j];
277         }
278 }
279
280 static void cubic_calc_speed(
281         const Cubic *cubic, const double t, const uint dims,
282         double r_v[])
283 {
284         CUBIC_VARS_CONST(cubic, dims, p0, p1, p2, p3);
285         const double s = 1.0 - t;
286         for (uint j = 0; j < dims; j++) {
287                 r_v[j] =  3.0 * ((p1[j] - p0[j]) * s * s + 2.0 *
288                                  (p2[j] - p0[j]) * s * t +
289                                  (p3[j] - p2[j]) * t * t);
290         }
291 }
292
293 static void cubic_calc_acceleration(
294         const Cubic *cubic, const double t, const uint dims,
295         double r_v[])
296 {
297         CUBIC_VARS_CONST(cubic, dims, p0, p1, p2, p3);
298         const double s = 1.0 - t;
299         for (uint j = 0; j < dims; j++) {
300                 r_v[j] = 6.0 * ((p2[j] - 2.0 * p1[j] + p0[j]) * s +
301                                 (p3[j] - 2.0 * p2[j] + p1[j]) * t);
302         }
303 }
304
305 /**
306  * Returns a 'measure' of the maximum distance (squared) of the points specified
307  * by points_offset from the corresponding cubic(u[]) points.
308  */
309 static double cubic_calc_error(
310         const Cubic *cubic,
311         const double *points_offset,
312         const uint points_offset_len,
313         const double *u,
314         const uint dims,
315
316         uint *r_error_index)
317 {
318         double error_max_sq = 0.0;
319         uint   error_index = 0;
320
321         const double *pt_real = points_offset + dims;
322 #ifdef USE_VLA
323         double        pt_eval[dims];
324 #else
325         double       *pt_eval = alloca(sizeof(double) * dims);
326 #endif
327
328         for (uint i = 1; i < points_offset_len - 1; i++, pt_real += dims) {
329                 cubic_evaluate(cubic, u[i], dims, pt_eval);
330
331                 const double err_sq = len_squared_vnvn(pt_real, pt_eval, dims);
332                 if (err_sq >= error_max_sq) {
333                         error_max_sq = err_sq;
334                         error_index = i;
335                 }
336         }
337
338         *r_error_index = error_index;
339         return error_max_sq;
340 }
341
342 /**
343  * Bezier multipliers
344  */
345
346 static double B1(double u)
347 {
348         double tmp = 1.0 - u;
349         return 3.0 * u * tmp * tmp;
350 }
351
352 static double B2(double u)
353 {
354         return 3.0 * u * u * (1.0 - u);
355 }
356
357 static double B0plusB1(double u)
358 {
359     double tmp = 1.0 - u;
360     return tmp * tmp * (1.0 + 2.0 * u);
361 }
362
363 static double B2plusB3(double u)
364 {
365     return u * u * (3.0 - 2.0 * u);
366 }
367
368 static void points_calc_center_weighted(
369         const double *points_offset,
370         const uint    points_offset_len,
371         const uint    dims,
372
373         double r_center[])
374 {
375         /*
376          * Calculate a center that compensates for point spacing.
377          */
378
379         const double *pt_prev = &points_offset[(points_offset_len - 2) * dims];
380         const double *pt_curr = pt_prev + dims;
381         const double *pt_next = points_offset;
382
383         double w_prev = len_vnvn(pt_prev, pt_curr, dims);
384
385         zero_vn(r_center, dims);
386         double w_tot = 0.0;
387
388         for (uint i_next = 0; i_next < points_offset_len; i_next++) {
389                 const double w_next = len_vnvn(pt_curr, pt_next, dims);
390                 const double w = w_prev + w_next;
391                 w_tot += w;
392
393                 miadd_vn_vn_fl(r_center, pt_curr, w, dims);
394
395                 w_prev = w_next;
396
397                 pt_prev = pt_curr;
398                 pt_curr = pt_next;
399                 pt_next += dims;
400         }
401
402         if (w_tot != 0.0) {
403                 imul_vn_fl(r_center, 1.0 / w_tot, dims);
404         }
405 }
406
407 #ifdef USE_CIRCULAR_FALLBACK
408
409 /**
410  * Return a scale value, used to calculate how much the curve handles should be increased,
411  *
412  * This works by placing each end-point on an imaginary circle,
413  * the placement on the circle is based on the tangent vectors,
414  * where larger differences in tangent angle cover a larger part of the circle.
415  *
416  * Return the scale representing how much larger the distance around the circle is.
417  */
418 static double points_calc_circumference_factor(
419         const double  tan_l[],
420         const double  tan_r[],
421         const uint dims)
422 {
423         const double dot = dot_vnvn(tan_l, tan_r, dims);
424         const double len_tangent = dot < 0.0 ? len_vnvn(tan_l, tan_r, dims) : len_negated_vnvn(tan_l, tan_r, dims);
425         if (len_tangent > DBL_EPSILON) {
426                 /* only clamp to avoid precision error */
427                 double angle = acos(max(-fabs(dot), -1.0));
428                 /* Angle may be less than the length when the tangents define >180 degrees of the circle,
429                  * (tangents that point away from each other).
430                  * We could try support this but will likely cause extreme >1 scales which could cause other issues. */
431                 // assert(angle >= len_tangent);
432                 double factor = (angle / len_tangent);
433                 assert(factor < (M_PI / 2) + DBL_EPSILON);
434                 return factor;
435         }
436         else {
437                 /* tangents are exactly aligned (think two opposite sides of a circle). */
438                 return (M_PI / 2);
439         }
440 }
441
442 /**
443  * Return the value which the distance between points will need to be scaled by,
444  * to define a handle, given both points are on a perfect circle.
445  *
446  * \note the return value will need to be multiplied by 1.3... for correct results.
447  */
448 static double points_calc_circle_tangent_factor(
449         const double  tan_l[],
450         const double  tan_r[],
451         const uint dims)
452 {
453         const double angle_sin = len_vnvn(tan_l, tan_r, dims) / 2.0;
454         if (angle_sin != 0.0) {
455                 const double tan_dot = dot_vnvn(tan_l, tan_r, dims);
456                 double scale;
457                 if (tan_dot > -1.0) {
458                         const double angle = acos(tan_dot) / 2.0;
459                         const double angle_cos = cos(angle);
460                         scale = (1.0 - angle_cos) / (angle_sin * 2.0);
461                 }
462                 else {
463                         scale = 1.0 / 2.0;
464                 }
465                 return (scale / angle_sin);
466         }
467         else {
468                 return (1.0 / 3.0) * 0.75;
469         }
470 }
471
472 /**
473  * Calculate the scale the handles, which serves as a best-guess
474  * used as a fallback when the least-square solution fails.
475  */
476 static double points_calc_cubic_scale(
477         const double v_l[], const double v_r[],
478         const double  tan_l[],
479         const double  tan_r[],
480         const double coords_length, uint dims)
481 {
482         const double len_direct = len_vnvn(v_l, v_r, dims);
483         const double len_circle_factor = points_calc_circle_tangent_factor(tan_l, tan_r, dims);
484
485         /* if this curve is a circle, this value doesn't need modification */
486         const double len_circle_handle = (len_direct * (len_circle_factor / 0.75));
487
488         /* scale by the difference from the circumference distance */
489         const double len_circle = len_direct * points_calc_circumference_factor(tan_l, tan_r, dims);
490         double scale_handle = (coords_length / len_circle);
491
492         /* Could investigate an accurate calculation here,
493          * though this gives close results */
494         scale_handle = ((scale_handle - 1.0) * 1.75) + 1.0;
495
496         return len_circle_handle * scale_handle;
497 }
498
499 static void cubic_from_points_fallback(
500         const double *points_offset,
501         const uint    points_offset_len,
502         const double  tan_l[],
503         const double  tan_r[],
504         const uint dims,
505
506         Cubic *r_cubic)
507 {
508         const double *p0 = &points_offset[0];
509         const double *p3 = &points_offset[(points_offset_len - 1) * dims];
510
511         double alpha = len_vnvn(p0, p3, dims) / 3.0;
512
513         double *p1 = CUBIC_PT(r_cubic, 1, dims);
514         double *p2 = CUBIC_PT(r_cubic, 2, dims);
515
516         copy_vnvn(CUBIC_PT(r_cubic, 0, dims), p0, dims);
517         copy_vnvn(CUBIC_PT(r_cubic, 3, dims), p3, dims);
518
519 #ifdef USE_ORIG_INDEX_DATA
520         r_cubic->orig_span = (points_offset_len - 1);
521 #endif
522
523         /* p1 = p0 - (tan_l * alpha_l);
524          * p2 = p3 + (tan_r * alpha_r);
525          */
526         msub_vn_vnvn_fl(p1, p0, tan_l, alpha, dims);
527         madd_vn_vnvn_fl(p2, p3, tan_r, alpha, dims);
528 }
529 #endif  /* USE_CIRCULAR_FALLBACK */
530
531 /**
532  * Use least-squares method to find Bezier control points for region.
533  */
534 static void cubic_from_points(
535         const double *points_offset,
536         const uint    points_offset_len,
537 #ifdef USE_CIRCULAR_FALLBACK
538         const double  points_offset_coords_length,
539 #endif
540         const double *u_prime,
541         const double  tan_l[],
542         const double  tan_r[],
543         const uint dims,
544
545         Cubic *r_cubic)
546 {
547
548         const double *p0 = &points_offset[0];
549         const double *p3 = &points_offset[(points_offset_len - 1) * dims];
550
551         /* Point Pairs */
552         double alpha_l, alpha_r;
553 #ifdef USE_VLA
554         double a[2][dims];
555         double tmp[dims];
556 #else
557         double *a[2] = {
558             alloca(sizeof(double) * dims),
559             alloca(sizeof(double) * dims),
560         };
561         double *tmp = alloca(sizeof(double) * dims);
562 #endif
563
564         {
565                 double x[2] = {0.0}, c[2][2] = {{0.0}};
566                 const double *pt = points_offset;
567
568                 for (uint i = 0; i < points_offset_len; i++, pt += dims) {
569                         mul_vnvn_fl(a[0], tan_l, B1(u_prime[i]), dims);
570                         mul_vnvn_fl(a[1], tan_r, B2(u_prime[i]), dims);
571
572                         c[0][0] += dot_vnvn(a[0], a[0], dims);
573                         c[0][1] += dot_vnvn(a[0], a[1], dims);
574                         c[1][1] += dot_vnvn(a[1], a[1], dims);
575
576                         c[1][0] = c[0][1];
577
578                         {
579                                 const double b0_plus_b1 = B0plusB1(u_prime[i]);
580                                 const double b2_plus_b3 = B2plusB3(u_prime[i]);
581                                 for (uint j = 0; j < dims; j++) {
582                                         tmp[j] = (pt[j] - (p0[j] * b0_plus_b1)) + (p3[j] * b2_plus_b3);
583                                 }
584
585                                 x[0] += dot_vnvn(a[0], tmp, dims);
586                                 x[1] += dot_vnvn(a[1], tmp, dims);
587                         }
588                 }
589
590                 double det_C0_C1 = c[0][0] * c[1][1] - c[0][1] * c[1][0];
591                 double det_C_0X  = x[1]    * c[0][0] - x[0]    * c[0][1];
592                 double det_X_C1  = x[0]    * c[1][1] - x[1]    * c[0][1];
593
594                 if (is_almost_zero(det_C0_C1)) {
595                         det_C0_C1 = c[0][0] * c[1][1] * 10e-12;
596                 }
597
598                 /* may still divide-by-zero, check below will catch nan values */
599                 alpha_l = det_X_C1 / det_C0_C1;
600                 alpha_r = det_C_0X / det_C0_C1;
601         }
602
603         /*
604          * The problem that the stupid values for alpha dare not put
605          * only when we realize that the sign and wrong,
606          * but even if the values are too high.
607          * But how do you evaluate it?
608          *
609          * Meanwhile, we should ensure that these values are sometimes
610          * so only problems absurd of approximation and not for bugs in the code.
611          */
612
613         bool use_clamp = true;
614
615         /* flip check to catch nan values */
616         if (!(alpha_l >= 0.0) ||
617             !(alpha_r >= 0.0))
618         {
619 #ifdef USE_CIRCULAR_FALLBACK
620                 alpha_l = alpha_r = points_calc_cubic_scale(p0, p3, tan_l, tan_r, points_offset_coords_length, dims);
621 #else
622                 alpha_l = alpha_r = len_vnvn(p0, p3, dims) / 3.0;
623 #endif
624
625                 /* skip clamping when we're using default handles */
626                 use_clamp = false;
627         }
628
629         double *p1 = CUBIC_PT(r_cubic, 1, dims);
630         double *p2 = CUBIC_PT(r_cubic, 2, dims);
631
632         copy_vnvn(CUBIC_PT(r_cubic, 0, dims), p0, dims);
633         copy_vnvn(CUBIC_PT(r_cubic, 3, dims), p3, dims);
634
635 #ifdef USE_ORIG_INDEX_DATA
636         r_cubic->orig_span = (points_offset_len - 1);
637 #endif
638
639         /* p1 = p0 - (tan_l * alpha_l);
640          * p2 = p3 + (tan_r * alpha_r);
641          */
642         msub_vn_vnvn_fl(p1, p0, tan_l, alpha_l, dims);
643         madd_vn_vnvn_fl(p2, p3, tan_r, alpha_r, dims);
644
645         /* ------------------------------------
646          * Clamping (we could make it optional)
647          */
648         if (use_clamp) {
649 #ifdef USE_VLA
650                 double center[dims];
651 #else
652                 double *center = alloca(sizeof(double) * dims);
653 #endif
654                 points_calc_center_weighted(points_offset, points_offset_len, dims, center);
655
656                 const double clamp_scale = 3.0;  /* clamp to 3x */
657                 double dist_sq_max = 0.0;
658
659                 {
660                         const double *pt = points_offset;
661                         for (uint i = 0; i < points_offset_len; i++, pt += dims) {
662 #if 0
663                                 double dist_sq_test = sq(len_vnvn(center, pt, dims) * clamp_scale);
664 #else
665                                 /* do inline */
666                                 double dist_sq_test = 0.0;
667                                 for (uint j = 0; j < dims; j++) {
668                                         dist_sq_test += sq((pt[j] - center[j]) * clamp_scale);
669                                 }
670 #endif
671                                 dist_sq_max = max(dist_sq_max, dist_sq_test);
672                         }
673                 }
674
675                 double p1_dist_sq = len_squared_vnvn(center, p1, dims);
676                 double p2_dist_sq = len_squared_vnvn(center, p2, dims);
677
678                 if (p1_dist_sq > dist_sq_max ||
679                     p2_dist_sq > dist_sq_max)
680                 {
681 #ifdef USE_CIRCULAR_FALLBACK
682                         alpha_l = alpha_r = points_calc_cubic_scale(p0, p3, tan_l, tan_r, points_offset_coords_length, dims);
683 #else
684                         alpha_l = alpha_r = len_vnvn(p0, p3, dims) / 3.0;
685 #endif
686
687                         /*
688                          * p1 = p0 - (tan_l * alpha_l);
689                          * p2 = p3 + (tan_r * alpha_r);
690                          */
691                         for (uint j = 0; j < dims; j++) {
692                                 p1[j] = p0[j] - (tan_l[j] * alpha_l);
693                                 p2[j] = p3[j] + (tan_r[j] * alpha_r);
694                         }
695
696                         p1_dist_sq = len_squared_vnvn(center, p1, dims);
697                         p2_dist_sq = len_squared_vnvn(center, p2, dims);
698                 }
699
700                 /* clamp within the 3x radius */
701                 if (p1_dist_sq > dist_sq_max) {
702                         isub_vnvn(p1, center, dims);
703                         imul_vn_fl(p1, sqrt(dist_sq_max) / sqrt(p1_dist_sq), dims);
704                         iadd_vnvn(p1, center, dims);
705                 }
706                 if (p2_dist_sq > dist_sq_max) {
707                         isub_vnvn(p2, center, dims);
708                         imul_vn_fl(p2, sqrt(dist_sq_max) / sqrt(p2_dist_sq), dims);
709                         iadd_vnvn(p2, center, dims);
710                 }
711         }
712         /* end clamping */
713 }
714
715 #ifdef USE_LENGTH_CACHE
716 static void points_calc_coord_length_cache(
717         const double *points_offset,
718         const uint    points_offset_len,
719         const uint    dims,
720
721         double     *r_points_length_cache)
722 {
723         const double *pt_prev = points_offset;
724         const double *pt = pt_prev + dims;
725         r_points_length_cache[0] = 0.0;
726         for (uint i = 1; i < points_offset_len; i++) {
727                 r_points_length_cache[i] = len_vnvn(pt, pt_prev, dims);
728                 pt_prev = pt;
729                 pt += dims;
730         }
731 }
732 #endif  /* USE_LENGTH_CACHE */
733
734 /**
735  * \return the accumulated length of \a points_offset.
736  */
737 static double points_calc_coord_length(
738         const double *points_offset,
739         const uint    points_offset_len,
740         const uint    dims,
741 #ifdef USE_LENGTH_CACHE
742         const double *points_length_cache,
743 #endif
744         double *r_u)
745 {
746         const double *pt_prev = points_offset;
747         const double *pt = pt_prev + dims;
748         r_u[0] = 0.0;
749         for (uint i = 1, i_prev = 0; i < points_offset_len; i++) {
750                 double length;
751
752 #ifdef USE_LENGTH_CACHE
753                 length = points_length_cache[i];
754
755                 assert(len_vnvn(pt, pt_prev, dims) == points_length_cache[i]);
756 #else
757                 length = len_vnvn(pt, pt_prev, dims);
758 #endif
759
760                 r_u[i] = r_u[i_prev] + length;
761                 i_prev = i;
762                 pt_prev = pt;
763                 pt += dims;
764         }
765         assert(!is_almost_zero(r_u[points_offset_len - 1]));
766         const double w = r_u[points_offset_len - 1];
767         for (uint i = 0; i < points_offset_len; i++) {
768                 r_u[i] /= w;
769         }
770         return w;
771 }
772
773 /**
774  * Use Newton-Raphson iteration to find better root.
775  *
776  * \param cubic: Current fitted curve.
777  * \param p: Point to test against.
778  * \param u: Parameter value for \a p.
779  *
780  * \note Return value may be `nan` caller must check for this.
781  */
782 static double cubic_find_root(
783         const Cubic *cubic,
784         const double p[],
785         const double u,
786         const uint dims)
787 {
788         /* Newton-Raphson Method. */
789         /* all vectors */
790 #ifdef USE_VLA
791         double q0_u[dims];
792         double q1_u[dims];
793         double q2_u[dims];
794 #else
795         double *q0_u = alloca(sizeof(double) * dims);
796         double *q1_u = alloca(sizeof(double) * dims);
797         double *q2_u = alloca(sizeof(double) * dims);
798 #endif
799
800         cubic_calc_point(cubic, u, dims, q0_u);
801         cubic_calc_speed(cubic, u, dims, q1_u);
802         cubic_calc_acceleration(cubic, u, dims, q2_u);
803
804         /* may divide-by-zero, caller must check for that case */
805         /* u - ((q0_u - p) * q1_u) / (q1_u.length_squared() + (q0_u - p) * q2_u) */
806         isub_vnvn(q0_u, p, dims);
807         return u - dot_vnvn(q0_u, q1_u, dims) /
808                (len_squared_vn(q1_u, dims) + dot_vnvn(q0_u, q2_u, dims));
809 }
810
811 static int compare_double_fn(const void *a_, const void *b_)
812 {
813         const double *a = a_;
814         const double *b = b_;
815         if      (*a > *b) return  1;
816         else if (*a < *b) return -1;
817         else              return  0;
818 }
819
820 /**
821  * Given set of points and their parameterization, try to find a better parameterization.
822  */
823 static bool cubic_reparameterize(
824         const Cubic *cubic,
825         const double *points_offset,
826         const uint    points_offset_len,
827         const double *u,
828         const uint    dims,
829
830         double       *r_u_prime)
831 {
832         /*
833          * Recalculate the values of u[] based on the Newton Raphson method
834          */
835
836         const double *pt = points_offset;
837         for (uint i = 0; i < points_offset_len; i++, pt += dims) {
838                 r_u_prime[i] = cubic_find_root(cubic, pt, u[i], dims);
839                 if (!isfinite(r_u_prime[i])) {
840                         return false;
841                 }
842         }
843
844         qsort(r_u_prime, points_offset_len, sizeof(double), compare_double_fn);
845
846         if ((r_u_prime[0] < 0.0) ||
847             (r_u_prime[points_offset_len - 1] > 1.0))
848         {
849                 return false;
850         }
851
852         assert(r_u_prime[0] >= 0.0);
853         assert(r_u_prime[points_offset_len - 1] <= 1.0);
854         return true;
855 }
856
857
858 static bool fit_cubic_to_points(
859         const double *points_offset,
860         const uint    points_offset_len,
861 #ifdef USE_LENGTH_CACHE
862         const double *points_length_cache,
863 #endif
864         const double  tan_l[],
865         const double  tan_r[],
866         const double  error_threshold_sq,
867         const uint    dims,
868
869         Cubic *r_cubic, double *r_error_max_sq, uint *r_split_index)
870 {
871         const uint iteration_max = 4;
872
873         if (points_offset_len == 2) {
874                 CUBIC_VARS(r_cubic, dims, p0, p1, p2, p3);
875
876                 copy_vnvn(p0, &points_offset[0 * dims], dims);
877                 copy_vnvn(p3, &points_offset[1 * dims], dims);
878
879                 const double dist = len_vnvn(p0, p3, dims) / 3.0;
880                 msub_vn_vnvn_fl(p1, p0, tan_l, dist, dims);
881                 madd_vn_vnvn_fl(p2, p3, tan_r, dist, dims);
882
883 #ifdef USE_ORIG_INDEX_DATA
884                 r_cubic->orig_span = 1;
885 #endif
886                 return true;
887         }
888
889         double *u = malloc(sizeof(double) * points_offset_len);
890
891 #ifdef USE_CIRCULAR_FALLBACK
892         const double points_offset_coords_length  =
893 #endif
894         points_calc_coord_length(
895                 points_offset, points_offset_len, dims,
896 #ifdef USE_LENGTH_CACHE
897                 points_length_cache,
898 #endif
899                 u);
900
901         double error_max_sq;
902         uint split_index;
903
904         /* Parameterize points, and attempt to fit curve */
905         cubic_from_points(
906                 points_offset, points_offset_len,
907 #ifdef USE_CIRCULAR_FALLBACK
908                 points_offset_coords_length,
909 #endif
910                 u, tan_l, tan_r, dims, r_cubic);
911
912         /* Find max deviation of points to fitted curve */
913         error_max_sq = cubic_calc_error(
914                 r_cubic, points_offset, points_offset_len, u, dims,
915                 &split_index);
916
917         Cubic *cubic_test = alloca(cubic_alloc_size(dims));
918
919 #ifdef USE_CIRCULAR_FALLBACK
920         if (!(error_max_sq < error_threshold_sq)) {
921                 /* Don't use the cubic calculated above, instead calculate a new fallback cubic,
922                  * since this tends to give more balanced split_index along the curve.
923                  * This is because the attempt to calcualte the cubic may contain spikes
924                  * along the curve which may give a lop-sided maximum distance. */
925                 cubic_from_points_fallback(
926                         points_offset, points_offset_len,
927                         tan_l, tan_r, dims, cubic_test);
928                 const double error_max_sq_test = cubic_calc_error(
929                         cubic_test, points_offset, points_offset_len, u, dims,
930                         &split_index);
931
932                 /* intentionally use the newly calculated 'split_index',
933                  * even if the 'error_max_sq_test' is worse. */
934                 if (error_max_sq > error_max_sq_test) {
935                         error_max_sq = error_max_sq_test;
936                         cubic_copy(r_cubic, cubic_test, dims);
937                 }
938         }
939 #endif
940
941         *r_error_max_sq = error_max_sq;
942         *r_split_index  = split_index;
943
944         if (error_max_sq < error_threshold_sq) {
945                 free(u);
946                 return true;
947         }
948         else {
949                 cubic_copy(cubic_test, r_cubic, dims);
950
951                 /* If error not too large, try some reparameterization and iteration */
952                 double *u_prime = malloc(sizeof(double) * points_offset_len);
953                 for (uint iter = 0; iter < iteration_max; iter++) {
954                         if (!cubic_reparameterize(
955                                 cubic_test, points_offset, points_offset_len, u, dims, u_prime))
956                         {
957                                 break;
958                         }
959
960                         cubic_from_points(
961                                 points_offset, points_offset_len,
962 #ifdef USE_CIRCULAR_FALLBACK
963                                 points_offset_coords_length,
964 #endif
965                                 u_prime, tan_l, tan_r, dims, cubic_test);
966                         error_max_sq = cubic_calc_error(
967                                 cubic_test, points_offset, points_offset_len, u_prime, dims,
968                                 &split_index);
969
970                         if (error_max_sq < error_threshold_sq) {
971                                 free(u_prime);
972                                 free(u);
973
974                                 cubic_copy(r_cubic, cubic_test, dims);
975                                 *r_error_max_sq = error_max_sq;
976                                 *r_split_index  = split_index;
977                                 return true;
978                         }
979                         else if (error_max_sq < *r_error_max_sq) {
980                                 cubic_copy(r_cubic, cubic_test, dims);
981                                 *r_error_max_sq = error_max_sq;
982                                 *r_split_index = split_index;
983                         }
984
985                         SWAP(double *, u, u_prime);
986                 }
987                 free(u_prime);
988                 free(u);
989
990                 return false;
991         }
992 }
993
994 static void fit_cubic_to_points_recursive(
995         const double *points_offset,
996         const uint    points_offset_len,
997 #ifdef USE_LENGTH_CACHE
998         const double *points_length_cache,
999 #endif
1000         const double  tan_l[],
1001         const double  tan_r[],
1002         const double  error_threshold_sq,
1003         const uint    dims,
1004         /* fill in the list */
1005         CubicList *clist)
1006 {
1007         Cubic *cubic = cubic_alloc(dims);
1008         uint split_index;
1009         double error_max_sq;
1010
1011         if (fit_cubic_to_points(
1012                 points_offset, points_offset_len,
1013 #ifdef USE_LENGTH_CACHE
1014                 points_length_cache,
1015 #endif
1016                 tan_l, tan_r, error_threshold_sq, dims,
1017                 cubic, &error_max_sq, &split_index))
1018         {
1019                 cubic_list_prepend(clist, cubic);
1020                 return;
1021         }
1022         cubic_free(cubic);
1023
1024
1025         /* Fitting failed -- split at max error point and fit recursively */
1026
1027         /* Check splinePoint is not an endpoint?
1028          *
1029          * This assert happens sometimes...
1030          * Look into it but disable for now. Campbell! */
1031
1032         // assert(split_index > 1)
1033 #ifdef USE_VLA
1034         double tan_center[dims];
1035 #else
1036         double *tan_center = alloca(sizeof(double) * dims);
1037 #endif
1038
1039         const double *pt_a = &points_offset[(split_index - 1) * dims];
1040         const double *pt_b = &points_offset[(split_index + 1) * dims];
1041
1042         assert(split_index < points_offset_len);
1043         if (equals_vnvn(pt_a, pt_b, dims)) {
1044                 pt_a += dims;
1045         }
1046
1047         {
1048 #ifdef USE_VLA
1049                 double tan_center_a[dims];
1050                 double tan_center_b[dims];
1051 #else
1052                 double *tan_center_a = alloca(sizeof(double) * dims);
1053                 double *tan_center_b = alloca(sizeof(double) * dims);
1054 #endif
1055                 const double *pt   = &points_offset[split_index * dims];
1056
1057                 /* tan_center = ((pt_a - pt).normalized() + (pt - pt_b).normalized()).normalized() */
1058                 normalize_vn_vnvn(tan_center_a, pt_a, pt, dims);
1059                 normalize_vn_vnvn(tan_center_b, pt, pt_b, dims);
1060                 add_vn_vnvn(tan_center, tan_center_a, tan_center_b, dims);
1061                 normalize_vn(tan_center, dims);
1062         }
1063
1064         fit_cubic_to_points_recursive(
1065                 points_offset, split_index + 1,
1066 #ifdef USE_LENGTH_CACHE
1067                 points_length_cache,
1068 #endif
1069                 tan_l, tan_center, error_threshold_sq, dims, clist);
1070         fit_cubic_to_points_recursive(
1071                 &points_offset[split_index * dims], points_offset_len - split_index,
1072 #ifdef USE_LENGTH_CACHE
1073                 points_length_cache + split_index,
1074 #endif
1075                 tan_center, tan_r, error_threshold_sq, dims, clist);
1076
1077 }
1078
1079 /** \} */
1080
1081
1082 /* -------------------------------------------------------------------- */
1083
1084 /** \name External API for Curve-Fitting
1085  * \{ */
1086
1087 /**
1088  * Main function:
1089  *
1090  * Take an array of 3d points.
1091  * return the cubic splines
1092  */
1093 int curve_fit_cubic_to_points_db(
1094         const double *points,
1095         const uint    points_len,
1096         const uint    dims,
1097         const double  error_threshold,
1098         const uint   *corners,
1099         uint          corners_len,
1100
1101         double **r_cubic_array, uint *r_cubic_array_len,
1102         uint **r_cubic_orig_index,
1103         uint **r_corner_index_array, uint *r_corner_index_len)
1104 {
1105         uint corners_buf[2];
1106         if (corners == NULL) {
1107                 assert(corners_len == 0);
1108                 corners_buf[0] = 0;
1109                 corners_buf[1] = points_len - 1;
1110                 corners = corners_buf;
1111                 corners_len = 2;
1112         }
1113
1114         CubicList clist = {0};
1115         clist.dims = dims;
1116
1117 #ifdef USE_VLA
1118         double tan_l[dims];
1119         double tan_r[dims];
1120 #else
1121         double *tan_l = alloca(sizeof(double) * dims);
1122         double *tan_r = alloca(sizeof(double) * dims);
1123 #endif
1124
1125 #ifdef USE_LENGTH_CACHE
1126         double *points_length_cache = NULL;
1127         uint    points_length_cache_len_alloc = 0;
1128 #endif
1129
1130         uint *corner_index_array = NULL;
1131         uint  corner_index = 0;
1132         if (r_corner_index_array && (corners != corners_buf)) {
1133                 corner_index_array = malloc(sizeof(uint) * corners_len);
1134                 corner_index_array[corner_index++] = corners[0];
1135         }
1136
1137         const double error_threshold_sq = sq(error_threshold);
1138
1139         for (uint i = 1; i < corners_len; i++) {
1140                 const uint points_offset_len = corners[i] - corners[i - 1] + 1;
1141                 const uint first_point = corners[i - 1];
1142
1143                 assert(points_offset_len >= 1);
1144                 if (points_offset_len > 1) {
1145                         const double *pt_l = &points[first_point * dims];
1146                         const double *pt_r = &points[(first_point + points_offset_len - 1) * dims];
1147                         const double *pt_l_next = pt_l + dims;
1148                         const double *pt_r_prev = pt_r - dims;
1149
1150                         /* tan_l = (pt_l - pt_l_next).normalized()
1151                          * tan_r = (pt_r_prev - pt_r).normalized()
1152                          */
1153                         normalize_vn_vnvn(tan_l, pt_l, pt_l_next, dims);
1154                         normalize_vn_vnvn(tan_r, pt_r_prev, pt_r, dims);
1155
1156 #ifdef USE_LENGTH_CACHE
1157                         if (points_length_cache_len_alloc < points_offset_len) {
1158                                 if (points_length_cache) {
1159                                         free(points_length_cache);
1160                                 }
1161                                 points_length_cache = malloc(sizeof(double) * points_offset_len);
1162                         }
1163                         points_calc_coord_length_cache(
1164                                 &points[first_point * dims], points_offset_len, dims,
1165                                 points_length_cache);
1166 #endif
1167
1168                         fit_cubic_to_points_recursive(
1169                                 &points[first_point * dims], points_offset_len,
1170 #ifdef USE_LENGTH_CACHE
1171                                 points_length_cache,
1172 #endif
1173                                 tan_l, tan_r, error_threshold_sq, dims, &clist);
1174                 }
1175                 else if (points_len == 1) {
1176                         assert(points_offset_len == 1);
1177                         assert(corners_len == 2);
1178                         assert(corners[0] == 0);
1179                         assert(corners[1] == 0);
1180                         const double *pt = &points[0];
1181                         Cubic *cubic = cubic_alloc(dims);
1182                         cubic_init(cubic, pt, pt, pt, pt, dims);
1183                         cubic_list_prepend(&clist, cubic);
1184                 }
1185
1186                 if (corner_index_array) {
1187                         corner_index_array[corner_index++] = clist.len;
1188                 }
1189         }
1190
1191 #ifdef USE_LENGTH_CACHE
1192         if (points_length_cache) {
1193                 free(points_length_cache);
1194         }
1195 #endif
1196
1197 #ifdef USE_ORIG_INDEX_DATA
1198         uint *cubic_orig_index = NULL;
1199         if (r_cubic_orig_index) {
1200                 cubic_orig_index = malloc(sizeof(uint) * (clist.len + 1));
1201         }
1202 #else
1203         *r_cubic_orig_index = NULL;
1204 #endif
1205
1206         /* allocate a contiguous array and free the linked list */
1207         *r_cubic_array = cubic_list_as_array(
1208                 &clist
1209 #ifdef USE_ORIG_INDEX_DATA
1210                 , corners[corners_len - 1], cubic_orig_index
1211 #endif
1212                 );
1213         *r_cubic_array_len = clist.len + 1;
1214
1215         cubic_list_clear(&clist);
1216
1217 #ifdef USE_ORIG_INDEX_DATA
1218         if (cubic_orig_index) {
1219                 *r_cubic_orig_index = cubic_orig_index;
1220         }
1221 #endif
1222
1223         if (corner_index_array) {
1224                 assert(corner_index == corners_len);
1225                 *r_corner_index_array = corner_index_array;
1226                 *r_corner_index_len = corner_index;
1227         }
1228
1229         return 0;
1230 }
1231
1232 /**
1233  * A version of #curve_fit_cubic_to_points_db to handle floats
1234  */
1235 int curve_fit_cubic_to_points_fl(
1236         const float  *points,
1237         const uint    points_len,
1238         const uint    dims,
1239         const float   error_threshold,
1240         const uint   *corners,
1241         const uint    corners_len,
1242
1243         float **r_cubic_array, uint *r_cubic_array_len,
1244         uint **r_cubic_orig_index,
1245         uint **r_corner_index_array, uint *r_corner_index_len)
1246 {
1247         const uint points_flat_len = points_len * dims;
1248         double *points_db = malloc(sizeof(double) * points_flat_len);
1249
1250         copy_vndb_vnfl(points_db, points, points_flat_len);
1251
1252         double *cubic_array_db = NULL;
1253         float  *cubic_array_fl = NULL;
1254         uint    cubic_array_len = 0;
1255
1256         int result = curve_fit_cubic_to_points_db(
1257                 points_db, points_len, dims, error_threshold, corners, corners_len,
1258                 &cubic_array_db, &cubic_array_len,
1259                 r_cubic_orig_index,
1260                 r_corner_index_array, r_corner_index_len);
1261         free(points_db);
1262
1263         if (!result) {
1264                 uint cubic_array_flat_len = cubic_array_len * 3 * dims;
1265                 cubic_array_fl = malloc(sizeof(float) * cubic_array_flat_len);
1266                 for (uint i = 0; i < cubic_array_flat_len; i++) {
1267                         cubic_array_fl[i] = (float)cubic_array_db[i];
1268                 }
1269                 free(cubic_array_db);
1270         }
1271
1272         *r_cubic_array = cubic_array_fl;
1273         *r_cubic_array_len = cubic_array_len;
1274
1275         return result;
1276 }
1277
1278 /**
1279  * Fit a single cubic to points.
1280  */
1281 int curve_fit_cubic_to_points_single_db(
1282         const double *points,
1283         const uint    points_len,
1284         const uint    dims,
1285         const double  error_threshold,
1286         const double tan_l[],
1287         const double tan_r[],
1288
1289         double  r_handle_l[],
1290         double  r_handle_r[],
1291         double  *r_error_max_sq)
1292 {
1293         Cubic *cubic = alloca(cubic_alloc_size(dims));
1294
1295         uint split_index;
1296
1297         /* in this instance theres no advantage in using length cache,
1298          * since we're not recursively calculating values. */
1299 #ifdef USE_LENGTH_CACHE
1300         double *points_length_cache = malloc(sizeof(double) * points_len);
1301         points_calc_coord_length_cache(
1302                 points, points_len, dims,
1303                 points_length_cache);
1304 #endif
1305
1306         fit_cubic_to_points(
1307                 points, points_len,
1308 #ifdef USE_LENGTH_CACHE
1309                 points_length_cache,
1310 #endif
1311                 tan_l, tan_r, error_threshold, dims,
1312
1313                 cubic, r_error_max_sq, &split_index);
1314
1315 #ifdef USE_LENGTH_CACHE
1316         free(points_length_cache);
1317 #endif
1318
1319         copy_vnvn(r_handle_l, CUBIC_PT(cubic, 1, dims), dims);
1320         copy_vnvn(r_handle_r, CUBIC_PT(cubic, 2, dims), dims);
1321
1322         return 0;
1323 }
1324
1325 int curve_fit_cubic_to_points_single_fl(
1326         const float  *points,
1327         const uint    points_len,
1328         const uint    dims,
1329         const float   error_threshold,
1330         const float   tan_l[],
1331         const float   tan_r[],
1332
1333         float   r_handle_l[],
1334         float   r_handle_r[],
1335         float  *r_error_sq)
1336 {
1337         const uint points_flat_len = points_len * dims;
1338         double *points_db = malloc(sizeof(double) * points_flat_len);
1339
1340         copy_vndb_vnfl(points_db, points, points_flat_len);
1341
1342 #ifdef USE_VLA
1343         double tan_l_db[dims];
1344         double tan_r_db[dims];
1345         double r_handle_l_db[dims];
1346         double r_handle_r_db[dims];
1347 #else
1348         double *tan_l_db = alloca(sizeof(double) * dims);
1349         double *tan_r_db = alloca(sizeof(double) * dims);
1350         double *r_handle_l_db = alloca(sizeof(double) * dims);
1351         double *r_handle_r_db = alloca(sizeof(double) * dims);
1352 #endif
1353         double r_error_sq_db;
1354
1355         copy_vndb_vnfl(tan_l_db, tan_l, dims);
1356         copy_vndb_vnfl(tan_r_db, tan_r, dims);
1357
1358         int result = curve_fit_cubic_to_points_single_db(
1359                 points_db, points_len, dims,
1360                 (double)error_threshold,
1361                 tan_l_db, tan_r_db,
1362                 r_handle_l_db, r_handle_r_db,
1363                 &r_error_sq_db);
1364
1365         free(points_db);
1366
1367         copy_vnfl_vndb(r_handle_l, r_handle_l_db, dims);
1368         copy_vnfl_vndb(r_handle_r, r_handle_r_db, dims);
1369         *r_error_sq = (float)r_error_sq_db;
1370
1371         return result;
1372 }
1373
1374 /** \} */