40 #define SOLID_ANGLE_TIME_PRECOMPUTE 0
42 #if SOLID_ANGLE_TIME_PRECOMPUTE
46 #define SOLID_ANGLE_DEBUG 0
51 #define TAYLOR_SERIES_ORDER 2
53 namespace HDK_Sample {
55 template<
typename T,
typename S>
61 memset(
this,0,
sizeof(*
this));
76 #if TAYLOR_SERIES_ORDER >= 1
86 #if TAYLOR_SERIES_ORDER >= 2
101 template<
typename T,
typename S>
108 , myTrianglePoints(nullptr)
110 , myPositions(nullptr)
113 template<
typename T,
typename S>
121 template<
typename T,
typename S>
123 const int ntriangles,
124 const int *
const triangle_points,
129 #if SOLID_ANGLE_DEBUG
132 UTdebugFormat(
"Building BVH for {} ntriangles on {} points:", ntriangles, npoints);
135 myNTriangles = ntriangles;
136 myTrianglePoints = triangle_points;
138 myPositions = positions;
140 #if SOLID_ANGLE_TIME_PRECOMPUTE
146 if (ntriangles < 16*1024)
148 const int *cur_triangle_points = triangle_points;
149 for (
int i = 0; i < ntriangles; ++i, cur_triangle_points += 3)
152 box.
initBounds(positions[cur_triangle_points[0]]);
161 const int *cur_triangle_points = triangle_points +
exint(r.begin())*3;
162 for (
int i = r.begin(),
end = r.end(); i <
end; ++i, cur_triangle_points += 3)
165 box.
initBounds(positions[cur_triangle_points[0]]);
171 #if SOLID_ANGLE_TIME_PRECOMPUTE
176 myTree.template init<UT::BVH_Heuristic::BOX_AREA,S,3>(triangle_boxes.
array(), ntriangles);
177 #if SOLID_ANGLE_TIME_PRECOMPUTE
179 UTdebugFormat(
"{} s to initialize UT_BVH structure. {} nodes", time, myTree.getNumNodes());
184 const int nnodes = myTree.getNumNodes();
188 myData.reset(box_data);
204 #if TAYLOR_SERIES_ORDER >= 1
212 #if TAYLOR_SERIES_ORDER >= 2
224 struct PrecomputeFunctors
228 const int *
const myTrianglePoints;
235 const int *triangle_points,
238 : myBoxData(box_data)
239 , myTriangleBoxes(triangle_boxes)
240 , myTrianglePoints(triangle_points)
241 , myPositions(positions)
244 constexpr
SYS_FORCE_INLINE bool pre(
const int nodei, LocalData *data_for_parent)
const
248 void item(
const int itemi,
const int parent_nodei, LocalData &data_for_parent)
const
251 const int *
const cur_triangle_points = myTrianglePoints + 3*itemi;
258 const UT::Box<S,3> &triangle_box = myTriangleBoxes[itemi];
259 data_for_parent.myBox.initBounds(triangle_box.
getMin(), triangle_box.
getMax());
264 const T area = SYSsqrt(area2);
266 data_for_parent.myAverageP = P;
267 data_for_parent.myAreaP = P*
area;
268 data_for_parent.myN =
N;
269 #if SOLID_ANGLE_DEBUG
271 UTdebugFormat(
"Triangle {}: P = {}; N = {}; area = {}", itemi, P, N, area);
275 data_for_parent.myArea =
area;
276 #if TAYLOR_SERIES_ORDER >= 1
283 data_for_parent.myNijDiag =
T(0);
284 data_for_parent.myNxy = 0; data_for_parent.myNyx = 0;
285 data_for_parent.myNyz = 0; data_for_parent.myNzy = 0;
286 data_for_parent.myNzx = 0; data_for_parent.myNxz = 0;
289 #if TAYLOR_SERIES_ORDER >= 2
296 data_for_parent.myNijkDiag =
T(0);
297 data_for_parent.mySumPermuteNxyz = 0;
298 data_for_parent.my2Nxxy_Nyxx = 0;
299 data_for_parent.my2Nxxz_Nzxx = 0;
300 data_for_parent.my2Nyyz_Nzyy = 0;
301 data_for_parent.my2Nyyx_Nxyy = 0;
302 data_for_parent.my2Nzzx_Nxzz = 0;
303 data_for_parent.my2Nzzy_Nyzz = 0;
314 int order_x[3] = {0,1,2};
317 if (values[order_x[0]].
x() > c.
x())
319 if (values[order_x[1]].
x() > values[order_x[2]].
x())
321 T dx = values[order_x[2]].
x() - values[order_x[0]].
x();
323 int order_y[3] = {0,1,2};
326 if (values[order_y[0]].
y() > c.
y())
328 if (values[order_y[1]].
y() > values[order_y[2]].
y())
330 T dy = values[order_y[2]].
y() - values[order_y[0]].
y();
332 int order_z[3] = {0,1,2};
335 if (values[order_z[0]].
z() > c.
z())
337 if (values[order_z[1]].
z() > values[order_z[2]].
z())
339 T dz = values[order_z[2]].
z() - values[order_z[0]].
z();
341 auto &&compute_integrals = [](
351 #if SOLID_ANGLE_DEBUG
352 UTdebugFormat(
" Splitting on {}; a = {}; b = {}; c = {}",
char(
'x'+i), a, b, c);
359 UT_ASSERT_MSG_P(oac[i] > 0,
"This should have been checked by the caller.");
360 const T t = oab[i]/oac[i];
361 UT_ASSERT_MSG_P(t >= 0 && t <= 1,
"Either sorting must have gone wrong, or there are input NaNs.");
363 const int j = (i==2) ? 0 : (i+1);
364 const int k = (j==2) ? 0 : (j+1);
365 const T jdiff = t*oac[j] - oab[j];
366 const T kdiff = t*oac[k] - oab[k];
367 const UT_Vector3T<T> cross_a(jdiff*oab[k] - kdiff*oab[j], kdiff*oab[i], jdiff*oab[i]);
368 const UT_Vector3T<T> cross_c(jdiff*ocb[k] - kdiff*ocb[j], kdiff*ocb[i], jdiff*ocb[i]);
369 const T area_scale_a = cross_a.
length();
370 const T area_scale_c = cross_c.
length();
371 const T Pai = a[i] - P[i];
372 const T Pci = c[i] - P[i];
377 const T int_ii_a = area_scale_a*(
T(0.5)*Pai*Pai +
T(2.0/3.0)*Pai*oab[i] +
T(0.25)*oab[i]*oab[i]);
378 const T int_ii_c = area_scale_c*(
T(0.5)*Pci*Pci +
T(2.0/3.0)*Pci*ocb[i] +
T(0.25)*ocb[i]*ocb[i]);
379 *integral_ii = int_ii_a + int_ii_c;
380 #if SOLID_ANGLE_DEBUG
381 UTdebugFormat(
" integral_{}{}_a = {}; integral_{}{}_c = {}",
char(
'x'+i),
char(
'x'+i), int_ii_a,
char(
'x'+i),
char(
'x'+i), int_ii_c);
385 T *integral = integral_ij;
391 T obmidj = b[jk] +
T(0.5)*diff;
392 T oabmidj = obmidj - a[jk];
393 T ocbmidj = obmidj - c[jk];
394 T Paj = a[jk] - P[jk];
395 T Pcj = c[jk] - P[jk];
397 const T int_ij_a = area_scale_a*(
T(0.5)*Pai*Paj +
T(1.0/3.0)*Pai*oabmidj +
T(1.0/3.0)*Paj*oab[i] +
T(0.25)*oab[i]*oabmidj);
398 const T int_ij_c = area_scale_c*(
T(0.5)*Pci*Pcj +
T(1.0/3.0)*Pci*ocbmidj +
T(1.0/3.0)*Pcj*ocb[i] +
T(0.25)*ocb[i]*ocbmidj);
399 *integral = int_ij_a + int_ij_c;
400 #if SOLID_ANGLE_DEBUG
401 UTdebugFormat(
" integral_{}{}_a = {}; integral_{}{}_c = {}",
char(
'x'+i),
char(
'x'+jk), int_ij_a,
char(
'x'+i),
char(
'x'+jk), int_ij_c);
407 integral = integral_ik;
424 values[order_x[0]], values[order_x[1]], values[order_x[2]], P,
425 &integral_xx, ((dx >= dy && dy > 0) ? &integral_xy :
nullptr), ((dx >= dz && dz > 0) ? &integral_zx :
nullptr), 0);
430 values[order_y[0]], values[order_y[1]], values[order_y[2]], P,
431 &integral_yy, ((dy >= dz && dz > 0) ? &integral_yz :
nullptr), ((dx < dy && dx > 0) ? &integral_xy :
nullptr), 1);
436 values[order_z[0]], values[order_z[1]], values[order_z[2]], P,
437 &integral_zz, ((dx < dz && dx > 0) ? &integral_zx :
nullptr), ((dy < dz && dy > 0) ? &integral_yz :
nullptr), 2);
441 data_for_parent.myNijkDiag = Niii;
442 data_for_parent.mySumPermuteNxyz = 2*
dot(n,
UT_Vector3T<T>(integral_yz, integral_zx, integral_xy));
443 T Nxxy = n.
x()*integral_xy;
444 T Nxxz = n.
x()*integral_zx;
445 T Nyyz = n.
y()*integral_yz;
446 T Nyyx = n.
y()*integral_xy;
447 T Nzzx = n.
z()*integral_zx;
448 T Nzzy = n.
z()*integral_yz;
449 data_for_parent.my2Nxxy_Nyxx = 2*Nxxy + n.
y()*integral_xx;
450 data_for_parent.my2Nxxz_Nzxx = 2*Nxxz + n.
z()*integral_xx;
451 data_for_parent.my2Nyyz_Nzyy = 2*Nyyz + n.
z()*integral_yy;
452 data_for_parent.my2Nyyx_Nxyy = 2*Nyyx + n.
x()*integral_yy;
453 data_for_parent.my2Nzzx_Nxzz = 2*Nzzx + n.
x()*integral_zz;
454 data_for_parent.my2Nzzy_Nyzz = 2*Nzzy + n.
y()*integral_zz;
455 #if SOLID_ANGLE_DEBUG
456 UTdebugFormat(
" integral_xx = {}; yy = {}; zz = {}", integral_xx, integral_yy, integral_zz);
457 UTdebugFormat(
" integral_xy = {}; yz = {}; zx = {}", integral_xy, integral_yz, integral_zx);
462 void post(
const int nodei,
const int parent_nodei, LocalData *data_for_parent,
const int nchildren,
const LocalData *child_data_array)
const
467 BoxData ¤t_box_data = myBoxData[nodei];
470 ((
T*)¤t_box_data.
myN[0])[0] = N[0];
471 ((
T*)¤t_box_data.
myN[1])[0] = N[1];
472 ((
T*)¤t_box_data.
myN[2])[0] = N[2];
474 T area = child_data_array[0].myArea;
476 ((
T*)¤t_box_data.
myAverageP[0])[0] = local_P[0];
477 ((
T*)¤t_box_data.
myAverageP[1])[0] = local_P[1];
478 ((
T*)¤t_box_data.
myAverageP[2])[0] = local_P[2];
479 for (
int i = 1; i < nchildren; ++i)
483 ((
T*)¤t_box_data.
myN[0])[i] = local_N[0];
484 ((
T*)¤t_box_data.
myN[1])[i] = local_N[1];
485 ((
T*)¤t_box_data.
myN[2])[i] = local_N[2];
486 areaP += child_data_array[i].myAreaP;
487 area += child_data_array[i].myArea;
489 ((
T*)¤t_box_data.
myAverageP[0])[i] = local_P[0];
490 ((
T*)¤t_box_data.
myAverageP[1])[i] = local_P[1];
491 ((
T*)¤t_box_data.
myAverageP[2])[i] = local_P[2];
493 for (
int i = nchildren; i < BVH_N; ++i)
496 ((
T*)¤t_box_data.
myN[0])[i] = 0;
497 ((
T*)¤t_box_data.
myN[1])[i] = 0;
498 ((
T*)¤t_box_data.
myN[2])[i] = 0;
503 data_for_parent->myN =
N;
504 data_for_parent->myAreaP = areaP;
505 data_for_parent->myArea =
area;
508 for (
int i = 1; i < nchildren; ++i)
514 averageP = areaP/
area;
517 data_for_parent->myAverageP = averageP;
519 data_for_parent->myBox =
box;
521 for (
int i = 0; i < nchildren; ++i)
528 for (
int i = nchildren; i < BVH_N; ++i)
532 ((
T*)¤t_box_data.
myMaxPDist2)[i] = std::numeric_limits<T>::infinity();
535 #if TAYLOR_SERIES_ORDER >= 1
540 data_for_parent->myNijDiag = child_data_array[0].myNijDiag;
541 data_for_parent->myNxy = 0;
542 data_for_parent->myNyx = 0;
543 data_for_parent->myNyz = 0;
544 data_for_parent->myNzy = 0;
545 data_for_parent->myNzx = 0;
546 data_for_parent->myNxz = 0;
547 #if TAYLOR_SERIES_ORDER >= 2
548 data_for_parent->myNijkDiag = child_data_array[0].myNijkDiag;
549 data_for_parent->mySumPermuteNxyz = child_data_array[0].mySumPermuteNxyz;
550 data_for_parent->my2Nxxy_Nyxx = child_data_array[0].my2Nxxy_Nyxx;
551 data_for_parent->my2Nxxz_Nzxx = child_data_array[0].my2Nxxz_Nzxx;
552 data_for_parent->my2Nyyz_Nzyy = child_data_array[0].my2Nyyz_Nzyy;
553 data_for_parent->my2Nyyx_Nxyy = child_data_array[0].my2Nyyx_Nxyy;
554 data_for_parent->my2Nzzx_Nxzz = child_data_array[0].my2Nzzx_Nxzz;
555 data_for_parent->my2Nzzy_Nyzz = child_data_array[0].my2Nzzy_Nyzz;
558 for (
int i = 1; i < nchildren; ++i)
560 data_for_parent->myNijDiag += child_data_array[i].myNijDiag;
561 #if TAYLOR_SERIES_ORDER >= 2
562 data_for_parent->myNijkDiag += child_data_array[i].myNijkDiag;
563 data_for_parent->mySumPermuteNxyz += child_data_array[i].mySumPermuteNxyz;
564 data_for_parent->my2Nxxy_Nyxx += child_data_array[i].my2Nxxy_Nyxx;
565 data_for_parent->my2Nxxz_Nzxx += child_data_array[i].my2Nxxz_Nzxx;
566 data_for_parent->my2Nyyz_Nzyy += child_data_array[i].my2Nyyz_Nzyy;
567 data_for_parent->my2Nyyx_Nxyy += child_data_array[i].my2Nyyx_Nxyy;
568 data_for_parent->my2Nzzx_Nxzz += child_data_array[i].my2Nzzx_Nxzz;
569 data_for_parent->my2Nzzy_Nyzz += child_data_array[i].my2Nzzy_Nyzz;
572 for (
int j = 0; j < 3; ++j)
573 ((
T*)¤t_box_data.
myNijDiag[j])[0] = child_data_array[0].myNijDiag[j];
574 ((
T*)¤t_box_data.
myNxy_Nyx)[0] = child_data_array[0].myNxy + child_data_array[0].myNyx;
575 ((
T*)¤t_box_data.
myNyz_Nzy)[0] = child_data_array[0].myNyz + child_data_array[0].myNzy;
576 ((
T*)¤t_box_data.
myNzx_Nxz)[0] = child_data_array[0].myNzx + child_data_array[0].myNxz;
577 for (
int j = 0; j < 3; ++j)
578 ((
T*)¤t_box_data.
myNijkDiag[j])[0] = child_data_array[0].myNijkDiag[j];
579 ((
T*)¤t_box_data.
mySumPermuteNxyz)[0] = child_data_array[0].mySumPermuteNxyz;
580 ((
T*)¤t_box_data.
my2Nxxy_Nyxx)[0] = child_data_array[0].my2Nxxy_Nyxx;
581 ((
T*)¤t_box_data.
my2Nxxz_Nzxx)[0] = child_data_array[0].my2Nxxz_Nzxx;
582 ((
T*)¤t_box_data.
my2Nyyz_Nzyy)[0] = child_data_array[0].my2Nyyz_Nzyy;
583 ((
T*)¤t_box_data.
my2Nyyx_Nxyy)[0] = child_data_array[0].my2Nyyx_Nxyy;
584 ((
T*)¤t_box_data.
my2Nzzx_Nxzz)[0] = child_data_array[0].my2Nzzx_Nxzz;
585 ((
T*)¤t_box_data.
my2Nzzy_Nyzz)[0] = child_data_array[0].my2Nzzy_Nyzz;
586 for (
int i = 1; i < nchildren; ++i)
588 for (
int j = 0; j < 3; ++j)
589 ((
T*)¤t_box_data.
myNijDiag[j])[i] = child_data_array[i].myNijDiag[j];
590 ((
T*)¤t_box_data.
myNxy_Nyx)[i] = child_data_array[i].myNxy + child_data_array[i].myNyx;
591 ((
T*)¤t_box_data.
myNyz_Nzy)[i] = child_data_array[i].myNyz + child_data_array[i].myNzy;
592 ((
T*)¤t_box_data.
myNzx_Nxz)[i] = child_data_array[i].myNzx + child_data_array[i].myNxz;
593 for (
int j = 0; j < 3; ++j)
594 ((
T*)¤t_box_data.
myNijkDiag[j])[i] = child_data_array[i].myNijkDiag[j];
595 ((
T*)¤t_box_data.
mySumPermuteNxyz)[i] = child_data_array[i].mySumPermuteNxyz;
596 ((
T*)¤t_box_data.
my2Nxxy_Nyxx)[i] = child_data_array[i].my2Nxxy_Nyxx;
597 ((
T*)¤t_box_data.
my2Nxxz_Nzxx)[i] = child_data_array[i].my2Nxxz_Nzxx;
598 ((
T*)¤t_box_data.
my2Nyyz_Nzyy)[i] = child_data_array[i].my2Nyyz_Nzyy;
599 ((
T*)¤t_box_data.
my2Nyyx_Nxyy)[i] = child_data_array[i].my2Nyyx_Nxyy;
600 ((
T*)¤t_box_data.
my2Nzzx_Nxzz)[i] = child_data_array[i].my2Nzzx_Nxzz;
601 ((
T*)¤t_box_data.
my2Nzzy_Nyzz)[i] = child_data_array[i].my2Nzzy_Nyzz;
603 for (
int i = nchildren; i < BVH_N; ++i)
606 for (
int j = 0; j < 3; ++j)
611 for (
int j = 0; j < 3; ++j)
622 for (
int i = 0; i < nchildren; ++i)
624 const LocalData &child_data = child_data_array[i];
629 data_for_parent->myNijDiag += N*displacement;
630 T Nxy = child_data.myNxy + N.
x()*displacement.y();
631 T Nyx = child_data.myNyx + N.
y()*displacement.x();
632 T Nyz = child_data.myNyz + N.
y()*displacement.z();
633 T Nzy = child_data.myNzy + N.
z()*displacement.y();
634 T Nzx = child_data.myNzx + N.
z()*displacement.x();
635 T Nxz = child_data.myNxz + N.
x()*displacement.z();
637 data_for_parent->myNxy += Nxy;
638 data_for_parent->myNyx += Nyx;
639 data_for_parent->myNyz += Nyz;
640 data_for_parent->myNzy += Nzy;
641 data_for_parent->myNzx += Nzx;
642 data_for_parent->myNxz += Nxz;
644 #if TAYLOR_SERIES_ORDER >= 2
648 data_for_parent->myNijkDiag +=
T(2)*displacement*child_data.myNijDiag + displacement*displacement*child_data.myN;
649 data_for_parent->mySumPermuteNxyz +=
dot(displacement,
UT_Vector3T<T>(Nyz+Nzy, Nzx+Nxz, Nxy+Nyx));
650 data_for_parent->my2Nxxy_Nyxx +=
651 2*(displacement.y()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxy + N.
x()*displacement.x()*displacement.y())
652 + 2*child_data.myNyx*displacement.x() + N.
y()*displacement.x()*displacement.x();
653 data_for_parent->my2Nxxz_Nzxx +=
654 2*(displacement.z()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxz + N.
x()*displacement.x()*displacement.z())
655 + 2*child_data.myNzx*displacement.x() + N.
z()*displacement.x()*displacement.x();
656 data_for_parent->my2Nyyz_Nzyy +=
657 2*(displacement.z()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyz + N.
y()*displacement.y()*displacement.z())
658 + 2*child_data.myNzy*displacement.y() + N.
z()*displacement.y()*displacement.y();
659 data_for_parent->my2Nyyx_Nxyy +=
660 2*(displacement.x()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyx + N.
y()*displacement.y()*displacement.x())
661 + 2*child_data.myNxy*displacement.y() + N.
x()*displacement.y()*displacement.y();
662 data_for_parent->my2Nzzx_Nxzz +=
663 2*(displacement.x()*child_data.myNijDiag.z() + displacement.z()*child_data.myNzx + N.
z()*displacement.z()*displacement.x())
664 + 2*child_data.myNxz*displacement.z() + N.
x()*displacement.z()*displacement.z();
665 data_for_parent->my2Nzzy_Nyzz +=
666 2*(displacement.y()*child_data.myNijDiag.z() + displacement.z()*child_data.myNzy + N.
z()*displacement.z()*displacement.y())
667 + 2*child_data.myNyz*displacement.z() + N.
y()*displacement.z()*displacement.z();
673 #if SOLID_ANGLE_DEBUG
677 #if TAYLOR_SERIES_ORDER >= 1
680 #if TAYLOR_SERIES_ORDER >= 2
691 #if SOLID_ANGLE_TIME_PRECOMPUTE
694 const PrecomputeFunctors functors(box_data, triangle_boxes.
array(), triangle_points, positions,
order);
696 LocalData local_data;
697 myTree.template traverseParallel<LocalData>(4096, functors, &local_data);
699 #if SOLID_ANGLE_TIME_PRECOMPUTE
705 template<
typename T,
typename S>
713 myTrianglePoints =
nullptr;
715 myPositions =
nullptr;
718 template<
typename T,
typename S>
721 const T accuracy_scale2 = accuracy_scale*accuracy_scale;
723 struct SolidAngleFunctors
725 const BoxData *
const myBoxData;
727 const T myAccuracyScale2;
729 const int *
const myTrianglePoints;
735 const T accuracy_scale2,
738 const int *
const triangle_points)
739 : myBoxData(box_data)
740 , myQueryPoint(query_point)
741 , myAccuracyScale2(accuracy_scale2)
743 , myPositions(positions)
744 , myTrianglePoints(triangle_points)
746 uint pre(
const int nodei,
T *data_for_parent)
const
755 const typename BoxData::Type qlength2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
761 v4uu descend_mask = (qlength2 <= maxP2*myAccuracyScale2);
762 uint descend_bitmask = _mm_movemask_ps(
V4SF(descend_mask.
vector));
763 constexpr
uint allchildbits = ((
uint(1)<<BVH_N)-1);
764 if (descend_bitmask == allchildbits)
766 *data_for_parent = 0;
773 const typename BoxData::Type qlength_m1 = SYSsqrt(qlength_m2);
780 #if TAYLOR_SERIES_ORDER >= 1
785 const typename BoxData::Type qlength_m3 = qlength_m2*qlength_m1;
792 Omega_approx += Omega_1;
793 #if TAYLOR_SERIES_ORDER >= 2
797 const typename BoxData::Type qlength_m4 = qlength_m2*qlength_m2;
811 Omega_approx += Omega_2;
819 const v4uu mask = Omega_approx.isFinite() & ~descend_mask;
820 Omega_approx = Omega_approx &
mask;
821 descend_bitmask = (~_mm_movemask_ps(
V4SF(mask.vector))) & allchildbits;
823 T sum = Omega_approx[0];
824 for (
int i = 1; i < BVH_N; ++i)
825 sum += Omega_approx[i];
826 *data_for_parent = sum;
828 return descend_bitmask;
830 void item(
const int itemi,
const int parent_nodei,
T &data_for_parent)
const
833 const int *
const cur_triangle_points = myTrianglePoints + 3*itemi;
840 SYS_FORCE_INLINE void post(
const int nodei,
const int parent_nodei,
T *data_for_parent,
const int nchildren,
const T *child_data_array,
const uint descend_bits)
const
842 T sum = (descend_bits&1) ? child_data_array[0] : 0;
843 for (
int i = 1; i < nchildren; ++i)
844 sum += ((descend_bits>>i)&1) ? child_data_array[i] : 0;
846 *data_for_parent += sum;
849 const SolidAngleFunctors functors(myData.get(), query_point, accuracy_scale2,
myOrder, myPositions, myTrianglePoints);
852 myTree.traverseVector(functors, &sum);
856 template<
typename T,
typename S>
862 memset(
this,0,
sizeof(*
this));
891 template<
typename T,
typename S>
898 , mySegmentPoints(nullptr)
900 , myPositions(nullptr)
903 template<
typename T,
typename S>
911 template<
typename T,
typename S>
914 const int *
const segment_points,
919 #if SOLID_ANGLE_DEBUG
922 UTdebugFormat(
"Building BVH for {} segments on {} points:", nsegments, npoints);
925 myNSegments = nsegments;
926 mySegmentPoints = segment_points;
928 myPositions = positions;
930 #if SOLID_ANGLE_TIME_PRECOMPUTE
936 if (nsegments < 16*1024)
938 const int *cur_segment_points = segment_points;
939 for (
int i = 0; i < nsegments; ++i, cur_segment_points += 2)
942 box.
initBounds(positions[cur_segment_points[0]]);
950 const int *cur_segment_points = segment_points +
exint(r.begin())*2;
951 for (
int i = r.begin(),
end = r.end(); i <
end; ++i, cur_segment_points += 2)
954 box.
initBounds(positions[cur_segment_points[0]]);
959 #if SOLID_ANGLE_TIME_PRECOMPUTE
964 myTree.template init<UT::BVH_Heuristic::BOX_AREA,S,2>(segment_boxes.
array(), nsegments);
965 #if SOLID_ANGLE_TIME_PRECOMPUTE
967 UTdebugFormat(
"{} s to initialize UT_BVH structure. {} nodes", time, myTree.getNumNodes());
972 const int nnodes = myTree.getNumNodes();
976 myData.reset(box_data);
1001 struct PrecomputeFunctors
1005 const int *
const mySegmentPoints;
1012 const int *segment_points,
1015 : myBoxData(box_data)
1016 , mySegmentBoxes(segment_boxes)
1017 , mySegmentPoints(segment_points)
1018 , myPositions(positions)
1021 constexpr
SYS_FORCE_INLINE bool pre(
const int nodei, LocalData *data_for_parent)
const
1025 void item(
const int itemi,
const int parent_nodei, LocalData &data_for_parent)
const
1028 const int *
const cur_segment_points = mySegmentPoints + 2*itemi;
1033 const UT::Box<S,2> &segment_box = mySegmentBoxes[itemi];
1034 data_for_parent.myBox = segment_box;
1039 const T length = SYSsqrt(length2);
1041 data_for_parent.myAverageP = P;
1042 data_for_parent.myLengthP = P*
length;
1043 data_for_parent.myN =
N;
1044 #if SOLID_ANGLE_DEBUG
1046 UTdebugFormat(
"Triangle {}: P = {}; N = {}; length = {}", itemi, P,
N, length);
1050 data_for_parent.myLength =
length;
1057 data_for_parent.myNijDiag =
T(0);
1058 data_for_parent.myNxy = 0; data_for_parent.myNyx = 0;
1066 data_for_parent.myNijkDiag =
T(0);
1067 data_for_parent.my2Nxxy_Nyxx = 0;
1068 data_for_parent.my2Nyyx_Nxyy = 0;
1072 T integral_xx = ab.
x()*ab.
x()/
T(12);
1073 T integral_xy = ab.
x()*ab.
y()/
T(12);
1074 T integral_yy = ab.
y()*ab.
y()/
T(12);
1075 data_for_parent.myNijkDiag.assign(integral_xx*
N.x(),integral_yy*
N.y());
1076 T Nxxy =
N.x()*integral_xy;
1077 T Nyxx =
N.y()*integral_xx;
1078 T Nyyx =
N.y()*integral_xy;
1079 T Nxyy =
N.x()*integral_yy;
1080 data_for_parent.my2Nxxy_Nyxx = 2*Nxxy + Nyxx;
1081 data_for_parent.my2Nyyx_Nxyy = 2*Nyyx + Nxyy;
1082 #if SOLID_ANGLE_DEBUG
1083 UTdebugFormat(
" integral_xx = {}; yy = {}", integral_xx, integral_yy);
1088 void post(
const int nodei,
const int parent_nodei, LocalData *data_for_parent,
const int nchildren,
const LocalData *child_data_array)
const
1093 BoxData ¤t_box_data = myBoxData[nodei];
1096 ((
T*)¤t_box_data.
myN[0])[0] = N[0];
1097 ((
T*)¤t_box_data.
myN[1])[0] = N[1];
1099 T length = child_data_array[0].myLength;
1101 ((
T*)¤t_box_data.
myAverageP[0])[0] = local_P[0];
1102 ((
T*)¤t_box_data.
myAverageP[1])[0] = local_P[1];
1103 for (
int i = 1; i < nchildren; ++i)
1107 ((
T*)¤t_box_data.
myN[0])[i] = local_N[0];
1108 ((
T*)¤t_box_data.
myN[1])[i] = local_N[1];
1109 lengthP += child_data_array[i].myLengthP;
1110 length += child_data_array[i].myLength;
1112 ((
T*)¤t_box_data.
myAverageP[0])[i] = local_P[0];
1113 ((
T*)¤t_box_data.
myAverageP[1])[i] = local_P[1];
1115 for (
int i = nchildren; i < BVH_N; ++i)
1118 ((
T*)¤t_box_data.
myN[0])[i] = 0;
1119 ((
T*)¤t_box_data.
myN[1])[i] = 0;
1123 data_for_parent->myN =
N;
1124 data_for_parent->myLengthP = lengthP;
1125 data_for_parent->myLength =
length;
1128 for (
int i = 1; i < nchildren; ++i)
1129 box.
combine(child_data_array[i].myBox);
1134 averageP = lengthP/
length;
1137 data_for_parent->myAverageP = averageP;
1139 data_for_parent->myBox =
box;
1141 for (
int i = 0; i < nchildren; ++i)
1143 const UT::Box<S,2> &local_box(child_data_array[i].myBox);
1148 for (
int i = nchildren; i < BVH_N; ++i)
1152 ((
T*)¤t_box_data.
myMaxPDist2)[i] = std::numeric_limits<T>::infinity();
1159 data_for_parent->myNijDiag = child_data_array[0].myNijDiag;
1160 data_for_parent->myNxy = 0;
1161 data_for_parent->myNyx = 0;
1162 data_for_parent->myNijkDiag = child_data_array[0].myNijkDiag;
1163 data_for_parent->my2Nxxy_Nyxx = child_data_array[0].my2Nxxy_Nyxx;
1164 data_for_parent->my2Nyyx_Nxyy = child_data_array[0].my2Nyyx_Nxyy;
1166 for (
int i = 1; i < nchildren; ++i)
1168 data_for_parent->myNijDiag += child_data_array[i].myNijDiag;
1169 data_for_parent->myNijkDiag += child_data_array[i].myNijkDiag;
1170 data_for_parent->my2Nxxy_Nyxx += child_data_array[i].my2Nxxy_Nyxx;
1171 data_for_parent->my2Nyyx_Nxyy += child_data_array[i].my2Nyyx_Nxyy;
1173 for (
int j = 0; j < 2; ++j)
1174 ((
T*)¤t_box_data.
myNijDiag[j])[0] = child_data_array[0].myNijDiag[j];
1175 ((
T*)¤t_box_data.
myNxy_Nyx)[0] = child_data_array[0].myNxy + child_data_array[0].myNyx;
1176 for (
int j = 0; j < 2; ++j)
1177 ((
T*)¤t_box_data.
myNijkDiag[j])[0] = child_data_array[0].myNijkDiag[j];
1178 ((
T*)¤t_box_data.
my2Nxxy_Nyxx)[0] = child_data_array[0].my2Nxxy_Nyxx;
1179 ((
T*)¤t_box_data.
my2Nyyx_Nxyy)[0] = child_data_array[0].my2Nyyx_Nxyy;
1180 for (
int i = 1; i < nchildren; ++i)
1182 for (
int j = 0; j < 2; ++j)
1183 ((
T*)¤t_box_data.
myNijDiag[j])[i] = child_data_array[i].myNijDiag[j];
1184 ((
T*)¤t_box_data.
myNxy_Nyx)[i] = child_data_array[i].myNxy + child_data_array[i].myNyx;
1185 for (
int j = 0; j < 2; ++j)
1186 ((
T*)¤t_box_data.
myNijkDiag[j])[i] = child_data_array[i].myNijkDiag[j];
1187 ((
T*)¤t_box_data.
my2Nxxy_Nyxx)[i] = child_data_array[i].my2Nxxy_Nyxx;
1188 ((
T*)¤t_box_data.
my2Nyyx_Nxyy)[i] = child_data_array[i].my2Nyyx_Nxyy;
1190 for (
int i = nchildren; i < BVH_N; ++i)
1193 for (
int j = 0; j < 2; ++j)
1196 for (
int j = 0; j < 2; ++j)
1202 for (
int i = 0; i < nchildren; ++i)
1204 const LocalData &child_data = child_data_array[i];
1209 data_for_parent->myNijDiag += N*displacement;
1210 T Nxy = child_data.myNxy + N.
x()*displacement.y();
1211 T Nyx = child_data.myNyx + N.
y()*displacement.x();
1213 data_for_parent->myNxy += Nxy;
1214 data_for_parent->myNyx += Nyx;
1219 data_for_parent->myNijkDiag +=
T(2)*displacement*child_data.myNijDiag + displacement*displacement*child_data.myN;
1220 data_for_parent->my2Nxxy_Nyxx +=
1221 2*(displacement.y()*child_data.myNijDiag.x() + displacement.x()*child_data.myNxy + N.
x()*displacement.x()*displacement.y())
1222 + 2*child_data.myNyx*displacement.x() + N.
y()*displacement.x()*displacement.x();
1223 data_for_parent->my2Nyyx_Nxyy +=
1224 2*(displacement.x()*child_data.myNijDiag.y() + displacement.y()*child_data.myNyx + N.
y()*displacement.y()*displacement.x())
1225 + 2*child_data.myNxy*displacement.y() + N.
x()*displacement.y()*displacement.y();
1229 #if SOLID_ANGLE_DEBUG
1241 #if SOLID_ANGLE_TIME_PRECOMPUTE
1244 const PrecomputeFunctors functors(box_data, segment_boxes.
array(), segment_points, positions,
order);
1246 LocalData local_data;
1247 myTree.template traverseParallel<LocalData>(4096, functors, &local_data);
1249 #if SOLID_ANGLE_TIME_PRECOMPUTE
1250 time = timer.
stop();
1255 template<
typename T,
typename S>
1263 mySegmentPoints =
nullptr;
1265 myPositions =
nullptr;
1268 template<
typename T,
typename S>
1271 const T accuracy_scale2 = accuracy_scale*accuracy_scale;
1273 struct AngleFunctors
1275 const BoxData *
const myBoxData;
1277 const T myAccuracyScale2;
1279 const int *
const mySegmentPoints;
1283 const BoxData *
const box_data,
1285 const T accuracy_scale2,
1288 const int *
const segment_points)
1289 : myBoxData(box_data)
1290 , myQueryPoint(query_point)
1291 , myAccuracyScale2(accuracy_scale2)
1293 , myPositions(positions)
1294 , mySegmentPoints(segment_points)
1296 uint pre(
const int nodei,
T *data_for_parent)
const
1304 const typename BoxData::Type qlength2 = q[0]*q[0] + q[1]*q[1];
1310 v4uu descend_mask = (qlength2 <= maxP2*myAccuracyScale2);
1311 uint descend_bitmask = _mm_movemask_ps(
V4SF(descend_mask.
vector));
1312 constexpr
uint allchildbits = ((
uint(1)<<BVH_N)-1);
1313 if (descend_bitmask == allchildbits)
1315 *data_for_parent = 0;
1316 return allchildbits;
1322 const typename BoxData::Type qlength_m1 = SYSsqrt(qlength_m2);
1337 Omega_approx += Omega_1;
1341 const typename BoxData::Type qlength_m3 = qlength_m2*qlength_m1;
1353 Omega_approx += Omega_2;
1359 const v4uu mask = Omega_approx.isFinite() & ~descend_mask;
1360 Omega_approx = Omega_approx &
mask;
1361 descend_bitmask = (~_mm_movemask_ps(
V4SF(mask.vector))) & allchildbits;
1363 T sum = Omega_approx[0];
1364 for (
int i = 1; i < BVH_N; ++i)
1365 sum += Omega_approx[i];
1366 *data_for_parent = sum;
1368 return descend_bitmask;
1370 void item(
const int itemi,
const int parent_nodei,
T &data_for_parent)
const
1373 const int *
const cur_segment_points = mySegmentPoints + 2*itemi;
1379 SYS_FORCE_INLINE void post(
const int nodei,
const int parent_nodei,
T *data_for_parent,
const int nchildren,
const T *child_data_array,
const uint descend_bits)
const
1381 T sum = (descend_bits&1) ? child_data_array[0] : 0;
1382 for (
int i = 1; i < nchildren; ++i)
1383 sum += ((descend_bits>>i)&1) ? child_data_array[i] : 0;
1385 *data_for_parent += sum;
1388 const AngleFunctors functors(myData.get(), query_point, accuracy_scale2,
myOrder, myPositions, mySegmentPoints);
1391 myTree.traverseVector(functors, &sum);
constexpr SYS_FORCE_INLINE T length2() const noexcept
UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
SYS_FORCE_INLINE UT_FixedVector< T, NAXES > getMax() const noexcept
GLboolean GLboolean GLboolean b
typename SYS_SelectType< UT_FixedVector< S, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< S, float >::value >::type SType
GLuint GLdouble GLdouble GLint GLint order
Axis-aligned bounding box (AABB).
#define SYS_STATIC_ASSERT_MSG(expr, msg)
SYS_FORCE_INLINE void initBounds() noexcept
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)
GT_API const UT_StringHolder time
void setSizeNoInit(exint newsize)
UT_FixedVector< Type, 2 > myNijDiag
constexpr SYS_FORCE_INLINE T & z() noexcept
UT_Vector3T< T > maxvec() const
UT_FixedVector< Type, 3 > myNijkDiag
typename SYS_SelectType< UT_FixedVector< T, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< T, float >::value >::type Type
constexpr SYS_FORCE_INLINE T length() const noexcept
#define UT_ASSERT_MSG_P(ZZ,...)
void UTparallelFor(const Range &range, const Body &body, const int subscribe_ratio=2, const int min_grain_size=1, const bool force_use_task_scope=false)
~UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
GLenum GLsizei GLsizei GLint * values
GLdouble GLdouble GLdouble GLdouble q
typename SYS_SelectType< UT_FixedVector< S, BVH_N >, v4uf, BVH_N==4 &&SYS_IsSame< S, float >::value >::type SType
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
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.
GLboolean GLboolean GLboolean GLboolean a
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...
GLdouble GLdouble GLdouble z
void enlargeBounds(const UT_Vector3T< T > &min, const UT_Vector3T< T > &max)
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
T computeAngle(const UT_Vector2T< T > &query_point, const T accuracy_scale=T(2.0)) const
SYS_FORCE_INLINE void enlargeBounds(const UT_FixedVector< T, NAXES > &pt) noexcept
T UTsignedAngleSegment(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, const UT_Vector2T< T > &query)
UT_FixedVector< Type, 2 > myNijkDiag
GLuint GLsizei GLsizei * length
constexpr SYS_FORCE_INLINE T length2() const noexcept
UT_SubtendedAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
GLsizei const GLint box[]
SYS_FORCE_INLINE UT_FixedVector< T, NAXES > getMin() const noexcept
UT_FixedVector< Type, 3 > myAverageP
Centre of mass of the mesh surface in this box.
UT_Vector3T< T > minvec() const
GA_API const UT_StringHolder N
void clear()
Frees myTree and myData, and clears the rest.
UT_FixedVector< Type, 3 > myNijDiag
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
SYS_FORCE_INLINE void combine(const Box< T, NAXES > &src) noexcept
T UTsignedSolidAngleTri(const UT_Vector3T< T > &a, const UT_Vector3T< T > &b, const UT_Vector3T< T > &c, const UT_Vector3T< T > &query)
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
UT_FixedVector< Type, 2 > myN
Unnormalized, area-weighted normal of the mesh in this box.
GA_API const UT_StringHolder area
#define UTdebugFormat(...)
constexpr SYS_FORCE_INLINE T & x() noexcept