OpenVDB 11.0.0
Loading...
Searching...
No Matches
PointRasterizeTrilinearImpl.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 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
12namespace openvdb {
14namespace OPENVDB_VERSION_NAME {
15namespace points {
16
17/// @cond OPENVDB_DOCS_INTERNAL
18
19namespace rasterize_trilinear_internal {
20
21template <typename TreeType,
22 typename PositionCodecT,
23 typename SourceValueT,
24 typename SourceCodecT>
25struct 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
85protected:
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
93template <typename TreeType,
94 typename PositionCodecT,
95 typename SourceValueT,
96 typename SourceCodecT>
97struct 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
104 static_assert(VecTraits<typename TreeType::ValueType>::IsVec,
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 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
184template <typename TreeType,
185 typename PositionCodecT,
186 typename SourceValueT,
187 typename SourceCodecT>
188struct 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 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 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 assert(centerw[0] >= 0.0f && centerw[0] <= 1.0f);
247 assert(centerw[1] >= 0.0f && centerw[1] <= 1.0f);
248 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
277template <bool Staggered,
278 typename ValueT,
279 typename CodecT,
280 typename PositionCodecT,
281 typename FilterT,
282 typename PointDataTreeT>
283inline
284#ifndef _MSC_VER
285typename TrilinearTraits<ValueT, Staggered>::template TreeT<PointDataTreeT>::Ptr
286#else
287auto
288#endif
289rasterizeTrilinear(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(),
313 tools::NN_FACE_EDGE_VERTEX, tools::EXPAND_TILES);
314
315 rasterize<PointDataTreeT, TransferT>(points, transfer, filter);
316 return tree;
317}
318
319} // namespace rasterize_trilinear_internal
320
321/// @endcond
322
323///////////////////////////////////////////////////
324
325template <bool Staggered,
326 typename ValueT,
327 typename FilterT,
328 typename PointDataTreeT>
329inline auto
330rasterizeTrilinear(const PointDataTreeT& points,
331 const std::string& attribute,
332 const FilterT& filter)
333{
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()) {
354 return rasterize_trilinear_internal::rasterizeTrilinear
355 <Staggered, ValueT, NullCodec, NullCodec>
356 (points, pidx, sidx, filter);
357 }
358 else {
359 return rasterize_trilinear_internal::rasterizeTrilinear
360 <Staggered, ValueT, UnknownCodec, NullCodec>
361 (points, pidx, sidx, filter);
362 }
363 }
364 else {
365 if (stype.second == NullCodec::name()) {
366 return rasterize_trilinear_internal::rasterizeTrilinear
367 <Staggered, ValueT, NullCodec, UnknownCodec>
368 (points, pidx, sidx, filter);
369 }
370 else {
371 return rasterize_trilinear_internal::rasterizeTrilinear
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
Definition Exceptions.h:63
BBox< Coord > CoordBBox
Definition NanoVDB.h:2535
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
Index32 Index
Definition Types.h:54
std::pair< Name, Name > NamePair
Definition AttributeArray.h:39
math::Vec3< float > Vec3f
Definition Types.h:74
OPENVDB_IMPORT void initialize()
Global registration of native Grid, Transform, Metadata and Point attribute types....
int32_t Int32
Definition Types.h:56
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Definition AttributeArray.h:470
Definition PointRasterizeTrilinear.h:36
Definition AttributeArray.h:466
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212