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