OpenVDB  12.0.0
PointRasterizeTrilinearImpl.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 //
4 /// @author Nick Avramoussis
5 ///
6 /// @file PointRasterizeTrilinearImpl.h
7 ///
8 
9 #ifndef OPENVDB_POINTS_RASTERIZE_TRILINEAR_IMPL_HAS_BEEN_INCLUDED
10 #define OPENVDB_POINTS_RASTERIZE_TRILINEAR_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 rasterize_trilinear_internal {
20 
21 template <typename TreeType,
22  typename PositionCodecT,
23  typename SourceValueT,
24  typename SourceCodecT>
25 struct TrilinearTransfer : public VolumeTransfer<TreeType>
26 {
27  using BaseT = VolumeTransfer<TreeType>;
28  using WeightT = typename TreeType::ValueType;
29  using PositionHandleT = points::AttributeHandle<Vec3f, PositionCodecT>;
30  using SourceHandleT = points::AttributeHandle<SourceValueT, SourceCodecT>;
31 
32  // precision of kernel arithmetic - aliased to the floating point
33  // element type of the source attribute.
34  using SourceElementT = typename ValueTraits<SourceValueT>::ElementType;
35  using RealT = SourceElementT;
36 
37  static_assert(std::is_floating_point<SourceElementT>::value,
38  "Trilinear rasterization only supports floating point values.");
39 
40  static const Index NUM_VALUES = TreeType::LeafNodeType::NUM_VALUES;
41 
42  TrilinearTransfer(const size_t pidx,
43  const size_t sidx, TreeType& tree)
44  : BaseT(tree)
45  , mPIdx(pidx)
46  , mSIdx(sidx)
47  , mPHandle()
48  , mSHandle()
49  , mWeights() {}
50 
51  TrilinearTransfer(const TrilinearTransfer& other)
52  : BaseT(other)
53  , mPIdx(other.mPIdx)
54  , mSIdx(other.mSIdx)
55  , mPHandle()
56  , mSHandle()
57  , mWeights() {}
58 
59  //// @note Kernel value evaluator
60  static inline RealT value(const RealT x)
61  {
62  const RealT abs_x = std::fabs(x);
63  if (abs_x < RealT(1.0)) return RealT(1.0) - abs_x;
64  return RealT(0.0);
65  }
66 
67  inline static Int32 range() { return 1; }
68  inline Int32 range(const Coord&, size_t) const { return this->range(); }
69 
70  inline void initialize(const Coord& origin, const size_t idx, const CoordBBox& bounds)
71  {
72  this->BaseT::initialize(origin, idx, bounds);
73  mWeights.fill(openvdb::zeroVal<WeightT>());
74  }
75 
76  inline bool startPointLeaf(const PointDataTree::LeafNodeType& leaf)
77  {
78  mPHandle.reset(new PositionHandleT(leaf.constAttributeArray(mPIdx)));
79  mSHandle.reset(new SourceHandleT(leaf.constAttributeArray(mSIdx)));
80  return true;
81  }
82 
83  inline bool endPointLeaf(const PointDataTree::LeafNodeType&) { return true; }
84 
85 protected:
86  const size_t mPIdx;
87  const size_t mSIdx;
88  typename PositionHandleT::UniquePtr mPHandle;
89  typename SourceHandleT::UniquePtr mSHandle;
90  std::array<WeightT, NUM_VALUES> mWeights;
91 };
92 
93 template <typename TreeType,
94  typename PositionCodecT,
95  typename SourceValueT,
96  typename SourceCodecT>
97 struct StaggeredTransfer :
98  public TrilinearTransfer<TreeType, PositionCodecT, SourceValueT, SourceCodecT>
99 {
100  using BaseT = TrilinearTransfer<TreeType, PositionCodecT, SourceValueT, SourceCodecT>;
101  using RealT = typename BaseT::RealT;
102  using BaseT::value;
103 
105  "Target Tree must be a vector tree for staggered rasterization");
106 
107  static const Index DIM = TreeType::LeafNodeType::DIM;
108  static const Index LOG2DIM = TreeType::LeafNodeType::LOG2DIM;
109 
110  StaggeredTransfer(const size_t pidx,
111  const size_t sidx, TreeType& tree)
112  : BaseT(pidx, sidx, tree) {}
113 
114  void rasterizePoint(const Coord& ijk,
115  const Index id,
116  const CoordBBox& bounds)
117  {
118  CoordBBox intersectBox(ijk.offsetBy(-1), ijk.offsetBy(1));
119  intersectBox.intersect(bounds);
120  if (intersectBox.empty()) return;
121 
122  auto* const data = this->buffer();
123  const auto& mask = *(this->mask());
124 
125  const math::Vec3<RealT> P(this->mPHandle->get(id));
126  const SourceValueT s(this->mSHandle->get(id));
127 
128  math::Vec3<RealT> centerw, macw;
129 
130  const Coord& a(intersectBox.min());
131  const Coord& b(intersectBox.max());
132  for (Coord c = a; c.x() <= b.x(); ++c.x()) {
133  // @todo can probably simplify the double call to value() in some way
134  const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM); // unsigned bit shift mult
135  const RealT x = static_cast<RealT>(c.x()-ijk.x()); // distance from ijk to c
136  centerw[0] = value(P.x() - x); // center dist
137  macw.x() = value(P.x() - (x-RealT(0.5))); // mac dist
138 
139  for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
140  const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
141  const RealT y = static_cast<RealT>(c.y()-ijk.y());
142  centerw[1] = value(P.y() - y);
143  macw.y() = value(P.y() - (y-RealT(0.5)));
144 
145  for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
146  OPENVDB_ASSERT(bounds.isInside(c));
147  const Index offset = ij + /*k*/(c.z() & (DIM-1u));
148  if (!mask.isOn(offset)) continue;
149  const RealT z = static_cast<RealT>(c.z()-ijk.z());
150  centerw[2] = value(P.z() - z);
151  macw.z() = value(P.z() - (z-RealT(0.5)));
152 
153  const math::Vec3<RealT> r {
154  (macw[0] * centerw[1] * centerw[2]),
155  (macw[1] * centerw[0] * centerw[2]),
156  (macw[2] * centerw[0] * centerw[1])
157  };
158 
159  data[offset] += s * r;
160  this->mWeights[offset] += r;
161  }
162  }
163  }
164  }
165 
166  inline bool finalize(const Coord&, const size_t)
167  {
168  auto* const data = this->buffer();
169  const auto& mask = *(this->mask());
170 
171  for (auto iter = mask.beginOn(); iter; ++iter) {
172  const Index offset = iter.pos();
173  const auto& w = this->mWeights[offset];
174  auto& v = data[offset];
175  if (!math::isZero(w[0])) v[0] /= w[0];
176  if (!math::isZero(w[1])) v[1] /= w[1];
177  if (!math::isZero(w[2])) v[2] /= w[2];
178  }
179 
180  return true;
181  }
182 };
183 
184 template <typename TreeType,
185  typename PositionCodecT,
186  typename SourceValueT,
187  typename SourceCodecT>
188 struct CellCenteredTransfer :
189  public TrilinearTransfer<TreeType, PositionCodecT, SourceValueT, SourceCodecT>
190 {
191  using BaseT = TrilinearTransfer<TreeType, PositionCodecT, SourceValueT, SourceCodecT>;
192  using RealT = typename BaseT::RealT;
193  using BaseT::value;
194 
195  static const Index DIM = TreeType::LeafNodeType::DIM;
196  static const Index LOG2DIM = TreeType::LeafNodeType::LOG2DIM;
197 
198  CellCenteredTransfer(const size_t pidx,
199  const size_t sidx, TreeType& tree)
200  : BaseT(pidx, sidx, tree) {}
201 
202  void rasterizePoint(const Coord& ijk,
203  const Index id,
204  const CoordBBox& bounds)
205  {
206  const Vec3f P(this->mPHandle->get(id));
207 
208  // build area of influence depending on point position
209  CoordBBox intersectBox(ijk, ijk);
210  if (P.x() < 0.0f) intersectBox.min().x() -= 1;
211  else intersectBox.max().x() += 1;
212  if (P.y() < 0.0f) intersectBox.min().y() -= 1;
213  else intersectBox.max().y() += 1;
214  if (P.z() < 0.0f) intersectBox.min().z() -= 1;
215  else intersectBox.max().z() += 1;
216  OPENVDB_ASSERT(intersectBox.volume() == 8);
217 
218  intersectBox.intersect(bounds);
219  if (intersectBox.empty()) return;
220 
221  auto* const data = this->buffer();
222  const auto& mask = *(this->mask());
223 
224  const SourceValueT s(this->mSHandle->get(id));
225  math::Vec3<RealT> centerw;
226 
227  const Coord& a(intersectBox.min());
228  const Coord& b(intersectBox.max());
229  for (Coord c = a; c.x() <= b.x(); ++c.x()) {
230  const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM); // unsigned bit shift mult
231  const RealT x = static_cast<RealT>(c.x()-ijk.x()); // distance from ijk to c
232  centerw[0] = value(P.x() - x); // center dist
233 
234  for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
235  const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
236  const RealT y = static_cast<RealT>(c.y()-ijk.y());
237  centerw[1] = value(P.y() - y);
238 
239  for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
240  OPENVDB_ASSERT(bounds.isInside(c));
241  const Index offset = ij + /*k*/(c.z() & (DIM-1u));
242  if (!mask.isOn(offset)) continue;
243  const RealT z = static_cast<RealT>(c.z()-ijk.z());
244  centerw[2] = value(P.z() - z);
245 
246  OPENVDB_ASSERT(centerw[0] >= 0.0f && centerw[0] <= 1.0f);
247  OPENVDB_ASSERT(centerw[1] >= 0.0f && centerw[1] <= 1.0f);
248  OPENVDB_ASSERT(centerw[2] >= 0.0f && centerw[2] <= 1.0f);
249 
250  const RealT weight = centerw.product();
251  data[offset] += s * weight;
252  this->mWeights[offset] += weight;
253  }
254  }
255  }
256  }
257 
258  inline bool finalize(const Coord&, const size_t)
259  {
260  auto* const data = this->buffer();
261  const auto& mask = *(this->mask());
262 
263  for (auto iter = mask.beginOn(); iter; ++iter) {
264  const Index offset = iter.pos();
265  const auto& w = this->mWeights[offset];
266  auto& v = data[offset];
267  if (!math::isZero(w)) v /= w;
268  }
269  return true;
270  }
271 };
272 
273 // @note If building with MSVC we have to use auto to deduce the return type
274 // due to a compiler bug. We can also use that for the public API - but
275 // we explicitly define it in non-msvc builds to ensure the API remains
276 // consistent
277 template <bool Staggered,
278  typename ValueT,
279  typename CodecT,
280  typename PositionCodecT,
281  typename FilterT,
282  typename PointDataTreeT>
283 inline
284 #ifndef _MSC_VER
285 typename TrilinearTraits<ValueT, Staggered>::template TreeT<PointDataTreeT>::Ptr
286 #else
287 auto
288 #endif
289 rasterizeTrilinear(const PointDataTreeT& points,
290  const size_t pidx,
291  const size_t sidx,
292  const FilterT& filter)
293 {
294  using TraitsT = TrilinearTraits<ValueT, Staggered>;
295  using TargetTreeT = typename TraitsT::template TreeT<PointDataTree>;
296  using TransferT = typename std::conditional<Staggered,
297  StaggeredTransfer<TargetTreeT, PositionCodecT, ValueT, CodecT>,
298  CellCenteredTransfer<TargetTreeT, PositionCodecT, ValueT, CodecT>
299  >::type;
300 
301  typename TargetTreeT::Ptr tree(new TargetTreeT);
302  if (std::is_same<FilterT, NullFilter>::value) {
303  tree->topologyUnion(points);
304  }
305  else {
306  using MaskTreeT = typename PointDataTreeT::template ValueConverter<ValueMask>::Type;
307  auto mask = convertPointsToMask<PointDataTreeT, MaskTreeT>(points, filter);
308  tree->topologyUnion(*mask);
309  }
310 
311  TransferT transfer(pidx, sidx, *tree);
312  tools::dilateActiveValues(*tree, transfer.range(),
314 
315  rasterize<PointDataTreeT, TransferT>(points, transfer, filter);
316  return tree;
317 }
318 
319 } // namespace rasterize_trilinear_internal
320 
321 /// @endcond
322 
323 ///////////////////////////////////////////////////
324 
325 template <bool Staggered,
326  typename ValueT,
327  typename FilterT,
328  typename PointDataTreeT>
329 inline auto
330 rasterizeTrilinear(const PointDataTreeT& points,
331  const std::string& attribute,
332  const FilterT& filter)
333 {
334  using TraitsT = TrilinearTraits<ValueT, Staggered>;
335  using TargetTreeT = typename TraitsT::template TreeT<PointDataTree>;
336 
337  const auto iter = points.cbeginLeaf();
338  if (!iter) return typename TargetTreeT::Ptr(new TargetTreeT);
339 
340  const AttributeSet::Descriptor& descriptor = iter->attributeSet().descriptor();
341  const size_t pidx = descriptor.find("P");
342  const size_t sidx = descriptor.find(attribute);
343  if (pidx == AttributeSet::INVALID_POS) {
344  OPENVDB_THROW(RuntimeError, "Failed to find position attribute");
345  }
346  if (sidx == AttributeSet::INVALID_POS) {
347  OPENVDB_THROW(RuntimeError, "Failed to find source attribute");
348  }
349 
350  const NamePair& ptype = descriptor.type(pidx);
351  const NamePair& stype = descriptor.type(sidx);
352  if (ptype.second == NullCodec::name()) {
353  if (stype.second == NullCodec::name()) {
355  <Staggered, ValueT, NullCodec, NullCodec>
356  (points, pidx, sidx, filter);
357  }
358  else {
360  <Staggered, ValueT, UnknownCodec, NullCodec>
361  (points, pidx, sidx, filter);
362  }
363  }
364  else {
365  if (stype.second == NullCodec::name()) {
367  <Staggered, ValueT, NullCodec, UnknownCodec>
368  (points, pidx, sidx, filter);
369  }
370  else {
372  <Staggered, ValueT, UnknownCodec, UnknownCodec>
373  (points, pidx, sidx, filter);
374  }
375  }
376 }
377 
378 
379 } // namespace points
380 } // namespace OPENVDB_VERSION_NAME
381 } // namespace openvdb
382 
383 #endif //OPENVDB_POINTS_RASTERIZE_TRILINEAR_IMPL_HAS_BEEN_INCLUDED
static const bool IsVec
Definition: Types.h:245
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Definition: PointRasterizeTrilinear.h:36
Definition: Morphology.h:82
int32_t Int32
Definition: Types.h:56
Index32 Index
Definition: Types.h:54
OPENVDB_IMPORT void initialize()
Global registration of native Grid, Transform, Metadata and Point attribute types. Also initializes blosc (if enabled).
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:337
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1057
Definition: Exceptions.h:13
math::Vec3< float > Vec3f
Definition: Types.h:74
std::pair< Name, Name > NamePair
Definition: AttributeArray.h:40
Definition: AttributeArray.h:433
Definition: AttributeArray.h:436
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
auto rasterizeTrilinear(const PointDataTreeT &points, const std::string &attribute, const FilterT &filter=NullFilter())
Perform weighted trilinear rasterization of all points within a voxel. This method takes and returns ...
Definition: PointRasterizeTrilinearImpl.h:330
typename T::ValueType ElementType
Definition: Types.h:302
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
Definition: Exceptions.h:63