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