HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GU_VDBFromParticleFluid.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: GU_VDBFromParticleFluid.h (GU Library, C++)
7  *
8  * COMMENTS:s
9  */
10 
11 #ifndef GU_VDBFROMPARTICLEFLUID_H_
12 #define GU_VDBFROMPARTICLEFLUID_H_
13 
14 #include <openvdb/Types.h>
15 #include <openvdb/Grid.h>
17 
18 #include <UT/UT_Array.h>
19 #include <UT/UT_ParallelUtil.h>
20 #include <GEO/GEO_PrimVDB.h>
21 
22 #include "GU_Detail.h"
23 #include "GU_API.h"
24 
25 #define USE_W2SUM_MASS 0
26 #define USE_DIST 0
27 #define USE_WEIGHTED_DIST 0
28 #define USE_VELOCITY 0
29 
30 
31 template<typename RT>
33 {
34 public:
35  typedef RT RealT;
38 
40  {
41  PHANTOM_PARTICLES = 0,
42  REFINED_AVERAGE_POSITION
43  };
44 
46  {
47  BUCKET_SAMPLE = 0,
49  AP_SAMPLE
50  };
51 
52  struct Parms
53  {
54  Parms() :
55  sampling_method(SMOOTH_SAMPLE),
56  surfacing_method(PHANTOM_PARTICLES),
57  particle_separation(0.1),
58  voxel_size(0.5),
59  influence_radius(3),
60  surface_distance(1),
61  resampling_iterations(-1),
62  resample_input(true),
63  tile_size(10),
64  additional_padding(0.5),
65  bandwidth_voxels(3),
66  memory_usage(-1),
67  time_step(0),
68  velocity_effect(0),
69  smoothing_iterations(0),
70  feature_preservation(0.0),
71  sheeting_scale(1.0),
72  feature_enhancement(0.3),
73  tile_grid_orig(Vector3R(0.0, 0.0, 0.0))
74  {}
75 
78 
79  // A base unit length, e.g. particle separation, based on which other
80  // parameters are meseared.
82 
83  // Voxel size measured in particle_separations
85 
86  // Common influence radius of all particles measured in particle_separations
88 
89  // Default droplet size measured in particle_separations if pscale attribute
90  // is not present or in (per-particle) pscales when there is a pscale point
91  // attribute
93 
94  // Maximum number of resampling iterations, -1 means unlimited
96 
97  // Half the width of narrow band in voxels
99 
100  // Total allowed work buffer memory usage in KB.
101  // Memory used to for hold the outputis excluded.
102  // -1 means unlimited.
104 
105  // Side length of a tile (excluding padding) measured as a multiple of
106  // influence radius
108 
109  // Extra padding (in addition of influence radius that is added by default)
110  // measured as a multiple of influence radius
112 
113  // The simulation time step (used in velocity based particle interation)
115 
116  // Degree of integration of velocity effect ranging between 0.0 (minimum) and 1.0 (maximum)
118 
119  // If true, the input particle set is resamped for the given influence radius.
120  // This substantially improves the performance and eliminates hideous artifacts
121  // of irregularities in particle samling
123 
124  // experimental, currently ignored
127 
130 
132  };
133 
134 
135 private:
136 
137  struct TileParms
138  {
139  TileParms() :
140  surfacing_method(PHANTOM_PARTICLES),
141  sampling_method(SMOOTH_SAMPLE),
142  world_orig(0.0, 0.0, 0.0),
143  world_orig_voxel(0, 0, 0),
144  voxel_size(-1.0),
145  tile_voxels(-1.0),
146  padding_voxels(-1.0),
147  influence_voxels(-1.0),
148  surface_distance(-1.0),
149  sampling_frequency(-1.0),
150  sampling_iterations(0),
151  time_step(-1.0),
152  smoothing_iterations(0),
153  feature_preservation(0.0),
154  sheeting_scale(1.0),
155  feature_enhancement(0.3)
156  {}
157 
158  SurfacingMethod surfacing_method;
159  SamplingMethod sampling_method;
160  Vector3R world_orig;
161  CoordT world_orig_voxel;
162  RealT voxel_size;
163  exint tile_voxels;
164  exint padding_voxels;
165  RealT influence_voxels;
166  RealT surface_distance;
167  RealT particle_separation;
168  RealT sampling_frequency;
169  int sampling_iterations;
170  RealT time_step;
171  int smoothing_iterations;
172  RealT feature_preservation;
173  RealT sheeting_scale;
174  RealT feature_enhancement;
175  };
176 
177  struct GridPoint
178  {
179  union {
180  RealT wsum; // sum of weights
181  RealT tmp; // reused as tmp in smoothing
182  };
183 
184  union {
185  RealT rad; // weighted radus (surface distance) sum
186  RealT final; // reused for final density values
187  };
188 
189  Vector3R pos; // wieighted position sum
190 
191  int16 iter; // last iteration that updated the point
192 
193 #if USE_VELOCITY
194  RealT vel;
195 #endif
196 
197 #if USE_DIST
198  RealT dist;
199 #endif
200 
201 #if USE_WEIGHTED_DIST
202  RealT distwm;
203 #endif
204 
205  inline Vector3R
206  center() { return pos / wsum; }
207  inline RealT
208  density() { return (pos.length() - rad) / wsum; }
209  inline RealT
210  height() { return pos.length() / wsum; }
211  inline RealT
212  radius() { return rad / wsum; }
213  };
214 
215 
216  struct SamplingGridPoint
217  {
218  int16 iter; // last iteration to have updated the point
219  RealT wsum; // sum of weights
220 
221  RealT mass;
222 
223 #if USE_VELOCITY
224  RealT vel;
225 #endif
226  Vector3R pos; // wieighted sum of positions
227  RealT rad; // weighted sum of radii (surface distances)
228 
229  RealT oldmass;
230  Vector3R oldpos;
231  };
232 
233  typedef typename openvdb::math::Vec2<int> Vec2T;
234  typedef typename openvdb::tree::Tree4<Vec2T, 5, 4, 3>::Type TileTreeT;
235  typedef typename TileTreeT::LeafNodeType TileLeafT;
236  typedef typename openvdb::Grid<TileTreeT> TileGridT;
237  typedef typename TileGridT::Ptr TileGridPtr;
238  typedef typename TileGridT::Accessor TileAccessorT;
239  typedef typename TileGridT::ConstAccessor TileConstAccessorT;
240  typedef typename openvdb::tree::LeafManager<TileTreeT> TileLeafManagerT;
241  typedef typename TileLeafManagerT::LeafRange TileRangeT;
243 
244 public:
246  typedef typename openvdb::Grid<RealTreeT> RealGridT;
247  typedef typename RealGridT::Ptr RealGridPtr;
248  typedef typename RealTreeT::LeafNodeType RealLeafT;
249  typedef typename RealGridT::Accessor RealAccessorT;
250  typedef typename RealGridT::ConstAccessor ConstRealAccessorT;
251  typedef typename openvdb::tree::LeafManager<RealTreeT> RealLeafManagerT;
252  typedef typename RealLeafManagerT::LeafRange RealRangeT;
253 
254 
255  GU_VDBFromParticleFluid(const GEO_Detail *gdp, const Parms &);
256 
258 
260 
261  RealGridPtr getGrid() { return myOutputVDBGrid; }
262 
264  {
265  return myErrorMessage.str();
266  }
267 
268  int evaluate(RealGridPtr result_grid, bool threaded = true);
269  void drawGuide(GU_Detail *gdp, RealT padding);
270  void drawGuide1(GU_Detail *gdp);
271  void drawGuide2(GU_Detail *gdp);
272 
273 
274 
275  void operator()(const UT_BlockedRange<exint>& range);
276  void join(GU_VDBFromParticleFluid& other);
277 
278 private:
279 
280  SYS_PACKED_STRUCT_HINT_BEGIN(GridIdxKey, 4)
281  {
282  exint grididx;
283  int key;
285 
286  class ut_isGridIdxLess
287  {
288  public:
289  inline bool operator()(const GridIdxKey &a,
290  const GridIdxKey &b) const
291  {
292  if (a.grididx == b.grididx)
293  return a.key < b.key;
294  return a.grididx < b.grididx;
295  }
296  };
297 
298 
299  struct SentinelVoxelsCleaner
300  {
301  void run(RealGridPtr grid);
302  void operator()(const RealRangeT& r) const;
303  private:
304  RealGridPtr myGrid;
305  };
306 
307  RealT bytesPerVoxel(RealT sampling_frequency_in_voxels);
308 
309  // convert the coordinates of a tile into a linear index
310  inline exint
311  tileCoordIndex(exint i, exint j, exint k)
312  {
313  return i * myMaxJK + j * myMaxK + k;
314  }
315 
316  // convert from linear index to a tile coordinate
317  inline void indexTileCoord(exint idx, exint &i, exint &j, exint &k)
318  {
319  i = idx / myMaxJK;
320  idx -= i * myMaxJK;
321  j = idx / myMaxK;
322  idx -= j * myMaxK;
323  k = idx;
324  }
325 
326 
327  struct Task;
328  struct ParticleInstancer;
329  struct InstanceMerger;
330 
331  void instanceParticles(RealT, RealT, UT_Array<GridIdxKey> &);
332  void buildTiles();
333 
334  class sop_IsTaskLess
335  {
336  public:
337  inline bool operator()(const TileLeafT::ValueOnIter &a,
338  const TileLeafT::ValueOnIter &b) const
339  {
340  Vec2T vala = a.getValue();
341  Vec2T valb = b.getValue();
342  return vala(1) - vala(0) < valb(1) - valb(0);
343  }
344  };
345 
346 
347  // some utility functions
348  inline Vector3R vecFloor(const Vector3R v)
349  {
350  return Vector3R(SYSfloor(v(0)), SYSfloor(v(1)), SYSfloor(v(2)));
351  }
352 
353  inline Vector3R vecCeil(const Vector3R v)
354  {
355  return Vector3R(SYSceil(v(0)), SYSceil(v(1)), SYSceil(v(2)));
356  }
357 
358  inline Vector3R vecRound(const Vector3R v)
359  {
360  return Vector3R(SYSfloor(v(0) + RealT(0.5)), SYSfloor(v(1) + RealT(0.5)), SYSfloor(v(2) + RealT(0.5)));
361  }
362 
363  inline RealT maxComp(const Vector3R v)
364  {
365  return SYSmax(v(0), v(1), v(2));
366  }
367 
368  inline RealT minComp(const Vector3R v)
369  {
370  return SYSmin(v(0), v(1), v(2));
371  }
372 
373  class Tile
374  {
375  public:
376  Tile(const GEO_Detail *gdp, const TileParms &parms, RealAccessorT& acc, UT_IntArray &keys);
377  void position(exint start_idx, exint stop_idx, Vector3R tile_orig, CoordT tile_orig_voxel);
378  void run(bool resample_input, int band_voxels);
379 
380  private:
381  void reset();
382 
383  inline exint
384  coordToIndex(CoordT v, exint grid_side_points, exint grid_side_points_sqr)
385  {
386  return grid_side_points_sqr * v(0) +
387  grid_side_points * v(1) + v(2);
388  }
389 
390  inline exint
391  coordToIndex(exint i, exint j, exint k, exint grid_side_points,
392  exint grid_side_points_sqr)
393  {
394  return grid_side_points_sqr * i + grid_side_points * j + k;
395  }
396 
397  inline CoordT
398  indexToCoord(exint index, exint grid_side_points, exint grid_side_points_sqr)
399  {
400  exint i = index / grid_side_points_sqr;
401  index -= i * grid_side_points_sqr;
402  exint j = index / grid_side_points;
403  index -= j * grid_side_points;
404  return CoordT(i, j, index);
405  }
406 
407  // p is in buffer coordinates (voxel size is unit)
408  // r is the particle surface distance value in voxels
409  inline void
410  pasteParticle(const Vector3R &pos, RealT rad,
411  RealT mass, RealT squeeze = RealT(1.0));
412 
413  inline void
414  pasteParticleToSamplingGrid(const Vector3R &pos, RealT rad,
415  RealT mass, RealT squeeze = RealT(1.0));
416 
417 
418  // Trilinear estimation of density at an arbitrary location
419  inline RealT evalDensity(Vector3R xyz);
420  inline RealT evalIsoDensity(Vector3R xyz, RealT *dist = NULL);
421 
422  inline RealT adjustedHeight(exint voxel_index, RealT srad);
423 
424  void expandIndexList(UT_IntArray &idx_list, int num_voxels, exint grid_side_points);
425 
426  void markBandVoxels(int band_radius, bool mark_all = false);
427 
428 #if CLIP_DEEP_PARTICLES
429  void separateDeepAndShallowSamples();
430  void floodIndexList(UT_IntArray &idx_list, bool ascending);
431  void floodDeepVoxels();
432 #endif
433 
434  void finalizeDensities();
435  void writeToVDB(int band_voxels);
436 
437  void collectNewSamples();
438 
439 
440  inline void recordSample(const Vector3R &pos, RealT rad, RealT mass = RealT(1.0),
441  RealT vel = RealT(0.0), bool clean_voxel_only = true);
442 
443  inline int findSamplingVoxelCorners(const Vector3R &pos, exint *offsets, RealT *weights);
444 
445  void resampleInputParticles();
446  exint addAllParticles();
447  exint pasteNewSamples();
448 
449 
450 #if EXACT_NN
451  void buildPointTree();
452  int nearestParticle(Vector3R &xyz);
453 #endif
454 
455 
456  inline RealT softClipAbove(RealT min, RealT max, RealT clip, RealT x)
457  {
458  if (x < min)
459  return x;
460  if (x > max)
461  return clip;
462  return ((max - x) * x + (x - min) * clip)/(max - min);
463  }
464 
465  inline RealT softClipBelow(RealT min, RealT max, RealT clip, RealT x)
466  {
467  if (x < min)
468  return clip;
469  if (x > max)
470  return x;
471  return ((max - x) * clip + (x - min) * x) / (max - min);
472  }
473 
474  void boxFilter(RealT rwidth, int iterations, RealT bias);
475  void boxFilterV2(int width, int iterations, RealT bias);
476 
477  private:
478 
479  const GEO_Detail *myGdp;
480  UT_IntArray &myKeys;
481  RealAccessorT myAcc;
482 
483  exint myStartIdx, myStopIdx;
484 
485  RealT myTimeStep;
486 
487  // Buffer Grid
488  RealT myBufferVoxelSize;
489  exint myBufferGridSidePoints;
490  exint myBufferGridSidePointsSqr;
491  exint myBufferGridSize;
492 
493  UT_Array<GridPoint> myBufferGrid;
494 
495  // Sampling Grid
496  RealT mySamplingVoxelSize;
497  exint mySamplingGridSidePoints;
498  exint mySamplingGridSidePointsSqr;
499  exint mySamplingGridSize;
500  UT_Array<SamplingGridPoint> mySamplingGrid;
501 
502  UT_IntArray myDirtySamplingIndices;
503 #if CLIP_DEEP_PARTICLES
504  UT_IntArray myDeepSamplingIndices;
505 #endif
506  UT_IntArray myLoadedBufferIndices;
507  UT_IntArray myDirtyBufferIndices;
508 
509  // Sampling buffer origin relative to the buffer (local) origin
510  Vector3R mySamplingOrig;
511 
512  // Origin of the tile in global position
513  Vector3R myTileOrig;
514 
515  CoordT myTileOrigVoxel;
516  Vector3R myWorldOrig;
517  CoordT myWorldOrigVoxel;
518  Vector3R myBufferOrig;
519  CoordT myBufferOrigVoxel;
520 
521  exint myPaddingVoxels;
522  RealT mySamplingFrequency;
523  int mySamplingIterations;
524  RealT mySurfaceDistance;
525  RealT myParticleSeparation;
526  RealT myInfluenceVoxels;
527  RealT myInfluenceVoxelsSqr;
528 
529  int mySmoothingIterations;
530  RealT myFeaturePreservation;
531 
532 #if EXACT_NN
533  GEO_PointTreeT<exint> myPointTree;
534  UT_Vector3Array myPoints;
535 #endif
536 
537  int16 myCurrentIteration;
538 
539  RealT mySheetingScale;
540  RealT myFeatureEnhancement;
541 
542  SamplingMethod mySamplingMethod;
543  SurfacingMethod mySurfacingMethod;
544 
545  UT_Interrupt *myBoss;
546  };
547 
548 private:
549 
550  const GEO_Detail *myGdp;
551 
552  RealT myTimeStep;
553  RealT myInfluenceRadius;
554  RealT myInfluenceVoxels;
555  RealT mySamplingFrequency;
556  int mySamplingIterations;
557 
558  UT_IntArray myHardKeys;
559  UT_IntArray &myKeys;
560 
561  exint myMaxJK, myMaxI, myMaxJ, myMaxK;
562  Vector3R myBoundingBoxOrig;
563  CoordT myBoundingBoxOrigVoxel;
564 
565  RealT myTilePadding;
566  exint myTilePaddingVoxels;
567  exint myTileSideVoxels;
568 
569  RealT myTileSideLength;
570  Vector3R myTileGridOrig;
571 
572  // VDB Grid storing for each tile the pair of its start and end indices
573  TileGridPtr myTileIdxGrid;
574 
575  RealT myVoxelSize;
576  RealGridPtr myOutputVDBGrid;
577 
578  bool myResampleInput;
579  int myBandVoxels;
580  RealT mySurfaceDistance;
581  RealT myParticleSeparation;
582 
583  TaskArrayT myHardTasks;
584  TaskArrayT &myTasks;
585  int myNumTiles;
586  int myMaxThreads;
587 
588  std::stringstream myErrorMessage;
589  int myError;
590 
591  int mySmoothingIterations;
592  RealT myFeaturePreservation;
593 
594  SamplingMethod mySamplingMethod;
595  SurfacingMethod mySurfacingMethod;
596 
597  bool myIsSlave;
598  SYS_AtomicCounter *myNextTileToDo;
599  SYS_AtomicCounter *myCompletedParticles;
600 
601  RealT mySheetingScale;
602  RealT myFeatureEnhancement;
603 };
604 
605 #endif /* GU_VDBFROMPARTICLEFLUID_H_ */
#define SYSmax(a, b)
Definition: SYS_Math.h:1570
GA_API const UT_StringHolder dist
#define SYS_PACKED_STRUCT_HINT_END
Definition: SYS_Types.h:469
openvdb::tree::Tree4< RealT, 5, 4, 3 >::Type RealTreeT
GLenum GLint * range
Definition: glcorearb.h:1925
RealGridT::Accessor RealAccessorT
#define SYS_PACKED_STRUCT_HINT_BEGIN(name, n)
Definition: SYS_Types.h:468
const GLdouble * v
Definition: glcorearb.h:837
GLsizei const GLchar *const * string
Definition: glcorearb.h:814
RealTreeT::LeafNodeType RealLeafT
int64 exint
Definition: SYS_Types.h:125
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
RealGridT::ConstAccessor ConstRealAccessorT
openvdb::tree::LeafManager< RealTreeT > RealLeafManagerT
const std::string getErrorMessage()
GLuint GLsizei const GLuint const GLintptr * offsets
Definition: glcorearb.h:2621
GLint GLsizei GLsizei height
Definition: glcorearb.h:103
GLboolean reset
Definition: glad.h:5138
UT_Vector3T< RealT > Vector3R
openvdb::Grid< RealTreeT > RealGridT
#define GU_API
Definition: GU_API.h:14
SYS_API fpreal32 SYSfloor(fpreal32 val)
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
GLint GLenum GLint x
Definition: glcorearb.h:409
GLint j
Definition: glad.h:2733
GA_API const UT_StringHolder mass
short int16
Definition: SYS_Types.h:37
GLuint index
Definition: glcorearb.h:786
RealLeafManagerT::LeafRange RealRangeT
SIM_API const UT_StringHolder position
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GLint GLsizei width
Definition: glcorearb.h:103
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
GLboolean r
Definition: glcorearb.h:1222
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
#define SYSmin(a, b)
Definition: SYS_Math.h:1571
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T clip(const T &p, const Box< T > &box) IMATH_NOEXCEPT
Definition: ImathBoxAlgo.h:29
Declare prior to use.
auto join(It begin, Sentinel end, string_view sep) -> join_view< It, Sentinel >
Definition: format.h:2559
SYS_API fpreal32 SYSceil(fpreal32 val)