HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_Spline.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_Spline.h ( UT Library, C++)
7  *
8  * COMMENTS: Simple spline class.
9  *
10  * The linear and catmull-rom splines expect a parametric evaluation coordinate
11  * between 0 and 1.
12  */
13 
14 #ifndef __UT_Spline__
15 #define __UT_Spline__
16 
17 #include "UT_API.h"
18 #include "UT_Array.h"
19 #include "UT_BoundingBox.h"
20 #include "UT_Color.h"
21 #include "UT_Matrix4.h"
22 #include "UT_RootFinder.h"
23 #include "UT_Vector3.h"
24 #include "UT_Vector4.h"
25 #include <VM/VM_SSEFunc.h>
26 #include <SYS/SYS_Floor.h>
27 #include <SYS/SYS_Inline.h>
28 #include <SYS/SYS_Math.h>
29 #include <SYS/SYS_Types.h>
30 
31 typedef enum {
32  // These splines all work on uniform keys.
33  UT_SPLINE_CONSTANT, // Constant between keys
34  UT_SPLINE_LINEAR, // Linear interpolation between keys
35  UT_SPLINE_CATMULL_ROM, // Catmull-Rom Cardinal Spline
36  UT_SPLINE_MONOTONECUBIC, // Monotone Cubic Hermite Spline
37 
38  // This interpolation works a little differently. It takes a set of scalar
39  // values, and "fits" the parametric coordinate into the keys. That is, it
40  // performs the binary search to find where the parametric coordinate maps,
41  // then, it performs a linear interpolation between the two nearest key
42  // values to figure out what the coordinate should be.
44 
45  UT_SPLINE_BEZIER, // Bezier Curve
46  UT_SPLINE_BSPLINE, // B-Spline
47  UT_SPLINE_HERMITE, // Hermite Spline
49 
50 UT_API extern const char * UTnameFromSplineBasis(UT_SPLINE_BASIS basis);
52 
53 
55 {
56 public:
57  /// Evaluates an open spline at a given location.
58  /// The given CV list must have 4 elements in it!
59  /// The cvs should be for the current segment, and t in [0, 1]
60  template <typename T>
61  static T evalOpen(const T *cvs, float t)
62  {
63  UT_Matrix4 weightmatrix = getOpenWeights();
64  float t2 = t*t;
65  float t3 = t2*t;
66  UT_Vector4 powers(1, t, t2, t3);
67 
68  UT_Vector4 weights = colVecMult(weightmatrix, powers);
69 
70  T value;
71 
72  value = cvs[0] * weights[0];
73  value += cvs[1] * weights[1];
74  value += cvs[2] * weights[2];
75  value += cvs[3] * weights[3];
76 
77  return value;
78  }
79 
80  /// Evaluates a range of t values in uniform increasing manner.
81  /// The cvs list should have 3 + nseg entries.
82  template <typename T>
83  static void evalRangeOpen(T *results, const T *cvs, float start_t, float step_t, int len_t, int nseg)
84  {
85  int curseg;
86  curseg = SYSfastFloor(start_t);
87  curseg = SYSclamp(curseg, 0, nseg-1);
88  float t = start_t - curseg;
89 
90  for (int i = 0; i < len_t; i++)
91  {
92  results[i] = evalOpen(&cvs[curseg], t);
93  t += step_t;
94  if (t > 1)
95  {
96  while (curseg < nseg-1)
97  {
98  curseg++;
99  t -= 1;
100  if (t <= 1)
101  break;
102  }
103  }
104  }
105  }
106 
107  /// Evaluates a closed spline at the given location. The given
108  /// cv list must have 4 elements. Whether we are end interpolated or
109  /// not depends on which segment this represents. The cvs list
110  /// should be the cvs for the current segment and t in [0, 1]
111  template <typename T>
112  static T evalClosed(const T *cvs, float t, int seg, int nseg, bool deriv = false)
113  {
114  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg, deriv);
115  float t2 = t*t;
116  float t3 = t2*t;
117  UT_Vector4 powers(1, t, t2, t3);
118 
119  UT_Vector4 weights = colVecMult(weightmatrix, powers);
120 
121  T value;
122 
123  value = cvs[0] * weights[0];
124  value += cvs[1] * weights[1];
125  value += cvs[2] * weights[2];
126  value += cvs[3] * weights[3];
127 
128  return value;
129  }
130 
131  /// Evaluates a range of t values in uniform increasing manner.
132  /// The cvs list should have 3 + nseg entries.
133  template <typename T>
134  static void evalRangeClosed(T *results, const T *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv = false)
135  {
136  int curseg;
137  curseg = SYSfastFloor(start_t);
138  curseg = SYSclamp(curseg, 0, nseg-1);
139  float t = start_t - curseg;
140 
141  for (int i = 0; i < len_t; i++)
142  {
143  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
144  t += step_t;
145  if (t > 1)
146  {
147  while (curseg < nseg-1)
148  {
149  curseg++;
150  t -= 1;
151  if (t <= 1)
152  break;
153  }
154  }
155  }
156  }
157 
158  template <typename T>
159  static T evalSubDStart(const T *cvs, float t)
160  {
161  // First segment is (1 - t + (1/6)t^3)*P0 + (t - (1/3)*t^3)*P1 + ((1/6)t^3)*P2
162  const float onesixth = 0.16666666666666667f;
163  float onesixtht3 = onesixth*t*t*t;
164  float w0 = 1 - t + onesixtht3;
165  float w1 = t - 2*onesixtht3;
166  float w2 = onesixtht3;
167 
168  T value = w0*cvs[0];
169  value += w1*cvs[1];
170  value += w2*cvs[2];
171  return value;
172  }
173 
174  template <typename T>
175  static T evalSubDEnd(const T *cvs, float t)
176  {
177  // Reverse t relative to evalSubDStart
178  t = 1.0f-t;
179 
180  const float onesixth = 0.16666666666666667f;
181  float onesixtht3 = onesixth*t*t*t;
182  float w0 = 1 - t + onesixtht3;
183  float w1 = t - 2*onesixtht3;
184  float w2 = onesixtht3;
185 
186  // Also reverse point order relative to evalSubDStart
187  T value = w0*cvs[2];
188  value += w1*cvs[1];
189  value += w2*cvs[0];
190  return value;
191  }
192 
193  template <float (func)(const float *,float)>
195  UT_BoundingBox &box, const UT_Vector3 *cvs,
196  const UT_Vector3 &a, const UT_Vector3 &b, const UT_Vector3 &c,
197  float rootmin, float rootmax)
198  {
199  // If the value of the quadratic has equal signs at zero
200  // and one, AND the derivative has equal signs at zero and one,
201  // it can't have crossed zero between zero and one, so we
202  // can skip the root find in that case. The other rejection
203  // case of a negative b^2-4ac is already checked by
204  // UT_RootFinder, because it doesn't depend on the range
205  // limits.
206 
207  // a+b+c = value of quadratic at one
208  // (a+b+c)*c > 0 iff signs of values are equal
209  UT_Vector3 abc = a + b + c;
210  abc *= c;
211  // 2a+b = derivative of quadratic at one
212  // (2a+b)*b > 0 iff signs of derivatives are equal
213  UT_Vector3 b2a = a * 2.0F + b;
214  b2a *= b;
215 
216  for (int DIM = 0; DIM < 3; DIM++)
217  {
218  // No chance of crossing zero in case descirbed above
219  // NOTE: The abc == 0 case can be rejected, because we
220  // already did enlargeBounds on both p values.
221  // The abc > 0 && b2a == 0 case can be rejected,
222  // because the peak of the quadratic has the same
223  // sign as the rest, so never crosses zero.
224  if (abc[DIM] >= 0 && b2a[DIM] >= 0)
225  continue;
226 
227  float t1, t2;
228  int nroots = UT_RootFinder::quadratic(a[DIM], b[DIM], c[DIM], t1, t2);
229  if (nroots == 0)
230  continue;
231 
232  float fcvs[4];
233  fcvs[0] = cvs[0][DIM];
234  fcvs[1] = cvs[1][DIM];
235  fcvs[2] = cvs[2][DIM];
236  fcvs[3] = cvs[3][DIM];
237 
238  // Add any minima/maxima to the bounding box
239  if (t1 > rootmin && t1 < rootmax)
240  {
241  float v = func(fcvs, t1);
242  box.vals[DIM][0] = SYSmin(box.vals[DIM][0], v);
243  box.vals[DIM][1] = SYSmax(box.vals[DIM][1], v);
244  }
245  if (nroots == 2 && t2 > rootmin && t2 < rootmax)
246  {
247  float v = func(fcvs, t2);
248  box.vals[DIM][0] = SYSmin(box.vals[DIM][0], v);
249  box.vals[DIM][1] = SYSmax(box.vals[DIM][1], v);
250  }
251  }
252  }
253 
254  /// Enlarges box by any minima/maxima of the cubic curve defined by 4 cvs, that lie between rootmin and rootmax.
255  /// NOTE: This must be defined below so that it doesn't instantiate evalOpen before its specialization below.
256  static inline void enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax);
257 
258  /// Enlarges box by any minima/maxima of the cubic curve defined by 3 cvs, that lie between rootmin and rootmax.
259  /// The curve in this case is the start segment of a subdivision curve.
260  static inline void enlargeBoundingBoxSubDStart(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax);
261 
262  /// Enlarges box by any minima/maxima of the cubic curve defined by cvs, that lie between rootmin and rootmax.
263  /// The curve in this case is the end segment of a subdivision curve.
264  static inline void enlargeBoundingBoxSubDEnd(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax);
265 
266  /// Returns the weights for a power-basis evaluation of a segment.
267  /// The t values should be normalized inside the segment.
268  /// The format is (1, t, t^2, t^3), and colVecMult.
269  /// Assumes uniform knots.
271  {
272  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
273  4/6., 0/6., -6/6., 3/6.,
274  1/6., 3/6., 3/6., -3/6.,
275  0/6., 0/6., 0/6., 1/6. );
276  }
278  {
279  return UT_Matrix4( 1/6., 4/6., 1/6., 0/6.,
280  -3/6., 0/6., 3/6., 0/6.,
281  3/6., -6/6., 3/6., 0/6.,
282  -1/6., 3/6., -3/6., 1/6. );
283  }
284 
285  template<typename T>
286  static T evalMatrix(const UT_Matrix4 &basis, const T cvs[4], float t)
287  {
288  float t2 = t*t;
289  UT_Vector4 tpow(1.0f, t, t2, t2*t);
290 
291  UT_Vector4 coeff = colVecMult(basis, tpow);
292 
293  T val = cvs[0]*coeff[0] + cvs[1]*coeff[1] + cvs[2]*coeff[2] + cvs[3]*coeff[3];
294 
295  return val;
296  }
297 
298  /// This function is for evaluating a subdivision curve that is open.
299  /// For simplicitly, the parameter range is [0,1].
300  /// It's implemented in a way that maximizes stability
301  /// and readability, not necessarily performance.
302  template<typename T>
303  static T evalSubDCurve(const T *cvs, float t, int npts, bool deriv=false)
304  {
305  T p0;
306  T diff; // p1-p0
307  T c0; // Average of neighbours of p0, minus p0
308  T c1; // Average of neighbours of p1, minus p1
309 
310  // npts-1 segments, since npts points in whole curve
311  t *= (npts-1);
312 
313  int i = int(t);
314 
315  if (i < 0)
316  i = 0;
317  else if (i > npts-1)
318  i = npts-1;
319 
320  t -= i;
321  p0 = cvs[i];
322  diff = cvs[i+1]-cvs[i];
323 
324  if (i > 0)
325  c0 = 0.5*(cvs[i-1]+cvs[i+1]) - cvs[i];
326  else
327  c0 = 0;
328 
329  if (i < npts-1)
330  c1 = 0.5*(cvs[i]+cvs[i+2]) - cvs[i+1];
331  else
332  c1 = 0;
333 
334  float ti = 1-t;
335  if (!deriv)
336  {
337  float t3 = t*t*t/3;
338  float ti3 = ti*ti*ti/3;
339  // Order of addition should reduce roundoff in common cases.
340  return p0 + (diff*t + (c0*ti3 + c1*t3));
341  }
342  else
343  {
344  float t2 = t*t;
345  float ti2 = ti*ti;
346  // Order of addition should reduce roundoff in common cases.
347  return diff + (c1*t2 - c0*ti2);
348  }
349  }
350 
351  /// Basis for first segment of subd curve. Evaluation is:
352  /// [p[0] p[1] p[2] p[3]] * theSubDFirstBasis * [1 t t^2 t^3]^T
353  /// The last segment can be evaluated as: (NOTE the reversed order and 1-t)
354  /// [p[n-1] p[n-2] p[n-3] p[n-4]] * theSubDFirstBasis * [1 (1-t) (1-t)^2 (1-t)^3]^T
355  /// FYI: The last row is all zero, since it only depends on 3 points.
357 
358  /// Basis for derivative of first segment of subd curve. Evaluation is:
359  /// [p[0] p[1] p[2] p[3]] * theSubDFirstDerivBasis * [1 t t^2 t^3]^T
360  /// The last segment derivative can be evaluated as: (NOTE the reversed order and 1-t)
361  /// [p[n-1] p[n-2] p[n-3] p[n-4]] * theSubDFirstDerivBasis * [1 (1-t) (1-t)^2 t^3]^T
362  /// FYI: The last row is all zero, since it only depends on 3 points.
363  /// The last column is all zero, since the derivative has no cubic component.
365 
366  /// Basis for middle segment of subd curve or uniform, open, cubic NURBS.
367  /// Evaluation is:
368  /// [p[-1] p[0] p[1] p[2]] * theOpenBasis * [1 t t^2 t^3]^T
369  static const UT_Matrix4 theOpenBasis;
370 
371  /// Basis for derivative of middle segment of subd curve or uniform, open, cubic NURBS.
372  /// Evaluation is:
373  /// [p[-1] p[0] p[1] p[2]] * theOpenDerivBasis * [1 t t^2 t^3]^T
374  /// FYI: The last column is all zero, since the derivative has no cubic component.
376 
377  /// Basis for first segment of interpolating curve. Evaluation is:
378  /// [p[0] p[1] p[2] p[3]] * theInterpFirstBasis * [1 t t^2 t^3]^T
379  /// The last segment can be evaluated as: (NOTE the reversed order and 1-t)
380  /// [p[n-1] p[n-2] p[n-3] p[n-4]] * theInterpFirstBasis * [1 (1-t) (1-t)^2 (1-t)^3]^T
381  /// FYI: The last row is all zero, since it only depends on 3 points.
383 
384  /// Basis for middle segment of interpolating curve. Evaluation is:
385  /// [p[-1] p[0] p[1] p[2]] * theInterpBasis * [1 t t^2 t^3]^T
387 
388  /// Basis for Hermite cubic curve using two values (p) and two derivatives (v):
389  /// Evaluation is:
390  /// [p[0] p[1] v[0] v[1]] * theHermiteBasis * [1 t t^2 t^3]^T
392 
393  /// Basis for derivative of Hermite cubic curve using two values (p) and two derivatives (v):
394  /// Evaluation is:
395  /// [p[0] p[1] v[0] v[1]] * theHermiteDerivBasis * [1 t t^2 t^3]^T
396  /// FYI: The last column is all zero, since the derivative has no cubic component.
398 
399  /// Uniform knots with closed end conditions. seg is which segment
400  /// is being evaluates, nseg is the total. nseg should be
401  /// number of vertices minus three as we have cubics.
402  static UT_Matrix4 getClosedWeights(int seg, int nseg, bool deriv = false)
403  {
404  // these matrices come from $GEO/support/computespline.py
405  // which computes the power basis form of end-interpolated
406  // uniform bsplines.
407 
408  if (deriv == false)
409  {
410  if (nseg <= 1)
411  {
412  // Bezier.
413  return UT_Matrix4( 1, -3, 3, -1,
414  0, 3, -6, 3,
415  0, 0, 3, -3,
416  0, 0, 0, 1 );
417  }
418  else if (nseg == 2)
419  {
420  // 0, 0, 0, 1, 2, 2, 2,
421  if (seg == 0)
422  return UT_Matrix4( 1, -3, 3, -1,
423  0, 3, -4.5, 1.75,
424  0, 0, 1.5, -1,
425  0, 0, 0, 0.25 );
426  else
427  return UT_Matrix4( .25, -.75, .75, -0.25,
428  0.5, 0, -1.5, 1,
429  0.25, 0.75, 0.75, -1.75,
430  0, 0, 0, 1 );
431  }
432  else if (nseg == 3)
433  {
434  // 0, 0, 0, 1, 2, 3, 3, 3
435  if (seg == 0)
436  return UT_Matrix4( 1, -3, 3, -1,
437  0, 3, -4.5, 1.75,
438  0, 0, 1.5, -11/12.,
439  0, 0, 0, 1/6.);
440  else if (seg == 1)
441  return UT_Matrix4( .25, -.75, .75, -0.25,
442  7/12., 0.25, -1.25, 7/12.,
443  1/6., 0.5, 0.5, -7/12.,
444  0, 0, 0, 0.25 );
445  else
446  return UT_Matrix4( 1/6., -.5, .5, -1/6.,
447  7/12., -0.25, -1.25, 11/12.,
448  0.25, 0.75, 0.75, -1.75,
449  0, 0, 0, 1 );
450  }
451  else
452  {
453  // Either on an end, or in the middle
454  if (seg >= 2 && seg < nseg-2)
455  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
456  4/6., 0/6., -6/6., 3/6.,
457  1/6., 3/6., 3/6., -3/6.,
458  0/6., 0/6., 0/6., 1/6. );
459  else if (seg == 0)
460  return UT_Matrix4( 1, -3, 3, -1,
461  0, 3, -4.5, 1.75,
462  0, 0, 1.5, -11/12.,
463  0, 0, 0, 1/6. );
464  else if (seg == 1)
465  return UT_Matrix4( 0.25, -0.75, 0.75, -.25,
466  7/12., 0.25, -1.25, 7/12.,
467  1/6., 0.5, 0.5, -0.5,
468  0, 0, 0, 1/6. );
469  else if (seg == nseg-2)
470  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
471  2/3., 0, -1, 0.5,
472  1/6., 0.5, 0.5, -7/12.,
473  0, 0, 0, 0.25 );
474  else // if (seg == nseg-1)
475  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
476  7/12., -.25, -1.25, 11/12.,
477  0.25, 0.75, 0.75, -1.75,
478  0, 0, 0, 1 );
479  }
480  }
481  else
482  {
483  if (nseg <= 1)
484  {
485  // Bezier.
486  return UT_Matrix4( -3, 6, -3, 0,
487  3, -12, 9, 0,
488  0, 6, -9, 0,
489  0, 0, 3, 0 );
490  }
491  else if (nseg == 2)
492  {
493  // 0, 0, 0, 1, 2, 2, 2,
494  if (seg == 0)
495  return UT_Matrix4( -3, 6, -3, 0,
496  3, -9, 5.25, 0,
497  0, 3, -3, 0,
498  0, 0, 0.75, 9 );
499  else
500  return UT_Matrix4( -.75, 1.5, -0.75, 0,
501  0, -3, 3, 0,
502  0.75, 1.5, -5.25, 0,
503  0, 0, 3, 0 );
504  }
505  else if (nseg == 3)
506  {
507  // 0, 0, 0, 1, 2, 3, 3, 3
508  if (seg == 0)
509  return UT_Matrix4( -3, 6, -3, 0,
510  3, -9, 5.25, 0,
511  0, 3, -11/4., 0,
512  0, 0, .5, 0);
513  else if (seg == 1)
514  return UT_Matrix4( -.75, 1.5, -0.75, 0,
515  0.25, -2.5, 7/4., 0,
516  0.5, 1, -7/4., 0,
517  0, 0, 0.75, 0);
518  else
519  return UT_Matrix4( -.5, 1, -0.5, 0,
520  -0.25, -2.5, 11/4., 0,
521  0.75, 1.5, -5.25, 0,
522  0, 0, 3, 0);
523  }
524  else
525  {
526  // Either on an end, or in the middle
527  if (seg >= 2 && seg < nseg-2)
528  return UT_Matrix4( -3/6., 1.0, -0.5, 0,
529  0/6., -2.0, 1.5, 0,
530  3/6., 1.0, -1.5, 0,
531  0/6., 0, 0.5, 0);
532  else if (seg == 0)
533  return UT_Matrix4( -3, 6, -3, 0,
534  3, -9, 5.25, 0,
535  0, 3, -11/4., 0,
536  0, 0, 0.5, 0);
537  else if (seg == 1)
538  return UT_Matrix4( -0.75, 1.5, -.75, 0,
539  0.25, -2.5, 7/4., 0,
540  0.5, 1, -1.5, 0,
541  0, 0, 0.5, 0 );
542  else if (seg == nseg-2)
543  return UT_Matrix4( -3/6., 1, -0.5, 0,
544  0, -2, 1.5, 0,
545  0.5, 1, -7/4., 0,
546  0, 0, 0.75, 0 );
547  else // if (seg == nseg-1)
548  return UT_Matrix4(-3/6., 1, -.5, 0,
549  -.25, -2.5, 11/4., 0,
550  0.75, 1.5, -5.25, 0,
551  0, 0, 3, 0);
552  }
553  }
554  }
555  static UT_Matrix4 getClosedWeightsTranspose(int seg, int nseg, bool deriv = false)
556  {
557  if (deriv == false)
558  {
559  // these matrices come from $GEO/support/computespline.py
560  // which computes the power basis form of end-interpolated
561  // uniform bsplines.
562  if (nseg <= 1)
563  {
564  // Bezier.
565  return UT_Matrix4( 1, 0, 0, 0,
566  -3, 3, 0, 0,
567  3, -6, 3, 0,
568  -1, 3, -3, 1 );
569  }
570  else if (nseg == 2)
571  {
572  // 0, 0, 0, 1, 2, 2, 2,
573  if (seg == 0)
574  return UT_Matrix4( 1, 0, 0, 0,
575  -3, 3, 0, 0,
576  3, -4.5, 1.5, 0,
577  -1, 1.75, -1, 0.25 );
578  else
579  return UT_Matrix4( .25, .5, .25, 0,
580  -.75, 0, .75, 0,
581  0.75, -1.5, 0.75, 0,
582  -0.25, 1, -1.75, 1 );
583  }
584  else if (nseg == 3)
585  {
586  // 0, 0, 0, 1, 2, 3, 3, 3
587  if (seg == 0)
588  return UT_Matrix4( 1, 0, 0, 0,
589  -3, 3, 0, 0,
590  3,-4.5,1.5, 0,
591  -1,1.75,-11/12.,1/6. );
592  else if (seg == 1)
593  return UT_Matrix4( 0.25, 7/12., 1/6., 0,
594  -.75, 0.25, 0.5, 0,
595  .75,-1.25, 0.5, 0,
596  -.25,7/12.,-7/12.,0.25 );
597  else
598  return UT_Matrix4( 1/6., 7/12., 0.25, 0,
599  -.5, -0.25, 0.75, 0,
600  .5, -1.25, 0.75, 0,
601  -1/6.,11/12.,-1.75, 1 );
602 
603  }
604  else
605  {
606  // Either on an end, or in the middle
607  if (seg >= 2 && seg < nseg-2)
608  return UT_Matrix4( 1/6., 4/6., 1/6., 0/6.,
609  -3/6., 0/6., 3/6., 0/6.,
610  3/6., -6/6., 3/6., 0/6.,
611  -1/6., 3/6., -3/6., 1/6. );
612  else if (seg == 0)
613  return UT_Matrix4( 1, 0, 0, 0,
614  -3, 3, 0, 0,
615  3,-4.5,1.5, 0,
616  -1,1.75,-11/12., 1/6. );
617  else if (seg == 1)
618  return UT_Matrix4( 0.25, 7/12., 1/6., 0,
619  -0.75, 0.25, 0.5, 0,
620  0.75,-1.25, 0.5, 0,
621  -0.25,7/12., -0.5, 1/6. );
622  else if (seg == nseg-2)
623  return UT_Matrix4( 1/6., 2/3., 1/6., 0,
624  -3/6., 0, 0.5, 0,
625  3/6., -1, 0.5, 0,
626  -1/6., 0.5,-7/12.,0.25 );
627  else // if (seg == nseg-1)
628  return UT_Matrix4( 1/6., 7/12., 0.25, 0,
629  -3/6., -.25, 0.75, 0,
630  3/6.,-1.25, 0.75, 0,
631  -1/6.,11/12.,-1.75, 1 );
632  }
633  }
634  else
635  {
636  if (nseg <= 1)
637  {
638  // Bezier.
639  return UT_Matrix4(-3, 3, 0, 0,
640  6, -12, 6, 0,
641  -3, 9, -9, 3,
642  0, 0, 0, 0);
643  }
644  else if (nseg == 2)
645  {
646  // 0, 0, 0, 1, 2, 2, 2,
647  if (seg == 0)
648  return UT_Matrix4(-3, 3, 0, 0,
649  6, -9, 3, 0,
650  -3, 5.25, -3, 0.75,
651  0, 0, 0, 0);
652  else
653  return UT_Matrix4(-.75, 0, .75, 0,
654  1.5, -3, 1.5, 0,
655  -0.75, 3, -5.25, 3,
656  0, 0, 0, 0);
657  }
658  else if (nseg == 3)
659  {
660  // 0, 0, 0, 1, 2, 3, 3, 3
661  if (seg == 0)
662  return UT_Matrix4(-3, 3, 0, 0,
663  6, -9, 3, 0,
664  -3, 5.25, -11/4., .5,
665  0, 0, 0, 0);
666  else if (seg == 1)
667  return UT_Matrix4(-.75, 0.25, 0.5, 0,
668  1.5,-2.5, 1, 0,
669  -.75,7/4.,-7/4.,0.75,
670  0, 0, 0, 0);
671  else
672  return UT_Matrix4(-.5, -0.25, 0.75, 0,
673  1, -2.5, 1.5, 0,
674  -.5, 11/4., -5.25, 3,
675  0, 0, 0, 0);
676 
677  }
678  else
679  {
680  // Either on an end, or in the middle
681  if (seg >= 2 && seg < nseg-2)
682  return UT_Matrix4(-3/6., 0/6., 3/6., 0/6.,
683  1, -2, 1, 0,
684  -0.5, 1.5, -1.5, 0.5,
685  0, 0, 0, 0);
686  else if (seg == 0)
687  return UT_Matrix4(-3, 3, 0, 0,
688  6, -9, 3, 0,
689  -3, 5.25, -11/4., .5,
690  0, 0, 0, 0);
691  else if (seg == 1)
692  return UT_Matrix4(-0.75, 0.25, 0.5, 0,
693  1.5, -2.5, 1, 0,
694  -0.75, 7/4., -1.5, .5,
695  0, 0, 0, 0);
696 
697  else if (seg == nseg-2)
698  return UT_Matrix4(-3/6., 0, 0.5, 0,
699  1, -2, 1, 0,
700  -0.5, 1.5, -7/4., 0.75,
701  0, 0, 0, 0);
702  else // if (seg == nseg-1)
703  return UT_Matrix4(-3/6., -.25, 0.75, 0,
704  1, -2.5, 1.5, 0,
705  -.5, 11/4.,-5.25, 3,
706  0, 0, 0, 0);
707  }
708  }
709  }
710 };
711 
712 
713 /// The Linear & Catmull-Rom splines expect a parametric coordinate for
714 /// evaluation between 0 and 1. The Catmull-Rom spline requires additional
715 /// key values at the beginning and end of the spline to evaluate the slopes
716 /// of the Hermite spline properly.
717 ///
718 /// The LinearSolve only works on scalar values. It will compute the
719 /// parametric coordinate associated with the value passed in. This can be
720 /// used to simulate non-uniform keys on the spline.
722 public:
723  UT_Spline();
724  ~UT_Spline();
725 
726  /// Return the amount of memory owned by this UT_Spline in bytes
727  int64 getMemoryUsage(bool inclusive) const;
728 
729  /// Query the basis or knot length of the spline
730  UT_SPLINE_BASIS getGlobalBasis() const { return myGlobalBasis; }
731  int getVectorSize() const { return myVectorSize; }
732  int getKnotLength() const { return myKnotLength; }
733  fpreal64 getTension() const { return myTension; }
734 
736  { myGlobalBasis = b; }
737 
738  /// Sets the spline's global basis to b, and also set the basis
739  /// of every knot to b.
740  void setAllBases(UT_SPLINE_BASIS b);
741 
742  /// Construction of the spline object. All values are initialized to 0.
743  /// Warning, calling setSize() will clear all existing values.
744  void setSize(int nkeys, int vector_size);
745  /// Cubic splines may have a "tension". The tension defaults to 0.5 which
746  /// results in Catmull-Rom splines.
747  void setTension(fpreal64 t);
748 
749  /// Once the spline has been constructed, the values need to be set.
750  /// It is possible to change values between evaluations.
751  // @{
752  void setValue(int key, const fpreal32 *value, int size);
753  void setValue(int key, const fpreal64 *value, int size);
754  // @}
755 
756  /// Set the basis for the given key index.
757  /// This will also set the global basis.
758  void setBasis(int key, UT_SPLINE_BASIS b);
759 
760  /// Evaluate the spline using the global basis.
761  /// When interp_space is not UT_RGB, then values are treated as UT_RGBA
762  /// and converted into the desired color space before being interpolated.
763  /// The result is always returned as UT_RGBA.
764  ///
765  /// If order == 0, the value is returned.
766  /// If order == 1, the derivative with respect to t is returned.
767  // @{
768  bool evaluate(fpreal t, fpreal32 *result, int size,
769  UT_ColorType interp_space, int order = 0) const;
770  bool evaluate(fpreal t, fpreal64 *result, int size,
771  UT_ColorType interp_space, int order = 0) const;
772  // @}
773 
774  /// Evaluate the spline using multiple basis types depending on t.
775  /// Also unlike evaluate(), evaluateMulti() doesn't require extra keys for
776  /// Catmull-Rom on the ends. It always evaluates using a 0 slope at the
777  /// ends.
778  ///
779  /// If order == 0, the value is returned.
780  /// If order == 1, the derivative with respect to t is returned.
781  // @{
782  bool evaluateMulti(fpreal t, fpreal32 *result, int n,
783  UT_ColorType interp_space,
784  int knot_segment_hint = -1,
785  int order = 0) const;
786  bool evaluateMulti(fpreal t, fpreal64 *result, int n,
787  UT_ColorType interp_space,
788  int knot_segment_hint = -1,
789  int order = 0) const;
790  // @}
791 
792  /// Return _monotone_ cubic Hermite slopes at the current knot given the
793  /// previous and next knots.
794  // @{
795  /// Fritsch-Carlson (local) method.
796  /// It gives relatively good looking C1 curves but might not give a C2
797  /// solution even if it exists.
798  /// 1. Fritsch, F.N., Carlson, R.E., Monotone piecewise cubic interpolant,
799  /// SIAM J. Numer. Anal., 17(1980), 238-246.
800  static fpreal64 getMonotoneSlopeFC(fpreal64 v_cur, fpreal64 t_cur,
801  fpreal64 v_prev, fpreal64 t_prev,
802  fpreal64 v_next, fpreal64 t_next);
803  /// Paul Kupan's method.
804  /// Similar to Fritsch-Carlson method except it gives more visually
805  /// pleasing results when the intervals next to the knot are uneven.
806  /// 2. Kupan, P.A., Monotone Interpolant Built with Slopes Obtained by
807  /// Linear Combination, Studia Universitatis Babes-Bolyai Mathematica,
808  /// 53(2008), 59-66.
809  static fpreal64 getMonotoneSlopePA(fpreal64 v_cur, fpreal64 t_cur,
810  fpreal64 v_prev, fpreal64 t_prev,
811  fpreal64 v_next, fpreal64 t_next);
812  // @}
813 
814  /// Given the position within the two knots and the first knot index
815  /// number, normalize the position from knot-length domain to unit domain.
816  fpreal normalizeParameter(fpreal parm, fpreal t0, fpreal t1,
817  int t0_index) const;
818 
819  // Given parm in [0,1] interval in which CVs are at 'knots' values,
820  // find the parameter t in [0,1] interval in which CVs are evenly spaced.
821  // The returned reparameterized value is calculated in such a way that it
822  // yields smooth curve at CVs (unlike normalizeation that uses linear
823  // interpolation that yields tangent disconinuity at CVs).
824  // If 'parm_knot_segment' is given, this function sets it to the index
825  // of the knot segment into which the parm falls; it can be used as a hint
826  // to the evaluateMulti() method, especially for b-splines that may remap
827  // parameters to an adjacent segment (for continuity at end control points).
828  // If 'order' == 1, the derivative of t with respect to parm is returned.
829  fpreal solveSpline(fpreal parm,
830  const UT_Array<fpreal> &knots,
831  int *parm_knot_segment = nullptr,
832  int order = 0) const;
833 
834  /// Given the keys surrounding a channel segment, evaluate it as cubic
835  /// hermite spline. This function assumes dt > 0.
836  /// kt is [0,1] which maps over the dt time interval.
837  /// The `order`-th derivative with respect to kt is returned.
838  /// In particular, for order == 0 (default), the value is returned.
839  template <typename T>
840  static T evalCubic(T kt, T dt, T iv, T im, T ov, T om, int order = 0)
841  {
842  T x0 = iv;
843  T x1 = im*dt;
844  T x3 = 2*(iv - ov) + (im + om)*dt;
845  T x2 = ov - iv - im*dt - x3;
846  switch (order)
847  {
848  case 0:
849  return x0 + kt*(x1 + kt*(x2 + kt*x3));
850  case 1:
851  return x1 + kt*(2*x2 + kt*3*x3);
852  case 2:
853  return 2*x2 + kt*6*x3;
854  case 3:
855  return 6*x3;
856  default:
857  return 0;
858  }
859  }
860 
861  /// Given the keys surrounding a channel segment, evaluate it as
862  /// a quintic hermite spline. See evalCubic above.
863  template <typename T>
864  static T evalQuintic(T kt, T dt, T iv, T im, T ia, T ov, T om, T oa, int order = 0)
865  {
866  fpreal x0 = iv;
867  fpreal x1 = im * dt;
868  fpreal x2 = (fpreal).5 * ia * dt;
869  fpreal b0 = ov - (x0 + x1 + x2);
870  fpreal b1 = om * dt - (x1 + ia * dt);
871  fpreal b2 = oa * dt - ia * dt;
872  fpreal x5 = (fpreal).5 * (b2 - 6 * b1 + 12 * b0);
873  fpreal x4 = (fpreal).25 * (b2 - 2 * b1 - 10 * x5);
874  fpreal x3 = b0 - x4 - x5;
875 
876  switch (order)
877  {
878  case 0:
879  return x0 + kt*(x1 + kt*(x2 + kt*(x3 + kt*(x4 + kt*x5))));
880  case 1:
881  return x1 + kt*(2*x2 + kt*(3*x3 + kt*(4*x4 + kt*5*x5)));
882  case 2:
883  return 2*x2 + kt*(6*x3 + kt*(12*x4 + kt*20*x5));
884  case 3:
885  return 6*x3 + kt*(24*x4 + kt*60*x5);
886  case 4:
887  return 24*x4 + kt*120*x5;
888  case 5:
889  return 120*x5;
890  default:
891  return 0;
892  }
893  }
894 
895 private:
896  UT_Spline(const UT_Spline &copy); // not implemented yet
897  UT_Spline &operator =(const UT_Spline &copy); // not implemented yet
898 
899 private:
900  void grow(int size);
901  int getInterp(fpreal t, int knots[], fpreal weights[]);
902 
903  int64 getSizeOfValues() const
904  { return myKnotLength*myVectorSize*sizeof(fpreal64); }
905  int64 getSizeOfBases() const
906  { return myKnotLength*sizeof(UT_SPLINE_BASIS); }
907 
908  template <typename T>
909  inline void combineKeys(T *result, int vector_size,
910  fpreal64 weights[], int indices[],
911  int num_indices,
912  UT_ColorType interp_space) const;
913 
914  template <typename T>
915  inline void evalMonotoneCubic(fpreal t, T *result, int n,
916  int order, UT_ColorType interp_space,
917  bool do_multi) const;
918 
919  template <typename T>
920  inline void setValueInternal(int key, const T *value, int size);
921 
922  template <typename T>
923  inline bool evaluateInternal(fpreal t, T *result, int size,
924  UT_ColorType interp_space, int order) const;
925 
926  template <typename T>
927  inline bool evaluateMultiInternal(fpreal t, T *result, int n,
928  UT_ColorType interp_space,
929  int knot_segment_hint, int order) const;
930 
931  fpreal64 *myValues;
932  UT_SPLINE_BASIS *myBases;
933  fpreal64 myTension;
934  int myVectorSize;
935  int myKnotLength;
936  UT_SPLINE_BASIS myGlobalBasis;
937 };
938 
939 #include <VM/VM_SIMD.h>
940 
941 template <>
944 {
945 #if defined(CPU_HAS_SIMD_INSTR)
946  v4uf row1(1/6., 4/6., 1/6., 0/6.);
947  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
948  v4uf row3(3/6., -6/6., 3/6., 0/6.);
949  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
950 
951  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
952  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
953  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
954  v4uf vt(t);
955  v4uf vt2 = vt*vt;
956  v4uf vt3 = vt2*vt;
957 
958  v4uf weights;
959 
960  weights = row1;
961  weights += row2 * vt;
962  weights += row3 * vt2;
963  weights += row4 * vt3;
964 
965  vcvsx *= weights;
966  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
967  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
968  vcvsy *= weights;
969  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
970  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
971  vcvsz *= weights;
972  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
973  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
974 
975  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
976 #else
977  UT_Matrix4 weightmatrix = getOpenWeights();
978  float t2 = t*t;
979  float t3 = t2*t;
980  UT_Vector4 powers(1, t, t2, t3);
981 
982  UT_Vector4 weights = colVecMult(weightmatrix, powers);
983 
985 
986  value = cvs[0] * weights[0];
987  value += cvs[1] * weights[1];
988  value += cvs[2] * weights[2];
989  value += cvs[3] * weights[3];
990 
991  return value;
992 #endif
993 }
994 
995 template <>
996 SYS_FORCE_INLINE float
997 UT_SplineCubic::evalOpen(const float *cvs, float t)
998 {
999 #if defined(CPU_HAS_SIMD_INSTR)
1000  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1001  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1002  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1003  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1004 
1005  v4uf vcvs(cvs);
1006  v4uf vt(t);
1007  v4uf vt2 = vt*vt;
1008  v4uf vt3 = vt2*vt;
1009 
1010  v4uf weights;
1011 
1012  weights = row1;
1013  weights += row2 * vt;
1014  weights += row3 * vt2;
1015  weights += row4 * vt3;
1016 
1017  vcvs *= weights;
1018  vcvs += vcvs.swizzle<1, 1, 3, 3>();
1019  vcvs += vcvs.swizzle<2, 2, 2, 2>();
1020 
1021  return vcvs[0];
1022 #else
1023  UT_Matrix4 weightmatrix = getOpenWeights();
1024  float t2 = t*t;
1025  float t3 = t2*t;
1026  UT_Vector4 powers(1, t, t2, t3);
1027 
1028  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1029 
1030  float value;
1031 
1032  value = cvs[0] * weights[0];
1033  value += cvs[1] * weights[1];
1034  value += cvs[2] * weights[2];
1035  value += cvs[3] * weights[3];
1036 
1037  return value;
1038 #endif
1039 }
1040 
1041 template <>
1042 inline void
1043 UT_SplineCubic::evalRangeOpen(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg)
1044 {
1045  int curseg;
1046  curseg = SYSfastFloor(start_t);
1047  curseg = SYSclamp(curseg, 0, nseg-1);
1048  float t = start_t - curseg;
1049 
1050 #if defined(CPU_HAS_SIMD_INSTR)
1051  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1052  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1053  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1054  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1055 
1056  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1057  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1058  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1059 
1060  for (int i = 0; i < len_t; i++)
1061  {
1062  {
1063  v4uf weights;
1064  float t2 = t*t;
1065  float t3 = t2*t;
1066 
1067  weights = row1;
1068  weights += row2 * t;
1069  weights += row3 * t2;
1070  weights += row4 * t3;
1071 
1072  v4uf vx = vcvsx * weights;
1073  vx += vx.swizzle<1, 1, 3, 3>();
1074  vx += vx.swizzle<2, 2, 2, 2>();
1075  v4uf vy = vcvsy * weights;
1076  vy += vy.swizzle<1, 1, 3, 3>();
1077  vy += vy.swizzle<2, 2, 2, 2>();
1078  v4uf vz = vcvsz * weights;
1079  vz += vz.swizzle<1, 1, 3, 3>();
1080  vz += vz.swizzle<2, 2, 2, 2>();
1081  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1082  }
1083 
1084  t += step_t;
1085  if (t > 1)
1086  {
1087  while (curseg < nseg-1)
1088  {
1089  curseg++;
1090  t -= 1;
1091  if (t <= 1)
1092  break;
1093  }
1094  if (i < len_t-1)
1095  {
1096  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1097  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1098  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1099  }
1100  }
1101  }
1102 #else
1103  for (int i = 0; i < len_t; i++)
1104  {
1105  results[i] = evalOpen(&cvs[curseg], t);
1106  t += step_t;
1107  if (t > 1)
1108  {
1109  while (curseg < nseg-1)
1110  {
1111  curseg++;
1112  t -= 1;
1113  if (t <= 1)
1114  break;
1115  }
1116  }
1117  }
1118 #endif
1119 }
1120 
1121 template <>
1122 inline void
1123 UT_SplineCubic::evalRangeOpen(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg)
1124 {
1125  int curseg;
1126  curseg = SYSfastFloor(start_t);
1127  curseg = SYSclamp(curseg, 0, nseg-1);
1128  float t = start_t - curseg;
1129 
1130 #if defined(CPU_HAS_SIMD_INSTR)
1131  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1132  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1133  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1134  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1135 
1136  v4uf vcvs(&cvs[curseg]);
1137 
1138  for (int i = 0; i < len_t; i++)
1139  {
1140  {
1141  v4uf weights;
1142  float t2 = t*t;
1143  float t3 = t2*t;
1144 
1145  weights = row1;
1146  weights += row2 * t;
1147  weights += row3 * t2;
1148  weights += row4 * t3;
1149 
1150  v4uf v = vcvs * weights;
1151  v += v.swizzle<1, 1, 3, 3>();
1152  v += v.swizzle<2, 2, 2, 2>();
1153  results[i] = v[0];
1154  }
1155 
1156  t += step_t;
1157  if (t > 1)
1158  {
1159  while (curseg < nseg-1)
1160  {
1161  curseg++;
1162  t -= 1;
1163  if (t <= 1)
1164  break;
1165  }
1166  if (i < len_t-1)
1167  {
1168  vcvs = v4uf(&cvs[curseg]);
1169  }
1170  }
1171  }
1172 #else
1173  for (int i = 0 ; i < len_t; i++)
1174  {
1175  results[i] = evalOpen(&cvs[curseg], t);
1176  t += step_t;
1177  if (t > 1)
1178  {
1179  while (curseg < nseg-1)
1180  {
1181  curseg++;
1182  t -= 1;
1183  if (t <= 1)
1184  break;
1185  }
1186  }
1187  }
1188 #endif
1189 }
1190 
1191 template <>
1193 UT_SplineCubic::evalClosed(const UT_Vector3 *cvs, float t, int seg, int nseg, bool deriv)
1194 {
1195 #if defined(CPU_HAS_SIMD_INSTR)
1196  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1197 
1198  v4uf row1(weightmatrix.data());
1199  v4uf row2(weightmatrix.data()+4);
1200  v4uf row3(weightmatrix.data()+8);
1201  v4uf row4(weightmatrix.data()+12);
1202 
1203  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
1204  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
1205  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
1206  v4uf vt(t);
1207  v4uf vt2 = vt*vt;
1208  v4uf vt3 = vt2*vt;
1209 
1210  v4uf weights;
1211 
1212  weights = row1;
1213  weights += row2 * vt;
1214  weights += row3 * vt2;
1215  weights += row4 * vt3;
1216 
1217  vcvsx *= weights;
1218  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
1219  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
1220  vcvsy *= weights;
1221  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
1222  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
1223  vcvsz *= weights;
1224  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
1225  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
1226 
1227  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
1228 #else
1229  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1230  float t2 = t*t;
1231  float t3 = t2*t;
1232  UT_Vector4 powers(1, t, t2, t3);
1233 
1234  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1235 
1236  UT_Vector3 value;
1237 
1238  value = cvs[0] * weights[0];
1239  value += cvs[1] * weights[1];
1240  value += cvs[2] * weights[2];
1241  value += cvs[3] * weights[3];
1242 
1243  return value;
1244 #endif
1245 }
1246 
1247 template <>
1248 SYS_FORCE_INLINE float
1249 UT_SplineCubic::evalClosed(const float *cvs, float t, int seg, int nseg, bool deriv)
1250 {
1251 #if defined(CPU_HAS_SIMD_INSTR)
1252  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1253 
1254  v4uf row1(weightmatrix.data());
1255  v4uf row2(weightmatrix.data()+4);
1256  v4uf row3(weightmatrix.data()+8);
1257  v4uf row4(weightmatrix.data()+12);
1258 
1259  v4uf vcvs(cvs);
1260  float t2 = t*t;
1261  float t3 = t2*t;
1262 
1263  v4uf weights;
1264 
1265  weights = row1;
1266  weights += row2 * t;
1267  weights += row3 * t2;
1268  weights += row4 * t3;
1269 
1270  vcvs *= weights;
1271 
1272  vcvs += vcvs.swizzle<1, 1, 3, 3>();
1273  vcvs += vcvs.swizzle<2, 2, 2, 2>();
1274 
1275  return vcvs[0];
1276 #else
1277  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1278  float t2 = t*t;
1279  float t3 = t2*t;
1280  UT_Vector4 powers(1, t, t2, t3);
1281 
1282  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1283 
1284  float value;
1285 
1286  value = cvs[0] * weights[0];
1287  value += cvs[1] * weights[1];
1288  value += cvs[2] * weights[2];
1289  value += cvs[3] * weights[3];
1290 
1291  return value;
1292 #endif
1293 }
1294 
1295 template <>
1296 inline void
1297 UT_SplineCubic::evalRangeClosed(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1298 {
1299  int curseg;
1300  curseg = SYSfastFloor(start_t);
1301  curseg = SYSclamp(curseg, 0, nseg-1);
1302  float t = start_t - curseg;
1303 
1304 #if defined(CPU_HAS_SIMD_INSTR)
1305  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1306 
1307  v4uf row1(weightmatrix.data());
1308  v4uf row2(weightmatrix.data()+4);
1309  v4uf row3(weightmatrix.data()+8);
1310  v4uf row4(weightmatrix.data()+12);
1311 
1312  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1313  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1314  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1315 
1316  for (int i = 0; i < len_t; i++)
1317  {
1318  {
1319  v4uf weights;
1320  float t2 = t*t;
1321  float t3 = t2*t;
1322 
1323  weights = row1;
1324  weights += row2 * t;
1325  weights += row3 * t2;
1326  weights += row4 * t3;
1327 
1328  v4uf vx = vcvsx * weights;
1329  vx += vx.swizzle<1, 1, 3, 3>();
1330  vx += vx.swizzle<2, 2, 2, 2>();
1331  v4uf vy = vcvsy * weights;
1332  vy += vy.swizzle<1, 1, 3, 3>();
1333  vy += vy.swizzle<2, 2, 2, 2>();
1334  v4uf vz = vcvsz * weights;
1335  vz += vz.swizzle<1, 1, 3, 3>();
1336  vz += vz.swizzle<2, 2, 2, 2>();
1337  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1338  }
1339 
1340  t += step_t;
1341  if (t > 1)
1342  {
1343  while (curseg < nseg-1)
1344  {
1345  curseg++;
1346  t -= 1;
1347  if (t <= 1)
1348  break;
1349  }
1350  if (i < len_t-1)
1351  {
1352  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1353 
1354  row1 = v4uf(weightmatrix.data());
1355  row2 = v4uf(weightmatrix.data()+4);
1356  row3 = v4uf(weightmatrix.data()+8);
1357  row4 = v4uf(weightmatrix.data()+12);
1358 
1359  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1360  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1361  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1362  }
1363  }
1364  }
1365 #else
1366  for (int i = 0; i < len_t; i++)
1367  {
1368  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1369  t += step_t;
1370  if (t > 1)
1371  {
1372  while (curseg < nseg-1)
1373  {
1374  curseg++;
1375  t -= 1;
1376  if (t <= 1)
1377  break;
1378  }
1379  }
1380  }
1381 #endif
1382 }
1383 
1384 template <>
1385 inline void
1386 UT_SplineCubic::evalRangeClosed(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1387 {
1388  int curseg;
1389  curseg = SYSfastFloor(start_t);
1390  curseg = SYSclamp(curseg, 0, nseg-1);
1391  float t = start_t - curseg;
1392 
1393 #if defined(CPU_HAS_SIMD_INSTR)
1394  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1395 
1396  v4uf row1(weightmatrix.data());
1397  v4uf row2(weightmatrix.data()+4);
1398  v4uf row3(weightmatrix.data()+8);
1399  v4uf row4(weightmatrix.data()+12);
1400 
1401  v4uf vcvs(&cvs[curseg]);
1402 
1403  for (int i = 0; i < len_t; i++)
1404  {
1405  {
1406  v4uf weights;
1407  float t2 = t*t;
1408  float t3 = t2*t;
1409 
1410  weights = row1;
1411  weights += row2 * t;
1412  weights += row3 * t2;
1413  weights += row4 * t3;
1414 
1415  v4uf v = vcvs * weights;
1416  v += v.swizzle<1, 1, 3, 3>();
1417  v += v.swizzle<2, 2, 2, 2>();
1418  results[i] = v[0];
1419  }
1420 
1421  t += step_t;
1422  if (t > 1)
1423  {
1424  while (curseg < nseg-1)
1425  {
1426  curseg++;
1427  t -= 1;
1428  if (t <= 1)
1429  break;
1430  }
1431  if (i < len_t-1)
1432  {
1433  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1434 
1435  row1 = v4uf(weightmatrix.data());
1436  row2 = v4uf(weightmatrix.data()+4);
1437  row3 = v4uf(weightmatrix.data()+8);
1438  row4 = v4uf(weightmatrix.data()+12);
1439 
1440  vcvs = v4uf(&cvs[curseg]);
1441  }
1442  }
1443  }
1444 #else
1445  for (int i = 0 ; i < len_t; i++)
1446  {
1447  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1448  t += step_t;
1449  if (t > 1)
1450  {
1451  while (curseg < nseg-1)
1452  {
1453  curseg++;
1454  t -= 1;
1455  if (t <= 1)
1456  break;
1457  }
1458  }
1459  }
1460 #endif
1461 }
1462 
1463 void
1464 UT_SplineCubic::enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1465 {
1466  // We need to find any minimum or maximum in each dimension
1467  // to enlarge the bounding box.
1468  // To do this, for each, dimension, we take the derivative
1469  // of the cubic, leaving a quadratic, and find the zeros of it.
1470  // The quadratic is such that its ith derivatives at zero are
1471  // the (i+1)th derivatives of the curve segment at zero.
1472  // a = (1/2) * 3rd derivative of curve segment at zero
1473  UT_Vector3 a = -cvs[0] + cvs[1] * 3.0F + cvs[2] * (-3.0F) + cvs[3];
1474  a *= 0.5F;
1475 
1476  // b = 2nd derivative of curve segment at zero
1477  // (this is equivalent to the 2nd difference)
1478  UT_Vector3 b = cvs[0] + cvs[1] * (-2.0F) + cvs[2];
1479  // c = 1st derivative of curve segment at zero
1480  // (this is equivalent to the central difference)
1481  UT_Vector3 c = cvs[2] - cvs[0];
1482  c *= 0.5F;
1483 
1484  enlargeBoundingBoxCommon<UT_SplineCubic::evalOpen<float> >(box, cvs, a, b, c, rootmin, rootmax);
1485 }
1486 
1487 void
1488 UT_SplineCubic::enlargeBoundingBoxSubDStart(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1489 {
1490  // We need to find any minimum or maximum in each dimension
1491  // to enlarge the bounding box.
1492  // To do this, for each, dimension, we take the derivative
1493  // of the cubic, leaving a quadratic, and find the zeros of it.
1494  // The quadratic is such that its ith derivatives at zero are
1495  // the (i+1)th derivatives of the curve segment at zero.
1496 
1497  // First segment is (1 - t + (1/6)t^3)*P0 + (t - (1/3)*t^3)*P1 + ((1/6)t^3)*P2
1498  // 1st derivative is (-1 + (1/2)t^2)*P0 + (1 - t^2)*P1 + ((1/2)t^2)*P2
1499  // 2nd derivative is (t)*P0 + (-2t)*P1 + (t)*P2
1500  // 3rd derivative is (1)*P0 + (-2)*P1 + (1)*P2
1501 
1502  // a = (1/2) * 3rd derivative of curve segment at zero
1503  UT_Vector3 a = cvs[0] - 2.0f*cvs[1] + cvs[2];
1504  a *= 0.5F;
1505 
1506  // b = 2nd derivative of curve segment at zero
1507  // (this is equivalent to the 2nd difference)
1508  UT_Vector3 b(0,0,0);
1509  // c = 1st derivative of curve segment at zero
1510  // (this is equivalent to the central difference)
1511  UT_Vector3 c = cvs[1] - cvs[0];
1512 
1513  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDStart<float> >(box, cvs, a, b, c, rootmin, rootmax);
1514 }
1515 
1516 void
1517 UT_SplineCubic::enlargeBoundingBoxSubDEnd(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1518 {
1519  // We need to find any minimum or maximum in each dimension
1520  // to enlarge the bounding box.
1521  // To do this, for each, dimension, we take the derivative
1522  // of the cubic, leaving a quadratic, and find the zeros of it.
1523  // The quadratic is such that its ith derivatives at zero are
1524  // the (i+1)th derivatives of the curve segment at zero.
1525 
1526  // First segment is ((1/6)(1-t)^3)*P0 + ((1-t) - (1/3)*(1-t)^3)*P1 + (1 - (1-t) + (1/6)(1-t)^3)*P2
1527  // 1st derivative is (-(1/2)(1-t)^2)*P0 + (-1 + (1-t)^2)*P1 + (1 - (1/2)(1-t)^2)*P2
1528  // 2nd derivative is (1-t)*P0 + (-2(1-t))*P1 + (1-t)*P2
1529  // 3rd derivative is (-1)*P0 + (2)*P1 + (-1)*P2
1530 
1531  // a = (1/2) * 3rd derivative of curve segment at zero
1532  // b = 2nd derivative of curve segment at zero
1533  // (this is equivalent to the 2nd difference)
1534  UT_Vector3 b = cvs[0] - 2.0f*cvs[1] + cvs[2];
1535  UT_Vector3 a = -0.5f*b;
1536 
1537  // c = 1st derivative of curve segment at zero
1538  // (this is equivalent to the central difference)
1539  UT_Vector3 c = cvs[2] - cvs[0];
1540  c *= 0.5f;
1541 
1542  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDEnd<float> >(box, cvs, a, b, c, rootmin, rootmax);
1543 }
1544 
1545 #endif
static const UT_Matrix4 theHermiteDerivBasis
Definition: UT_Spline.h:397
static UT_Matrix4 getClosedWeightsTranspose(int seg, int nseg, bool deriv=false)
Definition: UT_Spline.h:555
#define SYSmax(a, b)
Definition: SYS_Math.h:1570
UT_SPLINE_BASIS getGlobalBasis() const
Query the basis or knot length of the spline.
Definition: UT_Spline.h:730
static T evalClosed(const T *cvs, float t, int seg, int nseg, bool deriv=false)
Definition: UT_Spline.h:112
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
GLsizei GLenum const void * indices
Definition: glcorearb.h:406
OIIO_UTIL_API bool copy(string_view from, string_view to, std::string &err)
const GLdouble * v
Definition: glcorearb.h:837
int getKnotLength() const
Definition: UT_Spline.h:732
GLsizei const GLfloat * value
Definition: glcorearb.h:824
static UT_Matrix4 getOpenWeightsTranspose()
Definition: UT_Spline.h:277
static void evalRangeClosed(T *results, const T *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv=false)
Definition: UT_Spline.h:134
UT_Vector3T< float > UT_Vector3
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
#define UT_API
Definition: UT_API.h:14
void setGlobalBasis(UT_SPLINE_BASIS b)
Definition: UT_Spline.h:735
GLint y
Definition: glcorearb.h:103
**But if you need a result
Definition: thread.h:613
static T evalSubDStart(const T *cvs, float t)
Definition: UT_Spline.h:159
__hostdev__ void setValue(uint32_t offset, bool v)
Definition: NanoVDB.h:5750
UT_SPLINE_BASIS
Definition: UT_Spline.h:31
float fpreal32
Definition: SYS_Types.h:200
static UT_Matrix4 getClosedWeights(int seg, int nseg, bool deriv=false)
Definition: UT_Spline.h:402
static void evalRangeOpen(T *results, const T *cvs, float start_t, float step_t, int len_t, int nseg)
Definition: UT_Spline.h:83
fpreal64 getTension() const
Definition: UT_Spline.h:733
GLdouble GLdouble x2
Definition: glad.h:2349
UT_Matrix4T< float > UT_Matrix4
static const UT_Matrix4 theInterpFirstBasis
Definition: UT_Spline.h:382
double fpreal64
Definition: SYS_Types.h:201
static UT_Matrix4 getOpenWeights()
Definition: UT_Spline.h:270
static const UT_Matrix4 theHermiteBasis
Definition: UT_Spline.h:391
GLdouble n
Definition: glcorearb.h:2008
GLfloat f
Definition: glcorearb.h:1926
int getVectorSize() const
Definition: UT_Spline.h:731
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
UT_Vector3T< T > SYSclamp(const UT_Vector3T< T > &v, const UT_Vector3T< T > &min, const UT_Vector3T< T > &max)
Definition: UT_Vector3.h:1057
Definition: VM_SIMD.h:188
GLdouble GLdouble GLint GLint order
Definition: glad.h:2676
long long int64
Definition: SYS_Types.h:116
static void enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
Definition: UT_Spline.h:1464
#define SYS_STATIC_FORCE_INLINE
Definition: SYS_Inline.h:48
static const UT_Matrix4 theOpenDerivBasis
Definition: UT_Spline.h:375
static const UT_Matrix4 theSubDFirstBasis
Definition: UT_Spline.h:356
GLuint const GLchar * name
Definition: glcorearb.h:786
static const UT_Matrix4 theOpenBasis
Definition: UT_Spline.h:369
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
GLint GLenum GLint x
Definition: glcorearb.h:409
static T evalCubic(T kt, T dt, T iv, T im, T ov, T om, int order=0)
Definition: UT_Spline.h:840
static T evalQuintic(T kt, T dt, T iv, T im, T ia, T ov, T om, T oa, int order=0)
Definition: UT_Spline.h:864
GLdouble t
Definition: glad.h:2397
static void enlargeBoundingBoxSubDStart(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
Definition: UT_Spline.h:1488
static T evalSubDCurve(const T *cvs, float t, int npts, bool deriv=false)
Definition: UT_Spline.h:303
GLsizeiptr size
Definition: glcorearb.h:664
GLenum func
Definition: glcorearb.h:783
UT_ColorType
Definition: UT_Color.h:24
const T * data() const
Return the raw matrix data.
Definition: UT_Matrix4.h:1163
static T evalMatrix(const UT_Matrix4 &basis, const T cvs[4], float t)
Definition: UT_Spline.h:286
static void enlargeBoundingBoxSubDEnd(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
Definition: UT_Spline.h:1517
fpreal64 fpreal
Definition: SYS_Types.h:277
static T evalOpen(const T *cvs, float t)
Definition: UT_Spline.h:61
LeafData & operator=(const LeafData &)=delete
UT_Vector3T< T > colVecMult(const UT_Matrix3T< S > &m, const UT_Vector3T< T > &v)
Definition: UT_Matrix3.h:1534
GLuint GLfloat * val
Definition: glcorearb.h:1608
static T evalSubDEnd(const T *cvs, float t)
Definition: UT_Spline.h:175
SYS_STATIC_FORCE_INLINE void enlargeBoundingBoxCommon(UT_BoundingBox &box, const UT_Vector3 *cvs, const UT_Vector3 &a, const UT_Vector3 &b, const UT_Vector3 &c, float rootmin, float rootmax)
Definition: UT_Spline.h:194
SYS_FORCE_INLINE v4uf swizzle() const
Definition: VM_SIMD.h:335
static const UT_Matrix4 theInterpBasis
Definition: UT_Spline.h:386
Definition: core.h:1131
UT_API UT_SPLINE_BASIS UTsplineBasisFromName(const char *name)
static int quadratic(T a, T b, T c, T &v0, T &v1)
UT_API const char * UTnameFromSplineBasis(UT_SPLINE_BASIS basis)
#define SYSmin(a, b)
Definition: SYS_Math.h:1571
static const UT_Matrix4 theSubDFirstDerivBasis
Definition: UT_Spline.h:364