OpenVDB 10.0.1
Loading...
Searching...
No Matches
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
20namespace 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)
33template<typename ValueType, typename RealT = ValueType>
34__hostdev__ inline ValueType
35WENO5(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
60template <typename RealT>
61__hostdev__ inline RealT
62GodunovsNormSqrd(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
80template<typename RealT>
81__hostdev__ inline RealT
82GodunovsNormSqrd(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)
95template<typename DerivedType, int SIZE, typename GridT>
97{
98public:
99 using ValueType = typename GridT::ValueType;
100 using GridType = GridT;
101 using TreeType = typename GridT::TreeType;
102 using AccessorType = typename GridT::AccessorType;// ReadAccessor<ValueType>;
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>
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.
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;
229 __hostdev__ void set(int i) { bits |= (1 << i); }
230 __hostdev__ bool test(int i) const { return bits & (1 << i); }
231 __hostdev__ bool any() const { return bits > 0u; }
232 __hostdev__ bool all() const { return bits == 255u; }
233 __hostdev__ bool none() const { return bits == 0u; }
234 __hostdev__ 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).
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
268protected:
269 // Constructor is protected to prevent direct instantiation.
271 : mGrid(&grid)
272 , mAcc(grid)
273 , mCenter(Coord::max())
274 {
275 }
276
281
282}; // BaseStencil class
283
284
285// ---------------------------- BoxStencil ----------------------------
286
287
288namespace { // 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
303template<typename GridT>
304class BoxStencil: public BaseStencil<BoxStencil<GridT>, 8, GridT>
305{
306 using SelfT = BoxStencil<GridT>;
308public:
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
414private:
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
435namespace { // 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.
453template<typename GridT>
454class GradStencil : public BaseStencil<GradStencil<GridT>, 7, GridT>
455{
458public:
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
556private:
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
579namespace { // 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.
612template<typename GridT, typename RealT = typename GridT::ValueType>
613class WenoStencil: public BaseStencil<WenoStencil<GridT>, 19, GridT>
614{
617public:
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).
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
713private:
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
747namespace { // 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
778template<typename GridT, typename RealT = typename GridT::ValueType>
779class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT>, 19, GridT>
780{
783public:
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
930private:
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
ValueT value
Definition: GridBuilder.h:1290
#define __hostdev__
Definition: NanoVDB.h:192
#define NANOVDB_ASSERT(x)
Definition: NanoVDB.h:173
Definition: Stencils.h:97
ValueType min() const
Return the smallest value in the stencil buffer.
Definition: Stencils.h:189
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
ValueType mValues[SIZE]
Definition: Stencils.h:279
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:266
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:132
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:157
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
ValueType mean() const
Return the mean value of the current stencil.
Definition: Stencils.h:181
void setValue(const ValueType &value)
Set the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:172
void moveTo(const Coord &ijk, const ValueType &centerValue)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors....
Definition: Stencils.h:119
typename GridT::AccessorType AccessorType
Definition: Stencils.h:102
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:107
Coord mCenter
Definition: Stencils.h:280
const GridType * mGrid
Definition: Stencils.h:277
const ValueType & getCenterValue() const
Return the value at the center of the stencil.
Definition: Stencils.h:212
const ValueType & getValue() const
Return the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:165
const Coord & getCenterCoord() const
Return the coordinates of the center point of the stencil.
Definition: Stencils.h:209
typename GridT::TreeType TreeType
Definition: Stencils.h:101
typename GridT::ValueType ValueType
Definition: Stencils.h:99
BaseStencil(const GridType &grid)
Definition: Stencils.h:270
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:262
static int size()
Return the size of the stencil buffer.
Definition: Stencils.h:178
ValueType max() const
Return the largest value in the stencil buffer.
Definition: Stencils.h:199
AccessorType mAcc
Definition: Stencils.h:278
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:100
Definition: Stencils.h:305
static constexpr int SIZE
Definition: Stencils.h:313
Vec3< ValueType > gradient(const Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:374
ValueType interpolation(const Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:342
bool intersects(ValueType isoValue=ValueType(0)) const
Return true if the center of the stencil intersects the.
Definition: Stencils.h:323
BoxStencil(const GridType &grid)
Definition: Stencils.h:315
typename GridT::TreeType TreeType
Definition: Stencils.h:310
typename GridT::ValueType ValueType
Definition: Stencils.h:311
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:319
GridT GridType
Definition: Stencils.h:309
Signed (i, j, k) 32-bit integer coordinate class, similar to openvdb::math::Coord.
Definition: NanoVDB.h:967
Coord offsetBy(ValueType dx, ValueType dy, ValueType dz) const
Definition: NanoVDB.h:1116
Definition: Stencils.h:780
void curvaturesNormGrad(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:872
static constexpr int SIZE
Definition: Stencils.h:788
CurvatureStencil(const GridType &grid)
Definition: Stencils.h:790
ValueType meanCurvature() const
Return the mean curvature at the previously buffered location.
Definition: Stencils.h:808
Vec3< ValueType > gradient() const
Definition: Stencils.h:918
ValueType gaussianCurvatureNormGrad() const
Definition: Stencils.h:860
void principalCurvatures(ValueType &min, ValueType &max) const
Computes the minimum and maximum principal curvature at the previously buffered location.
Definition: Stencils.h:888
void curvatures(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:831
ValueType laplacian() const
Definition: Stencils.h:905
typename GridT::TreeType TreeType
Definition: Stencils.h:785
typename GridT::ValueType ValueType
Definition: Stencils.h:786
ValueType meanCurvatureNormGrad() const
Definition: Stencils.h:848
ValueType gaussianCurvature() const
Return the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:819
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:928
GridT GridType
Definition: Stencils.h:784
CurvatureStencil(const GridType &grid, double dx)
Definition: Stencils.h:797
Definition: Stencils.h:455
static constexpr int SIZE
Definition: Stencils.h:463
GradStencil(const GridType &grid, double dx)
Definition: Stencils.h:472
GradStencil(const GridType &grid)
Definition: Stencils.h:465
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
ValueType laplacian() const
Definition: Stencils.h:520
Vec3< ValueType > gradient(const Vec3< ValueType > &V) const
Return the first-order upwind gradient corresponding to the direction V.
Definition: Stencils.h:510
bool zeroCrossing() const
Definition: Stencils.h:529
typename GridT::TreeType TreeType
Definition: Stencils.h:460
typename GridT::ValueType ValueType
Definition: Stencils.h:461
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:554
ValueType normSqGrad() const
Return the norm square of the single-sided upwind gradient (computed via Godunov's scheme) at the pre...
Definition: Stencils.h:484
GridT GridType
Definition: Stencils.h:459
A simple vector class with three double components, similar to openvdb::math::Vec3.
Definition: NanoVDB.h:1158
This is a special 19-point stencil that supports optimal fifth-order WENO upwinding,...
Definition: Stencils.h:614
static constexpr int SIZE
Definition: Stencils.h:622
Vec3< ValueType > gradient() const
Definition: Stencils.h:680
WenoStencil(const GridType &grid)
Definition: Stencils.h:624
WenoStencil(const GridType &grid, double dx)
Definition: Stencils.h:632
ValueType laplacian() const
Definition: Stencils.h:692
Vec3< ValueType > gradient(const Vec3< ValueType > &V) const
Definition: Stencils.h:664
bool zeroCrossing() const
Definition: Stencils.h:702
typename GridT::TreeType TreeType
Definition: Stencils.h:619
typename GridT::ValueType ValueType
Definition: Stencils.h:620
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:711
GridT GridType
Definition: Stencils.h:618
ValueType normSqGrad(ValueType isoValue=ValueType(0)) const
Return the norm-square of the WENO upwind gradient (computed via WENO upwinding and Godunov's scheme)...
Definition: Stencils.h:645
Definition: NanoVDB.h:208
uint32_t CountOn(uint64_t v)
Definition: NanoVDB.h:1931
Type Min(Type a, Type b)
Definition: NanoVDB.h:758
T Pow4(T x)
Definition: NanoVDB.h:849
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: NanoVDB.h:902
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
T Pow3(T x)
Definition: NanoVDB.h:843
CoordT RoundDown(const Vec3T< RealT > &xyz)
Definition: NanoVDB.h:895
T Pow2(T x)
Definition: NanoVDB.h:837
Type Max(Type a, Type b)
Definition: NanoVDB.h:779
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
Definition: Stencils.h:226
int count() const
Definition: Stencils.h:234
bool test(int i) const
Definition: Stencils.h:230
bool none() const
Definition: Stencils.h:233
uint8_t bits
Definition: Stencils.h:227
bool all() const
Definition: Stencils.h:232
void set(int i)
Definition: Stencils.h:229
bool any() const
Definition: Stencils.h:231
Mask()
Definition: Stencils.h:228