HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_SolidAngle.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_SolidAngle.h (UT Library, C++)
7  *
8  * COMMENTS: Functions for computing solid angles.
9  */
10 
11 #ifndef __UT_SolidAngle_h__
12 #define __UT_SolidAngle_h__
13 
14 #include "UT_API.h"
15 #include "UT_UniquePtr.h"
16 #include "UT_Vector3.h"
17 #include <SYS/SYS_Math.h>
18 
19 #define SOLID_ANGLE_USE_BVH 1
20 #if SOLID_ANGLE_USE_BVH
21 #include "UT_BVH.h"
22 #endif
23 
24 /// Returns the signed solid angle subtended by triangle abc
25 /// from query point.
26 ///
27 /// WARNING: This uses the right-handed normal convention, whereas most of
28 /// Houdini uses the left-handed normal convention, so either
29 /// negate the output, or swap b and c if you want it to be
30 /// positive inside and negative outside.
31 template<typename T>
33  const UT_Vector3T<T> &a,
34  const UT_Vector3T<T> &b,
35  const UT_Vector3T<T> &c,
36  const UT_Vector3T<T> &query)
37 {
38  // Make a, b, and c relative to query
39  UT_Vector3T<T> qa = a-query;
40  UT_Vector3T<T> qb = b-query;
41  UT_Vector3T<T> qc = c-query;
42 
43  const T alength = qa.length();
44  const T blength = qb.length();
45  const T clength = qc.length();
46 
47  // If any triangle vertices are coincident with query,
48  // query is on the surface, which we treat as no solid angle.
49  if (alength == 0 || blength == 0 || clength == 0)
50  return T(0);
51 
52  // Normalize the vectors
53  qa /= alength;
54  qb /= blength;
55  qc /= clength;
56 
57  // The formula on Wikipedia has roughly dot(qa,cross(qb,qc)),
58  // but that's unstable when qa, qb, and qc are very close,
59  // (e.g. if the input triangle was very far away).
60  // This should be equivalent, but more stable.
61  const T numerator = dot(qa, cross(qb-qa, qc-qa));
62 
63  // If numerator is 0, regardless of denominator, query is on the
64  // surface, which we treat as no solid angle.
65  if (numerator == 0)
66  return T(0);
67 
68  const T denominator = T(1) + dot(qa,qb) + dot(qa,qc) + dot(qb,qc);
69 
70  return T(2)*SYSatan2(numerator, denominator);
71 }
72 
73 template<typename T>
75  const UT_Vector3T<T> &a,
76  const UT_Vector3T<T> &b,
77  const UT_Vector3T<T> &c,
78  const UT_Vector3T<T> &d,
79  const UT_Vector3T<T> &query)
80 {
81  // Make a, b, c, and d relative to query
82  UT_Vector3T<T> v[4] = {
83  a-query,
84  b-query,
85  c-query,
86  d-query
87  };
88 
89  const T lengths[4] = {
90  v[0].length(),
91  v[1].length(),
92  v[2].length(),
93  v[3].length()
94  };
95 
96  // If any quad vertices are coincident with query,
97  // query is on the surface, which we treat as no solid angle.
98  // We could add the contribution from the non-planar part,
99  // but in the context of a mesh, we'd still miss some, like
100  // we do in the triangle case.
101  if (lengths[0] == T(0) || lengths[1] == T(0) || lengths[2] == T(0) || lengths[3] == T(0))
102  return T(0);
103 
104  // Normalize the vectors
105  v[0] /= lengths[0];
106  v[1] /= lengths[1];
107  v[2] /= lengths[2];
108  v[3] /= lengths[3];
109 
110  // Compute (unnormalized, but consistently-scaled) barycentric coordinates
111  // for the query point inside the tetrahedron of points.
112  // If 0 or 4 of the coordinates are positive, (or slightly negative), the
113  // query is (approximately) inside, so the choice of triangulation matters.
114  // Otherwise, the triangulation doesn't matter.
115 
116  const UT_Vector3T<T> diag02 = v[2]-v[0];
117  const UT_Vector3T<T> diag13 = v[3]-v[1];
118  const UT_Vector3T<T> v01 = v[1]-v[0];
119  const UT_Vector3T<T> v23 = v[3]-v[2];
120 
121  T bary[4];
122  bary[0] = dot(v[3],cross(v23,diag13));
123  bary[1] = -dot(v[2],cross(v23,diag02));
124  bary[2] = -dot(v[1],cross(v01,diag13));
125  bary[3] = dot(v[0],cross(v01,diag02));
126 
127  const T dot01 = dot(v[0],v[1]);
128  const T dot12 = dot(v[1],v[2]);
129  const T dot23 = dot(v[2],v[3]);
130  const T dot30 = dot(v[3],v[0]);
131 
132  T omega = T(0);
133 
134  // Equation of a bilinear patch in barycentric coordinates of its
135  // tetrahedron is x0*x2 = x1*x3. Less is one side; greater is other.
136  if (bary[0]*bary[2] < bary[1]*bary[3])
137  {
138  // Split 0-2: triangles 0,1,2 and 0,2,3
139  const T numerator012 = bary[3];
140  const T numerator023 = bary[1];
141  const T dot02 = dot(v[0],v[2]);
142 
143  // If numerator is 0, regardless of denominator, query is on the
144  // surface, which we treat as no solid angle.
145  if (numerator012 != T(0))
146  {
147  const T denominator012 = T(1) + dot01 + dot12 + dot02;
148  omega = SYSatan2(numerator012, denominator012);
149  }
150  if (numerator023 != T(0))
151  {
152  const T denominator023 = T(1) + dot02 + dot23 + dot30;
153  omega += SYSatan2(numerator023, denominator023);
154  }
155  }
156  else
157  {
158  // Split 1-3: triangles 0,1,3 and 1,2,3
159  const T numerator013 = -bary[2];
160  const T numerator123 = -bary[0];
161  const T dot13 = dot(v[1],v[3]);
162 
163  // If numerator is 0, regardless of denominator, query is on the
164  // surface, which we treat as no solid angle.
165  if (numerator013 != T(0))
166  {
167  const T denominator013 = T(1) + dot01 + dot13 + dot30;
168  omega = SYSatan2(numerator013, denominator013);
169  }
170  if (numerator123 != T(0))
171  {
172  const T denominator123 = T(1) + dot12 + dot23 + dot13;
173  omega += SYSatan2(numerator123, denominator123);
174  }
175  }
176  return T(2)*omega;
177 }
178 
179 #if !SOLID_ANGLE_USE_BVH
180 // Forward declaration, to avoid having to include UT_RTree.h
181 template<int MAX_ORDER>
182 class UT_RTreeGeneric;
183 #endif
184 
185 /// Class for quickly approximating signed solid angle of a large mesh
186 /// from many query points. This is useful for computing the
187 /// generalized winding number at many points.
188 ///
189 /// NOTE: This is currently only instantiated for <float,float>.
190 template<typename T,typename S>
192 {
193 public:
194  /// This is outlined so that we don't need to include UT_BVHImpl.h
195  UT_SolidAngle();
196  /// This is outlined so that we don't need to include UT_BVHImpl.h
197  ~UT_SolidAngle();
198 
199  /// NOTE: This does not take ownership over triangle_points or positions,
200  /// but does keep pointers to them, so the caller must keep them in
201  /// scope for the lifetime of this structure.
203  const int ntriangles,
204  const int *const triangle_points,
205  const int npoints,
206  const UT_Vector3T<S> *const positions,
207  const int order = 2)
208  : UT_SolidAngle()
209  { init(ntriangles, triangle_points, npoints, positions, order); }
210 
211  /// Initialize the tree and data.
212  /// NOTE: It is safe to call init on a UT_SolidAngle that has had init
213  /// called on it before, to re-initialize it.
214  void init(
215  const int ntriangles,
216  const int *const triangle_points,
217  const int npoints,
218  const UT_Vector3T<S> *const positions,
219  const int order = 2);
220 
221  /// Frees myTree and myData, and clears the rest.
222  void clear();
223 
224  /// Returns true if this is clear
225  bool isClear() const
226  { return myNTriangles == 0; }
227 
228  /// Returns an approximation of the signed solid angle of the mesh from the specified query_point
229  /// accuracy_scale is the value of (maxP/q) beyond which the approximation of the box will be used.
230  T computeSolidAngle(const UT_Vector3T<T> &query_point, const T accuracy_scale = T(2.0)) const;
231 
232 private:
233  struct BoxData;
234 
235 #if SOLID_ANGLE_USE_BVH
236  static constexpr uint BVH_N = 4;
237  UT_BVH<BVH_N> myTree;
238 #else
240 #endif
241  int myNBoxes;
242  int myOrder;
244  int myNTriangles;
245  const int *myTrianglePoints;
246  int myNPoints;
247  const UT_Vector3T<S> *myPositions;
248 };
249 
250 #undef SOLID_ANGLE_USE_BVH
251 
252 template<typename T>
254  const UT_Vector2T<T> &a,
255  const UT_Vector2T<T> &b,
256  const UT_Vector2T<T> &query)
257 {
258  // Make a and b relative to query
259  UT_Vector2T<T> qa = a-query;
260  UT_Vector2T<T> qb = b-query;
261 
262  // If any segment vertices are coincident with query,
263  // query is on the segment, which we treat as no angle.
264  if (qa.isZero() || qb.isZero())
265  return T(0);
266 
267  // numerator = |qa||qb|sin(theta)
268  const T numerator = cross(qa, qb);
269 
270  // If numerator is 0, regardless of denominator, query is on the
271  // surface, which we treat as no solid angle.
272  if (numerator == 0)
273  return T(0);
274 
275  // denominator = |qa||qb|cos(theta)
276  const T denominator = dot(qa,qb);
277 
278  // numerator/denominator = tan(theta)
279  return SYSatan2(numerator, denominator);
280 }
281 
282 /// Class for quickly approximating signed subtended angle of a large curve
283 /// from many query points. This is useful for computing the
284 /// generalized winding number at many points.
285 ///
286 /// NOTE: This is currently only instantiated for <float,float>.
287 template<typename T,typename S>
289 {
290 public:
291  /// This is outlined so that we don't need to include UT_BVHImpl.h
293  /// This is outlined so that we don't need to include UT_BVHImpl.h
295 
296  /// NOTE: This does not take ownership over segment_points or positions,
297  /// but does keep pointers to them, so the caller must keep them in
298  /// scope for the lifetime of this structure.
300  const int nsegments,
301  const int *const segment_points,
302  const int npoints,
303  const UT_Vector2T<S> *const positions,
304  const int order = 2)
306  { init(nsegments, segment_points, npoints, positions, order); }
307 
308  /// Initialize the tree and data.
309  /// NOTE: It is safe to call init on a UT_SolidAngle that has had init
310  /// called on it before, to re-initialize it.
311  void init(
312  const int nsegments,
313  const int *const segment_points,
314  const int npoints,
315  const UT_Vector2T<S> *const positions,
316  const int order = 2);
317 
318  /// Frees myTree and myData, and clears the rest.
319  void clear();
320 
321  /// Returns true if this is clear
322  bool isClear() const
323  { return myNSegments == 0; }
324 
325  /// Returns an approximation of the signed solid angle of the mesh from the specified query_point
326  /// accuracy_scale is the value of (maxP/q) beyond which the approximation of the box will be used.
327  T computeAngle(const UT_Vector2T<T> &query_point, const T accuracy_scale = T(2.0)) const;
328 
329 private:
330  struct BoxData;
331 
332  static constexpr uint BVH_N = 4;
333  UT_BVH<BVH_N> myTree;
334  int myNBoxes;
335  int myOrder;
337  int myNSegments;
338  const int *mySegmentPoints;
339  int myNPoints;
340  const UT_Vector2T<S> *myPositions;
341 };
342 
343 #endif
GLenum query
Definition: glad.h:2772
int myOrder
Definition: GT_CurveEval.h:263
const GLdouble * v
Definition: glcorearb.h:837
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
#define UT_API
Definition: UT_API.h:14
constexpr SYS_FORCE_INLINE T length() const noexcept
Definition: UT_Vector3.h:361
3D Vector class.
2D Vector class.
Definition: UT_Vector2.h:159
UT_SubtendedAngle(const int nsegments, const int *const segment_points, const int npoints, const UT_Vector2T< S > *const positions, const int order=2)
std::unique_ptr< T, Deleter > UT_UniquePtr
A smart pointer for unique ownership of dynamically allocated objects.
Definition: UT_UniquePtr.h:39
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
UT_SolidAngle(const int ntriangles, const int *const triangle_points, const int npoints, const UT_Vector3T< S > *const positions, const int order=2)
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:130
GLdouble GLdouble GLint GLint order
Definition: glad.h:2676
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
constexpr SYS_FORCE_INLINE bool isZero() const noexcept
Definition: UT_Vector2.h:326
bool isClear() const
Returns true if this is clear.
bool isClear() const
Returns true if this is clear.
T UTsignedAngleSegment(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, const UT_Vector2T< T > &query)
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
unsigned int uint
Definition: SYS_Types.h:45
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
GLsizei GLenum GLenum GLuint GLenum GLsizei * lengths
Definition: glcorearb.h:2542