OpenVDB  9.1.1
Stencils.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @author Ken Museth
5 ///
6 /// @date April 9, 2021
7 ///
8 /// @file Stencils.h
9 ///
10 /// @brief Defines various finite-difference stencils that allow for the
11 /// computation of gradients of order 1 to 5, mean curvatures,
12 /// gaussian curvatures, principal curvatures, tri-linear interpolation,
13 /// zero-crossing, laplacian, and closest point transform.
14 
15 #ifndef NANOVDB_STENCILS_HAS_BEEN_INCLUDED
16 #define NANOVDB_STENCILS_HAS_BEEN_INCLUDED
17 
18 #include "../NanoVDB.h"// for __hosedev__, Vec3, Min, Max, Pow2, Pow3, Pow4
19 
20 namespace nanovdb {
21 
22 // ---------------------------- WENO5 ----------------------------
23 
24 /// @brief Implementation of nominally fifth-order finite-difference WENO
25 /// @details This function returns the numerical flux. See "High Order Finite Difference and
26 /// Finite Volume WENO Schemes and Discontinuous Galerkin Methods for CFD" - Chi-Wang Shu
27 /// ICASE Report No 2001-11 (page 6). Also see ICASE No 97-65 for a more complete reference
28 /// (Shu, 1997).
29 /// Given v1 = f(x-2dx), v2 = f(x-dx), v3 = f(x), v4 = f(x+dx) and v5 = f(x+2dx),
30 /// return an interpolated value f(x+dx/2) with the special property that
31 /// ( f(x+dx/2) - f(x-dx/2) ) / dx = df/dx (x) + error,
32 /// where the error is fifth-order in smooth regions: O(dx) <= error <=O(dx^5)
33 template<typename ValueType, typename RealT = ValueType>
34 __hostdev__ inline ValueType
35 WENO5(const ValueType& v1,
36  const ValueType& v2,
37  const ValueType& v3,
38  const ValueType& v4,
39  const ValueType& v5,
40  RealT scale2 = 1.0)// openvdb uses scale2 = 0.01
41 {
42  static const RealT C = 13.0 / 12.0;
43  // WENO is formulated for non-dimensional equations, here the optional scale2
44  // is a reference value (squared) for the function being interpolated. For
45  // example if 'v' is of order 1000, then scale2 = 10^6 is ok. But in practice
46  // leave scale2 = 1.
47  const RealT eps = RealT(1.0e-6) * scale2;
48  // {\tilde \omega_k} = \gamma_k / ( \beta_k + \epsilon)^2 in Shu's ICASE report)
49  const RealT A1 = RealT(0.1)/Pow2(C*Pow2(v1-2*v2+v3)+RealT(0.25)*Pow2(v1-4*v2+3*v3)+eps),
50  A2 = RealT(0.6)/Pow2(C*Pow2(v2-2*v3+v4)+RealT(0.25)*Pow2(v2-v4)+eps),
51  A3 = RealT(0.3)/Pow2(C*Pow2(v3-2*v4+v5)+RealT(0.25)*Pow2(3*v3-4*v4+v5)+eps);
52 
53  return static_cast<ValueType>((A1*(2*v1 - 7*v2 + 11*v3) +
54  A2*(5*v3 - v2 + 2*v4) +
55  A3*(2*v3 + 5*v4 - v5))/(6*(A1+A2+A3)));
56 }
57 
58 // ---------------------------- GodunovsNormSqrd ----------------------------
59 
60 template <typename RealT>
61 __hostdev__ inline RealT
62 GodunovsNormSqrd(bool isOutside,
63  RealT dP_xm, RealT dP_xp,
64  RealT dP_ym, RealT dP_yp,
65  RealT dP_zm, RealT dP_zp)
66 {
67  RealT dPLen2;
68  if (isOutside) { // outside
69  dPLen2 = Max(Pow2(Max(dP_xm, RealT(0))), Pow2(Min(dP_xp, RealT(0)))); // (dP/dx)2
70  dPLen2 += Max(Pow2(Max(dP_ym, RealT(0))), Pow2(Min(dP_yp, RealT(0)))); // (dP/dy)2
71  dPLen2 += Max(Pow2(Max(dP_zm, RealT(0))), Pow2(Min(dP_zp, RealT(0)))); // (dP/dz)2
72  } else { // inside
73  dPLen2 = Max(Pow2(Min(dP_xm, RealT(0))), Pow2(Max(dP_xp, RealT(0)))); // (dP/dx)2
74  dPLen2 += Max(Pow2(Min(dP_ym, RealT(0))), Pow2(Max(dP_yp, RealT(0)))); // (dP/dy)2
75  dPLen2 += Max(Pow2(Min(dP_zm, RealT(0))), Pow2(Max(dP_zp, RealT(0)))); // (dP/dz)2
76  }
77  return dPLen2; // |\nabla\phi|^2
78 }
79 
80 template<typename RealT>
81 __hostdev__ inline RealT
82 GodunovsNormSqrd(bool isOutside,
83  const Vec3<RealT>& gradient_m,
84  const Vec3<RealT>& gradient_p)
85 {
86  return GodunovsNormSqrd<RealT>(isOutside,
87  gradient_m[0], gradient_p[0],
88  gradient_m[1], gradient_p[1],
89  gradient_m[2], gradient_p[2]);
90 }
91 
92 // ---------------------------- BaseStencil ----------------------------
93 
94 // BaseStencil uses curiously recurring template pattern (CRTP)
95 template<typename DerivedType, int SIZE, typename GridT>
97 {
98 public:
99  using ValueType = typename GridT::ValueType;
100  using GridType = GridT;
101  using TreeType = typename GridT::TreeType;
103 
104  /// @brief Initialize the stencil buffer with the values of voxel (i, j, k)
105  /// and its neighbors.
106  /// @param ijk Index coordinates of stencil center
107  __hostdev__ inline void moveTo(const Coord& ijk)
108  {
109  mCenter = ijk;
110  mValues[0] = mAcc.getValue(ijk);
111  static_cast<DerivedType&>(*this).init(mCenter);
112  }
113 
114  /// @brief Initialize the stencil buffer with the values of voxel (i, j, k)
115  /// and its neighbors. The method also takes a value of the center
116  /// element of the stencil, assuming it is already known.
117  /// @param ijk Index coordinates of stencil center
118  /// @param centerValue Value of the center element of the stencil
119  __hostdev__ inline void moveTo(const Coord& ijk, const ValueType& centerValue)
120  {
121  mCenter = ijk;
122  mValues[0] = centerValue;
123  static_cast<DerivedType&>(*this).init(mCenter);
124  }
125 
126  /// @brief Initialize the stencil buffer with the values of voxel
127  /// (x, y, z) and its neighbors.
128  ///
129  /// @note This version is slightly faster than the one above, since
130  /// the center voxel's value is read directly from the iterator.
131  template<typename IterType>
132  __hostdev__ inline void moveTo(const IterType& iter)
133  {
134  mCenter = iter.getCoord();
135  mValues[0] = *iter;
136  static_cast<DerivedType&>(*this).init(mCenter);
137  }
138 
139  /// @brief Initialize the stencil buffer with the values of voxel (x, y, z)
140  /// and its neighbors.
141  /// @param xyz Floating point voxel coordinates of stencil center
142  /// @details This method will check to see if it is necessary to
143  /// update the stencil based on the cached index coordinates of
144  /// the center point.
145  template<typename RealType>
146  __hostdev__ inline void moveTo(const Vec3<RealType>& xyz)
147  {
148  Coord ijk = RoundDown(xyz);
149  if (ijk != mCenter) this->moveTo(ijk);
150  }
151 
152  /// @brief Return the value from the stencil buffer with linear
153  /// offset pos.
154  ///
155  /// @note The default (@a pos = 0) corresponds to the first element
156  /// which is typically the center point of the stencil.
157  __hostdev__ inline const ValueType& getValue(unsigned int pos = 0) const
158  {
159  NANOVDB_ASSERT(pos < SIZE);
160  return mValues[pos];
161  }
162 
163  /// @brief Return the value at the specified location relative to the center of the stencil
164  template<int i, int j, int k>
165  __hostdev__ inline const ValueType& getValue() const
166  {
167  return mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()];
168  }
169 
170  /// @brief Set the value at the specified location relative to the center of the stencil
171  template<int i, int j, int k>
172  __hostdev__ inline void setValue(const ValueType& value)
173  {
174  mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()] = value;
175  }
176 
177  /// @brief Return the size of the stencil buffer.
178  __hostdev__ static int size() { return SIZE; }
179 
180  /// @brief Return the mean value of the current stencil.
181  __hostdev__ inline ValueType mean() const
182  {
183  ValueType sum = 0.0;
184  for (int i = 0; i < SIZE; ++i) sum += mValues[i];
185  return sum / ValueType(SIZE);
186  }
187 
188  /// @brief Return the smallest value in the stencil buffer.
189  __hostdev__ inline ValueType min() const
190  {
191  ValueType v = mValues[0];
192  for (int i=1; i<SIZE; ++i) {
193  if (mValues[i] < v) v = mValues[i];
194  }
195  return v;
196  }
197 
198  /// @brief Return the largest value in the stencil buffer.
199  __hostdev__ inline ValueType max() const
200  {
201  ValueType v = mValues[0];
202  for (int i=1; i<SIZE; ++i) {
203  if (mValues[i] > v) v = mValues[i];
204  }
205  return v;
206  }
207 
208  /// @brief Return the coordinates of the center point of the stencil.
209  __hostdev__ inline const Coord& getCenterCoord() const { return mCenter; }
210 
211  /// @brief Return the value at the center of the stencil
212  __hostdev__ inline const ValueType& getCenterValue() const { return mValues[0]; }
213 
214  /// @brief Return true if the center of the stencil intersects the
215  /// iso-contour specified by the isoValue
216  __hostdev__ inline bool intersects(const ValueType &isoValue = ValueType(0) ) const
217  {
218  const bool less = this->getValue< 0, 0, 0>() < isoValue;
219  return (less ^ (this->getValue<-1, 0, 0>() < isoValue)) ||
220  (less ^ (this->getValue< 1, 0, 0>() < isoValue)) ||
221  (less ^ (this->getValue< 0,-1, 0>() < isoValue)) ||
222  (less ^ (this->getValue< 0, 1, 0>() < isoValue)) ||
223  (less ^ (this->getValue< 0, 0,-1>() < isoValue)) ||
224  (less ^ (this->getValue< 0, 0, 1>() < isoValue)) ;
225  }
226  struct Mask {
227  uint8_t bits;
228  Mask() : bits(0u) {}
229  void set(int i) { bits |= (1 << i); }
230  bool test(int i) const { return bits & (1 << i); }
231  bool any() const { return bits > 0u; }
232  bool all() const { return bits == 255u; }
233  bool none() const { return bits == 0u; }
234  int count() const { return CountOn(bits); }
235  };// Mask
236 
237  /// @brief Return true a bit-mask where the 6 lower bits indicates if the
238  /// center of the stencil intersects the iso-contour specified by the isoValue.
239  ///
240  /// @note There are 2^6 = 64 different possible cases, including no intersections!
241  ///
242  /// @details The ordering of bit mask is ( -x, +x, -y, +y, -z, +z ), so to
243  /// check if there is an intersection in -y use (mask & (1u<<2)) where mask is
244  /// ther return value from this function. To check if there are any
245  /// intersections use mask!=0u, and for no intersections use mask==0u.
246  /// To count the number of intersections use __builtin_popcount(mask).
247  inline Mask intersectionMask(ValueType isoValue = ValueType(0)) const
248  {
249  Mask mask;
250  const bool less = this->getValue< 0, 0, 0>() < isoValue;
251  if (less ^ (this->getValue<-1, 0, 0>() < isoValue)) mask.set(0);// |= 1u;
252  if (less ^ (this->getValue< 1, 0, 0>() < isoValue)) mask.set(1);// |= 2u;
253  if (less ^ (this->getValue< 0,-1, 0>() < isoValue)) mask.set(2);// |= 4u;
254  if (less ^ (this->getValue< 0, 1, 0>() < isoValue)) mask.set(3);// |= 8u;
255  if (less ^ (this->getValue< 0, 0,-1>() < isoValue)) mask.set(4);// |= 16u;
256  if (less ^ (this->getValue< 0, 0, 1>() < isoValue)) mask.set(5);// |= 32u;
257  return mask;
258  }
259 
260  /// @brief Return a const reference to the grid from which this
261  /// stencil was constructed.
262  __hostdev__ inline const GridType& grid() const { return *mGrid; }
263 
264  /// @brief Return a const reference to the ValueAccessor
265  /// associated with this Stencil.
266  __hostdev__ inline const AccessorType& accessor() const { return mAcc; }
267 
268 protected:
269  // Constructor is protected to prevent direct instantiation.
271  : mGrid(&grid)
272  , mAcc(grid.tree().root())
273  , mCenter(Coord::max())
274  {
275  }
276 
277  const GridType* mGrid;
281 
282 }; // BaseStencil class
283 
284 
285 // ---------------------------- BoxStencil ----------------------------
286 
287 
288 namespace { // anonymous namespace for stencil-layout map
289 
290  // the eight point box stencil
291  template<int i, int j, int k> struct BoxPt {};
292  template<> struct BoxPt< 0, 0, 0> { enum { idx = 0 }; };
293  template<> struct BoxPt< 0, 0, 1> { enum { idx = 1 }; };
294  template<> struct BoxPt< 0, 1, 1> { enum { idx = 2 }; };
295  template<> struct BoxPt< 0, 1, 0> { enum { idx = 3 }; };
296  template<> struct BoxPt< 1, 0, 0> { enum { idx = 4 }; };
297  template<> struct BoxPt< 1, 0, 1> { enum { idx = 5 }; };
298  template<> struct BoxPt< 1, 1, 1> { enum { idx = 6 }; };
299  template<> struct BoxPt< 1, 1, 0> { enum { idx = 7 }; };
300 
301 }
302 
303 template<typename GridT>
304 class BoxStencil: public BaseStencil<BoxStencil<GridT>, 8, GridT>
305 {
306  using SelfT = BoxStencil<GridT>;
308 public:
309  using GridType = GridT;
310  using TreeType = typename GridT::TreeType;
311  using ValueType = typename GridT::ValueType;
312 
313  static constexpr int SIZE = 8;
314 
316 
317  /// Return linear offset for the specified stencil point relative to its center
318  template<int i, int j, int k>
319  __hostdev__ unsigned int pos() const { return BoxPt<i,j,k>::idx; }
320 
321  /// @brief Return true if the center of the stencil intersects the
322  /// iso-contour specified by the isoValue
323  __hostdev__ inline bool intersects(ValueType isoValue = ValueType(0)) const
324  {
325  const bool less = mValues[0] < isoValue;
326  return (less ^ (mValues[1] < isoValue)) ||
327  (less ^ (mValues[2] < isoValue)) ||
328  (less ^ (mValues[3] < isoValue)) ||
329  (less ^ (mValues[4] < isoValue)) ||
330  (less ^ (mValues[5] < isoValue)) ||
331  (less ^ (mValues[6] < isoValue)) ||
332  (less ^ (mValues[7] < isoValue)) ;
333  }
334 
335  /// @brief Return the trilinear interpolation at the normalized position.
336  /// @param xyz Floating point coordinate position. Index space and NOT world space.
337  /// @warning It is assumed that the stencil has already been moved
338  /// to the relevant voxel position, e.g. using moveTo(xyz).
339  /// @note Trilinear interpolation kernal reads as:
340  /// v000 (1-u)(1-v)(1-w) + v001 (1-u)(1-v)w + v010 (1-u)v(1-w) + v011 (1-u)vw
341  /// + v100 u(1-v)(1-w) + v101 u(1-v)w + v110 uv(1-w) + v111 uvw
343  {
344  const ValueType u = xyz[0] - mCenter[0];
345  const ValueType v = xyz[1] - mCenter[1];
346  const ValueType w = xyz[2] - mCenter[2];
347 
348  NANOVDB_ASSERT(u>=0 && u<=1);
349  NANOVDB_ASSERT(v>=0 && v<=1);
350  NANOVDB_ASSERT(w>=0 && w<=1);
351 
352  ValueType V = BaseType::template getValue<0,0,0>();
353  ValueType A = V + (BaseType::template getValue<0,0,1>() - V) * w;
354  V = BaseType::template getValue< 0, 1, 0>();
355  ValueType B = V + (BaseType::template getValue<0,1,1>() - V) * w;
356  ValueType C = A + (B - A) * v;
357 
358  V = BaseType::template getValue<1,0,0>();
359  A = V + (BaseType::template getValue<1,0,1>() - V) * w;
360  V = BaseType::template getValue<1,1,0>();
361  B = V + (BaseType::template getValue<1,1,1>() - V) * w;
362  ValueType D = A + (B - A) * v;
363 
364  return C + (D - C) * u;
365  }
366 
367  /// @brief Return the gradient in world space of the trilinear interpolation kernel.
368  /// @param xyz Floating point coordinate position.
369  /// @warning It is assumed that the stencil has already been moved
370  /// to the relevant voxel position, e.g. using moveTo(xyz).
371  /// @note Computed as partial derivatives of the trilinear interpolation kernel:
372  /// v000 (1-u)(1-v)(1-w) + v001 (1-u)(1-v)w + v010 (1-u)v(1-w) + v011 (1-u)vw
373  /// + v100 u(1-v)(1-w) + v101 u(1-v)w + v110 uv(1-w) + v111 uvw
375  {
376  const ValueType u = xyz[0] - mCenter[0];
377  const ValueType v = xyz[1] - mCenter[1];
378  const ValueType w = xyz[2] - mCenter[2];
379 
380  NANOVDB_ASSERT(u>=0 && u<=1);
381  NANOVDB_ASSERT(v>=0 && v<=1);
382  NANOVDB_ASSERT(w>=0 && w<=1);
383 
384  ValueType D[4]={BaseType::template getValue<0,0,1>()-BaseType::template getValue<0,0,0>(),
385  BaseType::template getValue<0,1,1>()-BaseType::template getValue<0,1,0>(),
386  BaseType::template getValue<1,0,1>()-BaseType::template getValue<1,0,0>(),
387  BaseType::template getValue<1,1,1>()-BaseType::template getValue<1,1,0>()};
388 
389  // Z component
390  ValueType A = D[0] + (D[1]- D[0]) * v;
391  ValueType B = D[2] + (D[3]- D[2]) * v;
392  Vec3<ValueType> grad(0, 0, A + (B - A) * u);
393 
394  D[0] = BaseType::template getValue<0,0,0>() + D[0] * w;
395  D[1] = BaseType::template getValue<0,1,0>() + D[1] * w;
396  D[2] = BaseType::template getValue<1,0,0>() + D[2] * w;
397  D[3] = BaseType::template getValue<1,1,0>() + D[3] * w;
398 
399  // X component
400  A = D[0] + (D[1] - D[0]) * v;
401  B = D[2] + (D[3] - D[2]) * v;
402 
403  grad[0] = B - A;
404 
405  // Y component
406  A = D[1] - D[0];
407  B = D[3] - D[2];
408 
409  grad[1] = A + (B - A) * u;
410 
411  return BaseType::mGrid->map().applyIJT(grad);
412  }
413 
414 private:
415  __hostdev__ inline void init(const Coord& ijk)
416  {
417  mValues[ 1] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
418  mValues[ 2] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
419  mValues[ 3] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
420  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
421  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
422  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 1, 1, 1));
423  mValues[ 7] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
424  }
425 
426  template<typename, int, typename> friend class BaseStencil; // allow base class to call init()
427  using BaseType::mAcc;
428  using BaseType::mValues;
429  using BaseType::mCenter;
430 };// BoxStencil class
431 
432 
433 // ---------------------------- GradStencil ----------------------------
434 
435 namespace { // anonymous namespace for stencil-layout map
436 
437  template<int i, int j, int k> struct GradPt {};
438  template<> struct GradPt< 0, 0, 0> { enum { idx = 0 }; };
439  template<> struct GradPt< 1, 0, 0> { enum { idx = 2 }; };
440  template<> struct GradPt< 0, 1, 0> { enum { idx = 4 }; };
441  template<> struct GradPt< 0, 0, 1> { enum { idx = 6 }; };
442  template<> struct GradPt<-1, 0, 0> { enum { idx = 1 }; };
443  template<> struct GradPt< 0,-1, 0> { enum { idx = 3 }; };
444  template<> struct GradPt< 0, 0,-1> { enum { idx = 5 }; };
445 }
446 
447 /// This is a simple 7-point nearest neighbor stencil that supports
448 /// gradient by second-order central differencing, first-order upwinding,
449 /// Laplacian, closest-point transform and zero-crossing test.
450 ///
451 /// @note For optimal random access performance this class
452 /// includes its own grid accessor.
453 template<typename GridT>
454 class GradStencil : public BaseStencil<GradStencil<GridT>, 7, GridT>
455 {
456  using SelfT = GradStencil<GridT>;
458 public:
459  using GridType = GridT;
460  using TreeType = typename GridT::TreeType;
461  using ValueType = typename GridT::ValueType;
462 
463  static constexpr int SIZE = 7;
464 
466  : BaseType(grid)
467  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
468  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
469  {
470  }
471 
473  : BaseType(grid)
474  , mInv2Dx(ValueType(0.5 / dx))
475  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
476  {
477  }
478 
479  /// @brief Return the norm square of the single-sided upwind gradient
480  /// (computed via Godunov's scheme) at the previously buffered location.
481  ///
482  /// @note This method should not be called until the stencil
483  /// buffer has been populated via a call to moveTo(ijk).
485  {
486  return mInvDx2 * GodunovsNormSqrd(mValues[0] > ValueType(0),
487  mValues[0] - mValues[1],
488  mValues[2] - mValues[0],
489  mValues[0] - mValues[3],
490  mValues[4] - mValues[0],
491  mValues[0] - mValues[5],
492  mValues[6] - mValues[0]);
493  }
494 
495  /// @brief Return the gradient computed at the previously buffered
496  /// location by second order central differencing.
497  ///
498  /// @note This method should not be called until the stencil
499  /// buffer has been populated via a call to moveTo(ijk).
501  {
502  return Vec3<ValueType>(mValues[2] - mValues[1],
503  mValues[4] - mValues[3],
504  mValues[6] - mValues[5])*mInv2Dx;
505  }
506  /// @brief Return the first-order upwind gradient corresponding to the direction V.
507  ///
508  /// @note This method should not be called until the stencil
509  /// buffer has been populated via a call to moveTo(ijk).
511  {
512  return Vec3<ValueType>(
513  V[0]>0 ? mValues[0] - mValues[1] : mValues[2] - mValues[0],
514  V[1]>0 ? mValues[0] - mValues[3] : mValues[4] - mValues[0],
515  V[2]>0 ? mValues[0] - mValues[5] : mValues[6] - mValues[0])*2*mInv2Dx;
516  }
517 
518  /// Return the Laplacian computed at the previously buffered
519  /// location by second-order central differencing.
521  {
522  return mInvDx2 * (mValues[1] + mValues[2] +
523  mValues[3] + mValues[4] +
524  mValues[5] + mValues[6] - 6*mValues[0]);
525  }
526 
527  /// Return @c true if the sign of the value at the center point of the stencil
528  /// is different from the signs of any of its six nearest neighbors.
529  __hostdev__ inline bool zeroCrossing() const
530  {
531  return (mValues[0]>0 ? (mValues[1]<0 || mValues[2]<0 || mValues[3]<0 || mValues[4]<0 || mValues[5]<0 || mValues[6]<0)
532  : (mValues[1]>0 || mValues[2]>0 || mValues[3]>0 || mValues[4]>0 || mValues[5]>0 || mValues[6]>0));
533  }
534 
535  /// @brief Compute the closest-point transform to a level set.
536  /// @return the closest point in index space to the surface
537  /// from which the level set was derived.
538  ///
539  /// @note This method assumes that the grid represents a level set
540  /// with distances in world units and a simple affine transfrom
541  /// with uniform scaling.
543  {
544  const Coord& ijk = BaseType::getCenterCoord();
545  const ValueType d = ValueType(mValues[0] * 0.5 * mInvDx2); // distance in voxels / (2dx^2)
546  const auto value = Vec3<ValueType>(ijk[0] - d*(mValues[2] - mValues[1]),
547  ijk[1] - d*(mValues[4] - mValues[3]),
548  ijk[2] - d*(mValues[6] - mValues[5]));
549  return value;
550  }
551 
552  /// Return linear offset for the specified stencil point relative to its center
553  template<int i, int j, int k>
554  __hostdev__ unsigned int pos() const { return GradPt<i,j,k>::idx; }
555 
556 private:
557 
558  __hostdev__ inline void init(const Coord& ijk)
559  {
560  mValues[ 1] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
561  mValues[ 2] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
562 
563  mValues[ 3] = mAcc.getValue(ijk.offsetBy( 0,-1, 0));
564  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
565 
566  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 0, 0,-1));
567  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
568  }
569 
570  template<typename, int, typename> friend class BaseStencil; // allow base class to call init()
571  using BaseType::mAcc;
572  using BaseType::mValues;
573  const ValueType mInv2Dx, mInvDx2;
574 }; // GradStencil class
575 
576 
577 // ---------------------------- WenoStencil ----------------------------
578 
579 namespace { // anonymous namespace for stencil-layout map
580 
581  template<int i, int j, int k> struct WenoPt {};
582  template<> struct WenoPt< 0, 0, 0> { enum { idx = 0 }; };
583 
584  template<> struct WenoPt<-3, 0, 0> { enum { idx = 1 }; };
585  template<> struct WenoPt<-2, 0, 0> { enum { idx = 2 }; };
586  template<> struct WenoPt<-1, 0, 0> { enum { idx = 3 }; };
587  template<> struct WenoPt< 1, 0, 0> { enum { idx = 4 }; };
588  template<> struct WenoPt< 2, 0, 0> { enum { idx = 5 }; };
589  template<> struct WenoPt< 3, 0, 0> { enum { idx = 6 }; };
590 
591  template<> struct WenoPt< 0,-3, 0> { enum { idx = 7 }; };
592  template<> struct WenoPt< 0,-2, 0> { enum { idx = 8 }; };
593  template<> struct WenoPt< 0,-1, 0> { enum { idx = 9 }; };
594  template<> struct WenoPt< 0, 1, 0> { enum { idx =10 }; };
595  template<> struct WenoPt< 0, 2, 0> { enum { idx =11 }; };
596  template<> struct WenoPt< 0, 3, 0> { enum { idx =12 }; };
597 
598  template<> struct WenoPt< 0, 0,-3> { enum { idx =13 }; };
599  template<> struct WenoPt< 0, 0,-2> { enum { idx =14 }; };
600  template<> struct WenoPt< 0, 0,-1> { enum { idx =15 }; };
601  template<> struct WenoPt< 0, 0, 1> { enum { idx =16 }; };
602  template<> struct WenoPt< 0, 0, 2> { enum { idx =17 }; };
603  template<> struct WenoPt< 0, 0, 3> { enum { idx =18 }; };
604 
605 }
606 
607 /// @brief This is a special 19-point stencil that supports optimal fifth-order WENO
608 /// upwinding, second-order central differencing, Laplacian, and zero-crossing test.
609 ///
610 /// @note For optimal random access performance this class
611 /// includes its own grid accessor.
612 template<typename GridT, typename RealT = typename GridT::ValueType>
613 class WenoStencil: public BaseStencil<WenoStencil<GridT>, 19, GridT>
614 {
615  using SelfT = WenoStencil<GridT>;
617 public:
618  using GridType = GridT;
619  using TreeType = typename GridT::TreeType;
620  using ValueType = typename GridT::ValueType;
621 
622  static constexpr int SIZE = 19;
623 
625  : BaseType(grid)
626  , mDx2(ValueType(Pow2(grid.voxelSize()[0])))
627  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
628  , mInvDx2(ValueType(1.0 / mDx2))
629  {
630  }
631 
633  : BaseType(grid)
634  , mDx2(ValueType(dx * dx))
635  , mInv2Dx(ValueType(0.5 / dx))
636  , mInvDx2(ValueType(1.0 / mDx2))
637  {
638  }
639 
640  /// @brief Return the norm-square of the WENO upwind gradient (computed via
641  /// WENO upwinding and Godunov's scheme) at the previously buffered location.
642  ///
643  /// @note This method should not be called until the stencil
644  /// buffer has been populated via a call to moveTo(ijk).
645  __hostdev__ inline ValueType normSqGrad(ValueType isoValue = ValueType(0)) const
646  {
647  const ValueType* v = mValues;
648  const RealT
649  dP_xm = WENO5<RealT>(v[ 2]-v[ 1],v[ 3]-v[ 2],v[ 0]-v[ 3],v[ 4]-v[ 0],v[ 5]-v[ 4],mDx2),
650  dP_xp = WENO5<RealT>(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0],v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
651  dP_ym = WENO5<RealT>(v[ 8]-v[ 7],v[ 9]-v[ 8],v[ 0]-v[ 9],v[10]-v[ 0],v[11]-v[10],mDx2),
652  dP_yp = WENO5<RealT>(v[12]-v[11],v[11]-v[10],v[10]-v[ 0],v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
653  dP_zm = WENO5<RealT>(v[14]-v[13],v[15]-v[14],v[ 0]-v[15],v[16]-v[ 0],v[17]-v[16],mDx2),
654  dP_zp = WENO5<RealT>(v[18]-v[17],v[17]-v[16],v[16]-v[ 0],v[ 0]-v[15],v[15]-v[14],mDx2);
655  return mInvDx2*static_cast<ValueType>(
656  GodunovsNormSqrd(v[0]>isoValue, dP_xm, dP_xp, dP_ym, dP_yp, dP_zm, dP_zp));
657  }
658 
659  /// Return the optimal fifth-order upwind gradient corresponding to the
660  /// direction V.
661  ///
662  /// @note This method should not be called until the stencil
663  /// buffer has been populated via a call to moveTo(ijk).
665  {
666  const ValueType* v = mValues;
667  return 2*mInv2Dx * Vec3<ValueType>(
668  V[0]>0 ? WENO5<RealT>(v[ 2]-v[ 1],v[ 3]-v[ 2],v[ 0]-v[ 3], v[ 4]-v[ 0],v[ 5]-v[ 4],mDx2)
669  : WENO5<RealT>(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0], v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
670  V[1]>0 ? WENO5<RealT>(v[ 8]-v[ 7],v[ 9]-v[ 8],v[ 0]-v[ 9], v[10]-v[ 0],v[11]-v[10],mDx2)
671  : WENO5<RealT>(v[12]-v[11],v[11]-v[10],v[10]-v[ 0], v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
672  V[2]>0 ? WENO5<RealT>(v[14]-v[13],v[15]-v[14],v[ 0]-v[15], v[16]-v[ 0],v[17]-v[16],mDx2)
673  : WENO5<RealT>(v[18]-v[17],v[17]-v[16],v[16]-v[ 0], v[ 0]-v[15],v[15]-v[14],mDx2));
674  }
675  /// Return the gradient computed at the previously buffered
676  /// location by second-order central differencing.
677  ///
678  /// @note This method should not be called until the stencil
679  /// buffer has been populated via a call to moveTo(ijk).
681  {
682  return mInv2Dx * Vec3<ValueType>(mValues[ 4] - mValues[ 3],
683  mValues[10] - mValues[ 9],
684  mValues[16] - mValues[15]);
685  }
686 
687  /// Return the Laplacian computed at the previously buffered
688  /// location by second-order central differencing.
689  ///
690  /// @note This method should not be called until the stencil
691  /// buffer has been populated via a call to moveTo(ijk).
693  {
694  return mInvDx2 * (
695  mValues[ 3] + mValues[ 4] +
696  mValues[ 9] + mValues[10] +
697  mValues[15] + mValues[16] - 6*mValues[0]);
698  }
699 
700  /// Return @c true if the sign of the value at the center point of the stencil
701  /// differs from the sign of any of its six nearest neighbors
702  __hostdev__ inline bool zeroCrossing() const
703  {
704  const ValueType* v = mValues;
705  return (v[ 0]>0 ? (v[ 3]<0 || v[ 4]<0 || v[ 9]<0 || v[10]<0 || v[15]<0 || v[16]<0)
706  : (v[ 3]>0 || v[ 4]>0 || v[ 9]>0 || v[10]>0 || v[15]>0 || v[16]>0));
707  }
708 
709  /// Return linear offset for the specified stencil point relative to its center
710  template<int i, int j, int k>
711  __hostdev__ unsigned int pos() const { return WenoPt<i,j,k>::idx; }
712 
713 private:
714  __hostdev__ inline void init(const Coord& ijk)
715  {
716  mValues[ 1] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
717  mValues[ 2] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
718  mValues[ 3] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
719  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
720  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
721  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
722 
723  mValues[ 7] = mAcc.getValue(ijk.offsetBy( 0, -3, 0));
724  mValues[ 8] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
725  mValues[ 9] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
726  mValues[10] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
727  mValues[11] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
728  mValues[12] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
729 
730  mValues[13] = mAcc.getValue(ijk.offsetBy( 0, 0, -3));
731  mValues[14] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
732  mValues[15] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
733  mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
734  mValues[17] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
735  mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
736  }
737 
738  template<typename, int, typename> friend class BaseStencil; // allow base class to call init()
739  using BaseType::mAcc;
740  using BaseType::mValues;
741  const ValueType mDx2, mInv2Dx, mInvDx2;
742 }; // WenoStencil class
743 
744 
745 // ---------------------------- CurvatureStencil ----------------------------
746 
747 namespace { // anonymous namespace for stencil-layout map
748 
749  template<int i, int j, int k> struct CurvPt {};
750  template<> struct CurvPt< 0, 0, 0> { enum { idx = 0 }; };
751 
752  template<> struct CurvPt<-1, 0, 0> { enum { idx = 1 }; };
753  template<> struct CurvPt< 1, 0, 0> { enum { idx = 2 }; };
754 
755  template<> struct CurvPt< 0,-1, 0> { enum { idx = 3 }; };
756  template<> struct CurvPt< 0, 1, 0> { enum { idx = 4 }; };
757 
758  template<> struct CurvPt< 0, 0,-1> { enum { idx = 5 }; };
759  template<> struct CurvPt< 0, 0, 1> { enum { idx = 6 }; };
760 
761  template<> struct CurvPt<-1,-1, 0> { enum { idx = 7 }; };
762  template<> struct CurvPt< 1,-1, 0> { enum { idx = 8 }; };
763  template<> struct CurvPt<-1, 1, 0> { enum { idx = 9 }; };
764  template<> struct CurvPt< 1, 1, 0> { enum { idx =10 }; };
765 
766  template<> struct CurvPt<-1, 0,-1> { enum { idx =11 }; };
767  template<> struct CurvPt< 1, 0,-1> { enum { idx =12 }; };
768  template<> struct CurvPt<-1, 0, 1> { enum { idx =13 }; };
769  template<> struct CurvPt< 1, 0, 1> { enum { idx =14 }; };
770 
771  template<> struct CurvPt< 0,-1,-1> { enum { idx =15 }; };
772  template<> struct CurvPt< 0, 1,-1> { enum { idx =16 }; };
773  template<> struct CurvPt< 0,-1, 1> { enum { idx =17 }; };
774  template<> struct CurvPt< 0, 1, 1> { enum { idx =18 }; };
775 
776 }
777 
778 template<typename GridT, typename RealT = typename GridT::ValueType>
779 class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT>, 19, GridT>
780 {
783 public:
784  using GridType = GridT;
785  using TreeType = typename GridT::TreeType;
786  using ValueType = typename GridT::ValueType;
787 
788  static constexpr int SIZE = 19;
789 
791  : BaseType(grid)
792  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
793  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
794  {
795  }
796 
798  : BaseType(grid)
799  , mInv2Dx(ValueType(0.5 / dx))
800  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
801  {
802  }
803 
804  /// @brief Return the mean curvature at the previously buffered location.
805  ///
806  /// @note This method should not be called until the stencil
807  /// buffer has been populated via a call to moveTo(ijk).
809  {
810  RealT alpha, normGrad;
811  return this->meanCurvature(alpha, normGrad) ?
812  ValueType(alpha*mInv2Dx/Pow3(normGrad)) : 0;
813  }
814 
815  /// @brief Return the Gaussian curvature at the previously buffered location.
816  ///
817  /// @note This method should not be called until the stencil
818  /// buffer has been populated via a call to moveTo(ijk).
820  {
821  RealT alpha, normGrad;
822  return this->gaussianCurvature(alpha, normGrad) ?
823  ValueType(alpha*mInvDx2/Pow4(normGrad)) : 0;
824  }
825 
826  /// @brief Return both the mean and the Gaussian curvature at the
827  /// previously buffered location.
828  ///
829  /// @note This method should not be called until the stencil
830  /// buffer has been populated via a call to moveTo(ijk).
831  __hostdev__ inline void curvatures(ValueType &mean, ValueType& gauss) const
832  {
833  RealT alphaM, alphaG, normGrad;
834  if (this->curvatures(alphaM, alphaG, normGrad)) {
835  mean = ValueType(alphaM*mInv2Dx/Pow3(normGrad));
836  gauss = ValueType(alphaG*mInvDx2/Pow4(normGrad));
837  } else {
838  mean = gauss = 0;
839  }
840  }
841 
842  /// Return the mean curvature multiplied by the norm of the
843  /// central-difference gradient. This method is very useful for
844  /// mean-curvature flow of level sets!
845  ///
846  /// @note This method should not be called until the stencil
847  /// buffer has been populated via a call to moveTo(ijk).
849  {
850  RealT alpha, normGrad;
851  return this->meanCurvature(alpha, normGrad) ?
852  ValueType(alpha*mInvDx2/(2*Pow2(normGrad))) : 0;
853  }
854 
855  /// Return the mean Gaussian multiplied by the norm of the
856  /// central-difference gradient.
857  ///
858  /// @note This method should not be called until the stencil
859  /// buffer has been populated via a call to moveTo(ijk).
861  {
862  RealT alpha, normGrad;
863  return this->gaussianCurvature(alpha, normGrad) ?
864  ValueType(2*alpha*mInv2Dx*mInvDx2/Pow3(normGrad)) : 0;
865  }
866 
867  /// @brief Return both the mean and the Gaussian curvature at the
868  /// previously buffered location.
869  ///
870  /// @note This method should not be called until the stencil
871  /// buffer has been populated via a call to moveTo(ijk).
873  {
874  RealT alphaM, alphaG, normGrad;
875  if (this->curvatures(alphaM, alphaG, normGrad)) {
876  mean = ValueType(alphaM*mInvDx2/(2*Pow2(normGrad)));
877  gauss = ValueType(2*alphaG*mInv2Dx*mInvDx2/Pow3(normGrad));
878  } else {
879  mean = gauss = 0;
880  }
881  }
882 
883  /// @brief Computes the minimum and maximum principal curvature at the
884  /// previously buffered location.
885  ///
886  /// @note This method should not be called until the stencil
887  /// buffer has been populated via a call to moveTo(ijk).
889  {
890  min = max = 0;
891  RealT alphaM, alphaG, normGrad;
892  if (this->curvatures(alphaM, alphaG, normGrad)) {
893  const RealT mean = alphaM*mInv2Dx/Pow3(normGrad);
894  const RealT tmp = Sqrt(mean*mean - alphaG*mInvDx2/Pow4(normGrad));
895  min = ValueType(mean - tmp);
896  max = ValueType(mean + tmp);
897  }
898  }
899 
900  /// Return the Laplacian computed at the previously buffered
901  /// location by second-order central differencing.
902  ///
903  /// @note This method should not be called until the stencil
904  /// buffer has been populated via a call to moveTo(ijk).
906  {
907  return mInvDx2 * (
908  mValues[1] + mValues[2] +
909  mValues[3] + mValues[4] +
910  mValues[5] + mValues[6] - 6*mValues[0]);
911  }
912 
913  /// Return the gradient computed at the previously buffered
914  /// location by second-order central differencing.
915  ///
916  /// @note This method should not be called until the stencil
917  /// buffer has been populated via a call to moveTo(ijk).
919  {
920  return Vec3<ValueType>(
921  mValues[2] - mValues[1],
922  mValues[4] - mValues[3],
923  mValues[6] - mValues[5])*mInv2Dx;
924  }
925 
926  /// Return linear offset for the specified stencil point relative to its center
927  template<int i, int j, int k>
928  __hostdev__ unsigned int pos() const { return CurvPt<i,j,k>::idx; }
929 
930 private:
931  __hostdev__ inline void init(const Coord &ijk)
932  {
933  mValues[ 1] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
934  mValues[ 2] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
935 
936  mValues[ 3] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
937  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
938 
939  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
940  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
941 
942  mValues[ 7] = mAcc.getValue(ijk.offsetBy(-1, -1, 0));
943  mValues[ 8] = mAcc.getValue(ijk.offsetBy( 1, -1, 0));
944  mValues[ 9] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
945  mValues[10] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
946 
947  mValues[11] = mAcc.getValue(ijk.offsetBy(-1, 0, -1));
948  mValues[12] = mAcc.getValue(ijk.offsetBy( 1, 0, -1));
949  mValues[13] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
950  mValues[14] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
951 
952  mValues[15] = mAcc.getValue(ijk.offsetBy( 0, -1, -1));
953  mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 1, -1));
954  mValues[17] = mAcc.getValue(ijk.offsetBy( 0, -1, 1));
955  mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
956  }
957 
958  __hostdev__ inline RealT Dx() const { return 0.5*(mValues[2] - mValues[1]); }// * 1/dx
959  __hostdev__ inline RealT Dy() const { return 0.5*(mValues[4] - mValues[3]); }// * 1/dx
960  __hostdev__ inline RealT Dz() const { return 0.5*(mValues[6] - mValues[5]); }// * 1/dx
961  __hostdev__ inline RealT Dxx() const { return mValues[2] - 2 * mValues[0] + mValues[1]; }// * 1/dx2
962  __hostdev__ inline RealT Dyy() const { return mValues[4] - 2 * mValues[0] + mValues[3]; }// * 1/dx2}
963  __hostdev__ inline RealT Dzz() const { return mValues[6] - 2 * mValues[0] + mValues[5]; }// * 1/dx2
964  __hostdev__ inline RealT Dxy() const { return 0.25 * (mValues[10] - mValues[ 8] + mValues[ 7] - mValues[ 9]); }// * 1/dx2
965  __hostdev__ inline RealT Dxz() const { return 0.25 * (mValues[14] - mValues[12] + mValues[11] - mValues[13]); }// * 1/dx2
966  __hostdev__ inline RealT Dyz() const { return 0.25 * (mValues[18] - mValues[16] + mValues[15] - mValues[17]); }// * 1/dx2
967 
968  __hostdev__ inline bool meanCurvature(RealT& alpha, RealT& normGrad) const
969  {
970  // For performance all finite differences are unscaled wrt dx
971  const RealT Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
972  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
973  if (normGrad2 <= Tolerance<RealT>::value()) {
974  alpha = normGrad = 0;
975  return false;
976  }
977  const RealT Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz();
978  alpha = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
979  2*(Dx*(Dy*this->Dxy() + Dz*this->Dxz()) + Dy*Dz*this->Dyz());// * 1/dx^4
980  normGrad = Sqrt(normGrad2); // * 1/dx
981  return true;
982  }
983 
984  __hostdev__ inline bool gaussianCurvature(RealT& alpha, RealT& normGrad) const
985  {
986  // For performance all finite differences are unscaled wrt dx
987  const RealT Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
988  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
989  if (normGrad2 <= Tolerance<RealT>::value()) {
990  alpha = normGrad = 0;
991  return false;
992  }
993  const RealT Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
994  Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
995  alpha = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
996  2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// * 1/dx^6
997  normGrad = Sqrt(normGrad2); // * 1/dx
998  return true;
999  }
1000 
1001  __hostdev__ inline bool curvatures(RealT& alphaM, RealT& alphaG, RealT& normGrad) const
1002  {
1003  // For performance all finite differences are unscaled wrt dx
1004  const RealT Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1005  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1006  if (normGrad2 <= Tolerance<RealT>::value()) {
1007  alphaM = alphaG =normGrad = 0;
1008  return false;
1009  }
1010  const RealT Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
1011  Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
1012  alphaM = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
1013  2*(Dx*(Dy*Dxy + Dz*Dxz) + Dy*Dz*Dyz);// *1/dx^4
1014  alphaG = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
1015  2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// *1/dx^6
1016  normGrad = Sqrt(normGrad2); // * 1/dx
1017  return true;
1018  }
1019 
1020  template<typename, int, typename> friend class BaseStencil; // allow base class to call init()
1021  using BaseType::mAcc;
1022  using BaseType::mValues;
1023  const ValueType mInv2Dx, mInvDx2;
1024 }; // CurvatureStencil class
1025 
1026 } // end nanovdb namespace
1027 
1028 #endif // NANOVDB_STENCILS_HAS_BEEN_INCLUDED
ValueType gaussianCurvature() const
Return the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:819
bool intersects(ValueType isoValue=ValueType(0)) const
Return true if the center of the stencil intersects the.
Definition: Stencils.h:323
typename GridT::TreeType TreeType
Definition: Stencils.h:619
typename GridT::TreeType TreeType
Definition: Stencils.h:310
Definition: Stencils.h:779
Vec3< ValueType > gradient() const
Definition: Stencils.h:918
WenoStencil(const GridType &grid)
Definition: Stencils.h:624
ValueType meanCurvatureNormGrad() const
Definition: Stencils.h:848
bool zeroCrossing() const
Definition: Stencils.h:529
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:266
bool none() const
Definition: Stencils.h:233
WenoStencil(const GridType &grid, double dx)
Definition: Stencils.h:632
Mask intersectionMask(ValueType isoValue=ValueType(0)) const
Return true a bit-mask where the 6 lower bits indicates if the center of the stencil intersects the i...
Definition: Stencils.h:247
GridT GridType
Definition: Stencils.h:309
typename GridT::ValueType ValueType
Definition: Stencils.h:311
ValueType laplacian() const
Definition: Stencils.h:692
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:554
Coord mCenter
Definition: Stencils.h:280
bool intersects(const ValueType &isoValue=ValueType(0)) const
Return true if the center of the stencil intersects the iso-contour specified by the isoValue...
Definition: Stencils.h:216
CoordT RoundDown(const Vec3T< RealT > &xyz)
Definition: NanoVDB.h:788
GridT GridType
Definition: Stencils.h:784
RealT GodunovsNormSqrd(bool isOutside, RealT dP_xm, RealT dP_xp, RealT dP_ym, RealT dP_yp, RealT dP_zm, RealT dP_zp)
Definition: Stencils.h:62
ValueType laplacian() const
Definition: Stencils.h:905
typename GridT::TreeType TreeType
Definition: Stencils.h:101
GradStencil(const GridType &grid, double dx)
Definition: Stencils.h:472
GradStencil(const GridType &grid)
Definition: Stencils.h:465
CurvatureStencil(const GridType &grid)
Definition: Stencils.h:790
ValueType normSqGrad(ValueType isoValue=ValueType(0)) const
Return the norm-square of the WENO upwind gradient (computed via WENO upwinding and Godunov&#39;s scheme)...
Definition: Stencils.h:645
Type Max(Type a, Type b)
Definition: NanoVDB.h:672
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1034
Type Min(Type a, Type b)
Definition: NanoVDB.h:651
ValueType max() const
Return the largest value in the stencil buffer.
Definition: Stencils.h:199
typename GridT::ValueType ValueType
Definition: Stencils.h:786
void curvatures(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:831
void moveTo(const Coord &ijk, const ValueType &centerValue)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors. The method also takes a value of the center element of the stencil, assuming it is already known.
Definition: Stencils.h:119
Mask()
Definition: Stencils.h:228
ValueType min() const
Return the smallest value in the stencil buffer.
Definition: Stencils.h:189
This is a special 19-point stencil that supports optimal fifth-order WENO upwinding, second-order central differencing, Laplacian, and zero-crossing test.
Definition: Stencils.h:613
void setValue(const ValueType &value)
Set the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:172
Definition: Stencils.h:454
ValueType laplacian() const
Definition: Stencils.h:520
typename GridT::ValueType ValueType
Definition: Stencils.h:461
ValueType gaussianCurvatureNormGrad() const
Definition: Stencils.h:860
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:928
Definition: NanoVDB.h:184
const Coord & getCenterCoord() const
Return the coordinates of the center point of the stencil.
Definition: Stencils.h:209
Bit-mask to encode active states and facilitate sequential iterators and a fast codec for I/O compres...
Definition: NanoVDB.h:1795
T Pow4(T x)
Definition: NanoVDB.h:742
Definition: Stencils.h:304
A simple vector class with three double components, similar to openvdb::math::Vec3.
Definition: NanoVDB.h:856
Vec3< ValueType > gradient(const Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:374
GridT GridType
Definition: Stencils.h:459
AccessorType mAcc
Definition: Stencils.h:278
CurvatureStencil(const GridType &grid, double dx)
Definition: Stencils.h:797
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, RealT scale2=1.0)
Implementation of nominally fifth-order finite-difference WENO.
Definition: Stencils.h:35
void curvaturesNormGrad(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:872
void moveTo(const Vec3< RealType > &xyz)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:146
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:711
Vec3< ValueType > gradient(const Vec3< ValueType > &V) const
Return the first-order upwind gradient corresponding to the direction V.
Definition: Stencils.h:510
BaseStencil(const GridType &grid)
Definition: Stencils.h:270
static int size()
Return the size of the stencil buffer.
Definition: Stencils.h:178
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:262
typename GridT::TreeType TreeType
Definition: Stencils.h:785
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:157
uint32_t CountOn(uint64_t v)
Definition: NanoVDB.h:1774
Definition: Stencils.h:226
typename GridT::ValueType ValueType
Definition: Stencils.h:620
Vec3< ValueType > gradient() const
Return the gradient computed at the previously buffered location by second order central differencing...
Definition: Stencils.h:500
Vec3< ValueType > cpt()
Compute the closest-point transform to a level set.
Definition: Stencils.h:542
ValueT value
Definition: GridBuilder.h:1287
Coord offsetBy(ValueType dx, ValueType dy, ValueType dz) const
Definition: NanoVDB.h:1002
bool any() const
Definition: Stencils.h:231
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:132
ValueType mean() const
Return the mean value of the current stencil.
Definition: Stencils.h:181
typename GridT::TreeType TreeType
Definition: Stencils.h:460
int count() const
Definition: Stencils.h:234
void set(uint32_t n, bool On)
Definition: NanoVDB.h:1924
#define NANOVDB_ASSERT(x)
Definition: NanoVDB.h:149
ValueType meanCurvature() const
Return the mean curvature at the previously buffered location.
Definition: Stencils.h:808
bool zeroCrossing() const
Definition: Stencils.h:702
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:216
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: NanoVDB.h:795
const GridType * mGrid
Definition: Stencils.h:277
bool test(int i) const
Definition: Stencils.h:230
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:107
Definition: Stencils.h:96
ValueType mValues[SIZE]
Definition: Stencils.h:279
bool all() const
Definition: Stencils.h:232
Tolerance for floating-point comparison.
Definition: NanoVDB.h:581
T Pow2(T x)
Definition: NanoVDB.h:730
BoxStencil(const GridType &grid)
Definition: Stencils.h:315
typename GridT::ValueType ValueType
Definition: Stencils.h:99
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:319
Vec3< ValueType > gradient() const
Definition: Stencils.h:680
const ValueType & getCenterValue() const
Return the value at the center of the stencil.
Definition: Stencils.h:212
#define __hostdev__
Definition: NanoVDB.h:168
ValueType interpolation(const Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:342
T Pow3(T x)
Definition: NanoVDB.h:736
Signed (i, j, k) 32-bit integer coordinate class, similar to openvdb::math::Coord.
Definition: NanoVDB.h:859
GridT GridType
Definition: Stencils.h:618
void principalCurvatures(ValueType &min, ValueType &max) const
Computes the minimum and maximum principal curvature at the previously buffered location.
Definition: Stencils.h:888
const ValueType & getValue() const
Return the value at the specified location relative to the center of the stencil. ...
Definition: Stencils.h:165
ValueType normSqGrad() const
Return the norm square of the single-sided upwind gradient (computed via Godunov&#39;s scheme) at the pre...
Definition: Stencils.h:484
Vec3< ValueType > gradient(const Vec3< ValueType > &V) const
Definition: Stencils.h:664
uint8_t bits
Definition: Stencils.h:227