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