HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PointRasterizeFrustumImpl.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 Dan Bailey
5 ///
6 /// @file PointRasterizeFrustumImpl.h
7 ///
8 
9 #ifndef OPENVDB_POINTS_POINT_RASTERIZE_FRUSTUM_IMPL_HAS_BEEN_INCLUDED
10 #define OPENVDB_POINTS_POINT_RASTERIZE_FRUSTUM_IMPL_HAS_BEEN_INCLUDED
11 
12 namespace openvdb {
14 namespace OPENVDB_VERSION_NAME {
15 namespace points {
16 
17 /// @cond OPENVDB_DOCS_INTERNAL
18 
19 namespace point_rasterize_internal {
20 
21 
22 // By default this trait simply no-op copies the filter
23 template <typename FilterT>
24 struct RasterGroupTraits
25 {
26  using NewFilterT = FilterT;
27  static NewFilterT resolve(const FilterT& filter, const points::AttributeSet&)
28  {
29  return filter;
30  }
31 };
32 
33 // For RasterGroups objects, this returns a new MultiGroupFilter with names-to-indices resolved
34 template <>
35 struct RasterGroupTraits<RasterGroups>
36 {
37  using NewFilterT = points::MultiGroupFilter;
38  static NewFilterT resolve(const RasterGroups& groups, const points::AttributeSet& attributeSet)
39  {
40  return NewFilterT(groups.includeNames, groups.excludeNames, attributeSet);
41  }
42 };
43 
44 
45 // Traits to indicate bool and ValueMask are bool types
46 template <typename T> struct BoolTraits { static const bool IsBool = false; };
47 template <> struct BoolTraits<bool> { static const bool IsBool = true; };
48 template <> struct BoolTraits<ValueMask> { static const bool IsBool = true; };
49 
50 
51 struct TrueOp {
52  bool mOn;
53  explicit TrueOp(double scale) : mOn(scale > 0.0) { }
54  template<typename ValueType>
55  void operator()(ValueType& v) const { v = static_cast<ValueType>(mOn); }
56 };
57 
58 
59 template <typename ValueT>
60 typename std::enable_if<std::is_integral<typename ValueTraits<ValueT>::ElementType>::value, ValueT>::type
61 castValue(const double value)
62 {
63  return ValueT(math::Ceil(value));
64 }
65 
66 
67 template <typename ValueT>
68 typename std::enable_if<!std::is_integral<typename ValueTraits<ValueT>::ElementType>::value, ValueT>::type
69 castValue(const double value)
70 {
71  return static_cast<ValueT>(value);
72 }
73 
74 
75 template <typename ValueT>
76 typename std::enable_if<!ValueTraits<ValueT>::IsVec, bool>::type
77 greaterThan(const ValueT& value, const float threshold)
78 {
79  return value >= static_cast<ValueT>(threshold);
80 }
81 
82 
83 template <typename ValueT>
84 typename std::enable_if<ValueTraits<ValueT>::IsVec, bool>::type
85 greaterThan(const ValueT& value, const float threshold)
86 {
87  return static_cast<double>(value.lengthSqr()) >= threshold*threshold;
88 }
89 
90 
91 template <typename AttributeT, typename HandleT, typename StridedHandleT>
92 typename std::enable_if<!ValueTraits<AttributeT>::IsVec, AttributeT>::type
93 getAttributeScale(HandleT& handlePtr, StridedHandleT&, Index index)
94 {
95  if (handlePtr) {
96  return handlePtr->get(index);
97  }
98  return AttributeT(1);
99 }
100 
101 
102 template <typename AttributeT, typename HandleT, typename StridedHandleT>
103 typename std::enable_if<ValueTraits<AttributeT>::IsVec, AttributeT>::type
104 getAttributeScale(HandleT& handlePtr, StridedHandleT& stridedHandlePtr, Index index)
105 {
106  if (handlePtr) {
107  return handlePtr->get(index);
108  } else if (stridedHandlePtr) {
109  return AttributeT(
110  stridedHandlePtr->get(index, 0),
111  stridedHandlePtr->get(index, 1),
112  stridedHandlePtr->get(index, 2));
113  }
114  return AttributeT(1);
115 }
116 
117 
118 template <typename ValueT>
119 struct MultiplyOp
120 {
121  template <typename AttributeT>
122  static ValueT mul(const ValueT& a, const AttributeT& b)
123  {
124  return a * b;
125  }
126 };
127 
128 template <>
129 struct MultiplyOp<bool>
130 {
131  template <typename AttributeT>
132  static bool mul(const bool& a, const AttributeT& b)
133  {
134  return a && (b != zeroVal<AttributeT>());
135  }
136 };
137 
138 template <typename PointDataGridT, typename AttributeT, typename GridT,
139  typename FilterT>
140 struct RasterizeOp
141 {
142  using TreeT = typename GridT::TreeType;
143  using LeafT = typename TreeT::LeafNodeType;
144  using ValueT = typename TreeT::ValueType;
145  using PointLeafT = typename PointDataGridT::TreeType::LeafNodeType;
146  using PointIndexT = typename PointLeafT::ValueType;
147  using CombinableT = tbb::combinable<GridT>;
148  using SumOpT = tools::valxform::SumOp<ValueT>;
149  using MaxOpT = tools::valxform::MaxOp<ValueT>;
150  using PositionHandleT = openvdb::points::AttributeHandle<Vec3f>;
151  using VelocityHandleT = openvdb::points::AttributeHandle<Vec3f>;
152  using RadiusHandleT = openvdb::points::AttributeHandle<float>;
153 
154  // to prevent checking the interrupter too frequently, only check every 32 voxels cubed
155  static const int interruptThreshold = 32*32*32;
156 
157  RasterizeOp(
158  const PointDataGridT& grid,
159  const std::vector<Index64>& offsets,
160  const size_t attributeIndex,
161  const Name& velocityAttribute,
162  const Name& radiusAttribute,
163  CombinableT& combinable,
164  CombinableT* weightCombinable,
165  const bool dropBuffers,
166  const float scale,
167  const FilterT& filter,
168  const bool computeMax,
169  const bool alignedTransform,
170  const bool staticCamera,
171  const FrustumRasterizerSettings& settings,
172  const FrustumRasterizerMask& mask,
173  util::NullInterrupter* interrupt)
174  : mGrid(grid)
175  , mOffsets(offsets)
176  , mAttributeIndex(attributeIndex)
177  , mVelocityAttribute(velocityAttribute)
178  , mRadiusAttribute(radiusAttribute)
179  , mCombinable(combinable)
180  , mWeightCombinable(weightCombinable)
181  , mDropBuffers(dropBuffers)
182  , mScale(scale)
183  , mFilter(filter)
184  , mComputeMax(computeMax)
185  , mAlignedTransform(alignedTransform)
186  , mStaticCamera(staticCamera)
187  , mSettings(settings)
188  , mMask(mask)
189  , mInterrupter(interrupt) { }
190 
191  template <typename PointOpT>
192  static void rasterPoint(const Coord& ijk, const double scale,
193  const AttributeT& attributeScale, PointOpT& op)
194  {
195  op(ijk, scale, attributeScale);
196  }
197 
198  template <typename PointOpT>
199  static void rasterPoint(const Vec3d& position, const double scale,
200  const AttributeT& attributeScale, PointOpT& op)
201  {
202  Coord ijk = Coord::round(position);
203  op(ijk, scale, attributeScale);
204  }
205 
206  template <typename SphereOpT>
207  static void rasterVoxelSphere(const Vec3d& position, const double scale,
208  const AttributeT& attributeScale, const float radius, util::NullInterrupter* interrupter, SphereOpT& op)
209  {
210  assert(radius > 0.0f);
211  Coord ijk = Coord::round(position);
212  int &i = ijk[0], &j = ijk[1], &k = ijk[2];
213  const int imin=math::Floor(position[0]-radius), imax=math::Ceil(position[0]+radius);
214  const int jmin=math::Floor(position[1]-radius), jmax=math::Ceil(position[1]+radius);
215  const int kmin=math::Floor(position[2]-radius), kmax=math::Ceil(position[2]+radius);
216 
217  const bool interrupt = interrupter && (imax-imin)*(jmax-jmin)*(kmax-kmin) > interruptThreshold;
218 
219  for (i = imin; i <= imax; ++i) {
220  if (interrupt && interrupter->wasInterrupted()) break;
221  const auto x2 = math::Pow2(i - position[0]);
222  for (j = jmin; j <= jmax; ++j) {
223  if (interrupt && interrupter->wasInterrupted()) break;
224  const auto x2y2 = math::Pow2(j - position[1]) + x2;
225  for (k = kmin; k <= kmax; ++k) {
226  const auto x2y2z2 = x2y2 + math::Pow2(k - position[2]);
227  op(ijk, scale, attributeScale, x2y2z2, radius*radius);
228  }
229  }
230  }
231  }
232 
233  template <typename SphereOpT>
234  static void rasterApproximateFrustumSphere(const Vec3d& position, const double scale,
235  const AttributeT& attributeScale, const float radiusWS,
236  const math::Transform& frustum, const CoordBBox* clipBBox,
237  util::NullInterrupter* interrupter, SphereOpT& op)
238  {
239  Vec3d voxelSize = frustum.voxelSize(position);
240  Vec3d radius = Vec3d(radiusWS)/voxelSize;
241  CoordBBox frustumBBox(Coord::floor(position-radius), Coord::ceil(position+radius));
242 
243  // if clipping to frustum is enabled, clip the index bounding box
244  if (clipBBox) {
245  frustumBBox.intersect(*clipBBox);
246  }
247 
248  Vec3i outMin = frustumBBox.min().asVec3i();
249  Vec3i outMax = frustumBBox.max().asVec3i();
250 
251  const bool interrupt = interrupter && frustumBBox.volume() > interruptThreshold;
252 
253  // back-project each output voxel into the input tree.
254  Vec3R xyz;
255  Coord outXYZ;
256  int &x = outXYZ.x(), &y = outXYZ.y(), &z = outXYZ.z();
257  for (x = outMin.x(); x <= outMax.x(); ++x) {
258  if (interrupt && interrupter->wasInterrupted()) break;
259  xyz.x() = x;
260  for (y = outMin.y(); y <= outMax.y(); ++y) {
261  if (interrupt && interrupter->wasInterrupted()) break;
262  xyz.y() = y;
263  for (z = outMin.z(); z <= outMax.z(); ++z) {
264  xyz.z() = z;
265 
266  // approximate conversion to world-space
267 
268  Vec3d offset = (xyz - position) * voxelSize;
269 
270  double distanceSqr = offset.dot(offset);
271 
272  op(outXYZ, scale, attributeScale, distanceSqr, radiusWS*radiusWS);
273  }
274  }
275  }
276  }
277 
278  template <typename SphereOpT>
279  static void rasterFrustumSphere(const Vec3d& position, const double scale,
280  const AttributeT& attributeScale, const float radiusWS,
281  const math::Transform& frustum, const CoordBBox* clipBBox,
282  util::NullInterrupter* interrupter, SphereOpT& op)
283  {
284  const Vec3d positionWS = frustum.indexToWorld(position);
285 
286  BBoxd inputBBoxWS(positionWS-radiusWS, positionWS+radiusWS);
287  // Transform the corners of the input tree's bounding box
288  // and compute the enclosing bounding box in the output tree.
289  Vec3R frustumMin;
290  Vec3R frustumMax;
291  math::calculateBounds(frustum, inputBBoxWS.min(), inputBBoxWS.max(),
292  frustumMin, frustumMax);
293  CoordBBox frustumBBox(Coord::floor(frustumMin), Coord::ceil(frustumMax));
294 
295  // if clipping to frustum is enabled, clip the index bounding box
296  if (clipBBox) {
297  frustumBBox.intersect(*clipBBox);
298  }
299 
300  Vec3i outMin = frustumBBox.min().asVec3i();
301  Vec3i outMax = frustumBBox.max().asVec3i();
302 
303  const bool interrupt = interrupter && frustumBBox.volume() > interruptThreshold;
304 
305  // back-project each output voxel into the input tree.
306  Vec3R xyz;
307  Coord outXYZ;
308  int &x = outXYZ.x(), &y = outXYZ.y(), &z = outXYZ.z();
309  for (x = outMin.x(); x <= outMax.x(); ++x) {
310  if (interrupt && interrupter->wasInterrupted()) break;
311  xyz.x() = x;
312  for (y = outMin.y(); y <= outMax.y(); ++y) {
313  if (interrupt && interrupter->wasInterrupted()) break;
314  xyz.y() = y;
315  for (z = outMin.z(); z <= outMax.z(); ++z) {
316  xyz.z() = z;
317 
318  Vec3R xyzWS = frustum.indexToWorld(xyz);
319 
320  double xDist = xyzWS.x() - positionWS.x();
321  double yDist = xyzWS.y() - positionWS.y();
322  double zDist = xyzWS.z() - positionWS.z();
323 
324  double distanceSqr = xDist*xDist+yDist*yDist+zDist*zDist;
325  op(outXYZ, scale, attributeScale, distanceSqr, radiusWS*radiusWS);
326  }
327  }
328  }
329  }
330 
331  void operator()(const PointLeafT& leaf, size_t) const
332  {
333  if (mInterrupter && mInterrupter->wasInterrupted()) {
335  return;
336  }
337 
338  using AccessorT = tree::ValueAccessor<typename GridT::TreeType>;
339  using MaskAccessorT = tree::ValueAccessor<const MaskTree>;
340 
341  constexpr bool isBool = BoolTraits<ValueT>::IsBool;
342  const bool isFrustum = !mSettings.transform->isLinear();
343 
344  const bool useRaytrace = mSettings.velocityMotionBlur || !mStaticCamera;
345  const bool useRadius = mSettings.useRadius;
346  const float radiusScale = mSettings.radiusScale;
347  const float radiusThreshold = std::sqrt(3.0f);
348 
349  const float threshold = mSettings.threshold;
350  const bool scaleByVoxelVolume = !useRadius && mSettings.scaleByVoxelVolume;
351 
352  const float shutterStartDt = mSettings.camera.shutterStart()/mSettings.framesPerSecond;
353  const float shutterEndDt = mSettings.camera.shutterEnd()/mSettings.framesPerSecond;
354  const int motionSteps = std::max(1, mSettings.motionSamples-1);
355 
356  std::vector<Vec3d> motionPositions(motionSteps+1, Vec3d());
357  std::vector<Vec2s> frustumRadiusSizes(motionSteps+1, Vec2s());
358 
359  const auto& pointsTransform = mGrid.constTransform();
360 
361  const float voxelSize = static_cast<float>(mSettings.transform->voxelSize()[0]);
362 
363  auto& grid = mCombinable.local();
364  auto& tree = grid.tree();
365  const auto& transform = *(mSettings.transform);
366 
367  grid.setTransform(mSettings.transform->copy());
368 
369  AccessorT valueAccessor(tree);
370 
371  std::unique_ptr<MaskAccessorT> maskAccessor;
372  auto maskTree = mMask.getTreePtr();
373  if (maskTree) maskAccessor.reset(new MaskAccessorT(*maskTree));
374 
375  const CoordBBox* clipBBox = !mMask.clipBBox().empty() ? &mMask.clipBBox() : nullptr;
376 
377  std::unique_ptr<AccessorT> weightAccessor;
378  if (mWeightCombinable) {
379  auto& weightGrid = mWeightCombinable->local();
380  auto& weightTree = weightGrid.tree();
381  weightAccessor.reset(new AccessorT(weightTree));
382  }
383 
384  // create a temporary tree when rasterizing with radius and accumulate
385  // or weighted average methods
386 
387  std::unique_ptr<TreeT> tempTree;
388  std::unique_ptr<TreeT> tempWeightTree;
389  std::unique_ptr<AccessorT> tempValueAccessor;
390  std::unique_ptr<AccessorT> tempWeightAccessor;
391  if (useRadius && !mComputeMax) {
392  tempTree.reset(new TreeT(tree.background()));
393  tempValueAccessor.reset(new AccessorT(*tempTree));
394  if (weightAccessor) {
395  tempWeightTree.reset(new TreeT(weightAccessor->tree().background()));
396  tempWeightAccessor.reset(new AccessorT(*tempWeightTree));
397  }
398  }
399 
400  // point rasterization
401 
402  // impl - modify a single voxel by coord, handles temporary trees and all supported modes
403  auto doModifyVoxelOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale,
404  const bool isTemp, const bool forceSum)
405  {
406  if (mMask && !mMask.valid(ijk, maskAccessor.get())) return;
407  if (isBool) {
408  // only modify the voxel if the attributeScale is positive and non-zero
409  if (!math::isZero<AttributeT>(attributeScale) && !math::isNegative<AttributeT>(attributeScale)) {
410  valueAccessor.modifyValue(ijk, TrueOp(scale));
411  }
412  } else {
413  ValueT weightValue = castValue<ValueT>(scale);
414  ValueT newValue = MultiplyOp<ValueT>::mul(weightValue, attributeScale);
415  if (scaleByVoxelVolume) {
416  newValue /= static_cast<ValueT>(transform.voxelSize(ijk.asVec3d()).product());
417  }
418  if (point_rasterize_internal::greaterThan(newValue, threshold)) {
419  if (isTemp) {
420  tempValueAccessor->modifyValue(ijk, MaxOpT(newValue));
421  if (tempWeightAccessor) {
422  tempWeightAccessor->modifyValue(ijk, MaxOpT(weightValue));
423  }
424  } else {
425  if (mComputeMax && !forceSum) {
426  valueAccessor.modifyValue(ijk, MaxOpT(newValue));
427  } else {
428  valueAccessor.modifyValue(ijk, SumOpT(newValue));
429  }
430  if (weightAccessor) {
431  weightAccessor->modifyValue(ijk, SumOpT(weightValue));
432  }
433  }
434  }
435  }
436  };
437 
438  // modify a single voxel by coord, disable temporary trees
439  auto modifyVoxelOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale)
440  {
441  doModifyVoxelOp(ijk, scale, attributeScale, /*isTemp=*/false, /*forceSum=*/false);
442  };
443 
444  // sum a single voxel by coord, no temporary trees or maximum mode
445  auto sumVoxelOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale)
446  {
447  doModifyVoxelOp(ijk, scale, attributeScale, /*isTemp=*/false, /*forceSum=*/true);
448  };
449 
450  // sphere rasterization
451 
452  // impl - modify a single voxel by coord based on distance from sphere origin
453  auto doModifyVoxelByDistanceOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale,
454  const double distanceSqr, const double radiusSqr, const bool isTemp)
455  {
456  if (distanceSqr >= radiusSqr) return;
457  if (isBool) {
458  valueAccessor.modifyValue(ijk, TrueOp(scale));
459  } else {
460  double distance = std::sqrt(distanceSqr);
461  double radius = std::sqrt(radiusSqr);
462  double result = 1.0 - distance/radius;
463  doModifyVoxelOp(ijk, result * scale, attributeScale, isTemp, /*forceSum=*/false);
464  }
465  };
466 
467  // modify a single voxel by coord based on distance from sphere origin, disable temporary trees
468  auto modifyVoxelByDistanceOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale,
469  const double distanceSqr, const double radiusSqr)
470  {
471  doModifyVoxelByDistanceOp(ijk, scale, attributeScale, distanceSqr, radiusSqr, /*isTemp=*/false);
472  };
473 
474  // modify a single voxel by coord based on distance from sphere origin, enable temporary trees
475  auto modifyTempVoxelByDistanceOp = [&](const Coord& ijk, const double scale, const AttributeT& attributeScale,
476  const double distanceSqr, const double radiusSqr)
477  {
478  doModifyVoxelByDistanceOp(ijk, scale, attributeScale, distanceSqr, radiusSqr, /*isTemp=*/true);
479  };
480 
481  typename points::AttributeHandle<AttributeT>::Ptr attributeHandle;
482  using ElementT = typename ValueTraits<AttributeT>::ElementType;
483  typename points::AttributeHandle<ElementT>::Ptr stridedAttributeHandle;
484 
485  if (mAttributeIndex != points::AttributeSet::INVALID_POS) {
486  const auto& attributeArray = leaf.constAttributeArray(mAttributeIndex);
487  if (attributeArray.stride() == 3) {
488  stridedAttributeHandle = points::AttributeHandle<ElementT>::create(attributeArray);
489  } else {
490  attributeHandle = points::AttributeHandle<AttributeT>::create(attributeArray);
491  }
492  }
493 
494  size_t positionIndex = leaf.attributeSet().find("P");
495  size_t velocityIndex = leaf.attributeSet().find(mVelocityAttribute);
496  size_t radiusIndex = leaf.attributeSet().find(mRadiusAttribute);
497 
498  auto positionHandle = PositionHandleT::create(
499  leaf.constAttributeArray(positionIndex));
500  auto velocityHandle = (useRaytrace && leaf.hasAttribute(velocityIndex)) ?
501  VelocityHandleT::create(leaf.constAttributeArray(velocityIndex)) :
502  VelocityHandleT::Ptr();
503  auto radiusHandle = (useRadius && leaf.hasAttribute(radiusIndex)) ?
504  RadiusHandleT::create(leaf.constAttributeArray(radiusIndex)) :
505  RadiusHandleT::Ptr();
506 
507  for (auto iter = leaf.beginIndexOn(mFilter); iter; ++iter) {
508 
509  // read attribute value if it exists, attributes with stride 3 are composited
510  // into a vector grid such as float[3] => vec3s
511 
512  const AttributeT attributeScale = getAttributeScale<AttributeT>(
513  attributeHandle, stridedAttributeHandle, *iter);
514 
515  float radiusWS = 1.0f, radius = 1.0f;
516  if (useRadius) {
517  radiusWS *= radiusScale;
518  if (radiusHandle) radiusWS *= radiusHandle->get(*iter);
519  if (isFrustum) {
520  radius = radiusWS;
521  } else if (voxelSize > 0.0f) {
522  radius = radiusWS / voxelSize;
523  }
524  }
525 
526  // frustum thresholding is done later to factor in changing voxel sizes
527 
528  bool doRadius = useRadius;
529  if (!isFrustum && radius < radiusThreshold) {
530  doRadius = false;
531  }
532 
533  float increment = shutterEndDt - shutterStartDt;
534 
535  // disable ray-tracing if velocity is very small
536 
537  openvdb::Vec3f velocity(0.0f);
538  bool doRaytrace = useRaytrace;
539  if (doRaytrace) {
540  if (increment < openvdb::math::Delta<float>::value()) {
541  doRaytrace = false;
542  } else {
543  if (velocityHandle) velocity = velocityHandle->get(*iter);
544  if (mStaticCamera && velocity.lengthSqr() < openvdb::math::Delta<float>::value()) {
545  doRaytrace = false;
546  }
547  }
548  }
549 
550  if (motionSteps > 1) increment /= float(motionSteps);
551 
552  Vec3d position = positionHandle->get(*iter) + iter.getCoord().asVec3d();
553 
554  if (doRaytrace) {
555  // raytracing is done in index space
556  position = pointsTransform.indexToWorld(position);
557 
558  for (int motionStep = 0; motionStep <= motionSteps; motionStep++) {
559 
560  float offset = motionStep == motionSteps ? shutterEndDt :
561  (shutterStartDt + increment * static_cast<float>(motionStep));
562  Vec3d samplePosition = position + velocity * offset;
563 
564  const math::Transform* sampleTransform = &transform;
565  if (!mSettings.camera.isStatic()) {
566  sampleTransform = &mSettings.camera.transform(motionStep);
567  if (mSettings.camera.hasWeight(motionStep)) {
568  const float weight = mSettings.camera.weight(motionStep);
569  const Vec3d referencePosition = transform.worldToIndex(samplePosition);
570  const Vec3d adjustedPosition = sampleTransform->worldToIndex(samplePosition);
571  motionPositions[motionStep] = referencePosition + (adjustedPosition - referencePosition) * weight;
572  } else {
573  motionPositions[motionStep] = sampleTransform->worldToIndex(samplePosition);
574  }
575  } else {
576  motionPositions[motionStep] = sampleTransform->worldToIndex(samplePosition);
577  }
578  if (doRadius && isFrustum) {
579  Vec3d left = sampleTransform->worldToIndex(samplePosition - Vec3d(radiusWS, 0, 0));
580  Vec3d right = sampleTransform->worldToIndex(samplePosition + Vec3d(radiusWS, 0, 0));
581  float width = static_cast<float>((right - left).length());
582  Vec3d top = sampleTransform->worldToIndex(samplePosition + Vec3d(0, radiusWS, 0));
583  Vec3d bottom = sampleTransform->worldToIndex(samplePosition - Vec3d(0, radiusWS, 0));
584  float height = static_cast<float>((top - bottom).length());
585  frustumRadiusSizes[motionStep].x() = width;
586  frustumRadiusSizes[motionStep].y() = height;
587  }
588  }
589 
590  double totalDistance = 0.0;
591  for (size_t i = 0; i < motionPositions.size()-1; i++) {
592  Vec3d direction = motionPositions[i+1] - motionPositions[i];
593  double distance = direction.length();
594  totalDistance += distance;
595  }
596 
597  double distanceWeight = totalDistance > 0.0 ? 1.0 / totalDistance : 1.0;
598 
599  if (doRadius && !mComputeMax) {
600  // mark all voxels inactive
601  for (auto leaf = tempTree->beginLeaf(); leaf; ++leaf) {
602  leaf->setValuesOff();
603  }
604  if (tempWeightAccessor) {
605  for (auto leaf = tempWeightTree->beginLeaf(); leaf; ++leaf) {
606  leaf->setValuesOff();
607  }
608  }
609  }
610 
611  for (int motionStep = 0; motionStep < motionSteps; motionStep++) {
612 
613  Vec3d startPosition = motionPositions[motionStep];
614  Vec3d endPosition = motionPositions[motionStep+1];
615 
616  Vec3d direction(endPosition - startPosition);
617  double distance = direction.length();
618 
619  // if rasterizing into a frustum grid, compute the index-space radii for start
620  // and end positions and if below the radius threshold disable using radius
621 
622  float maxRadius = radius;
623 
624  if (doRadius && isFrustum) {
625  const Vec2s& startRadius = frustumRadiusSizes[motionStep];
626  const Vec2s& endRadius = frustumRadiusSizes[motionStep+1];
627 
628  if (startRadius[0] < radiusThreshold && startRadius[1] < radiusThreshold &&
629  endRadius[0] < radiusThreshold && endRadius[1] < radiusThreshold) {
630  doRadius = false;
631  } else {
632  // max radius is the largest index-space radius factoring in
633  // that in frustum space a sphere is not rasterized as spherical
634  maxRadius = std::max(startRadius[0], startRadius[1]);
635  maxRadius = std::max(maxRadius, endRadius[0]);
636  maxRadius = std::max(maxRadius, endRadius[1]);
637  }
638  }
639 
640  if (doRadius) {
641  distanceWeight = std::min(distanceWeight, 1.0);
642 
643  // these arbitrary constants are how tightly spheres are packed together
644  // irrespective of how large they are in index-space - if it is too low, the shape of
645  // the individual spheres becomes visible in the rasterized path,
646  // if it is too high, rasterization becomes less efficient with
647  // diminishing returns towards the accuracy of the rasterized path
648  double spherePacking = mSettings.accurateSphereMotionBlur ? 4.0 : 2.0;
649 
650  const int steps = std::max(2, math::Ceil(distance * spherePacking / double(maxRadius)) + 1);
651 
652  Vec3d sample(startPosition);
653  const Vec3d offset(direction / (steps-1));
654 
655  for (int step = 0; step < steps; step++) {
656  if (isFrustum) {
657  if (mComputeMax) {
658  if (mSettings.accurateFrustumRadius) {
659  this->rasterFrustumSphere(sample, mScale * distanceWeight,
660  attributeScale, radiusWS, transform, clipBBox, mInterrupter, modifyVoxelByDistanceOp);
661  } else {
662  this->rasterApproximateFrustumSphere(sample, mScale * distanceWeight,
663  attributeScale, radiusWS, transform, clipBBox, mInterrupter, modifyVoxelByDistanceOp);
664  }
665  } else {
666  if (mSettings.accurateFrustumRadius) {
667  this->rasterFrustumSphere(sample, mScale * distanceWeight,
668  attributeScale, radiusWS, transform, clipBBox, mInterrupter, modifyTempVoxelByDistanceOp);
669  } else {
670  this->rasterApproximateFrustumSphere(sample, mScale * distanceWeight,
671  attributeScale, radiusWS, transform, clipBBox, mInterrupter, modifyTempVoxelByDistanceOp);
672  }
673  }
674  } else {
675  if (mComputeMax) {
676  this->rasterVoxelSphere(sample, mScale * distanceWeight,
677  attributeScale, radius, mInterrupter, modifyVoxelByDistanceOp);
678  } else {
679  this->rasterVoxelSphere(sample, mScale * distanceWeight,
680  attributeScale, radius, mInterrupter, modifyTempVoxelByDistanceOp);
681  }
682  }
683  if (step < (steps-1)) sample += offset;
684  else sample = endPosition;
685  }
686  } else {
687  // compute direction and store vector length as max time
688  mDdaRay.setMinTime(math::Delta<double>::value());
689  mDdaRay.setMaxTime(std::max(distance, math::Delta<double>::value()*2.0));
690 
691  // DDA requires normalized directions
692  direction.normalize();
693  mDdaRay.setDir(direction);
694 
695  // dda assumes node-centered ray-tracing, so compensate by adding half a voxel first
696  mDdaRay.setEye(startPosition + Vec3d(0.5));
697 
698  // clip the ray to the frustum bounding box
699  if (clipBBox) {
700  mDdaRay.clip(*clipBBox);
701  }
702 
703  // first rasterization in a subsequent DDA traversal should always sum contributions
704  // in order to smoothly stitch the beginning and end of two rays together
705  bool forceSum = motionStep > 0;
706 
707  mDda.init(mDdaRay);
708  while (true) {
709  const Coord& voxel = mDda.voxel();
710  double delta = (mDda.next() - mDda.time()) * distanceWeight;
711  if (forceSum) {
712  this->rasterPoint(voxel, mScale * delta,
713  attributeScale, sumVoxelOp);
714  forceSum = false;
715  } else {
716  this->rasterPoint(voxel, mScale * delta,
717  attributeScale, modifyVoxelOp);
718  }
719  if (!mDda.step()) break;
720  }
721  }
722  }
723 
724  if (doRadius && !mComputeMax) {
725  // copy values into valueAccessor
726  for (auto iter = tempTree->cbeginValueOn(); iter; ++iter) {
727  valueAccessor.modifyValue(iter.getCoord(), SumOpT(*iter));
728  }
729 
730  // copy values into weightAccessor
731  if (weightAccessor) {
732  for (auto iter = tempWeightTree->cbeginValueOn(); iter; ++iter) {
733  weightAccessor->modifyValue(iter.getCoord(), SumOpT(*iter));
734  }
735  }
736  }
737 
738  } else {
739  if (!mAlignedTransform) {
740  position = transform.worldToIndex(
741  pointsTransform.indexToWorld(position));
742  }
743 
744  if (doRadius) {
745  if (isFrustum) {
746  if (mSettings.accurateFrustumRadius) {
747  this->rasterFrustumSphere(position, mScale, attributeScale,
748  radiusWS, transform, clipBBox, mInterrupter, modifyVoxelByDistanceOp);
749  } else {
750  this->rasterApproximateFrustumSphere(position, mScale, attributeScale,
751  radiusWS, transform, clipBBox, mInterrupter, modifyVoxelByDistanceOp);
752  }
753  } else {
754  this->rasterVoxelSphere(position, mScale, attributeScale,
755  radius, mInterrupter, modifyVoxelByDistanceOp);
756  }
757  } else {
758  this->rasterPoint(position, mScale, attributeScale, modifyVoxelOp);
759  }
760  }
761  }
762 
763  // if drop buffers is enabled, swap the leaf buffer with a partial (empty) buffer,
764  // so the voxel data is de-allocated to reduce memory
765 
766  if (mDropBuffers) {
767  typename PointLeafT::Buffer emptyBuffer(PartialCreate(), zeroVal<PointIndexT>());
768  (const_cast<PointLeafT&>(leaf)).swap(emptyBuffer);
769  }
770  }
771 
772 // TODO: would be better to move some of the immutable values into a shared struct
773 
774 private:
775  mutable math::Ray<double> mDdaRay;
776  mutable math::DDA<openvdb::math::Ray<double>> mDda;
777  const PointDataGridT& mGrid;
778  const std::vector<Index64>& mOffsets;
779  const size_t mAttributeIndex;
780  const Name mVelocityAttribute;
781  const Name mRadiusAttribute;
782  CombinableT& mCombinable;
783  CombinableT* mWeightCombinable;
784  const bool mDropBuffers;
785  const float mScale;
786  const FilterT& mFilter;
787  const bool mComputeMax;
788  const bool mAlignedTransform;
789  const bool mStaticCamera;
790  const FrustumRasterizerSettings& mSettings;
791  const FrustumRasterizerMask& mMask;
792  util::NullInterrupter* mInterrupter;
793 }; // struct RasterizeOp
794 
795 
796 /// @brief Combines multiple grids into one by stealing leaf nodes and summing voxel values
797 /// This class is designed to work with thread local storage containers such as tbb::combinable
798 template<typename GridT>
799 struct GridCombinerOp
800 {
801  using CombinableT = typename tbb::combinable<GridT>;
802 
803  using TreeT = typename GridT::TreeType;
804  using LeafT = typename TreeT::LeafNodeType;
805  using ValueType = typename TreeT::ValueType;
806  using SumOpT = tools::valxform::SumOp<typename TreeT::ValueType>;
807  using MaxOpT = tools::valxform::MaxOp<typename TreeT::ValueType>;
808 
809  GridCombinerOp(GridT& grid, bool sum)
810  : mTree(grid.tree())
811  , mSum(sum) {}
812 
813  void operator()(const GridT& grid)
814  {
815  for (auto leaf = grid.tree().beginLeaf(); leaf; ++leaf) {
816  auto* newLeaf = mTree.probeLeaf(leaf->origin());
817  if (!newLeaf) {
818  // if the leaf doesn't yet exist in the new tree, steal it
819  auto& tree = const_cast<GridT&>(grid).tree();
820  mTree.addLeaf(tree.template stealNode<LeafT>(leaf->origin(),
821  zeroVal<ValueType>(), false));
822  }
823  else {
824  // otherwise increment existing values
825  if (mSum) {
826  for (auto iter = leaf->cbeginValueOn(); iter; ++iter) {
827  newLeaf->modifyValue(iter.offset(), SumOpT(ValueType(*iter)));
828  }
829  } else {
830  for (auto iter = leaf->cbeginValueOn(); iter; ++iter) {
831  newLeaf->modifyValue(iter.offset(), MaxOpT(ValueType(*iter)));
832  }
833  }
834  }
835  }
836  }
837 
838 private:
839  TreeT& mTree;
840  bool mSum;
841 }; // struct GridCombinerOp
842 
843 
844 template<typename GridT>
845 struct CombinableTraits {
846  using OpT = GridCombinerOp<GridT>;
847  using T = typename OpT::CombinableT;
848 };
849 
850 
851 template <typename PointDataGridT>
852 class GridToRasterize
853 {
854 public:
855  using GridPtr = typename PointDataGridT::Ptr;
856  using GridConstPtr = typename PointDataGridT::ConstPtr;
857  using PointDataTreeT = typename PointDataGridT::TreeType;
858  using PointDataLeafT = typename PointDataTreeT::LeafNodeType;
859 
860  GridToRasterize(GridPtr& grid, bool stream,
861  const FrustumRasterizerSettings& settings,
862  const FrustumRasterizerMask& mask)
863  : mGrid(grid)
864  , mStream(stream)
865  , mSettings(settings)
866  , mMask(mask) { }
867 
868  GridToRasterize(GridConstPtr& grid, const FrustumRasterizerSettings& settings,
869  const FrustumRasterizerMask& mask)
870  : mGrid(ConstPtrCast<PointDataGridT>(grid))
871  , mStream(false)
872  , mSettings(settings)
873  , mMask(mask) { }
874 
875  const typename PointDataGridT::TreeType& tree() const
876  {
877  return mGrid->constTree();
878  }
879 
880  void setLeafPercentage(int percentage)
881  {
882  mLeafPercentage = math::Clamp(percentage, 0, 100);
883  }
884 
885  int leafPercentage() const
886  {
887  return mLeafPercentage;
888  }
889 
890  size_t memUsage() const
891  {
892  return mGrid->memUsage() + mLeafOffsets.capacity();
893  }
894 
895  template <typename AttributeT, typename GridT, typename FilterT>
896  void rasterize(const Name& attribute,
897  typename CombinableTraits<GridT>::T& combiner, typename CombinableTraits<GridT>::T* weightCombiner,
898  const float scale, const FilterT& groupFilter, const bool computeMax, const bool reduceMemory,
899  util::NullInterrupter* interrupter)
900  {
901  using point_rasterize_internal::RasterizeOp;
902 
903  const auto& velocityAttribute = mSettings.velocityMotionBlur ?
904  mSettings.velocityAttribute : "";
905  const auto& radiusAttribute = mSettings.useRadius ?
906  mSettings.radiusAttribute : "";
907 
908  bool isPositionAttribute = attribute == "P";
909  bool isVelocityAttribute = attribute == mSettings.velocityAttribute;
910  bool isRadiusAttribute = attribute == mSettings.radiusAttribute;
911 
912  // find the attribute index
913  const auto& attributeSet = mGrid->constTree().cbeginLeaf()->attributeSet();
914  const size_t attributeIndex = attributeSet.find(attribute);
915  const bool attributeExists = attributeIndex != points::AttributeSet::INVALID_POS;
916 
917  if (attributeExists)
918  {
919  // throw if attribute type doesn't match AttributeT
920  const auto* attributeArray = attributeSet.getConst(attributeIndex);
921  const Index stride = bool(attributeArray) ? attributeArray->stride() : Index(1);
922  const auto& actualValueType = attributeSet.descriptor().valueType(attributeIndex);
923  const auto requestedValueType = Name(typeNameAsString<AttributeT>());
924  const bool packVector = stride == 3 &&
925  ((actualValueType == typeNameAsString<float>() &&
926  requestedValueType == typeNameAsString<Vec3f>()) ||
927  (actualValueType == typeNameAsString<double>() &&
928  requestedValueType == typeNameAsString<Vec3d>()) ||
929  (actualValueType == typeNameAsString<int32_t>() &&
930  requestedValueType == typeNameAsString<Vec3I>()));
931  if (!packVector && actualValueType != requestedValueType) {
932  OPENVDB_THROW(TypeError,
933  "Attribute type requested for rasterization \"" << requestedValueType << "\""
934  " must match attribute value type \"" << actualValueType << "\"");
935  }
936  }
937 
938  tree::LeafManager<PointDataTreeT> leafManager(mGrid->tree());
939 
940  // de-allocate voxel buffers if the points are being streamed and the caches are being released
941  const bool dropBuffers = mStream && reduceMemory;
942 
943  // turn RasterGroups into a MultiGroupFilter for this VDB Grid (if applicable)
944  using ResolvedFilterT = typename point_rasterize_internal::RasterGroupTraits<FilterT>::NewFilterT;
945  auto resolvedFilter = point_rasterize_internal::RasterGroupTraits<FilterT>::resolve(
946  groupFilter, attributeSet);
947 
948  // generate leaf offsets (if necessary)
949  if (mLeafOffsets.empty()) {
950  openvdb::points::pointOffsets(mLeafOffsets, mGrid->constTree(), resolvedFilter,
951  /*inCoreOnly=*/false, mSettings.threaded);
952  }
953 
954  // set streaming arbitrary attribute array flags
955  if (mStream) {
956  if (attributeExists && !isPositionAttribute && !isVelocityAttribute && !isRadiusAttribute) {
957  leafManager.foreach(
958  [&](PointDataLeafT& leaf, size_t /*idx*/) {
959  leaf.attributeArray(attributeIndex).setStreaming(true);
960  },
961  mSettings.threaded);
962  }
963 
964  if (reduceMemory) {
965  size_t positionIndex = attributeSet.find("P");
966  leafManager.foreach(
967  [&](PointDataLeafT& leaf, size_t /*idx*/) {
968  leaf.attributeArray(positionIndex).setStreaming(true);
969  },
970  mSettings.threaded);
971  if (mSettings.velocityMotionBlur || !mSettings.camera.isStatic()) {
972  size_t velocityIndex = attributeSet.find(velocityAttribute);
973  if (velocityIndex != points::AttributeSet::INVALID_POS) {
974  leafManager.foreach(
975  [&](PointDataLeafT& leaf, size_t /*idx*/) {
976  leaf.attributeArray(velocityIndex).setStreaming(true);
977  },
978  mSettings.threaded);
979  }
980  }
981  if (mSettings.useRadius) {
982  size_t radiusIndex = attributeSet.find(radiusAttribute);
983  if (radiusIndex != points::AttributeSet::INVALID_POS) {
984  leafManager.foreach(
985  [&](PointDataLeafT& leaf, size_t /*idx*/) {
986  leaf.attributeArray(radiusIndex).setStreaming(true);
987  },
988  mSettings.threaded);
989  }
990  }
991  }
992  }
993 
994  const bool alignedTransform = *(mSettings.transform) == mGrid->constTransform();
995 
996  RasterizeOp<PointDataGridT, AttributeT, GridT, ResolvedFilterT> rasterizeOp(
997  *mGrid, mLeafOffsets, attributeIndex, velocityAttribute, radiusAttribute, combiner, weightCombiner,
998  dropBuffers, scale, resolvedFilter, computeMax, alignedTransform, mSettings.camera.isStatic(),
999  mSettings, mMask, interrupter);
1000  leafManager.foreach(rasterizeOp, mSettings.threaded);
1001 
1002  // clean up leaf offsets (if reduce memory is enabled)
1003  if (reduceMemory && !mLeafOffsets.empty()) {
1004  mLeafOffsets.clear();
1005  // de-allocate the vector data to ensure we keep memory footprint as low as possible
1006  mLeafOffsets.shrink_to_fit();
1007  }
1008  }
1009 
1010 private:
1011  GridPtr mGrid;
1012  const bool mStream;
1013  const FrustumRasterizerSettings& mSettings;
1014  const FrustumRasterizerMask& mMask;
1015  int mLeafPercentage = -1;
1016  std::vector<Index64> mLeafOffsets;
1017 }; // class GridToRasterize
1018 
1019 
1020 template <typename ValueT>
1021 typename std::enable_if<!ValueTraits<ValueT>::IsVec, ValueT>::type
1022 computeWeightedValue(const ValueT& value, const ValueT& weight)
1023 {
1024  constexpr bool isSignedInt = std::is_integral<ValueT>() && std::is_signed<ValueT>();
1025 
1026  if (!math::isFinite(weight) || math::isApproxZero(weight) ||
1027  (isSignedInt && math::isNegative(weight))) {
1028  return zeroVal<ValueT>();
1029  } else {
1030  return value / weight;
1031  }
1032 }
1033 
1034 
1035 template <typename ValueT>
1036 typename std::enable_if<ValueTraits<ValueT>::IsVec, ValueT>::type
1037 computeWeightedValue(const ValueT& value, const ValueT& weight)
1038 {
1039  using ElementT = typename ValueTraits<ValueT>::ElementType;
1040 
1041  constexpr bool isSignedInt = std::is_integral<ElementT>() && std::is_signed<ElementT>();
1042 
1043  ValueT result(value);
1044  for (int i=0; i<ValueTraits<ValueT>::Size; ++i) {
1045  if (!math::isFinite(weight[i]) || math::isApproxZero(weight[i]) ||
1046  (isSignedInt && math::isNegative(weight[i]))) {
1047  result[i] = zeroVal<ElementT>();
1048  } else {
1049  result[i] = value[i] / weight[i];
1050  }
1051  }
1052  return result;
1053 }
1054 
1055 
1056 } // namespace point_rasterize_internal
1057 
1058 
1059 ////////////////////////////////////////////////////////////////////////////
1060 
1061 
1062 RasterCamera::RasterCamera(const math::Transform& transform)
1063  : mTransforms{transform}
1064  , mWeights{1} { }
1065 
1066 bool RasterCamera::isStatic() const
1067 {
1068  return mTransforms.size() <= 1;
1069 }
1070 
1071 void RasterCamera::clear()
1072 {
1073  mTransforms.clear();
1074  mWeights.clear();
1075 }
1076 
1077 void RasterCamera::appendTransform(const math::Transform& transform, float weight)
1078 {
1079  mTransforms.push_back(transform);
1080  mWeights.push_back(weight);
1081 }
1082 
1083 size_t RasterCamera::size() const
1084 {
1085  return mTransforms.size();
1086 }
1087 
1089 {
1090  // if two or more identical transforms, only keep one
1091  if (mTransforms.size() >= 2) {
1092  const auto& transform = mTransforms.front();
1093  bool isStatic = true;
1094  for (const auto& testTransform : mTransforms) {
1095  if (transform != testTransform) {
1096  isStatic = false;
1097  }
1098  }
1099  if (isStatic) {
1100  while (mTransforms.size() > 1) {
1101  mTransforms.pop_back();
1102  }
1103  }
1104  }
1105  // if all weights are equal to one, delete the weights array
1106  if (!mWeights.empty()) {
1107  bool hasWeight = false;
1108  for (Index i = 0; i < mWeights.size(); i++) {
1109  if (this->hasWeight(i)) {
1110  hasWeight = true;
1111  break;
1112  }
1113  }
1114  if (!hasWeight) mWeights.clear();
1115  }
1116 }
1117 
1118 bool RasterCamera::hasWeight(Index i) const
1119 {
1120  if (mWeights.empty()) return false;
1121  assert(i < mWeights.size());
1122  return !openvdb::math::isApproxEqual(mWeights[i], 1.0f, 1e-3f);
1123 }
1124 
1125 float RasterCamera::weight(Index i) const
1126 {
1127  if (mWeights.empty()) {
1128  return 1.0f;
1129  } else {
1130  assert(i < mWeights.size());
1131  return mWeights[i];
1132  }
1133 }
1134 
1135 const math::Transform& RasterCamera::transform(Index i) const
1136 {
1137  if (mTransforms.size() == 1) {
1138  return mTransforms.front();
1139  } else {
1140  assert(i < mTransforms.size());
1141  return mTransforms[i];
1142  }
1143 }
1144 
1145 const math::Transform& RasterCamera::firstTransform() const
1146 {
1147  assert(!mTransforms.empty());
1148  return mTransforms.front();
1149 }
1150 
1151 const math::Transform& RasterCamera::lastTransform() const
1152 {
1153  assert(!mTransforms.empty());
1154  return mTransforms.back();
1155 }
1156 
1157 void RasterCamera::setShutter(float start, float end)
1158 {
1159  mShutterStart = start;
1160  mShutterEnd = end;
1161 }
1162 
1163 float RasterCamera::shutterStart() const
1164 {
1165  return mShutterStart;
1166 }
1167 
1168 float RasterCamera::shutterEnd() const
1169 {
1170  return mShutterEnd;
1171 }
1172 
1173 
1174 ////////////////////////////////////////////////////////////////////////////
1175 
1176 
1177 FrustumRasterizerMask::FrustumRasterizerMask( const math::Transform& transform,
1178  const MaskGrid* mask,
1179  const BBoxd& bbox,
1180  const bool clipToFrustum,
1181  const bool invertMask)
1182  : mMask()
1183  , mClipBBox()
1184  , mInvert(invertMask)
1185 {
1186  // TODO: the current OpenVDB implementation for resampling masks is particularly slow,
1187  // this is primarily because it uses a scatter reduction-style method that only samples
1188  // and generates leaf nodes and relies on pruning to generate tiles, this could be
1189  // significantly improved!
1190 
1191  // convert world-space clip mask to index-space
1192  if (mask) {
1193  // resample mask to index space
1194  mMask = MaskGrid::create();
1195  mMask->setTransform(transform.copy());
1196  tools::resampleToMatch<tools::PointSampler>(*mask, *mMask);
1197 
1198  // prune the clip mask
1199  tools::prune(mMask->tree());
1200  }
1201 
1202  // convert world-space clip bbox to index-space
1203  if (!bbox.empty()) {
1204  // create world-space mask (with identity linear transform)
1205  MaskGrid::Ptr tempMask = MaskGrid::create();
1206  CoordBBox coordBBox(Coord::floor(bbox.min()),
1207  Coord::ceil(bbox.max()));
1208  tempMask->sparseFill(coordBBox, true, true);
1209 
1210  // resample mask to index space
1211  MaskGrid::Ptr bboxMask = MaskGrid::create();
1212  bboxMask->setTransform(mMask ? mMask->transformPtr() : transform.copy());
1213  tools::resampleToMatch<tools::PointSampler>(*tempMask, *bboxMask);
1214 
1215  // replace or union the mask
1216  if (mMask) {
1217  mMask->topologyUnion(*bboxMask);
1218  } else {
1219  mMask = bboxMask;
1220  }
1221  }
1222 
1223  if (clipToFrustum) {
1224  auto frustumMap = transform.template constMap<math::NonlinearFrustumMap>();
1225  if (frustumMap) {
1226  const auto& frustumBBox = frustumMap->getBBox();
1227  mClipBBox.reset(Coord::floor(frustumBBox.min()),
1228  Coord::ceil(frustumBBox.max()));
1229  }
1230  }
1231 }
1232 
1233 FrustumRasterizerMask::operator bool() const
1234 {
1235  return mMask || !mClipBBox.empty();
1236 }
1237 
1240 {
1241  return mMask ? mMask->treePtr() : MaskTree::ConstPtr();
1242 }
1243 
1244 const CoordBBox&
1246 {
1247  return mClipBBox;
1248 }
1249 
1250 bool
1251 FrustumRasterizerMask::valid(const Coord& ijk, const tree::ValueAccessor<const MaskTree>* acc) const
1252 {
1253  const bool maskValue = acc ? acc->isValueOn(ijk) : true;
1254  const bool insideMask = mInvert ? !maskValue : maskValue;
1255  const bool insideFrustum = mClipBBox.empty() ? true : mClipBBox.isInside(ijk);
1256  return insideMask && insideFrustum;
1257 }
1258 
1259 
1260 ////////////////////////////////////////////////////////////////////////////
1261 
1262 
1263 template <typename PointDataGridT>
1264 FrustumRasterizer<PointDataGridT>::FrustumRasterizer(const FrustumRasterizerSettings& settings,
1265  const FrustumRasterizerMask& mask,
1266  util::NullInterrupter* interrupt)
1267  : mSettings(settings)
1268  , mMask(mask)
1269  , mInterrupter(interrupt)
1270 {
1271  if (mSettings.velocityAttribute.empty() && mSettings.velocityMotionBlur) {
1272  OPENVDB_THROW(ValueError,
1273  "Using velocity motion blur during rasterization requires a velocity attribute.");
1274  }
1275 }
1276 
1277 template <typename PointDataGridT>
1278 void
1280 {
1281  // skip any empty grids
1282  if (!grid || grid->tree().empty()) return;
1283 
1284  // note that streaming is not possible with a const grid
1285  mPointGrids.emplace_back(grid, mSettings, mMask);
1286 }
1287 
1288 template <typename PointDataGridT>
1289 void
1291 {
1292  // skip any empty grids
1293  if (!grid || grid->tree().empty()) return;
1294 
1295  mPointGrids.emplace_back(grid, stream, mSettings, mMask);
1296 }
1297 
1298 template <typename PointDataGridT>
1299 void
1301 {
1302  mPointGrids.clear();
1303 }
1304 
1305 template <typename PointDataGridT>
1306 size_t
1308 {
1309  return mPointGrids.size();
1310 }
1311 
1312 template <typename PointDataGridT>
1313 size_t
1315 {
1316  size_t mem = sizeof(*this) + sizeof(mPointGrids);
1317  for (const auto& grid : mPointGrids) {
1318  mem += grid.memUsage();
1319  }
1320  return mem;
1321 }
1322 
1323 template <typename PointDataGridT>
1324 template <typename FilterT>
1327  RasterMode mode, bool reduceMemory, float scale, const FilterT& filter)
1328 {
1329  // no attribute to rasterize, so just provide an empty string and default to float type
1330  auto density = rasterizeAttribute<FloatGrid, float>("", mode, reduceMemory, scale, filter);
1331  // hard-code grid name to density
1332  density->setName("density");
1333  return density;
1334 }
1335 
1336 template <typename PointDataGridT>
1337 template <typename FilterT>
1340  const openvdb::Name& attribute, RasterMode mode, bool reduceMemory, float scale, const FilterT& filter)
1341 {
1342  auto density = rasterizeAttribute<FloatGrid, float>(attribute, mode, reduceMemory, scale, filter);
1343  // hard-code grid name to density
1344  density->setName("density");
1345  return density;
1346 }
1347 
1348 template <typename PointDataGridT>
1349 template <typename FilterT>
1352  const Name& attribute, RasterMode mode, bool reduceMemory, float scale, const FilterT& filter)
1353 {
1354  // retrieve the source type of the attribute
1355 
1356  Name sourceType;
1357  Index stride(0);
1358  for (const auto& points : mPointGrids) {
1359  auto leaf = points.tree().cbeginLeaf();
1360  const auto& descriptor = leaf->attributeSet().descriptor();
1361  const size_t targetIndex = descriptor.find(attribute);
1362  // ignore grids which don't have the attribute
1363  if (targetIndex == points::AttributeSet::INVALID_POS) continue;
1364  const auto* attributeArray = leaf->attributeSet().getConst(attribute);
1365  if (!attributeArray) continue;
1366  stride = attributeArray->stride();
1367  sourceType = descriptor.valueType(targetIndex);
1368  if (!sourceType.empty()) break;
1369  }
1370 
1371  // no valid point attributes for rasterization
1372  // TODO: add a warning / error in the case that there are non-zero grids
1373  if (sourceType.empty()) return {};
1374 
1375  if (stride == 1 && sourceType == typeNameAsString<float>()) {
1376  using GridT = typename PointDataGridT::template ValueConverter<float>::Type;
1377  return rasterizeAttribute<GridT, float>(attribute, mode, reduceMemory, scale, filter);
1378 // this define is for lowering compilation time when debugging by instantiating only the float
1379 // code path (default is to instantiate all code paths)
1380 #ifndef ONLY_RASTER_FLOAT
1381  } else if ( sourceType == typeNameAsString<Vec3f>() ||
1382  (stride == 3 && sourceType == typeNameAsString<float>())) {
1383  using GridT = typename PointDataGridT::template ValueConverter<Vec3f>::Type;
1384  return rasterizeAttribute<GridT, Vec3f>(attribute, mode, reduceMemory, scale, filter);
1385  } else if ( sourceType == typeNameAsString<Vec3d>() ||
1386  (stride == 3 && sourceType == typeNameAsString<double>())) {
1387  using GridT = typename PointDataGridT::template ValueConverter<Vec3d>::Type;
1388  return rasterizeAttribute<GridT, Vec3d>(attribute, mode, reduceMemory, scale, filter);
1389  } else if ( sourceType == typeNameAsString<Vec3i>() ||
1390  (stride == 3 && sourceType == typeNameAsString<int32_t>())) {
1391  using GridT = typename PointDataGridT::template ValueConverter<Vec3i>::Type;
1392  return rasterizeAttribute<GridT, Vec3i>(attribute, mode, reduceMemory, scale, filter);
1393  } else if (stride == 1 && sourceType == typeNameAsString<int16_t>()) {
1394  using GridT = typename PointDataGridT::template ValueConverter<Int32>::Type;
1395  return rasterizeAttribute<GridT, int16_t>(attribute, mode, reduceMemory, scale, filter);
1396  } else if (stride == 1 && sourceType == typeNameAsString<int32_t>()) {
1397  using GridT = typename PointDataGridT::template ValueConverter<Int32>::Type;
1398  return rasterizeAttribute<GridT, int32_t>(attribute, mode, reduceMemory, scale, filter);
1399  } else if (stride == 1 && sourceType == typeNameAsString<int64_t>()) {
1400  using GridT = typename PointDataGridT::template ValueConverter<Int64>::Type;
1401  return rasterizeAttribute<GridT, int64_t>(attribute, mode, reduceMemory, scale, filter);
1402  } else if (stride == 1 && sourceType == typeNameAsString<double>()) {
1403  using GridT = typename PointDataGridT::template ValueConverter<double>::Type;
1404  return rasterizeAttribute<GridT, double>(attribute, mode, reduceMemory, scale, filter);
1405  } else if (stride == 1 && sourceType == typeNameAsString<bool>()) {
1406  using GridT = typename PointDataGridT::template ValueConverter<bool>::Type;
1407  return rasterizeAttribute<GridT, bool>(attribute, mode, reduceMemory, true, filter);
1408 #endif
1409  } else {
1410  std::ostringstream ostr;
1411  ostr << "Cannot rasterize attribute of type - " << sourceType;
1412  if (stride > 1) ostr << " x " << stride;
1413  OPENVDB_THROW(TypeError, ostr.str());
1414  }
1415 }
1416 
1417 template <typename PointDataGridT>
1418 template <typename GridT, typename AttributeT, typename FilterT>
1419 typename GridT::Ptr
1421  bool reduceMemory, float scale, const FilterT& filter)
1422 {
1423  if (attribute == "P") {
1424  OPENVDB_THROW(ValueError, "Cannot rasterize position attribute.")
1425  }
1426 
1427  auto grid = GridT::create();
1428  grid->setName(attribute);
1429  grid->setTransform(mSettings.transform->copy());
1430 
1431  this->performRasterization<AttributeT>(*grid, mode, attribute, reduceMemory, scale, filter);
1432 
1433  return grid;
1434 }
1435 
1436 template <typename PointDataGridT>
1437 template <typename GridT, typename FilterT>
1438 typename GridT::Ptr
1439 FrustumRasterizer<PointDataGridT>::rasterizeMask(bool reduceMemory, const FilterT& filter)
1440 {
1441  using ValueT = typename GridT::ValueType;
1442 
1443  static_assert(point_rasterize_internal::BoolTraits<ValueT>::IsBool,
1444  "Value type of mask to be rasterized must be bool or ValueMask.");
1445 
1446  auto grid = rasterizeAttribute<GridT, ValueT>("", RasterMode::ACCUMULATE, reduceMemory, true, filter);
1447  grid->setName("mask");
1448 
1449  return grid;
1450 }
1451 
1452 template <typename PointDataGridT>
1453 template <typename AttributeT, typename GridT, typename FilterT>
1454 void
1455 FrustumRasterizer<PointDataGridT>::performRasterization(
1456  GridT& grid, RasterMode mode, const openvdb::Name& attribute, bool reduceMemory,
1457  float scale, const FilterT& filter)
1458 {
1459  using openvdb::points::point_mask_internal::GridCombinerOp;
1460  using point_rasterize_internal::computeWeightedValue;
1461 
1462  using TreeT = typename GridT::TreeType;
1463  using LeafT = typename TreeT::LeafNodeType;
1464  using ValueT = typename TreeT::ValueType;
1465  using CombinerOpT = typename point_rasterize_internal::CombinableTraits<GridT>::OpT;
1466  using CombinableT = typename point_rasterize_internal::CombinableTraits<GridT>::T;
1467 
1468  // start interrupter with a descriptive name
1469  if (mInterrupter) {
1470  std::stringstream ss;
1471  ss << "Rasterizing Points ";
1472  if (!mSettings.camera.isStatic() || mSettings.velocityMotionBlur) {
1473  ss << "with Motion Blur ";
1474  }
1475  if (mSettings.transform->isLinear()) {
1476  ss << "to Linear Volume";
1477  } else {
1478  ss << "to Frustum Volume";
1479  }
1480  mInterrupter->start(ss.str().c_str());
1481  }
1482 
1483  const bool useMaximum = mode == RasterMode::MAXIMUM;
1484  const bool useWeight = mode == RasterMode::AVERAGE;
1485 
1486  if (useMaximum && VecTraits<ValueT>::IsVec) {
1487  OPENVDB_THROW(ValueError,
1488  "Cannot use maximum mode when rasterizing vector attributes.");
1489  }
1490 
1491  if ((useMaximum || useWeight) && point_rasterize_internal::BoolTraits<ValueT>::IsBool) {
1492  OPENVDB_THROW(ValueError,
1493  "Cannot use maximum or average modes when rasterizing bool attributes.");
1494  }
1495 
1496  CombinableT combiner;
1497  CombinableT weightCombiner;
1498  CombinableT* weightCombinerPtr = useWeight ? &weightCombiner : nullptr;
1499 
1500  // use leaf count as an approximate indicator for progress as it
1501  // doesn't even require loading the topology
1502 
1503  if (mInterrupter) {
1504  if (mPointGrids.size() == 1) {
1505  mPointGrids[0].setLeafPercentage(100);
1506  }
1507  else {
1508  // compute cumulative leaf counts per grid
1509  Index64 leafCount(0);
1510  std::vector<Index64> leafCounts;
1511  leafCounts.reserve(mPointGrids.size());
1512  for (auto& points : mPointGrids) {
1513  leafCount += points.tree().leafCount();
1514  leafCounts.push_back(leafCount);
1515  }
1516 
1517  // avoid dividing by zero
1518  if (leafCount == Index64(0)) leafCount++;
1519 
1520  // compute grid percentages based on leaf count
1521  for (size_t i = 0; i < leafCounts.size(); i++) {
1522  int percentage = static_cast<int>(math::Round((static_cast<float>(leafCounts[i]))/
1523  static_cast<float>(leafCount)));
1524  mPointGrids[i].setLeafPercentage(percentage);
1525  }
1526  }
1527  }
1528 
1529  // rasterize each point grid into each grid in turn
1530 
1531  for (auto& points : mPointGrids) {
1532 
1533  points.template rasterize<AttributeT, GridT>(
1534  attribute, combiner, weightCombinerPtr, scale, filter, mode == RasterMode::MAXIMUM, reduceMemory, mInterrupter);
1535 
1536  // interrupt if requested and update the progress percentage
1537  // note that even when interrupting, the operation to combine the local grids
1538  // is completed so that the user receives a partially rasterized volume
1539 
1540  if (mInterrupter &&
1541  mInterrupter->wasInterrupted(points.leafPercentage())) {
1542  break;
1543  }
1544  }
1545 
1546  // combine the value grids into one
1547 
1548  CombinerOpT combineOp(grid, mode != RasterMode::MAXIMUM);
1549  combiner.combine_each(combineOp);
1550 
1551  if (useWeight) {
1552 
1553  // combine the weight grids into one
1554 
1555  auto weightGrid = GridT::create(ValueT(1));
1556  CombinerOpT weightCombineOp(*weightGrid, /*sum=*/true);
1557  weightCombiner.combine_each(weightCombineOp);
1558 
1559  tree::LeafManager<TreeT> leafManager(grid.tree());
1560  leafManager.foreach(
1561  [&](LeafT& leaf, size_t /*idx*/) {
1562  auto weightAccessor = weightGrid->getConstAccessor();
1563  for (auto iter = leaf.beginValueOn(); iter; ++iter) {
1564  auto weight = weightAccessor.getValue(iter.getCoord());
1565  iter.setValue(computeWeightedValue(iter.getValue(), weight));
1566  }
1567  },
1568  mSettings.threaded);
1569  }
1570 
1571  if (mInterrupter) mInterrupter->end();
1572 }
1573 
1574 
1575 } // namespace points
1576 } // namespace OPENVDB_VERSION_NAME
1577 } // namespace openvdb
1578 
1579 #endif // OPENVDB_POINTS_POINT_RASTERIZE_FRUSTUM_IMPL_HAS_BEEN_INCLUDED
GLfloat GLfloat GLfloat top
Definition: glew.h:15525
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
Definition: ImfName.h:30
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:859
GLenum GLenum GLenum GLenum GLenum scale
Definition: glew.h:14163
math::Vec3< Real > Vec3R
Definition: Types.h:72
SharedPtr< GridBase > Ptr
Definition: Grid.h:80
void rasterize(const PointDataTreeOrGridT &points, TransferT &transfer, const FilterT &filter=NullFilter(), InterrupterT *interrupter=nullptr)
Perform potentially complex rasterization from a user defined transfer scheme.
RasterMode
How to composite points into a volume.
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
IMATH_HOSTDEVICE constexpr int floor(T x) IMATH_NOEXCEPT
Definition: ImathFun.h:112
GLint left
Definition: glcorearb.h:2005
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
const char * typeNameAsString< int32_t >()
Definition: Types.h:487
GLuint start
Definition: glcorearb.h:475
const char * typeNameAsString< Vec3i >()
Definition: Types.h:495
Index64 memUsage(const TreeT &tree, bool threaded=true)
Return the total amount of memory in bytes occupied by this tree.
Definition: Count.h:493
GLuint GLsizei const GLuint const GLintptr * offsets
Definition: glcorearb.h:2621
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7481
const char * typeNameAsString< int16_t >()
Definition: Types.h:485
GLuint GLuint stream
Definition: glcorearb.h:1832
size_t size() const
Return number of PointDataGrids in the rasterizer.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:207
FrustumRasterizer(const FrustumRasterizerSettings &settings, const FrustumRasterizerMask &mask=FrustumRasterizerMask(), util::NullInterrupter *interrupt=nullptr)
main constructor
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:195
bool isNegative(const Type &x)
Return true if x is less than zero.
Definition: Math.h:368
GLfloat right
Definition: glew.h:15525
const char * typeNameAsString< Vec3d >()
Definition: Types.h:497
FloatGrid::Ptr rasterizeUniformDensity(RasterMode mode=RasterMode::MAXIMUM, bool reduceMemory=false, float scale=1.0f, const FilterT &filter=FilterT())
SharedPtr< T > ConstPtrCast(const SharedPtr< U > &ptr)
Return a new shared pointer that points to the same object as the given pointer but with possibly dif...
Definition: Types.h:126
GLuint GLuint GLfloat weight
Definition: glew.h:13892
GLint GLenum GLint x
Definition: glcorearb.h:409
GLuint GLenum GLenum transform
Definition: glew.h:15055
const char * typeNameAsString< Vec3f >()
Definition: Types.h:496
GridBase::Ptr rasterizeAttribute(const Name &attribute, RasterMode mode=RasterMode::ACCUMULATE, bool reduceMemory=false, float scale=1.0f, const FilterT &filter=FilterT())
GLuint64EXT * result
Definition: glew.h:14311
bool valid(const Coord &ijk, AccessorT *acc) const
void appendTransform(const math::Transform &, float weight=1.0f)
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition: Math.h:407
static Ptr create()
Return a new grid with background value zero.
Definition: Grid.h:1313
GLsizei GLsizei GLfloat distance
Definition: glew.h:13923
GLint GLuint mask
Definition: glcorearb.h:124
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:822
const char * typeNameAsString< float >()
Definition: Types.h:481
size_t memUsage() const
Return memory usage of the rasterizer.
const GLdouble * v
Definition: glcorearb.h:837
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
Grid< MaskTree > MaskGrid
Definition: openvdb.h:78
GLuint GLuint end
Definition: glcorearb.h:475
GLint GLint bottom
Definition: glcorearb.h:2005
RasterCamera(const math::Transform &transform)
OPENVDB_API void calculateBounds(const Transform &t, const Vec3d &minWS, const Vec3d &maxWS, Vec3d &minIS, Vec3d &maxIS)
Calculate an axis-aligned bounding box in index space from an axis-aligned bounding box in world spac...
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
void addPoints(GridConstPtr &points)
Append a PointDataGrid to the rasterizer (but don't rasterize yet).
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:198
const char * typeNameAsString< int64_t >()
Definition: Types.h:489
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition: glcorearb.h:1297
SharedPtr< const Tree > ConstPtr
Definition: Tree.h:181
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:350
vfloat4 round(const vfloat4 &a)
Definition: simd.h:7436
GLint GLsizei width
Definition: glcorearb.h:103
GLint GLsizei GLsizei height
Definition: glcorearb.h:103
GridT::Ptr rasterizeMask(bool reduceMemory=false, const FilterT &filter=FilterT())
const char * typeNameAsString< double >()
Definition: Types.h:482
GLenum mode
Definition: glcorearb.h:99
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:260
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glew.h:3460
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
IMATH_HOSTDEVICE constexpr int ceil(T x) IMATH_NOEXCEPT
Definition: ImathFun.h:119
static Ptr create(const AttributeArray &array, const bool collapseOnDestruction=true)
const char * typeNameAsString< bool >()
Definition: Types.h:478
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:89
GLuint index
Definition: glcorearb.h:786
const math::Transform & firstTransform() const
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
math::BBox< Vec3d > BBoxd
Definition: Types.h:84
Index64 pointOffsets(std::vector< Index64 > &pointOffsets, const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Populate an array of cumulative point offsets per leaf node.
Definition: PointCount.h:122
GLsizei const GLfloat * value
Definition: glcorearb.h:824
OIIO_API bool attribute(string_view name, TypeDesc type, const void *val)
GLfloat f
Definition: glcorearb.h:1926
int Floor(float x)
Return the floor of x.
Definition: Math.h:851
Definition: core.h:1131
GLintptr offset
Definition: glcorearb.h:665
GLint GLenum GLboolean GLsizei stride
Definition: glcorearb.h:872
void prune(TreeT &tree, typename TreeT::ValueType tolerance=zeroVal< typename TreeT::ValueType >(), bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with tiles any nodes whose values are all the same...
Definition: Prune.h:335
type
Definition: core.h:1059
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
GLsizei GLuint * groups
Definition: glew.h:2749
const math::Transform & transform(Index i) const
GLint y
Definition: glcorearb.h:103
void clear()
Clear all PointDataGrids in the rasterizer.
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
bool isFinite(const float x)
Return true if x is finite.
Definition: Math.h:376
const math::Transform & lastTransform() const
ImageBuf OIIO_API mul(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
FloatGrid::Ptr rasterizeDensity(const openvdb::Name &attribute, RasterMode mode=RasterMode::MAXIMUM, bool reduceMemory=false, float scale=1.0f, const FilterT &filter=FilterT())
math::Vec3< float > Vec3f
Definition: Types.h:74