HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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  /// Uniform knots with closed end conditions. seg is which segment
383  /// is being evaluates, nseg is the total. nseg should be
384  /// number of vertices minus three as we have cubics.
385  static UT_Matrix4 getClosedWeights(int seg, int nseg, bool deriv = false)
386  {
387  // these matrices come from $GEO/support/computespline.py
388  // which computes the power basis form of end-interpolated
389  // uniform bsplines.
390 
391  if (deriv == false)
392  {
393  if (nseg <= 1)
394  {
395  // Bezier.
396  return UT_Matrix4( 1, -3, 3, -1,
397  0, 3, -6, 3,
398  0, 0, 3, -3,
399  0, 0, 0, 1 );
400  }
401  else if (nseg == 2)
402  {
403  // 0, 0, 0, 1, 2, 2, 2,
404  if (seg == 0)
405  return UT_Matrix4( 1, -3, 3, -1,
406  0, 3, -4.5, 1.75,
407  0, 0, 1.5, -1,
408  0, 0, 0, 0.25 );
409  else
410  return UT_Matrix4( .25, -.75, .75, -0.25,
411  0.5, 0, -1.5, 1,
412  0.25, 0.75, 0.75, -1.75,
413  0, 0, 0, 1 );
414  }
415  else if (nseg == 3)
416  {
417  // 0, 0, 0, 1, 2, 3, 3, 3
418  if (seg == 0)
419  return UT_Matrix4( 1, -3, 3, -1,
420  0, 3, -4.5, 1.75,
421  0, 0, 1.5, -11/12.,
422  0, 0, 0, 1/6.);
423  else if (seg == 1)
424  return UT_Matrix4( .25, -.75, .75, -0.25,
425  7/12., 0.25, -1.25, 7/12.,
426  1/6., 0.5, 0.5, -7/12.,
427  0, 0, 0, 0.25 );
428  else
429  return UT_Matrix4( 1/6., -.5, .5, -1/6.,
430  7/12., -0.25, -1.25, 11/12.,
431  0.25, 0.75, 0.75, -1.75,
432  0, 0, 0, 1 );
433  }
434  else
435  {
436  // Either on an end, or in the middle
437  if (seg >= 2 && seg < nseg-2)
438  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
439  4/6., 0/6., -6/6., 3/6.,
440  1/6., 3/6., 3/6., -3/6.,
441  0/6., 0/6., 0/6., 1/6. );
442  else if (seg == 0)
443  return UT_Matrix4( 1, -3, 3, -1,
444  0, 3, -4.5, 1.75,
445  0, 0, 1.5, -11/12.,
446  0, 0, 0, 1/6. );
447  else if (seg == 1)
448  return UT_Matrix4( 0.25, -0.75, 0.75, -.25,
449  7/12., 0.25, -1.25, 7/12.,
450  1/6., 0.5, 0.5, -0.5,
451  0, 0, 0, 1/6. );
452  else if (seg == nseg-2)
453  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
454  2/3., 0, -1, 0.5,
455  1/6., 0.5, 0.5, -7/12.,
456  0, 0, 0, 0.25 );
457  else // if (seg == nseg-1)
458  return UT_Matrix4( 1/6., -3/6., 3/6., -1/6.,
459  7/12., -.25, -1.25, 11/12.,
460  0.25, 0.75, 0.75, -1.75,
461  0, 0, 0, 1 );
462  }
463  }
464  else
465  {
466  if (nseg <= 1)
467  {
468  // Bezier.
469  return UT_Matrix4( -3, 6, -3, 0,
470  3, -12, 9, 0,
471  0, 6, -9, 0,
472  0, 0, 3, 0 );
473  }
474  else if (nseg == 2)
475  {
476  // 0, 0, 0, 1, 2, 2, 2,
477  if (seg == 0)
478  return UT_Matrix4( -3, 6, -3, 0,
479  3, -9, 5.25, 0,
480  0, 3, -3, 0,
481  0, 0, 0.75, 9 );
482  else
483  return UT_Matrix4( -.75, 1.5, -0.75, 0,
484  0, -3, 3, 0,
485  0.75, 1.5, -5.25, 0,
486  0, 0, 3, 0 );
487  }
488  else if (nseg == 3)
489  {
490  // 0, 0, 0, 1, 2, 3, 3, 3
491  if (seg == 0)
492  return UT_Matrix4( -3, 6, -3, 0,
493  3, -9, 5.25, 0,
494  0, 3, -11/4., 0,
495  0, 0, .5, 0);
496  else if (seg == 1)
497  return UT_Matrix4( -.75, 1.5, -0.75, 0,
498  0.25, -2.5, 7/4., 0,
499  0.5, 1, -7/4., 0,
500  0, 0, 0.75, 0);
501  else
502  return UT_Matrix4( -.5, 1, -0.5, 0,
503  -0.25, -2.5, 11/4., 0,
504  0.75, 1.5, -5.25, 0,
505  0, 0, 3, 0);
506  }
507  else
508  {
509  // Either on an end, or in the middle
510  if (seg >= 2 && seg < nseg-2)
511  return UT_Matrix4( -3/6., 1.0, -0.5, 0,
512  0/6., -2.0, 1.5, 0,
513  3/6., 1.0, -1.5, 0,
514  0/6., 0, 0.5, 0);
515  else if (seg == 0)
516  return UT_Matrix4( -3, 6, -3, 0,
517  3, -9, 5.25, 0,
518  0, 3, -11/4., 0,
519  0, 0, 0.5, 0);
520  else if (seg == 1)
521  return UT_Matrix4( -0.75, 1.5, -.75, 0,
522  0.25, -2.5, 7/4., 0,
523  0.5, 1, -1.5, 0,
524  0, 0, 0.5, 0 );
525  else if (seg == nseg-2)
526  return UT_Matrix4( -3/6., 1, -0.5, 0,
527  0, -2, 1.5, 0,
528  0.5, 1, -7/4., 0,
529  0, 0, 0.75, 0 );
530  else // if (seg == nseg-1)
531  return UT_Matrix4(-3/6., 1, -.5, 0,
532  -.25, -2.5, 11/4., 0,
533  0.75, 1.5, -5.25, 0,
534  0, 0, 3, 0);
535  }
536  }
537  }
538  static UT_Matrix4 getClosedWeightsTranspose(int seg, int nseg, bool deriv = false)
539  {
540  if (deriv == false)
541  {
542  // these matrices come from $GEO/support/computespline.py
543  // which computes the power basis form of end-interpolated
544  // uniform bsplines.
545  if (nseg <= 1)
546  {
547  // Bezier.
548  return UT_Matrix4( 1, 0, 0, 0,
549  -3, 3, 0, 0,
550  3, -6, 3, 0,
551  -1, 3, -3, 1 );
552  }
553  else if (nseg == 2)
554  {
555  // 0, 0, 0, 1, 2, 2, 2,
556  if (seg == 0)
557  return UT_Matrix4( 1, 0, 0, 0,
558  -3, 3, 0, 0,
559  3, -4.5, 1.5, 0,
560  -1, 1.75, -1, 0.25 );
561  else
562  return UT_Matrix4( .25, .5, .25, 0,
563  -.75, 0, .75, 0,
564  0.75, -1.5, 0.75, 0,
565  -0.25, 1, -1.75, 1 );
566  }
567  else if (nseg == 3)
568  {
569  // 0, 0, 0, 1, 2, 3, 3, 3
570  if (seg == 0)
571  return UT_Matrix4( 1, 0, 0, 0,
572  -3, 3, 0, 0,
573  3,-4.5,1.5, 0,
574  -1,1.75,-11/12.,1/6. );
575  else if (seg == 1)
576  return UT_Matrix4( 0.25, 7/12., 1/6., 0,
577  -.75, 0.25, 0.5, 0,
578  .75,-1.25, 0.5, 0,
579  -.25,7/12.,-7/12.,0.25 );
580  else
581  return UT_Matrix4( 1/6., 7/12., 0.25, 0,
582  -.5, -0.25, 0.75, 0,
583  .5, -1.25, 0.75, 0,
584  -1/6.,11/12.,-1.75, 1 );
585 
586  }
587  else
588  {
589  // Either on an end, or in the middle
590  if (seg >= 2 && seg < nseg-2)
591  return UT_Matrix4( 1/6., 4/6., 1/6., 0/6.,
592  -3/6., 0/6., 3/6., 0/6.,
593  3/6., -6/6., 3/6., 0/6.,
594  -1/6., 3/6., -3/6., 1/6. );
595  else if (seg == 0)
596  return UT_Matrix4( 1, 0, 0, 0,
597  -3, 3, 0, 0,
598  3,-4.5,1.5, 0,
599  -1,1.75,-11/12., 1/6. );
600  else if (seg == 1)
601  return UT_Matrix4( 0.25, 7/12., 1/6., 0,
602  -0.75, 0.25, 0.5, 0,
603  0.75,-1.25, 0.5, 0,
604  -0.25,7/12., -0.5, 1/6. );
605  else if (seg == nseg-2)
606  return UT_Matrix4( 1/6., 2/3., 1/6., 0,
607  -3/6., 0, 0.5, 0,
608  3/6., -1, 0.5, 0,
609  -1/6., 0.5,-7/12.,0.25 );
610  else // if (seg == nseg-1)
611  return UT_Matrix4( 1/6., 7/12., 0.25, 0,
612  -3/6., -.25, 0.75, 0,
613  3/6.,-1.25, 0.75, 0,
614  -1/6.,11/12.,-1.75, 1 );
615  }
616  }
617  else
618  {
619  if (nseg <= 1)
620  {
621  // Bezier.
622  return UT_Matrix4(-3, 3, 0, 0,
623  6, -12, 6, 0,
624  -3, 9, -9, 3,
625  0, 0, 0, 0);
626  }
627  else if (nseg == 2)
628  {
629  // 0, 0, 0, 1, 2, 2, 2,
630  if (seg == 0)
631  return UT_Matrix4(-3, 3, 0, 0,
632  6, -9, 3, 0,
633  -3, 5.25, -3, 0.75,
634  0, 0, 0, 0);
635  else
636  return UT_Matrix4(-.75, 0, .75, 0,
637  1.5, -3, 1.5, 0,
638  -0.75, 3, -5.25, 3,
639  0, 0, 0, 0);
640  }
641  else if (nseg == 3)
642  {
643  // 0, 0, 0, 1, 2, 3, 3, 3
644  if (seg == 0)
645  return UT_Matrix4(-3, 3, 0, 0,
646  6, -9, 3, 0,
647  -3, 5.25, -11/4., .5,
648  0, 0, 0, 0);
649  else if (seg == 1)
650  return UT_Matrix4(-.75, 0.25, 0.5, 0,
651  1.5,-2.5, 1, 0,
652  -.75,7/4.,-7/4.,0.75,
653  0, 0, 0, 0);
654  else
655  return UT_Matrix4(-.5, -0.25, 0.75, 0,
656  1, -2.5, 1.5, 0,
657  -.5, 11/4., -5.25, 3,
658  0, 0, 0, 0);
659 
660  }
661  else
662  {
663  // Either on an end, or in the middle
664  if (seg >= 2 && seg < nseg-2)
665  return UT_Matrix4(-3/6., 0/6., 3/6., 0/6.,
666  1, -2, 1, 0,
667  -0.5, 1.5, -1.5, 0.5,
668  0, 0, 0, 0);
669  else if (seg == 0)
670  return UT_Matrix4(-3, 3, 0, 0,
671  6, -9, 3, 0,
672  -3, 5.25, -11/4., .5,
673  0, 0, 0, 0);
674  else if (seg == 1)
675  return UT_Matrix4(-0.75, 0.25, 0.5, 0,
676  1.5, -2.5, 1, 0,
677  -0.75, 7/4., -1.5, .5,
678  0, 0, 0, 0);
679 
680  else if (seg == nseg-2)
681  return UT_Matrix4(-3/6., 0, 0.5, 0,
682  1, -2, 1, 0,
683  -0.5, 1.5, -7/4., 0.75,
684  0, 0, 0, 0);
685  else // if (seg == nseg-1)
686  return UT_Matrix4(-3/6., -.25, 0.75, 0,
687  1, -2.5, 1.5, 0,
688  -.5, 11/4.,-5.25, 3,
689  0, 0, 0, 0);
690  }
691  }
692  }
693 };
694 
695 
696 /// The Linear & Catmull-Rom splines expect a parametric coordinate for
697 /// evaluation between 0 and 1. The Catmull-Rom spline requires additional
698 /// key values at the beginning and end of the spline to evaluate the slopes
699 /// of the Hermite spline properly.
700 ///
701 /// The LinearSolve only works on scalar values. It will compute the
702 /// parametric coordinate associated with the value passed in. This can be
703 /// used to simulate non-uniform keys on the spline.
705 public:
706  UT_Spline();
707  ~UT_Spline();
708 
709  /// Return the amount of memory owned by this UT_Spline in bytes
710  int64 getMemoryUsage(bool inclusive) const;
711 
712  /// Query the basis or knot length of the spline
713  UT_SPLINE_BASIS getGlobalBasis() const { return myGlobalBasis; }
714  int getVectorSize() const { return myVectorSize; }
715  int getKnotLength() const { return myKnotLength; }
716  fpreal64 getTension() const { return myTension; }
717 
719  { myGlobalBasis = b; }
720 
721  /// Construction of the spline object. All values are initialized to 0.
722  /// Warning, calling setSize() will clear all existing values.
723  void setSize(int nkeys, int vector_size);
724  /// Cubic splines may have a "tension". The tension defaults to 0.5 which
725  /// results in Catmull-Rom splines.
726  void setTension(fpreal64 t);
727 
728  /// Once the spline has been constructed, the values need to be set.
729  /// It is possible to change values between evaluations.
730  // @{
731  void setValue(int key, const fpreal32 *value, int size);
732  void setValue(int key, const fpreal64 *value, int size);
733  // @}
734 
735  /// Set the basis for the given key index.
736  /// This will also set the global basis.
737  void setBasis(int key, UT_SPLINE_BASIS b);
738 
739  /// Evaluate the spline using the global basis.
740  /// When interp_space is not UT_RGB, then values are treated as UT_RGBA
741  /// and converted into the desired color space before being interpolated.
742  /// The result is always returned as UT_RGBA.
743  // @{
744  bool evaluate(fpreal t, fpreal32 *result, int size,
745  UT_ColorType interp_space) const;
746  bool evaluate(fpreal t, fpreal64 *result, int size,
747  UT_ColorType interp_space) const;
748  // @}
749 
750  /// Evaluate the spline using multiple basis types depending on t.
751  /// Also unlike evaluate(), evaluateMulti() doesn't require extra keys for
752  /// Catmull-Rom on the ends. It always evaluates using a 0 slope at the
753  /// ends.
754  // @{
755  bool evaluateMulti(fpreal t, fpreal32 *result, int n,
756  UT_ColorType interp_space,
757  int knot_segment_hint = -1) const;
758  bool evaluateMulti(fpreal t, fpreal64 *result, int n,
759  UT_ColorType interp_space,
760  int knot_segment_hint = -1) const;
761  // @}
762 
763  /// Return _monotone_ cubic Hermite slopes at the current knot given the
764  /// previous and next knots.
765  // @{
766  /// Fritsch-Carlson (local) method.
767  /// It gives relatively good looking C1 curves but might not give a C2
768  /// solution even if it exists.
769  /// 1. Fritsch, F.N., Carlson, R.E., Monotone piecewise cubic interpolant,
770  /// SIAM J. Numer. Anal., 17(1980), 238-246.
771  static fpreal64 getMonotoneSlopeFC(fpreal64 v_cur, fpreal64 t_cur,
772  fpreal64 v_prev, fpreal64 t_prev,
773  fpreal64 v_next, fpreal64 t_next);
774  /// Paul Kupan's method.
775  /// Similar to Fritsch-Carlson method except it gives more visually
776  /// pleasing results when the intervals next to the knot are uneven.
777  /// 2. Kupan, P.A., Monotone Interpolant Built with Slopes Obtained by
778  /// Linear Combination, Studia Universitatis Babes-Bolyai Mathematica,
779  /// 53(2008), 59-66.
780  static fpreal64 getMonotoneSlopePA(fpreal64 v_cur, fpreal64 t_cur,
781  fpreal64 v_prev, fpreal64 t_prev,
782  fpreal64 v_next, fpreal64 t_next);
783  // @}
784 
785  /// Given the position within the two knots and the first knot index
786  /// number, normalize the position from knot-length domain to unit domain.
787  fpreal normalizeParameter(fpreal parm, fpreal t0, fpreal t1,
788  int t0_index) const;
789 
790  // Given parm in [0,1] interval in which CVs are at 'knots' values,
791  // find the parameter t in [0,1] interval in which CVs are evenly spaced.
792  // The returned reparameterized value is calculated in such a way that it
793  // yields smooth curve at CVs (unlike normalizeation that uses linear
794  // interpolation that yields tangent disconinuity at CVs).
795  // If 'parm_knot_segment' is given, this function sets it to the index
796  // of the knot segment into which the parm falls; it can be used as a hint
797  // to the evaluateMulti() method, especially for b-splines that may remap
798  // parameters to an adjacent segment (for continuity at end control points).
799  fpreal solveSpline(fpreal parm,
800  const UT_FprealArray &knots,
801  int *parm_knot_segment = NULL) const;
802 
803  /// Given the keys surrounding a channel segment, evaluate it as cubic
804  /// hermite spline. This function assumes dt > 0.
805  /// kt is [0,1] which maps over the dt time interval.
806  template <typename T>
807  static T evalCubic(T kt, T dt, T iv, T im, T ov, T om)
808  {
809  T x0 = iv;
810  T x1 = im*dt;
811  T x3 = 2*(iv - ov) + (im + om)*dt;
812  T x2 = ov - iv - im*dt - x3;
813  return x0 + kt*(x1 + kt*(x2 + kt*x3));
814  }
815 
816 private:
817  UT_Spline(const UT_Spline &copy); // not implemented yet
818  UT_Spline &operator =(const UT_Spline &copy); // not implemented yet
819 
820 private:
821  void grow(int size);
822  int getInterp(fpreal t, int knots[], fpreal weights[]);
823 
824  int64 getSizeOfValues() const
825  { return myKnotLength*myVectorSize*sizeof(fpreal64); }
826  int64 getSizeOfBases() const
827  { return myKnotLength*sizeof(UT_SPLINE_BASIS); }
828 
829  template <typename T>
830  inline void combineKeys(T *result, int vector_size,
831  fpreal64 weights[], int indices[],
832  int num_indices,
833  UT_ColorType interp_space) const;
834 
835  template <typename T>
836  inline void evalMonotoneCubic(fpreal t, T *result,
837  int n, UT_ColorType interp_space,
838  bool do_multi) const;
839 
840  template <typename T>
841  inline void setValueInternal(int key, const T *value, int size);
842 
843  template <typename T>
844  inline bool evaluateInternal(fpreal t, T *result, int size,
845  UT_ColorType interp_space) const;
846 
847  template <typename T>
848  inline bool evaluateMultiInternal(fpreal t, T *result, int n,
849  UT_ColorType interp_space,
850  int knot_segment_hint) const;
851 
852  fpreal64 *myValues;
853  UT_SPLINE_BASIS *myBases;
854  fpreal64 myTension;
855  int myVectorSize;
856  int myKnotLength;
857  UT_SPLINE_BASIS myGlobalBasis;
858 };
859 
860 #include <VM/VM_SIMD.h>
861 
862 template <>
865 {
866 #if defined(CPU_HAS_SIMD_INSTR)
867  v4uf row1(1/6., 4/6., 1/6., 0/6.);
868  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
869  v4uf row3(3/6., -6/6., 3/6., 0/6.);
870  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
871 
872  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
873  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
874  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
875  v4uf vt(t);
876  v4uf vt2 = vt*vt;
877  v4uf vt3 = vt2*vt;
878 
879  v4uf weights;
880 
881  weights = row1;
882  weights += row2 * vt;
883  weights += row3 * vt2;
884  weights += row4 * vt3;
885 
886  vcvsx *= weights;
887  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
888  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
889  vcvsy *= weights;
890  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
891  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
892  vcvsz *= weights;
893  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
894  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
895 
896  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
897 #else
898  UT_Matrix4 weightmatrix = getOpenWeights();
899  float t2 = t*t;
900  float t3 = t2*t;
901  UT_Vector4 powers(1, t, t2, t3);
902 
903  UT_Vector4 weights = colVecMult(weightmatrix, powers);
904 
906 
907  value = cvs[0] * weights[0];
908  value += cvs[1] * weights[1];
909  value += cvs[2] * weights[2];
910  value += cvs[3] * weights[3];
911 
912  return value;
913 #endif
914 }
915 
916 template <>
917 SYS_FORCE_INLINE float
918 UT_SplineCubic::evalOpen(const float *cvs, float t)
919 {
920 #if defined(CPU_HAS_SIMD_INSTR)
921  v4uf row1(1/6., 4/6., 1/6., 0/6.);
922  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
923  v4uf row3(3/6., -6/6., 3/6., 0/6.);
924  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
925 
926  v4uf vcvs(cvs);
927  v4uf vt(t);
928  v4uf vt2 = vt*vt;
929  v4uf vt3 = vt2*vt;
930 
931  v4uf weights;
932 
933  weights = row1;
934  weights += row2 * vt;
935  weights += row3 * vt2;
936  weights += row4 * vt3;
937 
938  vcvs *= weights;
939  vcvs += vcvs.swizzle<1, 1, 3, 3>();
940  vcvs += vcvs.swizzle<2, 2, 2, 2>();
941 
942  return vcvs[0];
943 #else
944  UT_Matrix4 weightmatrix = getOpenWeights();
945  float t2 = t*t;
946  float t3 = t2*t;
947  UT_Vector4 powers(1, t, t2, t3);
948 
949  UT_Vector4 weights = colVecMult(weightmatrix, powers);
950 
951  float value;
952 
953  value = cvs[0] * weights[0];
954  value += cvs[1] * weights[1];
955  value += cvs[2] * weights[2];
956  value += cvs[3] * weights[3];
957 
958  return value;
959 #endif
960 }
961 
962 template <>
963 inline void
964 UT_SplineCubic::evalRangeOpen(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg)
965 {
966  int curseg;
967  curseg = SYSfastFloor(start_t);
968  curseg = SYSclamp(curseg, 0, nseg-1);
969  float t = start_t - curseg;
970 
971 #if defined(CPU_HAS_SIMD_INSTR)
972  v4uf row1(1/6., 4/6., 1/6., 0/6.);
973  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
974  v4uf row3(3/6., -6/6., 3/6., 0/6.);
975  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
976 
977  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
978  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
979  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
980 
981  for (int i = 0; i < len_t; i++)
982  {
983  {
984  v4uf weights;
985  float t2 = t*t;
986  float t3 = t2*t;
987 
988  weights = row1;
989  weights += row2 * t;
990  weights += row3 * t2;
991  weights += row4 * t3;
992 
993  v4uf vx = vcvsx * weights;
994  vx += vx.swizzle<1, 1, 3, 3>();
995  vx += vx.swizzle<2, 2, 2, 2>();
996  v4uf vy = vcvsy * weights;
997  vy += vy.swizzle<1, 1, 3, 3>();
998  vy += vy.swizzle<2, 2, 2, 2>();
999  v4uf vz = vcvsz * weights;
1000  vz += vz.swizzle<1, 1, 3, 3>();
1001  vz += vz.swizzle<2, 2, 2, 2>();
1002  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1003  }
1004 
1005  t += step_t;
1006  if (t > 1)
1007  {
1008  while (curseg < nseg-1)
1009  {
1010  curseg++;
1011  t -= 1;
1012  if (t <= 1)
1013  break;
1014  }
1015  if (i < len_t-1)
1016  {
1017  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1018  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1019  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1020  }
1021  }
1022  }
1023 #else
1024  for (int i = 0; i < len_t; i++)
1025  {
1026  results[i] = evalOpen(&cvs[curseg], t);
1027  t += step_t;
1028  if (t > 1)
1029  {
1030  while (curseg < nseg-1)
1031  {
1032  curseg++;
1033  t -= 1;
1034  if (t <= 1)
1035  break;
1036  }
1037  }
1038  }
1039 #endif
1040 }
1041 
1042 template <>
1043 inline void
1044 UT_SplineCubic::evalRangeOpen(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg)
1045 {
1046  int curseg;
1047  curseg = SYSfastFloor(start_t);
1048  curseg = SYSclamp(curseg, 0, nseg-1);
1049  float t = start_t - curseg;
1050 
1051 #if defined(CPU_HAS_SIMD_INSTR)
1052  v4uf row1(1/6., 4/6., 1/6., 0/6.);
1053  v4uf row2(-3/6., 0/6., 3/6., 0/6.);
1054  v4uf row3(3/6., -6/6., 3/6., 0/6.);
1055  v4uf row4(-1/6., 3/6., -3/6., 1/6. );
1056 
1057  v4uf vcvs(&cvs[curseg]);
1058 
1059  for (int i = 0; i < len_t; i++)
1060  {
1061  {
1062  v4uf weights;
1063  float t2 = t*t;
1064  float t3 = t2*t;
1065 
1066  weights = row1;
1067  weights += row2 * t;
1068  weights += row3 * t2;
1069  weights += row4 * t3;
1070 
1071  v4uf v = vcvs * weights;
1072  v += v.swizzle<1, 1, 3, 3>();
1073  v += v.swizzle<2, 2, 2, 2>();
1074  results[i] = v[0];
1075  }
1076 
1077  t += step_t;
1078  if (t > 1)
1079  {
1080  while (curseg < nseg-1)
1081  {
1082  curseg++;
1083  t -= 1;
1084  if (t <= 1)
1085  break;
1086  }
1087  if (i < len_t-1)
1088  {
1089  vcvs = v4uf(&cvs[curseg]);
1090  }
1091  }
1092  }
1093 #else
1094  for (int i = 0 ; i < len_t; i++)
1095  {
1096  results[i] = evalOpen(&cvs[curseg], t);
1097  t += step_t;
1098  if (t > 1)
1099  {
1100  while (curseg < nseg-1)
1101  {
1102  curseg++;
1103  t -= 1;
1104  if (t <= 1)
1105  break;
1106  }
1107  }
1108  }
1109 #endif
1110 }
1111 
1112 template <>
1114 UT_SplineCubic::evalClosed(const UT_Vector3 *cvs, float t, int seg, int nseg, bool deriv)
1115 {
1116 #if defined(CPU_HAS_SIMD_INSTR)
1117  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1118 
1119  v4uf row1(weightmatrix.data());
1120  v4uf row2(weightmatrix.data()+4);
1121  v4uf row3(weightmatrix.data()+8);
1122  v4uf row4(weightmatrix.data()+12);
1123 
1124  v4uf vcvsx(cvs[0].x(), cvs[1].x(), cvs[2].x(), cvs[3].x());
1125  v4uf vcvsy(cvs[0].y(), cvs[1].y(), cvs[2].y(), cvs[3].y());
1126  v4uf vcvsz(cvs[0].z(), cvs[1].z(), cvs[2].z(), cvs[3].z());
1127  v4uf vt(t);
1128  v4uf vt2 = vt*vt;
1129  v4uf vt3 = vt2*vt;
1130 
1131  v4uf weights;
1132 
1133  weights = row1;
1134  weights += row2 * vt;
1135  weights += row3 * vt2;
1136  weights += row4 * vt3;
1137 
1138  vcvsx *= weights;
1139  vcvsx += vcvsx.swizzle<1, 1, 3, 3>();
1140  vcvsx += vcvsx.swizzle<2, 2, 2, 2>();
1141  vcvsy *= weights;
1142  vcvsy += vcvsy.swizzle<1, 1, 3, 3>();
1143  vcvsy += vcvsy.swizzle<2, 2, 2, 2>();
1144  vcvsz *= weights;
1145  vcvsz += vcvsz.swizzle<1, 1, 3, 3>();
1146  vcvsz += vcvsz.swizzle<2, 2, 2, 2>();
1147 
1148  return UT_Vector3( vcvsx[0], vcvsy[0], vcvsz[0] );
1149 #else
1150  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1151  float t2 = t*t;
1152  float t3 = t2*t;
1153  UT_Vector4 powers(1, t, t2, t3);
1154 
1155  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1156 
1157  UT_Vector3 value;
1158 
1159  value = cvs[0] * weights[0];
1160  value += cvs[1] * weights[1];
1161  value += cvs[2] * weights[2];
1162  value += cvs[3] * weights[3];
1163 
1164  return value;
1165 #endif
1166 }
1167 
1168 template <>
1169 SYS_FORCE_INLINE float
1170 UT_SplineCubic::evalClosed(const float *cvs, float t, int seg, int nseg, bool deriv)
1171 {
1172 #if defined(CPU_HAS_SIMD_INSTR)
1173  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(seg, nseg, deriv);
1174 
1175  v4uf row1(weightmatrix.data());
1176  v4uf row2(weightmatrix.data()+4);
1177  v4uf row3(weightmatrix.data()+8);
1178  v4uf row4(weightmatrix.data()+12);
1179 
1180  v4uf vcvs(cvs);
1181  float t2 = t*t;
1182  float t3 = t2*t;
1183 
1184  v4uf weights;
1185 
1186  weights = row1;
1187  weights += row2 * t;
1188  weights += row3 * t2;
1189  weights += row4 * t3;
1190 
1191  vcvs *= weights;
1192 
1193  vcvs += vcvs.swizzle<1, 1, 3, 3>();
1194  vcvs += vcvs.swizzle<2, 2, 2, 2>();
1195 
1196  return vcvs[0];
1197 #else
1198  UT_Matrix4 weightmatrix = getClosedWeights(seg, nseg);
1199  float t2 = t*t;
1200  float t3 = t2*t;
1201  UT_Vector4 powers(1, t, t2, t3);
1202 
1203  UT_Vector4 weights = colVecMult(weightmatrix, powers);
1204 
1205  float value;
1206 
1207  value = cvs[0] * weights[0];
1208  value += cvs[1] * weights[1];
1209  value += cvs[2] * weights[2];
1210  value += cvs[3] * weights[3];
1211 
1212  return value;
1213 #endif
1214 }
1215 
1216 template <>
1217 inline void
1218 UT_SplineCubic::evalRangeClosed(UT_Vector3 *results, const UT_Vector3 *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1219 {
1220  int curseg;
1221  curseg = SYSfastFloor(start_t);
1222  curseg = SYSclamp(curseg, 0, nseg-1);
1223  float t = start_t - curseg;
1224 
1225 #if defined(CPU_HAS_SIMD_INSTR)
1226  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1227 
1228  v4uf row1(weightmatrix.data());
1229  v4uf row2(weightmatrix.data()+4);
1230  v4uf row3(weightmatrix.data()+8);
1231  v4uf row4(weightmatrix.data()+12);
1232 
1233  v4uf vcvsx(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1234  v4uf vcvsy(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1235  v4uf vcvsz(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1236 
1237  for (int i = 0; i < len_t; i++)
1238  {
1239  {
1240  v4uf weights;
1241  float t2 = t*t;
1242  float t3 = t2*t;
1243 
1244  weights = row1;
1245  weights += row2 * t;
1246  weights += row3 * t2;
1247  weights += row4 * t3;
1248 
1249  v4uf vx = vcvsx * weights;
1250  vx += vx.swizzle<1, 1, 3, 3>();
1251  vx += vx.swizzle<2, 2, 2, 2>();
1252  v4uf vy = vcvsy * weights;
1253  vy += vy.swizzle<1, 1, 3, 3>();
1254  vy += vy.swizzle<2, 2, 2, 2>();
1255  v4uf vz = vcvsz * weights;
1256  vz += vz.swizzle<1, 1, 3, 3>();
1257  vz += vz.swizzle<2, 2, 2, 2>();
1258  results[i] = UT_Vector3( vx[0], vy[0], vz[0] );
1259  }
1260 
1261  t += step_t;
1262  if (t > 1)
1263  {
1264  while (curseg < nseg-1)
1265  {
1266  curseg++;
1267  t -= 1;
1268  if (t <= 1)
1269  break;
1270  }
1271  if (i < len_t-1)
1272  {
1273  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1274 
1275  row1 = v4uf(weightmatrix.data());
1276  row2 = v4uf(weightmatrix.data()+4);
1277  row3 = v4uf(weightmatrix.data()+8);
1278  row4 = v4uf(weightmatrix.data()+12);
1279 
1280  vcvsx = v4uf(cvs[curseg].x(), cvs[curseg+1].x(), cvs[curseg+2].x(), cvs[curseg+3].x());
1281  vcvsy = v4uf(cvs[curseg].y(), cvs[curseg+1].y(), cvs[curseg+2].y(), cvs[curseg+3].y());
1282  vcvsz = v4uf(cvs[curseg].z(), cvs[curseg+1].z(), cvs[curseg+2].z(), cvs[curseg+3].z());
1283  }
1284  }
1285  }
1286 #else
1287  for (int i = 0; i < len_t; i++)
1288  {
1289  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1290  t += step_t;
1291  if (t > 1)
1292  {
1293  while (curseg < nseg-1)
1294  {
1295  curseg++;
1296  t -= 1;
1297  if (t <= 1)
1298  break;
1299  }
1300  }
1301  }
1302 #endif
1303 }
1304 
1305 template <>
1306 inline void
1307 UT_SplineCubic::evalRangeClosed(float *results, const float *cvs, float start_t, float step_t, int len_t, int nseg, bool deriv)
1308 {
1309  int curseg;
1310  curseg = SYSfastFloor(start_t);
1311  curseg = SYSclamp(curseg, 0, nseg-1);
1312  float t = start_t - curseg;
1313 
1314 #if defined(CPU_HAS_SIMD_INSTR)
1315  UT_Matrix4 weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1316 
1317  v4uf row1(weightmatrix.data());
1318  v4uf row2(weightmatrix.data()+4);
1319  v4uf row3(weightmatrix.data()+8);
1320  v4uf row4(weightmatrix.data()+12);
1321 
1322  v4uf vcvs(&cvs[curseg]);
1323 
1324  for (int i = 0; i < len_t; i++)
1325  {
1326  {
1327  v4uf weights;
1328  float t2 = t*t;
1329  float t3 = t2*t;
1330 
1331  weights = row1;
1332  weights += row2 * t;
1333  weights += row3 * t2;
1334  weights += row4 * t3;
1335 
1336  v4uf v = vcvs * weights;
1337  v += v.swizzle<1, 1, 3, 3>();
1338  v += v.swizzle<2, 2, 2, 2>();
1339  results[i] = v[0];
1340  }
1341 
1342  t += step_t;
1343  if (t > 1)
1344  {
1345  while (curseg < nseg-1)
1346  {
1347  curseg++;
1348  t -= 1;
1349  if (t <= 1)
1350  break;
1351  }
1352  if (i < len_t-1)
1353  {
1354  weightmatrix = getClosedWeightsTranspose(curseg, nseg, deriv);
1355 
1356  row1 = v4uf(weightmatrix.data());
1357  row2 = v4uf(weightmatrix.data()+4);
1358  row3 = v4uf(weightmatrix.data()+8);
1359  row4 = v4uf(weightmatrix.data()+12);
1360 
1361  vcvs = v4uf(&cvs[curseg]);
1362  }
1363  }
1364  }
1365 #else
1366  for (int i = 0 ; i < len_t; i++)
1367  {
1368  results[i] = evalClosed(&cvs[curseg], t, curseg, nseg, deriv);
1369  t += step_t;
1370  if (t > 1)
1371  {
1372  while (curseg < nseg-1)
1373  {
1374  curseg++;
1375  t -= 1;
1376  if (t <= 1)
1377  break;
1378  }
1379  }
1380  }
1381 #endif
1382 }
1383 
1384 void
1385 UT_SplineCubic::enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1386 {
1387  // We need to find any minimum or maximum in each dimension
1388  // to enlarge the bounding box.
1389  // To do this, for each, dimension, we take the derivative
1390  // of the cubic, leaving a quadratic, and find the zeros of it.
1391  // The quadratic is such that its ith derivatives at zero are
1392  // the (i+1)th derivatives of the curve segment at zero.
1393  // a = (1/2) * 3rd derivative of curve segment at zero
1394  UT_Vector3 a = -cvs[0] + cvs[1] * 3.0F + cvs[2] * (-3.0F) + cvs[3];
1395  a *= 0.5F;
1396 
1397  // b = 2nd derivative of curve segment at zero
1398  // (this is equivalent to the 2nd difference)
1399  UT_Vector3 b = cvs[0] + cvs[1] * (-2.0F) + cvs[2];
1400  // c = 1st derivative of curve segment at zero
1401  // (this is equivalent to the central difference)
1402  UT_Vector3 c = cvs[2] - cvs[0];
1403  c *= 0.5F;
1404 
1405  enlargeBoundingBoxCommon<UT_SplineCubic::evalOpen<float> >(box, cvs, a, b, c, rootmin, rootmax);
1406 }
1407 
1408 void
1409 UT_SplineCubic::enlargeBoundingBoxSubDStart(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1410 {
1411  // We need to find any minimum or maximum in each dimension
1412  // to enlarge the bounding box.
1413  // To do this, for each, dimension, we take the derivative
1414  // of the cubic, leaving a quadratic, and find the zeros of it.
1415  // The quadratic is such that its ith derivatives at zero are
1416  // the (i+1)th derivatives of the curve segment at zero.
1417 
1418  // First segment is (1 - t + (1/6)t^3)*P0 + (t - (1/3)*t^3)*P1 + ((1/6)t^3)*P2
1419  // 1st derivative is (-1 + (1/2)t^2)*P0 + (1 - t^2)*P1 + ((1/2)t^2)*P2
1420  // 2nd derivative is (t)*P0 + (-2t)*P1 + (t)*P2
1421  // 3rd derivative is (1)*P0 + (-2)*P1 + (1)*P2
1422 
1423  // a = (1/2) * 3rd derivative of curve segment at zero
1424  UT_Vector3 a = cvs[0] - 2.0f*cvs[1] + cvs[2];
1425  a *= 0.5F;
1426 
1427  // b = 2nd derivative of curve segment at zero
1428  // (this is equivalent to the 2nd difference)
1429  UT_Vector3 b(0,0,0);
1430  // c = 1st derivative of curve segment at zero
1431  // (this is equivalent to the central difference)
1432  UT_Vector3 c = cvs[1] - cvs[0];
1433 
1434  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDStart<float> >(box, cvs, a, b, c, rootmin, rootmax);
1435 }
1436 
1437 void
1438 UT_SplineCubic::enlargeBoundingBoxSubDEnd(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
1439 {
1440  // We need to find any minimum or maximum in each dimension
1441  // to enlarge the bounding box.
1442  // To do this, for each, dimension, we take the derivative
1443  // of the cubic, leaving a quadratic, and find the zeros of it.
1444  // The quadratic is such that its ith derivatives at zero are
1445  // the (i+1)th derivatives of the curve segment at zero.
1446 
1447  // 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
1448  // 1st derivative is (-(1/2)(1-t)^2)*P0 + (-1 + (1-t)^2)*P1 + (1 - (1/2)(1-t)^2)*P2
1449  // 2nd derivative is (1-t)*P0 + (-2(1-t))*P1 + (1-t)*P2
1450  // 3rd derivative is (-1)*P0 + (2)*P1 + (-1)*P2
1451 
1452  // a = (1/2) * 3rd derivative of curve segment at zero
1453  // b = 2nd derivative of curve segment at zero
1454  // (this is equivalent to the 2nd difference)
1455  UT_Vector3 b = cvs[0] - 2.0f*cvs[1] + cvs[2];
1456  UT_Vector3 a = -0.5f*b;
1457 
1458  // c = 1st derivative of curve segment at zero
1459  // (this is equivalent to the central difference)
1460  UT_Vector3 c = cvs[2] - cvs[0];
1461  c *= 0.5f;
1462 
1463  enlargeBoundingBoxCommon<UT_SplineCubic::evalSubDEnd<float> >(box, cvs, a, b, c, rootmin, rootmax);
1464 }
1465 
1466 #endif
static UT_Matrix4 getClosedWeightsTranspose(int seg, int nseg, bool deriv=false)
Definition: UT_Spline.h:538
#define SYSmax(a, b)
Definition: SYS_Math.h:1365
UT_SPLINE_BASIS getGlobalBasis() const
Query the basis or knot length of the spline.
Definition: UT_Spline.h:713
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:715
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:12
void setGlobalBasis(UT_SPLINE_BASIS b)
Definition: UT_Spline.h:718
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:385
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:716
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:959
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:100
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:807
int getVectorSize() const
Definition: UT_Spline.h:714
double fpreal64
Definition: SYS_Types.h:185
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
Definition: VM_SIMD.h:126
static void enlargeBoundingBoxOpen(UT_BoundingBox &box, const UT_Vector3 *cvs, float rootmin, float rootmax)
Definition: UT_Spline.h:1385
#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:1409
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:263
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:1438
static T evalOpen(const T *cvs, float t)
Definition: UT_Spline.h:55
GLint GLenum GLint x
Definition: glcorearb.h:408
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:258
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)
UT_Vector2T< T > colVecMult(const UT_Matrix2T< S > &m, const UT_Vector2T< T > &v)
Definition: UT_Vector2.h:481
#define SYSmin(a, b)
Definition: SYS_Math.h:1366
static const UT_Matrix4 theSubDFirstDerivBasis
Definition: UT_Spline.h:358
float fpreal32
Definition: SYS_Types.h:184