HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GU_CurveFrame.C
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2024
3  * Side Effects Software Inc. All rights reserved.
4  *
5  * Redistribution and use of Houdini Development Kit samples in source and
6  * binary forms, with or without modification, are permitted provided that the
7  * following conditions are met:
8  * 1. Redistributions of source code must retain the above copyright notice,
9  * this list of conditions and the following disclaimer.
10  * 2. The name of Side Effects Software may not be used to endorse or
11  * promote products derived from this software without specific prior
12  * written permission.
13  *
14  * THIS SOFTWARE IS PROVIDED BY SIDE EFFECTS SOFTWARE `AS IS' AND ANY EXPRESS
15  * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
17  * NO EVENT SHALL SIDE EFFECTS SOFTWARE BE LIABLE FOR ANY DIRECT, INDIRECT,
18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
19  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
20  * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
21  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
22  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
23  * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  *
25  *----------------------------------------------------------------------------
26  * Definitions of functions for computing reference frames
27  * for curve vertices.
28  */
29 
30 #include "GU_CurveFrame.h"
31 
32 #include <GEO/GEO_Curve.h>
34 #include <GA/GA_Iterator.h>
35 #include <GA/GA_OffsetList.h>
36 #include <GA/GA_ElementGroup.h>
37 #include <GA/GA_Range.h>
38 #include <GA/GA_SplittableRange.h>
39 #include <GA/GA_Types.h>
40 #include <UT/UT_Assert.h>
41 #include <UT/UT_Interrupt.h>
42 #include <UT/UT_Matrix3.h>
43 #include <UT/UT_Matrix4.h>
44 #include <UT/UT_ParallelUtil.h>
45 #include <UT/UT_Ramp.h>
46 #include <UT/UT_SmallArray.h>
47 #include <UT/UT_StringHolder.h>
48 #include <UT/UT_Vector3.h>
49 #include <SYS/SYS_Math.h>
50 #include <SYS/SYS_Types.h>
51 
52 using namespace UT::Literal;
53 
54 namespace HDK_Sample {
55 
56 namespace GU_CurveFrame {
57 
58 static constexpr float theExtremelySmallLength2 = 1e-37f;
59 static constexpr float theQuiteSmallRelativeLength = 1e-5f;
60 static constexpr float theQuiteSmallRelativeLength2 = 1e-10f;
61 
62 template<typename T>
63 static void
64 interpolateTangent(
65  const TangentType tangent_type,
66  const bool stretch_using_backbone,
67  const T max_stretch_scale,
68  const T max_stretch_length_threshold, // 2/max_stretch_length
69  const UT_Vector3T<T> &prev_edge,
70  const UT_Vector3T<T> &next_edge,
71  const T prev_length2,
72  const T next_length2,
73  const bool is_last, // false for first, true for last
74  UT_Vector3T<T> *tangent,
75  UT_Vector3T<T> *pend_stretch_dir,
76  T *pend_stretch_scale,
77  bool normalize)
78 {
79  UT_ASSERT_MSG_P(prev_length2 > 0 && next_length2 > 0, "This function doesn't handle zero-length edges. The caller should.");
80  UT_Vector3T<T> next_dir = next_edge;
81  if (next_length2 != 1) {
82  next_dir /= SYSsqrt(next_length2);
83  }
84  // first_dir and next_dir are unit vectors.
85  // We want to project back what prev_dir would be
86  // and compute the tangent from first_dir
87  // and prev_dir.
88  // Projection back the same angle as
89  // the angle between first_dir and next_dir.
90  UT_Vector3T<T> prev_dir = prev_edge;
91  if (prev_length2 != 1) {
92  prev_dir /= SYSsqrt(prev_length2);
93  }
94  UT_Vector3T<T> mid_dir;
95  T mid_length;
96  if (tangent_type == TangentType::CIRCULAR || stretch_using_backbone)
97  {
98  mid_dir = prev_dir + next_dir;
99  T mid_length2 = mid_dir.length2();
100  mid_length = SYSsqrt(mid_length2);
101  if (mid_length < theQuiteSmallRelativeLength) {
102  // Almost perfectly backtracking, so just pick the inner one.
103  mid_dir = is_last ? prev_dir : next_dir;
104  }
105  else {
106  mid_dir /= mid_length;
107  }
108  }
109  if (stretch_using_backbone) {
110  T stretch_scale;
111  UT_Vector3T<T> stretch_dir;
112  if (mid_length < max_stretch_length_threshold) {
113  stretch_scale = max_stretch_scale;
114 
115  if (mid_length < theQuiteSmallRelativeLength) {
116  // Almost perfectly backtracking, so just pick the inner one.
117  stretch_dir = is_last ? prev_dir : next_dir;
118  }
119  else {
120  stretch_dir = next_dir - prev_dir;
121  stretch_dir.normalize();
122  }
123  }
124  else {
125  // It takes a bit of trig to confirm, but the
126  // stretch scale is 1/cos(angle/2), which is 2/L_mid.
127  stretch_scale = SYSmin(2.0f / mid_length, max_stretch_scale);
128  stretch_dir = next_dir - prev_dir;
129  stretch_dir.normalize();
130  }
131  *pend_stretch_dir = stretch_dir;
132  *pend_stretch_scale = stretch_scale;
133  }
134  if (!tangent)
135  return;
136 
137  switch (tangent_type) {
139  *tangent = mid_dir;
140  break;
141  case TangentType::PREV:
142  *tangent = normalize ? prev_dir : prev_edge;
143  break;
144  case TangentType::NEXT:
145  *tangent = normalize ? next_dir : next_edge;
146  break;
147  case TangentType::SUBD:
148  // We need to use the unnormalized directions to interpolate in the subd case,
149  // even if we're normalizing after.
150  mid_dir = prev_edge + next_edge;
151  if (normalize)
152  {
153  T mid_length2 = mid_dir.length2();
154  if (mid_length2 < theQuiteSmallRelativeLength2*next_length2) {
155  // Almost perfectly stopping, so just pick the inner one (normalized).
156  mid_dir = is_last ? prev_dir : next_dir;
157  }
158  else {
159  mid_dir /= SYSsqrt(mid_length2);
160  }
161  }
162  else
163  mid_dir *= T(0.5);
164  *tangent = mid_dir;
165  break;
166  case TangentType::NONE:
167  UT_ASSERT_MSG(0, "This case should have been excluded earlier.");
168  break;
169  }
170 }
171 
172 template<typename T>
174 extrapolateVectorCircular(const UT_Vector3T<T> &a, const UT_Vector3T<T> &b) {
175  // If 'b' is a unit vector, this is a reflection of a through the line that goes through 'b' and the origin,
176  // resulting in a vector that's 'a' rotated by twice the rotation that takes 'a' to 'b'.
177  return (2*dot(a,b))*b - a;
178 }
179 
180 template<typename T>
182 extrapolateVectorLinear(const UT_Vector3T<T> &a, const UT_Vector3T<T> &b) {
183  // This is a linear extrapolation as far from 'b' as 'b' is from 'a'.
184  return b + (b-a);
185 }
186 
187 template<typename T>
188 static void
189 extrapolateEndTangent(
190  const TangentType tangent_type,
191  const bool stretch_using_backbone,
192  const T max_stretch_scale,
193  const T max_stretch_length_threshold, // 2/max_stretch_length
194  const UT_Vector3T<T> &outer_edge,
195  const UT_Vector3T<T> &inner_edge,
196  const UT_Vector3T<T> &outer_dir, // outer_edge normalized
197  const bool is_last, // false for first, true for last
198  UT_Vector3T<T> *end_tangent,
199  UT_Vector3T<T> *pend_stretch_dir,
200  T *pend_stretch_scale,
201  const bool outer_edge_is_tangent)
202 {
203  UT_Vector3T<T> inner_dir = inner_edge;
204  T inner_length2 = inner_dir.length2();
205  if (inner_length2 == 0) {
206  if (end_tangent)
207  *end_tangent = outer_dir;
208  return;
209  }
210 
211  inner_dir /= SYSsqrt(inner_length2);
212 
213  UT_Vector3T<T> extrapolated_edge;
214  T extrapolated_length2;
215  UT_Vector3T<T> end_edge_or_dir;
216  T end_length2;
217  if (tangent_type != TangentType::SUBD) {
218  extrapolated_edge = extrapolateVectorCircular(inner_dir, outer_dir);
219  extrapolated_length2 = 1;
220  end_edge_or_dir = outer_edge_is_tangent ? inner_dir : outer_dir;
221  end_length2 = 1;
222  }
223  else {
224  // We need to use the unnormalized directions to extrapolate in the subd case.
225  extrapolated_edge = extrapolateVectorLinear(inner_edge, outer_edge);
226  extrapolated_length2 = extrapolated_edge.length2();
227  end_edge_or_dir = outer_edge_is_tangent ? inner_edge : outer_edge;
228  end_length2 = outer_edge_is_tangent ? inner_length2 : outer_edge.length2();
229  }
230 
231  interpolateTangent(tangent_type,
232  stretch_using_backbone,
233  max_stretch_scale,
234  max_stretch_length_threshold,
235  is_last ? end_edge_or_dir : extrapolated_edge,
236  is_last ? extrapolated_edge : end_edge_or_dir,
237  is_last ? end_length2 : extrapolated_length2,
238  is_last ? extrapolated_length2 : end_length2,
239  is_last,
240  end_tangent, pend_stretch_dir, pend_stretch_scale,
241  true); // FIXME: Support non-normalized scales!!!
242 }
243 
244 template<typename T,typename T2>
245 static void
246 computeSingleBackboneFrames(
247  bool &rotate_using_backbone,
248  bool &stretch_using_backbone,
249  bool need_directions,
250  UT_Array<UT_Vector3T<T>> &directions,
251  UT_Array<UT_Vector3T<T>> &tangents,
252  UT_Array<UT_Vector3T<T>> &up_vectors,
253  UT_Array<UT_Vector3T<T>> &stretch_dirs,
254  UT_Array<T> &stretch_scales,
255  T &total_twist_around_loop,
256  const GEO_Detail *geo,
257  GA_Size nedges,
258  GA_Size npoints,
259  GA_Size nverts,
260  bool closed,
261  TangentType tangent_type,
262  bool extrapolate_end_tangents,
263  const GA_ROHandleV3D &instance_N,
264  const GA_Offset primoff,
265  const GA_OffsetListRef &vertices,
266  const bool use_nurbs_tangents,
267  T max_stretch_scale,
268  T max_stretch_length_threshold,
269  const GA_ROHandleV3D &instance_up,
270  UT_Vector3T<T> target_up_vector,
271  const bool use_normal_vector_up,
272  const bool target_up_vector_at_start,
273  const bool continuous_closed_curves,
274  const bool use_end_target_up_vector,
275  UT_Vector3T<T> end_target_up_vector,
276  RotationPer twist_per,
277  const GA_ROHandleT<T2> &roll_attrib,
278  const int ucomponent)
279 {
280  rotate_using_backbone = (rotate_using_backbone && nedges >= 1);
281  stretch_using_backbone = (stretch_using_backbone && npoints >= 3);
282 
283  if (!rotate_using_backbone && !stretch_using_backbone && !need_directions)
284  return;
285 
286  // Compute the backbone frames for all vertices of the curve,
287  // as well as the stretch vectors, if applicable.
288 
289  directions.setSize(nedges);
290 
291  const UT_Vector3T<T> pos0 = geo->getPos3(geo->vertexPoint(vertices(0)));
292  UT_Vector3T<T> prev_pos = pos0;
293  T max_length_squared = 0;
294  for (GA_Size i = 1; i < nverts; ++i)
295  {
296  const UT_Vector3T<T> cur_pos = geo->getPos3(geo->vertexPoint(vertices(i)));
297  UT_Vector3T<T> dir = cur_pos - prev_pos;
298  T length2 = dir.length2();
299  // Force very short edges to exactly zero, so that length2()
300  // will always produce exactly zero for them, even if the
301  // compiler pulls some inconsistent shenanigans.
302  if (length2 < theExtremelySmallLength2)
303  {
304  length2 = 0;
305  dir.assign(0,0,0);
306  }
307  directions(i-1) = dir;
308  max_length_squared = SYSmax(max_length_squared, length2);
309  prev_pos = cur_pos;
310  }
311  if (closed && nverts == nedges)
312  {
313  // One last edge if closed and not unrolled
314  UT_Vector3T<T> dir = pos0 - prev_pos;
315  T length2 = dir.length2();
316  // Force very short edges to exactly zero, as above.
317  if (length2 < theExtremelySmallLength2)
318  {
319  length2 = 0;
320  dir.assign(0,0,0);
321  }
322  directions(nverts-1) = dir;
323  max_length_squared = SYSmax(max_length_squared, length2);
324  }
325 
326  if (!rotate_using_backbone && !stretch_using_backbone)
327  return;
328 
329  if (max_length_squared == 0)
330  {
331  // If all edges are effectively zero-length, stick
332  // with no rotation or stretching.
333  rotate_using_backbone = false;
334  stretch_using_backbone = false;
335  return;
336  }
337 
338  // We need to specially handle any zero-length edges at the
339  // beginning or end of the curve.
340  // These loops rely on length2() being consistently
341  // zero or non-zero, because we know that there's
342  // at least one non-zero-length edge; we just don't
343  // want to incorrectly loop all the way past the end.
344  GA_Size first_nonzero = 0;
345  GA_Size last_nonzero = directions.size()-1;
346  while (first_nonzero < last_nonzero && directions(first_nonzero).length2() == 0)
347  {
348  ++first_nonzero;
349  }
350  while (first_nonzero < last_nonzero && directions(last_nonzero).length2() == 0)
351  {
352  --last_nonzero;
353  }
354 
355  tangents.setSize(npoints);
356  up_vectors.setSize(npoints);
357 
358  if (instance_N.isValid())
359  {
360  UT_ASSERT(!instance_up.isValid());
361  for (exint i = 0; i < npoints; ++i)
362  {
363  GA_Offset ptoff = geo->vertexPoint(vertices(i));
364  UT_Vector3D normal = instance_N.get(ptoff);
365  normal.normalize();
366  tangents[i] = normal;
367  }
368  }
369  if (instance_up.isValid())
370  {
371  UT_ASSERT(!instance_N.isValid());
372  for (exint i = 0; i < npoints; ++i)
373  {
374  GA_Offset ptoff = geo->vertexPoint(vertices(i));
375  UT_Vector3D up = instance_up.get(ptoff);
376  up.normalize();
377  up_vectors[i] = up;
378  }
379  }
380 
381  if (first_nonzero == last_nonzero)
382  {
383  // Exactly one non-zero-length edge, (usually just one total).
384  if (!instance_N.isValid())
385  {
386  // The single edge length must have length sqrt(max_edge_length_squared).
387  UT_Vector3T<T> dir = directions(first_nonzero) / SYSsqrt(max_length_squared);
388  tangents.constant(dir);
389 
390  if (!instance_up.isValid())
391  {
392  // No N or up, so same up for all points
393  UT_Vector3T<T> other = target_up_vector;
394  other -= dot(other,dir)*dir;
395  T length2 = other.length2();
396  if (length2 >= theQuiteSmallRelativeLength2)
397  other /= SYSsqrt(length2);
398  else
399  other.arbitraryPerp(dir);
400  up_vectors.constant(other);
401  }
402  else
403  {
404  // up attribute, so get closest to up perpendicular to dir
405  for (exint i = 0; i < npoints; ++i)
406  {
407  UT_Vector3T<T> other = up_vectors[i];
408  other -= dot(other,dir)*dir;
409  T length2 = other.length2();
410  if (length2 >= theQuiteSmallRelativeLength2)
411  other /= SYSsqrt(length2);
412  else
413  other.arbitraryPerp(dir);
414  up_vectors[i] = other;
415  }
416  }
417  }
418  else
419  {
420  // N attribute, but no up, so get closest to target up perpendicular to N
421  UT_ASSERT(!instance_up.isValid());
422 
423  for (exint i = 0; i < npoints; ++i)
424  {
425  UT_Vector3D dir = tangents[i];
426  UT_Vector3D other = target_up_vector;
427  other -= dot(other,dir)*dir;
428  T length2 = other.length2();
429  if (length2 >= theQuiteSmallRelativeLength2)
430  other /= SYSsqrt(length2);
431  else
432  other.arbitraryPerp(dir);
433  up_vectors[i] = other;
434  }
435  }
436 
437  // No stretching if only one non-zero-length edge.
438  stretch_using_backbone = false;
439 
440  // Nothing more to do in this case.
441  return;
442  }
443 
444  // More than one non-zero-length edge.
445 
446  if (stretch_using_backbone)
447  {
448  stretch_dirs.setSizeNoInit(npoints);
449  // setSize would initialize to zero, but let's make it explicitly clear
450  // that we want to initialize these to a default of zero.
451  // A stretch vector being zero will result in no stretch.
452  stretch_dirs.constant(UT_Vector3T<T>(0,0,0));
453  stretch_scales.setSizeNoInit(npoints);
454  stretch_scales.constant(T(1));
455  }
456 
457  if (instance_N.isValid() && stretch_using_backbone)
458  {
459  // If instance_N is valid, don't compute tangents, but do still
460  // compute stretch_dirs and stretch_scales if stretch_using_backbone,
461  // including extrapolating ends as needed.
462  if (!closed)
463  {
464  // Always act as if extrapolating, but based on tangent,
465  // since end tangents may not line up with end directions.
466  if (first_nonzero == 0)
467  {
468  extrapolateEndTangent(tangent_type,
469  stretch_using_backbone, max_stretch_scale,
470  max_stretch_length_threshold,
471  tangents[0],
472  directions[0],
473  tangents[0], // Already normalized
474  false, // Not last edge (since first edge)
475  (UT_Vector3T<T>*)nullptr,
476  &stretch_dirs(0),
477  &stretch_scales(0),
478  true);// outer_edge is a tangent
479  }
480  if (last_nonzero == directions.size()-1)
481  {
482  extrapolateEndTangent(tangent_type,
483  stretch_using_backbone, max_stretch_scale,
484  max_stretch_length_threshold,
485  tangents.last(),
486  directions.last(),
487  tangents.last(), // Already normalized
488  true, // Last edge
489  (UT_Vector3T<T>*)nullptr,
490  &stretch_dirs(npoints-1),
491  &stretch_scales(npoints-1),
492  true);// outer_edge is a tangent
493  }
494  }
495  else if (first_nonzero == 0 && last_nonzero == directions.size()-1)
496  {
497  // Closed backbone and no zero-length edges at beginning or end, so
498  // it's like a normal case, except with non-adjacent edge indices.
499 
500  interpolateTangent(
501  tangent_type,
502  stretch_using_backbone,
503  max_stretch_scale,
504  max_stretch_length_threshold,
505  directions.last(),
506  directions(0),
507  directions.last().length2(),
508  directions(0).length2(),
509  false, // Arbitrarily choose not last
510  (UT_Vector3T<T>*)nullptr,
511  &stretch_dirs(0),
512  &stretch_scales(0),
513  true); // FIXME: Support non-normalized scales!!!
514  }
515 
516  // Middle stretch cases
517  GA_Size prev_nonzero = first_nonzero;
518  T prev_length2 = directions[first_nonzero].length2();
519  for (GA_Size next_edgei = first_nonzero+1; next_edgei <= last_nonzero; ++next_edgei) {
520  T next_length2 = directions[next_edgei].length2();
521  if (next_length2 != 0)
522  {
523  // Common case: Two adjacent, non-zero-length edges.
524  interpolateTangent(
525  tangent_type,
526  stretch_using_backbone,
527  max_stretch_scale,
528  max_stretch_length_threshold,
529  directions[prev_nonzero],
530  directions[next_edgei],
531  prev_length2,
532  next_length2,
533  false, // Arbitrarily choose not last
534  (UT_Vector3T<T>*)nullptr,
535  &stretch_dirs[next_edgei],
536  &stretch_scales[next_edgei],
537  true); // FIXME: Support non-normalized scales!!!
538  }
539  else
540  {
541  // Find next non-zero-length edge.
542  do {
543  ++next_edgei;
544  next_length2 = directions(next_edgei).length2();
545  } while (next_edgei != last_nonzero && next_length2 == 0);
546  }
547  prev_nonzero = next_edgei;
548  prev_length2 = next_length2;
549  }
550  }
551 
552  if (!instance_N.isValid() && use_nurbs_tangents)
553  {
554  int primtype = geo->getPrimitiveTypeId(primoff);
555  if (primtype == GA_PRIMNURBCURVE || primtype == GA_PRIMBEZCURVE)
556  {
557  const GEO_Curve *curve = UTverify_cast<const GEO_Curve *>(geo->getPrimitive(primoff));
558  const bool curve_closed = curve->isClosed();
559  const GA_Basis *basis = curve->getBasis();
560 
561  UT_Array<float> greville_us;
562  greville_us.setSizeNoInit(nverts);
563  for (exint i = 0; i < nverts; ++i)
564  {
565  float u = basis->getGreville(i, true, curve_closed);
566  greville_us(i) = u;
567  }
568 
569  // Request tangents from curve.
570  for (exint i = 0; i < nverts; ++i)
571  {
572  UT_Vector4 cur_tangent;
573  float u = greville_us(i);
574  curve->evaluate(u, cur_tangent, /*du=*/1);
575 
576 
577 
578 
579  }
580  }
581  }
582 
583  // Handle the beginning and end of the backbone first,
584  // so that everything else is just a common middle case.
585  UT_Vector3T<T> first_dir = directions(first_nonzero);
586  UT_Vector3T<T> last_dir = directions(last_nonzero);
587  first_dir.normalize();
588  last_dir.normalize();
589  if (!instance_N.isValid())
590  {
591  if (!closed) {
592  if (first_nonzero > 0) {
593  // If there are zero-length edges at the beginning,
594  // set tangents to the first non-zero-length edge direction,
595  // up to and including the first point of that edge.
596  for (GA_Size i = 0; i <= first_nonzero; ++i) {
597  tangents(i) = first_dir;
598  }
599  }
600  else if (extrapolate_end_tangents) {
601  extrapolateEndTangent(tangent_type,
602  stretch_using_backbone, max_stretch_scale,
603  max_stretch_length_threshold,
604  directions(0),
605  directions(1), // We know there are at least two, so this is safe.
606  first_dir,
607  false, // Not last edge (since first edge)
608  &tangents(0),
609  stretch_using_backbone ? &stretch_dirs(0) : nullptr,
610  stretch_using_backbone ? &stretch_scales(0) : nullptr,
611  false);// outer_edge is not a tangent
612  }
613  else {
614  // Not extrapolating end tangents, so just copy.
615  tangents(0) = first_dir;
616  }
617 
618  if (last_nonzero < directions.size()-1) {
619  // Same as above for zero-length edges at the end.
620  for (GA_Size i = last_nonzero+1; i < npoints; ++i) {
621  tangents(i) = last_dir;
622  }
623  }
624  else if (extrapolate_end_tangents) {
625  extrapolateEndTangent(tangent_type,
626  stretch_using_backbone, max_stretch_scale,
627  max_stretch_length_threshold,
628  directions.last(),
629  directions(directions.size()-2), // We know there are at least two, so this is safe.
630  last_dir,
631  true, // Last edge
632  &tangents(npoints-1),
633  stretch_using_backbone ? &stretch_dirs(npoints-1) : nullptr,
634  stretch_using_backbone ? &stretch_scales(npoints-1) : nullptr,
635  false);// outer_edge is not a tangent
636  }
637  else {
638  // Not extrapolating end tangents, so just copy.
639  tangents(npoints-1) = last_dir;
640  }
641  }
642  else if (first_nonzero > 0 || last_nonzero < directions.size()-1) {
643  // Closed backbone and at least one zero-length edge at beginning or end.
644 
645  switch (tangent_type) {
647  {
648  if (first_nonzero == 1 && last_nonzero == directions.size()-1) {
649  tangents(0) = last_dir;
650  tangents(1) = first_dir;
651  }
652  else if (first_nonzero == 0 && last_nonzero == directions.size()-2) {
653  tangents.last() = last_dir;
654  tangents(0) = first_dir;
655  }
656  else {
657  // 3+ coincident, so circularly interpolate tangent
658  T cos_angle = dot(first_dir,last_dir);
659  UT_Vector3T<T> sin_axis = first_dir - cos_angle*last_dir;
660  sin_axis.normalize();
661  // TODO: atan2(sin_angle,cos_angle) is more stable for small angles or half turns.
662  T full_angle = SYSacos(cos_angle);
663  GA_Size nsteps = first_nonzero + (directions.size()-1 - last_nonzero);
664  GA_Size j = 0;
665  for (GA_Size i = last_nonzero+1; i < npoints; ++i, ++j) {
666  T part_angle = j*full_angle/nsteps;
667  tangents(i) = SYScos(part_angle)*last_dir + SYSsin(part_angle)*sin_axis;
668  }
669  for (GA_Size i = 0; i <= first_nonzero; ++i, ++j) {
670  T part_angle = j*full_angle/nsteps;
671  tangents(i) = SYScos(part_angle)*last_dir + SYSsin(part_angle)*sin_axis;
672  }
673  }
674  break;
675  }
676  case TangentType::SUBD:
677  {
678  if (first_nonzero == 1 && last_nonzero == directions.size()-1) {
679  tangents(0) = last_dir;
680  tangents(1) = first_dir;
681  }
682  else if (first_nonzero == 0 && last_nonzero == directions.size()-2) {
683  tangents.last() = last_dir;
684  tangents(0) = first_dir;
685  }
686  else {
687  // 3+ coincident, so linearly interpolate tangent
688  GA_Size nsteps = first_nonzero + (directions.size()-1 - last_nonzero);
689  GA_Size j = 0;
690  // NOTE: This value 0,0,1 should never end up being used. It should be overwritten first.
691  UT_Vector3T<T> prev_tangent(0,0,1);
692  for (GA_Size i = last_nonzero+1; i < npoints; ++i, ++j) {
693  UT_Vector3T<T> tangent = SYSlerp(directions(last_nonzero), directions(first_nonzero), T(j)/nsteps);
694  tangent.normalize();
695  if (tangent.length2() < 0.5f) {
696  // If didn't normalize, fall back to previous tangent,
697  // (shouldn't happen on first iteration, since tangent should be last_dir).
698  tangent = prev_tangent;
699  }
700  else {
701  prev_tangent = tangent;
702  }
703  tangents(i) = tangent;
704  }
705  for (GA_Size i = 0; i <= first_nonzero; ++i, ++j) {
706  UT_Vector3T<T> tangent = SYSlerp(directions(last_nonzero), directions(first_nonzero), T(j)/nsteps);
707  tangent.normalize();
708  if (tangent.length2() < 0.5f) {
709  // If didn't normalize, fall back to previous tangent,
710  // (shouldn't happen on first iteration, since tangent should be last_dir).
711  tangent = prev_tangent;
712  }
713  else {
714  prev_tangent = tangent;
715  }
716  tangents(i) = tangent;
717  }
718  }
719  break;
720  }
721  case TangentType::PREV:
722  case TangentType::NEXT:
723  {
724  UT_Vector3T<T> tangent = (tangent_type == TangentType::PREV) ? last_dir : first_dir;
725  for (GA_Size i = last_nonzero+1; i < npoints; ++i) {
726  tangents(i) = tangent;
727  }
728  for (GA_Size i = 0; i <= first_nonzero; ++i) {
729  tangents(i) = tangent;
730  }
731  break;
732  }
733  case TangentType::NONE:
734  UT_ASSERT_MSG(0, "This case should have been excluded earlier.");
735  break;
736  }
737  }
738  else {
739  // Closed backbone and no zero-length edges at beginning or end, so
740  // it's like a normal case, except with non-adjacent edge indices.
741 
742  interpolateTangent(
743  tangent_type,
744  stretch_using_backbone,
745  max_stretch_scale,
746  max_stretch_length_threshold,
747  directions.last(),
748  directions(0),
749  directions.last().length2(),
750  directions(0).length2(),
751  false, // Arbitrarily choose not last
752  &tangents(0),
753  stretch_using_backbone ? &stretch_dirs(0) : nullptr,
754  stretch_using_backbone ? &stretch_scales(0) : nullptr,
755  true); // FIXME: Support non-normalized scales!!!
756  }
757 
758  // The beginning and end have been handled, and we have two non-zero-length
759  // edges bounding the remaining edges.
760 
761 
762  GA_Size prev_nonzero = first_nonzero;
763  T prev_length2 = directions(first_nonzero).length2();
764  for (GA_Size next_edgei = first_nonzero+1; next_edgei <= last_nonzero; ++next_edgei) {
765  T next_length2 = directions(next_edgei).length2();
766  if (next_length2 != 0) {
767  // Common case: Two adjacent, non-zero-length edges.
768  interpolateTangent(
769  tangent_type,
770  stretch_using_backbone,
771  max_stretch_scale,
772  max_stretch_length_threshold,
773  directions(prev_nonzero),
774  directions(next_edgei),
775  prev_length2,
776  next_length2,
777  false, // Arbitrarily choose not last
778  &tangents(next_edgei),
779  stretch_using_backbone ? &stretch_dirs(next_edgei) : nullptr,
780  stretch_using_backbone ? &stretch_scales(next_edgei) : nullptr,
781  true); // FIXME: Support non-normalized scales!!!
782  prev_nonzero = next_edgei;
783  prev_length2 = next_length2;
784  continue;
785  }
786 
787  // Less-common case: zero-length edge(s)
788  // Find next non-zero-length edge.
789  do {
790  ++next_edgei;
791  next_length2 = directions(next_edgei).length2();
792  } while (next_edgei != last_nonzero && next_length2 == 0);
793 
794  GA_Size nsteps = next_edgei - prev_nonzero - 1;
795 
796  UT_Vector3T<T> prev_dir = directions(prev_nonzero) / SYSsqrt(prev_length2);
797  UT_Vector3T<T> next_dir = directions(next_edgei) / SYSsqrt(next_length2);
798  switch (tangent_type) {
800  {
801  if (nsteps == 1) {
802  tangents(next_edgei-1) = prev_dir;
803  tangents(next_edgei) = next_dir;
804  }
805  else {
806  // 3+ coincident, so circularly interpolate tangent
807  T cos_angle = dot(next_dir,prev_dir);
808  UT_Vector3T<T> sin_axis = next_dir - cos_angle*prev_dir;
809  sin_axis.normalize();
810  // TODO: atan2(sin_angle,cos_angle) is more stable for small angles or half turns.
811  T full_angle = SYSacos(cos_angle);
812  GA_Size j = 0;
813  for (GA_Size i = prev_nonzero+1; i <= next_edgei; ++i, ++j) {
814  T part_angle = j*full_angle/nsteps;
815  tangents(i) = SYScos(part_angle)*prev_dir + SYSsin(part_angle)*sin_axis;
816  }
817  }
818  break;
819  }
820  case TangentType::SUBD:
821  {
822  if (nsteps == 1) {
823  tangents(next_edgei-1) = prev_dir;
824  tangents(next_edgei) = next_dir;
825  }
826  else {
827  // 3+ coincident, so linearly interpolate tangent
828  GA_Size j = 0;
829  // NOTE: This value 0,0,1 should never end up being used. It should be overwritten first.
830  UT_Vector3T<T> prev_tangent(0,0,1);
831  for (GA_Size i = prev_nonzero+1; i <= next_edgei; ++i, ++j) {
832  UT_Vector3T<T> tangent = SYSlerp(directions(prev_nonzero), directions(next_edgei), T(j)/nsteps);
833  tangent.normalize();
834  if (tangent.length2() < 0.5f) {
835  // If didn't normalize, fall back to previous tangent,
836  // (shouldn't happen on first iteration, since tangent should be last_dir).
837  tangent = prev_tangent;
838  }
839  else {
840  prev_tangent = tangent;
841  }
842  tangents(i) = tangent;
843  }
844  }
845  break;
846  }
847  case TangentType::PREV:
848  case TangentType::NEXT:
849  {
850  UT_Vector3T<T> tangent = (tangent_type == TangentType::PREV) ? prev_dir : next_dir;
851  for (GA_Size i = prev_nonzero+1; i <= next_edgei; ++i) {
852  tangents(i) = tangent;
853  }
854  break;
855  }
856  case TangentType::NONE:
857  UT_ASSERT_MSG(0, "This case should have been excluded earlier.");
858  break;
859  }
860 
861  prev_nonzero = next_edgei;
862  prev_length2 = next_length2;
863  }
864  }
865 
866  // We now have all of the tangents.
867  // stretch_dirs and stretch_scales have also been initialized
868  // if stretch_using_backbone is true.
869 
870  if (instance_up.isValid())
871  {
872  // up attribute, but we still need to force up_vectors
873  // to be perpendicular to tangents.
874 
875  for (exint i = 0; i < npoints; ++i)
876  {
877  UT_Vector3D dir = tangents[i];
878  UT_Vector3D other = up_vectors[i];
879  other -= dot(other,dir)*dir;
880  T length2 = other.length2();
881  if (length2 >= theQuiteSmallRelativeLength2)
882  other /= SYSsqrt(length2);
883  else
884  other.arbitraryPerp(dir);
885  up_vectors[i] = other;
886  }
887 
888  return;
889  }
890 
891  if (use_normal_vector_up) {
892  // Compute a face normal for target_up_vector, if possible.
893  // This isn't the most efficient way to compute a normal, but we want
894  // to try to fall back on the largest normal along the way, so we
895  // go one triangle at a time. We know that we have at least two
896  // non-zero-length edges here, so there is at least one triangle,
897  // even though it might have zero area.
898 
899  // Double-precision for the sum to reduce the risk of catastrophic
900  // roundoff error, (e.g. if a curve had more than 16,777,216 vertices,
901  // everything after that would effectively be lost in single-precision.)
902  UT_Vector3D normal_sum(0,0,0);
903  UT_Vector3D max_normal;
904  double length2_max_normal = 0;
905  T length2_max_chord = 0;
906  for (GA_Size i = 1; i < nverts-1; ++i) {
907  UT_Vector3T<T> posi = geo->getPos3(geo->vertexPoint(vertices(i)));
908  UT_Vector3T<T> chord = posi-pos0;
909 
910  // The length of the cross product is 2x the area of the triangle
911  // from 0 to i to i+1, but with the direction being the triangle's
912  // normal. When summed, they give a vector whose length is the
913  // area of the polygon projected into a plane perpendicular to the
914  // direction of the vector.
915  normal_sum += cross(directions(i), chord);
916 
917  // Record the longest normal along the way, for a fallback.
918  double length2 = normal_sum.length2();
919  if (length2 > length2_max_normal) {
920  length2_max_normal = length2;
921  max_normal = normal_sum;
922  }
923  length2_max_chord = SYSmax(length2_max_chord, chord.length2());
924  }
925  UT_Vector3T<T> posi = geo->getPos3(geo->vertexPoint(vertices(nverts-1)));
926  length2_max_chord = SYSmax(length2_max_chord, (posi-pos0).length2());
927 
928  // Double-precision to reduce risk of underflow/overflow below
929  double small_rel_length2 = theQuiteSmallRelativeLength2*length2_max_chord;
930  double small_rel_length4 = small_rel_length2*small_rel_length2;
931 
932  // length2_max_normal is the *square* of 2x the largest *area* along the way,
933  // so compare against length-squared *squared*.
934  if (length2_max_normal < small_rel_length4) {
935  // Even the max length normal along the way was quite small,
936  // which happens when the points are all along a line,
937  // so stick with default target_up_vector.
938  }
939  else {
940  double length2_normal = normal_sum.length2();
941  if (length2_normal < small_rel_length4) {
942  // normal_sum ended up being small, but it was large enough
943  // partway through, so use that. This can happen if a
944  // curve backtracks along itself to the beginning.
945  target_up_vector = max_normal / SYSsqrt(length2_max_normal);
946  }
947  else {
948  // The normal's good, so use that.
949  target_up_vector = normal_sum / SYSsqrt(length2_normal);
950  }
951  }
952  }
953 
954  // We have all tangents, but we need to compute consistent up vectors.
955  // First, we pick an arbitrary starting bitangent; it must be
956  // perpendicular to tangents(0) and it might as well be as close to
957  // perpendicular to target_up_vector as it can be.
958  UT_Vector3T<T> tangent0 = tangents(0);
959  UT_Vector3T<T> up_vectors0 = target_up_vector;
960  up_vectors0 -= dot(up_vectors0,tangent0)*tangent0;
961  T length2 = up_vectors0.length2();
962  if (length2 >= theQuiteSmallRelativeLength2) {
963  up_vectors0 /= SYSsqrt(length2);
964  }
965  else {
966  up_vectors0.arbitraryPerp(tangent0);
967  }
968  up_vectors(0) = up_vectors0;
969 
970  for (GA_Size i = 1; i < npoints; ++i) {
971  // Compute the minimal rotation matrix that takes tangents(i-1) to tangents(i).
972  // Vectors in tangents are already normalized, so no need to normalize inside (false).
973  // If this is a bottleneck, it's possible to compute the resulting vector
974  // without constructing the matrix itself, via:
975  // <v|(I-2|a><a|)(I-2|c><c|), where c is (a+b) normalized, and a and b are the tangents.
976  // v -= (2*dot(v,a))*a;
977  // v -= (2*dot(v,c))*c;
979  bool failure = rotation.dihedral(tangents(i-1), tangents(i), false);
980  if (failure) {
981  // The tangent flipped a half turn, so flip the bitangent a half turn,
982  // to keep the up vector the same.
983  up_vectors(i) = -up_vectors(i-1);
984  }
985  else {
986  // Rotate bitangent by the rotation matrix.
987  up_vectors(i) = up_vectors(i-1) * rotation;
988  }
989 
990  // Just to prevent roundoff error from getting bad, force up_vectors(i) to be orthogonal to tangents(i).
991  up_vectors(i) -= dot(up_vectors(i),tangents(i))*tangents(i);
992  // It seems implausible that up_vectors(i) could possibly be problematically small
993  // after the subtraction, since it should already have been approximately
994  // orthogonal to tangents(i) and unit length, so let's take the risk and just normalize.
995  up_vectors(i).normalize();
996  }
997 
998  // Closed backbone curves need an adjustment to make sure that beginning and end
999  // up_vectors line up, without a sudden roll rotation. In other words,
1000  // smooth the excess roll out over the whole curve.
1001  // Open curves with an end target up vector need the same type of adjustment.
1002  if ((closed && continuous_closed_curves) || use_end_target_up_vector)
1003  {
1004  UT_Vector3T<T> final_up_vector;
1005  UT_Vector3T<T> comparison_tangent;
1006  bool skip_twist = false;
1007  if (closed)
1008  {
1009  // Same as above regarding not really needing to compute this matrix
1010  // if it's a bottleneck.
1012  bool failure = rotation.dihedral(tangents[npoints-1], tangents[0], false);
1013  if (failure)
1014  {
1015  // The tangent flipped a half turn, so keep the up vector the same.
1016  // The bitangent will implicitly flip a half turn.
1017  final_up_vector = up_vectors[npoints-1];
1018  }
1019  else
1020  {
1021  // Rotate up vectors by the rotation matrix.
1022  final_up_vector = up_vectors[npoints-1] * rotation;
1023  }
1024  // Just to prevent roundoff error from getting bad, force final_up_vector to be orthogonal to tangents[0].
1025  final_up_vector -= dot(final_up_vector,tangents[0])*tangents[0];
1026  // It seems implausible that final_up_vector could possibly be problematically small
1027  // after the subtraction, since it should already have been approximately
1028  // orthogonal to tangents[0] and unit length, so let's take the risk and just normalize.
1029  final_up_vector.normalize();
1030 
1031  comparison_tangent = tangents[0];
1032  }
1033  else
1034  {
1035  // Open curve, so final up vector is already in the array.
1036  final_up_vector = up_vectors[npoints-1];
1037  comparison_tangent = tangents[npoints-1];
1038  }
1039 
1040  UT_Vector3T<T> comparison_up_vector;
1041  if (closed && continuous_closed_curves)
1042  {
1043  comparison_up_vector = up_vectors[0];
1044  }
1045  else
1046  {
1047  // Using an end target up vector.
1048  if (use_normal_vector_up)
1049  {
1050  // Normal vector was written into target_up_vector above.
1051  comparison_up_vector = target_up_vector;
1052  }
1053  else
1054  comparison_up_vector = end_target_up_vector;
1055  comparison_up_vector -= comparison_up_vector*dot(comparison_up_vector,comparison_tangent);
1056  T length2 = comparison_up_vector.length2();
1057  // If the target up vector is extremely close to the tangent, skip applying a twist.
1058  if (length2 >= theQuiteSmallRelativeLength2)
1059  comparison_up_vector /= SYSsqrt(length2);
1060  else
1061  skip_twist = true;
1062  }
1063 
1064  T net_twist = T(0);
1065  if (!skip_twist)
1066  {
1067  // NOTE: The cosine between the start and end up_vectors is
1068  // not sufficient to get the direction of the net twist.
1069  // We need to consider the direction turned relative to
1070  // the tangent. Thus, we need to use the cross product
1071  // between the two dotted against the tangent to get
1072  // the correct sine of the angle.
1073  UT_Vector3T<T> net_twist_vector = cross(final_up_vector, comparison_up_vector);
1074  T sin_net_twist = dot(net_twist_vector, comparison_tangent);
1075  T cos_net_twist = dot(final_up_vector, comparison_up_vector);
1076  net_twist = SYSatan2(sin_net_twist, cos_net_twist); // atan2(y,x) order
1077  }
1078 
1079  if (!SYSequalZero(net_twist, theQuiteSmallRelativeLength))
1080  {
1081  // Apply twist around the loop
1082  if (twist_per == RotationPer::EDGE || twist_per == RotationPer::FULLEDGES)
1083  {
1084  for (GA_Size i = 1; i < npoints; ++i)
1085  {
1086  T reverse_twist = -i*net_twist/nedges;
1087  T s = SYSsin(reverse_twist);
1088  T c = SYScos(reverse_twist);
1089  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1090  up_vectors(i) = s*bitangent + c*up_vectors(i);
1091  }
1092  }
1093  else if (twist_per == RotationPer::DISTANCE || twist_per == RotationPer::FULLDISTANCE ||
1094  (twist_per == RotationPer::ATTRIB && (
1095  !roll_attrib.isValid() ||
1096  roll_attrib->getOwner() == GA_ATTRIB_PRIMITIVE ||
1097  roll_attrib->getOwner() == GA_ATTRIB_DETAIL)))
1098  {
1099  double total_length = 0;
1100  for (GA_Size i = 0; i < nedges; ++i)
1101  {
1102  total_length += directions(i).length();
1103  }
1104 
1105  // NOTE: total_length shouldn't be zero, since we already ensured
1106  // that we have at least two non-zero-length edges.
1107  double cur_length = 0;
1108  for (GA_Size i = 1; i < npoints; ++i)
1109  {
1110  cur_length += directions(i-1).length();
1111  T reverse_twist = -cur_length*net_twist/total_length;
1112  T s = SYSsin(reverse_twist);
1113  T c = SYScos(reverse_twist);
1114  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1115  up_vectors(i) = s*bitangent + c*up_vectors(i);
1116  }
1117  }
1118  else
1119  { // twist_per == RotationPer::ATTRIB
1120  if (roll_attrib->getOwner() == GA_ATTRIB_VERTEX)
1121  {
1122 #if 0
1123  T2 ustart = 0;
1124  T2 uend = 1;
1125  if (nverts == npoints+1)
1126  {
1127  // We can really only use start and end uvs reliably
1128  // if the curve is unrolled, (1 extra vertex).
1129  // Otherwise, we default to assuming 0 to 1.
1130  ustart = roll_attrib.get(vertices(0), ucomponent);
1131  uend = roll_attrib.get(vertices.last(), ucomponent);
1132  }
1133  T2 uspan = (uend-ustart);
1134 #endif
1135 
1136  //if (uspan != T2(0))
1137  //{
1138  for (GA_Size i = 1; i < npoints; ++i)
1139  {
1140  T2 t = roll_attrib.get(vertices(i), ucomponent);
1141  //T2 t = (u-ustart)/uspan;
1142  T reverse_twist = -t*net_twist;
1143  T s = SYSsin(reverse_twist);
1144  T c = SYScos(reverse_twist);
1145  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1146  up_vectors(i) = s*bitangent + c*up_vectors(i);
1147  }
1148  //}
1149  }
1150  else
1151  {
1152  UT_ASSERT(roll_attrib->getOwner() == GA_ATTRIB_POINT);
1153  // Assume 0 to 1, since it won't end correctly and we have to guess.
1154 
1155  for (GA_Size i = 1; i < npoints; ++i)
1156  {
1157  T2 u = roll_attrib.get(geo->vertexPoint(vertices(i)), ucomponent);
1158  T reverse_twist = -u*net_twist;
1159  T s = SYSsin(reverse_twist);
1160  T c = SYScos(reverse_twist);
1161  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1162  up_vectors(i) = s*bitangent + c*up_vectors(i);
1163  }
1164  }
1165  }
1166 
1167  if (closed && continuous_closed_curves)
1168  {
1169  // The twist applied so far is the negative of the original net twist.
1170  total_twist_around_loop = -net_twist;
1171  }
1172  }
1173  }
1174 
1175  if (target_up_vector_at_start)
1176  {
1177  // Nothing else to do if target up vector is okay being at the start,
1178  // since that's what we used as our initial guess.
1179  return;
1180  }
1181 
1182  // Now that we have consistent up_vectors, we make them more stable by
1183  // making the average closer to perpendicular to the backbone polygon normal,
1184  // by introducing a global roll rotation.
1185 
1186  // This sounds pretty nebulous, but the math simplifies down nicely.
1187  // We want to find the z-axis rotation Rz such that
1188  // [0 1 0] * Rz * M
1189  // is in the direction of target_up_vector,
1190  // where M is the current average rotation matrix.
1191  // (Note that M itself likely isn't a rotation matrix, may not
1192  // have orthogonal rows, and may even be singular or near-singular.)
1193  // Depending on which direction we choose as the angle, this is:
1194  // [sin(theta) cos(theta) 0] * M
1195  // = [s c 0] * M
1196  // = s*Mx + c*My
1197  // If we solve:
1198  // [Mx.Mx Mx.My][x] [Mx.up]
1199  // [Mx.My My.My][y] = [My.up]
1200  // we'll get the vector that is least-squares closest to up
1201  // that can be made using Mx and My, as a linear combination of Mx and My.
1202  // If we normalize [x y], that gives us a valid [s c]. Then, we can apply
1203  // Rz = [ c -s 0]
1204  // [ s c 0]
1205  // [ 0 0 1]
1206  // on the left of each frame.
1207  // Of course, if Mx and My are parallel, or one of the two is much shorter
1208  // than the other, (or equivalently, if the determinant of the 2x2 matrix above
1209  // is small compared to the length-squared-squared of Mx and My), we need to handle
1210  // things differently.
1211 
1212  // Find Mx and My, (Mz isn't needed), in double-precision to avoid
1213  // catastrophic roundoff error in the accumulation.
1214  UT_Vector3D Mx(0,0,0);
1215  UT_Vector3D My(0,0,0);
1216  for (GA_Size i = 0; i < npoints; ++i) {
1217  My += up_vectors(i);
1218  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1219  Mx += bitangent;
1220  }
1221  // We don't strictly need to divide by npoints for the average,
1222  // since we'll be normalizing later, and there's little risk of overflow
1223  // using double-precision.
1224 
1225  double MxMx = dot(Mx,Mx);
1226  double MxMy = dot(Mx,My);
1227  double MyMy = dot(My,My);
1228  double Mxup = dot(Mx,UT_Vector3D(target_up_vector));
1229  double Myup = dot(My,UT_Vector3D(target_up_vector));
1230 
1231  // Part of this is effectively the same as UT_Matrix2T::solve
1232  // Side note: via awesome properties of cross & dot products,
1233  // determinant is equal to (Mx x My).(Mx x My),
1234  // so it's the length-squared of the cross product.
1235  double determinant = MxMx*MyMy - MxMy*MxMy;
1236  double scale2 = SYSmax(MxMx*MxMx, MyMy*MyMy);
1237  T s; T c;
1238  // 1e-6 is because if Mx is 1000x longer than My and they're in the ballpark
1239  // of orthogonal, My is almost certainly due to slightly imbalanced vectors
1240  // that should really cancel out exactly. If they're in the ballpark of
1241  // parallel and the non-parallel component of My is 1000x smaller than Mx,
1242  // they're effectively the same.
1243  // In either case, that yields determinant/scale2 ~ 1e-6
1244  if (!SYSequalZero(determinant, 1e-6*scale2)) {
1245  // Not near singular, so we can solve it.
1246  // We'd normally divide by determinant to solve a 2x2 system,
1247  // but we're normalizing after anyway.
1248  s = (MyMy*Mxup - MxMy*Myup);
1249  c = (MxMx*Myup - MxMy*Mxup);
1250  T normalization = s*s + c*c;
1251  if (SYSequalZero(normalization, theExtremelySmallLength2)) {
1252  s = 0;
1253  c = 1;
1254  }
1255  else {
1256  normalization = T(1) / SYSsqrt(normalization);
1257  s *= normalization;
1258  c *= normalization;
1259  }
1260  }
1261  else if (MxMx > MyMy) {
1262  s = (Mxup >= 0) ? 1 : -1;
1263  c = 0;
1264  }
1265  else {
1266  s = 0;
1267  c = (Myup >= 0) ? 1 : -1;
1268  }
1269 
1270  // New up vector =
1271  // [ c -s 0] [<--Mx-->] [cMx-sMy ]
1272  // [0 1 0] [ s c 0] [<--My-->] = [0 1 0] [sMx+cMy ] = [sMx+cMy ]
1273  // [ 0 0 1] [<--Mz-->] [<--Mz-->]
1274 
1275  for (GA_Size i = 0; i < npoints; ++i) {
1276  UT_Vector3T<T> bitangent = cross(up_vectors(i), tangents(i));
1277  up_vectors(i) = s*bitangent + c*up_vectors(i);
1278  }
1279 
1280  // Finally done computing up_vectors! Yay! :)
1281 }
1282 
1283 bool
1285  const GEO_Detail *geometry,
1286  const GA_OffsetListRef &vertices,
1287  exint &nedges,
1288  bool &closed,
1289  bool &unrolled)
1290 {
1291  if (vertices.size() == 0) {
1292  // Skip cross-sections with zero vertices
1293  return false;
1294  }
1295  if (vertices.size() == 1) {
1296  // Always treat single-vertex polygons as open,
1297  // to try to reduce confusion.
1298  closed = false;
1299  unrolled = false;
1300  nedges = 0;
1301  return true;
1302  }
1303 
1304  closed = vertices.getExtraFlag();
1305  unrolled = false;
1306  nedges = vertices.size();
1307 
1308  if (!closed) {
1309  --nedges;
1310 
1311  GA_Offset pt0 = geometry->vertexPoint(vertices(0));
1312  GA_Offset lastpt = geometry->vertexPoint(vertices.last());
1313  // If first and last point are the same, treat it as if it's closed
1314  // in most cases, but we want to record that it's an unrolled curve,
1315  // so that the row and col cases treat it as open with a shared point.
1316  unrolled = (pt0 == lastpt);
1317  closed = unrolled;
1318  }
1319  return true;
1320 }
1321 
1322 template<typename T>
1323 static void
1324 createRotationMatrix(
1326  const UT_Vector3T<T> &cur_angles,
1327  const UT_Axis3::axis order[3])
1328 {
1329  transform.identity();
1330 
1331  for (int i = 0; i < 3; ++i)
1332  {
1333  float angle = cur_angles[i];
1334  if (angle == 0)
1335  continue;
1336  transform.rotate(order[i], angle);
1337  }
1338 }
1339 
1340 template<typename T>
1341 void
1343  const GEO_Detail *const geo,
1344  const GA_PrimitiveGroup *curve_group,
1345  const GA_RWHandleT<UT_Matrix4T<T>> &transform_attrib,
1346  const CurveFrameParms<T> &parms)
1347 {
1348  UT_AutoInterrupt interrupt("Computing curve transforms");
1349  if (interrupt.wasInterrupted())
1350  return;
1351 
1352  UT_ASSERT(transform_attrib.isValid());
1353  if (!transform_attrib.isValid())
1354  return;
1355 
1356  UT_ASSERT(transform_attrib->getOwner() == GA_ATTRIB_VERTEX);
1357 
1358  // We're going to be writing to tangent_attrib in parallel, and although
1359  // we're using GA_SplittableRange, we're splitting over primitive pages,
1360  // not vertex pages, so we have to harden in advance.
1361  transform_attrib->hardenAllPages();
1362 
1363  bool rotate_using_backbone = parms.myTangentType != TangentType::NONE;
1364  UT_Vector3T<T> target_up_vector = parms.myTargetUpVector;
1365  UT_Vector3T<T> end_target_up_vector = parms.myEndTargetUpVector;
1366  if (rotate_using_backbone)
1367  {
1368  // Normalize up_vector
1369  T length2 = target_up_vector.length2();
1370  if (length2 < theExtremelySmallLength2)
1371  {
1372  // Fall back to 0,1,0 if need be.
1373  target_up_vector.assign(0,1,0);
1374  }
1375  else
1376  {
1377  target_up_vector /= SYSsqrt(length2);
1378  }
1379  if (parms.myUseEndTargetUpVector)
1380  {
1381  // Normalize end up_vector
1382  length2 = end_target_up_vector.length2();
1383  if (length2 < theExtremelySmallLength2)
1384  {
1385  // Fall back to 0,1,0 if need be.
1386  end_target_up_vector.assign(0,1,0);
1387  }
1388  else
1389  {
1390  end_target_up_vector /= SYSsqrt(length2);
1391  }
1392  }
1393  }
1394 
1395  // If we're transforming by instance transform attributes,
1396  // e.g. N, up, pscale, scale, rot, orient, pivot, trans, xform,
1397  // set up the cache of those attributes.
1398  GA_AttributeInstanceMatrix instance_attribs;
1399  GA_ROHandleV3D instance_N;
1400  GA_ROHandleV3D instance_up;
1401  bool transform_by_instance_attribs = parms.myTransformByInstanceAttribs;
1402  if (transform_by_instance_attribs)
1403  {
1404  instance_attribs.initialize(geo->pointAttribs());
1405  if (parms.myNormalizeScales)
1406  instance_attribs.resetScales();
1407 
1408  // If there's only one of either N or up, we want to handle those separately.
1409  bool has_N = instance_attribs.myN.isValid();
1410  bool has_up = instance_attribs.myUp.isValid();
1411  if (has_N && !has_up)
1412  {
1413  instance_N = instance_attribs.myN;
1414  instance_attribs.myN.clear();
1415  }
1416  else if (has_up && !has_N)
1417  {
1418  instance_up = instance_attribs.myUp;
1419  instance_attribs.myUp.clear();
1420  }
1421 
1422  if (!instance_attribs.hasAnyAttribs())
1423  transform_by_instance_attribs = false;
1424 
1425  // We don't want to rotate both based on the backbone curves
1426  // and based on attributes, but we may still want to use pscale
1427  // with the backbone curve rotations.
1428  rotate_using_backbone &= !instance_attribs.hasNonScales();
1429  }
1430 
1431  const bool use_rotation_attrib[3] =
1432  {
1433  parms.myIncAnglePer[0] == RotationPer::ATTRIB && parms.myRotAttribs[0].isValid() && parms.myIncAngles[0] != 0.0,
1434  parms.myIncAnglePer[1] == RotationPer::ATTRIB && parms.myRotAttribs[1].isValid() && parms.myIncAngles[1] != 0.0,
1435  parms.myIncAnglePer[2] == RotationPer::ATTRIB && parms.myRotAttribs[2].isValid() && parms.myIncAngles[2] != 0.0
1436  };
1437  const GA_AttributeOwner rot_attrib_owner[3] =
1438  {
1439  use_rotation_attrib[0] ? parms.myRotAttribs[0]->getOwner() : GA_ATTRIB_INVALID,
1440  use_rotation_attrib[1] ? parms.myRotAttribs[1]->getOwner() : GA_ATTRIB_INVALID,
1441  use_rotation_attrib[2] ? parms.myRotAttribs[2]->getOwner() : GA_ATTRIB_INVALID
1442  };
1443  const bool varying_attrib[3] =
1444  {
1445  rot_attrib_owner[0] == GA_ATTRIB_POINT || rot_attrib_owner[0] == GA_ATTRIB_VERTEX,
1446  rot_attrib_owner[1] == GA_ATTRIB_POINT || rot_attrib_owner[1] == GA_ATTRIB_VERTEX,
1447  rot_attrib_owner[2] == GA_ATTRIB_POINT || rot_attrib_owner[2] == GA_ATTRIB_VERTEX
1448  };
1449 
1450  UT_Axis3::axis order[3];
1451  switch (parms.myRotationOrder)
1452  {
1453  case UT_XformOrder::XYZ: order[0] = UT_Axis3::XAXIS; order[1] = UT_Axis3::YAXIS; order[2] = UT_Axis3::ZAXIS; break;
1454  case UT_XformOrder::XZY: order[0] = UT_Axis3::XAXIS; order[1] = UT_Axis3::ZAXIS; order[2] = UT_Axis3::YAXIS; break;
1455  case UT_XformOrder::YXZ: order[0] = UT_Axis3::YAXIS; order[1] = UT_Axis3::XAXIS; order[2] = UT_Axis3::ZAXIS; break;
1456  case UT_XformOrder::YZX: order[0] = UT_Axis3::YAXIS; order[1] = UT_Axis3::ZAXIS; order[2] = UT_Axis3::XAXIS; break;
1457  case UT_XformOrder::ZXY: order[0] = UT_Axis3::ZAXIS; order[1] = UT_Axis3::XAXIS; order[2] = UT_Axis3::YAXIS; break;
1458  case UT_XformOrder::ZYX: order[0] = UT_Axis3::ZAXIS; order[1] = UT_Axis3::YAXIS; order[2] = UT_Axis3::XAXIS; break;
1459  }
1460 
1461  // Loop over primitives in parallel, using a splittable range and a lambda functor.
1463  [geo,&parms,transform_by_instance_attribs,&instance_N,&instance_up,
1464  &instance_attribs,rotate_using_backbone,target_up_vector,end_target_up_vector,
1465  &transform_attrib,&interrupt,use_rotation_attrib,rot_attrib_owner,varying_attrib,
1466  &order](const GA_SplittableRange &r)
1467  {
1468  // An array to help with computing backbone rotations
1469  UT_SmallArray<UT_Vector3T<T>, 8*sizeof(UT_Vector3T<T>)> directions;
1470  UT_SmallArray<UT_Vector3T<T>, 8*sizeof(UT_Vector3T<T>)> tangents;
1471  UT_SmallArray<UT_Vector3T<T>, 8*sizeof(UT_Vector3T<T>)> up_vectors;
1472  UT_SmallArray<UT_Vector3T<T>, 8*sizeof(UT_Vector3T<T>)> stretch_dirs;
1473  UT_SmallArray<T, 8*sizeof(T)> stretch_scales;
1474 
1475  // Shortest length of edgedira+edgedirb before max_stretch_scale must be applied
1476  const T max_stretch_length_threshold = 2.0f / parms.myMaxStretchScale;
1477 
1478  const bool rotate_using_parameters = !parms.myAngles.isZero() || !parms.myIncAngles.isZero();
1479 
1480  const int curve_u_component(parms.myRotAttribComponent);
1481 
1482  // Inside the functor, we have a sub-range of primitives, so loop over that range.
1483  // We use blockAdvance, instead of !it.atEnd() and ++it, for less looping overhead.
1485  for (GA_Iterator it(r); it.blockAdvance(start,end); )
1486  {
1487  // We probably don't need to check for interruption on every curve,
1488  // (unless people put in a curve with millions of vertices), so we
1489  // can check it once for every contiguous block of up to GA_PAGE_SIZE
1490  // primitive offsets.
1491  if (interrupt.wasInterrupted())
1492  return;
1493 
1494  // Loop over all primitives in this contiguous block of offsets.
1495  for (GA_Offset primoff = start; primoff < end; ++primoff)
1496  {
1497  const GA_OffsetListRef vertices = geo->getPrimitiveVertexList(primoff);
1498  GA_Size nedges;
1499  bool closed;
1500  bool unrolled;
1501  bool nonempty = getPolyProperties(geo, vertices, nedges, closed, unrolled);
1502  const GA_Size npoints = nedges + !closed;
1503  const GA_Size nverts = vertices.size();
1504 
1505  // Nothing to do if no vertices
1506  if (!nonempty)
1507  continue;
1508 
1509  UT_Vector3T<T> local_target_up_vector = target_up_vector;
1510  UT_Vector3T<T> local_end_target_up_vector = end_target_up_vector;
1511  if (rotate_using_backbone && parms.myTargetUpVectorAttrib.isValid())
1512  {
1514  local_target_up_vector = parms.myTargetUpVectorAttrib.get(primoff);
1515 
1516  // Normalize up_vector
1517  T length2 = local_target_up_vector.length2();
1518  if (length2 < theExtremelySmallLength2)
1519  {
1520  // Fall back to 0,1,0 if need be.
1521  local_target_up_vector.assign(0,1,0);
1522  }
1523  else
1524  {
1525  local_target_up_vector /= SYSsqrt(length2);
1526  }
1527  }
1528  if (rotate_using_backbone && parms.myUseEndTargetUpVector && parms.myEndTargetUpVectorAttrib.isValid())
1529  {
1531  local_end_target_up_vector = parms.myEndTargetUpVectorAttrib.get(primoff);
1532 
1533  // Normalize end up_vector
1534  T length2 = local_end_target_up_vector.length2();
1535  if (length2 < theExtremelySmallLength2)
1536  {
1537  // Fall back to 0,1,0 if need be.
1538  local_end_target_up_vector.assign(0,1,0);
1539  }
1540  else
1541  {
1542  local_end_target_up_vector /= SYSsqrt(length2);
1543  }
1544  }
1545 
1546  bool local_rotate_using_backbone = rotate_using_backbone;
1547  bool local_stretch_using_backbone = parms.myStretchAroundTurns;
1548 
1549  // If the curve is closed, we may have to introduce a partial twist
1550  // around the loop, instead of having the twist all at one point.
1551  // This is used to adjust how much additional twist to apply.
1552  T total_twist_around_loop = 0;
1553 
1554  const bool continuous_closed = closed && parms.myContinuousClosedCurves;
1555 
1556  // FIXME: Add parameter for use_nurbs_tangents once implemented!!!
1557  bool use_nurbs_tangents = false;
1558 
1559  RotationPer local_dangle_per[3] =
1560  {
1561  parms.myIncAnglePer[0],
1562  parms.myIncAnglePer[1],
1563  parms.myIncAnglePer[2]
1564  };
1565  UT_Vector3T<T> local_dangles = parms.myIncAngles;
1566 
1567  bool needs_length = false;
1568  bool needs_total_length = false;
1569  if (rotate_using_parameters)
1570  {
1571  needs_length =
1572  (local_dangle_per[0] == RotationPer::FULLDISTANCE || local_dangle_per[0] == RotationPer::DISTANCE) ||
1573  (local_dangle_per[1] == RotationPer::FULLDISTANCE || local_dangle_per[1] == RotationPer::DISTANCE) ||
1574  (local_dangle_per[2] == RotationPer::FULLDISTANCE || local_dangle_per[2] == RotationPer::DISTANCE);
1575  needs_total_length =
1576  (local_dangle_per[0] == RotationPer::FULLDISTANCE || (continuous_closed && local_dangle_per[0] == RotationPer::DISTANCE)) ||
1577  (local_dangle_per[1] == RotationPer::FULLDISTANCE || (continuous_closed && local_dangle_per[1] == RotationPer::DISTANCE)) ||
1578  (local_dangle_per[2] == RotationPer::FULLDISTANCE || (continuous_closed && local_dangle_per[2] == RotationPer::DISTANCE));
1579  }
1580  if (parms.myScaleRamp)
1581  {
1582  needs_length = true;
1583  needs_total_length = true;
1584  }
1585 
1586  computeSingleBackboneFrames(
1587  local_rotate_using_backbone,
1588  local_stretch_using_backbone,
1589  needs_length,
1590  directions, tangents, up_vectors,
1591  stretch_dirs, stretch_scales,
1592  total_twist_around_loop,
1593  geo, nedges, npoints, nverts, closed,
1595  instance_N,
1596  primoff, vertices, use_nurbs_tangents,
1597  parms.myMaxStretchScale, max_stretch_length_threshold,
1598  instance_up,
1599  local_target_up_vector, parms.myUseCurveNormalAsTargetUp,
1600  parms.myTargetUpVectorAtStart, continuous_closed,
1601  parms.myUseEndTargetUpVector, local_end_target_up_vector,
1602  parms.myIncAnglePer[2],
1603  parms.myRotAttribs[2], curve_u_component);
1604 
1605  UT_Vector3T<T> local_angles(0,0,0);
1606  if (use_rotation_attrib[0] && !varying_attrib[0])
1607  {
1608  GA_Offset offset = (rot_attrib_owner[0] == GA_ATTRIB_PRIMITIVE) ? primoff : GA_DETAIL_OFFSET;
1609  local_angles[0] = local_dangles[0]*parms.myRotAttribs[0].get(offset);
1610  }
1611  if (use_rotation_attrib[1] && !varying_attrib[1])
1612  {
1613  GA_Offset offset = (rot_attrib_owner[1] == GA_ATTRIB_PRIMITIVE) ? primoff : GA_DETAIL_OFFSET;
1614  local_angles[1] = local_dangles[1]*parms.myRotAttribs[1].get(offset);
1615  }
1616  if (use_rotation_attrib[2] && !varying_attrib[2])
1617  {
1618  GA_Offset offset = (rot_attrib_owner[2] == GA_ATTRIB_PRIMITIVE) ? primoff : GA_DETAIL_OFFSET;
1619  local_angles[2] = local_dangles[2]*parms.myRotAttribs[2].get(offset);
1620  }
1621 
1622  // Double-precision accumulators to avoid catastrophic roundoff error
1623  // for curves with > 16,777,216 vertices, (i.e. 2^24).
1624  double total_length = 0;
1625  double cur_length = 0;
1626  if (needs_total_length)
1627  {
1628  for (GA_Size i = 0; i < nedges; ++i)
1629  {
1630  total_length += directions(i).length();
1631  }
1632  }
1633 
1634  // For closed curves, the total twist must be a multiple of a full turn,
1635  // so we have to round it, but carefully.
1636  const bool rotate_using_dangles = !local_dangles.isZero();
1637  if (rotate_using_parameters && continuous_closed && rotate_using_dangles)
1638  {
1639  for (int axis = 0; axis < 3; ++axis)
1640  {
1641  if (local_dangle_per[axis] == RotationPer::EDGE)
1642  {
1643  local_dangles[axis] *= nedges;
1644  local_dangle_per[axis] = RotationPer::FULLEDGES;
1645  }
1646  else if (local_dangle_per[axis] == RotationPer::DISTANCE)
1647  {
1648  local_dangles[axis] *= total_length;
1649  local_dangle_per[axis] = RotationPer::FULLDISTANCE;
1650  }
1651  }
1652 
1653  // We have a separate twist due to the base frame needing to line up
1654  // with itself after going around the loop, so we want to take into
1655  // account that it'll be applied anyway.
1656  local_dangles.z() -= total_twist_around_loop;
1657  total_twist_around_loop = 0;
1658 
1659  // Round angles to multiples of a full turn
1660  local_dangles /= (2*M_PI);
1661  local_dangles.x() = SYSrint(local_dangles.x());
1662  local_dangles.y() = SYSrint(local_dangles.y());
1663  local_dangles.z() = SYSrint(local_dangles.z());
1664  local_dangles *= (2*M_PI);
1665  }
1666 
1667  // NOTE: Although we're really iterating over vertices,
1668  // we don't need to write to the last vertex of an
1669  // unrolled curve (open but last point same as first),
1670  // and the code above figured out backbone rotations
1671  // and stretching for npoints locations.
1672  for (GA_Size i = 0; i < npoints; ++i)
1673  {
1674  GA_Offset vtxoff = vertices(i);
1675  GA_Offset ptoff = geo->vertexPoint(vtxoff);
1676 
1677  UT_Vector3T<T> cur_angles = parms.myAngles+local_angles;
1678 
1679  // Apply the parameter rotations before any backbone or
1680  // attribute rotations, because they're rotations relative
1681  // to those transformed bases, as if the cross section was
1682  // rotated before applying this node.
1683  if (rotate_using_dangles)
1684  {
1685  for (int axis = 0; axis < 3; ++axis)
1686  {
1687  if (axis != 2 && local_dangles[axis] == 0)
1688  continue;
1689 
1690  const RotationPer cur_dangle_per = local_dangle_per[axis];
1691  if (cur_dangle_per == RotationPer::EDGE)
1692  {
1693  cur_angles[axis] += local_dangles[axis]*i;
1694 
1695  // Remove portion of twist added to make closed curve cycle smoothly.
1696  if (axis == 2 && total_twist_around_loop != 0)
1697  {
1698  T t = T(i)/nedges;
1699  cur_angles.z() -= total_twist_around_loop*t;
1700  }
1701  }
1702  else if (cur_dangle_per == RotationPer::FULLEDGES)
1703  {
1704  T t = T(i)/nedges;
1705  cur_angles[axis] += local_dangles[axis]*t;
1706 
1707  // Remove portion of twist added to make closed curve cycle smoothly.
1708  if (axis == 2)
1709  cur_angles.z() -= total_twist_around_loop*t;
1710  }
1711  else if (cur_dangle_per == RotationPer::DISTANCE)
1712  {
1713  cur_angles[axis] += local_dangles[axis]*cur_length;
1714 
1715  if (axis == 2 && total_length != 0 && total_twist_around_loop != 0)
1716  {
1717  // Remove portion of twist added to make closed curve cycle smoothly.
1718  T t = (cur_length/total_length);
1719  cur_angles.z() -= total_twist_around_loop*t;
1720  }
1721  }
1722  else if (cur_dangle_per == RotationPer::FULLDISTANCE)
1723  {
1724  // NOTE: total_length could be zero, but cur_length can't be more
1725  // than total_length, so checking for zero exactly suffices.
1726  if (total_length != 0)
1727  {
1728  T t = (cur_length/total_length);
1729  cur_angles[axis] += t*local_dangles[axis];
1730 
1731  // Remove portion of twist added to make closed curve cycle smoothly.
1732  if (axis == 2)
1733  cur_angles.z() -= total_twist_around_loop*t;
1734  }
1735  }
1736  else if (use_rotation_attrib[axis] && varying_attrib[axis])
1737  {
1738  UT_ASSERT_P(cur_dangle_per == RotationPer::ATTRIB);
1739  if (rot_attrib_owner[axis] == GA_ATTRIB_VERTEX)
1740  {
1741 #if 0
1742  T ustart = 0;
1743  T uend = 1;
1744  if (!closed || unrolled) {
1745  // We can really only use start and end uvs reliably if the
1746  // curve is open or unrolled, (separate start and end vertices).
1747  // Otherwise, we default to assuming 0 to 1.
1748  ustart = parms.myRotAttribs[axis].get(vertices(0), curve_u_component);
1749  uend = parms.myRotAttribs[axis].get(vertices.last(), curve_u_component);
1750  }
1751  T uspan = (uend-ustart);
1752 #endif
1753 
1754  T t = parms.myRotAttribs[axis].get(vtxoff, curve_u_component);
1755 
1756  //if (uspan != 0) {
1757  // T t = (u-ustart)/uspan;
1758  cur_angles[axis] += t*local_dangles[axis];
1759 
1760  // Remove portion of twist added to make closed curve cycle smoothly.
1761  if (axis == 2)
1762  cur_angles.z() -= total_twist_around_loop*t;
1763  //}
1764  }
1765  else
1766  {
1767  UT_ASSERT(rot_attrib_owner[axis] == GA_ATTRIB_POINT);
1768  // Assume 0 to 1, since it won't end correctly and we have to guess.
1769 
1770  T u = parms.myRotAttribs[axis].get(ptoff, curve_u_component);
1771  cur_angles[axis] += u*local_dangles[axis];
1772 
1773  // Remove portion of twist added to make closed curve cycle smoothly.
1774  if (axis == 2)
1775  cur_angles.z() -= total_twist_around_loop*u;
1776  }
1777  }
1778  }
1779  }
1780 
1782  if (rotate_using_parameters)
1783  {
1784  createRotationMatrix(transform, cur_angles, order);
1785  }
1786  else
1787  {
1788  transform.identity();
1789  }
1790 
1791  if (transform_by_instance_attribs)
1792  {
1793  UT_Matrix4T<T> instance_transform;
1794  // Save post-translate P for until after backbone rotation
1795  // or scale might be applied.
1796  instance_attribs.getMatrix(instance_transform, UT_Vector3T<T>(0,0,0), ptoff);
1797 
1798  transform *= instance_transform;
1799  }
1800 
1801  if (local_rotate_using_backbone)
1802  {
1803  UT_Vector3T<T> zaxis = tangents(i);
1804  UT_Vector3T<T> yaxis = up_vectors(i);
1805  UT_Vector3T<T> xaxis = cross(yaxis, zaxis);
1806  UT_Matrix3T<T> rotation(xaxis, yaxis, zaxis);
1807  transform *= rotation;
1808  }
1809 
1810  // Apply the stretch vector (to avoid collapsing on turns)
1811  // after any rotations, since it should always be applied
1812  // in the same final direction.
1813  if (local_stretch_using_backbone)
1814  {
1815  UT_Matrix3T<T> stretch_matrix;
1816  stretch_matrix.identity();
1817  stretch_matrix.outerproductUpdateT(stretch_scales(i)-1, stretch_dirs(i), stretch_dirs(i));
1818 
1819  transform *= stretch_matrix;
1820  }
1821 
1822  fpreal scale = parms.myUniformScale;
1823  if (parms.myScaleRamp)
1824  {
1825  const T t = (total_length != 0) ? (cur_length/total_length) : T(i)/nedges;
1826  float vals[4];
1827  parms.myScaleRamp->rampLookup(t, vals);
1828  scale *= vals[0];
1829  }
1830 
1831  // Apply the overall scale
1832  if (scale != 1)
1833  {
1834  // NOTE: We don't want to just do "transform *= scale",
1835  // because we don't want to scale the w components
1836  // of any rows, else translate below will be skewed.
1837  transform.scale(scale);
1838  }
1839 
1840  transform.translate(geo->getPos3(ptoff));
1841 
1842  transform_attrib.set(vtxoff, transform);
1843 
1844  if (needs_length && i+1 < npoints)
1845  cur_length += directions(i).length();
1846  }
1847  }
1848  }
1849  });
1850 }
1851 
1852 #define TEMPLATE_INST2(T) \
1853 template void computeCurveTransforms<T>( \
1854  const GEO_Detail *const geo, \
1855  const GA_PrimitiveGroup *curve_group, \
1856  const GA_RWHandleT<UT_Matrix4T<T>> &transform_attrib, \
1857  const CurveFrameParms<T> &parms); \
1858 /* end of macro */
1859 
1862 
1863 }
1864 
1865 } // End of HDK_Sample namespace
constexpr SYS_FORCE_INLINE T length2() const noexcept
Definition: UT_Vector3.h:356
#define SYSmax(a, b)
Definition: SYS_Math.h:1570
void UTparallelFor(const Range &range, const Body &body, const int subscribe_ratio=2, const int min_grain_size=1, const bool force_use_task_scope=true)
SYS_FORCE_INLINE GA_Primitive * getPrimitive(GA_Offset prim_off)
Definition: GA_Detail.h:429
SYS_FORCE_INLINE const GA_AttributeDict & pointAttribs() const
Definition: GEO_Detail.h:1939
Apply angle increment to each edge on top of previous edge's rotation.
Definition: GU_CurveFrame.h:69
void resetScales()
Resets only the scale attributes.
virtual fpreal getGreville(int idx, bool clamp=true, bool wrap=false) const =0
Iteration over a range of elements.
Definition: GA_Iterator.h:29
SIM_API const UT_StringHolder angle
SYS_FORCE_INLINE int getPrimitiveTypeId(GA_Offset primoff) const
Definition: GA_Primitive.h:912
GA_ROHandleT< UT_Vector3T< T > > myTargetUpVectorAttrib
#define M_PI
Definition: fmath.h:90
bool blockAdvance(GA_Offset &start, GA_Offset &end)
GLuint start
Definition: glcorearb.h:475
void setSizeNoInit(exint newsize)
Definition: UT_Array.h:695
SYS_FORCE_INLINE bool getExtraFlag() const
Synonym for isClosed()
constexpr SYS_FORCE_INLINE T & z() noexcept
Definition: UT_Vector3.h:667
int64 exint
Definition: SYS_Types.h:125
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
GLdouble s
Definition: glad.h:3009
SYS_FORCE_INLINE TO_T UTverify_cast(FROM_T from)
Definition: UT_Assert.h:229
#define UT_ASSERT_MSG_P(ZZ,...)
Definition: UT_Assert.h:158
SYS_FORCE_INLINE UT_Vector3 getPos3(GA_Offset ptoff) const
The ptoff passed is the point offset.
Definition: GA_Detail.h:185
void arbitraryPerp(const UT_Vector3T< T > &v)
Finds an arbitrary perpendicular to v, and sets this to it.
UT_Matrix2T< T > SYSlerp(const UT_Matrix2T< T > &v1, const UT_Matrix2T< T > &v2, S t)
Definition: UT_Matrix2.h:675
3D Vector class.
SYS_FORCE_INLINE bool isClosed() const
Definition: GEO_Face.h:248
float fpreal32
Definition: SYS_Types.h:200
exint GA_Size
Defines the bit width for index and offset types in GA.
Definition: GA_Types.h:236
int myRotAttribComponent
Component of myRotAttribs being read. (default 0)
GA_ROHandleT< UT_Vector3T< T > > myEndTargetUpVectorAttrib
bool hasAnyAttribs() const
Returns true if there are any attributes bound.
double fpreal64
Definition: SYS_Types.h:201
#define UT_ASSERT_MSG(ZZ,...)
Definition: UT_Assert.h:159
GA_Size GA_Offset
Definition: GA_Types.h:646
void rampLookup(fpreal pos, float values[4], int order=0) const
GA_API const UT_StringHolder scale
GLfloat f
Definition: glcorearb.h:1926
GLintptr offset
Definition: glcorearb.h:665
void identity()
Set the matrix to identity.
Definition: UT_Matrix3.h:1128
void clear()
Definition: GA_Handle.h:180
#define TEMPLATE_INST2(T)
#define UT_ASSERT_P(ZZ)
Definition: UT_Assert.h:155
void rotate(UT_Vector3T< S > &axis, T theta, int norm=1)
static UT_Matrix3T< T > dihedral(UT_Vector3T< S > &a, UT_Vector3T< S > &b, UT_Vector3T< S > &c, int norm=1)
SYS_FORCE_INLINE GA_OffsetListRef getPrimitiveVertexList(GA_Offset primoff) const
Definition: GA_Primitive.h:891
void getMatrix(UT_Matrix4 &xform, const UT_Vector3 &P, GA_Offset offset, float default_pscale=1) const
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:140
void outerproductUpdateT(T b, const UT_Vector3T< S > &v1, const UT_Vector3T< S > &v2)
Definition: UT_Matrix3.h:1431
GLuint GLuint end
Definition: glcorearb.h:475
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
Bezier or NURBS basis classes which maintain knot vectors.
Definition: GA_Basis.h:49
GLdouble GLdouble GLint GLint order
Definition: glad.h:2676
SIM_API const UT_StringHolder rotation
UT_Vector3T< fpreal64 > UT_Vector3D
SYS_FORCE_INLINE T get(GA_Offset off, int comp=0) const
Definition: GA_Handle.h:203
SYS_FORCE_INLINE GA_Offset vertexPoint(GA_Offset vertex) const
Given a vertex, return the point it references.
Definition: GA_Detail.h:529
void identity()
Set the matrix to identity.
Definition: UT_Matrix4.h:1128
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
GA_API const UT_StringHolder transform
bool evaluate(fpreal u, GEO_Vertex result, GEO_AttributeHandleList &gah, int du=0, int uoffset=-1) const
fpreal32 SYSrint(fpreal32 val)
Definition: SYS_Floor.h:163
GLdouble t
Definition: glad.h:2397
SYS_FORCE_INLINE bool isValid() const
Definition: GA_Handle.h:187
const GA_Basis * getBasis() const
Definition: GEO_Curve.h:310
bool SYSequalZero(const UT_Vector3T< T > &v)
Definition: UT_Vector3.h:1069
ToType last() const
Return the value of the last element.
void scale(T sx, T sy, T sz, T sw=1)
Definition: UT_Matrix4.h:701
void computeCurveTransforms(const GEO_Detail *const geo, const GA_PrimitiveGroup *curve_group, const GA_RWHandleT< UT_Matrix4T< T >> &transform_attrib, const CurveFrameParms< T > &parms)
GLint j
Definition: glad.h:2733
void assign(T xx=0.0f, T yy=0.0f, T zz=0.0f)
Set the values of the vector components.
Definition: UT_Vector3.h:694
GA_AttributeOwner
Definition: GA_Types.h:35
bool myTransformByInstanceAttribs
Use incoming N, up, rot, orient, pscale, scale, pivot, trans, transform.
void translate(T dx, T dy, T dz=0)
Definition: UT_Matrix4.h:773
fpreal64 fpreal
Definition: SYS_Types.h:277
constexpr SYS_FORCE_INLINE bool isZero() const noexcept
Definition: UT_Vector3.h:393
bool getPolyProperties(const GEO_Detail *geometry, const GA_OffsetListRef &vertices, exint &nedges, bool &closed, bool &unrolled)
Compute an instance transform given a set of attributes.
void constant(const T &v)
Quickly set the array to a single value.
SYS_FORCE_INLINE GA_AttributeOwner getOwner() const
Definition: GA_Attribute.h:210
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:156
SYS_FORCE_INLINE UT_StorageMathFloat_t< T > normalize() noexcept
Definition: UT_Vector3.h:376
GLboolean r
Definition: glcorearb.h:1222
#define GA_DETAIL_OFFSET
Definition: GA_Types.h:691
constexpr SYS_FORCE_INLINE T & y() noexcept
Definition: UT_Vector3.h:665
GA_Range getPrimitiveRange(const GA_PrimitiveGroup *group=0) const
Get a range of all primitives in the detail.
Definition: GA_Detail.h:1741
#define SYSmin(a, b)
Definition: SYS_Math.h:1571
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
void initialize(const GA_AttributeDict &dict, const UT_StringRef &N_name=GA_Names::N, const UT_StringRef &v_name=GA_Names::v)
constexpr T normalize(UT_FixedVector< T, D > &a) noexcept
RotationPer myIncAnglePer[3]
NOTE: myIncAnglePer[2] will also be used for ensuring closed curve continuity.
SYS_FORCE_INLINE FromType size() const
Returns the number of used elements in the list (always <= capacity())
constexpr SYS_FORCE_INLINE T & x() noexcept
Definition: UT_Vector3.h:663