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  /// Construction of the spline object. All values are initialized to 0.
739  /// Warning, calling setSize() will clear all existing values.
740  void setSize(int nkeys, int vector_size);
741  /// Cubic splines may have a "tension". The tension defaults to 0.5 which
742  /// results in Catmull-Rom splines.
743  void setTension(fpreal64 t);
744 
745  /// Once the spline has been constructed, the values need to be set.
746  /// It is possible to change values between evaluations.
747  // @{
748  void setValue(int key, const fpreal32 *value, int size);
749  void setValue(int key, const fpreal64 *value, int size);
750  // @}
751 
752  /// Set the basis for the given key index.
753  /// This will also set the global basis.
754  void setBasis(int key, UT_SPLINE_BASIS b);
755 
756  /// Evaluate the spline using the global basis.
757  /// When interp_space is not UT_RGB, then values are treated as UT_RGBA
758  /// and converted into the desired color space before being interpolated.
759  /// The result is always returned as UT_RGBA.
760  ///
761  /// If order == 0, the value is returned.
762  /// If order == 1, the derivative with respect to t is returned.
763  // @{
764  bool evaluate(fpreal t, fpreal32 *result, int size,
765  UT_ColorType interp_space, int order = 0) const;
766  bool evaluate(fpreal t, fpreal64 *result, int size,
767  UT_ColorType interp_space, int order = 0) const;
768  // @}
769 
770  /// Evaluate the spline using multiple basis types depending on t.
771  /// Also unlike evaluate(), evaluateMulti() doesn't require extra keys for
772  /// Catmull-Rom on the ends. It always evaluates using a 0 slope at the
773  /// ends.
774  ///
775  /// If order == 0, the value is returned.
776  /// If order == 1, the derivative with respect to t is returned.
777  // @{
778  bool evaluateMulti(fpreal t, fpreal32 *result, int n,
779  UT_ColorType interp_space,
780  int knot_segment_hint = -1,
781  int order = 0) const;
782  bool evaluateMulti(fpreal t, fpreal64 *result, int n,
783  UT_ColorType interp_space,
784  int knot_segment_hint = -1,
785  int order = 0) const;
786  // @}
787 
788  /// Return _monotone_ cubic Hermite slopes at the current knot given the
789  /// previous and next knots.
790  // @{
791  /// Fritsch-Carlson (local) method.
792  /// It gives relatively good looking C1 curves but might not give a C2
793  /// solution even if it exists.
794  /// 1. Fritsch, F.N., Carlson, R.E., Monotone piecewise cubic interpolant,
795  /// SIAM J. Numer. Anal., 17(1980), 238-246.
796  static fpreal64 getMonotoneSlopeFC(fpreal64 v_cur, fpreal64 t_cur,
797  fpreal64 v_prev, fpreal64 t_prev,
798  fpreal64 v_next, fpreal64 t_next);
799  /// Paul Kupan's method.
800  /// Similar to Fritsch-Carlson method except it gives more visually
801  /// pleasing results when the intervals next to the knot are uneven.
802  /// 2. Kupan, P.A., Monotone Interpolant Built with Slopes Obtained by
803  /// Linear Combination, Studia Universitatis Babes-Bolyai Mathematica,
804  /// 53(2008), 59-66.
805  static fpreal64 getMonotoneSlopePA(fpreal64 v_cur, fpreal64 t_cur,
806  fpreal64 v_prev, fpreal64 t_prev,
807  fpreal64 v_next, fpreal64 t_next);
808  // @}
809 
810  /// Given the position within the two knots and the first knot index
811  /// number, normalize the position from knot-length domain to unit domain.
812  fpreal normalizeParameter(fpreal parm, fpreal t0, fpreal t1,
813  int t0_index) const;
814 
815  // Given parm in [0,1] interval in which CVs are at 'knots' values,
816  // find the parameter t in [0,1] interval in which CVs are evenly spaced.
817  // The returned reparameterized value is calculated in such a way that it
818  // yields smooth curve at CVs (unlike normalizeation that uses linear
819  // interpolation that yields tangent disconinuity at CVs).
820  // If 'parm_knot_segment' is given, this function sets it to the index
821  // of the knot segment into which the parm falls; it can be used as a hint
822  // to the evaluateMulti() method, especially for b-splines that may remap
823  // parameters to an adjacent segment (for continuity at end control points).
824  // If 'order' == 1, the derivative of t with respect to parm is returned.
825  fpreal solveSpline(fpreal parm,
826  const UT_Array<fpreal> &knots,
827  int *parm_knot_segment = nullptr,
828  int order = 0) const;
829 
830  /// Given the keys surrounding a channel segment, evaluate it as cubic
831  /// hermite spline. This function assumes dt > 0.
832  /// kt is [0,1] which maps over the dt time interval.
833  /// The `order`-th derivative with respect to kt is returned.
834  /// In particular, for order == 0 (default), the value is returned.
835  template <typename T>
836  static T evalCubic(T kt, T dt, T iv, T im, T ov, T om, int order = 0)
837  {
838  T x0 = iv;
839  T x1 = im*dt;
840  T x3 = 2*(iv - ov) + (im + om)*dt;
841  T x2 = ov - iv - im*dt - x3;
842  switch (order)
843  {
844  case 0:
845  return x0 + kt*(x1 + kt*(x2 + kt*x3));
846  case 1:
847  return x1 + kt*(2*x2 + kt*3*x3);
848  case 2:
849  return 2*x2 + kt*6*x3;
850  case 3:
851  return 6*x3;
852  default:
853  return 0;
854  }
855  }
856 
857 private:
858  UT_Spline(const UT_Spline &copy); // not implemented yet
859  UT_Spline &operator =(const UT_Spline &copy); // not implemented yet
860 
861 private:
862  void grow(int size);
863  int getInterp(fpreal t, int knots[], fpreal weights[]);
864 
865  int64 getSizeOfValues() const
866  { return myKnotLength*myVectorSize*sizeof(fpreal64); }
867  int64 getSizeOfBases() const
868  { return myKnotLength*sizeof(UT_SPLINE_BASIS); }
869 
870  template <typename T>
871  inline void combineKeys(T *result, int vector_size,
872  fpreal64 weights[], int indices[],
873  int num_indices,
874  UT_ColorType interp_space) const;
875 
876  template <typename T>
877  inline void evalMonotoneCubic(fpreal t, T *result, int n,
878  int order, UT_ColorType interp_space,
879  bool do_multi) const;
880 
881  template <typename T>
882  inline void setValueInternal(int key, const T *value, int size);
883 
884  template <typename T>
885  inline bool evaluateInternal(fpreal t, T *result, int size,
886  UT_ColorType interp_space, int order) const;
887 
888  template <typename T>
889  inline bool evaluateMultiInternal(fpreal t, T *result, int n,
890  UT_ColorType interp_space,
891  int knot_segment_hint, int order) const;
892 
893  fpreal64 *myValues;
894  UT_SPLINE_BASIS *myBases;
895  fpreal64 myTension;
896  int myVectorSize;
897  int myKnotLength;
898  UT_SPLINE_BASIS myGlobalBasis;
899 };
900 
901 #include <VM/VM_SIMD.h>
902 
903 template <>
906 {
907 #if defined(CPU_HAS_SIMD_INSTR)
908  v4uf row1(1/6., 4/6., 1/6., 0/6.);
909  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
910  v4uf row3(3/6., -6/6., 3/6., 0/6.);
911  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
912 
913  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
914  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
915  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
916  v4uf vt(t);
917  v4uf vt2 = vt*vt;
918  v4uf vt3 = vt2*vt;
919 
920  v4uf weights;
921 
922  weights = row1;
923  weights += row2 * vt;
924  weights += row3 * vt2;
925  weights += row4 * vt3;
926 
927  vcvsx *= weights;
928  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
929  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
930  vcvsy *= weights;
931  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
932  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
933  vcvsz *= weights;
934  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
935  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
936 
937  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
938 #else
939  UT_Matrix4 weightmatrix = getOpenWeights();
940  float t2 = t*t;
941  float t3 = t2*t;
942  UT_Vector4 powers(1, t, t2, t3);
943 
944  UT_Vector4 weights = colVecMult(weightmatrix, powers);
945 
947 
948  value = cvs[0] * weights[0];
949  value += cvs[1] * weights[1];
950  value += cvs[2] * weights[2];
951  value += cvs[3] * weights[3];
952 
953  return value;
954 #endif
955 }
956 
957 template <>
958 SYS_FORCE_INLINE float
959 UT_SplineCubic::evalOpen(const float *cvs, float t)
960 {
961 #if defined(CPU_HAS_SIMD_INSTR)
962  v4uf row1(1/6., 4/6., 1/6., 0/6.);
963  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
964  v4uf row3(3/6., -6/6., 3/6., 0/6.);
965  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
966 
967  v4uf vcvs(cvs);
968  v4uf vt(t);
969  v4uf vt2 = vt*vt;
970  v4uf vt3 = vt2*vt;
971 
972  v4uf weights;
973 
974  weights = row1;
975  weights += row2 * vt;
976  weights += row3 * vt2;
977  weights += row4 * vt3;
978 
979  vcvs *= weights;
980  vcvs += vcvs.swizzle<1, 1, 3, 3>();
981  vcvs += vcvs.swizzle<2, 2, 2, 2>();
982 
983  return vcvs[0];
984 #else
985  UT_Matrix4 weightmatrix = getOpenWeights();
986  float t2 = t*t;
987  float t3 = t2*t;
988  UT_Vector4 powers(1, t, t2, t3);
989 
990  UT_Vector4 weights = colVecMult(weightmatrix, powers);
991 
992  float value;
993 
994  value = cvs[0] * weights[0];
995  value += cvs[1] * weights[1];
996  value += cvs[2] * weights[2];
997  value += cvs[3] * weights[3];
998 
999  return value;
1000 #endif
1001 }
1002 
1003 template <>
1004 inline void
1005 UT_SplineCubic::evalRangeOpen(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg)
1006 {
1007  int curseg;
1008  curseg = SYSfastFloor(start_t);
1009  curseg = SYSclamp(curseg, 0, nseg-1);
1010  float t = start_t - curseg;
1011 
1012 #if defined(CPU_HAS_SIMD_INSTR)
1013  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1014  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1015  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1016  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1017 
1018  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1019  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1020  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1021 
1022  for (int i = 0; i < len_t; i++)
1023  {
1024  {
1025  v4uf weights;
1026  float t2 = t*t;
1027  float t3 = t2*t;
1028 
1029  weights = row1;
1030  weights += row2 * t;
1031  weights += row3 * t2;
1032  weights += row4 * t3;
1033 
1034  v4uf vx = vcvsx * weights;
1035  vx += vx.swizzle<1, 1, 3, 3>();
1036  vx += vx.swizzle<2, 2, 2, 2>();
1037  v4uf vy = vcvsy * weights;
1038  vy += vy.swizzle<1, 1, 3, 3>();
1039  vy += vy.swizzle<2, 2, 2, 2>();
1040  v4uf vz = vcvsz * weights;
1041  vz += vz.swizzle<1, 1, 3, 3>();
1042  vz += vz.swizzle<2, 2, 2, 2>();
1043  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1044  }
1045 
1046  t += step_t;
1047  if (t > 1)
1048  {
1049  while (curseg < nseg-1)
1050  {
1051  curseg++;
1052  t -= 1;
1053  if (t <= 1)
1054  break;
1055  }
1056  if (i < len_t-1)
1057  {
1058  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1059  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1060  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1061  }
1062  }
1063  }
1064 #else
1065  for (int i = 0; i < len_t; i++)
1066  {
1067  results[i] = evalOpen(&cvs[curseg], t);
1068  t += step_t;
1069  if (t > 1)
1070  {
1071  while (curseg < nseg-1)
1072  {
1073  curseg++;
1074  t -= 1;
1075  if (t <= 1)
1076  break;
1077  }
1078  }
1079  }
1080 #endif
1081 }
1082 
1083 template <>
1084 inline void
1085 UT_SplineCubic::evalRangeOpen(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg)
1086 {
1087  int curseg;
1088  curseg = SYSfastFloor(start_t);
1089  curseg = SYSclamp(curseg, 0, nseg-1);
1090  float t = start_t - curseg;
1091 
1092 #if defined(CPU_HAS_SIMD_INSTR)
1093  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1094  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1095  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1096  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1097 
1098  v4uf vcvs(&cvs[curseg]);
1099 
1100  for (int i = 0; i < len_t; i++)
1101  {
1102  {
1103  v4uf weights;
1104  float t2 = t*t;
1105  float t3 = t2*t;
1106 
1107  weights = row1;
1108  weights += row2 * t;
1109  weights += row3 * t2;
1110  weights += row4 * t3;
1111 
1112  v4uf v = vcvs * weights;
1113  v += v.swizzle<1, 1, 3, 3>();
1114  v += v.swizzle<2, 2, 2, 2>();
1115  results[i] = v[0];
1116  }
1117 
1118  t += step_t;
1119  if (t > 1)
1120  {
1121  while (curseg < nseg-1)
1122  {
1123  curseg++;
1124  t -= 1;
1125  if (t <= 1)
1126  break;
1127  }
1128  if (i < len_t-1)
1129  {
1130  vcvs = v4uf(&cvs[curseg]);
1131  }
1132  }
1133  }
1134 #else
1135  for (int i = 0 ; i < len_t; i++)
1136  {
1137  results[i] = evalOpen(&cvs[curseg], t);
1138  t += step_t;
1139  if (t > 1)
1140  {
1141  while (curseg < nseg-1)
1142  {
1143  curseg++;
1144  t -= 1;
1145  if (t <= 1)
1146  break;
1147  }
1148  }
1149  }
1150 #endif
1151 }
1152 
1153 template <>
1155 UT_SplineCubic::evalClosed(const UT_Vector3 *cvs, float t, int seg, int nseg, bool deriv)
1156 {
1157 #if defined(CPU_HAS_SIMD_INSTR)
1158  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1159 
1160  v4uf row1(weightmatrix.data());
1161  v4uf row2(weightmatrix.data()+4);
1162  v4uf row3(weightmatrix.data()+8);
1163  v4uf row4(weightmatrix.data()+12);
1164 
1165  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
1166  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
1167  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
1168  v4uf vt(t);
1169  v4uf vt2 = vt*vt;
1170  v4uf vt3 = vt2*vt;
1171 
1172  v4uf weights;
1173 
1174  weights = row1;
1175  weights += row2 * vt;
1176  weights += row3 * vt2;
1177  weights += row4 * vt3;
1178 
1179  vcvsx *= weights;
1180  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
1181  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
1182  vcvsy *= weights;
1183  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
1184  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
1185  vcvsz *= weights;
1186  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
1187  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
1188 
1189  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
1190 #else
1191  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1192  float t2 = t*t;
1193  float t3 = t2*t;
1194  UT_Vector4 powers(1, t, t2, t3);
1195 
1196  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1197 
1198  UT_Vector3 value;
1199 
1200  value = cvs[0] * weights[0];
1201  value += cvs[1] * weights[1];
1202  value += cvs[2] * weights[2];
1203  value += cvs[3] * weights[3];
1204 
1205  return value;
1206 #endif
1207 }
1208 
1209 template <>
1210 SYS_FORCE_INLINE float
1211 UT_SplineCubic::evalClosed(const float *cvs, float t, int seg, int nseg, bool deriv)
1212 {
1213 #if defined(CPU_HAS_SIMD_INSTR)
1214  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1215 
1216  v4uf row1(weightmatrix.data());
1217  v4uf row2(weightmatrix.data()+4);
1218  v4uf row3(weightmatrix.data()+8);
1219  v4uf row4(weightmatrix.data()+12);
1220 
1221  v4uf vcvs(cvs);
1222  float t2 = t*t;
1223  float t3 = t2*t;
1224 
1225  v4uf weights;
1226 
1227  weights = row1;
1228  weights += row2 * t;
1229  weights += row3 * t2;
1230  weights += row4 * t3;
1231 
1232  vcvs *= weights;
1233 
1234  vcvs += vcvs.swizzle<1, 1, 3, 3>();
1235  vcvs += vcvs.swizzle<2, 2, 2, 2>();
1236 
1237  return vcvs[0];
1238 #else
1239  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1240  float t2 = t*t;
1241  float t3 = t2*t;
1242  UT_Vector4 powers(1, t, t2, t3);
1243 
1244  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1245 
1246  float value;
1247 
1248  value = cvs[0] * weights[0];
1249  value += cvs[1] * weights[1];
1250  value += cvs[2] * weights[2];
1251  value += cvs[3] * weights[3];
1252 
1253  return value;
1254 #endif
1255 }
1256 
1257 template <>
1258 inline void
1259 UT_SplineCubic::evalRangeClosed(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1260 {
1261  int curseg;
1262  curseg = SYSfastFloor(start_t);
1263  curseg = SYSclamp(curseg, 0, nseg-1);
1264  float t = start_t - curseg;
1265 
1266 #if defined(CPU_HAS_SIMD_INSTR)
1267  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1268 
1269  v4uf row1(weightmatrix.data());
1270  v4uf row2(weightmatrix.data()+4);
1271  v4uf row3(weightmatrix.data()+8);
1272  v4uf row4(weightmatrix.data()+12);
1273 
1274  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1275  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1276  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1277 
1278  for (int i = 0; i < len_t; i++)
1279  {
1280  {
1281  v4uf weights;
1282  float t2 = t*t;
1283  float t3 = t2*t;
1284 
1285  weights = row1;
1286  weights += row2 * t;
1287  weights += row3 * t2;
1288  weights += row4 * t3;
1289 
1290  v4uf vx = vcvsx * weights;
1291  vx += vx.swizzle<1, 1, 3, 3>();
1292  vx += vx.swizzle<2, 2, 2, 2>();
1293  v4uf vy = vcvsy * weights;
1294  vy += vy.swizzle<1, 1, 3, 3>();
1295  vy += vy.swizzle<2, 2, 2, 2>();
1296  v4uf vz = vcvsz * weights;
1297  vz += vz.swizzle<1, 1, 3, 3>();
1298  vz += vz.swizzle<2, 2, 2, 2>();
1299  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1300  }
1301 
1302  t += step_t;
1303  if (t > 1)
1304  {
1305  while (curseg < nseg-1)
1306  {
1307  curseg++;
1308  t -= 1;
1309  if (t <= 1)
1310  break;
1311  }
1312  if (i < len_t-1)
1313  {
1314  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1315 
1316  row1 = v4uf(weightmatrix.data());
1317  row2 = v4uf(weightmatrix.data()+4);
1318  row3 = v4uf(weightmatrix.data()+8);
1319  row4 = v4uf(weightmatrix.data()+12);
1320 
1321  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1322  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1323  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1324  }
1325  }
1326  }
1327 #else
1328  for (int i = 0; i < len_t; i++)
1329  {
1330  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1331  t += step_t;
1332  if (t > 1)
1333  {
1334  while (curseg < nseg-1)
1335  {
1336  curseg++;
1337  t -= 1;
1338  if (t <= 1)
1339  break;
1340  }
1341  }
1342  }
1343 #endif
1344 }
1345 
1346 template <>
1347 inline void
1348 UT_SplineCubic::evalRangeClosed(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1349 {
1350  int curseg;
1351  curseg = SYSfastFloor(start_t);
1352  curseg = SYSclamp(curseg, 0, nseg-1);
1353  float t = start_t - curseg;
1354 
1355 #if defined(CPU_HAS_SIMD_INSTR)
1356  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1357 
1358  v4uf row1(weightmatrix.data());
1359  v4uf row2(weightmatrix.data()+4);
1360  v4uf row3(weightmatrix.data()+8);
1361  v4uf row4(weightmatrix.data()+12);
1362 
1363  v4uf vcvs(&cvs[curseg]);
1364 
1365  for (int i = 0; i < len_t; i++)
1366  {
1367  {
1368  v4uf weights;
1369  float t2 = t*t;
1370  float t3 = t2*t;
1371 
1372  weights = row1;
1373  weights += row2 * t;
1374  weights += row3 * t2;
1375  weights += row4 * t3;
1376 
1377  v4uf v = vcvs * weights;
1378  v += v.swizzle<1, 1, 3, 3>();
1379  v += v.swizzle<2, 2, 2, 2>();
1380  results[i] = v[0];
1381  }
1382 
1383  t += step_t;
1384  if (t > 1)
1385  {
1386  while (curseg < nseg-1)
1387  {
1388  curseg++;
1389  t -= 1;
1390  if (t <= 1)
1391  break;
1392  }
1393  if (i < len_t-1)
1394  {
1395  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1396 
1397  row1 = v4uf(weightmatrix.data());
1398  row2 = v4uf(weightmatrix.data()+4);
1399  row3 = v4uf(weightmatrix.data()+8);
1400  row4 = v4uf(weightmatrix.data()+12);
1401 
1402  vcvs = v4uf(&cvs[curseg]);
1403  }
1404  }
1405  }
1406 #else
1407  for (int i = 0 ; i < len_t; i++)
1408  {
1409  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1410  t += step_t;
1411  if (t > 1)
1412  {
1413  while (curseg < nseg-1)
1414  {
1415  curseg++;
1416  t -= 1;
1417  if (t <= 1)
1418  break;
1419  }
1420  }
1421  }
1422 #endif
1423 }
1424 
1425 void
1426 UT_SplineCubic::enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1427 {
1428  // We need to find any minimum or maximum in each dimension
1429  // to enlarge the bounding box.
1430  // To do this, for each, dimension, we take the derivative
1431  // of the cubic, leaving a quadratic, and find the zeros of it.
1432  // The quadratic is such that its ith derivatives at zero are
1433  // the (i+1)th derivatives of the curve segment at zero.
1434  // a = (1/2) * 3rd derivative of curve segment at zero
1435  UT_Vector3 a = -cvs[0] + cvs[1] * 3.0F + cvs[2] * (-3.0F) + cvs[3];
1436  a *= 0.5F;
1437 
1438  // b = 2nd derivative of curve segment at zero
1439  // (this is equivalent to the 2nd difference)
1440  UT_Vector3 b = cvs[0] + cvs[1] * (-2.0F) + cvs[2];
1441  // c = 1st derivative of curve segment at zero
1442  // (this is equivalent to the central difference)
1443  UT_Vector3 c = cvs[2] - cvs[0];
1444  c *= 0.5F;
1445 
1446  enlargeBoundingBoxCommon<UT_SplineCubic::evalOpen<float> >(box, cvs, a, b, c, rootmin, rootmax);
1447 }
1448 
1449 void
1451 {
1452  // We need to find any minimum or maximum in each dimension
1453  // to enlarge the bounding box.
1454  // To do this, for each, dimension, we take the derivative
1455  // of the cubic, leaving a quadratic, and find the zeros of it.
1456  // The quadratic is such that its ith derivatives at zero are
1457  // the (i+1)th derivatives of the curve segment at zero.
1458 
1459  // First segment is (1 - t + (1/6)t^3)*P0 + (t - (1/3)*t^3)*P1 + ((1/6)t^3)*P2
1460  // 1st derivative is (-1 + (1/2)t^2)*P0 + (1 - t^2)*P1 + ((1/2)t^2)*P2
1461  // 2nd derivative is (t)*P0 + (-2t)*P1 + (t)*P2
1462  // 3rd derivative is (1)*P0 + (-2)*P1 + (1)*P2
1463 
1464  // a = (1/2) * 3rd derivative of curve segment at zero
1465  UT_Vector3 a = cvs[0] - 2.0f*cvs[1] + cvs[2];
1466  a *= 0.5F;
1467 
1468  // b = 2nd derivative of curve segment at zero
1469  // (this is equivalent to the 2nd difference)
1470  UT_Vector3 b(0,0,0);
1471  // c = 1st derivative of curve segment at zero
1472  // (this is equivalent to the central difference)
1473  UT_Vector3 c = cvs[1] - cvs[0];
1474 
1475  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDStart<float> >(box, cvs, a, b, c, rootmin, rootmax);
1476 }
1477 
1478 void
1479 UT_SplineCubic::enlargeBoundingBoxSubDEnd(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1480 {
1481  // We need to find any minimum or maximum in each dimension
1482  // to enlarge the bounding box.
1483  // To do this, for each, dimension, we take the derivative
1484  // of the cubic, leaving a quadratic, and find the zeros of it.
1485  // The quadratic is such that its ith derivatives at zero are
1486  // the (i+1)th derivatives of the curve segment at zero.
1487 
1488  // 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
1489  // 1st derivative is (-(1/2)(1-t)^2)*P0 + (-1 + (1-t)^2)*P1 + (1 - (1/2)(1-t)^2)*P2
1490  // 2nd derivative is (1-t)*P0 + (-2(1-t))*P1 + (1-t)*P2
1491  // 3rd derivative is (-1)*P0 + (2)*P1 + (-1)*P2
1492 
1493  // a = (1/2) * 3rd derivative of curve segment at zero
1494  // b = 2nd derivative of curve segment at zero
1495  // (this is equivalent to the 2nd difference)
1496  UT_Vector3 b = cvs[0] - 2.0f*cvs[1] + cvs[2];
1497  UT_Vector3 a = -0.5f*b;
1498 
1499  // c = 1st derivative of curve segment at zero
1500  // (this is equivalent to the central difference)
1501  UT_Vector3 c = cvs[2] - cvs[0];
1502  c *= 0.5f;
1503 
1504  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDEnd<float> >(box, cvs, a, b, c, rootmin, rootmax);
1505 }
1506 
1507 #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:1535
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
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
GLuint GLdouble GLdouble GLint GLint order
Definition: glew.h:3460
int getKnotLength() const
Definition: UT_Spline.h:732
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
const GLfloat * c
Definition: glew.h:16631
UT_Vector3T< float > UT_Vector3
#define UT_API
Definition: UT_API.h:14
void setGlobalBasis(UT_SPLINE_BASIS b)
Definition: UT_Spline.h:735
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: glew.h:12900
GLuint const GLchar * name
Definition: glcorearb.h:785
static T evalSubDStart(const T *cvs, float t)
Definition: UT_Spline.h:159
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
GLdouble GLdouble t
Definition: glew.h:1403
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
UT_Matrix4T< float > UT_Matrix4
GLint GLenum GLint x
Definition: glcorearb.h:408
static const UT_Matrix4 theInterpFirstBasis
Definition: UT_Spline.h:382
GLsizeiptr size
Definition: glcorearb.h:663
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
GLuint64EXT * result
Definition: glew.h:14311
int getVectorSize() const
Definition: UT_Spline.h:731
GLuint GLfloat GLfloat GLfloat x1
Definition: glew.h:12900
const GLdouble * v
Definition: glcorearb.h:836
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
GLsizei GLenum const void * indices
Definition: glcorearb.h:405
UT_Vector3T< T > SYSclamp(const UT_Vector3T< T > &v, const UT_Vector3T< T > &min, const UT_Vector3T< T > &max)
Definition: UT_Vector3.h:832
Definition: VM_SIMD.h:188
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
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:1426
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
#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
static const UT_Matrix4 theOpenBasis
Definition: UT_Spline.h:369
static T evalCubic(T kt, T dt, T iv, T im, T ov, T om, int order=0)
Definition: UT_Spline.h:836
static void enlargeBoundingBoxSubDStart(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
Definition: UT_Spline.h:1450
GLdouble n
Definition: glcorearb.h:2007
static T evalSubDCurve(const T *cvs, float t, int npts, bool deriv=false)
Definition: UT_Spline.h:303
GLuint GLfloat * val
Definition: glcorearb.h:1607
GLenum func
Definition: glcorearb.h:782
GLsizei const GLint box[]
Definition: glew.h:11654
UT_ColorType
Definition: UT_Color.h:24
const T * data() const
Return the raw matrix data.
Definition: UT_Matrix4.h:1100
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:1479
GLuint GLfloat x0
Definition: glew.h:12900
fpreal64 fpreal
Definition: SYS_Types.h:277
static T evalOpen(const T *cvs, float t)
Definition: UT_Spline.h:61
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition: glew.h:12900
OIIO_API bool copy(string_view from, string_view to, std::string &err)
GLbyte * weights
Definition: glew.h:7581
UT_Vector3T< T > colVecMult(const UT_Matrix3T< S > &m, const UT_Vector3T< T > &v)
Definition: UT_Matrix3.h:1495
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
GLsizei const GLfloat * value
Definition: glcorearb.h:823
GLfloat f
Definition: glcorearb.h:1925
SYS_FORCE_INLINE v4uf swizzle() const
Definition: VM_SIMD.h:335
static const UT_Matrix4 theInterpBasis
Definition: UT_Spline.h:386
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:1536
static const UT_Matrix4 theSubDFirstDerivBasis
Definition: UT_Spline.h:364
GLint y
Definition: glcorearb.h:102