HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SOP_WindingNumber.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  * This SOP computes the winding number induced by a mesh at the specified points.
27  */
28 
29 #include "SOP_WindingNumber.h"
30 
31 // This is an automatically generated header file based on theDsFile, below,
32 // to provide SOP_WindingNumberParms, an easy way to access parameter values from
33 // SOP_WindingNumberVerb::cook with the correct type.
35 
36 #include "UT_SolidAngle.h"
37 
38 #include <SOP/SOP_NodeVerb.h>
39 #include <GU/GU_Detail.h>
40 #include <GEO/GEO_Curve.h>
41 #include <GEO/GEO_PrimMesh.h>
42 #include <GEO/GEO_PrimPolySoup.h>
43 #include <GEO/GEO_PrimSphere.h>
45 #include <GEO/GEO_TPSurf.h>
46 #include <GA/GA_Handle.h>
47 #include <GA/GA_SplittableRange.h>
48 #include <GA/GA_Types.h>
49 #include <OP/OP_Operator.h>
50 #include <OP/OP_OperatorTable.h>
52 #include <UT/UT_DSOVersion.h>
53 #include <UT/UT_Interrupt.h>
54 #include <UT/UT_ParallelUtil.h>
55 #include <UT/UT_StringHolder.h>
56 #include <UT/UT_UniquePtr.h>
57 #include <SYS/SYS_Math.h>
58 
59 // UT_SolidAngle doesn't currently have a special case to handle quads correctly
60 // as bilinear patches; it always treats them as two triangles, which can result
61 // in a difference of about 1 between points that are just inside or just
62 // outside. If you want the accurate version to treat them as triangles too,
63 // (e.g. for testing the approximation error of UT_SolidAngle), set this to 0.
64 #define QUADS_AS_BILINEAR_PATCHES 1
65 
66 
67 using namespace UT::Literal;
68 using namespace HDK_Sample;
69 
70 //
71 // Help is stored in a "wiki" style text file. This text file should be copied
72 // to $HOUDINI_PATH/help/nodes/sop/hdk_windingnumber.txt
73 //
74 // See the sample_install.sh file for an example.
75 //
76 
77 // Forward declarations of file-scope functions used before they're defined
78 
79 static void
80 sopAccumulateTriangles(
81  const GA_Detail *const mesh_geo,
82  const GA_Offset primoff,
83  const UT_Array<int> &ptmap,
84  UT_Array<int> &triangle_points);
85 
86 static void
87 sopAccumulateSegments(
88  const GA_Detail *const mesh_geo,
89  const GA_Offset primoff,
90  const UT_Array<int> &ptmap,
91  UT_Array<int> &segment_points);
92 
93 // This gives easy access to the enum types for this operator's parameter
94 using namespace SOP_WindingNumberEnums;
95 
96 //******************************************************************************
97 //* Setup *
98 //******************************************************************************
99 
100 typedef int64 GA_DataId;
101 #define GA_INVALID_DATAID GA_DataId(-1)
102 
103 /// This class is for caching data between cooks, e.g. so that we don't have to
104 /// rebuild the UT_SolidAngle tree on every cook if the input hasn't changed.
106 {
107 public:
109  , mySolidAngleTree()
110  , mySubtendedAngleTree()
111  , myTopologyDataID(GA_INVALID_DATAID)
112  , myPrimitiveListDataID(GA_INVALID_DATAID)
113  , myPDataID(GA_INVALID_DATAID)
114  , myApproxOrder(-1)
115  , myAxis0(-1)
116  , myHadGroup(false)
117  , myGroupString()
118  , myUniqueId(-1)
119  , myMetaCacheCount(-1)
120  {}
122 
123  void update3D(const GA_Detail &mesh_geo, const GA_PrimitiveGroup *prim_group, const UT_StringHolder &group_string,const int approx_order)
124  {
125  const GA_DataId topology_data_id = mesh_geo.getTopology().getDataId();
126  const GA_DataId primitive_list_data_id = mesh_geo.getPrimitiveList().getDataId();
127  const GA_DataId P_data_id = mesh_geo.getP()->getDataId();
128  const bool has_group = (prim_group != nullptr);
129  if (mySubtendedAngleTree.isClear() &&
130  topology_data_id == myTopologyDataID &&
131  primitive_list_data_id == myPrimitiveListDataID &&
132  P_data_id == myPDataID &&
133  approx_order == myApproxOrder &&
134  has_group == myHadGroup &&
135  (!has_group || (
136  group_string == myGroupString &&
137  mesh_geo.getUniqueId() == myUniqueId &&
138  mesh_geo.getMetaCacheCount() == myMetaCacheCount)))
139  {
140  return;
141  }
142  mySubtendedAngleTree.clear();
143  myPositions2D.clear();
144 
145  UT_AutoInterrupt boss("Constructing Solid Angle Tree");
146 
147  myTopologyDataID = topology_data_id;
148  myPrimitiveListDataID = primitive_list_data_id;
149  myPDataID = P_data_id;
150  myApproxOrder = approx_order;
151  myHadGroup = has_group;
152  myGroupString = group_string;
153  myUniqueId = mesh_geo.getUniqueId();
154  myMetaCacheCount = mesh_geo.getMetaCacheCount();
155 
156  GA_Size nprimitives = prim_group ? prim_group->entries() : mesh_geo.getNumPrimitives();
157  myTrianglePoints.setSize(0);
158  myTrianglePoints.setCapacity(3*nprimitives);
159  myPositions3D.setSize(0);
160  myPositions3D.setCapacity(0);
161  UT_Array<int> ptmap;
162  if (!prim_group)
163  {
164  // Copy all point positions
165  mesh_geo.getAttributeAsArray(mesh_geo.getP(), mesh_geo.getPointRange(), myPositions3D);
166  }
167  else
168  {
169  // Find the used points
170  GA_PointGroup mesh_ptgroup(mesh_geo);
171  mesh_ptgroup.combine(prim_group);
172  myPositions3D.setSizeNoInit(mesh_ptgroup.entries());
173  ptmap.setSizeNoInit(mesh_geo.getNumPointOffsets());
174  int ptnum = 0;
175  mesh_geo.forEachPoint(&mesh_ptgroup, [&mesh_geo,this,&ptmap,&ptnum](const GA_Offset ptoff)
176  {
177  ptmap[ptoff] = ptnum;
178  myPositions3D[ptnum] = mesh_geo.getPos3(ptoff);
179  ++ptnum;
180  });
181  }
182 
184  GA_Offset end;
185  for (GA_Iterator it(mesh_geo.getPrimitiveRange(prim_group)); it.blockAdvance(start,end); )
186  {
187  for (GA_Offset primoff = start; primoff < end; ++primoff)
188  {
189  sopAccumulateTriangles(&mesh_geo, primoff, ptmap, myTrianglePoints);
190  }
191  }
192 
193  mySolidAngleTree.init(myTrianglePoints.size()/3, myTrianglePoints.array(), myPositions3D.size(), myPositions3D.array(), approx_order);
194  }
195 
196  void update2D(
197  const GA_Detail &mesh_geo,
198  const GA_PrimitiveGroup *prim_group,
199  const UT_StringHolder &group_string,
200  const int approx_order,
201  const int axis0, const int axis1)
202  {
203  const GA_DataId topology_data_id = mesh_geo.getTopology().getDataId();
204  const GA_DataId primitive_list_data_id = mesh_geo.getPrimitiveList().getDataId();
205  const GA_DataId P_data_id = mesh_geo.getP()->getDataId();
206  const bool has_group = (prim_group != nullptr);
207  if (mySolidAngleTree.isClear() &&
208  topology_data_id == myTopologyDataID &&
209  primitive_list_data_id == myPrimitiveListDataID &&
210  P_data_id == myPDataID &&
211  approx_order == myApproxOrder &&
212  axis0 == myAxis0 &&
213  has_group == myHadGroup &&
214  (!has_group || (
215  group_string == myGroupString &&
216  mesh_geo.getUniqueId() == myUniqueId &&
217  mesh_geo.getMetaCacheCount() == myMetaCacheCount)))
218  {
219  return;
220  }
221  mySolidAngleTree.clear();
222  myPositions3D.clear();
223 
224  UT_AutoInterrupt boss("Constructing Solid Angle Tree");
225 
226  myTopologyDataID = topology_data_id;
227  myPrimitiveListDataID = primitive_list_data_id;
228  myPDataID = P_data_id;
229  myApproxOrder = approx_order;
230  myAxis0 = axis0;
231  myHadGroup = has_group;
232  myGroupString = group_string;
233  myUniqueId = mesh_geo.getUniqueId();
234  myMetaCacheCount = mesh_geo.getMetaCacheCount();
235 
236  GA_Size nprimitives = prim_group ? prim_group->entries() : mesh_geo.getNumPrimitives();
237  // NOTE: myTrianglePoints is actually segment points
238  myTrianglePoints.setSize(0);
239  myTrianglePoints.setCapacity(2*nprimitives);
240  myPositions2D.setSize(0);
241  myPositions2D.setCapacity(0);
242  UT_Array<int> ptmap;
243  if (!prim_group)
244  {
245  // Copy all point positions
246  myPositions2D.setSizeNoInit(mesh_geo.getNumPoints());
247  int ptnum = 0;
248  mesh_geo.forEachPoint([&mesh_geo,this,&ptnum,axis0,axis1](const GA_Offset ptoff)
249  {
250  UT_Vector3 pos = mesh_geo.getPos3(ptoff);
251  myPositions2D[ptnum] = UT_Vector2(pos[axis0],pos[axis1]);
252  ++ptnum;
253  });
254  }
255  else
256  {
257  // Find the used points
258  GA_PointGroup mesh_ptgroup(mesh_geo);
259  mesh_ptgroup.combine(prim_group);
260  myPositions2D.setSizeNoInit(mesh_ptgroup.entries());
261  ptmap.setSizeNoInit(mesh_geo.getNumPointOffsets());
262  int ptnum = 0;
263  mesh_geo.forEachPoint(&mesh_ptgroup, [&mesh_geo,this,&ptmap,&ptnum,axis0,axis1](const GA_Offset ptoff)
264  {
265  ptmap[ptoff] = ptnum;
266  UT_Vector3 pos = mesh_geo.getPos3(ptoff);
267  myPositions2D[ptnum] = UT_Vector2(pos[axis0],pos[axis1]);
268  ++ptnum;
269  });
270  }
271 
273  GA_Offset end;
274  for (GA_Iterator it(mesh_geo.getPrimitiveRange(prim_group)); it.blockAdvance(start,end); )
275  {
276  for (GA_Offset primoff = start; primoff < end; ++primoff)
277  {
278  sopAccumulateSegments(&mesh_geo, primoff, ptmap, myTrianglePoints);
279  }
280  }
281 
282  mySubtendedAngleTree.init(myTrianglePoints.size()/2, myTrianglePoints.array(), myPositions2D.size(), myPositions2D.array(), approx_order);
283  }
284 
285  void clear()
286  {
287  mySolidAngleTree.clear();
288  mySubtendedAngleTree.clear();
289  myTrianglePoints.setCapacity(0);
290  myPositions2D.setCapacity(0);
291  myPositions3D.setCapacity(0);
292  myTopologyDataID = GA_INVALID_DATAID;
293  myPrimitiveListDataID = GA_INVALID_DATAID;
294  myPDataID = GA_INVALID_DATAID;
295  myApproxOrder = -1;
296  myHadGroup = false;
297  myGroupString.clear();
298  myUniqueId = -1;
299  myMetaCacheCount = -1;
300  }
301 
311  int myAxis0;
316 };
317 
318 
320 {
321 public:
322  virtual SOP_NodeParms *allocParms() const { return new SOP_WindingNumberParms(); }
323  virtual SOP_NodeCache *allocCache() const { return new SOP_WindingNumberCache(); }
324  virtual UT_StringHolder name() const { return theSOPTypeName; }
325 
326  /// This SOP wouldn't get any benefit from the results of the previous cook,
327  /// (except for what's stored in SOP_WindingNumberCache), and it would
328  /// always duplicate its input at the start of the cook anyway, so it might
329  /// as well use COOK_INPLACE, (which always either steals the first input's
330  /// detail for the output detail or duplicates it into the output detail),
331  /// instead of COOK_GENERIC, (which would let us have an output detail
332  /// that's separate from the input detail and might be the same output
333  /// detail as on a previous cook).
334  virtual CookMode cookMode(const SOP_NodeParms *parms) const { return COOK_INPLACE; }
335 
336  virtual void cook(const CookParms &cookparms) const;
337 
338  /// This is the internal name of the SOP type.
339  /// It isn't allowed to be the same as any other SOP's type name.
341 
342  /// This static data member automatically registers
343  /// this verb class at library load time.
345 
346  /// This is the parameter interface string, below.
347  static const char *const theDsFile;
348 };
349 
350 // The static member variable definitions have to be outside the class definition.
351 // The declarations are inside the class.
352 const UT_StringHolder SOP_WindingNumberVerb::theSOPTypeName("hdk_windingnumber"_sh);
354 
355 PRM_Template *SOP_WindingNumber::buildTemplates()
356 {
357  static PRM_TemplateBuilder templ("SOP_WindingNumber.C"_sh, SOP_WindingNumberVerb::theDsFile);
358  if (templ.justBuilt())
359  {
360  templ.setChoiceListPtr("querypoints"_sh, &SOP_Node::pointGroupMenu);
361  templ.setChoiceListPtr("meshprims"_sh, &SOP_Node::primGroupMenu);
362  }
363  return templ.templates();
364 }
365 
366 const SOP_NodeVerb *SOP_WindingNumber::cookVerb() const
367 {
368  return SOP_WindingNumberVerb::theVerb.get();
369 }
370 
371 /// newSopOperator is the hook that Houdini grabs from this dll
372 /// and invokes to register the SOP. In this case, we add ourselves
373 /// to the specified operator table.
375 {
376  table->addOperator(new OP_Operator(
377  SOP_WindingNumberVerb::theSOPTypeName, // Internal name
378  "Winding Number", // UI name
379  SOP_WindingNumber::myConstructor, // How to build the SOP
380  SOP_WindingNumber::buildTemplates(), // My parameters
381  2, // Min # of sources
382  2, // Max # of sources
383  nullptr,// Custom local variables (none)
384  0)); // No flags: not a generator, no merge input, not an output
385 }
386 
387 //******************************************************************************
388 //* Parameters *
389 //******************************************************************************
390 
391 /// This is a multi-line raw string specifying the parameter interface for this SOP.
392 const char *const SOP_WindingNumberVerb::theDsFile = R"THEDSFILE(
393 {
394  name parameters
395  parm {
396  name "querypoints"
397  cppname "QueryPoints"
398  label "Query Points"
399  type string
400  default { "" }
401  parmtag { "script_action" "import soputils\nkwargs['geometrytype'] = (hou.geometryType.Points,)\nkwargs['inputindex'] = 0\nsoputils.selectGroupParm(kwargs)" }
402  parmtag { "script_action_help" "Select geometry from an available viewport.\nShift-click to turn on Select Groups." }
403  parmtag { "script_action_icon" "BUTTONS_reselect" }
404  }
405  parm {
406  name "meshprims"
407  cppname "MeshPrims"
408  label "Mesh Primitives"
409  type string
410  default { "" }
411  parmtag { "script_action" "import soputils\nkwargs['geometrytype'] = (hou.geometryType.Primitives,)\nkwargs['inputindex'] = 1\nsoputils.selectGroupParm(kwargs)" }
412  parmtag { "script_action_help" "Select geometry from an available viewport.\nShift-click to turn on Select Groups." }
413  parmtag { "script_action_icon" "BUTTONS_reselect" }
414  parmtag { "sop_input" "1" }
415  }
416  parm {
417  name "type"
418  label "Winding Number Type"
419  type ordinal
420  default { "0" }
421  menu {
422  "xyz" "3D"
423  "xy" "2D in XY Plane"
424  "yz" "2D in YZ Plane"
425  "zx" "2D in ZX Plane"
426  }
427  }
428  parm {
429  name "attrib"
430  label "Attribute Name"
431  type string
432  default { "winding_number" }
433  }
434  parm {
435  name "assolidangle"
436  cppname "AsSolidAngle"
437  label "Scale to Solid Angle"
438  type toggle
439  default { "0" }
440  joinnext
441  }
442  parm {
443  name "negate"
444  cppname "Negate"
445  label "Negate Value (Reverse)"
446  type toggle
447  default { "0" }
448  }
449  parm {
450  name "fullaccuracy"
451  cppname "FullAccuracy"
452  label "Full Accuracy (Slow)"
453  type toggle
454  default { "0" }
455  }
456  parm {
457  name "accuracyscale"
458  cppname "AccuracyScale"
459  label "Accuracy Scale"
460  type float
461  default { "2" }
462  range { 1! 20 }
463  disablewhen "{ fullaccuracy == 1 }"
464  }
465 }
466 )THEDSFILE";
467 
468 //******************************************************************************
469 //* Cooking *
470 //******************************************************************************
471 
472 static void
473 sopSumContributions3D(
474  double *sum,
475  const UT_Vector3D &query_point,
476  const GEO_Detail *const mesh_geo,
477  const GA_OffsetList &primoffs,
478  const exint start, const exint end)
479 {
480  if (end-start > 1024)
481  {
482  exint mid = (start+end)>>1;
483  double sum0;
484  double sum1;
485  UTparallelInvoke(true, [&] {
486  sopSumContributions3D(&sum0, query_point, mesh_geo, primoffs, start, mid);
487  }, [&] {
488  sopSumContributions3D(&sum1, query_point, mesh_geo, primoffs, mid, end);
489  });
490  *sum = sum0 + sum1;
491  return;
492  }
493 
494  double local_sum = 0;
495  for (exint arrayi = start; arrayi < end; ++arrayi)
496  {
497  GA_Offset primoff = primoffs(arrayi);
498  int primtype = mesh_geo->getPrimitiveTypeId(primoff);
499  if (primtype == GA_PRIMPOLY)
500  {
501  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
502  const GA_Size n = vertices.size();
503  const bool closed = vertices.getExtraFlag();
504  if (n < 3 || !closed)
505  continue;
506 
507  if (n == 3)
508  {
509  const GA_Offset pt0 = mesh_geo->vertexPoint(vertices(0));
510  const GA_Offset pt1 = mesh_geo->vertexPoint(vertices(1));
511  const GA_Offset pt2 = mesh_geo->vertexPoint(vertices(2));
512  const UT_Vector3D p0 = mesh_geo->getPos3(pt0);
513  const UT_Vector3D p1 = mesh_geo->getPos3(pt1);
514  const UT_Vector3D p2 = mesh_geo->getPos3(pt2);
515  const double poly_sum = UTsignedSolidAngleTri(p0, p1, p2, query_point);
516  local_sum += poly_sum;
517  continue;
518  }
519  if (n == 4)
520  {
521  // Special case for quads, to divide the correct diagonal
522  // in the case where the point is inside the convex hull.
523  const GA_Offset pt0 = mesh_geo->vertexPoint(vertices(0));
524  const GA_Offset pt1 = mesh_geo->vertexPoint(vertices(1));
525  const GA_Offset pt2 = mesh_geo->vertexPoint(vertices(2));
526  const GA_Offset pt3 = mesh_geo->vertexPoint(vertices(3));
527  const UT_Vector3D p0 = mesh_geo->getPos3(pt0);
528  const UT_Vector3D p1 = mesh_geo->getPos3(pt1);
529  const UT_Vector3D p2 = mesh_geo->getPos3(pt2);
530  const UT_Vector3D p3 = mesh_geo->getPos3(pt3);
531 #if QUADS_AS_BILINEAR_PATCHES
532  const double poly_sum = UTsignedSolidAngleQuad(p0, p1, p2, p3, query_point);
533 #else
534  const double poly_sum =
535  UTsignedSolidAngleTri(p0, p1, p2, query_point) +
536  UTsignedSolidAngleTri(p0, p2, p3, query_point);
537 #endif
538  local_sum += poly_sum;
539  continue;
540  }
541  // A triangle fan suffices, even if the polygon is non-convex,
542  // because the contributions in the opposite direction will
543  // partly cancel out the ones in the other direction,
544  // in just the right amount.
545  double poly_sum = 0;
546  const UT_Vector3D p0 = mesh_geo->getPos3(mesh_geo->vertexPoint(vertices(0)));
547  UT_Vector3D prev = mesh_geo->getPos3(mesh_geo->vertexPoint(vertices(1)));
548  for (GA_Size i = 2; i < n; ++i)
549  {
550  const UT_Vector3D next = mesh_geo->getPos3(mesh_geo->vertexPoint(vertices(i)));
551  poly_sum += UTsignedSolidAngleTri(p0, prev, next, query_point);
552  prev = next;
553  }
554  local_sum += poly_sum;
555  }
556  else if (primtype == GA_PRIMTETRAHEDRON)
557  {
558  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
559  const GEO_PrimTetrahedron tet(SYSconst_cast(mesh_geo), primoff, vertices);
560 
561  double tet_sum = 0;
562  for (int i = 0; i < 4; ++i)
563  {
564  // Ignore shared tet faces. They would contribute exactly opposite amounts.
565  if (tet.isFaceShared(i))
566  continue;
567 
568  const int *face_indices = GEO_PrimTetrahedron::fastFaceIndices(i);
569  const GA_Offset pt0 = mesh_geo->vertexPoint(vertices(face_indices[0]));
570  const GA_Offset pt1 = mesh_geo->vertexPoint(vertices(face_indices[1]));
571  const GA_Offset pt2 = mesh_geo->vertexPoint(vertices(face_indices[2]));
572  const UT_Vector3D a = mesh_geo->getPos3(pt0);
573  const UT_Vector3D b = mesh_geo->getPos3(pt1);
574  const UT_Vector3D c = mesh_geo->getPos3(pt2);
575 
576  tet_sum += UTsignedSolidAngleTri(a, b, c, query_point);
577  }
578  local_sum += tet_sum;
579  }
580  else if (primtype == GA_PRIMPOLYSOUP)
581  {
582  double soup_sum = 0;
583  const GEO_PrimPolySoup *soup = UTverify_cast<const GEO_PrimPolySoup *>(mesh_geo->getPrimitive(primoff));
584  for (GEO_PrimPolySoup::PolygonIterator poly(*soup); !poly.atEnd(); ++poly)
585  {
586  GA_Size n = poly.nvertices();
587  if (n < 3)
588  continue;
589 
590  if (n == 3)
591  {
592  const UT_Vector3D p0 = poly.getPos3(0);
593  const UT_Vector3D p1 = poly.getPos3(1);
594  const UT_Vector3D p2 = poly.getPos3(2);
595  const double poly_sum = UTsignedSolidAngleTri(p0, p1, p2, query_point);
596  local_sum += poly_sum;
597  continue;
598  }
599  if (n == 4)
600  {
601  // Special case for quads, to divide the correct diagonal
602  // in the case where the point is inside the convex hull.
603  const UT_Vector3D p0 = poly.getPos3(0);
604  const UT_Vector3D p1 = poly.getPos3(1);
605  const UT_Vector3D p2 = poly.getPos3(2);
606  const UT_Vector3D p3 = poly.getPos3(3);
607 #if QUADS_AS_BILINEAR_PATCHES
608  const double poly_sum = UTsignedSolidAngleQuad(p0, p1, p2, p3, query_point);
609 #else
610  const double poly_sum =
611  UTsignedSolidAngleTri(p0, p1, p2, query_point) +
612  UTsignedSolidAngleTri(p0, p2, p3, query_point);
613 #endif
614  local_sum += poly_sum;
615  continue;
616  }
617  // A triangle fan suffices, even if the polygon is non-convex,
618  // because the contributions in the opposite direction will
619  // partly cancel out the ones in the other direction,
620  // in just the right amount.
621  double poly_sum = 0;
622  const UT_Vector3D p0 = poly.getPos3(0);
623  UT_Vector3D prev = poly.getPos3(1);
624  for (GA_Size i = 2; i < n; ++i)
625  {
626  const UT_Vector3D next = poly.getPos3(i);
627  poly_sum += UTsignedSolidAngleTri(p0, prev, next, query_point);
628  prev = next;
629  }
630  soup_sum += poly_sum;
631  }
632  local_sum += soup_sum;
633  }
634  else if (primtype == GEO_PRIMMESH)
635  {
636  double mesh_sum = 0;
637  const GEO_PrimMesh *mesh = UTverify_cast<const GEO_PrimMesh *>(mesh_geo->getPrimitive(primoff));
638  const int nquadrows = mesh->getNumRows() - !mesh->isWrappedV();
639  const int nquadcols = mesh->getNumCols() - !mesh->isWrappedU();
640  for (int row = 0; row < nquadrows; ++row)
641  {
642  for (int col = 0; col < nquadcols; ++col)
643  {
644  GEO_Hull::Poly poly(*mesh, row, col);
645  const UT_Vector3D p0 = poly.getPos3(0);
646  const UT_Vector3D p1 = poly.getPos3(1);
647  const UT_Vector3D p2 = poly.getPos3(2);
648  const UT_Vector3D p3 = poly.getPos3(3);
649 #if QUADS_AS_BILINEAR_PATCHES
650  const double poly_sum = UTsignedSolidAngleQuad(p0, p1, p2, p3, query_point);
651 #else
652  const double poly_sum =
653  UTsignedSolidAngleTri(p0, p1, p2, query_point) +
654  UTsignedSolidAngleTri(p0, p2, p3, query_point);
655 #endif
656  mesh_sum += poly_sum;
657  }
658  }
659  local_sum += mesh_sum;
660  }
661  else if (primtype == GEO_PRIMNURBSURF || primtype == GEO_PRIMBEZSURF)
662  {
663  const GEO_TPSurf *surf = UTverify_cast<const GEO_TPSurf *>(mesh_geo->getPrimitive(primoff));
664  const int nquadrows = surf->getNumRows() - !surf->isWrappedV();
665  const int nquadcols = surf->getNumCols() - !surf->isWrappedU();
666  if (nquadcols <= 0 || nquadrows <= 0)
667  continue;
668 
669  // For slightly better accuracy, we use the greville points
670  // instead of the hull points.
671  double surf_sum = 0;
672  UT_FloatArray ugrevilles;
673  UT_FloatArray vgrevilles;
674  surf->getUGrevilles(ugrevilles);
675  if (surf->isWrappedU())
676  ugrevilles.append(0);
677  surf->getVGrevilles(vgrevilles);
678  if (surf->isWrappedV())
679  vgrevilles.append(0);
680  UT_ASSERT(ugrevilles.size() >= nquadcols+1);
681  UT_ASSERT(vgrevilles.size() >= nquadrows+1);
682  for (int row = 0; row < nquadrows; ++row)
683  {
684  for (int col = 0; col < nquadcols; ++col)
685  {
686  UT_Vector4 a4;
687  UT_Vector4 b4;
688  UT_Vector4 c4;
689  UT_Vector4 d4;
690  surf->evaluateInteriorPoint(a4, ugrevilles(col), vgrevilles(row));
691  surf->evaluateInteriorPoint(b4, ugrevilles(col+1), vgrevilles(row));
692  surf->evaluateInteriorPoint(c4, ugrevilles(col+1), vgrevilles(row+1));
693  surf->evaluateInteriorPoint(d4, ugrevilles(col), vgrevilles(row+1));
694  const UT_Vector3D p0(a4);
695  const UT_Vector3D p1(b4);
696  const UT_Vector3D p2(c4);
697  const UT_Vector3D p3(d4);
698 #if QUADS_AS_BILINEAR_PATCHES
699  const double poly_sum = UTsignedSolidAngleQuad(p0, p1, p2, p3, query_point);
700 #else
701  const double poly_sum =
702  UTsignedSolidAngleTri(p0, p1, p2, query_point) +
703  UTsignedSolidAngleTri(p0, p2, p3, query_point);
704 #endif
705  surf_sum += poly_sum;
706  }
707  }
708  local_sum += surf_sum;
709  }
710  else if (primtype == GEO_PRIMSPHERE)
711  {
712  const GEO_PrimSphere *sphere = UTverify_cast<const GEO_PrimSphere *>(mesh_geo->getPrimitive(primoff));
713  UT_Vector3D query_minus_centre = query_point - UT_Vector3D(sphere->getPos3());
715  sphere->getLocalTransform(transform);
716 
717  // If we're outside, no contribution.
718  // If we're inside and the determinant is negative, +4pi
719  // If we're inside and the determinant is positive, -4pi
720  // The negation is for consistency with the winding order of polygons
721  // that would be created if the sphere were converted to polygons.
722  double determinant = transform.determinant();
723  if (determinant == 0)
724  continue;
725 
726  int failed = transform.invert();
727  if (failed)
728  continue;
729 
730  query_minus_centre.rowVecMult(transform);
731  float length2 = query_minus_centre.length2();
732 
733  double sphere_sum = 0;
734  if (length2 < 1)
735  sphere_sum = 4*M_PI;
736  else if (length2 == 1)
737  sphere_sum = 2*M_PI;
738 
739  if (determinant > 0)
740  sphere_sum = -sphere_sum;
741 
742  local_sum += sphere_sum;
743  }
744  }
745  *sum = local_sum;
746 }
747 
748 static void
749 sopSumContributions2D(
750  double *sum,
751  const UT_Vector2D &query_point,
752  const GEO_Detail *const mesh_geo,
753  const GA_OffsetList &primoffs,
754  const exint start, const exint end,
755  const int axis0, const int axis1)
756 {
757  if (end-start > 1024)
758  {
759  exint mid = (start+end)>>1;
760  double sum0;
761  double sum1;
762  UTparallelInvoke(true, [&] {
763  sopSumContributions2D(&sum0, query_point, mesh_geo, primoffs, start, mid, axis0, axis1);
764  }, [&] {
765  sopSumContributions2D(&sum1, query_point, mesh_geo, primoffs, mid, end, axis0, axis1);
766  });
767  *sum = sum0 + sum1;
768  return;
769  }
770 
771  double local_sum = 0;
772  for (exint arrayi = start; arrayi < end; ++arrayi)
773  {
774  GA_Offset primoff = primoffs(arrayi);
775  int primtype = mesh_geo->getPrimitiveTypeId(primoff);
776  if (primtype == GA_PRIMPOLY)
777  {
778  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
779  const GA_Size n = vertices.size();
780  const bool closed = vertices.getExtraFlag();
781  if (n < 2+int(closed))
782  continue;
783 
784  double poly_sum = 0;
785  UT_Vector3 pos = mesh_geo->getPos3(mesh_geo->vertexPoint(vertices(0)));
786  UT_Vector2D prev(pos[axis0], pos[axis1]);
787  const UT_Vector2D p0 = prev;
788  for (GA_Size i = 1; i < n; ++i)
789  {
790  pos = mesh_geo->getPos3(mesh_geo->vertexPoint(vertices(i)));
791  const UT_Vector2D next(pos[axis0], pos[axis1]);
792  poly_sum += UTsignedAngleSegment(prev, next, query_point);
793  prev = next;
794  }
795  if (closed)
796  poly_sum += UTsignedAngleSegment(prev, p0, query_point);
797 
798  local_sum += poly_sum;
799  }
800  else if (primtype == GA_PRIMTETRAHEDRON)
801  {
802  // NOTE: Tetrahedra never contribute to the 2D winding number, since
803  // every point is contained in 1 forward triangle and 1 backward triangle.
804  }
805  else if (primtype == GA_PRIMPOLYSOUP)
806  {
807  double soup_sum = 0;
808  const GEO_PrimPolySoup *soup = UTverify_cast<const GEO_PrimPolySoup *>(mesh_geo->getPrimitive(primoff));
809  for (GEO_PrimPolySoup::PolygonIterator poly(*soup); !poly.atEnd(); ++poly)
810  {
811  GA_Size n = poly.nvertices();
812  if (n < 3)
813  continue;
814 
815  double poly_sum = 0;
816  UT_Vector3 pos = poly.getPos3(0);
817  UT_Vector2D prev(pos[axis0], pos[axis1]);
818  const UT_Vector2D p0 = prev;
819  for (GA_Size i = 1; i < n; ++i)
820  {
821  pos = poly.getPos3(i);
822  const UT_Vector2D next(pos[axis0], pos[axis1]);
823  poly_sum += UTsignedAngleSegment(prev, next, query_point);
824  prev = next;
825  }
826  poly_sum += UTsignedAngleSegment(prev, p0, query_point);
827  soup_sum += poly_sum;
828  }
829  local_sum += soup_sum;
830  }
831  else if (primtype == GEO_PRIMMESH)
832  {
833  double mesh_sum = 0;
834  const GEO_PrimMesh *mesh = UTverify_cast<const GEO_PrimMesh *>(mesh_geo->getPrimitive(primoff));
835  const int nquadrows = mesh->getNumRows() - !mesh->isWrappedV();
836  const int nquadcols = mesh->getNumCols() - !mesh->isWrappedU();
837  for (int row = 0; row < nquadrows; ++row)
838  {
839  for (int col = 0; col < nquadcols; ++col)
840  {
841  GEO_Hull::Poly poly(*mesh, row, col);
842  double poly_sum = 0;
843  UT_Vector3 pos = poly.getPos3(0);
844  UT_Vector2D prev(pos[axis0], pos[axis1]);
845  const UT_Vector2D p0 = prev;
846  for (GA_Size i = 1; i < 4; ++i)
847  {
848  pos = poly.getPos3(i);
849  const UT_Vector2D next(pos[axis0], pos[axis1]);
850  poly_sum += UTsignedAngleSegment(prev, next, query_point);
851  prev = next;
852  }
853  poly_sum += UTsignedAngleSegment(prev, p0, query_point);
854  mesh_sum += poly_sum;
855  }
856  }
857  local_sum += mesh_sum;
858  }
859  else if (primtype == GEO_PRIMNURBSURF || primtype == GEO_PRIMBEZSURF)
860  {
861  const GEO_TPSurf *surf = UTverify_cast<const GEO_TPSurf *>(mesh_geo->getPrimitive(primoff));
862  const int nquadrows = surf->getNumRows() - !surf->isWrappedV();
863  const int nquadcols = surf->getNumCols() - !surf->isWrappedU();
864  if (nquadcols <= 0 || nquadrows <= 0)
865  continue;
866 
867  // For slightly better accuracy, we use the greville points
868  // instead of the hull points.
869  double surf_sum = 0;
870  UT_FloatArray ugrevilles;
871  UT_FloatArray vgrevilles;
872  surf->getUGrevilles(ugrevilles);
873  if (surf->isWrappedU())
874  ugrevilles.append(0);
875  surf->getVGrevilles(vgrevilles);
876  if (surf->isWrappedV())
877  vgrevilles.append(0);
878  UT_ASSERT(ugrevilles.size() >= nquadcols+1);
879  UT_ASSERT(vgrevilles.size() >= nquadrows+1);
880  for (int row = 0; row < nquadrows; ++row)
881  {
882  for (int col = 0; col < nquadcols; ++col)
883  {
884  UT_Vector4 a4;
885  UT_Vector4 b4;
886  UT_Vector4 c4;
887  UT_Vector4 d4;
888  surf->evaluateInteriorPoint(a4, ugrevilles(col), vgrevilles(row));
889  surf->evaluateInteriorPoint(b4, ugrevilles(col+1), vgrevilles(row));
890  surf->evaluateInteriorPoint(c4, ugrevilles(col+1), vgrevilles(row+1));
891  surf->evaluateInteriorPoint(d4, ugrevilles(col), vgrevilles(row+1));
892  const UT_Vector2D p0(a4[axis0],a4[axis1]);
893  const UT_Vector2D p1(b4[axis0],b4[axis1]);
894  const UT_Vector2D p2(c4[axis0],c4[axis1]);
895  const UT_Vector2D p3(d4[axis0],d4[axis1]);
896  double poly_sum;
897  poly_sum = UTsignedAngleSegment(p0, p1, query_point);
898  poly_sum += UTsignedAngleSegment(p1, p2, query_point);
899  poly_sum += UTsignedAngleSegment(p2, p3, query_point);
900  poly_sum += UTsignedAngleSegment(p3, p0, query_point);
901  surf_sum += poly_sum;
902  }
903  }
904  local_sum += surf_sum;
905  }
906  else if (primtype == GEO_PRIMNURBCURVE || primtype == GEO_PRIMBEZCURVE)
907  {
908  const GEO_Curve *curve = UTverify_cast<const GEO_Curve *>(mesh_geo->getPrimitive(primoff));
909  const int nedges = curve->getVertexCount() - !curve->isClosed();
910  if (nedges <= 0)
911  continue;
912 
913  // For slightly better accuracy, we use the greville points
914  // instead of the hull points.
915  double curve_sum = 0;
916  UT_Array<GEO_Greville> grevilles;
917  grevilles.setCapacity(nedges+1);
918  curve->buildGrevilles(grevilles);
919  if (curve->isClosed())
920  grevilles.append(grevilles[0]);
921  UT_ASSERT(grevilles.size() >= nedges+1);
922  UT_Vector2D prev(grevilles[0].getCVImage()[axis0], grevilles[0].getCVImage()[axis1]);
923  for (int i = 0; i < nedges; ++i)
924  {
925  const UT_Vector2D next(grevilles[i].getCVImage()[axis0], grevilles[i].getCVImage()[axis1]);
926  curve_sum += UTsignedAngleSegment(prev, next, query_point);
927  }
928  local_sum += curve_sum;
929  }
930  }
931  *sum = local_sum;
932 }
933 
934 static void
935 sopAccumulateTriangles(
936  const GA_Detail *const mesh_geo,
937  const GA_Offset primoff,
938  const UT_Array<int> &ptmap,
939  UT_Array<int> &triangle_points)
940 {
941  const bool has_ptmap = !ptmap.isEmpty();
942  int primtype = mesh_geo->getPrimitiveTypeId(primoff);
943  if (primtype == GA_PRIMPOLY)
944  {
945  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
946  const GA_Size n = vertices.size();
947  const bool closed = vertices.getExtraFlag();
948  if (n < 3 || !closed)
949  return;
950 
951  // A triangle fan suffices, even if the polygon is non-convex,
952  // because the contributions in the opposite direction will
953  // partly cancel out the ones in the other direction,
954  // in just the right amount.
955  const GA_Offset p0 = mesh_geo->vertexPoint(vertices(0));
956  const int p0i = has_ptmap ? ptmap[p0] : int(mesh_geo->pointIndex(p0));
957  GA_Offset prev = mesh_geo->vertexPoint(vertices(1));
958  int previ = has_ptmap ? ptmap[prev] : int(mesh_geo->pointIndex(prev));
959  for (GA_Size i = 2; i < n; ++i)
960  {
961  const GA_Offset next = mesh_geo->vertexPoint(vertices(i));
962  const int nexti = has_ptmap ? ptmap[next] : int(mesh_geo->pointIndex(next));
963  triangle_points.append(p0i);
964  triangle_points.append(previ);
965  triangle_points.append(nexti);
966  previ = nexti;
967  }
968  }
969  else if (primtype == GA_PRIMTETRAHEDRON)
970  {
971  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
972  const GEO_PrimTetrahedron tet(SYSconst_cast(mesh_geo), primoff, vertices);
973 
974  for (int i = 0; i < 4; ++i)
975  {
976  // Ignore shared tet faces. They would contribute exactly opposite amounts.
977  if (tet.isFaceShared(i))
978  continue;
979 
980  const int *face_indices = GEO_PrimTetrahedron::fastFaceIndices(i);
981  const GA_Offset a = mesh_geo->vertexPoint(vertices(face_indices[0]));
982  const GA_Offset b = mesh_geo->vertexPoint(vertices(face_indices[1]));
983  const GA_Offset c = mesh_geo->vertexPoint(vertices(face_indices[2]));
984  const int ai = has_ptmap ? ptmap[a] : int(mesh_geo->pointIndex(a));
985  const int bi = has_ptmap ? ptmap[b] : int(mesh_geo->pointIndex(b));
986  const int ci = has_ptmap ? ptmap[c] : int(mesh_geo->pointIndex(c));
987  triangle_points.append(ai);
988  triangle_points.append(bi);
989  triangle_points.append(ci);
990  }
991  }
992  else if (primtype == GA_PRIMPOLYSOUP)
993  {
994  const GEO_PrimPolySoup *soup = UTverify_cast<const GEO_PrimPolySoup *>(mesh_geo->getPrimitive(primoff));
995  for (GEO_PrimPolySoup::PolygonIterator poly(*soup); !poly.atEnd(); ++poly)
996  {
997  GA_Size n = poly.nvertices();
998  if (n < 3)
999  continue;
1000 
1001  // A triangle fan suffices, even if the polygon is non-convex,
1002  // because the contributions in the opposite direction will
1003  // partly cancel out the ones in the other direction,
1004  // in just the right amount.
1005  const GA_Offset p0 = poly.getPointOffset(0);
1006  const int p0i = has_ptmap ? ptmap[p0] : int(mesh_geo->pointIndex(p0));
1007  GA_Offset prev = poly.getPointOffset(1);
1008  int previ = has_ptmap ? ptmap[prev] : int(mesh_geo->pointIndex(prev));
1009  for (GA_Size i = 2; i < n; ++i)
1010  {
1011  const GA_Offset next = poly.getPointOffset(i);
1012  const int nexti = has_ptmap ? ptmap[next] : int(mesh_geo->pointIndex(next));
1013  triangle_points.append(p0i);
1014  triangle_points.append(previ);
1015  triangle_points.append(nexti);
1016  previ = nexti;
1017  }
1018  }
1019  }
1020  else if (primtype == GEO_PRIMMESH || primtype == GEO_PRIMNURBSURF || primtype == GEO_PRIMBEZSURF)
1021  {
1022  // In this mode, we're only using points in the detail, so no grevilles.
1023  const GEO_Hull *mesh = UTverify_cast<const GEO_Hull *>(mesh_geo->getPrimitive(primoff));
1024  const int nquadrows = mesh->getNumRows() - !mesh->isWrappedV();
1025  const int nquadcols = mesh->getNumCols() - !mesh->isWrappedU();
1026  for (int row = 0; row < nquadrows; ++row)
1027  {
1028  for (int col = 0; col < nquadcols; ++col)
1029  {
1030  GEO_Hull::Poly poly(*mesh, row, col);
1031  const GA_Offset a = poly.getPointOffset(0);
1032  const GA_Offset b = poly.getPointOffset(1);
1033  const GA_Offset c = poly.getPointOffset(2);
1034  const GA_Offset d = poly.getPointOffset(3);
1035  const int ai = has_ptmap ? ptmap[a] : int(mesh_geo->pointIndex(a));
1036  const int bi = has_ptmap ? ptmap[b] : int(mesh_geo->pointIndex(b));
1037  const int ci = has_ptmap ? ptmap[c] : int(mesh_geo->pointIndex(c));
1038  const int di = has_ptmap ? ptmap[d] : int(mesh_geo->pointIndex(d));
1039  triangle_points.append(ai);
1040  triangle_points.append(bi);
1041  triangle_points.append(ci);
1042  triangle_points.append(ai);
1043  triangle_points.append(ci);
1044  triangle_points.append(di);
1045  }
1046  }
1047  }
1048 }
1049 
1050 static void
1051 sopAccumulateSegments(
1052  const GA_Detail *const mesh_geo,
1053  const GA_Offset primoff,
1054  const UT_Array<int> &ptmap,
1055  UT_Array<int> &segment_points)
1056 {
1057  const bool has_ptmap = !ptmap.isEmpty();
1058  int primtype = mesh_geo->getPrimitiveTypeId(primoff);
1059  if (primtype == GA_PRIMPOLY || primtype == GA_PRIMBEZCURVE || primtype == GA_PRIMNURBCURVE)
1060  {
1061  const GA_OffsetListRef vertices = mesh_geo->getPrimitiveVertexList(primoff);
1062  const GA_Size n = vertices.size();
1063  const bool closed = vertices.getExtraFlag();
1064  if (n < 2+int(closed))
1065  return;
1066 
1067  GA_Offset prev = mesh_geo->vertexPoint(vertices(0));
1068  int previ = has_ptmap ? ptmap[prev] : int(mesh_geo->pointIndex(prev));
1069  const int pt0i = previ;
1070  for (GA_Size i = 1; i < n; ++i)
1071  {
1072  const GA_Offset next = mesh_geo->vertexPoint(vertices(i));
1073  const int nexti = has_ptmap ? ptmap[next] : int(mesh_geo->pointIndex(next));
1074  segment_points.append(previ);
1075  segment_points.append(nexti);
1076  previ = nexti;
1077  }
1078  if (closed)
1079  {
1080  segment_points.append(previ);
1081  segment_points.append(pt0i);
1082  }
1083  }
1084  else if (primtype == GA_PRIMTETRAHEDRON)
1085  {
1086  // NOTE: Tetrahedra never contribute to the 2D winding number, since
1087  // every point is contained in 1 forward triangle and 1 backward triangle.
1088  }
1089  else if (primtype == GA_PRIMPOLYSOUP)
1090  {
1091  const GEO_PrimPolySoup *soup = UTverify_cast<const GEO_PrimPolySoup *>(mesh_geo->getPrimitive(primoff));
1092  for (GEO_PrimPolySoup::PolygonIterator poly(*soup); !poly.atEnd(); ++poly)
1093  {
1094  GA_Size n = poly.nvertices();
1095  if (n < 3)
1096  continue;
1097 
1098  GA_Offset prev = poly.getPointOffset(0);
1099  int previ = has_ptmap ? ptmap[prev] : int(mesh_geo->pointIndex(prev));
1100  const int pt0i = previ;
1101  for (GA_Size i = 1; i < n; ++i)
1102  {
1103  const GA_Offset next = poly.getPointOffset(i);
1104  const int nexti = has_ptmap ? ptmap[next] : int(mesh_geo->pointIndex(next));
1105  segment_points.append(previ);
1106  segment_points.append(nexti);
1107  previ = nexti;
1108  }
1109  segment_points.append(previ);
1110  segment_points.append(pt0i);
1111  }
1112  }
1113  else if (primtype == GEO_PRIMMESH || primtype == GEO_PRIMNURBSURF || primtype == GEO_PRIMBEZSURF)
1114  {
1115  // In this mode, we're only using points in the detail, so no grevilles.
1116  const GEO_Hull *mesh = UTverify_cast<const GEO_Hull *>(mesh_geo->getPrimitive(primoff));
1117  const int nquadrows = mesh->getNumRows() - !mesh->isWrappedV();
1118  const int nquadcols = mesh->getNumCols() - !mesh->isWrappedU();
1119  for (int row = 0; row < nquadrows; ++row)
1120  {
1121  for (int col = 0; col < nquadcols; ++col)
1122  {
1123  GEO_Hull::Poly poly(*mesh, row, col);
1124  const GA_Offset a = poly.getPointOffset(0);
1125  const GA_Offset b = poly.getPointOffset(1);
1126  const GA_Offset c = poly.getPointOffset(2);
1127  const GA_Offset d = poly.getPointOffset(3);
1128  const int ai = has_ptmap ? ptmap[a] : int(mesh_geo->pointIndex(a));
1129  const int bi = has_ptmap ? ptmap[b] : int(mesh_geo->pointIndex(b));
1130  const int ci = has_ptmap ? ptmap[c] : int(mesh_geo->pointIndex(c));
1131  const int di = has_ptmap ? ptmap[d] : int(mesh_geo->pointIndex(d));
1132  segment_points.append(ai);
1133  segment_points.append(bi);
1134  segment_points.append(bi);
1135  segment_points.append(ci);
1136  segment_points.append(ci);
1137  segment_points.append(di);
1138  segment_points.append(di);
1139  segment_points.append(ai);
1140  }
1141  }
1142  }
1143 }
1144 
1145 static void
1146 sop3DFullAccuracy(
1147  const GEO_Detail *const query_points,
1148  const GA_SplittableRange &point_range,
1149  const GEO_Detail *const mesh_geo,
1150  const GA_PrimitiveGroup *const mesh_prim_group,
1151  const GA_RWHandleF &winding_number_attrib,
1152  const bool as_solid_angle,
1153  const bool negate)
1154 {
1155  UT_AutoInterrupt boss("Computing Winding Numbers");
1156 
1157  // See comment below for why we create this GA_Offsetlist of mesh primitive offsets.
1158  GA_OffsetList primoffs;
1159  if (!mesh_prim_group && mesh_geo->getPrimitiveMap().isTrivialMap())
1160  {
1161  primoffs.setTrivial(GA_Offset(0), mesh_geo->getNumPrimitives());
1162  }
1163  else
1164  {
1165  GA_Offset start;
1166  GA_Offset end;
1167  for (GA_Iterator it(mesh_geo->getPrimitiveRange(mesh_prim_group)); it.fullBlockAdvance(start, end); )
1168  {
1169  primoffs.setTrivialRange(primoffs.size(), start, end-start);
1170  }
1171  }
1172 
1173  UTparallelFor(point_range, [query_points,mesh_geo,winding_number_attrib,&primoffs,as_solid_angle,negate,&boss](const GA_SplittableRange &r)
1174  {
1175  GA_Offset start;
1176  GA_Offset end;
1177  for (GA_Iterator it(r); it.blockAdvance(start, end); )
1178  {
1179  for (GA_Offset ptoff = start; ptoff != end; ++ptoff)
1180  {
1181  if (boss.wasInterrupted())
1182  return;
1183 
1184  const UT_Vector3 query_point = query_points->getPos3(ptoff);
1185 
1186  // NOTE: We can't use UTparallelReduce, because that would have
1187  // nondeterministic roundoff error due to floating-point
1188  // addition being non-associative. We can't just use
1189  // UTparallelInvoke with GA_SplittableRange either, because
1190  // the roundoff error would change with defragmentation,
1191  // e.g. from locking the mesh input and reloading the HIP file.
1192  // Instead, we always split in the middle of the list of offsets.
1193  double sum;
1194  sopSumContributions3D(&sum, query_point, mesh_geo, primoffs, 0, primoffs.size());
1195 
1196  if (!as_solid_angle)
1197  sum *= (0.25*M_1_PI); // Divide by 4pi (solid angle of full sphere)
1198  if (negate)
1199  sum = -sum;
1200 
1201  winding_number_attrib.set(ptoff, sum);
1202  }
1203  }
1204  });
1205 }
1206 
1207 static void
1208 sop2DFullAccuracy(
1209  const GEO_Detail *const query_points,
1210  const GA_SplittableRange &point_range,
1211  const GEO_Detail *const mesh_geo,
1212  const GA_PrimitiveGroup *const mesh_prim_group,
1213  const GA_RWHandleF &winding_number_attrib,
1214  const bool as_solid_angle,
1215  const bool negate,
1216  const int axis0,
1217  const int axis1)
1218 {
1219  UT_AutoInterrupt boss("Computing Winding Numbers");
1220 
1221  // See comment below for why we create this GA_Offsetlist of mesh primitive offsets.
1222  GA_OffsetList primoffs;
1223  if (!mesh_prim_group && mesh_geo->getPrimitiveMap().isTrivialMap())
1224  {
1225  primoffs.setTrivial(GA_Offset(0), mesh_geo->getNumPrimitives());
1226  }
1227  else
1228  {
1229  GA_Offset start;
1230  GA_Offset end;
1231  for (GA_Iterator it(mesh_geo->getPrimitiveRange(mesh_prim_group)); it.fullBlockAdvance(start, end); )
1232  {
1233  primoffs.setTrivialRange(primoffs.size(), start, end-start);
1234  }
1235  }
1236 
1237  UTparallelFor(point_range, [query_points,mesh_geo,winding_number_attrib,&primoffs,as_solid_angle,negate,&boss,axis0,axis1](const GA_SplittableRange &r)
1238  {
1239  GA_Offset start;
1240  GA_Offset end;
1241  for (GA_Iterator it(r); it.blockAdvance(start, end); )
1242  {
1243  for (GA_Offset ptoff = start; ptoff != end; ++ptoff)
1244  {
1245  if (boss.wasInterrupted())
1246  return;
1247 
1248  const UT_Vector3D query_point_3d = query_points->getPos3(ptoff);
1249  const UT_Vector2D query_point(query_point_3d[axis0], query_point_3d[axis1]);
1250 
1251  // NOTE: We can't use UTparallelReduce, because that would have
1252  // nondeterministic roundoff error due to floating-point
1253  // addition being non-associative. We can't just use
1254  // UTparallelInvoke with GA_SplittableRange either, because
1255  // the roundoff error would change with defragmentation,
1256  // e.g. from locking the mesh input and reloading the HIP file.
1257  // Instead, we always split in the middle of the list of offsets.
1258  double sum;
1259  sopSumContributions2D(&sum, query_point, mesh_geo, primoffs, 0, primoffs.size(), axis0, axis1);
1260 
1261  if (!as_solid_angle)
1262  sum *= (0.5*M_1_PI); // Divide by 2pi (angle of full circle)
1263  if (negate)
1264  sum = -sum;
1265 
1266  winding_number_attrib.set(ptoff, sum);
1267  }
1268  }
1269  });
1270 }
1271 
1272 static void
1273 sop3DApproximate(
1274  const GEO_Detail *const query_points,
1275  const GA_SplittableRange &point_range,
1276  const UT_SolidAngle<float,float> &solid_angle_tree,
1277  const double accuracy_scale,
1278  const GA_RWHandleF &winding_number_attrib,
1279  const bool as_solid_angle,
1280  const bool negate)
1281 {
1282  UT_AutoInterrupt boss("Computing Winding Numbers");
1283 
1284  UTparallelFor(point_range, [query_points,&winding_number_attrib,&solid_angle_tree,as_solid_angle,negate,accuracy_scale,&boss](const GA_SplittableRange &r)
1285  {
1286  GA_Offset start;
1287  GA_Offset end;
1288  for (GA_Iterator it(r); it.blockAdvance(start, end); )
1289  {
1290  if (boss.wasInterrupted())
1291  return;
1292 
1293  for (GA_Offset ptoff = start; ptoff != end; ++ptoff)
1294  {
1295  const UT_Vector3 query_point = query_points->getPos3(ptoff);
1296 
1297  double sum = solid_angle_tree.computeSolidAngle(query_point, accuracy_scale);
1298 
1299  if (!as_solid_angle)
1300  sum *= (0.25*M_1_PI); // Divide by 4pi (solid angle of full sphere)
1301  if (negate)
1302  sum = -sum;
1303 
1304  winding_number_attrib.set(ptoff, sum);
1305  }
1306  }
1307  }, 10); // Large subscribe ratio, because expensive points are often clustered
1308 }
1309 
1310 static void
1311 sop2DApproximate(
1312  const GEO_Detail *const query_points,
1313  const GA_SplittableRange &point_range,
1314  const UT_SubtendedAngle<float,float> &subtended_angle_tree,
1315  const double accuracy_scale,
1316  const GA_RWHandleF &winding_number_attrib,
1317  const bool as_angle,
1318  const bool negate,
1319  const int axis0,
1320  const int axis1)
1321 {
1322  UT_AutoInterrupt boss("Computing Winding Numbers");
1323 
1324  UTparallelFor(point_range, [query_points,&winding_number_attrib,&subtended_angle_tree,as_angle,negate,accuracy_scale,&boss,axis0,axis1](const GA_SplittableRange &r)
1325  {
1326  GA_Offset start;
1327  GA_Offset end;
1328  for (GA_Iterator it(r); it.blockAdvance(start, end); )
1329  {
1330  if (boss.wasInterrupted())
1331  return;
1332 
1333  for (GA_Offset ptoff = start; ptoff != end; ++ptoff)
1334  {
1335  const UT_Vector3 query_point_3d = query_points->getPos3(ptoff);
1336  const UT_Vector2D query_point(query_point_3d[axis0], query_point_3d[axis1]);
1337 
1338  double sum = subtended_angle_tree.computeAngle(query_point, accuracy_scale);
1339 
1340  if (!as_angle)
1341  sum *= (0.5*M_1_PI); // Divide by 2pi (angle of full circle)
1342  if (negate)
1343  sum = -sum;
1344 
1345  winding_number_attrib.set(ptoff, sum);
1346  }
1347  }
1348  }, 10); // Large subscribe ratio, because expensive points are often clustered
1349 }
1350 
1351 /// This is the function that does the actual work.
1352 void SOP_WindingNumberVerb::cook(const CookParms &cookparms) const
1353 {
1354  // This gives easy access to all of the current parameter values
1355  auto &&sopparms = cookparms.parms<SOP_WindingNumberParms>();
1356  auto sopcache = (SOP_WindingNumberCache *)cookparms.cache();
1357 
1358  // The output detail
1359  GEO_Detail *const query_points = cookparms.gdh().gdpNC();
1360 
1361  const GEO_Detail *const mesh_geo = cookparms.inputGeo(1);
1362 
1363  GA_RWHandleF winding_number_attrib(query_points->findFloatTuple(GA_ATTRIB_POINT, sopparms.getAttrib()));
1364  if (!winding_number_attrib.isValid())
1365  winding_number_attrib.bind(query_points->addFloatTuple(GA_ATTRIB_POINT, sopparms.getAttrib(), 1));
1366 
1367  if (!winding_number_attrib.isValid())
1368  {
1369  cookparms.sopAddError(SOP_ATTRIBUTE_INVALID, sopparms.getAttrib());
1370  return;
1371  }
1372 
1373  GA_Size npoints = query_points->getNumPoints();
1374  if (npoints == 0)
1375  return;
1376 
1377  GOP_Manager group_manager;
1378 
1379  // Parse group of mesh primitives
1380  const UT_StringHolder &mesh_prim_group_string = sopparms.getMeshPrims();
1381  const GA_PrimitiveGroup *mesh_prim_group = nullptr;
1382  if (mesh_prim_group_string.isstring())
1383  {
1384  bool success;
1385  mesh_prim_group = group_manager.parsePrimitiveDetached(
1386  mesh_prim_group_string.c_str(),
1387  mesh_geo, true, success);
1388  }
1389 
1390  // Parse group of query points
1391  const UT_StringHolder &query_point_group_string = sopparms.getQueryPoints();
1392  const GA_PointGroup *query_point_group = nullptr;
1393  if (query_point_group_string.isstring())
1394  {
1395  bool success;
1396  query_point_group = group_manager.parsePointDetached(
1397  query_point_group_string.c_str(),
1398  query_points, true, success);
1399  }
1400 
1401  const GA_SplittableRange point_range(query_points->getPointRange(query_point_group));
1402 
1403  const bool full_accuracy = sopparms.getFullAccuracy();
1404  const bool as_solid_angle = sopparms.getAsSolidAngle();
1405 
1406  Type winding_number_type = sopparms.getType();
1407 
1408  if (winding_number_type == Type::XYZ)
1409  {
1410  // NOTE: The negation is because Houdini's normals are
1411  // left-handed with respect to the winding order,
1412  // whereas most use the right-handed convension.
1413  const bool negate = !sopparms.getNegate();
1414  if (full_accuracy)
1415  {
1416  sopcache->clear();
1417 
1418  sop3DFullAccuracy(
1419  query_points, point_range,
1420  mesh_geo, mesh_prim_group,
1421  winding_number_attrib,
1422  as_solid_angle, negate);
1423  }
1424  else
1425  {
1426  sopcache->update3D(*mesh_geo, mesh_prim_group, mesh_prim_group_string, 2);
1427 
1428  const UT_SolidAngle<float,float> &solid_angle_tree = sopcache->mySolidAngleTree;
1429  const double accuracy_scale = sopparms.getAccuracyScale();
1430 
1431  sop3DApproximate(
1432  query_points, point_range,
1433  solid_angle_tree, accuracy_scale,
1434  winding_number_attrib,
1435  as_solid_angle, negate);
1436  }
1437  }
1438  else
1439  {
1440  int axis0 = int(winding_number_type)-1;
1441  int axis1 = (axis0 == 2) ? 0 : (axis0+1);
1442 
1443  const bool negate = sopparms.getNegate();
1444  if (full_accuracy)
1445  {
1446  sopcache->clear();
1447 
1448  sop2DFullAccuracy(
1449  query_points, point_range,
1450  mesh_geo, mesh_prim_group,
1451  winding_number_attrib,
1452  as_solid_angle, negate,
1453  axis0, axis1);
1454  }
1455  else
1456  {
1457  sopcache->update2D(*mesh_geo, mesh_prim_group, mesh_prim_group_string, 2, axis0, axis1);
1458 
1459  const UT_SubtendedAngle<float,float> &subtended_angle_tree = sopcache->mySubtendedAngleTree;
1460  const double accuracy_scale = sopparms.getAccuracyScale();
1461 
1462  sop2DApproximate(
1463  query_points, point_range,
1464  subtended_angle_tree, accuracy_scale,
1465  winding_number_attrib,
1466  as_solid_angle, negate,
1467  axis0, axis1);
1468  }
1469  }
1470 
1471  winding_number_attrib->bumpDataId();
1472 }
bool isWrappedU() const
Definition: GEO_Hull.h:488
static PRM_ChoiceList primGroupMenu
Definition: SOP_Node.h:1189
SYS_FORCE_INLINE void bumpDataId()
Definition: GA_Attribute.h:306
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
virtual void cook(const CookParms &cookparms) const
This is the function that does the actual work.
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)
UT_Array< UT_Vector3 > myPositions3D
SYS_FORCE_INLINE GA_Offset getPointOffset(GA_Size i) const
Definition: GA_Primitive.h:254
SYS_FORCE_INLINE GA_Primitive * getPrimitive(GA_Offset prim_off)
Definition: GA_Detail.h:429
fpreal64 getAccuracyScale() const
bool isWrappedV() const
Definition: GEO_Hull.h:489
virtual UT_StringHolder name() const
void getVGrevilles(UT_Array< float > &vgrevilles) const
UT_StringHolder myGroupString
const GA_IndexMap & getPrimitiveMap() const
Definition: GA_Detail.h:743
SYS_FORCE_INLINE GA_Size getVertexCount() const
Return the number of vertices used by this primitive.
Definition: GA_Primitive.h:232
virtual CookMode cookMode(const SOP_NodeParms *parms) const
GA_DataId getDataId() const
Iteration over a range of elements.
Definition: GA_Iterator.h:29
void setChoiceListPtr(const UT_StringRef &name, PRM_ChoiceList *list)
SYS_FORCE_INLINE int getPrimitiveTypeId(GA_Offset primoff) const
Definition: GA_Primitive.h:907
UT_Vector2T< float > UT_Vector2
GA_Attribute * addFloatTuple(GA_AttributeOwner owner, GA_AttributeScope scope, const UT_StringHolder &name, int tuple_size, const GA_Defaults &defaults=GA_Defaults(0.0), const UT_Options *creation_args=0, const GA_AttributeOptions *attribute_options=0, GA_Storage storage=GA_STORE_REAL32, const GA_ReuseStrategy &reuse=GA_ReuseStrategy())
void update3D(const GA_Detail &mesh_geo, const GA_PrimitiveGroup *prim_group, const UT_StringHolder &group_string, const int approx_order)
int64 GA_DataId
Definition: GA_Types.h:687
#define GA_INVALID_DATAID
#define M_PI
Definition: fmath.h:90
bool blockAdvance(GA_Offset &start, GA_Offset &end)
GLuint start
Definition: glcorearb.h:475
GA_Size entries() const overridefinal
Will return the number of primary elements.
void setSizeNoInit(exint newsize)
Definition: UT_Array.h:695
bool getAttributeAsArray(const GA_Attribute *atr, const GA_Range &range, UT_Array< T > &result) const
Get/set all the point attribute data from/to a contiguous array.
bool evaluateInteriorPoint(GA_Offset result_vtx, GA_AttributeRefMap &map, fpreal u, fpreal v, fpreal w=0) const
SYS_FORCE_INLINE T * SYSconst_cast(const T *foo)
Definition: SYS_Types.h:136
SYS_FORCE_INLINE bool getExtraFlag() const
Synonym for isClosed()
GA_DataId getDataId() const
Return the data ID for the topology attributes.
int64 exint
Definition: SYS_Types.h:125
A soup of polygons.
GA_Attribute * getP()
Convenience method to access the P attribute.
Definition: GA_Detail.h:164
void setTrivial(ToType startvalue, GA_Size size)
Makes the list a trivial list with the specified start value and size.
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
void setTrivialRange(FromType startindex, ToType startvalue, GA_Size nelements)
void setCapacity(exint new_capacity)
static const SOP_NodeVerb::Register< SOP_WindingNumberVerb > theVerb
SYS_FORCE_INLINE TO_T UTverify_cast(FROM_T from)
Definition: UT_Assert.h:229
bool addOperator(OP_Operator *op, std::ostream *err=nullptr)
SYS_FORCE_INLINE UT_Vector3 getPos3(GA_Offset ptoff) const
The ptoff passed is the point offset.
Definition: GA_Detail.h:185
static const UT_StringHolder theSOPTypeName
SYS_FORCE_INLINE bool isClosed() const
Definition: GEO_Face.h:248
static PRM_ChoiceList pointGroupMenu
Definition: SOP_Node.h:1190
exint size() const
Definition: UT_Array.h:646
exint GA_Size
Defines the bit width for index and offset types in GA.
Definition: GA_Types.h:235
UT_ErrorSeverity sopAddError(int code, const char *msg=0, const UT_SourceLocation *loc=0) const
Definition: SOP_NodeVerb.h:505
bool combine(const GA_Group *src) overridefinal
UT_SolidAngle< float, float > mySolidAngleTree
GA_Size GA_Offset
Definition: GA_Types.h:641
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:32
void getLocalTransform(UT_Matrix3D &x) const override
const T & parms() const
Definition: SOP_NodeVerb.h:412
GLdouble n
Definition: glcorearb.h:2008
GA_Range getPointRange(const GA_PointGroup *group=0) const
Get a range of all points in the detail.
Definition: GA_Detail.h:1730
Constructs a PRM_Template list from an embedded .ds file or an istream.
virtual SOP_NodeParms * allocParms() const
SYS_FORCE_INLINE void forEachPoint(FUNCTOR &&functor) const
Definition: GA_Detail.h:1499
PRM_Template * templates() const
SYS_FORCE_INLINE GA_OffsetListRef getPrimitiveVertexList(GA_Offset primoff) const
Definition: GA_Primitive.h:886
SYS_FORCE_INLINE int getNumRows() const
Definition: GEO_Hull.h:368
GLuint GLuint end
Definition: glcorearb.h:475
T computeSolidAngle(const UT_Vector3T< T > &query_point, const T accuracy_scale=T(2.0)) const
SYS_FORCE_INLINE GA_Offset getNumPointOffsets() const
Definition: GA_Detail.h:341
long long int64
Definition: SYS_Types.h:116
virtual SOP_NodeCache * allocCache() const
SYS_FORCE_INLINE const char * c_str() const
UT_Vector3T< fpreal64 > UT_Vector3D
SYS_FORCE_INLINE GA_Offset vertexPoint(GA_Offset vertex) const
Given a vertex, return the point it references.
Definition: GA_Detail.h:529
#define M_1_PI
Definition: fmath.h:102
SYS_FORCE_INLINE GA_DataId getDataId() const
Definition: GA_Attribute.h:299
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
GA_API const UT_StringHolder transform
GU_Detail * gdpNC()
GLenum GLenum GLsizei void * table
Definition: glad.h:5129
SYS_FORCE_INLINE bool atEnd() const
Definition: GA_PolyCounts.h:96
exint append()
Definition: UT_Array.h:142
exint getUniqueId() const
Definition: GA_Detail.h:117
GA_Topology & getTopology()
Definition: GA_Detail.h:798
SYS_FORCE_INLINE UT_Vector3 getPos3() const
The fast point position accessor.
Definition: GEO_Quadric.h:111
T computeAngle(const UT_Vector2T< T > &query_point, const T accuracy_scale=T(2.0)) const
SYS_FORCE_INLINE bool isValid() const
Definition: GA_Handle.h:187
SYS_FORCE_INLINE GA_Index pointIndex(GA_Offset offset) const
Given a point's data offset, return its index.
Definition: GA_Detail.h:349
SYS_FORCE_INLINE bool isTrivialMap() const
Definition: GA_IndexMap.h:290
const GA_Attribute * findFloatTuple(GA_AttributeOwner owner, GA_AttributeScope scope, const UT_StringRef &name, int min_size=1, int max_size=-1) const
const GA_PointGroup * parsePointDetached(const char *pat, const GEO_Detail *pgdp, bool forceexistence, bool &success)
UT_Array< UT_Vector2 > myPositions2D
void newSopOperator(OP_OperatorTable *table)
void UTparallelInvoke(bool parallel, F1 &&f1, F2 &&f2)
bool buildGrevilles(UT_Array< GEO_Greville > &dest) const
SYS_FORCE_INLINE void set(GA_Offset off, const T &val) const
Definition: GA_Handle.h:354
T UTsignedAngleSegment(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, const UT_Vector2T< T > &query)
SYS_FORCE_INLINE GA_Size getNumPrimitives() const
Return the number of primitives.
Definition: GA_Detail.h:408
int64 getMetaCacheCount() const
Definition: GA_Detail.h:2370
void update2D(const GA_Detail &mesh_geo, const GA_PrimitiveGroup *prim_group, const UT_StringHolder &group_string, const int approx_order, const int axis0, const int axis1)
const GU_Detail * inputGeo(exint idx) const
Definition: SOP_NodeVerb.h:376
const GA_PrimitiveList & getPrimitiveList() const
Definition: GA_Detail.h:792
const GA_PrimitiveGroup * parsePrimitiveDetached(const char *pat, const GEO_Detail *pgdp, bool forceexistence, bool &success)
SOP_NodeCache * cache() const
Definition: SOP_NodeVerb.h:417
Container class for all geometry.
Definition: GA_Detail.h:96
static const char *const theDsFile
This is the parameter interface string, below.
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:156
SYS_FORCE_INLINE int getNumCols() const
Definition: GEO_Hull.h:375
GLenum GLenum GLsizei void * row
Definition: glad.h:5135
void bind(GA_Detail *gdp, GA_AttributeOwner owner, const UT_StringRef &name, int minsize=1)
int64 GA_DataId
GLboolean r
Definition: glcorearb.h:1222
T UTsignedSolidAngleQuad(const UT_Vector3T< T > &a, const UT_Vector3T< T > &b, const UT_Vector3T< T > &c, const UT_Vector3T< T > &d, const UT_Vector3T< T > &query)
Definition: UT_SolidAngle.h:74
UT_SubtendedAngle< float, float > mySubtendedAngleTree
GA_Range getPrimitiveRange(const GA_PrimitiveGroup *group=0) const
Get a range of all primitives in the detail.
Definition: GA_Detail.h:1733
GU_DetailHandle & gdh() const
The initial state of gdh depends on the cookMode()
Definition: SOP_NodeVerb.h:341
bool fullBlockAdvance(GA_Offset &start, GA_Offset &end)
SYS_FORCE_INLINE bool isstring() const
UT_Array< int > myTrianglePoints
SYS_FORCE_INLINE FromType size() const
Returns the number of used elements in the list (always <= capacity())
void getUGrevilles(UT_Array< float > &ugrevilles) const
static SYS_FORCE_INLINE const int * fastFaceIndices(GA_Size faceno)
SYS_FORCE_INLINE GA_Size getNumPoints() const
Return the number of points.
Definition: GA_Detail.h:334
bool isEmpty() const
Returns true iff there are no occupied elements in the array.
Definition: UT_Array.h:650