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