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  * Copyright (c) 2021
3  * Side Effects Software Inc. All rights reserved.
4  *
5  * Redistribution and use of Houdini Development Kit samples in source and
6  * binary forms, with or without modification, are permitted provided that the
7  * following conditions are met:
8  * 1. Redistributions of source code must retain the above copyright notice,
9  * this list of conditions and the following disclaimer.
10  * 2. The name of Side Effects Software may not be used to endorse or
11  * promote products derived from this software without specific prior
12  * written permission.
13  *
14  * THIS SOFTWARE IS PROVIDED BY SIDE EFFECTS SOFTWARE `AS IS' AND ANY EXPRESS
15  * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
17  * NO EVENT SHALL SIDE EFFECTS SOFTWARE BE LIABLE FOR ANY DIRECT, INDIRECT,
18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
19  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
20  * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
21  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
22  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
23  * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  *
25  *----------------------------------------------------------------------------
26  * Functions for computing solid angles.
27  */
28 
29 #pragma once
30 
31 #ifndef __HDK_UT_SolidAngle_h__
32 #define __HDK_UT_SolidAngle_h__
33 
34 #include "UT_BVH.h"
35 
36 #include <UT/UT_UniquePtr.h>
37 #include <UT/UT_Vector3.h>
38 #include <SYS/SYS_Math.h>
39 
40 namespace HDK_Sample {
41 
42 /// Returns the signed solid angle subtended by triangle abc
43 /// from query point.
44 ///
45 /// WARNING: This uses the right-handed normal convention, whereas most of
46 /// Houdini uses the left-handed normal convention, so either
47 /// negate the output, or swap b and c if you want it to be
48 /// positive inside and negative outside.
49 template<typename T>
51  const UT_Vector3T<T> &a,
52  const UT_Vector3T<T> &b,
53  const UT_Vector3T<T> &c,
54  const UT_Vector3T<T> &query)
55 {
56  // Make a, b, and c relative to query
57  UT_Vector3T<T> qa = a-query;
58  UT_Vector3T<T> qb = b-query;
59  UT_Vector3T<T> qc = c-query;
60 
61  const T alength = qa.length();
62  const T blength = qb.length();
63  const T clength = qc.length();
64 
65  // If any triangle vertices are coincident with query,
66  // query is on the surface, which we treat as no solid angle.
67  if (alength == 0 || blength == 0 || clength == 0)
68  return T(0);
69 
70  // Normalize the vectors
71  qa /= alength;
72  qb /= blength;
73  qc /= clength;
74 
75  // The formula on Wikipedia has roughly dot(qa,cross(qb,qc)),
76  // but that's unstable when qa, qb, and qc are very close,
77  // (e.g. if the input triangle was very far away).
78  // This should be equivalent, but more stable.
79  const T numerator = dot(qa, cross(qb-qa, qc-qa));
80 
81  // If numerator is 0, regardless of denominator, query is on the
82  // surface, which we treat as no solid angle.
83  if (numerator == 0)
84  return T(0);
85 
86  const T denominator = T(1) + dot(qa,qb) + dot(qa,qc) + dot(qb,qc);
87 
88  return T(2)*SYSatan2(numerator, denominator);
89 }
90 
91 template<typename T>
93  const UT_Vector3T<T> &a,
94  const UT_Vector3T<T> &b,
95  const UT_Vector3T<T> &c,
96  const UT_Vector3T<T> &d,
97  const UT_Vector3T<T> &query)
98 {
99  // Make a, b, c, and d relative to query
100  UT_Vector3T<T> v[4] = {
101  a-query,
102  b-query,
103  c-query,
104  d-query
105  };
106 
107  const T lengths[4] = {
108  v[0].length(),
109  v[1].length(),
110  v[2].length(),
111  v[3].length()
112  };
113 
114  // If any quad vertices are coincident with query,
115  // query is on the surface, which we treat as no solid angle.
116  // We could add the contribution from the non-planar part,
117  // but in the context of a mesh, we'd still miss some, like
118  // we do in the triangle case.
119  if (lengths[0] == T(0) || lengths[1] == T(0) || lengths[2] == T(0) || lengths[3] == T(0))
120  return T(0);
121 
122  // Normalize the vectors
123  v[0] /= lengths[0];
124  v[1] /= lengths[1];
125  v[2] /= lengths[2];
126  v[3] /= lengths[3];
127 
128  // Compute (unnormalized, but consistently-scaled) barycentric coordinates
129  // for the query point inside the tetrahedron of points.
130  // If 0 or 4 of the coordinates are positive, (or slightly negative), the
131  // query is (approximately) inside, so the choice of triangulation matters.
132  // Otherwise, the triangulation doesn't matter.
133 
134  const UT_Vector3T<T> diag02 = v[2]-v[0];
135  const UT_Vector3T<T> diag13 = v[3]-v[1];
136  const UT_Vector3T<T> v01 = v[1]-v[0];
137  const UT_Vector3T<T> v23 = v[3]-v[2];
138 
139  T bary[4];
140  bary[0] = dot(v[3],cross(v23,diag13));
141  bary[1] = -dot(v[2],cross(v23,diag02));
142  bary[2] = -dot(v[1],cross(v01,diag13));
143  bary[3] = dot(v[0],cross(v01,diag02));
144 
145  const T dot01 = dot(v[0],v[1]);
146  const T dot12 = dot(v[1],v[2]);
147  const T dot23 = dot(v[2],v[3]);
148  const T dot30 = dot(v[3],v[0]);
149 
150  T omega = T(0);
151 
152  // Equation of a bilinear patch in barycentric coordinates of its
153  // tetrahedron is x0*x2 = x1*x3. Less is one side; greater is other.
154  if (bary[0]*bary[2] < bary[1]*bary[3])
155  {
156  // Split 0-2: triangles 0,1,2 and 0,2,3
157  const T numerator012 = bary[3];
158  const T numerator023 = bary[1];
159  const T dot02 = dot(v[0],v[2]);
160 
161  // If numerator is 0, regardless of denominator, query is on the
162  // surface, which we treat as no solid angle.
163  if (numerator012 != T(0))
164  {
165  const T denominator012 = T(1) + dot01 + dot12 + dot02;
166  omega = SYSatan2(numerator012, denominator012);
167  }
168  if (numerator023 != T(0))
169  {
170  const T denominator023 = T(1) + dot02 + dot23 + dot30;
171  omega += SYSatan2(numerator023, denominator023);
172  }
173  }
174  else
175  {
176  // Split 1-3: triangles 0,1,3 and 1,2,3
177  const T numerator013 = -bary[2];
178  const T numerator123 = -bary[0];
179  const T dot13 = dot(v[1],v[3]);
180 
181  // If numerator is 0, regardless of denominator, query is on the
182  // surface, which we treat as no solid angle.
183  if (numerator013 != T(0))
184  {
185  const T denominator013 = T(1) + dot01 + dot13 + dot30;
186  omega = SYSatan2(numerator013, denominator013);
187  }
188  if (numerator123 != T(0))
189  {
190  const T denominator123 = T(1) + dot12 + dot23 + dot13;
191  omega += SYSatan2(numerator123, denominator123);
192  }
193  }
194  return T(2)*omega;
195 }
196 
197 /// Class for quickly approximating signed solid angle of a large mesh
198 /// from many query points. This is useful for computing the
199 /// generalized winding number at many points.
200 ///
201 /// NOTE: This is currently only instantiated for <float,float>.
202 template<typename T,typename S>
204 {
205 public:
206  /// This is outlined so that we don't need to include UT_BVHImpl.h
207  UT_SolidAngle();
208  /// This is outlined so that we don't need to include UT_BVHImpl.h
209  ~UT_SolidAngle();
210 
211  /// NOTE: This does not take ownership over triangle_points or positions,
212  /// but does keep pointers to them, so the caller must keep them in
213  /// scope for the lifetime of this structure.
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  : UT_SolidAngle()
221  { init(ntriangles, triangle_points, npoints, positions, order); }
222 
223  /// Initialize the tree and data.
224  /// NOTE: It is safe to call init on a UT_SolidAngle that has had init
225  /// called on it before, to re-initialize it.
226  void init(
227  const int ntriangles,
228  const int *const triangle_points,
229  const int npoints,
230  const UT_Vector3T<S> *const positions,
231  const int order = 2);
232 
233  /// Frees myTree and myData, and clears the rest.
234  void clear();
235 
236  /// Returns true if this is clear
237  bool isClear() const
238  { return myNTriangles == 0; }
239 
240  /// Returns an approximation of the signed solid angle of the mesh from the specified query_point
241  /// accuracy_scale is the value of (maxP/q) beyond which the approximation of the box will be used.
242  T computeSolidAngle(const UT_Vector3T<T> &query_point, const T accuracy_scale = T(2.0)) const;
243 
244 private:
245  struct BoxData;
246 
247  static constexpr uint BVH_N = 4;
248  UT_BVH<BVH_N> myTree;
249  int myNBoxes;
250  int myOrder;
252  int myNTriangles;
253  const int *myTrianglePoints;
254  int myNPoints;
255  const UT_Vector3T<S> *myPositions;
256 };
257 
258 template<typename T>
260  const UT_Vector2T<T> &a,
261  const UT_Vector2T<T> &b,
262  const UT_Vector2T<T> &query)
263 {
264  // Make a and b relative to query
265  UT_Vector2T<T> qa = a-query;
266  UT_Vector2T<T> qb = b-query;
267 
268  // If any segment vertices are coincident with query,
269  // query is on the segment, which we treat as no angle.
270  if (qa.isZero() || qb.isZero())
271  return T(0);
272 
273  // numerator = |qa||qb|sin(theta)
274  const T numerator = cross(qa, qb);
275 
276  // If numerator is 0, regardless of denominator, query is on the
277  // surface, which we treat as no solid angle.
278  if (numerator == 0)
279  return T(0);
280 
281  // denominator = |qa||qb|cos(theta)
282  const T denominator = dot(qa,qb);
283 
284  // numerator/denominator = tan(theta)
285  return SYSatan2(numerator, denominator);
286 }
287 
288 /// Class for quickly approximating signed subtended angle of a large curve
289 /// from many query points. This is useful for computing the
290 /// generalized winding number at many points.
291 ///
292 /// NOTE: This is currently only instantiated for <float,float>.
293 template<typename T,typename S>
295 {
296 public:
297  /// This is outlined so that we don't need to include UT_BVHImpl.h
299  /// This is outlined so that we don't need to include UT_BVHImpl.h
301 
302  /// NOTE: This does not take ownership over segment_points or positions,
303  /// but does keep pointers to them, so the caller must keep them in
304  /// scope for the lifetime of this structure.
306  const int nsegments,
307  const int *const segment_points,
308  const int npoints,
309  const UT_Vector2T<S> *const positions,
310  const int order = 2)
312  { init(nsegments, segment_points, npoints, positions, order); }
313 
314  /// Initialize the tree and data.
315  /// NOTE: It is safe to call init on a UT_SolidAngle that has had init
316  /// called on it before, to re-initialize it.
317  void init(
318  const int nsegments,
319  const int *const segment_points,
320  const int npoints,
321  const UT_Vector2T<S> *const positions,
322  const int order = 2);
323 
324  /// Frees myTree and myData, and clears the rest.
325  void clear();
326 
327  /// Returns true if this is clear
328  bool isClear() const
329  { return myNSegments == 0; }
330 
331  /// Returns an approximation of the signed solid angle of the mesh from the specified query_point
332  /// accuracy_scale is the value of (maxP/q) beyond which the approximation of the box will be used.
333  T computeAngle(const UT_Vector2T<T> &query_point, const T accuracy_scale = T(2.0)) const;
334 
335 private:
336  struct BoxData;
337 
338  static constexpr uint BVH_N = 4;
339  UT_BVH<BVH_N> myTree;
340  int myNBoxes;
341  int myOrder;
343  int myNSegments;
344  const int *mySegmentPoints;
345  int myNPoints;
346  const UT_Vector2T<S> *myPositions;
347 };
348 
349 } // End HDK_Sample namespace
350 #endif
UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
UT_SolidAngle(const int ntriangles, const int *const triangle_points, const int npoints, const UT_Vector3T< S > *const positions, const int order=2)
int myOrder
Definition: GT_CurveEval.h:263
SYS_FORCE_INLINE bool isZero() const noexcept
GLboolean GLboolean GLboolean GLboolean a
Definition: glew.h:9477
GLenum query
Definition: glew.h:5704
const GLdouble * v
Definition: glew.h:1391
3D Vector class.
2D Vector class.
Definition: UT_Vector2.h:149
~UT_SolidAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
bool isClear() const
Returns true if this is clear.
GLuint GLdouble GLdouble GLint GLint order
Definition: glew.h:3446
std::unique_ptr< T, Deleter > UT_UniquePtr
A smart pointer for unique ownership of dynamically allocated objects.
Definition: UT_UniquePtr.h:33
INT32 * numerator
Definition: wglew.h:1180
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
INT32 INT32 * denominator
Definition: wglew.h:1180
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:127
void clear()
Frees myTree and myData, and clears the rest.
void init(const int nsegments, const int *const segment_points, const int npoints, const UT_Vector2T< S > *const positions, const int order=2)
T computeSolidAngle(const UT_Vector3T< T > &query_point, const T accuracy_scale=T(2.0)) const
const GLfloat * c
Definition: glew.h:16296
void init(const int ntriangles, const int *const triangle_points, const int npoints, const UT_Vector3T< S > *const positions, const int order=2)
GLsizei GLenum GLuint GLuint GLsizei * lengths
Definition: glew.h:2581
T computeAngle(const UT_Vector2T< T > &query_point, const T accuracy_scale=T(2.0)) const
T UTsignedAngleSegment(const UT_Vector2T< T > &a, const UT_Vector2T< T > &b, const UT_Vector2T< T > &query)
GLdouble GLdouble GLdouble b
Definition: glew.h:9122
bool isClear() const
Returns true if this is clear.
UT_SubtendedAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
UT_SubtendedAngle(const int nsegments, const int *const segment_points, const int npoints, const UT_Vector2T< S > *const positions, const int order=2)
void clear()
Frees myTree and myData, and clears the rest.
~UT_SubtendedAngle()
This is outlined so that we don't need to include UT_BVHImpl.h.
SYS_FORCE_INLINE Storage::MathFloat length() const
unsigned int uint
Definition: SYS_Types.h:45
T UTsignedSolidAngleTri(const UT_Vector3T< T > &a, const UT_Vector3T< T > &b, const UT_Vector3T< T > &c, const UT_Vector3T< T > &query)
Definition: UT_SolidAngle.h:50
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)