HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_SolidAngle.C
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2022
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  * Functions for computing solid angles.
27  */
28 
29 #include "UT_SolidAngle.h"
30 #include "UT_BVHImpl.h"
31 
32 #include <UT/UT_Assert.h>
33 #include <UT/UT_BoundingBox.h>
34 #include <UT/UT_SmallArray.h>
35 #include <UT/UT_Vector3.h>
36 #include <VM/VM_SIMD.h>
37 #include <SYS/SYS_TypeTraits.h>
38 #include <utility>
39 
40 #define SOLID_ANGLE_TIME_PRECOMPUTE 0
41 
42 #if SOLID_ANGLE_TIME_PRECOMPUTE
43 #include <UT/UT_StopWatch.h>
44 #endif
45 
46 #define SOLID_ANGLE_DEBUG 0
47 #if SOLID_ANGLE_DEBUG
48 #include <UT/UT_Debug.h>
49 #endif
50 
51 #define TAYLOR_SERIES_ORDER 2
52 
53 namespace HDK_Sample {
54 
55 template<typename T,typename S>
56 struct UT_SolidAngle<T,S>::BoxData
57 {
58  void clear()
59  {
60  // Set everything to zero
61  memset(this,0,sizeof(*this));
62  }
63 
64  using Type = typename SYS_SelectType<UT_FixedVector<T,BVH_N>, v4uf, BVH_N==4 && SYS_IsSame<T,float>::value>::type;
65  using SType = typename SYS_SelectType<UT_FixedVector<S,BVH_N>, v4uf, BVH_N==4 && SYS_IsSame<S,float>::value>::type;
66 
67  /// An upper bound on the squared distance from myAverageP to the farthest point in the box.
68  SType myMaxPDist2;
69 
70  /// Centre of mass of the mesh surface in this box
71  UT_FixedVector<Type,3> myAverageP;
72 
73  /// Unnormalized, area-weighted normal of the mesh in this box
74  UT_FixedVector<Type,3> myN;
75 
76 #if TAYLOR_SERIES_ORDER >= 1
77  /// Values for Omega_1
78  /// @{
79  UT_FixedVector<Type,3> myNijDiag; // Nxx, Nyy, Nzz
80  Type myNxy_Nyx; // Nxy+Nyx
81  Type myNyz_Nzy; // Nyz+Nzy
82  Type myNzx_Nxz; // Nzx+Nxz
83  /// @}
84 #endif
85 
86 #if TAYLOR_SERIES_ORDER >= 2
87  /// Values for Omega_2
88  /// @{
89  UT_FixedVector<Type,3> myNijkDiag; // Nxxx, Nyyy, Nzzz
90  Type mySumPermuteNxyz; // (Nxyz+Nxzy+Nyzx+Nyxz+Nzxy+Nzyx) = 2*(Nxyz+Nyzx+Nzxy)
91  Type my2Nxxy_Nyxx; // Nxxy+Nxyx+Nyxx = 2Nxxy+Nyxx
92  Type my2Nxxz_Nzxx; // Nxxz+Nxzx+Nzxx = 2Nxxz+Nzxx
93  Type my2Nyyz_Nzyy; // Nyyz+Nyzy+Nzyy = 2Nyyz+Nzyy
94  Type my2Nyyx_Nxyy; // Nyyx+Nyxy+Nxyy = 2Nyyx+Nxyy
95  Type my2Nzzx_Nxzz; // Nzzx+Nzxz+Nxzz = 2Nzzx+Nxzz
96  Type my2Nzzy_Nyzz; // Nzzy+Nzyz+Nyzz = 2Nzzy+Nyzz
97  /// @}
98 #endif
99 };
100 
101 template<typename T,typename S>
102 UT_SolidAngle<T,S>::UT_SolidAngle()
103  : myTree()
104  , myNBoxes(0)
105  , myOrder(2)
106  , myData(nullptr)
107  , myNTriangles(0)
108  , myTrianglePoints(nullptr)
109  , myNPoints(0)
110  , myPositions(nullptr)
111 {}
112 
113 template<typename T,typename S>
114 UT_SolidAngle<T,S>::~UT_SolidAngle()
115 {
116  // Default destruction works, but this needs to be outlined
117  // to avoid having to include UT_BVHImpl.h in the header,
118  // (for the UT_UniquePtr destructor.)
119 }
120 
121 template<typename T,typename S>
122 void UT_SolidAngle<T,S>::init(
123  const int ntriangles,
124  const int *const triangle_points,
125  const int npoints,
126  const UT_Vector3T<S> *const positions,
127  const int order)
128 {
129 #if SOLID_ANGLE_DEBUG
130  UTdebugFormat("");
131  UTdebugFormat("");
132  UTdebugFormat("Building BVH for {} ntriangles on {} points:", ntriangles, npoints);
133 #endif
134  myOrder = order;
135  myNTriangles = ntriangles;
136  myTrianglePoints = triangle_points;
137  myNPoints = npoints;
138  myPositions = positions;
139 
140 #if SOLID_ANGLE_TIME_PRECOMPUTE
141  UT_StopWatch timer;
142  timer.start();
143 #endif
144  UT_SmallArray<UT::Box<S,3>> triangle_boxes;
145  triangle_boxes.setSizeNoInit(ntriangles);
146  if (ntriangles < 16*1024)
147  {
148  const int *cur_triangle_points = triangle_points;
149  for (int i = 0; i < ntriangles; ++i, cur_triangle_points += 3)
150  {
151  UT::Box<S,3> &box = triangle_boxes[i];
152  box.initBounds(positions[cur_triangle_points[0]]);
153  box.enlargeBounds(positions[cur_triangle_points[1]]);
154  box.enlargeBounds(positions[cur_triangle_points[2]]);
155  }
156  }
157  else
158  {
159  UTparallelFor(UT_BlockedRange<int>(0,ntriangles), [triangle_points,&triangle_boxes,positions](const UT_BlockedRange<int> &r)
160  {
161  const int *cur_triangle_points = triangle_points + exint(r.begin())*3;
162  for (int i = r.begin(), end = r.end(); i < end; ++i, cur_triangle_points += 3)
163  {
164  UT::Box<S,3> &box = triangle_boxes[i];
165  box.initBounds(positions[cur_triangle_points[0]]);
166  box.enlargeBounds(positions[cur_triangle_points[1]]);
167  box.enlargeBounds(positions[cur_triangle_points[2]]);
168  }
169  });
170  }
171 #if SOLID_ANGLE_TIME_PRECOMPUTE
172  double time = timer.stop();
173  UTdebugFormat("{} s to create bounding boxes.", time);
174  timer.start();
175 #endif
176  myTree.template init<UT::BVH_Heuristic::BOX_AREA,S,3>(triangle_boxes.array(), ntriangles);
177 #if SOLID_ANGLE_TIME_PRECOMPUTE
178  time = timer.stop();
179  UTdebugFormat("{} s to initialize UT_BVH structure. {} nodes", time, myTree.getNumNodes());
180 #endif
181 
182  //myTree.debugDump();
183 
184  const int nnodes = myTree.getNumNodes();
185 
186  myNBoxes = nnodes;
187  BoxData *box_data = new BoxData[nnodes];
188  myData.reset(box_data);
189 
190  // Some data are only needed during initialization.
191  struct LocalData
192  {
193  // Bounding box
194  UT_BoundingBoxT<S> myBox;
195 
196  // P and N are needed from each child for computing Nij.
197  UT_Vector3T<T> myAverageP;
198  UT_Vector3T<T> myAreaP;
199  UT_Vector3T<T> myN;
200 
201  // Unsigned area is needed for computing the average position.
202  T myArea;
203 
204 #if TAYLOR_SERIES_ORDER >= 1
205  // These are needed for computing Nijk.
206  UT_Vector3T<T> myNijDiag;
207  T myNxy; T myNyx;
208  T myNyz; T myNzy;
209  T myNzx; T myNxz;
210 #endif
211 
212 #if TAYLOR_SERIES_ORDER >= 2
213  UT_Vector3T<T> myNijkDiag; // Nxxx, Nyyy, Nzzz
214  T mySumPermuteNxyz; // (Nxyz+Nxzy+Nyzx+Nyxz+Nzxy+Nzyx) = 2*(Nxyz+Nyzx+Nzxy)
215  T my2Nxxy_Nyxx; // Nxxy+Nxyx+Nyxx = 2Nxxy+Nyxx
216  T my2Nxxz_Nzxx; // Nxxz+Nxzx+Nzxx = 2Nxxz+Nzxx
217  T my2Nyyz_Nzyy; // Nyyz+Nyzy+Nzyy = 2Nyyz+Nzyy
218  T my2Nyyx_Nxyy; // Nyyx+Nyxy+Nxyy = 2Nyyx+Nxyy
219  T my2Nzzx_Nxzz; // Nzzx+Nzxz+Nxzz = 2Nzzx+Nxzz
220  T my2Nzzy_Nyzz; // Nzzy+Nzyz+Nyzz = 2Nzzy+Nyzz
221 #endif
222  };
223 
224  struct PrecomputeFunctors
225  {
226  BoxData *const myBoxData;
227  const UT::Box<S,3> *const myTriangleBoxes;
228  const int *const myTrianglePoints;
229  const UT_Vector3T<S> *const myPositions;
230  const int myOrder;
231 
232  PrecomputeFunctors(
233  BoxData *box_data,
234  const UT::Box<S,3> *triangle_boxes,
235  const int *triangle_points,
236  const UT_Vector3T<S> *positions,
237  const int order)
238  : myBoxData(box_data)
239  , myTriangleBoxes(triangle_boxes)
240  , myTrianglePoints(triangle_points)
241  , myPositions(positions)
242  , myOrder(order)
243  {}
244  constexpr SYS_FORCE_INLINE bool pre(const int nodei, LocalData *data_for_parent) const
245  {
246  return true;
247  }
248  void item(const int itemi, const int parent_nodei, LocalData &data_for_parent) const
249  {
250  const UT_Vector3T<S> *const positions = myPositions;
251  const int *const cur_triangle_points = myTrianglePoints + 3*itemi;
252  const UT_Vector3T<T> a = positions[cur_triangle_points[0]];
253  const UT_Vector3T<T> b = positions[cur_triangle_points[1]];
254  const UT_Vector3T<T> c = positions[cur_triangle_points[2]];
255  const UT_Vector3T<T> ab = b-a;
256  const UT_Vector3T<T> ac = c-a;
257 
258  const UT::Box<S,3> &triangle_box = myTriangleBoxes[itemi];
259  data_for_parent.myBox.initBounds(triangle_box.getMin(), triangle_box.getMax());
260 
261  // Area-weighted normal (unnormalized)
262  const UT_Vector3T<T> N = T(0.5)*cross(ab,ac);
263  const T area2 = N.length2();
264  const T area = SYSsqrt(area2);
265  const UT_Vector3T<T> P = (a+b+c)/3;
266  data_for_parent.myAverageP = P;
267  data_for_parent.myAreaP = P*area;
268  data_for_parent.myN = N;
269 #if SOLID_ANGLE_DEBUG
270  UTdebugFormat("");
271  UTdebugFormat("Triangle {}: P = {}; N = {}; area = {}", itemi, P, N, area);
272  UTdebugFormat(" box = {}", data_for_parent.myBox);
273 #endif
274 
275  data_for_parent.myArea = area;
276 #if TAYLOR_SERIES_ORDER >= 1
277  const int order = myOrder;
278  if (order < 1)
279  return;
280 
281  // NOTE: Due to P being at the centroid, triangles have Nij = 0
282  // contributions to Nij.
283  data_for_parent.myNijDiag = T(0);
284  data_for_parent.myNxy = 0; data_for_parent.myNyx = 0;
285  data_for_parent.myNyz = 0; data_for_parent.myNzy = 0;
286  data_for_parent.myNzx = 0; data_for_parent.myNxz = 0;
287 #endif
288 
289 #if TAYLOR_SERIES_ORDER >= 2
290  if (order < 2)
291  return;
292 
293  // If it's zero-length, the results are zero, so we can skip.
294  if (area == 0)
295  {
296  data_for_parent.myNijkDiag = T(0);
297  data_for_parent.mySumPermuteNxyz = 0;
298  data_for_parent.my2Nxxy_Nyxx = 0;
299  data_for_parent.my2Nxxz_Nzxx = 0;
300  data_for_parent.my2Nyyz_Nzyy = 0;
301  data_for_parent.my2Nyyx_Nxyy = 0;
302  data_for_parent.my2Nzzx_Nxzz = 0;
303  data_for_parent.my2Nzzy_Nyzz = 0;
304  return;
305  }
306 
307  // We need to use the NORMALIZED normal to multiply the integrals by.
308  UT_Vector3T<T> n = N/area;
309 
310  // Figure out the order of a, b, and c in x, y, and z
311  // for use in computing the integrals for Nijk.
312  UT_Vector3T<T> values[3] = {a, b, c};
313 
314  int order_x[3] = {0,1,2};
315  if (a.x() > b.x())
316  std::swap(order_x[0],order_x[1]);
317  if (values[order_x[0]].x() > c.x())
318  std::swap(order_x[0],order_x[2]);
319  if (values[order_x[1]].x() > values[order_x[2]].x())
320  std::swap(order_x[1],order_x[2]);
321  T dx = values[order_x[2]].x() - values[order_x[0]].x();
322 
323  int order_y[3] = {0,1,2};
324  if (a.y() > b.y())
325  std::swap(order_y[0],order_y[1]);
326  if (values[order_y[0]].y() > c.y())
327  std::swap(order_y[0],order_y[2]);
328  if (values[order_y[1]].y() > values[order_y[2]].y())
329  std::swap(order_y[1],order_y[2]);
330  T dy = values[order_y[2]].y() - values[order_y[0]].y();
331 
332  int order_z[3] = {0,1,2};
333  if (a.z() > b.z())
334  std::swap(order_z[0],order_z[1]);
335  if (values[order_z[0]].z() > c.z())
336  std::swap(order_z[0],order_z[2]);
337  if (values[order_z[1]].z() > values[order_z[2]].z())
338  std::swap(order_z[1],order_z[2]);
339  T dz = values[order_z[2]].z() - values[order_z[0]].z();
340 
341  auto &&compute_integrals = [](
342  const UT_Vector3T<T> &a,
343  const UT_Vector3T<T> &b,
344  const UT_Vector3T<T> &c,
345  const UT_Vector3T<T> &P,
346  T *integral_ii,
347  T *integral_ij,
348  T *integral_ik,
349  const int i)
350  {
351 #if SOLID_ANGLE_DEBUG
352  UTdebugFormat(" Splitting on {}; a = {}; b = {}; c = {}", char('x'+i), a, b, c);
353 #endif
354  // NOTE: a, b, and c must be in order of the i axis.
355  // We're splitting the triangle at the middle i coordinate.
356  const UT_Vector3T<T> oab = b - a;
357  const UT_Vector3T<T> oac = c - a;
358  const UT_Vector3T<T> ocb = b - c;
359  UT_ASSERT_MSG_P(oac[i] > 0, "This should have been checked by the caller.");
360  const T t = oab[i]/oac[i];
361  UT_ASSERT_MSG_P(t >= 0 && t <= 1, "Either sorting must have gone wrong, or there are input NaNs.");
362 
363  const int j = (i==2) ? 0 : (i+1);
364  const int k = (j==2) ? 0 : (j+1);
365  const T jdiff = t*oac[j] - oab[j];
366  const T kdiff = t*oac[k] - oab[k];
367  const UT_Vector3T<T> cross_a(jdiff*oab[k] - kdiff*oab[j], kdiff*oab[i], jdiff*oab[i]);
368  const UT_Vector3T<T> cross_c(jdiff*ocb[k] - kdiff*ocb[j], kdiff*ocb[i], jdiff*ocb[i]);
369  const T area_scale_a = cross_a.length();
370  const T area_scale_c = cross_c.length();
371  const T Pai = a[i] - P[i];
372  const T Pci = c[i] - P[i];
373 
374  // Integral over the area of the triangle of (pi^2)dA,
375  // by splitting the triangle into two at b, the a side
376  // and the c side.
377  const T int_ii_a = area_scale_a*(T(0.5)*Pai*Pai + T(2.0/3.0)*Pai*oab[i] + T(0.25)*oab[i]*oab[i]);
378  const T int_ii_c = area_scale_c*(T(0.5)*Pci*Pci + T(2.0/3.0)*Pci*ocb[i] + T(0.25)*ocb[i]*ocb[i]);
379  *integral_ii = int_ii_a + int_ii_c;
380 #if SOLID_ANGLE_DEBUG
381  UTdebugFormat(" integral_{}{}_a = {}; integral_{}{}_c = {}", char('x'+i), char('x'+i), int_ii_a, char('x'+i), char('x'+i), int_ii_c);
382 #endif
383 
384  int jk = j;
385  T *integral = integral_ij;
386  T diff = jdiff;
387  while (true) // This only does 2 iterations, one for j and one for k
388  {
389  if (integral)
390  {
391  T obmidj = b[jk] + T(0.5)*diff;
392  T oabmidj = obmidj - a[jk];
393  T ocbmidj = obmidj - c[jk];
394  T Paj = a[jk] - P[jk];
395  T Pcj = c[jk] - P[jk];
396  // Integral over the area of the triangle of (pi*pj)dA
397  const T int_ij_a = area_scale_a*(T(0.5)*Pai*Paj + T(1.0/3.0)*Pai*oabmidj + T(1.0/3.0)*Paj*oab[i] + T(0.25)*oab[i]*oabmidj);
398  const T int_ij_c = area_scale_c*(T(0.5)*Pci*Pcj + T(1.0/3.0)*Pci*ocbmidj + T(1.0/3.0)*Pcj*ocb[i] + T(0.25)*ocb[i]*ocbmidj);
399  *integral = int_ij_a + int_ij_c;
400 #if SOLID_ANGLE_DEBUG
401  UTdebugFormat(" integral_{}{}_a = {}; integral_{}{}_c = {}", char('x'+i), char('x'+jk), int_ij_a, char('x'+i), char('x'+jk), int_ij_c);
402 #endif
403  }
404  if (jk == k)
405  break;
406  jk = k;
407  integral = integral_ik;
408  diff = kdiff;
409  }
410  };
411 
412  T integral_xx = 0;
413  T integral_xy = 0;
414  T integral_yy = 0;
415  T integral_yz = 0;
416  T integral_zz = 0;
417  T integral_zx = 0;
418  // Note that if the span of any axis is zero, the integral must be zero,
419  // since there's a factor of (p_i-P_i), i.e. value minus average,
420  // and every value must be equal to the average, giving zero.
421  if (dx > 0)
422  {
423  compute_integrals(
424  values[order_x[0]], values[order_x[1]], values[order_x[2]], P,
425  &integral_xx, ((dx >= dy && dy > 0) ? &integral_xy : nullptr), ((dx >= dz && dz > 0) ? &integral_zx : nullptr), 0);
426  }
427  if (dy > 0)
428  {
429  compute_integrals(
430  values[order_y[0]], values[order_y[1]], values[order_y[2]], P,
431  &integral_yy, ((dy >= dz && dz > 0) ? &integral_yz : nullptr), ((dx < dy && dx > 0) ? &integral_xy : nullptr), 1);
432  }
433  if (dz > 0)
434  {
435  compute_integrals(
436  values[order_z[0]], values[order_z[1]], values[order_z[2]], P,
437  &integral_zz, ((dx < dz && dx > 0) ? &integral_zx : nullptr), ((dy < dz && dy > 0) ? &integral_yz : nullptr), 2);
438  }
439 
440  UT_Vector3T<T> Niii = n*UT_Vector3T<T>(integral_xx, integral_yy, integral_zz);
441  data_for_parent.myNijkDiag = Niii;
442  data_for_parent.mySumPermuteNxyz = 2*dot(n,UT_Vector3T<T>(integral_yz, integral_zx, integral_xy));
443  T Nxxy = n.x()*integral_xy;
444  T Nxxz = n.x()*integral_zx;
445  T Nyyz = n.y()*integral_yz;
446  T Nyyx = n.y()*integral_xy;
447  T Nzzx = n.z()*integral_zx;
448  T Nzzy = n.z()*integral_yz;
449  data_for_parent.my2Nxxy_Nyxx = 2*Nxxy + n.y()*integral_xx;
450  data_for_parent.my2Nxxz_Nzxx = 2*Nxxz + n.z()*integral_xx;
451  data_for_parent.my2Nyyz_Nzyy = 2*Nyyz + n.z()*integral_yy;
452  data_for_parent.my2Nyyx_Nxyy = 2*Nyyx + n.x()*integral_yy;
453  data_for_parent.my2Nzzx_Nxzz = 2*Nzzx + n.x()*integral_zz;
454  data_for_parent.my2Nzzy_Nyzz = 2*Nzzy + n.y()*integral_zz;
455 #if SOLID_ANGLE_DEBUG
456  UTdebugFormat(" integral_xx = {}; yy = {}; zz = {}", integral_xx, integral_yy, integral_zz);
457  UTdebugFormat(" integral_xy = {}; yz = {}; zx = {}", integral_xy, integral_yz, integral_zx);
458 #endif
459 #endif
460  }
461 
462  void post(const int nodei, const int parent_nodei, LocalData *data_for_parent, const int nchildren, const LocalData *child_data_array) const
463  {
464  // NOTE: Although in the general case, data_for_parent may be null for the root call,
465  // this functor assumes that it's non-null, so the call below must pass a non-null pointer.
466 
467  BoxData &current_box_data = myBoxData[nodei];
468 
469  UT_Vector3T<T> N = child_data_array[0].myN;
470  ((T*)&current_box_data.myN[0])[0] = N[0];
471  ((T*)&current_box_data.myN[1])[0] = N[1];
472  ((T*)&current_box_data.myN[2])[0] = N[2];
473  UT_Vector3T<T> areaP = child_data_array[0].myAreaP;
474  T area = child_data_array[0].myArea;
475  UT_Vector3T<T> local_P = child_data_array[0].myAverageP;
476  ((T*)&current_box_data.myAverageP[0])[0] = local_P[0];
477  ((T*)&current_box_data.myAverageP[1])[0] = local_P[1];
478  ((T*)&current_box_data.myAverageP[2])[0] = local_P[2];
479  for (int i = 1; i < nchildren; ++i)
480  {
481  const UT_Vector3T<T> local_N = child_data_array[i].myN;
482  N += local_N;
483  ((T*)&current_box_data.myN[0])[i] = local_N[0];
484  ((T*)&current_box_data.myN[1])[i] = local_N[1];
485  ((T*)&current_box_data.myN[2])[i] = local_N[2];
486  areaP += child_data_array[i].myAreaP;
487  area += child_data_array[i].myArea;
488  const UT_Vector3T<T> local_P = child_data_array[i].myAverageP;
489  ((T*)&current_box_data.myAverageP[0])[i] = local_P[0];
490  ((T*)&current_box_data.myAverageP[1])[i] = local_P[1];
491  ((T*)&current_box_data.myAverageP[2])[i] = local_P[2];
492  }
493  for (int i = nchildren; i < BVH_N; ++i)
494  {
495  // Set to zero, just to avoid false positives for uses of uninitialized memory.
496  ((T*)&current_box_data.myN[0])[i] = 0;
497  ((T*)&current_box_data.myN[1])[i] = 0;
498  ((T*)&current_box_data.myN[2])[i] = 0;
499  ((T*)&current_box_data.myAverageP[0])[i] = 0;
500  ((T*)&current_box_data.myAverageP[1])[i] = 0;
501  ((T*)&current_box_data.myAverageP[2])[i] = 0;
502  }
503  data_for_parent->myN = N;
504  data_for_parent->myAreaP = areaP;
505  data_for_parent->myArea = area;
506 
507  UT_BoundingBoxT<S> box(child_data_array[0].myBox);
508  for (int i = 1; i < nchildren; ++i)
509  box.enlargeBounds(child_data_array[i].myBox);
510 
511  // Normalize P
512  UT_Vector3T<T> averageP;
513  if (area > 0)
514  averageP = areaP/area;
515  else
516  averageP = T(0.5)*(box.minvec() + box.maxvec());
517  data_for_parent->myAverageP = averageP;
518 
519  data_for_parent->myBox = box;
520 
521  for (int i = 0; i < nchildren; ++i)
522  {
523  const UT_BoundingBoxT<S> &local_box(child_data_array[i].myBox);
524  const UT_Vector3T<T> &local_P = child_data_array[i].myAverageP;
525  const UT_Vector3T<T> maxPDiff = SYSmax(local_P-UT_Vector3T<T>(local_box.minvec()), UT_Vector3T<T>(local_box.maxvec())-local_P);
526  ((T*)&current_box_data.myMaxPDist2)[i] = maxPDiff.length2();
527  }
528  for (int i = nchildren; i < BVH_N; ++i)
529  {
530  // This child is non-existent. If we set myMaxPDist2 to infinity, it will never
531  // use the approximation, and the traverseVector function can check for EMPTY.
532  ((T*)&current_box_data.myMaxPDist2)[i] = std::numeric_limits<T>::infinity();
533  }
534 
535 #if TAYLOR_SERIES_ORDER >= 1
536  const int order = myOrder;
537  if (order >= 1)
538  {
539  // We now have the current box's P, so we can adjust Nij and Nijk
540  data_for_parent->myNijDiag = child_data_array[0].myNijDiag;
541  data_for_parent->myNxy = 0;
542  data_for_parent->myNyx = 0;
543  data_for_parent->myNyz = 0;
544  data_for_parent->myNzy = 0;
545  data_for_parent->myNzx = 0;
546  data_for_parent->myNxz = 0;
547 #if TAYLOR_SERIES_ORDER >= 2
548  data_for_parent->myNijkDiag = child_data_array[0].myNijkDiag;
549  data_for_parent->mySumPermuteNxyz = child_data_array[0].mySumPermuteNxyz;
550  data_for_parent->my2Nxxy_Nyxx = child_data_array[0].my2Nxxy_Nyxx;
551  data_for_parent->my2Nxxz_Nzxx = child_data_array[0].my2Nxxz_Nzxx;
552  data_for_parent->my2Nyyz_Nzyy = child_data_array[0].my2Nyyz_Nzyy;
553  data_for_parent->my2Nyyx_Nxyy = child_data_array[0].my2Nyyx_Nxyy;
554  data_for_parent->my2Nzzx_Nxzz = child_data_array[0].my2Nzzx_Nxzz;
555  data_for_parent->my2Nzzy_Nyzz = child_data_array[0].my2Nzzy_Nyzz;
556 #endif
557 
558  for (int i = 1; i < nchildren; ++i)
559  {
560  data_for_parent->myNijDiag += child_data_array[i].myNijDiag;
561 #if TAYLOR_SERIES_ORDER >= 2
562  data_for_parent->myNijkDiag += child_data_array[i].myNijkDiag;
563  data_for_parent->mySumPermuteNxyz += child_data_array[i].mySumPermuteNxyz;
564  data_for_parent->my2Nxxy_Nyxx += child_data_array[i].my2Nxxy_Nyxx;
565  data_for_parent->my2Nxxz_Nzxx += child_data_array[i].my2Nxxz_Nzxx;
566  data_for_parent->my2Nyyz_Nzyy += child_data_array[i].my2Nyyz_Nzyy;
567  data_for_parent->my2Nyyx_Nxyy += child_data_array[i].my2Nyyx_Nxyy;
568  data_for_parent->my2Nzzx_Nxzz += child_data_array[i].my2Nzzx_Nxzz;
569  data_for_parent->my2Nzzy_Nyzz += child_data_array[i].my2Nzzy_Nyzz;
570 #endif
571  }
572  for (int j = 0; j < 3; ++j)
573  ((T*)&current_box_data.myNijDiag[j])[0] = child_data_array[0].myNijDiag[j];
574  ((T*)&current_box_data.myNxy_Nyx)[0] = child_data_array[0].myNxy + child_data_array[0].myNyx;
575  ((T*)&current_box_data.myNyz_Nzy)[0] = child_data_array[0].myNyz + child_data_array[0].myNzy;
576  ((T*)&current_box_data.myNzx_Nxz)[0] = child_data_array[0].myNzx + child_data_array[0].myNxz;
577  for (int j = 0; j < 3; ++j)
578  ((T*)&current_box_data.myNijkDiag[j])[0] = child_data_array[0].myNijkDiag[j];
579  ((T*)&current_box_data.mySumPermuteNxyz)[0] = child_data_array[0].mySumPermuteNxyz;
580  ((T*)&current_box_data.my2Nxxy_Nyxx)[0] = child_data_array[0].my2Nxxy_Nyxx;
581  ((T*)&current_box_data.my2Nxxz_Nzxx)[0] = child_data_array[0].my2Nxxz_Nzxx;
582  ((T*)&current_box_data.my2Nyyz_Nzyy)[0] = child_data_array[0].my2Nyyz_Nzyy;
583  ((T*)&current_box_data.my2Nyyx_Nxyy)[0] = child_data_array[0].my2Nyyx_Nxyy;
584  ((T*)&current_box_data.my2Nzzx_Nxzz)[0] = child_data_array[0].my2Nzzx_Nxzz;
585  ((T*)&current_box_data.my2Nzzy_Nyzz)[0] = child_data_array[0].my2Nzzy_Nyzz;
586  for (int i = 1; i < nchildren; ++i)
587  {
588  for (int j = 0; j < 3; ++j)
589  ((T*)&current_box_data.myNijDiag[j])[i] = child_data_array[i].myNijDiag[j];
590  ((T*)&current_box_data.myNxy_Nyx)[i] = child_data_array[i].myNxy + child_data_array[i].myNyx;
591  ((T*)&current_box_data.myNyz_Nzy)[i] = child_data_array[i].myNyz + child_data_array[i].myNzy;
592  ((T*)&current_box_data.myNzx_Nxz)[i] = child_data_array[i].myNzx + child_data_array[i].myNxz;
593  for (int j = 0; j < 3; ++j)
594  ((T*)&current_box_data.myNijkDiag[j])[i] = child_data_array[i].myNijkDiag[j];
595  ((T*)&current_box_data.mySumPermuteNxyz)[i] = child_data_array[i].mySumPermuteNxyz;
596  ((T*)&current_box_data.my2Nxxy_Nyxx)[i] = child_data_array[i].my2Nxxy_Nyxx;
597  ((T*)&current_box_data.my2Nxxz_Nzxx)[i] = child_data_array[i].my2Nxxz_Nzxx;
598  ((T*)&current_box_data.my2Nyyz_Nzyy)[i] = child_data_array[i].my2Nyyz_Nzyy;
599  ((T*)&current_box_data.my2Nyyx_Nxyy)[i] = child_data_array[i].my2Nyyx_Nxyy;
600  ((T*)&current_box_data.my2Nzzx_Nxzz)[i] = child_data_array[i].my2Nzzx_Nxzz;
601  ((T*)&current_box_data.my2Nzzy_Nyzz)[i] = child_data_array[i].my2Nzzy_Nyzz;
602  }
603  for (int i = nchildren; i < BVH_N; ++i)
604  {
605  // Set to zero, just to avoid false positives for uses of uninitialized memory.
606  for (int j = 0; j < 3; ++j)
607  ((T*)&current_box_data.myNijDiag[j])[i] = 0;
608  ((T*)&current_box_data.myNxy_Nyx)[i] = 0;
609  ((T*)&current_box_data.myNyz_Nzy)[i] = 0;
610  ((T*)&current_box_data.myNzx_Nxz)[i] = 0;
611  for (int j = 0; j < 3; ++j)
612  ((T*)&current_box_data.myNijkDiag[j])[i] = 0;
613  ((T*)&current_box_data.mySumPermuteNxyz)[i] = 0;
614  ((T*)&current_box_data.my2Nxxy_Nyxx)[i] = 0;
615  ((T*)&current_box_data.my2Nxxz_Nzxx)[i] = 0;
616  ((T*)&current_box_data.my2Nyyz_Nzyy)[i] = 0;
617  ((T*)&current_box_data.my2Nyyx_Nxyy)[i] = 0;
618  ((T*)&current_box_data.my2Nzzx_Nxzz)[i] = 0;
619  ((T*)&current_box_data.my2Nzzy_Nyzz)[i] = 0;
620  }
621 
622  for (int i = 0; i < nchildren; ++i)
623  {
624  const LocalData &child_data = child_data_array[i];
625  UT_Vector3T<T> displacement = child_data.myAverageP - UT_Vector3T<T>(data_for_parent->myAverageP);
626  UT_Vector3T<T> N = child_data.myN;
627 
628  // Adjust Nij for the change in centre P
629  data_for_parent->myNijDiag += N*displacement;
630  T Nxy = child_data.myNxy + N.x()*displacement.y();
631  T Nyx = child_data.myNyx + N.y()*displacement.x();
632  T Nyz = child_data.myNyz + N.y()*displacement.z();
633  T Nzy = child_data.myNzy + N.z()*displacement.y();
634  T Nzx = child_data.myNzx + N.z()*displacement.x();
635  T Nxz = child_data.myNxz + N.x()*displacement.z();
636 
637  data_for_parent->myNxy += Nxy;
638  data_for_parent->myNyx += Nyx;
639  data_for_parent->myNyz += Nyz;
640  data_for_parent->myNzy += Nzy;
641  data_for_parent->myNzx += Nzx;
642  data_for_parent->myNxz += Nxz;
643 
644 #if TAYLOR_SERIES_ORDER >= 2
645  if (order >= 2)
646  {
647  // Adjust Nijk for the change in centre P
648  data_for_parent->myNijkDiag += T(2)*displacement*child_data.myNijDiag + displacement*displacement*child_data.myN;
649  data_for_parent->mySumPermuteNxyz += dot(displacement, UT_Vector3T<T>(Nyz+Nzy, Nzx+Nxz, Nxy+Nyx));
650  data_for_parent->my2Nxxy_Nyxx +=
651  2*(displacement.y()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxy + N.x()*displacement.x()*displacement.y())
652  + 2*child_data.myNyx*displacement.x() + N.y()*displacement.x()*displacement.x();
653  data_for_parent->my2Nxxz_Nzxx +=
654  2*(displacement.z()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxz + N.x()*displacement.x()*displacement.z())
655  + 2*child_data.myNzx*displacement.x() + N.z()*displacement.x()*displacement.x();
656  data_for_parent->my2Nyyz_Nzyy +=
657  2*(displacement.z()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyz + N.y()*displacement.y()*displacement.z())
658  + 2*child_data.myNzy*displacement.y() + N.z()*displacement.y()*displacement.y();
659  data_for_parent->my2Nyyx_Nxyy +=
660  2*(displacement.x()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyx + N.y()*displacement.y()*displacement.x())
661  + 2*child_data.myNxy*displacement.y() + N.x()*displacement.y()*displacement.y();
662  data_for_parent->my2Nzzx_Nxzz +=
663  2*(displacement.x()*child_data.myNijDiag.z() + displacement.z()*child_data.myNzx + N.z()*displacement.z()*displacement.x())
664  + 2*child_data.myNxz*displacement.z() + N.x()*displacement.z()*displacement.z();
665  data_for_parent->my2Nzzy_Nyzz +=
666  2*(displacement.y()*child_data.myNijDiag.z() + displacement.z()*child_data.myNzy + N.z()*displacement.z()*displacement.y())
667  + 2*child_data.myNyz*displacement.z() + N.y()*displacement.z()*displacement.z();
668  }
669 #endif
670  }
671  }
672 #endif
673 #if SOLID_ANGLE_DEBUG
674  UTdebugFormat("");
675  UTdebugFormat("Node {}: nchildren = {}; maxP = {}", nodei, nchildren, SYSsqrt(current_box_data.myMaxPDist2));
676  UTdebugFormat(" P = {}; N = {}", current_box_data.myAverageP, current_box_data.myN);
677 #if TAYLOR_SERIES_ORDER >= 1
678  UTdebugFormat(" Nii = {}", current_box_data.myNijDiag);
679  UTdebugFormat(" Nxy+Nyx = {}; Nyz+Nzy = {}; Nyz+Nzy = {}", current_box_data.myNxy_Nyx, current_box_data.myNyz_Nzy, current_box_data.myNzx_Nxz);
680 #if TAYLOR_SERIES_ORDER >= 2
681  UTdebugFormat(" Niii = {}; 2(Nxyz+Nyzx+Nzxy) = {}", current_box_data.myNijkDiag, current_box_data.mySumPermuteNxyz);
682  UTdebugFormat(" 2Nxxy+Nyxx = {}; 2Nxxz+Nzxx = {}", current_box_data.my2Nxxy_Nyxx, current_box_data.my2Nxxz_Nzxx);
683  UTdebugFormat(" 2Nyyz+Nzyy = {}; 2Nyyx+Nxyy = {}", current_box_data.my2Nyyz_Nzyy, current_box_data.my2Nyyx_Nxyy);
684  UTdebugFormat(" 2Nzzx+Nxzz = {}; 2Nzzy+Nyzz = {}", current_box_data.my2Nzzx_Nxzz, current_box_data.my2Nzzy_Nyzz);
685 #endif
686 #endif
687 #endif
688  }
689  };
690 
691 #if SOLID_ANGLE_TIME_PRECOMPUTE
692  timer.start();
693 #endif
694  const PrecomputeFunctors functors(box_data, triangle_boxes.array(), triangle_points, positions, order);
695  // NOTE: post-functor relies on non-null data_for_parent, so we have to pass one.
696  LocalData local_data;
697  myTree.template traverseParallel<LocalData>(4096, functors, &local_data);
698  //myTree.template traverse<LocalData>(functors);
699 #if SOLID_ANGLE_TIME_PRECOMPUTE
700  time = timer.stop();
701  UTdebugFormat("{} s to precompute coefficients.", time);
702 #endif
703 }
704 
705 template<typename T,typename S>
706 void UT_SolidAngle<T, S>::clear()
707 {
708  myTree.clear();
709  myNBoxes = 0;
710  myOrder = 2;
711  myData.reset();
712  myNTriangles = 0;
713  myTrianglePoints = nullptr;
714  myNPoints = 0;
715  myPositions = nullptr;
716 }
717 
718 template<typename T,typename S>
719 T UT_SolidAngle<T, S>::computeSolidAngle(const UT_Vector3T<T> &query_point, const T accuracy_scale) const
720 {
721  const T accuracy_scale2 = accuracy_scale*accuracy_scale;
722 
723  struct SolidAngleFunctors
724  {
725  const BoxData *const myBoxData;
726  const UT_Vector3T<T> myQueryPoint;
727  const T myAccuracyScale2;
728  const UT_Vector3T<S> *const myPositions;
729  const int *const myTrianglePoints;
730  const int myOrder;
731 
732  SolidAngleFunctors(
733  const BoxData *const box_data,
734  const UT_Vector3T<T> &query_point,
735  const T accuracy_scale2,
736  const int order,
737  const UT_Vector3T<S> *const positions,
738  const int *const triangle_points)
739  : myBoxData(box_data)
740  , myQueryPoint(query_point)
741  , myAccuracyScale2(accuracy_scale2)
742  , myOrder(order)
743  , myPositions(positions)
744  , myTrianglePoints(triangle_points)
745  {}
746  uint pre(const int nodei, T *data_for_parent) const
747  {
748  const BoxData &data = myBoxData[nodei];
749  const typename BoxData::Type maxP2 = data.myMaxPDist2;
750  UT_FixedVector<typename BoxData::Type,3> q;
751  q[0] = typename BoxData::Type(myQueryPoint.x());
752  q[1] = typename BoxData::Type(myQueryPoint.y());
753  q[2] = typename BoxData::Type(myQueryPoint.z());
754  q -= data.myAverageP;
755  const typename BoxData::Type qlength2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
756 
757  // If the query point is within a factor of accuracy_scale of the box radius,
758  // it's assumed to be not a good enough approximation, so it needs to descend.
759  // TODO: Is there a way to estimate the error?
760  SYS_STATIC_ASSERT_MSG((SYS_IsSame<typename BoxData::Type,v4uf>::value), "FIXME: Implement support for other tuple types!");
761  v4uu descend_mask = (qlength2 <= maxP2*myAccuracyScale2);
762  uint descend_bitmask = _mm_movemask_ps(V4SF(descend_mask.vector));
763  constexpr uint allchildbits = ((uint(1)<<BVH_N)-1);
764  if (descend_bitmask == allchildbits)
765  {
766  *data_for_parent = 0;
767  return allchildbits;
768  }
769 
770  // qlength2 must be non-zero, since it's strictly greater than something.
771  // We still need to be careful for NaNs, though, because the 4th power might cause problems.
772  const typename BoxData::Type qlength_m2 = typename BoxData::Type(1.0)/qlength2;
773  const typename BoxData::Type qlength_m1 = SYSsqrt(qlength_m2);
774 
775  // Normalize q to reduce issues with overflow/underflow, since we'd need the 7th power
776  // if we didn't normalize, and (1e-6)^-7 = 1e42, which overflows single-precision.
777  q *= qlength_m1;
778 
779  typename BoxData::Type Omega_approx = -qlength_m2*dot(q,data.myN);
780 #if TAYLOR_SERIES_ORDER >= 1
781  const int order = myOrder;
782  if (order >= 1)
783  {
784  const UT_FixedVector<typename BoxData::Type,3> q2 = q*q;
785  const typename BoxData::Type qlength_m3 = qlength_m2*qlength_m1;
786  const typename BoxData::Type Omega_1 =
787  qlength_m3*(data.myNijDiag[0] + data.myNijDiag[1] + data.myNijDiag[2]
788  -typename BoxData::Type(3.0)*(dot(q2,data.myNijDiag) +
789  q[0]*q[1]*data.myNxy_Nyx +
790  q[0]*q[2]*data.myNzx_Nxz +
791  q[1]*q[2]*data.myNyz_Nzy));
792  Omega_approx += Omega_1;
793 #if TAYLOR_SERIES_ORDER >= 2
794  if (order >= 2)
795  {
796  const UT_FixedVector<typename BoxData::Type,3> q3 = q2*q;
797  const typename BoxData::Type qlength_m4 = qlength_m2*qlength_m2;
798  typename BoxData::Type temp0[3] = {
799  data.my2Nyyx_Nxyy+data.my2Nzzx_Nxzz,
800  data.my2Nzzy_Nyzz+data.my2Nxxy_Nyxx,
801  data.my2Nxxz_Nzxx+data.my2Nyyz_Nzyy
802  };
803  typename BoxData::Type temp1[3] = {
804  q[1]*data.my2Nxxy_Nyxx + q[2]*data.my2Nxxz_Nzxx,
805  q[2]*data.my2Nyyz_Nzyy + q[0]*data.my2Nyyx_Nxyy,
806  q[0]*data.my2Nzzx_Nxzz + q[1]*data.my2Nzzy_Nyzz
807  };
808  const typename BoxData::Type Omega_2 =
809  qlength_m4*(typename BoxData::Type(1.5)*dot(q, typename BoxData::Type(3)*data.myNijkDiag + UT_FixedVector<typename BoxData::Type,3>(temp0))
810  -typename BoxData::Type(7.5)*(dot(q3,data.myNijkDiag) + q[0]*q[1]*q[2]*data.mySumPermuteNxyz + dot(q2, UT_FixedVector<typename BoxData::Type,3>(temp1))));
811  Omega_approx += Omega_2;
812  }
813 #endif
814  }
815 #endif
816 
817  // If q is so small that we got NaNs and we just have a
818  // small bounding box, it needs to descend.
819  const v4uu mask = Omega_approx.isFinite() & ~descend_mask;
820  Omega_approx = Omega_approx & mask;
821  descend_bitmask = (~_mm_movemask_ps(V4SF(mask.vector))) & allchildbits;
822 
823  T sum = Omega_approx[0];
824  for (int i = 1; i < BVH_N; ++i)
825  sum += Omega_approx[i];
826  *data_for_parent = sum;
827 
828  return descend_bitmask;
829  }
830  void item(const int itemi, const int parent_nodei, T &data_for_parent) const
831  {
832  const UT_Vector3T<S> *const positions = myPositions;
833  const int *const cur_triangle_points = myTrianglePoints + 3*itemi;
834  const UT_Vector3T<T> a = positions[cur_triangle_points[0]];
835  const UT_Vector3T<T> b = positions[cur_triangle_points[1]];
836  const UT_Vector3T<T> c = positions[cur_triangle_points[2]];
837 
838  data_for_parent = UTsignedSolidAngleTri(a, b, c, myQueryPoint);
839  }
840  SYS_FORCE_INLINE void post(const int nodei, const int parent_nodei, T *data_for_parent, const int nchildren, const T *child_data_array, const uint descend_bits) const
841  {
842  T sum = (descend_bits&1) ? child_data_array[0] : 0;
843  for (int i = 1; i < nchildren; ++i)
844  sum += ((descend_bits>>i)&1) ? child_data_array[i] : 0;
845 
846  *data_for_parent += sum;
847  }
848  };
849  const SolidAngleFunctors functors(myData.get(), query_point, accuracy_scale2, myOrder, myPositions, myTrianglePoints);
850 
851  T sum;
852  myTree.traverseVector(functors, &sum);
853  return sum;
854 }
855 
856 template<typename T,typename S>
857 struct UT_SubtendedAngle<T,S>::BoxData
858 {
859  void clear()
860  {
861  // Set everything to zero
862  memset(this,0,sizeof(*this));
863  }
864 
865  using Type = typename SYS_SelectType<UT_FixedVector<T,BVH_N>, v4uf, BVH_N==4 && SYS_IsSame<T,float>::value>::type;
866  using SType = typename SYS_SelectType<UT_FixedVector<S,BVH_N>, v4uf, BVH_N==4 && SYS_IsSame<S,float>::value>::type;
867 
868  /// An upper bound on the squared distance from myAverageP to the farthest point in the box.
869  SType myMaxPDist2;
870 
871  /// Centre of mass of the mesh surface in this box
872  UT_FixedVector<Type,2> myAverageP;
873 
874  /// Unnormalized, area-weighted normal of the mesh in this box
875  UT_FixedVector<Type,2> myN;
876 
877  /// Values for Omega_1
878  /// @{
879  UT_FixedVector<Type,2> myNijDiag; // Nxx, Nyy
880  Type myNxy_Nyx; // Nxy+Nyx
881  /// @}
882 
883  /// Values for Omega_2
884  /// @{
885  UT_FixedVector<Type,2> myNijkDiag; // Nxxx, Nyyy
886  Type my2Nxxy_Nyxx; // Nxxy+Nxyx+Nyxx = 2Nxxy+Nyxx
887  Type my2Nyyx_Nxyy; // Nyyx+Nyxy+Nxyy = 2Nyyx+Nxyy
888  /// @}
889 };
890 
891 template<typename T,typename S>
892 UT_SubtendedAngle<T,S>::UT_SubtendedAngle()
893  : myTree()
894  , myNBoxes(0)
895  , myOrder(2)
896  , myData(nullptr)
897  , myNSegments(0)
898  , mySegmentPoints(nullptr)
899  , myNPoints(0)
900  , myPositions(nullptr)
901 {}
902 
903 template<typename T,typename S>
904 UT_SubtendedAngle<T,S>::~UT_SubtendedAngle()
905 {
906  // Default destruction works, but this needs to be outlined
907  // to avoid having to include UT_BVHImpl.h in the header,
908  // (for the UT_UniquePtr destructor.)
909 }
910 
911 template<typename T,typename S>
912 void UT_SubtendedAngle<T,S>::init(
913  const int nsegments,
914  const int *const segment_points,
915  const int npoints,
916  const UT_Vector2T<S> *const positions,
917  const int order)
918 {
919 #if SOLID_ANGLE_DEBUG
920  UTdebugFormat("");
921  UTdebugFormat("");
922  UTdebugFormat("Building BVH for {} segments on {} points:", nsegments, npoints);
923 #endif
924  myOrder = order;
925  myNSegments = nsegments;
926  mySegmentPoints = segment_points;
927  myNPoints = npoints;
928  myPositions = positions;
929 
930 #if SOLID_ANGLE_TIME_PRECOMPUTE
931  UT_StopWatch timer;
932  timer.start();
933 #endif
934  UT_SmallArray<UT::Box<S,2>> segment_boxes;
935  segment_boxes.setSizeNoInit(nsegments);
936  if (nsegments < 16*1024)
937  {
938  const int *cur_segment_points = segment_points;
939  for (int i = 0; i < nsegments; ++i, cur_segment_points += 2)
940  {
941  UT::Box<S,2> &box = segment_boxes[i];
942  box.initBounds(positions[cur_segment_points[0]]);
943  box.enlargeBounds(positions[cur_segment_points[1]]);
944  }
945  }
946  else
947  {
948  UTparallelFor(UT_BlockedRange<int>(0,nsegments), [segment_points,&segment_boxes,positions](const UT_BlockedRange<int> &r)
949  {
950  const int *cur_segment_points = segment_points + exint(r.begin())*2;
951  for (int i = r.begin(), end = r.end(); i < end; ++i, cur_segment_points += 2)
952  {
953  UT::Box<S,2> &box = segment_boxes[i];
954  box.initBounds(positions[cur_segment_points[0]]);
955  box.enlargeBounds(positions[cur_segment_points[1]]);
956  }
957  });
958  }
959 #if SOLID_ANGLE_TIME_PRECOMPUTE
960  double time = timer.stop();
961  UTdebugFormat("{} s to create bounding boxes.", time);
962  timer.start();
963 #endif
964  myTree.template init<UT::BVH_Heuristic::BOX_AREA,S,2>(segment_boxes.array(), nsegments);
965 #if SOLID_ANGLE_TIME_PRECOMPUTE
966  time = timer.stop();
967  UTdebugFormat("{} s to initialize UT_BVH structure. {} nodes", time, myTree.getNumNodes());
968 #endif
969 
970  //myTree.debugDump();
971 
972  const int nnodes = myTree.getNumNodes();
973 
974  myNBoxes = nnodes;
975  BoxData *box_data = new BoxData[nnodes];
976  myData.reset(box_data);
977 
978  // Some data are only needed during initialization.
979  struct LocalData
980  {
981  // Bounding box
982  UT::Box<S,2> myBox;
983 
984  // P and N are needed from each child for computing Nij.
985  UT_Vector2T<T> myAverageP;
986  UT_Vector2T<T> myLengthP;
987  UT_Vector2T<T> myN;
988 
989  // Unsigned length is needed for computing the average position.
990  T myLength;
991 
992  // These are needed for computing Nijk.
993  UT_Vector2T<T> myNijDiag;
994  T myNxy; T myNyx;
995 
996  UT_Vector2T<T> myNijkDiag; // Nxxx, Nyyy
997  T my2Nxxy_Nyxx; // Nxxy+Nxyx+Nyxx = 2Nxxy+Nyxx
998  T my2Nyyx_Nxyy; // Nyyx+Nyxy+Nxyy = 2Nyyx+Nxyy
999  };
1000 
1001  struct PrecomputeFunctors
1002  {
1003  BoxData *const myBoxData;
1004  const UT::Box<S,2> *const mySegmentBoxes;
1005  const int *const mySegmentPoints;
1006  const UT_Vector2T<S> *const myPositions;
1007  const int myOrder;
1008 
1009  PrecomputeFunctors(
1010  BoxData *box_data,
1011  const UT::Box<S,2> *segment_boxes,
1012  const int *segment_points,
1013  const UT_Vector2T<S> *positions,
1014  const int order)
1015  : myBoxData(box_data)
1016  , mySegmentBoxes(segment_boxes)
1017  , mySegmentPoints(segment_points)
1018  , myPositions(positions)
1019  , myOrder(order)
1020  {}
1021  constexpr SYS_FORCE_INLINE bool pre(const int nodei, LocalData *data_for_parent) const
1022  {
1023  return true;
1024  }
1025  void item(const int itemi, const int parent_nodei, LocalData &data_for_parent) const
1026  {
1027  const UT_Vector2T<S> *const positions = myPositions;
1028  const int *const cur_segment_points = mySegmentPoints + 2*itemi;
1029  const UT_Vector2T<T> a = positions[cur_segment_points[0]];
1030  const UT_Vector2T<T> b = positions[cur_segment_points[1]];
1031  const UT_Vector2T<T> ab = b-a;
1032 
1033  const UT::Box<S,2> &segment_box = mySegmentBoxes[itemi];
1034  data_for_parent.myBox = segment_box;
1035 
1036  // Length-weighted normal (unnormalized)
1037  const UT_Vector2T<T> N(ab.y(),-ab.x());
1038  const T length2 = ab.length2();
1039  const T length = SYSsqrt(length2);
1040  const UT_Vector2T<T> P = T(0.5)*(a+b);
1041  data_for_parent.myAverageP = P;
1042  data_for_parent.myLengthP = P*length;
1043  data_for_parent.myN = N;
1044 #if SOLID_ANGLE_DEBUG
1045  UTdebugFormat("");
1046  UTdebugFormat("Triangle {}: P = {}; N = {}; length = {}", itemi, P, N, length);
1047  UTdebugFormat(" box = {}", data_for_parent.myBox);
1048 #endif
1049 
1050  data_for_parent.myLength = length;
1051  const int order = myOrder;
1052  if (order < 1)
1053  return;
1054 
1055  // NOTE: Due to P being at the centroid, segments have Nij = 0
1056  // contributions to Nij.
1057  data_for_parent.myNijDiag = T(0);
1058  data_for_parent.myNxy = 0; data_for_parent.myNyx = 0;
1059 
1060  if (order < 2)
1061  return;
1062 
1063  // If it's zero-length, the results are zero, so we can skip.
1064  if (length == 0)
1065  {
1066  data_for_parent.myNijkDiag = T(0);
1067  data_for_parent.my2Nxxy_Nyxx = 0;
1068  data_for_parent.my2Nyyx_Nxyy = 0;
1069  return;
1070  }
1071 
1072  T integral_xx = ab.x()*ab.x()/T(12);
1073  T integral_xy = ab.x()*ab.y()/T(12);
1074  T integral_yy = ab.y()*ab.y()/T(12);
1075  data_for_parent.myNijkDiag.assign(integral_xx*N.x(),integral_yy*N.y());
1076  T Nxxy = N.x()*integral_xy;
1077  T Nyxx = N.y()*integral_xx;
1078  T Nyyx = N.y()*integral_xy;
1079  T Nxyy = N.x()*integral_yy;
1080  data_for_parent.my2Nxxy_Nyxx = 2*Nxxy + Nyxx;
1081  data_for_parent.my2Nyyx_Nxyy = 2*Nyyx + Nxyy;
1082 #if SOLID_ANGLE_DEBUG
1083  UTdebugFormat(" integral_xx = {}; yy = {}", integral_xx, integral_yy);
1084  UTdebugFormat(" integral_xy = {}", integral_xy);
1085 #endif
1086  }
1087 
1088  void post(const int nodei, const int parent_nodei, LocalData *data_for_parent, const int nchildren, const LocalData *child_data_array) const
1089  {
1090  // NOTE: Although in the general case, data_for_parent may be null for the root call,
1091  // this functor assumes that it's non-null, so the call below must pass a non-null pointer.
1092 
1093  BoxData &current_box_data = myBoxData[nodei];
1094 
1095  UT_Vector2T<T> N = child_data_array[0].myN;
1096  ((T*)&current_box_data.myN[0])[0] = N[0];
1097  ((T*)&current_box_data.myN[1])[0] = N[1];
1098  UT_Vector2T<T> lengthP = child_data_array[0].myLengthP;
1099  T length = child_data_array[0].myLength;
1100  const UT_Vector2T<T> local_P = child_data_array[0].myAverageP;
1101  ((T*)&current_box_data.myAverageP[0])[0] = local_P[0];
1102  ((T*)&current_box_data.myAverageP[1])[0] = local_P[1];
1103  for (int i = 1; i < nchildren; ++i)
1104  {
1105  const UT_Vector2T<T> local_N = child_data_array[i].myN;
1106  N += local_N;
1107  ((T*)&current_box_data.myN[0])[i] = local_N[0];
1108  ((T*)&current_box_data.myN[1])[i] = local_N[1];
1109  lengthP += child_data_array[i].myLengthP;
1110  length += child_data_array[i].myLength;
1111  const UT_Vector2T<T> local_P = child_data_array[i].myAverageP;
1112  ((T*)&current_box_data.myAverageP[0])[i] = local_P[0];
1113  ((T*)&current_box_data.myAverageP[1])[i] = local_P[1];
1114  }
1115  for (int i = nchildren; i < BVH_N; ++i)
1116  {
1117  // Set to zero, just to avoid false positives for uses of uninitialized memory.
1118  ((T*)&current_box_data.myN[0])[i] = 0;
1119  ((T*)&current_box_data.myN[1])[i] = 0;
1120  ((T*)&current_box_data.myAverageP[0])[i] = 0;
1121  ((T*)&current_box_data.myAverageP[1])[i] = 0;
1122  }
1123  data_for_parent->myN = N;
1124  data_for_parent->myLengthP = lengthP;
1125  data_for_parent->myLength = length;
1126 
1127  UT::Box<S,2> box(child_data_array[0].myBox);
1128  for (int i = 1; i < nchildren; ++i)
1129  box.combine(child_data_array[i].myBox);
1130 
1131  // Normalize P
1132  UT_Vector2T<T> averageP;
1133  if (length > 0)
1134  averageP = lengthP/length;
1135  else
1136  averageP = T(0.5)*(box.getMin() + box.getMax());
1137  data_for_parent->myAverageP = averageP;
1138 
1139  data_for_parent->myBox = box;
1140 
1141  for (int i = 0; i < nchildren; ++i)
1142  {
1143  const UT::Box<S,2> &local_box(child_data_array[i].myBox);
1144  const UT_Vector2T<T> &local_P = child_data_array[i].myAverageP;
1145  const UT_Vector2T<T> maxPDiff = SYSmax(local_P-UT_Vector2T<T>(local_box.getMin()), UT_Vector2T<T>(local_box.getMax())-local_P);
1146  ((T*)&current_box_data.myMaxPDist2)[i] = maxPDiff.length2();
1147  }
1148  for (int i = nchildren; i < BVH_N; ++i)
1149  {
1150  // This child is non-existent. If we set myMaxPDist2 to infinity, it will never
1151  // use the approximation, and the traverseVector function can check for EMPTY.
1152  ((T*)&current_box_data.myMaxPDist2)[i] = std::numeric_limits<T>::infinity();
1153  }
1154 
1155  const int order = myOrder;
1156  if (order >= 1)
1157  {
1158  // We now have the current box's P, so we can adjust Nij and Nijk
1159  data_for_parent->myNijDiag = child_data_array[0].myNijDiag;
1160  data_for_parent->myNxy = 0;
1161  data_for_parent->myNyx = 0;
1162  data_for_parent->myNijkDiag = child_data_array[0].myNijkDiag;
1163  data_for_parent->my2Nxxy_Nyxx = child_data_array[0].my2Nxxy_Nyxx;
1164  data_for_parent->my2Nyyx_Nxyy = child_data_array[0].my2Nyyx_Nxyy;
1165 
1166  for (int i = 1; i < nchildren; ++i)
1167  {
1168  data_for_parent->myNijDiag += child_data_array[i].myNijDiag;
1169  data_for_parent->myNijkDiag += child_data_array[i].myNijkDiag;
1170  data_for_parent->my2Nxxy_Nyxx += child_data_array[i].my2Nxxy_Nyxx;
1171  data_for_parent->my2Nyyx_Nxyy += child_data_array[i].my2Nyyx_Nxyy;
1172  }
1173  for (int j = 0; j < 2; ++j)
1174  ((T*)&current_box_data.myNijDiag[j])[0] = child_data_array[0].myNijDiag[j];
1175  ((T*)&current_box_data.myNxy_Nyx)[0] = child_data_array[0].myNxy + child_data_array[0].myNyx;
1176  for (int j = 0; j < 2; ++j)
1177  ((T*)&current_box_data.myNijkDiag[j])[0] = child_data_array[0].myNijkDiag[j];
1178  ((T*)&current_box_data.my2Nxxy_Nyxx)[0] = child_data_array[0].my2Nxxy_Nyxx;
1179  ((T*)&current_box_data.my2Nyyx_Nxyy)[0] = child_data_array[0].my2Nyyx_Nxyy;
1180  for (int i = 1; i < nchildren; ++i)
1181  {
1182  for (int j = 0; j < 2; ++j)
1183  ((T*)&current_box_data.myNijDiag[j])[i] = child_data_array[i].myNijDiag[j];
1184  ((T*)&current_box_data.myNxy_Nyx)[i] = child_data_array[i].myNxy + child_data_array[i].myNyx;
1185  for (int j = 0; j < 2; ++j)
1186  ((T*)&current_box_data.myNijkDiag[j])[i] = child_data_array[i].myNijkDiag[j];
1187  ((T*)&current_box_data.my2Nxxy_Nyxx)[i] = child_data_array[i].my2Nxxy_Nyxx;
1188  ((T*)&current_box_data.my2Nyyx_Nxyy)[i] = child_data_array[i].my2Nyyx_Nxyy;
1189  }
1190  for (int i = nchildren; i < BVH_N; ++i)
1191  {
1192  // Set to zero, just to avoid false positives for uses of uninitialized memory.
1193  for (int j = 0; j < 2; ++j)
1194  ((T*)&current_box_data.myNijDiag[j])[i] = 0;
1195  ((T*)&current_box_data.myNxy_Nyx)[i] = 0;
1196  for (int j = 0; j < 2; ++j)
1197  ((T*)&current_box_data.myNijkDiag[j])[i] = 0;
1198  ((T*)&current_box_data.my2Nxxy_Nyxx)[i] = 0;
1199  ((T*)&current_box_data.my2Nyyx_Nxyy)[i] = 0;
1200  }
1201 
1202  for (int i = 0; i < nchildren; ++i)
1203  {
1204  const LocalData &child_data = child_data_array[i];
1205  UT_Vector2T<T> displacement = child_data.myAverageP - UT_Vector2T<T>(data_for_parent->myAverageP);
1206  UT_Vector2T<T> N = child_data.myN;
1207 
1208  // Adjust Nij for the change in centre P
1209  data_for_parent->myNijDiag += N*displacement;
1210  T Nxy = child_data.myNxy + N.x()*displacement.y();
1211  T Nyx = child_data.myNyx + N.y()*displacement.x();
1212 
1213  data_for_parent->myNxy += Nxy;
1214  data_for_parent->myNyx += Nyx;
1215 
1216  if (order >= 2)
1217  {
1218  // Adjust Nijk for the change in centre P
1219  data_for_parent->myNijkDiag += T(2)*displacement*child_data.myNijDiag + displacement*displacement*child_data.myN;
1220  data_for_parent->my2Nxxy_Nyxx +=
1221  2*(displacement.y()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxy + N.x()*displacement.x()*displacement.y())
1222  + 2*child_data.myNyx*displacement.x() + N.y()*displacement.x()*displacement.x();
1223  data_for_parent->my2Nyyx_Nxyy +=
1224  2*(displacement.x()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyx + N.y()*displacement.y()*displacement.x())
1225  + 2*child_data.myNxy*displacement.y() + N.x()*displacement.y()*displacement.y();
1226  }
1227  }
1228  }
1229 #if SOLID_ANGLE_DEBUG
1230  UTdebugFormat("");
1231  UTdebugFormat("Node {}: nchildren = {}; maxP = {}", nodei, nchildren, SYSsqrt(current_box_data.myMaxPDist2));
1232  UTdebugFormat(" P = {}; N = {}", current_box_data.myAverageP, current_box_data.myN);
1233  UTdebugFormat(" Nii = {}", current_box_data.myNijDiag);
1234  UTdebugFormat(" Nxy+Nyx = {}", current_box_data.myNxy_Nyx);
1235  UTdebugFormat(" Niii = {}", current_box_data.myNijkDiag);
1236  UTdebugFormat(" 2Nxxy+Nyxx = {}; 2Nyyx+Nxyy = {}", current_box_data.my2Nxxy_Nyxx, current_box_data.my2Nyyx_Nxyy);
1237 #endif
1238  }
1239  };
1240 
1241 #if SOLID_ANGLE_TIME_PRECOMPUTE
1242  timer.start();
1243 #endif
1244  const PrecomputeFunctors functors(box_data, segment_boxes.array(), segment_points, positions, order);
1245  // NOTE: post-functor relies on non-null data_for_parent, so we have to pass one.
1246  LocalData local_data;
1247  myTree.template traverseParallel<LocalData>(4096, functors, &local_data);
1248  //myTree.template traverse<LocalData>(functors);
1249 #if SOLID_ANGLE_TIME_PRECOMPUTE
1250  time = timer.stop();
1251  UTdebugFormat("{} s to precompute coefficients.", time);
1252 #endif
1253 }
1254 
1255 template<typename T,typename S>
1256 void UT_SubtendedAngle<T, S>::clear()
1257 {
1258  myTree.clear();
1259  myNBoxes = 0;
1260  myOrder = 2;
1261  myData.reset();
1262  myNSegments = 0;
1263  mySegmentPoints = nullptr;
1264  myNPoints = 0;
1265  myPositions = nullptr;
1266 }
1267 
1268 template<typename T,typename S>
1269 T UT_SubtendedAngle<T, S>::computeAngle(const UT_Vector2T<T> &query_point, const T accuracy_scale) const
1270 {
1271  const T accuracy_scale2 = accuracy_scale*accuracy_scale;
1272 
1273  struct AngleFunctors
1274  {
1275  const BoxData *const myBoxData;
1276  const UT_Vector2T<T> myQueryPoint;
1277  const T myAccuracyScale2;
1278  const UT_Vector2T<S> *const myPositions;
1279  const int *const mySegmentPoints;
1280  const int myOrder;
1281 
1282  AngleFunctors(
1283  const BoxData *const box_data,
1284  const UT_Vector2T<T> &query_point,
1285  const T accuracy_scale2,
1286  const int order,
1287  const UT_Vector2T<S> *const positions,
1288  const int *const segment_points)
1289  : myBoxData(box_data)
1290  , myQueryPoint(query_point)
1291  , myAccuracyScale2(accuracy_scale2)
1292  , myOrder(order)
1293  , myPositions(positions)
1294  , mySegmentPoints(segment_points)
1295  {}
1296  uint pre(const int nodei, T *data_for_parent) const
1297  {
1298  const BoxData &data = myBoxData[nodei];
1299  const typename BoxData::Type maxP2 = data.myMaxPDist2;
1300  UT_FixedVector<typename BoxData::Type,2> q;
1301  q[0] = typename BoxData::Type(myQueryPoint.x());
1302  q[1] = typename BoxData::Type(myQueryPoint.y());
1303  q -= data.myAverageP;
1304  const typename BoxData::Type qlength2 = q[0]*q[0] + q[1]*q[1];
1305 
1306  // If the query point is within a factor of accuracy_scale of the box radius,
1307  // it's assumed to be not a good enough approximation, so it needs to descend.
1308  // TODO: Is there a way to estimate the error?
1309  SYS_STATIC_ASSERT_MSG((SYS_IsSame<typename BoxData::Type,v4uf>::value), "FIXME: Implement support for other tuple types!");
1310  v4uu descend_mask = (qlength2 <= maxP2*myAccuracyScale2);
1311  uint descend_bitmask = _mm_movemask_ps(V4SF(descend_mask.vector));
1312  constexpr uint allchildbits = ((uint(1)<<BVH_N)-1);
1313  if (descend_bitmask == allchildbits)
1314  {
1315  *data_for_parent = 0;
1316  return allchildbits;
1317  }
1318 
1319  // qlength2 must be non-zero, since it's strictly greater than something.
1320  // We still need to be careful for NaNs, though, because the 4th power might cause problems.
1321  const typename BoxData::Type qlength_m2 = typename BoxData::Type(1.0)/qlength2;
1322  const typename BoxData::Type qlength_m1 = SYSsqrt(qlength_m2);
1323 
1324  // Normalize q to reduce issues with overflow/underflow, since we'd need the 6th power
1325  // if we didn't normalize, and (1e-7)^-6 = 1e42, which overflows single-precision.
1326  q *= qlength_m1;
1327 
1328  typename BoxData::Type Omega_approx = -qlength_m1*dot(q,data.myN);
1329  const int order = myOrder;
1330  if (order >= 1)
1331  {
1332  const UT_FixedVector<typename BoxData::Type,2> q2 = q*q;
1333  const typename BoxData::Type Omega_1 =
1334  qlength_m2*(data.myNijDiag[0] + data.myNijDiag[1]
1335  -typename BoxData::Type(2.0)*(dot(q2,data.myNijDiag) +
1336  q[0]*q[1]*data.myNxy_Nyx));
1337  Omega_approx += Omega_1;
1338  if (order >= 2)
1339  {
1340  const UT_FixedVector<typename BoxData::Type,2> q3 = q2*q;
1341  const typename BoxData::Type qlength_m3 = qlength_m2*qlength_m1;
1342  typename BoxData::Type temp0[2] = {
1343  data.my2Nyyx_Nxyy,
1344  data.my2Nxxy_Nyxx
1345  };
1346  typename BoxData::Type temp1[2] = {
1347  q[1]*data.my2Nxxy_Nyxx,
1348  q[0]*data.my2Nyyx_Nxyy
1349  };
1350  const typename BoxData::Type Omega_2 =
1351  qlength_m3*(dot(q, typename BoxData::Type(3)*data.myNijkDiag + UT_FixedVector<typename BoxData::Type,2>(temp0))
1352  -typename BoxData::Type(4.0)*(dot(q3,data.myNijkDiag) + dot(q2, UT_FixedVector<typename BoxData::Type,2>(temp1))));
1353  Omega_approx += Omega_2;
1354  }
1355  }
1356 
1357  // If q is so small that we got NaNs and we just have a
1358  // small bounding box, it needs to descend.
1359  const v4uu mask = Omega_approx.isFinite() & ~descend_mask;
1360  Omega_approx = Omega_approx & mask;
1361  descend_bitmask = (~_mm_movemask_ps(V4SF(mask.vector))) & allchildbits;
1362 
1363  T sum = Omega_approx[0];
1364  for (int i = 1; i < BVH_N; ++i)
1365  sum += Omega_approx[i];
1366  *data_for_parent = sum;
1367 
1368  return descend_bitmask;
1369  }
1370  void item(const int itemi, const int parent_nodei, T &data_for_parent) const
1371  {
1372  const UT_Vector2T<S> *const positions = myPositions;
1373  const int *const cur_segment_points = mySegmentPoints + 2*itemi;
1374  const UT_Vector2T<T> a = positions[cur_segment_points[0]];
1375  const UT_Vector2T<T> b = positions[cur_segment_points[1]];
1376 
1377  data_for_parent = UTsignedAngleSegment(a, b, myQueryPoint);
1378  }
1379  SYS_FORCE_INLINE void post(const int nodei, const int parent_nodei, T *data_for_parent, const int nchildren, const T *child_data_array, const uint descend_bits) const
1380  {
1381  T sum = (descend_bits&1) ? child_data_array[0] : 0;
1382  for (int i = 1; i < nchildren; ++i)
1383  sum += ((descend_bits>>i)&1) ? child_data_array[i] : 0;
1384 
1385  *data_for_parent += sum;
1386  }
1387  };
1388  const AngleFunctors functors(myData.get(), query_point, accuracy_scale2, myOrder, myPositions, mySegmentPoints);
1389 
1390  T sum;
1391  myTree.traverseVector(functors, &sum);
1392  return sum;
1393 }
1394 
1395 // Instantiate our templates.
1396 template class UT_SolidAngle<fpreal32,fpreal32>;
1397 // FIXME: The SIMD parts will need to be handled differently in order to support fpreal64.
1398 //template class UT_SolidAngle<fpreal64,fpreal32>;
1399 //template class UT_SolidAngle<fpreal64,fpreal64>;
1400 template class UT_SubtendedAngle<fpreal32,fpreal32>;
1401 //template class UT_SubtendedAngle<fpreal64,fpreal32>;
1402 //template class UT_SubtendedAngle<fpreal64,fpreal64>;
1403 
1404 }
UT_Vector3T::length2
constexpr SYS_FORCE_INLINE T length2() const noexcept
Definition: UT_Vector3.h:346
SYSmax
#define SYSmax(a, b)
Definition: SYS_Math.h:1513
HDK_Sample::UT_SolidAngle::UT_SolidAngle
UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
Definition: UT_SolidAngle.C:102
HDK_Sample::UT::Box::getMax
SYS_FORCE_INLINE UT_FixedVector< T, NAXES > getMax() const noexcept
Definition: UT_BVH.h:172
b
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
HDK_Sample::UT_SubtendedAngle::BoxData::SType
typename SYS_SelectType< UT_FixedVector< S, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< S, float >::value >::type SType
Definition: UT_SolidAngle.C:866
order
GLuint GLdouble GLdouble GLint GLint order
Definition: glew.h:3460
UT_BoundingBoxT
Axis-aligned bounding box (AABB).
Definition: GEO_Detail.h:43
SYS_STATIC_ASSERT_MSG
#define SYS_STATIC_ASSERT_MSG(expr, msg)
Definition: SYS_StaticAssert.h:26
myOrder
int myOrder
Definition: GT_CurveEval.h:263
HDK_Sample::UT::Box::initBounds
SYS_FORCE_INLINE void initBounds() noexcept
Definition: UT_BVH.h:97
std::swap
void swap(UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &a, UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &b)
Definition: UT_ArraySet.h:1629
GT_Names::time
GT_API const UT_StringHolder time
SOP_AttribCreate_2_0Enums::Type
Type
Definition: SOP_AttribCreate-2.0.proto.h:49
UT_Array::setSizeNoInit
void setSizeNoInit(exint newsize)
Definition: UT_Array.h:658
HDK_Sample::UT_SolidAngle::BoxData::myNyz_Nzy
Type myNyz_Nzy
Definition: UT_SolidAngle.C:81
c
const GLfloat * c
Definition: glew.h:16631
HDK_Sample::UT_SubtendedAngle::BoxData::myNijDiag
UT_FixedVector< Type, 2 > myNijDiag
Definition: UT_SolidAngle.C:879
HDK_Sample::UT_SolidAngle::BoxData::mySumPermuteNxyz
Type mySumPermuteNxyz
Definition: UT_SolidAngle.C:90
UT_Vector3.h
UT_Vector3T::z
constexpr SYS_FORCE_INLINE T & z() noexcept
Definition: UT_Vector3.h:656
exint
int64 exint
Definition: SYS_Types.h:125
UT_BoundingBoxT::maxvec
UT_Vector3T< T > maxvec() const
Definition: UT_BoundingBox.h:120
UT_Vector2T::y
T & y()
Definition: UT_Vector2.h:428
HDK_Sample::UT_SubtendedAngle::BoxData
Definition: UT_SolidAngle.C:857
UT_Array::array
T * array()
Definition: UT_Array.h:782
UT_SmallArray.h
HDK_Sample::UT_SolidAngle::BoxData::myNijkDiag
UT_FixedVector< Type, 3 > myNijkDiag
Definition: UT_SolidAngle.C:89
HDK_Sample::UT_SolidAngle::BoxData::Type
typename SYS_SelectType< UT_FixedVector< T, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< T, float >::value >::type Type
Definition: UT_SolidAngle.C:64
UT_Vector3T::length
constexpr SYS_FORCE_INLINE T length() const noexcept
Definition: UT_Vector3.h:351
UT_ASSERT_MSG_P
#define UT_ASSERT_MSG_P(ZZ,...)
Definition: UT_Assert.h:155
UTparallelFor
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=false)
Definition: UT_ParallelUtil.h:289
SYS_IsSame
Definition: SYS_TypeTraits.h:109
UT_Vector3T< S >
UT_Vector2T
2D Vector class.
Definition: UT_Vector2.h:157
HDK_Sample::UT_SolidAngle::~UT_SolidAngle
~UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
Definition: UT_SolidAngle.C:114
t
GLdouble GLdouble t
Definition: glew.h:1403
UT_Assert.h
UT_Vector2T::x
T & x()
Definition: UT_Vector2.h:426
VM_SIMD.h
values
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1602
x
GLint GLenum GLint x
Definition: glcorearb.h:409
q
GLdouble GLdouble GLdouble GLdouble q
Definition: glew.h:1419
HDK_Sample::UT_SolidAngle::BoxData::my2Nxxz_Nzxx
Type my2Nxxz_Nzxx
Definition: UT_SolidAngle.C:92
mask
GLint GLuint mask
Definition: glcorearb.h:124
HDK_Sample::UT_SolidAngle::BoxData::SType
typename SYS_SelectType< UT_FixedVector< S, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< S, float >::value >::type SType
Definition: UT_SolidAngle.C:65
v4uu
Definition: VM_SIMD.h:48
dot
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:130
HDK_Sample::UT_SolidAngle::clear
void clear()
Frees myTree and myData, and clears the rest.
Definition: UT_SolidAngle.C:706
HDK_Sample::UT_SolidAngle::BoxData::myN
UT_FixedVector< Type, 3 > myN
Unnormalized, area-weighted normal of the mesh in this box.
Definition: UT_SolidAngle.C:74
a
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
end
GLuint GLuint end
Definition: glcorearb.h:475
SYS_FORCE_INLINE
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
HDK_Sample::UT_SubtendedAngle::BoxData::clear
void clear()
Definition: UT_SolidAngle.C:859
HDK_Sample::UT_SubtendedAngle::init
void init(const int nsegments, const int *const segment_points, const int npoints, const UT_Vector2T< S > *const positions, const int order=2)
Definition: UT_SolidAngle.C:912
HDK_Sample::UT_SolidAngle::computeSolidAngle
T computeSolidAngle(const UT_Vector3T< T > &query_point, const T accuracy_scale=T(2.0)) const
Definition: UT_SolidAngle.C:719
HDK_Sample::UT_SolidAngle::BoxData::myMaxPDist2
SType myMaxPDist2
An upper bound on the squared distance from myAverageP to the farthest point in the box...
Definition: UT_SolidAngle.C:68
v4uf
Definition: VM_SIMD.h:188
z
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
HDK_Sample::UT_SolidAngle::BoxData::my2Nyyz_Nzyy
Type my2Nyyz_Nzyy
Definition: UT_SolidAngle.C:93
UT_Debug.h
UT_BoundingBoxT::enlargeBounds
void enlargeBounds(const UT_Vector3T< T > &min, const UT_Vector3T< T > &max)
Definition: UT_BoundingBox.h:750
HDK_Sample::UT_SolidAngle::BoxData
Definition: UT_SolidAngle.C:56
OBJ_MatchTransform::T
HDK_Sample::UT_SubtendedAngle::BoxData::myAverageP
UT_FixedVector< Type, 2 > myAverageP
Centre of mass of the mesh surface in this box.
Definition: UT_SolidAngle.C:872
HDK_Sample::UT_SolidAngle::init
void init(const int ntriangles, const int *const triangle_points, const int npoints, const UT_Vector3T< S > *const positions, const int order=2)
Definition: UT_SolidAngle.C:122
HDK_Sample::UT_SubtendedAngle::BoxData::Type
typename SYS_SelectType< UT_FixedVector< T, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< T, float >::value >::type Type
Definition: UT_SolidAngle.C:865
HDK_Sample::UT_SubtendedAngle::computeAngle
T computeAngle(const UT_Vector2T< T > &query_point, const T accuracy_scale=T(2.0)) const
Definition: UT_SolidAngle.C:1269
n
GLdouble n
Definition: glcorearb.h:2008
UT_StopWatch
Definition: UT_StopWatch.h:23
HDK_Sample::UT_SubtendedAngle
Definition: UT_SolidAngle.h:294
HDK_Sample::UT_SolidAngle::BoxData::my2Nxxy_Nyxx
Type my2Nxxy_Nyxx
Definition: UT_SolidAngle.C:91
HDK_Sample::UT::Box::enlargeBounds
SYS_FORCE_INLINE void enlargeBounds(const UT_FixedVector< T, NAXES > &pt) noexcept
Definition: UT_BVH.h:155
HDK_Sample::UTsignedAngleSegment
T UTsignedAngleSegment(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, const UT_Vector2T< T > &query)
Definition: UT_SolidAngle.h:259
HDK_Sample::UT_SubtendedAngle::BoxData::myNijkDiag
UT_FixedVector< Type, 2 > myNijkDiag
Definition: UT_SolidAngle.C:885
HDK_Sample::UT_SolidAngle::BoxData::my2Nzzy_Nyzz
Type my2Nzzy_Nyzz
Definition: UT_SolidAngle.C:96
length
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
UT_Vector2T::length2
constexpr SYS_FORCE_INLINE T length2() const noexcept
Definition: UT_Vector2.h:292
HDK_Sample::UT::Box
Definition: UT_BVH.h:53
HDK_Sample::UT_SubtendedAngle::BoxData::myNxy_Nyx
Type myNxy_Nyx
Definition: UT_SolidAngle.C:880
UT_StopWatch::start
void start()
HDK_Sample::UT_SubtendedAngle::BoxData::my2Nyyx_Nxyy
Type my2Nyyx_Nxyy
Definition: UT_SolidAngle.C:887
HDK_Sample::UT_SubtendedAngle::UT_SubtendedAngle
UT_SubtendedAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
Definition: UT_SolidAngle.C:892
box
GLsizei const GLint box[]
Definition: glew.h:11654
HDK_Sample::UT_SolidAngle::BoxData::myNzx_Nxz
Type myNzx_Nxz
Definition: UT_SolidAngle.C:82
OBJ_MatchTransform::S
HDK_Sample::UT::Box::getMin
SYS_FORCE_INLINE UT_FixedVector< T, NAXES > getMin() const noexcept
Definition: UT_BVH.h:163
HDK_Sample::UT_SolidAngle::BoxData::myAverageP
UT_FixedVector< Type, 3 > myAverageP
Centre of mass of the mesh surface in this box.
Definition: UT_SolidAngle.C:71
UT_BoundingBoxT::minvec
UT_Vector3T< T > minvec() const
Definition: UT_BoundingBox.h:118
v4uu::vector
v4si vector
Definition: VM_SIMD.h:185
HDK_Sample::UT_SolidAngle::BoxData::myNxy_Nyx
Type myNxy_Nyx
Definition: UT_SolidAngle.C:80
UT_SmallArray
Definition: UT_SmallArray.h:28
HDK_Sample::UT_SubtendedAngle::BoxData::my2Nxxy_Nyxx
Type my2Nxxy_Nyxx
Definition: UT_SolidAngle.C:886
GA_Names::N
GA_API const UT_StringHolder N
UT_FixedVector< Type, 3 >
HDK_Sample::UT_SubtendedAngle::clear
void clear()
Frees myTree and myData, and clears the rest.
Definition: UT_SolidAngle.C:1256
SYS_SelectType
Definition: SYS_TypeTraits.h:205
UT_BoundingBox.h
HDK_Sample::UT_SolidAngle::BoxData::myNijDiag
UT_FixedVector< Type, 3 > myNijDiag
Definition: UT_SolidAngle.C:79
SYS_TypeTraits.h
UT_StopWatch::stop
fpreal64 stop()
UT_BVHImpl.h
V4SF
#define V4SF(A)
Definition: VM_BasicFunc.h:68
HDK_Sample::UT_SubtendedAngle::BoxData::myMaxPDist2
SType myMaxPDist2
An upper bound on the squared distance from myAverageP to the farthest point in the box...
Definition: UT_SolidAngle.C:869
HDK_Sample::UT_SolidAngle::BoxData::clear
void clear()
Definition: UT_SolidAngle.C:58
HDK_Sample::UT_SubtendedAngle::~UT_SubtendedAngle
~UT_SubtendedAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
Definition: UT_SolidAngle.C:904
type
type
Definition: core.h:1059
UT_Vector3T::y
constexpr SYS_FORCE_INLINE T & y() noexcept
Definition: UT_Vector3.h:654
r
GLboolean r
Definition: glcorearb.h:1222
UT_SolidAngle.h
HDK_Sample::UT_SolidAngle::BoxData::my2Nzzx_Nxzz
Type my2Nzzx_Nxzz
Definition: UT_SolidAngle.C:95
HDK_Sample::UT::Box::combine
SYS_FORCE_INLINE void combine(const Box< T, NAXES > &src) noexcept
Definition: UT_BVH.h:121
UT_BlockedRange
Declare prior to use.
Definition: UT_ParallelUtil.h:56
uint
unsigned int uint
Definition: SYS_Types.h:45
HDK_Sample::UTsignedSolidAngleTri
T UTsignedSolidAngleTri(const UT_Vector3T< T > &a, const UT_Vector3T< T > &b, const UT_Vector3T< T > &c, const UT_Vector3T< T > &query)
Definition: UT_SolidAngle.h:50
cross
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
Definition: SIM_DerVector3.h:304
y
GLint y
Definition: glcorearb.h:103
HDK_Sample::UT_SolidAngle
Definition: UT_SolidAngle.h:203
HDK_Sample::UT_SubtendedAngle::BoxData::myN
UT_FixedVector< Type, 2 > myN
Unnormalized, area-weighted normal of the mesh in this box.
Definition: UT_SolidAngle.C:875
data
Definition: format.h:895
UT_StopWatch.h
GA_Names::area
GA_API const UT_StringHolder area
HDK_Sample::UT_SolidAngle::BoxData::my2Nyyx_Nxyy
Type my2Nyyx_Nxyy
Definition: UT_SolidAngle.C:94
UTdebugFormat
#define UTdebugFormat(...)
Definition: UT_Debug.h:137
UT_Vector3T::x
constexpr SYS_FORCE_INLINE T & x() noexcept
Definition: UT_Vector3.h:652