1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /*
5 * Copyright (c) Side Effects Software Inc.
6 *
7 * Produced by:
8 * Side Effects Software Inc
9 * 477 Richmond Street West
10 * Toronto, Ontario
11 * Canada M5V 3E7
12 * 416-504-9876
13 *
14 * NAME: GEO_PrimVDB.C ( GEO Library, C++)
15 *
16 * COMMENTS: Base class for vdbs.
17 */
18
19 #if defined(SESI_OPENVDB) || defined(SESI_OPENVDB_PRIM)
20
21 #include "GEO_PrimVDB.h"
22
23 #include <SYS/SYS_AtomicPtr.h>
24 #include <SYS/SYS_Math.h>
25 #include <SYS/SYS_MathCbrt.h>
26 #include <SYS/SYS_SharedMemory.h>
27
28 #include <UT/UT_Debug.h>
29 #include <UT/UT_Defines.h>
30 #include <UT/UT_EnvControl.h>
31 #include <UT/UT_FSATable.h>
32 #include <UT/UT_IStream.h>
33 #include <UT/UT_JSONParser.h>
34 #include <UT/UT_JSONValue.h>
35 #include <UT/UT_JSONWriter.h>
36 #include <UT/UT_Matrix.h>
37 #include <UT/UT_MatrixSolver.h>
38 #include <UT/UT_MemoryCounter.h>
39 #include <UT/UT_SharedMemoryManager.h>
40 #include <UT/UT_SharedPtr.h>
41 #include <UT/UT_SparseArray.h>
42 #include <UT/UT_StackTrace.h>
43 #include <UT/UT_SysClone.h>
44 #include <UT/UT_UniquePtr.h>
45 #include "UT_VDBUtils.h"
46 #include <UT/UT_Vector.h>
47 #include <UT/UT_XformOrder.h>
48
49 #include <CE/CE_Context.h>
50 #include <CE/CE_VDBGrid.h>
51
52 #include <GA/GA_AttributeRefMap.h>
53 #include <GA/GA_AttributeRefMapDestHandle.h>
54 #include <GA/GA_Defragment.h>
55 #include <GA/GA_ElementWrangler.h>
56 #include <GA/GA_IntrinsicMacros.h>
57 #include <GA/GA_LoadMap.h>
58 #include <GA/GA_MergeMap.h>
59 #include <GA/GA_PrimitiveJSON.h>
60 #include <GA/GA_RangeMemberQuery.h>
61 #include <GA/GA_SaveMap.h>
62 #include <GA/GA_WeightedSum.h>
63 #include <GA/GA_WorkVertexBuffer.h>
64
65 #include <GEO/GEO_AttributeHandleList.h>
66 #include <GEO/GEO_Detail.h>
67 #include <GEO/GEO_PrimType.h>
68 #include <GEO/GEO_PrimVolume.h>
69 #include <GEO/GEO_VolumeOptions.h>
70 #include <GEO/GEO_WorkVertexBuffer.h>
71
72 #include <openvdb/io/Stream.h>
73 #include <openvdb/math/Maps.h>
74 #include <openvdb/tools/Composite.h>
75 #include <openvdb/tools/Interpolation.h>
76 #include <openvdb/tools/LevelSetMeasure.h>
77 #include <openvdb/tools/Statistics.h>
78 #include <openvdb/tools/VectorTransformer.h>
79
80 #include <iostream>
81 #include <stdexcept>
82
83 using namespace UT::Literal;
84
85 static const UT_StringHolder theKWVertex = "vertex"_sh;
86 static const UT_StringHolder theKWVDB = "vdb"_sh;
87 static const UT_StringHolder theKWVDBShm = "vdbshm"_sh;
88 static const UT_StringHolder theKWVDBVis = "vdbvis"_sh;
89
90
91 GEO_PrimVDB::UniqueId
nextUniqueId()92 GEO_PrimVDB::nextUniqueId()
93 {
94 static AtomicUniqueId theUniqueId;
95 return static_cast<UniqueId>(theUniqueId.add(1));
96 }
97
98
GEO_PrimVDB(GEO_Detail * d,GA_Offset offset)99 GEO_PrimVDB::GEO_PrimVDB(GEO_Detail *d, GA_Offset offset)
100 : GEO_Primitive(d, offset)
101 , myGridAccessor()
102 , myVis()
103 , myUniqueId(GEO_PrimVDB::nextUniqueId())
104 , myTreeUniqueId(GEO_PrimVDB::nextUniqueId())
105 , myMetadataUniqueId(GEO_PrimVDB::nextUniqueId())
106 , myTransformUniqueId(GEO_PrimVDB::nextUniqueId())
107 , myCEGrid(0)
108 , myCEGridAuthorative(false)
109 , myCEGridIsOwned(true)
110 {
111 }
112
~GEO_PrimVDB()113 GEO_PrimVDB::~GEO_PrimVDB()
114 {
115 if (myCEGridIsOwned)
116 delete myCEGrid;
117 }
118
119 void
stashed(bool beingstashed,GA_Offset offset)120 GEO_PrimVDB::stashed(bool beingstashed, GA_Offset offset)
121 {
122 // NB: Base class must be unstashed before we can call allocateVertex().
123 GEO_Primitive::stashed(beingstashed, offset);
124 if (!beingstashed)
125 {
126 // Reset to state as if freshly constructed
127 myVis = GEO_VolumeOptions();
128 myUniqueId.relaxedStore(GEO_PrimVDB::nextUniqueId());
129 myTreeUniqueId.relaxedStore(GEO_PrimVDB::nextUniqueId());
130 myMetadataUniqueId.relaxedStore(GEO_PrimVDB::nextUniqueId());
131 myTransformUniqueId.relaxedStore(GEO_PrimVDB::nextUniqueId());
132 }
133 else
134 {
135 // Free the grid and transform.
136 // This also makes sure that myGridAccessor will be
137 // as if freshly constructed when unstashing.
138 myGridAccessor.clear();
139
140 // Throw away any GPU cache.
141 if (myCEGridIsOwned)
142 delete myCEGrid;
143 myCEGrid = nullptr;
144 myCEGridAuthorative = false;
145 myCEGridIsOwned = true;
146 }
147 // Set our internal state to default
148 myVis = GEO_VolumeOptions(GEO_VOLUMEVIS_SMOKE, /*iso*/0.0, /*density*/1.0,
149 GEO_VOLUMEVISLOD_FULL);
150 }
151
152 bool
evaluatePointRefMap(GA_Offset result_vtx,GA_AttributeRefMap & map,fpreal,fpreal,uint,uint) const153 GEO_PrimVDB::evaluatePointRefMap(GA_Offset result_vtx,
154 GA_AttributeRefMap &map,
155 fpreal /*u*/, fpreal /*v*/,
156 uint /*du*/, uint /*dv*/) const
157 {
158 map.copyValue(GA_ATTRIB_VERTEX, result_vtx,
159 GA_ATTRIB_VERTEX, getVertexOffset(0));
160 return true;
161 }
162
163 // Houdini assumes that the depth scaling of the frustum is being done in the
164 // linear part of the NonlinearFrustumMap. This method ensures that if the
165 // grid has a frustum depth not equal to 1, then it returns an equivalent map
166 // which does.
167 static openvdb::math::NonlinearFrustumMap::ConstPtr
geoStandardFrustumMapPtr(const GEO_PrimVDB & vdb)168 geoStandardFrustumMapPtr(const GEO_PrimVDB &vdb)
169 {
170 using namespace openvdb::math;
171 using openvdb::Vec3d;
172
173 const Transform &transform = vdb.getGrid().transform();
174 UT_ASSERT(transform.baseMap()->isType<NonlinearFrustumMap>());
175 NonlinearFrustumMap::ConstPtr
176 frustum_map = transform.constMap<NonlinearFrustumMap>();
177
178 // If the depth is already 1, then just return the original
179 OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
180 if (frustum_map->getDepth() == 1.0)
181 return frustum_map;
182 OPENVDB_NO_FP_EQUALITY_WARNING_END
183
184 AffineMap secondMap = frustum_map->secondMap();
185 secondMap.accumPreScale(Vec3d(1, 1, frustum_map->getDepth()));
186 return NonlinearFrustumMap::ConstPtr(
187 new NonlinearFrustumMap(frustum_map->getBBox(),
188 frustum_map->getTaper(),
189 /*depth*/1.0, secondMap.copy()));
190 }
191
192 // The returned space's fromVoxelSpace() method will convert 0-1
193 // coordinates over the bbox extents to world space (and vice versa for
194 // toVoxelSpace()).
195 GEO_PrimVolumeXform
getSpaceTransform(const UT_BoundingBoxD & bbox) const196 GEO_PrimVDB::getSpaceTransform(const UT_BoundingBoxD &bbox) const
197 {
198 using namespace openvdb;
199 using namespace openvdb::math;
200 using openvdb::Vec3d;
201 using openvdb::Mat4d;
202
203 MapBase::ConstPtr base_map = getGrid().transform().baseMap();
204 BBoxd active_bbox(UTvdbConvert(bbox.minvec()),
205 UTvdbConvert(bbox.maxvec()));
206 UT_Matrix4D transform(1.0); // identity
207 fpreal new_taper(1.0); // no taper default
208
209 // If the active_bbox is empty(), we do not want to produce a singular
210 // matrix.
211 if (active_bbox.empty())
212 {
213 if (base_map->isType<NonlinearFrustumMap>())
214 {
215 // Ideally, we use getFrustumBounds() here but it returns the
216 // wrong type.
217 const NonlinearFrustumMap& frustum_map =
218 *getGrid().transform().constMap<NonlinearFrustumMap>();
219 active_bbox = frustum_map.getBBox();
220 active_bbox.translate(Vec3d(+0.5));
221 }
222 else
223 {
224 // Use a single voxel from index origin to act as a pass-through
225 active_bbox = BBoxd(Vec3d(0.0), 1.0);
226 }
227 }
228
229 // Shift the active_bbox by half a voxel to account for the fact that
230 // UT_VoxelArray's index coordinates are cell-centred while grid indices
231 // are cell edge-aligned.
232 active_bbox.translate(Vec3d(-0.5));
233
234 if (base_map->isType<NonlinearFrustumMap>())
235 {
236 // NOTES:
237 // - VDB's nonlinear frustum goes from index-space to world-space in
238 // two steps:
239 // 1. From index-space to NDC space where we have [-0.5,+0.5] XY
240 // on the Z=0 plane, tapered outwards to to the Z=1 plane
241 // (where depth=1).
242 // 2. Then the base_map->secondMap() is applied to convert it
243 // into world-space.
244 // - Our goal is to come up with an equivalent transform which goes
245 // from [-1,+1] unit-space to world-space, matching GEO_PrimVDB's
246 // node-centred samples with GEO_PrimVolume's cell-centered samples.
247
248 NonlinearFrustumMap::ConstPtr map_ptr = geoStandardFrustumMapPtr(*this);
249 const NonlinearFrustumMap& frustum_map = *map_ptr;
250
251 // Math below only handles NonlinearFrustumMap's with a depth of 1.
252 UT_ASSERT(frustum_map.getDepth() == 1);
253
254 BBoxd frustum_bbox = frustum_map.getBBox();
255 UT_Vector3D frustum_size(UTvdbConvert(frustum_bbox.extents()));
256 UT_Vector3D inv_frustum_size = 1.0 / frustum_size;
257
258 // Find active_bbox in the 1 by 1 -> taper by taper frustum
259 UT_Vector3D active_size(UTvdbConvert(active_bbox.extents()));
260 UT_Vector3D offset_uniform =
261 UTvdbConvert(active_bbox.min() - frustum_bbox.min())
262 * inv_frustum_size;
263 UT_Vector3 scale = active_size * inv_frustum_size;
264
265 // Compute the z coordinates of 'active_bbox' in the
266 // 0-1 space (normed to the frustum size)
267 fpreal z_min = offset_uniform.z();
268 fpreal z_max = offset_uniform.z() + scale.z();
269
270 // We need a new taper since active_bbox might have a different
271 // near/far plane ratio. The mag values are calculated using VDB's
272 // taper function but then we divide min by max because Houdini's taper
273 // ratio is inversed.
274 fpreal frustum_taper = frustum_map.getTaper();
275 fpreal gamma = 1.0/frustum_taper - 1.0;
276 fpreal mag_min = 1 + gamma * z_min;
277 fpreal mag_max = 1 + gamma * z_max;
278 new_taper = mag_min / mag_max;
279
280 // xform will go from 0-1 voxel space to the tapered, near=1x1 frustum
281 UT_Matrix4D xform(1);
282 xform.scale(mag_min, mag_min, 1.0);
283
284 xform.scale(0.5, 0.5, 0.5);
285 xform.scale(scale.x(), scale.y(), scale.z());
286
287 // Scale our correctly tapered box
288 // Offset the correctly scaled and tapered box into the right xy
289 // position.
290 OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
291 if (gamma != 0.0)
292 OPENVDB_NO_FP_EQUALITY_WARNING_END
293 {
294 // Scale by the inverse of the taper since NonlinearFrustumMap
295 // creates tapers in the -z direction (a positive taper will
296 // increase the ratio of the near / far) but only scales the far
297 // face (effectively, the +z face is scaled by 1 / taper and
298 // the -z face is kept at 1.0)
299 xform.scale(1.0 / new_taper, 1.0 / new_taper, 1.0);
300
301 // The distance from the near plane that the tapered frustum sides
302 // meet ie. the position of the z-plane that gets mapped to 0 in
303 // the taper map
304 fpreal z_i = 1.0 / gamma;
305 xform.translate(0, 0, z_i + z_min + 0.5 * scale.z());
306
307 // Compute the shear: it is the offset on the near plane of the
308 // current frustum to where we want it to be
309 UT_Vector3D frustum_center(0.5*frustum_size);
310 UT_Vector3D active_center(0.5*active_size);
311 // The current active_bbox position
312 UT_Vector3D bbox_offset = frustum_center - active_center;
313 // Compute the offset to the real position. We add back half-voxel
314 // so that the reference base is at the old min
315 UT_Vector3D shear = UTvdbConvert(active_bbox.min() + Vec3d(0.5))
316 - bbox_offset;
317 shear *= inv_frustum_size;
318 shear /= z_i;
319 xform.shear(0, shear.x(), shear.y());
320
321 // Translate ourselves back so that the tapered plane of
322 // frustum_bbox is at 0
323 xform.translate(0, 0, -z_i);
324 }
325 else
326 {
327 // Translate bottom corner of box to (0,0,0)
328 xform.translate(-0.5, -0.5, 0.0);
329 xform.translate(0.5*scale.x(), 0.5*scale.y(), 0.5*scale.z());
330 xform.translate(offset_uniform.x(),
331 offset_uniform.y(),
332 offset_uniform.z());
333 }
334
335 // `transform` now brings us from a tapered (1*y/x) box to world space,
336 // We want to go from a tapered, near=1x1 frustum to world space, so
337 // prescale by the aspect
338 fpreal aspect = frustum_bbox.extents().y() / frustum_bbox.extents().x();
339 xform.scale(1.0, aspect, 1.0);
340
341 Mat4d mat4 = frustum_map.secondMap().getMat4();
342 transform = xform * UTvdbConvert(mat4);
343 }
344 else
345 {
346 // NOTES:
347 // - VDB's grid transform goes from index-space to world-space.
348 // - Our goal is to come up with an equivalent transform which goes
349 // from [-1,+1] unit-space to world-space, matching GEO_PrimVDB's
350 // node-centred samples with GEO_PrimVolume's cell-centered samples.
351 // (NOTE: fromVoxelSpace() converts from [0,+1] to [-1,+1])
352
353 // Create transform which converts [-1,+1] unit-space to [0,+1]
354 transform.translate(1.0, 1.0, 1.0);
355 transform.scale(0.5, 0.5, 0.5);
356
357 // Go from the [0,1] volume space, to [min,max] where
358 // min and max are in the vdb index space. Note that UT_VoxelArray
359 // doesn't support negative coordinates which is why we need to shift
360 // the index origin to the bbox min.
361 transform.scale(active_bbox.extents().x(),
362 active_bbox.extents().y(),
363 active_bbox.extents().z());
364 transform.translate(active_bbox.min().x(),
365 active_bbox.min().y(),
366 active_bbox.min().z());
367
368 // Post-multiply by the affine map which converts index to world space
369 transform = transform * UTvdbConvert(base_map->getAffineMap()->getMat4());
370 }
371
372 UT_Vector3 translate;
373 transform.getTranslates(translate);
374
375 GEO_PrimVolumeXform result;
376 result.myXform = transform;
377 result.myCenter = translate;
378
379 OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
380 result.myHasTaper = (new_taper != 1.0);
381 OPENVDB_NO_FP_EQUALITY_WARNING_END
382
383 transform.invert();
384 result.myInverseXform = transform;
385
386 result.myTaperX = new_taper;
387 result.myTaperY = new_taper;
388
389 return result;
390 }
391
392 // Return a GEO_PrimVolumeXform which maps [-0.5,+0.5] Houdini voxel space
393 // coordinates over the VDB's active voxel bbox into world space.
394 GEO_PrimVolumeXform
getSpaceTransform() const395 GEO_PrimVDB::getSpaceTransform() const
396 {
397 const openvdb::CoordBBox bbox = getGrid().evalActiveVoxelBoundingBox();
398 return getSpaceTransform(UTvdbConvert(bbox));
399 }
400
401 bool
conditionMatrix(UT_Matrix4D & mat4)402 GEO_PrimVDB::conditionMatrix(UT_Matrix4D &mat4)
403 {
404 // This tolerance is just one factor larger than what
405 // AffineMap::updateAcceleration() uses to ensure that we never trigger the
406 // exception.
407 const double tol = 4.0 * openvdb::math::Tolerance<double>::value();
408 const double min_diag = SYScbrt(tol);
409 if (!SYSequalZero(mat4.determinant3(), tol))
410 {
411 // openvdb::math::simplify uses openvdb::math::isApproxEqual to detect
412 // uniform scaling, which has a more stringent tolerance than SYSisEqual.
413 // As a result we often have uniform voxel / axis aligned Volumes being
414 // converted to VDBs with Maps that are simplified to ScaleTranslate
415 // rather than UniformScaleTranslate. The latter is more correct, plus
416 // some operations (e.g. LevelSetMorph) don't work with non-Uniform
417 // scales. So we unify the diagonal if the diagonals are SYSisEqual,
418 // but not exactly equal.
419 if (SYSisEqual(mat4(0, 0), mat4(1, 1)) &&
420 SYSisEqual(mat4(0, 0), mat4(2, 2)) &&
421 !(mat4(0, 0) == mat4(1,1) && mat4(0, 0) == mat4(2,2)))
422 {
423 // Unify to mat(0, 0) like openvdb::math::simplify.
424 mat4(1, 1) = mat4(2, 2) = mat4(0, 0);
425 return true;
426 }
427 if (SYSalmostEqual((float)mat4(0, 0), (float)mat4(1, 1)) &&
428 SYSalmostEqual((float)mat4(0, 0), (float)mat4(2, 2)) &&
429 !(mat4(0, 0) == mat4(1,1) && mat4(0, 0) == mat4(2,2)))
430 {
431 // Unify to mat(0, 0) like openvdb::math::simplify.
432 mat4(1, 1) = mat4(2, 2) = mat4(0, 0);
433 return true;
434 }
435 return false;
436 }
437 UT_MatrixSolverD solver;
438 UT_Matrix3D mat3(mat4);
439 const int N = 3;
440 UT_MatrixD m(1,N, 1,N), v(1,N, 1,N), diag(1,N, 1,N), tmp(1,N, 1,N);
441 UT_VectorD w(1,N);
442
443 m.setSubmatrix3(1, 1, mat3);
444 if (!solver.SVDDecomp(m, w, v))
445 {
446 // Give up and return a scale matrix as small as possible
447 mat4.identity();
448 mat4.scale(min_diag, min_diag, min_diag);
449 }
450 else
451 {
452 v.transpose(tmp);
453 v = tmp;
454 diag.makeIdentity();
455 for (int i = 1; i <= N; i++)
456 diag(i,i) = SYSmax(min_diag, w(i));
457 m.postMult(diag, tmp);
458 tmp.postMult(v, m);
459 m.getSubmatrix3(mat3, 1, 1);
460 mat4 = mat3;
461 }
462 return true;
463 }
464
465 // All AffineMap creation must to through this to avoid crashes when passing
466 // singular matrices into OpenVDB
467 template<typename T>
468 static openvdb::SharedPtr<T>
geoCreateAffineMap(const UT_Matrix4D & mat4)469 geoCreateAffineMap(const UT_Matrix4D& mat4)
470 {
471 using namespace openvdb::math;
472
473 openvdb::SharedPtr<T> transform;
474 UT_Matrix4D new_mat4(mat4);
475 (void) GEO_PrimVDB::conditionMatrix(new_mat4);
476 try
477 {
478 transform.reset(new AffineMap(UTvdbConvert(new_mat4)));
479 }
480 catch (openvdb::ArithmeticError &)
481 {
482 // Fall back to trying to clear the last column first, since
483 // VDB seems to not like that, instead of falling back to identity.
484 new_mat4 = mat4;
485 new_mat4[0][3] = 0;
486 new_mat4[1][3] = 0;
487 new_mat4[2][3] = 0;
488 new_mat4[3][3] = 1;
489 (void) GEO_PrimVDB::conditionMatrix(new_mat4);
490 try
491 {
492 transform.reset(new AffineMap(UTvdbConvert(new_mat4)));
493 }
494 catch (openvdb::ArithmeticError &)
495 {
496 UT_ASSERT(!"Failed to create affine map");
497 transform.reset(new AffineMap());
498 }
499 }
500 return transform;
501 }
502
503 // All calls to createLinearTransform with a matrix4 must to through this to
504 // avoid crashes when passing singular matrices into OpenVDB
505 static openvdb::math::Transform::Ptr
geoCreateLinearTransform(const UT_Matrix4D & mat4)506 geoCreateLinearTransform(const UT_Matrix4D& mat4)
507 {
508 using namespace openvdb::math;
509 return Transform::Ptr(new Transform(geoCreateAffineMap<MapBase>(mat4)));
510 }
511
512 void
setSpaceTransform(const GEO_PrimVolumeXform & space,const UT_Vector3R & resolution,bool force_taper)513 GEO_PrimVDB::setSpaceTransform(
514 const GEO_PrimVolumeXform &space,
515 const UT_Vector3R &resolution,
516 bool force_taper)
517 {
518 using namespace openvdb;
519 using namespace openvdb::math;
520 using openvdb::Vec3d;
521
522 // VDB's nonlinear frustum goes from index-space to world-space in
523 // two steps:
524 // 1. From index-space to NDC space where we have [-0.5,+0.5] XY
525 // on the Z=0 plane, tapered outwards to to the Z=1 plane
526 // (where depth=1).
527 // 2. Then the base_map->secondMap() is applied to convert it
528 // into world-space.
529 // On the other hand, 'space' converts from [-1,+1] space over the given
530 // resolution into world-space.
531 //
532 // Our goal is to come up with an equivalent NonlinearTransformMap of
533 // 'space' which converts from index space to world-space, matching
534 // GEO_PrimVDB's node-centred samples with GEO_PrimVolume's cell-centered
535 // samples.
536
537 Transform::Ptr xform;
538
539 if (force_taper || space.myHasTaper)
540 {
541 // VDB only supports a single taper value so use average of the two
542 fpreal taper = 0.5*(space.myTaperX + space.myTaperY);
543
544 // Create a matrix which goes from vdb's post-taper XY(-0.5,+0.5) space
545 // to [-1,1] space so that we can post-multiply by space's transform to
546 // get into world-space.
547 //
548 // NonlinearFrustumMap use's 1/taper as its taper value, going from
549 // Z=0 to Z=1. So we first scale it by the taper to undo this.
550 UT_Matrix4D transform(1.0);
551 transform.scale(taper, taper, 1.0);
552 // Account for aspect ratio
553 transform.scale(1.0, resolution.x() / resolution.y(), 1.0);
554 // We now go from XY(-0.5,+0.5)/Z(0,1) to XY(-1,+1)/Z(2)
555 transform.scale(2.0, 2.0, 2.0);
556 // Translate into [-1,+1] on all axes by centering the Z axis
557 transform.translate(0.0, 0.0, -1.0);
558
559 // Now apply the space's linear transform
560 UT_Matrix4D linear;
561 linear = space.myXform; // convert UT_Matrix3 to UT_Matrix4
562 transform *= linear;
563 transform.translate(
564 space.myCenter.x(), space.myCenter.y(), space.myCenter.z());
565
566 // In order to get VDB to match Houdini, we create a frustum map using
567 // Houdini's bbox, so we offset by -0.5 voxel in order taper the
568 // Houdini bbox in VDB index space. This allows us to align Houdini's
569 // cell-centered samples with VDB node-centered ones.
570 BBoxd bbox(Vec3d(0.0), UTvdbConvert(resolution));
571 bbox.translate(Vec3d(-0.5));
572
573 // Build a NonlinearFrustumMap from this
574 MapBase::Ptr affine_map(geoCreateAffineMap<MapBase>(transform));
575 xform.reset(new Transform(MapBase::Ptr(
576 new NonlinearFrustumMap(bbox, taper, /*depth*/1.0, affine_map))));
577 }
578 else // optimize for a linear transform if we have no taper
579 {
580 // NOTES:
581 // - Houdini's transform goes from [-1,+1] unit-space to world-space.
582 // - VDB's grid transform goes from index-space to world-space.
583
584 UT_Matrix4D matx(/*identity*/1.0);
585 UT_Matrix4D mult;
586
587 // UT_VoxelArray's index coordinates are cell-centred while grid
588 // indices are cell edge-aligned. So first offset the VDB indices by
589 // +0.5 voxel to convert them into Houdini indices.
590 matx.translate(0.5, 0.5, 0.5);
591
592 // Now convert the indices from [0,dim] to [-1,+1]
593 matx.scale(2.0/resolution(0), 2.0/resolution(1), 2.0/resolution(2));
594 matx.translate(-1.0, -1.0, -1.0);
595
596 // Post-multiply with Houdini transform to get result that converts
597 // into world-space.
598 mult = space.myXform;
599 matx *= mult;
600 matx.translate(space.myCenter(0), space.myCenter(1), space.myCenter(2));
601
602 // Create a linear transform using the matrix
603 xform = geoCreateLinearTransform(matx);
604 }
605
606 myGridAccessor.setTransform(*xform, *this);
607 }
608
609 // This will give you the a GEO_PrimVolumeXform. Given an index, it will
610 // compute the world space position of that index.
611 GEO_PrimVolumeXform
getIndexSpaceTransform() const612 GEO_PrimVDB::getIndexSpaceTransform() const
613 {
614 using namespace openvdb;
615 using namespace openvdb::math;
616 using openvdb::Vec3d;
617 using openvdb::Mat4d;
618
619 // This taper function follows from the conversion code in
620 // GEO_PrimVolume::fromVoxelSpace() until until myXform/myCenter is
621 // applied. It has been simplified somewhat, and uses the definition that
622 // gamma = taper - 1.
623 struct Local
624 {
625 static fpreal taper(fpreal x, fpreal z, fpreal gamma)
626 {
627 return 2 * (x - 0.5) * (1 + gamma * (1 - z));
628 }
629 };
630
631 fpreal new_taper = 1.0;
632 UT_Matrix4D transform(1.0); // identity
633
634 MapBase::ConstPtr base_map = getGrid().transform().baseMap();
635 if (base_map->isType<NonlinearFrustumMap>())
636 {
637 // NOTES:
638 // - VDB's nonlinear frustum goes from index-space to world-space in
639 // two steps:
640 // 1. From index-space to NDC space where we have [-0.5,+0.5] XY
641 // on the Z=0 plane, tapered outwards to to the Z=1 plane
642 // (where depth=1).
643 // 2. Then the base_map->secondMap() is applied to convert it
644 // into world-space.
645 // - Our goal is to come up with an equivalent GEO_PrimVolumeXform
646 // which goes from index-space to world-space, matching GEO_PrimVDB's
647 // node-centred samples with GEO_PrimVolume's cell-centered samples.
648 // - The gotcha here is that callers use fromVoxelSpace() which will
649 // first do a mapping of [0,1] to [-1,+1] which we need to undo.
650
651 NonlinearFrustumMap::ConstPtr map_ptr = geoStandardFrustumMapPtr(*this);
652 const NonlinearFrustumMap& frustum_map = *map_ptr;
653
654 // Math below only handles NonlinearFrustumMap's with a depth of 1.
655 UT_ASSERT(frustum_map.getDepth() == 1);
656 fpreal taper = frustum_map.getTaper();
657
658 // We need to create a taper value for fromVoxelSpace()'s incoming
659 // Houdini index space coordinate, so the bbox we want to taper with is
660 // actually the Houdini index bbox.
661 UT_BoundingBox bbox;
662 getFrustumBounds(bbox);
663
664 fpreal x = bbox.xsize();
665 fpreal y = bbox.ysize();
666 fpreal z = bbox.zsize();
667
668 Vec3d real_min(bbox.xmin(), bbox.ymin(), bbox.zmin());
669 Vec3d real_max(bbox.xmax(), bbox.ymax(), bbox.zmax());
670
671 // Compute a new taper based on the expected ratio of these two
672 // z positions
673 fpreal z_min = real_min.z();
674 fpreal z_max = real_max.z();
675
676 //
677 // If t = (1+g(1-a))/(1+g(1-b)) then g = (1-t)/(t(1-b) - (1-a))
678 // where t = taper; g = new_taper - 1, a = z_min, b = z_max;
679 //
680 fpreal new_gamma = (1 - taper) / (taper * (1 - z_max) - (1 - z_min));
681 new_taper = new_gamma + 1;
682
683 // Since we are tapering the index space, the taper map adds a
684 // scale and a shear so we find these and invert them
685
686 fpreal x_max_pos = Local::taper(real_max.x(), z_max, new_gamma);
687 fpreal x_min_pos = Local::taper(real_min.x(), z_max, new_gamma);
688 // Now, move x_max_pos = -x_min_pos with a shear
689 fpreal x_scale = x_max_pos - x_min_pos;
690 fpreal shear_x = 0.5 * x_scale - x_max_pos;
691
692 // Do the same for y
693 fpreal y_max_pos = Local::taper(real_max.y(), z_max, new_gamma);
694 fpreal y_min_pos = Local::taper(real_min.y(), z_max, new_gamma);
695 fpreal y_scale = y_max_pos - y_min_pos;
696 fpreal shear_y = 0.5 * y_scale - y_max_pos;
697
698 transform.translate(0, 0, -2*(z_min - 0.5));
699
700 // Scale z so that our frustum depth range is 0-1
701 transform.scale(1, 1, 0.5/z);
702
703 // Apply the shear
704 OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
705 if (taper != 1.0)
706 OPENVDB_NO_FP_EQUALITY_WARNING_END
707 {
708 fpreal z_i = 1.0 / (taper - 1);
709 transform.translate(0, 0, -z_i-1.0);
710 transform.shear(0, -shear_x / z_i, -shear_y / z_i);
711 transform.translate(0, 0, z_i+1.0);
712 }
713 else
714 {
715 transform.translate(shear_x, shear_y, 0.0);
716 }
717
718 transform.scale(1.0/x_scale, 1.0/y_scale, 1.0);
719
720 // Scale by 1/taper to convert taper definitions
721 transform.scale(1.0 / taper, 1.0 / taper, 1.0);
722
723 // Account for aspect ratio
724 transform.scale(1, y/x, 1);
725
726 Mat4d mat4 = frustum_map.secondMap().getMat4();
727 transform *= UTvdbConvert(mat4);
728 }
729 else
730 {
731 // We only deal with nonlinear maps that are frustum maps
732 UT_ASSERT(base_map->isLinear()
733 && "Found unexpected nonlinear MapBase.");
734
735 // Since VDB's transform is already from index-space to world-space, we
736 // just need to undo the [0,1] -> [-1,+1] mapping that fromVoxelSpace()
737 // does before transforming by myXform/myCenter. The math is thus:
738 // scale(1/2)*translate(0.5)
739 // But we also want to shift VDB's node-centred samples to match
740 // GEO_PrimVolume's cell-centered ones so we want:
741 // scale(1/2)*translate(0.5)*translate(-0.5)
742 // This reduces down to just scale(1/2)
743 //
744 transform.scale(0.5, 0.5, 0.5);
745
746 transform *= UTvdbConvert(base_map->getAffineMap()->getMat4());
747 }
748
749 GEO_PrimVolumeXform result;
750 result.myXform = transform;
751 transform.getTranslates(result.myCenter);
752
753 OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
754 result.myHasTaper = (new_taper != 1.0);
755 OPENVDB_NO_FP_EQUALITY_WARNING_END
756
757 transform.invert();
758 result.myInverseXform = transform;
759
760 result.myTaperX = new_taper;
761 result.myTaperY = new_taper;
762
763 return result;
764 }
765
766 bool
isSDF() const767 GEO_PrimVDB::isSDF() const
768 {
769 if (getGrid().getGridClass() == openvdb::GRID_LEVEL_SET)
770 return true;
771
772 return false;
773 }
774
775 fpreal
getTaper() const776 GEO_PrimVDB::getTaper() const
777 {
778 return getSpaceTransform().myTaperX;
779 }
780
781 void
reverse()782 GEO_PrimVDB::reverse()
783 {
784 }
785
786 UT_Vector3
computeNormal() const787 GEO_PrimVDB::computeNormal() const
788 {
789 return UT_Vector3(0, 0, 0);
790 }
791
792
793
794 template <typename GridType>
795 static void
geo_calcVolume(const GridType & grid,fpreal & volume)796 geo_calcVolume(const GridType &grid, fpreal &volume)
797 {
798 bool calculated = false;
799 if (grid.getGridClass() == openvdb::GRID_LEVEL_SET) {
800 try {
801 volume = openvdb::tools::levelSetVolume(grid);
802 calculated = true;
803 } catch (std::exception& /*e*/) {
804 // do nothing
805 }
806 }
807
808 // Simply account for the total number of active voxels
809 if (!calculated) {
810 const openvdb::Vec3d size = grid.voxelSize();
811 volume = size[0] * size[1] * size[2] * grid.activeVoxelCount();
812 }
813 }
814
815 fpreal
calcVolume(const UT_Vector3 &) const816 GEO_PrimVDB::calcVolume(const UT_Vector3 &) const
817 {
818 fpreal volume = 0;
819 UTvdbCallAllTopology(getStorageType(), geo_calcVolume, getGrid(), volume);
820 return volume;
821 }
822
823 template <typename GridType>
824 static void
geo_calcArea(const GridType & grid,fpreal & area)825 geo_calcArea(const GridType &grid, fpreal &area)
826 {
827 bool calculated = false;
828 if (grid.getGridClass() == openvdb::GRID_LEVEL_SET) {
829 try {
830 area = openvdb::tools::levelSetArea(grid);
831 calculated = true;
832 } catch (std::exception& /*e*/) {
833 // do nothing
834 }
835 }
836
837 if (!calculated) {
838 using LeafIter = typename GridType::TreeType::LeafCIter;
839 using VoxelIter = typename GridType::TreeType::LeafNodeType::ValueOnCIter;
840 using openvdb::Coord;
841 const Coord normals[] = {Coord(0,0,-1), Coord(0,0,1), Coord(-1,0,0),
842 Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0)};
843 // NOTE: we assume rectangular prism voxels
844 openvdb::Vec3d voxel_size = grid.voxelSize();
845 const fpreal areas[] = {fpreal(voxel_size.x() * voxel_size.y()),
846 fpreal(voxel_size.x() * voxel_size.y()),
847 fpreal(voxel_size.y() * voxel_size.z()),
848 fpreal(voxel_size.y() * voxel_size.z()),
849 fpreal(voxel_size.z() * voxel_size.x()),
850 fpreal(voxel_size.z() * voxel_size.x())};
851 area = 0;
852 for (LeafIter leaf = grid.tree().cbeginLeaf(); leaf; ++leaf) {
853 // Visit all the active voxels in this leaf node.
854 for (VoxelIter iter = leaf->cbeginValueOn(); iter; ++iter) {
855 // Iterate through all the neighboring voxels
856 for (int i=0; i<6; i++)
857 if (!grid.tree().isValueOn(iter.getCoord() + normals[i])) area += areas[i];
858 }
859 }
860 }
861 }
862
863 fpreal
calcArea() const864 GEO_PrimVDB::calcArea() const
865 {
866 // Calculate the surface area of all the exterior voxels.
867 fpreal area = 0;
868 UTvdbCallAllTopology(getStorageType(), geo_calcArea, getGrid(), area);
869 return area;
870 }
871
872 void
enlargePointBounds(UT_BoundingBox & box) const873 GEO_PrimVDB::enlargePointBounds(UT_BoundingBox &box) const
874 {
875 UT_BoundingBox qbox;
876 if (getBBox(&qbox))
877 box.enlargeBounds(qbox);
878 }
879
880 bool
enlargeBoundingBox(UT_BoundingRect & box,const GA_Attribute * P) const881 GEO_PrimVDB::enlargeBoundingBox(UT_BoundingRect &box,
882 const GA_Attribute *P) const
883 {
884 const GA_Detail &gdp = getDetail();
885
886 if (!P)
887 P = gdp.getP();
888 else if (P != gdp.getP())
889 return GEO_Primitive::enlargeBoundingBox(box, P);
890
891 UT_BoundingBox my_bbox;
892 if (getBBox(&my_bbox))
893 {
894 box.enlargeBounds(my_bbox.xmin(), my_bbox.ymin());
895 box.enlargeBounds(my_bbox.xmax(), my_bbox.ymax());
896 return true;
897 }
898 return true;
899 }
900
901 bool
enlargeBoundingBox(UT_BoundingBox & box,const GA_Attribute * P) const902 GEO_PrimVDB::enlargeBoundingBox(UT_BoundingBox &box,
903 const GA_Attribute *P) const
904 {
905 const GA_Detail &gdp = getDetail();
906
907 if (!P)
908 P = gdp.getP();
909 else if (P != gdp.getP())
910 return GEO_Primitive::enlargeBoundingBox(box, P);
911
912 UT_BoundingBox my_bbox;
913 if (getBBox(&my_bbox))
914 {
915 box.enlargeBounds(my_bbox);
916 return true;
917 }
918 return true;
919 }
920
921 bool
enlargeBoundingSphere(UT_BoundingSphere & sphere,const GA_Attribute * P) const922 GEO_PrimVDB::enlargeBoundingSphere(UT_BoundingSphere &sphere,
923 const GA_Attribute *P) const
924 {
925 const GA_Detail &gdp = getDetail();
926
927 if (!P)
928 P = gdp.getP();
929 else if (P != gdp.getP())
930 return GEO_Primitive::enlargeBoundingSphere(sphere, P);
931
932 addToBSphere(&sphere);
933 return true;
934 }
935
936 int64
getBaseMemoryUsage() const937 GEO_PrimVDB::getBaseMemoryUsage() const
938 {
939 exint mem = 0;
940 if (hasGrid())
941 mem += getGrid().memUsage();
942 return mem;
943 }
944
945 void
countBaseMemory(UT_MemoryCounter & counter) const946 GEO_PrimVDB::countBaseMemory(UT_MemoryCounter &counter) const
947 {
948 if (hasGrid())
949 {
950 // NOTE: We don't share the grid object or its transform
951 // We don't know what type of Grid we have, but apart from what's
952 // in GridBase, it just has a shared pointer to the tree extra,
953 // so just add that in separately.
954 counter.countUnshared(sizeof(openvdb::GridBase) + sizeof(openvdb::TreeBase::Ptr));
955 // We don't know what type of MapBase the Transform uses,
956 // so just guess the largest one
957 counter.countUnshared(sizeof(openvdb::math::Transform) + sizeof(openvdb::math::NonlinearFrustumMap));
958
959 // The grid's tree is shared. In order to get the reference count,
960 // we need to get our own shared pointer to it, then use one less
961 // than the ref count (ours counts as one).
962 exint refcount;
963 exint size;
964 const openvdb::TreeBase *ptr;
965 {
966 openvdb::TreeBase::ConstPtr ref = getGrid().constBaseTreePtr();
967 refcount = ref.use_count() - 1;
968 size = ref->memUsage();
969 ptr = ref.get();
970 }
971 counter.countShared(size, refcount, ptr);
972 }
973 }
974
975
976 template <typename GridType>
977 static inline typename GridType::ValueType
geo_doubleToGridValue(double val)978 geo_doubleToGridValue(double val)
979 {
980 using ValueT = typename GridType::ValueType;
981 // This ugly construction avoids compiler warnings when,
982 // for example, initializing an openvdb::Vec3i with a double.
983 return ValueT(openvdb::zeroVal<ValueT>() + val);
984 }
985
986
987 template <typename GridType>
988 static fpreal
geo_sampleGrid(const GridType & grid,const UT_Vector3 & pos)989 geo_sampleGrid(const GridType &grid, const UT_Vector3 &pos)
990 {
991 const openvdb::math::Transform & xform = grid.transform();
992 openvdb::math::Vec3d vpos;
993 typename GridType::ValueType value;
994
995 vpos = openvdb::math::Vec3d(pos.x(), pos.y(), pos.z());
996 vpos = xform.worldToIndex(vpos);
997
998 openvdb::tools::BoxSampler::sample(grid.tree(), vpos, value);
999
1000 fpreal result = value;
1001
1002 return result;
1003 }
1004
1005 template <typename GridType>
1006 static fpreal
geo_sampleBoolGrid(const GridType & grid,const UT_Vector3 & pos)1007 geo_sampleBoolGrid(const GridType &grid, const UT_Vector3 &pos)
1008 {
1009 const openvdb::math::Transform & xform = grid.transform();
1010 openvdb::math::Vec3d vpos;
1011 typename GridType::ValueType value;
1012
1013 vpos = openvdb::math::Vec3d(pos.x(), pos.y(), pos.z());
1014 vpos = xform.worldToIndex(vpos);
1015
1016 openvdb::tools::PointSampler::sample(grid.tree(), vpos, value);
1017
1018 fpreal result = value;
1019
1020 return result;
1021 }
1022
1023 template <typename GridType>
1024 static UT_Vector3D
geo_sampleGridV3(const GridType & grid,const UT_Vector3 & pos)1025 geo_sampleGridV3(const GridType &grid, const UT_Vector3 &pos)
1026 {
1027 const openvdb::math::Transform & xform = grid.transform();
1028 openvdb::math::Vec3d vpos;
1029 typename GridType::ValueType value;
1030
1031 vpos = openvdb::math::Vec3d(pos.x(), pos.y(), pos.z());
1032 vpos = xform.worldToIndex(vpos);
1033
1034 openvdb::tools::BoxSampler::sample(grid.tree(), vpos, value);
1035
1036 UT_Vector3D result;
1037 result.x() = double(value[0]);
1038 result.y() = double(value[1]);
1039 result.z() = double(value[2]);
1040
1041 return result;
1042 }
1043
1044 template <typename GridType, typename T, typename IDX>
1045 static void
geo_sampleGridMany(const GridType & grid,T * f,int stride,const IDX * pos,int num)1046 geo_sampleGridMany(const GridType &grid,
1047 T *f, int stride,
1048 const IDX *pos,
1049 int num)
1050 {
1051 typename GridType::ConstAccessor accessor = grid.getAccessor();
1052
1053 const openvdb::math::Transform & xform = grid.transform();
1054 openvdb::math::Vec3d vpos;
1055 typename GridType::ValueType value;
1056
1057
1058 for (int i = 0; i < num; i++)
1059 {
1060 vpos = openvdb::math::Vec3d(pos[i].x(), pos[i].y(), pos[i].z());
1061 vpos = xform.worldToIndex(vpos);
1062
1063 openvdb::tools::BoxSampler::sample(accessor, vpos, value);
1064
1065 *f = T(value);
1066 f += stride;
1067 }
1068 }
1069
1070 template <typename GridType, typename T, typename IDX>
1071 static void
geo_sampleBoolGridMany(const GridType & grid,T * f,int stride,const IDX * pos,int num)1072 geo_sampleBoolGridMany(const GridType &grid,
1073 T *f, int stride,
1074 const IDX *pos,
1075 int num)
1076 {
1077 typename GridType::ConstAccessor accessor = grid.getAccessor();
1078
1079 const openvdb::math::Transform & xform = grid.transform();
1080 openvdb::math::Vec3d vpos;
1081 typename GridType::ValueType value;
1082
1083
1084 for (int i = 0; i < num; i++)
1085 {
1086 vpos = openvdb::math::Vec3d(pos[i].x(), pos[i].y(), pos[i].z());
1087 vpos = xform.worldToIndex(vpos);
1088
1089 openvdb::tools::PointSampler::sample(accessor, vpos, value);
1090
1091 *f = T(value);
1092 f += stride;
1093 }
1094 }
1095
1096 template <typename GridType, typename T, typename IDX>
1097 static void
geo_sampleVecGridMany(const GridType & grid,T * f,int stride,const IDX * pos,int num)1098 geo_sampleVecGridMany(const GridType &grid,
1099 T *f, int stride,
1100 const IDX *pos,
1101 int num)
1102 {
1103 typename GridType::ConstAccessor accessor = grid.getAccessor();
1104
1105 const openvdb::math::Transform & xform = grid.transform();
1106 openvdb::math::Vec3d vpos;
1107 typename GridType::ValueType value;
1108
1109
1110 for (int i = 0; i < num; i++)
1111 {
1112 vpos = openvdb::math::Vec3d(pos[i].x(), pos[i].y(), pos[i].z());
1113 vpos = xform.worldToIndex(vpos);
1114
1115 openvdb::tools::BoxSampler::sample(accessor, vpos, value);
1116
1117 f->x() = value[0];
1118 f->y() = value[1];
1119 f->z() = value[2];
1120 f += stride;
1121 }
1122 }
1123
1124 static fpreal
geoEvaluateVDB(const GEO_PrimVDB * vdb,const UT_Vector3 & pos)1125 geoEvaluateVDB(const GEO_PrimVDB *vdb, const UT_Vector3 &pos)
1126 {
1127 UTvdbReturnScalarType(vdb->getStorageType(), geo_sampleGrid, vdb->getGrid(), pos);
1128 UTvdbReturnBoolType(vdb->getStorageType(), geo_sampleBoolGrid, vdb->getGrid(), pos);
1129 return 0;
1130 }
1131
1132 static UT_Vector3D
geoEvaluateVDB_V3(const GEO_PrimVDB * vdb,const UT_Vector3 & pos)1133 geoEvaluateVDB_V3(const GEO_PrimVDB *vdb, const UT_Vector3 &pos)
1134 {
1135 UTvdbReturnVec3Type(vdb->getStorageType(),
1136 geo_sampleGridV3, vdb->getGrid(), pos);
1137 return UT_Vector3D(0, 0, 0);
1138 }
1139
1140 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,float * f,int stride,const UT_Vector3 * pos,int num)1141 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, float *f, int stride, const UT_Vector3 *pos, int num)
1142 {
1143 UTvdbReturnScalarType(vdb->getStorageType(),
1144 geo_sampleGridMany, vdb->getGrid(), f, stride, pos, num);
1145 UTvdbReturnBoolType(vdb->getStorageType(),
1146 geo_sampleBoolGridMany, vdb->getGrid(), f, stride, pos, num);
1147 for (int i = 0; i < num; i++)
1148 {
1149 *f = 0;
1150 f += stride;
1151 }
1152 }
1153
1154 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,int * f,int stride,const UT_Vector3 * pos,int num)1155 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, int *f, int stride, const UT_Vector3 *pos, int num)
1156 {
1157 UTvdbReturnScalarType(vdb->getStorageType(),
1158 geo_sampleGridMany, vdb->getGrid(), f, stride, pos, num);
1159 UTvdbReturnBoolType(vdb->getStorageType(),
1160 geo_sampleBoolGridMany, vdb->getGrid(), f, stride, pos, num);
1161 for (int i = 0; i < num; i++)
1162 {
1163 *f = 0;
1164 f += stride;
1165 }
1166 }
1167
1168 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,UT_Vector3 * f,int stride,const UT_Vector3 * pos,int num)1169 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, UT_Vector3 *f, int stride, const UT_Vector3 *pos, int num)
1170 {
1171 UTvdbReturnVec3Type(vdb->getStorageType(),
1172 geo_sampleVecGridMany, vdb->getGrid(), f, stride, pos, num);
1173 for (int i = 0; i < num; i++)
1174 {
1175 *f = 0;
1176 f += stride;
1177 }
1178 }
1179
1180 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,double * f,int stride,const UT_Vector3D * pos,int num)1181 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, double *f, int stride, const UT_Vector3D *pos, int num)
1182 {
1183 UTvdbReturnScalarType(vdb->getStorageType(),
1184 geo_sampleGridMany, vdb->getGrid(), f, stride, pos, num);
1185 UTvdbReturnBoolType(vdb->getStorageType(),
1186 geo_sampleBoolGridMany, vdb->getGrid(), f, stride, pos, num);
1187 for (int i = 0; i < num; i++)
1188 {
1189 *f = 0;
1190 f += stride;
1191 }
1192 }
1193
1194 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,exint * f,int stride,const UT_Vector3D * pos,int num)1195 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, exint *f, int stride, const UT_Vector3D *pos, int num)
1196 {
1197 UTvdbReturnScalarType(vdb->getStorageType(),
1198 geo_sampleGridMany, vdb->getGrid(), f, stride, pos, num);
1199 UTvdbReturnBoolType(vdb->getStorageType(),
1200 geo_sampleBoolGridMany, vdb->getGrid(), f, stride, pos, num);
1201 for (int i = 0; i < num; i++)
1202 {
1203 *f = 0;
1204 f += stride;
1205 }
1206 }
1207
1208 static void
geoEvaluateVDBMany(const GEO_PrimVDB * vdb,UT_Vector3D * f,int stride,const UT_Vector3D * pos,int num)1209 geoEvaluateVDBMany(const GEO_PrimVDB *vdb, UT_Vector3D *f, int stride, const UT_Vector3D *pos, int num)
1210 {
1211 UTvdbReturnVec3Type(vdb->getStorageType(),
1212 geo_sampleVecGridMany, vdb->getGrid(), f, stride, pos, num);
1213 for (int i = 0; i < num; i++)
1214 {
1215 *f = 0;
1216 f += stride;
1217 }
1218 }
1219
1220 fpreal
getValueF(const UT_Vector3 & pos) const1221 GEO_PrimVDB::getValueF(const UT_Vector3 &pos) const
1222 {
1223 return geoEvaluateVDB(this, pos);
1224 }
1225
1226 UT_Vector3D
getValueV3(const UT_Vector3 & pos) const1227 GEO_PrimVDB::getValueV3(const UT_Vector3 &pos) const
1228 {
1229 return geoEvaluateVDB_V3(this, pos);
1230 }
1231
1232 void
getValues(float * f,int stride,const UT_Vector3 * pos,int num) const1233 GEO_PrimVDB::getValues(float *f, int stride, const UT_Vector3 *pos, int num) const
1234 {
1235 return geoEvaluateVDBMany(this, f, stride, pos, num);
1236 }
1237
1238 void
getValues(int * f,int stride,const UT_Vector3 * pos,int num) const1239 GEO_PrimVDB::getValues(int *f, int stride, const UT_Vector3 *pos, int num) const
1240 {
1241 return geoEvaluateVDBMany(this, f, stride, pos, num);
1242 }
1243
1244 void
getValues(UT_Vector3 * f,int stride,const UT_Vector3 * pos,int num) const1245 GEO_PrimVDB::getValues(UT_Vector3 *f, int stride, const UT_Vector3 *pos, int num) const
1246 {
1247 return geoEvaluateVDBMany(this, f, stride, pos, num);
1248 }
1249
1250 void
getValues(double * f,int stride,const UT_Vector3D * pos,int num) const1251 GEO_PrimVDB::getValues(double *f, int stride, const UT_Vector3D *pos, int num) const
1252 {
1253 return geoEvaluateVDBMany(this, f, stride, pos, num);
1254 }
1255
1256 void
getValues(exint * f,int stride,const UT_Vector3D * pos,int num) const1257 GEO_PrimVDB::getValues(exint *f, int stride, const UT_Vector3D *pos, int num) const
1258 {
1259 return geoEvaluateVDBMany(this, f, stride, pos, num);
1260 }
1261
1262 void
getValues(UT_Vector3D * f,int stride,const UT_Vector3D * pos,int num) const1263 GEO_PrimVDB::getValues(UT_Vector3D *f, int stride, const UT_Vector3D *pos, int num) const
1264 {
1265 return geoEvaluateVDBMany(this, f, stride, pos, num);
1266 }
1267
1268 namespace // anonymous
1269 {
1270
1271 template <bool NORMALIZE>
1272 class geo_EvalGradients
1273 {
1274 public:
geo_EvalGradients(UT_Vector3 * gradients,int stride,const UT_Vector3 * positions,int num_positions)1275 geo_EvalGradients(
1276 UT_Vector3 *gradients,
1277 int stride,
1278 const UT_Vector3 *positions,
1279 int num_positions)
1280 : myGradients(gradients)
1281 , myStride(stride)
1282 , myPos(positions)
1283 , myNumPos(num_positions)
1284 {
1285 }
1286
1287 template<typename GridT>
operator ()(const GridT & grid)1288 void operator()(const GridT &grid)
1289 {
1290 using namespace openvdb;
1291 using AccessorT = typename GridT::ConstAccessor;
1292 using ValueT = typename GridT::ValueType;
1293
1294 const math::Transform & xform = grid.transform();
1295 const math::Vec3d dim = grid.voxelSize();
1296 const double vox_size = SYSmin(dim[0], dim[1], dim[2]);
1297 const double h = 0.5 * vox_size;
1298 const math::Vec3d mask[] =
1299 { math::Vec3d(h, 0, 0)
1300 , math::Vec3d(0, h, 0)
1301 , math::Vec3d(0, 0, h)
1302 };
1303 AccessorT accessor = grid.getConstAccessor();
1304 UT_Vector3 * gradient = myGradients;
1305
1306 for (int i = 0; i < myNumPos; i++, gradient += myStride)
1307 {
1308 const math::Vec3d pos(myPos[i].x(), myPos[i].y(), myPos[i].z());
1309
1310 for (int j = 0; j < 3; j++)
1311 {
1312 const math::Vec3d vpos0 = xform.worldToIndex(pos - mask[j]);
1313 const math::Vec3d vpos1 = xform.worldToIndex(pos + mask[j]);
1314 ValueT v0, v1;
1315
1316 tools::BoxSampler::sample<AccessorT>(accessor, vpos0, v0);
1317 tools::BoxSampler::sample<AccessorT>(accessor, vpos1, v1);
1318 if (NORMALIZE)
1319 (*gradient)(j) = (v1 - v0);
1320 else
1321 (*gradient)(j) = (v1 - v0) / vox_size;
1322 }
1323 if (NORMALIZE)
1324 gradient->normalize();
1325 }
1326 }
1327
1328 private:
1329 UT_Vector3 * myGradients;
1330 int myStride;
1331 const UT_Vector3 * myPos;
1332 int myNumPos;
1333 };
1334
1335 } // namespace anonymous
1336
1337 bool
evalGradients(UT_Vector3 * gradients,int stride,const UT_Vector3 * pos,int num_pos,bool normalize) const1338 GEO_PrimVDB::evalGradients(
1339 UT_Vector3 *gradients,
1340 int stride,
1341 const UT_Vector3 *pos,
1342 int num_pos,
1343 bool normalize) const
1344 {
1345 if (normalize)
1346 {
1347 geo_EvalGradients<true> eval(gradients, stride, pos, num_pos);
1348 return UTvdbProcessTypedGridScalar(getStorageType(), getGrid(), eval);
1349 }
1350 else
1351 {
1352 geo_EvalGradients<false> eval(gradients, stride, pos, num_pos);
1353 return UTvdbProcessTypedGridScalar(getStorageType(), getGrid(), eval);
1354 }
1355 }
1356
1357 bool
isAligned(const GEO_PrimVDB * vdb) const1358 GEO_PrimVDB::isAligned(const GEO_PrimVDB *vdb) const
1359 {
1360 if (getGrid().transform() == vdb->getGrid().transform())
1361 return true;
1362 return false;
1363 }
1364
1365 bool
isWorldAxisAligned() const1366 GEO_PrimVDB::isWorldAxisAligned() const
1367 {
1368 // Tapered are trivially not aligned.
1369 if (!SYSisEqual(getTaper(), 1))
1370 return false;
1371 UT_Matrix4D x;
1372
1373 x = getTransform4();
1374
1375 for (int i = 0; i < 3; i++)
1376 {
1377 for (int j = 0; j < 3; j++)
1378 {
1379 if (i == j)
1380 {
1381 if (x(i, j) <= 0)
1382 return false;
1383 }
1384 else
1385 {
1386 if (x(i,j) != 0)
1387 return false;
1388 }
1389 }
1390 }
1391
1392 return true;
1393 }
1394
1395 bool
isActiveRegionMatched(const GEO_PrimVDB * vdb) const1396 GEO_PrimVDB::isActiveRegionMatched(const GEO_PrimVDB *vdb) const
1397 {
1398 if (!isAligned(vdb))
1399 return false;
1400 // Ideally we'd invoke hasSameTopology?
1401 return vdb->getGrid().baseTreePtr() == getGrid().baseTreePtr();
1402 }
1403
1404 CE_VDBGrid *
getCEGrid(bool read,bool write) const1405 GEO_PrimVDB::getCEGrid(bool read, bool write) const
1406 {
1407 UT_ASSERT(!write);
1408 UT_ASSERT(read);
1409
1410 if (myCEGrid)
1411 return myCEGrid;
1412
1413 CE_VDBGrid *cegrid = new CE_VDBGrid();
1414
1415 try
1416 {
1417 cegrid->initFromVDB(getGrid());
1418 myCEGridIsOwned = true;
1419 }
1420 catch (cl::Error &err)
1421 {
1422 CE_Context::reportError(err);
1423 delete cegrid;
1424 cegrid = 0;
1425 }
1426
1427 myCEGrid = cegrid;
1428 myCEGridAuthorative = write;
1429
1430 return myCEGrid;
1431 }
1432
1433 void
flushCEWriteCaches()1434 GEO_PrimVDB::flushCEWriteCaches()
1435 {
1436 if (!myCEGrid)
1437 return;
1438 if (myCEGridAuthorative)
1439 {
1440 // Write back.
1441 UT_ASSERT(!"Not implemented");
1442 myCEGridAuthorative = false;
1443 }
1444 }
1445
1446 void
flushCECaches()1447 GEO_PrimVDB::flushCECaches()
1448 {
1449 // Write back if needed
1450 flushCEWriteCaches();
1451 // Destroy the grid.
1452 if (myCEGridIsOwned)
1453 delete myCEGrid;
1454 myCEGrid = 0;
1455 myCEGridAuthorative = false;
1456 myCEGridIsOwned = true;
1457 }
1458
1459 void
stealCEBuffers(const GA_Primitive * psrc)1460 GEO_PrimVDB::stealCEBuffers(const GA_Primitive *psrc)
1461 {
1462 UT_ASSERT(psrc && getTypeId() == psrc->getTypeId());
1463 if (!psrc || getTypeId() != psrc->getTypeId())
1464 return;
1465 const GEO_PrimVDB *src = UTverify_cast<const GEO_PrimVDB *>(psrc);
1466
1467 myCEGrid = src->myCEGrid;
1468 myCEGridAuthorative = src->myCEGridAuthorative;
1469 myCEGridIsOwned = src->myCEGridIsOwned;
1470
1471 src->myCEGrid = 0;
1472 src->myCEGridAuthorative = false;
1473 src->myCEGridIsOwned = true;
1474 }
1475
1476
1477 void
indexToPos(int x,int y,int z,UT_Vector3 & pos) const1478 GEO_PrimVDB::indexToPos(int x, int y, int z, UT_Vector3 &pos) const
1479 {
1480 openvdb::math::Vec3d vpos;
1481
1482 vpos = openvdb::math::Vec3d(x, y, z);
1483 vpos = getGrid().indexToWorld(vpos);
1484 pos = UT_Vector3(vpos[0], vpos[1], vpos[2]);
1485 }
1486
1487 void
findexToPos(UT_Vector3 idx,UT_Vector3 & pos) const1488 GEO_PrimVDB::findexToPos(UT_Vector3 idx, UT_Vector3 &pos) const
1489 {
1490 openvdb::math::Vec3d vpos;
1491
1492 vpos = openvdb::math::Vec3d(idx.x(), idx.y(), idx.z());
1493 vpos = getGrid().indexToWorld(vpos);
1494 pos = UT_Vector3(vpos[0], vpos[1], vpos[2]);
1495 }
1496
1497 void
posToIndex(UT_Vector3 pos,int & x,int & y,int & z) const1498 GEO_PrimVDB::posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const
1499 {
1500 openvdb::math::Vec3d vpos(pos.data());
1501 openvdb::math::Coord
1502 coord = getGrid().transform().worldToIndexCellCentered(vpos);
1503 x = coord.x();
1504 y = coord.y();
1505 z = coord.z();
1506 }
1507
1508 void
posToIndex(UT_Vector3 pos,UT_Vector3 & index) const1509 GEO_PrimVDB::posToIndex(UT_Vector3 pos, UT_Vector3 &index) const
1510 {
1511 openvdb::math::Vec3d vpos;
1512
1513 vpos = openvdb::math::Vec3d(pos.x(), pos.y(), pos.z());
1514 vpos = getGrid().worldToIndex(vpos);
1515
1516 index = UTvdbConvert(vpos);
1517 }
1518
1519 void
indexToPos(exint x,exint y,exint z,UT_Vector3D & pos) const1520 GEO_PrimVDB::indexToPos(exint x, exint y, exint z, UT_Vector3D &pos) const
1521 {
1522 openvdb::math::Vec3d vpos;
1523
1524 vpos = openvdb::math::Vec3d(x, y, z);
1525 vpos = getGrid().indexToWorld(vpos);
1526 pos = UT_Vector3D(vpos[0], vpos[1], vpos[2]);
1527 }
1528
1529 void
findexToPos(UT_Vector3D idx,UT_Vector3D & pos) const1530 GEO_PrimVDB::findexToPos(UT_Vector3D idx, UT_Vector3D &pos) const
1531 {
1532 openvdb::math::Vec3d vpos;
1533
1534 vpos = openvdb::math::Vec3d(idx.x(), idx.y(), idx.z());
1535 vpos = getGrid().indexToWorld(vpos);
1536 pos = UT_Vector3D(vpos[0], vpos[1], vpos[2]);
1537 }
1538
1539 void
posToIndex(UT_Vector3D pos,exint & x,exint & y,exint & z) const1540 GEO_PrimVDB::posToIndex(UT_Vector3D pos, exint &x, exint &y, exint &z) const
1541 {
1542 openvdb::math::Vec3d vpos(pos.data());
1543 openvdb::math::Coord
1544 coord = getGrid().transform().worldToIndexCellCentered(vpos);
1545 x = coord.x();
1546 y = coord.y();
1547 z = coord.z();
1548 }
1549
1550 void
posToIndex(UT_Vector3D pos,UT_Vector3D & index) const1551 GEO_PrimVDB::posToIndex(UT_Vector3D pos, UT_Vector3D &index) const
1552 {
1553 openvdb::math::Vec3d vpos;
1554
1555 vpos = openvdb::math::Vec3d(pos.x(), pos.y(), pos.z());
1556 vpos = getGrid().worldToIndex(vpos);
1557
1558 index = UTvdbConvert(vpos);
1559 }
1560
1561 template <typename GridType>
1562 static fpreal
geo_sampleIndex(const GridType & grid,int ix,int iy,int iz)1563 geo_sampleIndex(const GridType &grid, int ix, int iy, int iz)
1564 {
1565 openvdb::math::Coord xyz;
1566 typename GridType::ValueType value;
1567
1568 xyz = openvdb::math::Coord(ix, iy, iz);
1569
1570 value = grid.tree().getValue(xyz);
1571
1572 fpreal result = value;
1573
1574 return result;
1575 }
1576
1577 template <typename GridType>
1578 static UT_Vector3D
geo_sampleIndexV3(const GridType & grid,int ix,int iy,int iz)1579 geo_sampleIndexV3(const GridType &grid, int ix, int iy, int iz)
1580 {
1581 openvdb::math::Coord xyz;
1582 typename GridType::ValueType value;
1583
1584 xyz = openvdb::math::Coord(ix, iy, iz);
1585
1586 value = grid.tree().getValue(xyz);
1587
1588 UT_Vector3D result;
1589
1590 result.x() = double(value[0]);
1591 result.y() = double(value[1]);
1592 result.z() = double(value[2]);
1593
1594 return result;
1595 }
1596
1597 template <typename GridType, typename T, typename IDX>
1598 static void
geo_sampleIndexMany(const GridType & grid,T * f,int stride,const IDX * ix,const IDX * iy,const IDX * iz,int num)1599 geo_sampleIndexMany(const GridType &grid,
1600 T *f, int stride,
1601 const IDX *ix, const IDX *iy, const IDX *iz,
1602 int num)
1603 {
1604 typename GridType::ConstAccessor accessor = grid.getAccessor();
1605
1606 openvdb::math::Coord xyz;
1607 typename GridType::ValueType value;
1608
1609 for (int i = 0; i < num; i++)
1610 {
1611 xyz = openvdb::math::Coord(ix[i], iy[i], iz[i]);
1612
1613 value = accessor.getValue(xyz);
1614
1615 *f = T(value);
1616 f += stride;
1617 }
1618 }
1619
1620 template <typename GridType, typename T, typename IDX>
1621 static void
geo_sampleVecIndexMany(const GridType & grid,T * f,int stride,const IDX * ix,const IDX * iy,const IDX * iz,int num)1622 geo_sampleVecIndexMany(const GridType &grid,
1623 T *f, int stride,
1624 const IDX *ix, const IDX *iy, const IDX *iz,
1625 int num)
1626 {
1627 typename GridType::ConstAccessor accessor = grid.getAccessor();
1628
1629 openvdb::math::Coord xyz;
1630 typename GridType::ValueType value;
1631
1632 for (int i = 0; i < num; i++)
1633 {
1634 xyz = openvdb::math::Coord(ix[i], iy[i], iz[i]);
1635
1636 value = accessor.getValue(xyz);
1637
1638 f->x() = value[0];
1639 f->y() = value[1];
1640 f->z() = value[2];
1641 f += stride;
1642 }
1643 }
1644
1645 static fpreal
geoEvaluateIndexVDB(const GEO_PrimVDB * vdb,int ix,int iy,int iz)1646 geoEvaluateIndexVDB(const GEO_PrimVDB *vdb,
1647 int ix, int iy, int iz)
1648 {
1649 UTvdbReturnScalarType(vdb->getStorageType(),
1650 geo_sampleIndex, vdb->getGrid(), ix, iy, iz);
1651
1652 return 0.0;
1653 }
1654
1655 static UT_Vector3D
geoEvaluateIndexVDB_V3(const GEO_PrimVDB * vdb,int ix,int iy,int iz)1656 geoEvaluateIndexVDB_V3(const GEO_PrimVDB *vdb,
1657 int ix, int iy, int iz)
1658 {
1659 UTvdbReturnVec3Type(vdb->getStorageType(),
1660 geo_sampleIndexV3, vdb->getGrid(), ix, iy, iz);
1661
1662 return UT_Vector3D(0.0, 0, 0);
1663 }
1664
1665 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,float * f,int stride,const int * ix,const int * iy,const int * iz,int num)1666 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1667 float *f, int stride,
1668 const int *ix, const int *iy, const int *iz,
1669 int num)
1670 {
1671 UTvdbReturnScalarType(vdb->getStorageType(),
1672 geo_sampleIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1673 for (int i = 0; i < num; i++)
1674 {
1675 *f = 0;
1676 f += stride;
1677 }
1678 }
1679
1680 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,int * f,int stride,const int * ix,const int * iy,const int * iz,int num)1681 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1682 int *f, int stride,
1683 const int *ix, const int *iy, const int *iz,
1684 int num)
1685 {
1686 UTvdbReturnScalarType(vdb->getStorageType(),
1687 geo_sampleIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1688 for (int i = 0; i < num; i++)
1689 {
1690 *f = 0;
1691 f += stride;
1692 }
1693 }
1694
1695 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,UT_Vector3 * f,int stride,const int * ix,const int * iy,const int * iz,int num)1696 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1697 UT_Vector3 *f, int stride,
1698 const int *ix, const int *iy, const int *iz,
1699 int num)
1700 {
1701 UTvdbReturnVec3Type(vdb->getStorageType(),
1702 geo_sampleVecIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1703 for (int i = 0; i < num; i++)
1704 {
1705 *f = 0;
1706 f += stride;
1707 }
1708 }
1709
1710 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,double * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num)1711 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1712 double *f, int stride,
1713 const exint *ix, const exint *iy, const exint *iz,
1714 int num)
1715 {
1716 UTvdbReturnScalarType(vdb->getStorageType(),
1717 geo_sampleIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1718 for (int i = 0; i < num; i++)
1719 {
1720 *f = 0;
1721 f += stride;
1722 }
1723 }
1724
1725 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,exint * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num)1726 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1727 exint *f, int stride,
1728 const exint *ix, const exint *iy, const exint *iz,
1729 int num)
1730 {
1731 UTvdbReturnScalarType(vdb->getStorageType(),
1732 geo_sampleIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1733 for (int i = 0; i < num; i++)
1734 {
1735 *f = 0;
1736 f += stride;
1737 }
1738 }
1739
1740 static void
geoEvaluateIndexVDBMany(const GEO_PrimVDB * vdb,UT_Vector3D * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num)1741 geoEvaluateIndexVDBMany(const GEO_PrimVDB *vdb,
1742 UT_Vector3D *f, int stride,
1743 const exint *ix, const exint *iy, const exint *iz,
1744 int num)
1745 {
1746 UTvdbReturnVec3Type(vdb->getStorageType(),
1747 geo_sampleVecIndexMany, vdb->getGrid(), f, stride, ix, iy, iz, num);
1748 for (int i = 0; i < num; i++)
1749 {
1750 *f = 0;
1751 f += stride;
1752 }
1753 }
1754 fpreal
getValueAtIndexF(int ix,int iy,int iz) const1755 GEO_PrimVDB::getValueAtIndexF(int ix, int iy, int iz) const
1756 {
1757 return geoEvaluateIndexVDB(this, ix, iy, iz);
1758 }
1759
1760 UT_Vector3D
getValueAtIndexV3(int ix,int iy,int iz) const1761 GEO_PrimVDB::getValueAtIndexV3(int ix, int iy, int iz) const
1762 {
1763 return geoEvaluateIndexVDB_V3(this, ix, iy, iz);
1764 }
1765
1766 void
getValuesAtIndices(float * f,int stride,const int * ix,const int * iy,const int * iz,int num) const1767 GEO_PrimVDB::getValuesAtIndices(float *f, int stride, const int *ix, const int *iy, const int *iz, int num) const
1768 {
1769 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1770 }
1771
1772 void
getValuesAtIndices(int * f,int stride,const int * ix,const int * iy,const int * iz,int num) const1773 GEO_PrimVDB::getValuesAtIndices(int *f, int stride, const int *ix, const int *iy, const int *iz, int num) const
1774 {
1775 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1776 }
1777
1778 void
getValuesAtIndices(UT_Vector3 * f,int stride,const int * ix,const int * iy,const int * iz,int num) const1779 GEO_PrimVDB::getValuesAtIndices(UT_Vector3 *f, int stride, const int *ix, const int *iy, const int *iz, int num) const
1780 {
1781 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1782 }
1783
1784 void
getValuesAtIndices(double * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num) const1785 GEO_PrimVDB::getValuesAtIndices(double *f, int stride, const exint *ix, const exint *iy, const exint *iz, int num) const
1786 {
1787 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1788 }
1789
1790 void
getValuesAtIndices(exint * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num) const1791 GEO_PrimVDB::getValuesAtIndices(exint *f, int stride, const exint *ix, const exint *iy, const exint *iz, int num) const
1792 {
1793 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1794 }
1795
1796 void
getValuesAtIndices(UT_Vector3D * f,int stride,const exint * ix,const exint * iy,const exint * iz,int num) const1797 GEO_PrimVDB::getValuesAtIndices(UT_Vector3D *f, int stride, const exint *ix, const exint *iy, const exint *iz, int num) const
1798 {
1799 geoEvaluateIndexVDBMany(this, f, stride, ix, iy, iz, num);
1800 }
1801
1802 UT_Vector3
getGradient(const UT_Vector3 & pos) const1803 GEO_PrimVDB::getGradient(const UT_Vector3 &pos) const
1804 {
1805 UT_Vector3 grad;
1806
1807 grad = 0;
1808
1809 evalGradients(&grad, 1, &pos, 1, false);
1810
1811 return grad;
1812 }
1813
1814
1815 ////////////////////////////////////////
1816
1817
1818 namespace {
1819
1820 // Functor for use with UTvdbProcessTypedGridVec3() to apply a transform
1821 // to the voxel values of vector-valued grids
1822 struct gu_VecXformOp
1823 {
1824 openvdb::Mat4d mat;
gu_VecXformOp__anon16e4255b0211::gu_VecXformOp1825 gu_VecXformOp(const openvdb::Mat4d& _mat): mat(_mat) {}
operator ()__anon16e4255b0211::gu_VecXformOp1826 template<typename GridT> void operator()(GridT& grid) const
1827 {
1828 openvdb::tools::transformVectors(grid, mat);
1829 }
1830 };
1831
1832 } // unnamed namespace
1833
1834
1835 void
transform(const UT_Matrix4 & mat)1836 GEO_PrimVDB::transform(const UT_Matrix4 &mat)
1837 {
1838 if (!hasGrid()) return;
1839
1840 try {
1841 using openvdb::GridBase;
1842 using namespace openvdb::math;
1843
1844 // Get the transform
1845 const GridBase& const_grid = getConstGrid();
1846 MapBase::ConstPtr base_map = const_grid.transform().baseMap();
1847 Mat4d base_mat4 = base_map->getAffineMap()->getMat4();
1848
1849 // Get the 3x3 subcomponent of the matrix
1850 Vec3d translation = base_mat4.getTranslation();
1851 Mat3d vdbmatrix = base_mat4.getMat3();
1852
1853 // Multiply our mat with the mat3
1854 UT_Matrix3D transformed(mat);
1855 transformed = UTvdbConvert(vdbmatrix) * transformed;
1856
1857 // Put it into a mat4 and translate it
1858 UT_Matrix4D final;
1859 final = transformed;
1860 final.setTranslates(UTvdbConvert(translation));
1861
1862 // Make an affine matrix out of it
1863 AffineMap::Ptr map(geoCreateAffineMap<AffineMap>(final));
1864
1865 // Set the affine matrix from our base_map into this map
1866 MapBase::Ptr result = simplify(map);
1867 if (base_map->isType<NonlinearFrustumMap>()) {
1868 const NonlinearFrustumMap& frustum_map =
1869 *const_grid.transform().constMap<NonlinearFrustumMap>();
1870
1871 MapBase::Ptr new_frustum_map (new NonlinearFrustumMap(
1872 frustum_map.getBBox(),
1873 frustum_map.getTaper(),
1874 frustum_map.getDepth(),
1875 result));
1876
1877 result = new_frustum_map;
1878 }
1879
1880 // This sets the vertex position to `translation` as well
1881 myGridAccessor.setTransform(Transform(result), *this);
1882
1883 // If (and only if) the grid is vector-valued, apply the transform to
1884 // each voxel's value.
1885 if (const_grid.getVectorType() != openvdb::VEC_INVARIANT) {
1886 gu_VecXformOp op(UTvdbConvert(UT_Matrix4D(mat)));
1887 GEOvdbProcessTypedGridVec3(*this, op, /*make_unique*/true);
1888 }
1889
1890 } catch (std::exception& /*e*/) {
1891 UT_ASSERT(!"Failed to apply transform");
1892 }
1893 }
1894
1895
1896 void
copyGridFrom(const GEO_PrimVDB & src_prim,bool copyPosition)1897 GEO_PrimVDB::copyGridFrom(const GEO_PrimVDB& src_prim, bool copyPosition)
1898 {
1899 setGrid(src_prim.getGrid(), copyPosition); // makes a shallow copy
1900
1901 // Copy the source primitive's grid serial numbers.
1902 myTreeUniqueId.exchange(src_prim.getTreeUniqueId());
1903 myMetadataUniqueId.exchange(src_prim.getMetadataUniqueId());
1904 myTransformUniqueId.exchange(src_prim.getTransformUniqueId());
1905 }
1906
1907
1908 // If myGrid's tree is shared, replace the tree with a deep copy of itself.
1909 // Note: myGrid's metadata and transform are assumed to never be shared
1910 // (see setGrid()).
1911 void
makeGridUnique()1912 GEO_PrimVDB::GridAccessor::makeGridUnique()
1913 {
1914 if (myGrid) {
1915 UT_ASSERT(myGrid.unique());
1916 openvdb::TreeBase::Ptr localTreePtr = myGrid->baseTreePtr();
1917 if (localTreePtr.use_count() > 2) { // myGrid + localTreePtr = 2
1918 myGrid->setTree(myGrid->constBaseTree().copy());
1919 }
1920 }
1921 }
1922
1923 bool
isGridUnique() const1924 GEO_PrimVDB::GridAccessor::isGridUnique() const
1925 {
1926 if (myGrid) {
1927 // We require the grid to always be unique, it is the tree
1928 // that is allowed to be shared.
1929 UT_ASSERT(myGrid.unique());
1930 openvdb::TreeBase::Ptr localTreePtr = myGrid->baseTreePtr();
1931 if (localTreePtr.use_count() > 2) { // myGrid + localTreePtr = 2
1932 return false;
1933 }
1934 return true;
1935 }
1936 // Empty grids are trivially unique
1937 return true;
1938 }
1939
1940 void
setTransform4(const UT_Matrix4 & xform4)1941 GEO_PrimVDB::setTransform4(const UT_Matrix4 &xform4)
1942 {
1943 setTransform4(static_cast<UT_DMatrix4>(xform4));
1944 }
1945
1946 void
setTransform4(const UT_DMatrix4 & xform4)1947 GEO_PrimVDB::setTransform4(const UT_DMatrix4 &xform4)
1948 {
1949 using namespace openvdb::math;
1950 myGridAccessor.setTransform(*geoCreateLinearTransform(xform4), *this);
1951 }
1952
1953 void
getRes(int & rx,int & ry,int & rz) const1954 GEO_PrimVDB::getRes(int &rx, int &ry, int &rz) const
1955 {
1956 using namespace openvdb;
1957
1958 const GridBase & grid = getGrid();
1959 const math::Coord dim = grid.evalActiveVoxelDim();
1960
1961 // We could overflow with very large grids and get a negative
1962 // coordinate back.
1963 UT_ASSERT(dim[0] >= 0);
1964 UT_ASSERT(dim[1] >= 0);
1965 UT_ASSERT(dim[2] >= 0);
1966
1967 rx = SYSmax(dim[0], 0);
1968 ry = SYSmax(dim[1], 0);
1969 rz = SYSmax(dim[2], 0);
1970 }
1971
1972 fpreal
getVoxelDiameter() const1973 GEO_PrimVDB::getVoxelDiameter() const
1974 {
1975 UT_Vector3 p1, p2;
1976
1977 indexToPos(0, 0, 0, p1);
1978 indexToPos(1, 1, 1, p2);
1979
1980 p2 -= p1;
1981
1982 return p2.length();
1983 }
1984
1985 UT_Vector3
getVoxelSize() const1986 GEO_PrimVDB::getVoxelSize() const
1987 {
1988 UT_Vector3 p1, p2;
1989 UT_Vector3 vsize;
1990
1991 indexToPos(0, 0, 0, p1);
1992
1993 indexToPos(1, 0, 0, p2);
1994 p2 -= p1;
1995 vsize.x() = p2.length();
1996
1997 indexToPos(0, 1, 0, p2);
1998 p2 -= p1;
1999 vsize.y() = p2.length();
2000
2001 indexToPos(0, 0, 1, p2);
2002 p2 -= p1;
2003 vsize.z() = p2.length();
2004
2005 return vsize;
2006 }
2007
2008 template <typename GridType>
2009 static void
geo_calcMinVDB(GridType & grid,fpreal & result)2010 geo_calcMinVDB( GridType &grid,
2011 fpreal &result)
2012 {
2013 auto val = openvdb::tools::extrema(grid.cbeginValueOn());
2014 result = val.min();
2015 }
2016
2017 fpreal
calcMinimum() const2018 GEO_PrimVDB::calcMinimum() const
2019 {
2020 fpreal value = SYS_FPREAL_MAX;
2021 UTvdbCallScalarType( getStorageType(),
2022 geo_calcMinVDB,
2023 SYSconst_cast(getConstGrid()),
2024 value);
2025 return value;
2026 }
2027
2028 template <typename GridType>
2029 static void
geo_calcMaxVDB(GridType & grid,fpreal & result)2030 geo_calcMaxVDB( GridType &grid,
2031 fpreal &result)
2032 {
2033 auto val = openvdb::tools::extrema(grid.cbeginValueOn());
2034 result = val.max();
2035 }
2036
2037 fpreal
calcMaximum() const2038 GEO_PrimVDB::calcMaximum() const
2039 {
2040 fpreal value = -SYS_FPREAL_MAX;
2041 UTvdbCallScalarType( getStorageType(),
2042 geo_calcMaxVDB,
2043 SYSconst_cast(getConstGrid()),
2044 value);
2045 return value;
2046 }
2047
2048 template <typename GridType>
2049 static void
geo_calcAvgVDB(GridType & grid,fpreal & result)2050 geo_calcAvgVDB( GridType &grid,
2051 fpreal &result)
2052 {
2053 auto val = openvdb::tools::statistics(grid.cbeginValueOn());
2054 result = val.avg();
2055 }
2056
2057 fpreal
calcAverage() const2058 GEO_PrimVDB::calcAverage() const
2059 {
2060 fpreal value = 0;
2061 UTvdbCallScalarType( getStorageType(),
2062 geo_calcAvgVDB,
2063 SYSconst_cast(getConstGrid()),
2064 value);
2065 return value;
2066 }
2067
2068
2069 bool
getFrustumBounds(UT_BoundingBox & idxbox) const2070 GEO_PrimVDB::getFrustumBounds(UT_BoundingBox &idxbox) const
2071 {
2072 using namespace openvdb;
2073 using namespace openvdb::math;
2074 using openvdb::CoordBBox;
2075 using openvdb::Vec3d;
2076
2077 idxbox.makeInvalid();
2078
2079 // See if we have a non-linear map, this is the sign
2080 // we want to be bounded.
2081 MapBase::ConstPtr base_map = getGrid().transform().baseMap();
2082 if (base_map->isType<NonlinearFrustumMap>())
2083 {
2084 const NonlinearFrustumMap& frustum_map =
2085 *getGrid().transform().constMap<NonlinearFrustumMap>();
2086
2087 // The returned idxbox is intended to be used with
2088 // getIndexSpaceTransform() which will shift it by -0.5 voxel. So we
2089 // need to add +0.5 to compensate.
2090 BBoxd bbox = frustum_map.getBBox();
2091 bbox.translate(Vec3d(+0.5));
2092
2093 idxbox.initBounds( UTvdbConvert(bbox.min()) );
2094 idxbox.enlargeBounds( UTvdbConvert(bbox.max()) );
2095
2096 return true;
2097 }
2098
2099 return false;
2100 }
2101
2102 static bool
geoGetFrustumBoundsFromVDB(const GEO_PrimVDB * vdb,openvdb::CoordBBox & box)2103 geoGetFrustumBoundsFromVDB(const GEO_PrimVDB *vdb, openvdb::CoordBBox &box)
2104 {
2105 using namespace openvdb;
2106
2107 UT_BoundingBox clip;
2108 bool doclip;
2109
2110 doclip = vdb->getFrustumBounds(clip);
2111 if (doclip)
2112 {
2113 box = CoordBBox( Coord( (int)SYSrint(clip.xmin()), (int)SYSrint(clip.ymin()), (int)SYSrint(clip.zmin()) ),
2114 Coord( (int)SYSrint(clip.xmax()), (int)SYSrint(clip.ymax()), (int)SYSrint(clip.zmax()) ) );
2115 }
2116 return doclip;
2117 }
2118
2119 // The result of the intersection of active regions goes into grid_a
2120 template <typename GridTypeA, typename GridTypeB>
2121 static void
geoIntersect(GridTypeA & grid_a,const GridTypeB & grid_b)2122 geoIntersect(GridTypeA& grid_a, const GridTypeB &grid_b)
2123 {
2124 typename GridTypeA::Accessor access_a = grid_a.getAccessor();
2125 typename GridTypeB::ConstAccessor access_b = grid_b.getAccessor();
2126
2127 // For each on value in a, set it off if b is also off
2128 for (typename GridTypeA::ValueOnCIter
2129 iter = grid_a.cbeginValueOn(); iter; ++iter)
2130 {
2131 openvdb::CoordBBox bbox = iter.getBoundingBox();
2132 for (int k=bbox.min().z(); k<=bbox.max().z(); k++)
2133 {
2134 for (int j=bbox.min().y(); j<=bbox.max().y(); j++)
2135 {
2136 for (int i=bbox.min().x(); i<=bbox.max().x(); i++)
2137 {
2138 openvdb::Coord coord(i, j, k);
2139 if (!access_b.isValueOn(coord))
2140 {
2141 access_a.setValue(coord, grid_a.background());
2142 access_a.setValueOff(coord);
2143 }
2144 }
2145 }
2146 }
2147 }
2148 }
2149
2150 /// This class is used as a functor to set inactive voxels to the background
2151 /// value.
2152 template<typename GridType>
2153 class geoInactiveToBackground
2154 {
2155 public:
2156 typedef typename GridType::ValueOffIter Iterator;
2157 typedef typename GridType::ValueType ValueType;
2158
geoInactiveToBackground(const GridType & grid)2159 geoInactiveToBackground(const GridType& grid)
2160 {
2161 background = grid.background();
2162 }
2163
operator ()(const Iterator & iter) const2164 inline void operator()(const Iterator& iter) const
2165 {
2166 iter.setValue(background);
2167 }
2168
2169 private:
2170 ValueType background;
2171 };
2172
2173 template <typename GridType>
2174 static void
geoActivateBBox(GridType & grid,const openvdb::CoordBBox & bbox,bool setvalue,double value,GEO_PrimVDB::ActivateOperation operation,bool doclip,const openvdb::CoordBBox & clipbox)2175 geoActivateBBox(GridType& grid,
2176 const openvdb::CoordBBox &bbox,
2177 bool setvalue,
2178 double value,
2179 GEO_PrimVDB::ActivateOperation operation,
2180 bool doclip,
2181 const openvdb::CoordBBox &clipbox)
2182 {
2183 typename GridType::Accessor access = grid.getAccessor();
2184
2185 switch (operation)
2186 {
2187 case GEO_PrimVDB::ACTIVATE_UNION: // Union
2188 if (doclip)
2189 {
2190 openvdb::CoordBBox clipped = bbox;
2191 clipped = bbox;
2192 clipped.min().maxComponent(clipbox.min());
2193 clipped.max().minComponent(clipbox.max());
2194
2195 geoActivateBBox(grid, clipped, setvalue, value,
2196 operation,
2197 false,
2198 clipped);
2199 break;
2200 }
2201 if (setvalue)
2202 {
2203 grid.fill(bbox, geo_doubleToGridValue<GridType>(value), /*active*/true);
2204 }
2205 else
2206 {
2207 openvdb::MaskGrid mask(false);
2208 mask.denseFill(bbox, true, true);
2209 grid.topologyUnion(mask);
2210 }
2211 break;
2212 case GEO_PrimVDB::ACTIVATE_INTERSECT: // Intersect
2213 {
2214 openvdb::MaskGrid mask(false);
2215 mask.fill(bbox, true, true);
2216 grid.topologyIntersection(mask);
2217 geoInactiveToBackground<GridType> bgop(grid);
2218 openvdb::tools::foreach(grid.beginValueOff(), bgop);
2219 }
2220 break;
2221 case GEO_PrimVDB::ACTIVATE_SUBTRACT: // Difference
2222 // No matter what, we clear the background colour
2223 // for inactive.
2224 grid.fill(bbox, grid.background(), /*active*/false);
2225 break;
2226 case GEO_PrimVDB::ACTIVATE_COPY: // Copy
2227 // intersect
2228 geoActivateBBox(grid, bbox, setvalue, value, GEO_PrimVDB::ACTIVATE_INTERSECT, doclip, clipbox);
2229 // and union
2230 geoActivateBBox(grid, bbox, setvalue, value, GEO_PrimVDB::ACTIVATE_UNION, doclip, clipbox);
2231 break;
2232 }
2233 }
2234
2235 void
activateIndexBBoxAdapter(const void * bboxPtr,ActivateOperation operation,bool setvalue,fpreal value)2236 GEO_PrimVDB::activateIndexBBoxAdapter(const void* bboxPtr,
2237 ActivateOperation operation,
2238 bool setvalue,
2239 fpreal value)
2240 {
2241 using namespace openvdb;
2242
2243 // bboxPtr is assumed to point to an openvdb::vX_Y_Z::CoordBBox, for some
2244 // version X.Y.Z of OpenVDB that may be newer than the one with which
2245 // libHoudiniGEO.so was built. This is safe provided that CoordBBox and
2246 // its member objects are ABI-compatible between the two OpenVDB versions.
2247 const CoordBBox& bbox = *static_cast<const CoordBBox*>(bboxPtr);
2248
2249 bool doclip;
2250 CoordBBox clipbox;
2251 doclip = geoGetFrustumBoundsFromVDB(this, clipbox);
2252
2253 // Activate based on the parameters and inputs
2254 UTvdbCallAllTopology(this->getStorageType(), geoActivateBBox,
2255 this->getGrid(),
2256 bbox,
2257 setvalue,
2258 value,
2259 operation,
2260 doclip, clipbox);
2261 }
2262
2263 // Gets a conservative bounding box that maps to a coordinate
2264 // in index space.
2265 openvdb::CoordBBox
geoMapCoord(const openvdb::CoordBBox & bbox_b,GEO_PrimVolumeXform xform_a,GEO_PrimVolumeXform xform_b)2266 geoMapCoord(const openvdb::CoordBBox& bbox_b,
2267 GEO_PrimVolumeXform xform_a,
2268 GEO_PrimVolumeXform xform_b)
2269 {
2270 using openvdb::Coord;
2271 using openvdb::CoordBBox;
2272 // Get the eight corners of the voxel
2273 Coord x = Coord(bbox_b.extents().x(), 0, 0);
2274 Coord y = Coord(0, bbox_b.extents().y(), 0);
2275 Coord z = Coord(0, 0, bbox_b.extents().z());
2276 Coord m = bbox_b.min();
2277
2278 const Coord corners[] = {
2279 m, m+z, m+y, m+y+z, m+x, m+x+z, m+x+y, m+x+y+z,
2280 };
2281
2282 CoordBBox index_bbox;
2283 for (int i=0; i<8; i++)
2284 {
2285 UT_Vector3 corner = UT_Vector3(corners[i].x(), corners[i].y(), corners[i].z());
2286 UT_Vector3 index = xform_a.toVoxelSpace(xform_b.fromVoxelSpace(corner));
2287 Coord coord(int32(index.x()), int32(index.y()), int32(index.z()));
2288 if (i == 0)
2289 index_bbox = CoordBBox(coord, coord);
2290 else
2291 index_bbox.expand(coord);
2292 }
2293 return index_bbox;
2294 }
2295
2296 openvdb::CoordBBox
geoMapCoord(const openvdb::Coord & coord_b,GEO_PrimVolumeXform xform_a,GEO_PrimVolumeXform xform_b)2297 geoMapCoord(const openvdb::Coord& coord_b,
2298 GEO_PrimVolumeXform xform_a,
2299 GEO_PrimVolumeXform xform_b)
2300 {
2301 const openvdb::CoordBBox bbox_b(coord_b, coord_b + openvdb::Coord(1,1,1));
2302 return geoMapCoord(bbox_b, xform_a, xform_b);
2303 }
2304
2305 /// This class is used as a functor to create a mask for a grid's active
2306 /// region.
2307 template<typename GridType>
2308 class geoMaskTopology
2309 {
2310 public:
2311 typedef typename GridType::ValueOnCIter Iterator;
2312 typedef typename openvdb::MaskGrid::Accessor Accessor;
2313
geoMaskTopology(const GEO_PrimVolumeXform & a,const GEO_PrimVolumeXform & b)2314 geoMaskTopology(const GEO_PrimVolumeXform& a, const GEO_PrimVolumeXform& b)
2315 : xform_a(a), xform_b(b)
2316 {
2317 }
2318
operator ()(const Iterator & iter,Accessor & accessor) const2319 inline void operator()(const Iterator& iter, Accessor& accessor) const
2320 {
2321 openvdb::CoordBBox bbox = geoMapCoord(iter.getBoundingBox(), xform_a,
2322 xform_b);
2323 accessor.getTree()->fill(bbox, true, true);
2324 }
2325
2326 private:
2327 const GEO_PrimVolumeXform& xform_a;
2328 const GEO_PrimVolumeXform& xform_b;
2329 };
2330
2331 /// This class is used as a functor to create a mask for the intersection
2332 /// of two grids.
2333 template<typename GridTypeA, typename GridTypeB>
2334 class geoMaskIntersect
2335 {
2336 public:
2337 typedef typename GridTypeA::ValueOnCIter IteratorA;
2338 typedef typename GridTypeB::ConstAccessor AccessorB;
2339 typedef typename openvdb::MaskGrid::Accessor Accessor;
2340
geoMaskIntersect(const GridTypeB & source,const GEO_PrimVolumeXform & a,const GEO_PrimVolumeXform & b)2341 geoMaskIntersect(const GridTypeB& source,
2342 const GEO_PrimVolumeXform& a,
2343 const GEO_PrimVolumeXform& b)
2344 : myAccessor(source.getAccessor()),
2345 myXformA(a),
2346 myXformB(b)
2347 {
2348 }
2349
operator ()(const IteratorA & iter,Accessor & accessor) const2350 inline void operator()(const IteratorA& iter, Accessor& accessor) const
2351 {
2352 openvdb::CoordBBox bbox = iter.getBoundingBox();
2353 for(int k = bbox.min().z(); k <= bbox.max().z(); k++) {
2354 for (int j = bbox.min().y(); j <= bbox.max().y(); j++) {
2355 for (int i = bbox.min().x(); i <= bbox.max().x(); i++) {
2356 openvdb::Coord coord(i, j, k);
2357 accessor.setActiveState(coord,
2358 containsActiveVoxels(geoMapCoord(coord, myXformB, myXformA)));
2359 }
2360 }
2361 }
2362 }
2363
2364 private:
2365 AccessorB myAccessor;
2366 const GEO_PrimVolumeXform& myXformA;
2367 const GEO_PrimVolumeXform& myXformB;
2368
2369 // Returns true if there is at least one voxel in the source grid that is active
2370 // within the specified bounding box.
containsActiveVoxels(const openvdb::CoordBBox & bbox) const2371 inline bool containsActiveVoxels(const openvdb::CoordBBox& bbox) const
2372 {
2373 for(int k = bbox.min().z(); k <= bbox.max().z(); k++)
2374 {
2375 for(int j = bbox.min().y(); j <= bbox.max().y(); j++)
2376 {
2377 for(int i = bbox.min().x(); i <= bbox.max().x(); i++)
2378 {
2379 if(myAccessor.isValueOn(openvdb::Coord(i, j, k)))
2380 return true;
2381 }
2382 }
2383 }
2384 return false;
2385 }
2386 };
2387
2388 template <typename GridTypeA, typename GridTypeB>
2389 void
geoUnalignedUnion(GridTypeA & grid_a,const GridTypeB & grid_b,GEO_PrimVolumeXform xform_a,GEO_PrimVolumeXform xform_b,bool setvalue,double value,bool doclip,const openvdb::CoordBBox & clipbox)2390 geoUnalignedUnion(GridTypeA &grid_a,
2391 const GridTypeB &grid_b,
2392 GEO_PrimVolumeXform xform_a,
2393 GEO_PrimVolumeXform xform_b,
2394 bool setvalue, double value,
2395 bool doclip, const openvdb::CoordBBox &clipbox)
2396 {
2397 openvdb::MaskGrid mask(false);
2398 geoMaskTopology<GridTypeB> maskop(xform_a, xform_b);
2399 openvdb::tools::transformValues(grid_b.cbeginValueOn(), mask, maskop);
2400 if(doclip)
2401 mask.clip(clipbox);
2402
2403 if(setvalue)
2404 {
2405 typename GridTypeA::TreeType newTree(mask.tree(),
2406 geo_doubleToGridValue<GridTypeA>(value), openvdb::TopologyCopy());
2407 openvdb::tools::compReplace(grid_a.tree(), newTree);
2408 }
2409 else
2410 grid_a.tree().topologyUnion(mask.tree());
2411 }
2412
2413 template <typename GridTypeA, typename GridTypeB>
2414 void
geoUnalignedDifference(GridTypeA & grid_a,const GridTypeB & grid_b,GEO_PrimVolumeXform xform_a,GEO_PrimVolumeXform xform_b)2415 geoUnalignedDifference(GridTypeA &grid_a,
2416 const GridTypeB &grid_b,
2417 GEO_PrimVolumeXform xform_a,
2418 GEO_PrimVolumeXform xform_b)
2419 {
2420 openvdb::MaskGrid mask(false);
2421 geoMaskIntersect<GridTypeA, GridTypeB> maskop(grid_b, xform_a, xform_b);
2422 openvdb::tools::transformValues(grid_a.cbeginValueOn(), mask, maskop, true,
2423 // DO NOT SHARE THE OPERATOR, since grid_b's accessor does caching...
2424 false);
2425
2426 grid_a.tree().topologyDifference(mask.tree());
2427
2428 geoInactiveToBackground<GridTypeA> bgop(grid_a);
2429 openvdb::tools::foreach(grid_a.beginValueOff(), bgop);
2430 }
2431
2432 template <typename GridTypeA, typename GridTypeB>
2433 static void
geoUnalignedIntersect(GridTypeA & grid_a,const GridTypeB & grid_b,GEO_PrimVolumeXform xform_a,GEO_PrimVolumeXform xform_b)2434 geoUnalignedIntersect(GridTypeA &grid_a,
2435 const GridTypeB &grid_b,
2436 GEO_PrimVolumeXform xform_a,
2437 GEO_PrimVolumeXform xform_b)
2438 {
2439 openvdb::MaskGrid mask(false);
2440 geoMaskIntersect<GridTypeA, GridTypeB> maskop(grid_b, xform_a, xform_b);
2441 openvdb::tools::transformValues(grid_a.cbeginValueOn(), mask, maskop, true,
2442 // DO NOT SHARE THE OPERATOR, since grid_b's accessor does caching...
2443 false);
2444
2445 grid_a.tree().topologyIntersection(mask.tree());
2446
2447 geoInactiveToBackground<GridTypeA> bgop(grid_a);
2448 openvdb::tools::foreach(grid_a.beginValueOff(), bgop);
2449 }
2450
2451 // The result of the union of active regions goes into grid_a
2452 template <typename GridTypeA, typename GridTypeB>
2453 static void
geoUnion(GridTypeA & grid_a,const GridTypeB & grid_b,bool setvalue,double value,bool doclip,const openvdb::CoordBBox & clipbox)2454 geoUnion(GridTypeA& grid_a, const GridTypeB &grid_b, bool setvalue, double value, bool doclip, const openvdb::CoordBBox &clipbox)
2455 {
2456 typename GridTypeA::Accessor access_a = grid_a.getAccessor();
2457 typename GridTypeB::ConstAccessor access_b = grid_b.getAccessor();
2458
2459 if (!doclip && !setvalue) {
2460 grid_a.tree().topologyUnion(grid_b.tree());
2461 return;
2462 }
2463
2464 // For each on value in b, set a on
2465 for (typename GridTypeB::ValueOnCIter iter = grid_b.cbeginValueOn(); iter; ++iter) {
2466 openvdb::CoordBBox bbox = iter.getBoundingBox();
2467 // Intersect with our destination
2468 if (doclip) {
2469 bbox.min().maxComponent(clipbox.min());
2470 bbox.max().minComponent(clipbox.max());
2471 }
2472
2473 for (int k=bbox.min().z(); k<=bbox.max().z(); k++) {
2474 for (int j=bbox.min().y(); j<=bbox.max().y(); j++) {
2475 for (int i=bbox.min().x(); i<=bbox.max().x(); i++) {
2476 openvdb::Coord coord(i, j, k);
2477 if (setvalue) {
2478 access_a.setValue(coord, geo_doubleToGridValue<GridTypeA>(value));
2479 } else {
2480 access_a.setValueOn(coord);
2481 }
2482 }
2483 }
2484 }
2485 }
2486 }
2487
2488 // The result of the union of active regions goes into grid_a
2489 template <typename GridTypeA, typename GridTypeB>
2490 static void
geoDifference(GridTypeA & grid_a,const GridTypeB & grid_b)2491 geoDifference(GridTypeA& grid_a, const GridTypeB &grid_b)
2492 {
2493 typename GridTypeA::Accessor access_a = grid_a.getAccessor();
2494 typename GridTypeB::ConstAccessor access_b = grid_b.getAccessor();
2495
2496 // For each on value in a, set it off if b is on
2497 for (typename GridTypeA::ValueOnCIter
2498 iter = grid_a.cbeginValueOn(); iter; ++iter)
2499 {
2500 openvdb::CoordBBox bbox = iter.getBoundingBox();
2501 for (int k=bbox.min().z(); k<=bbox.max().z(); k++)
2502 {
2503 for (int j=bbox.min().y(); j<=bbox.max().y(); j++)
2504 {
2505 for (int i=bbox.min().x(); i<=bbox.max().x(); i++)
2506 {
2507 openvdb::Coord coord(i, j, k);
2508 // TODO: conditional needed? Profile please.
2509 if (access_b.isValueOn(coord))
2510 {
2511 access_a.setValue(coord, grid_a.background());
2512 access_a.setValueOff(coord);
2513 }
2514 }
2515 }
2516 }
2517 }
2518 }
2519
2520 template <typename GridTypeB>
2521 static void
geoDoUnion(const GridTypeB & grid_b,GEO_PrimVolumeXform xform_b,GEO_PrimVDB & vdb_a,bool setvalue,double value,bool doclip,const openvdb::CoordBBox & clipbox,bool ignore_transform)2522 geoDoUnion(const GridTypeB &grid_b,
2523 GEO_PrimVolumeXform xform_b,
2524 GEO_PrimVDB &vdb_a,
2525 bool setvalue, double value,
2526 bool doclip, const openvdb::CoordBBox &clipbox,
2527 bool ignore_transform)
2528 {
2529 // If the transforms are equal, we can do an aligned union
2530 if (ignore_transform || grid_b.transform() == vdb_a.getGrid().transform())
2531 {
2532 UTvdbCallAllTopology(vdb_a.getStorageType(), geoUnion,
2533 vdb_a.getGrid(), grid_b, setvalue, value,
2534 doclip, clipbox);
2535 }
2536 else
2537 {
2538 UTvdbCallAllTopology(vdb_a.getStorageType(), geoUnalignedUnion,
2539 vdb_a.getGrid(), grid_b,
2540 vdb_a.getIndexSpaceTransform(),
2541 xform_b, setvalue, value,
2542 doclip, clipbox);
2543 }
2544 }
2545
2546 template <typename GridTypeB>
2547 static void
geoDoIntersect(const GridTypeB & grid_b,GEO_PrimVolumeXform xform_b,GEO_PrimVDB & vdb_a,bool ignore_transform)2548 geoDoIntersect(
2549 const GridTypeB &grid_b,
2550 GEO_PrimVolumeXform xform_b,
2551 GEO_PrimVDB &vdb_a,
2552 bool ignore_transform)
2553 {
2554 if (ignore_transform || grid_b.transform() == vdb_a.getGrid().transform())
2555 {
2556 UTvdbCallAllTopology(vdb_a.getStorageType(),
2557 geoIntersect, vdb_a.getGrid(), grid_b);
2558 }
2559 else
2560 {
2561 UTvdbCallAllTopology(vdb_a.getStorageType(),
2562 geoUnalignedIntersect, vdb_a.getGrid(),
2563 grid_b, vdb_a.getIndexSpaceTransform(),
2564 xform_b);
2565 }
2566 }
2567
2568 template <typename GridTypeB>
2569 static void
geoDoDifference(const GridTypeB & grid_b,GEO_PrimVolumeXform xform_b,GEO_PrimVDB & vdb_a,bool ignore_transform)2570 geoDoDifference(
2571 const GridTypeB &grid_b,
2572 GEO_PrimVolumeXform xform_b,
2573 GEO_PrimVDB &vdb_a,
2574 bool ignore_transform)
2575 {
2576 if (ignore_transform || grid_b.transform() == vdb_a.getGrid().transform())
2577 {
2578 UTvdbCallAllTopology(vdb_a.getStorageType(), geoDifference,
2579 vdb_a.getGrid(), grid_b);
2580 }
2581 else
2582 {
2583 UTvdbCallAllTopology(vdb_a.getStorageType(), geoUnalignedDifference,
2584 vdb_a.getGrid(), grid_b,
2585 vdb_a.getIndexSpaceTransform(),
2586 xform_b);
2587 }
2588 }
2589
2590
2591 void
activateByVDB(const GEO_PrimVDB * input_vdb,ActivateOperation operation,bool setvalue,fpreal value,bool ignore_transform)2592 GEO_PrimVDB::activateByVDB(
2593 const GEO_PrimVDB *input_vdb,
2594 ActivateOperation operation,
2595 bool setvalue, fpreal value,
2596 bool ignore_transform)
2597 {
2598 const openvdb::GridBase& input_grid = input_vdb->getGrid();
2599
2600 bool doclip;
2601 openvdb::CoordBBox clipbox;
2602
2603 doclip = geoGetFrustumBoundsFromVDB(this, clipbox);
2604
2605 switch (operation)
2606 {
2607 case GEO_PrimVDB::ACTIVATE_UNION: // Union
2608 UTvdbCallAllTopology(input_vdb->getStorageType(),
2609 geoDoUnion, input_grid,
2610 input_vdb->getIndexSpaceTransform(),
2611 *this,
2612 setvalue,
2613 value,
2614 doclip, clipbox,
2615 ignore_transform);
2616 break;
2617 case GEO_PrimVDB::ACTIVATE_INTERSECT: // Intersect
2618 UTvdbCallAllTopology(input_vdb->getStorageType(),
2619 geoDoIntersect, input_grid,
2620 input_vdb->getIndexSpaceTransform(),
2621 *this,
2622 ignore_transform);
2623 break;
2624 case GEO_PrimVDB::ACTIVATE_SUBTRACT: // Difference
2625 UTvdbCallAllTopology(input_vdb->getStorageType(),
2626 geoDoDifference, input_grid,
2627 input_vdb->getIndexSpaceTransform(),
2628 *this,
2629 ignore_transform);
2630 break;
2631 case GEO_PrimVDB::ACTIVATE_COPY: // Copy
2632 UTvdbCallAllTopology(input_vdb->getStorageType(),
2633 geoDoIntersect, input_grid,
2634 input_vdb->getIndexSpaceTransform(),
2635 *this,
2636 ignore_transform);
2637 UTvdbCallAllTopology(input_vdb->getStorageType(),
2638 geoDoUnion, input_grid,
2639 input_vdb->getIndexSpaceTransform(),
2640 *this,
2641 setvalue,
2642 value,
2643 doclip, clipbox,
2644 ignore_transform);
2645 break;
2646 }
2647 }
2648
2649 UT_Matrix4D
getTransform4() const2650 GEO_PrimVDB::getTransform4() const
2651 {
2652 using namespace openvdb;
2653 using namespace openvdb::math;
2654
2655 UT_Matrix4D mat4;
2656 const Transform &gxform = getGrid().transform();
2657 NonlinearFrustumMap::ConstPtr fmap = gxform.map<NonlinearFrustumMap>();
2658 if (fmap)
2659 {
2660 const openvdb::BBoxd &bbox = fmap->getBBox();
2661 const openvdb::Vec3d center = bbox.getCenter();
2662 const openvdb::Vec3d size = bbox.extents();
2663
2664 // TODO: Use fmap->linearMap() once that actually works
2665 mat4.identity();
2666 mat4.translate(-center.x(), -center.y(), -bbox.min().z());
2667 // NOTE: We scale both XY axes by size.x() because the secondMap()
2668 // has the aspect ratio baked in
2669 mat4.scale(1.0/size.x(), 1.0/size.x(), 1.0/size.z());
2670 mat4 *= UTvdbConvert(fmap->secondMap().getMat4());
2671 }
2672 else
2673 {
2674 mat4 = UTvdbConvert(gxform.baseMap()->getAffineMap()->getMat4());
2675 }
2676 return mat4;
2677 }
2678
2679 void
getLocalTransform(UT_Matrix3D & result) const2680 GEO_PrimVDB::getLocalTransform(UT_Matrix3D &result) const
2681 {
2682 result = getTransform4();
2683 }
2684
2685 void
setLocalTransform(const UT_Matrix3D & new_mat3)2686 GEO_PrimVDB::setLocalTransform(const UT_Matrix3D &new_mat3)
2687 {
2688 using namespace openvdb;
2689 using namespace openvdb::math;
2690
2691 Transform::Ptr xform;
2692 UT_Matrix4D new_mat4;
2693 new_mat4 = new_mat3;
2694 new_mat4.setTranslates(getDetail().getPos3(vertexPoint(0)));
2695
2696 const Transform & gxform = getGrid().transform();
2697 NonlinearFrustumMap::ConstPtr fmap = gxform.map<NonlinearFrustumMap>();
2698 if (fmap)
2699 {
2700 fmap = geoStandardFrustumMapPtr(*this);
2701 const openvdb::BBoxd &bbox = fmap->getBBox();
2702 const openvdb::Vec3d center = bbox.getCenter();
2703 const openvdb::Vec3d size = bbox.extents();
2704
2705 // TODO: Use fmap->linearMap() once that actually works
2706 UT_Matrix4D second;
2707 second.identity();
2708 second.translate(-0.5, -0.5, 0.0); // adjust for frustum map center
2709 // NOTE: We scale both XY axes by size.x() because the secondMap()
2710 // has the aspect ratio baked in
2711 second.scale(size.x(), size.x(), size.z());
2712 second.translate(center.x(), center.y(), bbox.min().z());
2713 second *= new_mat4;
2714 xform.reset(new Transform(MapBase::Ptr(
2715 new NonlinearFrustumMap(fmap->getBBox(), fmap->getTaper(),
2716 /*depth*/1.0,
2717 geoCreateAffineMap<MapBase>(second)))));
2718 }
2719 else
2720 {
2721 xform = geoCreateLinearTransform(new_mat4);
2722 }
2723 myGridAccessor.setTransform(*xform, *this);
2724 }
2725
2726 int
detachPoints(GA_PointGroup & grp)2727 GEO_PrimVDB::detachPoints(GA_PointGroup &grp)
2728 {
2729 int count = 0;
2730
2731 if (grp.containsOffset(vertexPoint(0)))
2732 count++;
2733
2734 if (count == 0)
2735 return 0;
2736
2737 if (count == 1)
2738 return -2;
2739
2740 return -1;
2741 }
2742
2743 GA_Primitive::GA_DereferenceStatus
dereferencePoint(GA_Offset point,bool)2744 GEO_PrimVDB::dereferencePoint(GA_Offset point, bool)
2745 {
2746 return vertexPoint(0) == point
2747 ? GA_DEREFERENCE_DESTROY
2748 : GA_DEREFERENCE_OK;
2749 }
2750
2751 GA_Primitive::GA_DereferenceStatus
dereferencePoints(const GA_RangeMemberQuery & point_query,bool)2752 GEO_PrimVDB::dereferencePoints(const GA_RangeMemberQuery &point_query, bool)
2753 {
2754 return point_query.contains(vertexPoint(0))
2755 ? GA_DEREFERENCE_DESTROY
2756 : GA_DEREFERENCE_OK;
2757 }
2758
2759 ///
2760 /// JSON methods
2761 ///
2762
2763 namespace { // unnamed
2764
2765 class geo_PrimVDBJSON : public GA_PrimitiveJSON
2766 {
2767 public:
2768 static const char *theSharedMemKey;
2769
2770 public:
geo_PrimVDBJSON()2771 geo_PrimVDBJSON() {}
~geo_PrimVDBJSON()2772 ~geo_PrimVDBJSON() override {}
2773
2774 enum
2775 {
2776 geo_TBJ_VERTEX,
2777 geo_TBJ_VDB,
2778 geo_TBJ_VDB_SHMEM,
2779 geo_TBJ_VDB_VISUALIZATION,
2780 geo_TBJ_ENTRIES
2781 };
2782
vdb(const GA_Primitive * p) const2783 const GEO_PrimVDB *vdb(const GA_Primitive *p) const
2784 { return static_cast<const GEO_PrimVDB *>(p); }
vdb(GA_Primitive * p) const2785 GEO_PrimVDB *vdb(GA_Primitive *p) const
2786 { return static_cast<GEO_PrimVDB *>(p); }
2787
getEntries() const2788 int getEntries() const override
2789 { return geo_TBJ_ENTRIES; }
2790
2791 const UT_StringHolder &
getKeyword(int i) const2792 getKeyword(int i) const override
2793 {
2794 switch (i)
2795 {
2796 case geo_TBJ_VERTEX: return theKWVertex;
2797 case geo_TBJ_VDB: return theKWVDB;
2798 case geo_TBJ_VDB_SHMEM: return theKWVDBShm;
2799 case geo_TBJ_VDB_VISUALIZATION: return theKWVDBVis;
2800 case geo_TBJ_ENTRIES: break;
2801 }
2802 UT_ASSERT(0);
2803 return UT_StringHolder::theEmptyString;
2804 }
2805
2806 bool
shouldSaveField(const GA_Primitive * prim,int i,const GA_SaveMap & sm) const2807 shouldSaveField(const GA_Primitive *prim, int i,
2808 const GA_SaveMap &sm) const override
2809 {
2810 bool is_shmem = sm.getOptions().hasOption("geo:sharedmemowner");
2811 switch (i)
2812 {
2813 case geo_TBJ_VERTEX: return true;
2814 case geo_TBJ_VDB: return !is_shmem;
2815 case geo_TBJ_VDB_SHMEM: return is_shmem;
2816 case geo_TBJ_VDB_VISUALIZATION: return true;
2817 case geo_TBJ_ENTRIES: break;
2818 }
2819 UT_ASSERT(0);
2820 return false;
2821 }
2822
2823 bool
saveField(const GA_Primitive * pr,int i,UT_JSONWriter & w,const GA_SaveMap & map) const2824 saveField(const GA_Primitive *pr, int i, UT_JSONWriter &w,
2825 const GA_SaveMap &map) const override
2826 {
2827 switch (i)
2828 {
2829 case geo_TBJ_VERTEX:
2830 {
2831 GA_Offset vtx = vdb(pr)->getVertexOffset(0);
2832 return w.jsonInt(int64(map.getVertexIndex(vtx)));
2833 }
2834 case geo_TBJ_VDB:
2835 return vdb(pr)->saveVDB(w, map);
2836 case geo_TBJ_VDB_SHMEM:
2837 return vdb(pr)->saveVDB(w, map, true);
2838 case geo_TBJ_VDB_VISUALIZATION:
2839 return vdb(pr)->saveVisualization(w, map);
2840
2841 case geo_TBJ_ENTRIES:
2842 break;
2843 }
2844 return false;
2845 }
2846 bool
loadField(GA_Primitive * pr,int i,UT_JSONParser & p,const GA_LoadMap & map) const2847 loadField(GA_Primitive *pr, int i, UT_JSONParser &p,
2848 const GA_LoadMap &map) const override
2849 {
2850 switch (i)
2851 {
2852 case geo_TBJ_VERTEX:
2853 {
2854 int64 vidx;
2855 if (!p.parseInt(vidx))
2856 return false;
2857 GA_Offset voff = map.getVertexOffset(GA_Index(vidx));
2858 // Assign the preallocated vertex, but
2859 // do not bother updating the topology,
2860 // which will be done at the end of the
2861 // load anyway.
2862 vdb(pr)->assignVertex(voff, false);
2863 return true;
2864 }
2865 case geo_TBJ_VDB:
2866 return vdb(pr)->loadVDB(p);
2867 case geo_TBJ_VDB_SHMEM:
2868 return vdb(pr)->loadVDB(p, true);
2869 case geo_TBJ_VDB_VISUALIZATION:
2870 return vdb(pr)->loadVisualization(p, map);
2871
2872 case geo_TBJ_ENTRIES:
2873 break;
2874 }
2875 UT_ASSERT(0);
2876 return false;
2877 }
2878
2879 bool
saveField(const GA_Primitive * pr,int i,UT_JSONValue & val,const GA_SaveMap & map) const2880 saveField(const GA_Primitive *pr, int i, UT_JSONValue &val,
2881 const GA_SaveMap &map) const override
2882 {
2883 UT_AutoJSONWriter w(val);
2884 return saveField(pr, i, *w, map);
2885 }
2886 // Re-implement the H12.5 base class version, note that this was pure
2887 // virtual in H12.1.
2888 bool
loadField(GA_Primitive * pr,int i,UT_JSONParser & p,const UT_JSONValue & jval,const GA_LoadMap & map) const2889 loadField(GA_Primitive *pr, int i, UT_JSONParser &p,
2890 const UT_JSONValue &jval, const GA_LoadMap &map) const override
2891 {
2892 UT_AutoJSONParser parser(jval);
2893 bool ok = loadField(pr, i, *parser, map);
2894 p.stealErrors(*parser);
2895 return ok;
2896 }
2897
2898 bool
isEqual(int i,const GA_Primitive * p0,const GA_Primitive * p1) const2899 isEqual(int i,
2900 const GA_Primitive *p0, const GA_Primitive *p1) const override
2901 {
2902 switch (i)
2903 {
2904 case geo_TBJ_VERTEX:
2905 return (p0->getVertexOffset(0) == p1->getVertexOffset(0));
2906 case geo_TBJ_VDB_SHMEM:
2907 case geo_TBJ_VDB:
2908 return false; // never save these tags as uniform
2909 case geo_TBJ_VDB_VISUALIZATION:
2910 return (vdb(p0)->getVisOptions() == vdb(p1)->getVisOptions());
2911 case geo_TBJ_ENTRIES:
2912 break;
2913 }
2914 UT_ASSERT(0);
2915 return false;
2916 }
2917
2918 private:
2919 };
2920
2921 const char *geo_PrimVDBJSON::theSharedMemKey = "sharedkey";
2922
2923 } // namespace unnamed
2924
2925
2926 static const GA_PrimitiveJSON *
vdbJSON()2927 vdbJSON()
2928 {
2929 static SYS_AtomicPtr<GA_PrimitiveJSON> theJSON;
2930
2931 if (!theJSON) {
2932 GA_PrimitiveJSON* json = new geo_PrimVDBJSON;
2933 if (nullptr != theJSON.compare_swap(nullptr, json)) {
2934 delete json;
2935 json = nullptr;
2936 }
2937 }
2938 return theJSON;
2939 }
2940
2941 const GA_PrimitiveJSON *
getJSON() const2942 GEO_PrimVDB::getJSON() const
2943 {
2944 return vdbJSON();
2945 }
2946
2947 // This method is called by multiple places internally in Houdini.
2948 static void
geoSetVDBStreamCompression(openvdb::io::Stream & vos,bool backwards_compatible)2949 geoSetVDBStreamCompression(openvdb::io::Stream& vos, bool backwards_compatible)
2950 {
2951 // Always enable full compression, since it is fast and compresses level
2952 // sets and fog volumes well.
2953 uint32_t compression = openvdb::io::COMPRESS_ACTIVE_MASK;
2954 // Enable blosc compression unless we want it to be backwards compatible.
2955 if (vos.hasBloscCompression() && !backwards_compatible) {
2956 compression |= openvdb::io::COMPRESS_BLOSC;
2957 }
2958 vos.setCompression(compression);
2959 }
2960
2961 bool
saveVDB(UT_JSONWriter & w,const GA_SaveMap & sm,bool as_shmem) const2962 GEO_PrimVDB::saveVDB(UT_JSONWriter &w, const GA_SaveMap &sm,
2963 bool as_shmem) const
2964 {
2965 bool ok = true;
2966
2967 try
2968 {
2969 openvdb::GridCPtrVec grids;
2970 grids.push_back(getConstGridPtr());
2971
2972 if (as_shmem)
2973 {
2974 openvdb::MetaMap meta;
2975
2976 UT_String shmem_owner;
2977 sm.getOptions().importOption("geo:sharedmemowner", shmem_owner);
2978 if (!shmem_owner.isstring())
2979 return false;
2980
2981 // First do a pass to collect the final size
2982 SYS_SharedMemoryOutputStream os_count(NULL);
2983 {
2984 openvdb::io::Stream vos(os_count);
2985 geoSetVDBStreamCompression(vos, /*backwards_compatible*/false);
2986 vos.write(grids, meta);
2987 }
2988
2989 // Create the shmem segment
2990 UT_WorkBuffer shmem_key;
2991 shmem_key.sprintf("%s:%p", shmem_owner.buffer(), this);
2992 UT_SharedMemoryManager &shmgr = UT_SharedMemoryManager::get();
2993 SYS_SharedMemory *shmem = shmgr.get(shmem_key.buffer());
2994 if (shmem->size() != os_count.size())
2995 shmem->reset(os_count.size());
2996
2997 // Save the vdb stream to the shmem segment
2998 SYS_SharedMemoryOutputStream os_shm(shmem);
2999 {
3000 openvdb::io::Stream vos(os_shm);
3001 geoSetVDBStreamCompression(vos, /*backwards_compatible*/false);
3002 vos.write(grids, meta);
3003 }
3004
3005 // In the main json stream, just tag it with the shmem key
3006 ok = ok && w.jsonBeginArray();
3007 ok = ok && w.jsonKeyToken(geo_PrimVDBJSON::theSharedMemKey);
3008 ok = ok && w.jsonString(shmem->id());
3009 ok = ok && w.jsonEndArray();
3010 }
3011 else
3012 {
3013 UT_JSONWriter::TiledStream os(w);
3014
3015 openvdb::io::Stream vos(os);
3016 openvdb::MetaMap meta;
3017
3018 geoSetVDBStreamCompression(
3019 vos, UT_EnvControl::getInt(ENV_HOUDINI13_VOLUME_COMPATIBILITY));
3020
3021 // Visual C++ requires a default meta object declared on the stack
3022 vos.write(grids, meta);
3023 }
3024 }
3025 catch (std::exception &e)
3026 {
3027 std::cerr << "Save failure: " << e.what() << "\n";
3028 ok = false;
3029 }
3030 return ok;
3031 }
3032
3033 bool
loadVDB(UT_JSONParser & p,bool as_shmem)3034 GEO_PrimVDB::loadVDB(UT_JSONParser &p, bool as_shmem)
3035 {
3036 if (as_shmem)
3037 {
3038 bool array_error = false;
3039 UT_WorkBuffer key;
3040
3041 if (!p.parseBeginArray(array_error) || array_error)
3042 return false;
3043
3044 if (!p.parseString(key))
3045 return false;
3046
3047 if (key != geo_PrimVDBJSON::theSharedMemKey)
3048 return false;
3049
3050 UT_WorkBuffer shmem_key;
3051 if (!p.parseString(shmem_key))
3052 return false;
3053
3054 SYS_SharedMemory *shmem = new SYS_SharedMemory(shmem_key.buffer(),
3055 /*read_only*/true);
3056
3057 if (shmem->size())
3058 {
3059 try
3060 {
3061 SYS_SharedMemoryInputStream is_shm(*shmem);
3062 openvdb::io::Stream vis(is_shm, /*delayLoad*/false);
3063 openvdb::GridPtrVecPtr grids = vis.getGrids();
3064
3065 int count = (grids ? grids->size() : 0);
3066 if (count != 1)
3067 {
3068 UT_String mesg;
3069 mesg.sprintf("expected to read 1 grid, got %d grid%s",
3070 count, count == 1 ? "" : "s");
3071 throw std::runtime_error(mesg.nonNullBuffer());
3072 }
3073
3074 openvdb::GridBase::Ptr grid = (*grids)[0];
3075 UT_ASSERT(grid);
3076 if (grid) setGrid(*grid);
3077 }
3078 catch (std::exception &e)
3079 {
3080 std::cerr << "Shared memory load failure: " << e.what() << "\n";
3081 return false;
3082 }
3083 }
3084 else
3085 {
3086 // If the shared memory was set to zero, it probably died while
3087 // the IFD stream was in transit. Create a dummy grid so that
3088 // mantra doesn't flip out like a ninja.
3089 openvdb::GridBase::Ptr grid = openvdb::FloatGrid::create(0);
3090 setGrid(*grid);
3091 }
3092
3093 if (!p.parseEndArray(array_error) || array_error)
3094 return false;
3095 }
3096 else
3097 {
3098 try
3099 {
3100 UT_JSONParser::TiledStream is(p);
3101
3102 openvdb::io::Stream vis(is, /*delayLoad*/false);
3103
3104 openvdb::GridPtrVecPtr grids = vis.getGrids();
3105
3106 int count = (grids ? grids->size() : 0);
3107 if (count != 1)
3108 {
3109 UT_String mesg;
3110 mesg.sprintf("expected to read 1 grid, got %d grid%s",
3111 count, count == 1 ? "" : "s");
3112 throw std::runtime_error(mesg.nonNullBuffer());
3113 }
3114
3115 openvdb::GridBase::Ptr grid = (*grids)[0];
3116 UT_ASSERT(grid);
3117 if (grid)
3118 {
3119 // When we saved the grid, we auto-added metadata
3120 // which isn't reflected by our primitive attributes.
3121 // if any later node tries to sync the metadata from
3122 // the vdb primitive, we'll gain extra data such as
3123 // file_bbox
3124 const char *file_metadata[] =
3125 {
3126 "file_bbox_min",
3127 "file_bbox_max",
3128 "file_compression",
3129 "file_mem_bytes",
3130 "file_voxel_count",
3131 "file_delayed_load",
3132 0
3133 };
3134 for (int i = 0; file_metadata[i]; i++)
3135 {
3136 grid->removeMeta(file_metadata[i]);
3137 }
3138 setGrid(*grid);
3139 }
3140 }
3141 catch (std::exception &e)
3142 {
3143 std::cerr << "Load failure: " << e.what() << "\n";
3144 return false;
3145 }
3146 }
3147 return true;
3148 }
3149
3150 namespace // anonymous
3151 {
3152
3153 enum
3154 {
3155 geo_JVOL_VISMODE,
3156 geo_JVOL_VISISO,
3157 geo_JVOL_VISDENSITY,
3158 geo_JVOL_VISLOD,
3159 };
3160 UT_FSATable theJVolumeViz(
3161 geo_JVOL_VISMODE, "mode",
3162 geo_JVOL_VISISO, "iso",
3163 geo_JVOL_VISDENSITY, "density",
3164 geo_JVOL_VISLOD, "lod",
3165 -1, nullptr
3166 );
3167
3168 } // namespace anonymous
3169
3170 bool
saveVisualization(UT_JSONWriter & w,const GA_SaveMap &) const3171 GEO_PrimVDB::saveVisualization(UT_JSONWriter &w, const GA_SaveMap &) const
3172 {
3173 bool ok = true;
3174 ok = ok && w.jsonBeginMap();
3175
3176 ok = ok && w.jsonKeyToken(theJVolumeViz.getToken(geo_JVOL_VISMODE));
3177 ok = ok && w.jsonString(GEOgetVolumeVisToken(myVis.myMode));
3178
3179 ok = ok && w.jsonKeyToken(theJVolumeViz.getToken(geo_JVOL_VISISO));
3180 ok = ok && w.jsonReal(myVis.myIso);
3181
3182 ok = ok && w.jsonKeyToken(theJVolumeViz.getToken(geo_JVOL_VISDENSITY));
3183 ok = ok && w.jsonReal(myVis.myDensity);
3184
3185 // Only save myLod when non-default so that it loads in older builds
3186 if (myVis.myLod != GEO_VOLUMEVISLOD_FULL)
3187 {
3188 ok = ok && w.jsonKeyToken(theJVolumeViz.getToken(geo_JVOL_VISLOD));
3189 ok = ok && w.jsonString(GEOgetVolumeVisLodToken(myVis.myLod));
3190 }
3191
3192 return ok && w.jsonEndMap();
3193 }
3194
3195 bool
loadVisualization(UT_JSONParser & p,const GA_LoadMap &)3196 GEO_PrimVDB::loadVisualization(UT_JSONParser &p, const GA_LoadMap &)
3197 {
3198 UT_JSONParser::traverser it;
3199 GEO_VolumeVis mode = myVis.myMode;
3200 fpreal iso = myVis.myIso;
3201 fpreal density = myVis.myDensity;
3202 GEO_VolumeVisLod lod = myVis.myLod;
3203 UT_WorkBuffer key;
3204 fpreal64 fval;
3205 bool foundmap=false, ok = true;
3206
3207 for (it = p.beginMap(); ok && !it.atEnd(); ++it)
3208 {
3209 foundmap = true;
3210 if (!it.getLowerKey(key))
3211 {
3212 ok = false;
3213 break;
3214 }
3215 switch (theJVolumeViz.findSymbol(key.buffer()))
3216 {
3217 case geo_JVOL_VISMODE:
3218 if ((ok = p.parseString(key)))
3219 mode = GEOgetVolumeVisEnum(
3220 key.buffer(), GEO_VOLUMEVIS_SMOKE);
3221 break;
3222 case geo_JVOL_VISISO:
3223 if ((ok = p.parseReal(fval)))
3224 iso = fval;
3225 break;
3226 case geo_JVOL_VISDENSITY:
3227 if ((ok = p.parseReal(fval)))
3228 density = fval;
3229 break;
3230 case geo_JVOL_VISLOD:
3231 if ((ok = p.parseString(key)))
3232 lod = GEOgetVolumeVisLodEnum(
3233 key.buffer(), GEO_VOLUMEVISLOD_FULL);
3234 break;
3235 default:
3236 p.addWarning("Unexpected key for volume visualization: %s",
3237 key.buffer());
3238 ok = p.skipNextObject();
3239 break;
3240 }
3241 }
3242 if (!foundmap)
3243 {
3244 p.addFatal("Expected a JSON map for volume visualization data");
3245 ok = false;
3246 }
3247 if (ok)
3248 setVisualization(mode, iso, density, lod);
3249 return ok;
3250 }
3251
3252 template <typename GridType>
3253 static void
geo_sumPosDensity(const GridType & grid,fpreal64 & sum)3254 geo_sumPosDensity(const GridType &grid, fpreal64 &sum)
3255 {
3256 sum = 0;
3257 for (typename GridType::ValueOnCIter iter = grid.cbeginValueOn(); iter; ++iter) {
3258 fpreal value = *iter;
3259 if (value > 0) {
3260 if (iter.isTileValue()) sum += value * iter.getVoxelCount();
3261 else sum += value;
3262 }
3263 }
3264 }
3265
3266 fpreal
calcPositiveDensity() const3267 GEO_PrimVDB::calcPositiveDensity() const
3268 {
3269 fpreal64 density = 0;
3270
3271 UT_IF_ASSERT(UT_VDBType type = getStorageType();)
3272 UT_ASSERT(type == UT_VDB_FLOAT || type == UT_VDB_DOUBLE);
3273
3274 UTvdbCallRealType(getStorageType(), geo_sumPosDensity, getGrid(), density);
3275 UTvdbCallBoolType(getStorageType(), geo_sumPosDensity, getGrid(), density);
3276
3277 int numvoxel = getGrid().activeVoxelCount();
3278 if (numvoxel)
3279 density /= numvoxel;
3280
3281 UT_Vector3 zero(0, 0, 0);
3282 density *= calcVolume(zero);
3283
3284 return density;
3285 }
3286
3287 int
getBBox(UT_BoundingBox * bbox) const3288 GEO_PrimVDB::getBBox(UT_BoundingBox *bbox) const
3289 {
3290 if (hasGrid())
3291 {
3292 using namespace openvdb;
3293
3294 CoordBBox vbox;
3295
3296 const openvdb::GridBase &grid = getGrid();
3297 // NOTE: We use evalActiveVoxelBoundingBox() so that it matches
3298 // getRes() which calls evalActiveVoxelDim().
3299 if (!grid.baseTree().evalActiveVoxelBoundingBox(vbox))
3300 {
3301 bbox->makeInvalid();
3302 return false;
3303 }
3304 // Currently VDB may return true even if the final bounds
3305 // were zero, so to avoid generating a massive bound, detect
3306 // invalid bounding boxes and return false.
3307 if (vbox.min()[0] > vbox.max()[0])
3308 {
3309 bbox->makeInvalid();
3310 return false;
3311 }
3312
3313 const math::Transform &xform = grid.transform();
3314
3315 for (int i = 0; i < 8; i++)
3316 {
3317 math::Vec3d vpos(
3318 (i&1) ? vbox.min()[0] - 0.5 : vbox.max()[0] + 0.5,
3319 (i&2) ? vbox.min()[1] - 0.5 : vbox.max()[1] + 0.5,
3320 (i&4) ? vbox.min()[2] - 0.5 : vbox.max()[2] + 0.5);
3321 vpos = xform.indexToWorld(vpos);
3322
3323 UT_Vector3 worldpos(vpos.x(), vpos.y(), vpos.z());
3324 if (i == 0)
3325 bbox->initBounds(worldpos);
3326 else
3327 bbox->enlargeBounds(worldpos);
3328 }
3329 return true;
3330 }
3331
3332 bbox->initBounds(getDetail().getPos3(vertexPoint(0)));
3333
3334 return true;
3335 }
3336
3337 UT_Vector3
baryCenter() const3338 GEO_PrimVDB::baryCenter() const
3339 {
3340 // Return the center of the index space
3341 if (!hasGrid())
3342 return UT_Vector3(0, 0, 0);
3343
3344 const openvdb::GridBase &grid = getGrid();
3345 openvdb::CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
3346 UT_Vector3 pos;
3347 findexToPos(UTvdbConvert(bbox.getCenter()), pos);
3348 return pos;
3349 }
3350
3351 bool
isDegenerate() const3352 GEO_PrimVDB::isDegenerate() const
3353 {
3354 return false;
3355 }
3356
3357 //
3358 // Methods to handle vertex attributes for the attribute dictionary
3359 //
3360 void
copyPrimitive(const GEO_Primitive * psrc)3361 GEO_PrimVDB::copyPrimitive(const GEO_Primitive *psrc)
3362 {
3363 if (psrc == this) return;
3364
3365 // Ensure the vertices are copied.
3366 GEO_Primitive::copyPrimitive(psrc);
3367
3368 const GEO_PrimVDB *src = (const GEO_PrimVDB *)psrc;
3369
3370 copyGridFrom(*src); // makes a shallow copy
3371
3372 myVis = src->myVis;
3373 }
3374
3375 static inline
3376 openvdb::math::Vec3d
vdbTranslation(const openvdb::math::Transform & xform)3377 vdbTranslation(const openvdb::math::Transform &xform)
3378 {
3379 return xform.baseMap()->getAffineMap()->getMat4().getTranslation();
3380 }
3381
3382 // Replace the grid's translation with the prim's vertex position
3383 void
updateGridTranslates(const GEO_PrimVDB & prim) const3384 GEO_PrimVDB::GridAccessor::updateGridTranslates(const GEO_PrimVDB &prim) const
3385 {
3386 using namespace openvdb::math;
3387 const GA_Detail & geo = prim.getDetail();
3388
3389 // It is possible our vertex offset is invalid, such as us
3390 // being a stashed primitive.
3391 if (!GAisValid(prim.getVertexOffset(0)))
3392 return;
3393
3394 GA_Offset ptoff = prim.vertexPoint(0);
3395 Vec3d newpos = UTvdbConvert(geo.getPos3(ptoff));
3396 Vec3d oldpos = vdbTranslation(myGrid->transform());
3397 MapBase::ConstPtr map = myGrid->transform().baseMap();
3398
3399 if (isApproxEqual(oldpos, newpos))
3400 return;
3401
3402 const_cast<GEO_PrimVDB&>(prim).incrTransformUniqueId();
3403 Vec3d delta = newpos - oldpos;
3404 const_cast<GEO_PrimVDB::GridAccessor *>(this)->makeGridUnique();
3405 myGrid->setTransform(
3406 std::make_shared<Transform>(map->postTranslate(delta)));
3407 }
3408
3409 // Copy the translation from xform and set into our vertex position
3410 void
setVertexPositionAdapter(const void * xformPtr,GEO_PrimVDB & prim)3411 GEO_PrimVDB::GridAccessor::setVertexPositionAdapter(
3412 const void* xformPtr,
3413 GEO_PrimVDB &prim)
3414 {
3415 // xformPtr is assumed to point to an openvdb::vX_Y_Z::math::Transform,
3416 // for some version X.Y.Z of OpenVDB that may be newer than the one
3417 // with which libHoudiniGEO.so was built. This is safe provided that
3418 // math::Transform and its member objects are ABI-compatible between
3419 // the two OpenVDB versions.
3420 const openvdb::math::Transform& xform =
3421 *static_cast<const openvdb::math::Transform*>(xformPtr);
3422 if (myGrid && &myGrid->transform() == &xform)
3423 return;
3424 prim.incrTransformUniqueId();
3425 prim.getDetail().setPos3(
3426 prim.vertexPoint(0), UTvdbConvert(vdbTranslation(xform)));
3427 }
3428
3429 void
setTransformAdapter(const void * xformPtr,GEO_PrimVDB & prim)3430 GEO_PrimVDB::GridAccessor::setTransformAdapter(
3431 const void* xformPtr,
3432 GEO_PrimVDB &prim)
3433 {
3434 if (!myGrid)
3435 return;
3436 // xformPtr is assumed to point to an openvdb::vX_Y_Z::math::Transform,
3437 // for some version X.Y.Z of OpenVDB that may be newer than the one
3438 // with which libHoudiniGEO.so was built. This is safe provided that
3439 // math::Transform and its member objects are ABI-compatible between
3440 // the two OpenVDB versions.
3441 const openvdb::math::Transform& xform =
3442 *static_cast<const openvdb::math::Transform*>(xformPtr);
3443 setVertexPosition(xform, prim);
3444 myGrid->setTransform(xform.copy());
3445 }
3446
3447
3448 void
setGridAdapter(const void * gridPtr,GEO_PrimVDB & prim,bool copyPosition)3449 GEO_PrimVDB::GridAccessor::setGridAdapter(
3450 const void* gridPtr,
3451 GEO_PrimVDB &prim,
3452 bool copyPosition)
3453 {
3454 // gridPtr is assumed to point to an openvdb::vX_Y_Z::GridBase, for some
3455 // version X.Y.Z of OpenVDB that may be newer than the one with which
3456 // libHoudiniGEO.so was built. This is safe provided that GridBase and
3457 // its member objects are ABI-compatible between the two OpenVDB versions.
3458 const openvdb::GridBase& grid =
3459 *static_cast<const openvdb::GridBase*>(gridPtr);
3460 if (myGrid.get() == &grid)
3461 return;
3462 if (copyPosition)
3463 setVertexPosition(grid.transform(), prim);
3464 myGrid = openvdb::ConstPtrCast<openvdb::GridBase>(
3465 grid.copyGrid()); // always shallow-copy the source grid
3466 myStorageType = UTvdbGetGridType(*myGrid);
3467 }
3468
3469
3470 GEO_Primitive *
copy(int preserve_shared_pts) const3471 GEO_PrimVDB::copy(int preserve_shared_pts) const
3472 {
3473 GEO_Primitive *clone = GEO_Primitive::copy(preserve_shared_pts);
3474 if (!clone)
3475 return nullptr;
3476
3477 GEO_PrimVDB* vdb = static_cast<GEO_PrimVDB*>(clone);
3478
3479 // Give the clone the same serial number as this primitive.
3480 vdb->myUniqueId.exchange(this->getUniqueId());
3481
3482 // Give the clone a shallow copy of this primitive's grid.
3483 vdb->copyGridFrom(*this);
3484
3485 vdb->myVis = myVis;
3486
3487 return clone;
3488 }
3489
3490 void
copySubclassData(const GA_Primitive * source)3491 GEO_PrimVDB::copySubclassData(const GA_Primitive *source)
3492 {
3493 UT_ASSERT(source != this);
3494
3495 const GEO_PrimVDB* src = static_cast<const GEO_PrimVDB*>(source);
3496
3497 GEO_Primitive::copySubclassData(source);
3498
3499 // DO NOT copy P from the source, since copySubclassData
3500 // should be independent of any attributes!
3501 copyGridFrom(*src, false); // makes a shallow copy
3502
3503 myVis = src->myVis;
3504 }
3505
3506 void
assignVertex(GA_Offset new_vtx,bool update_topology)3507 GEO_PrimVDB::assignVertex(GA_Offset new_vtx, bool update_topology)
3508 {
3509 if (getVertexCount() == 1)
3510 {
3511 GA_Offset orig_vtx = getVertexOffset();
3512 if (orig_vtx == new_vtx)
3513 return;
3514 UT_ASSERT_P(GAisValid(orig_vtx));
3515 destroyVertex(orig_vtx);
3516 myVertexList.set(0, new_vtx);
3517 }
3518 else
3519 {
3520 myVertexList.setTrivial(new_vtx, 1);
3521 }
3522 if (update_topology)
3523 registerVertex(new_vtx);
3524 }
3525
3526 const char *
getGridName() const3527 GEO_PrimVDB::getGridName() const
3528 {
3529 GA_ROHandleS nameAttr(getParent(), GA_ATTRIB_PRIMITIVE, "name");
3530 return nameAttr.isValid()
3531 ? static_cast<const char *>(nameAttr.get(getMapOffset())) : "";
3532 }
3533
3534
3535 namespace // anonymous
3536 {
3537 using geo_Size = GA_Size;
3538
3539 // Intrinsic attributes
3540 enum geo_Intrinsic
3541 {
3542 geo_INTRINSIC_BACKGROUND,
3543 geo_INTRINSIC_VOXELSIZE,
3544 geo_INTRINSIC_ACTIVEVOXELDIM,
3545 geo_INTRINSIC_ACTIVEVOXELCOUNT,
3546 geo_INTRINSIC_TRANSFORM,
3547 geo_INTRINSIC_VOLUMEVISUALMODE,
3548 geo_INTRINSIC_VOLUMEVISUALDENSITY,
3549 geo_INTRINSIC_VOLUMEVISUALISO,
3550 geo_INTRINSIC_VOLUMEVISUALLOD,
3551
3552 geo_INTRINSIC_VOLUMEMINVALUE,
3553 geo_INTRINSIC_VOLUMEMAXVALUE,
3554 geo_INTRINSIC_VOLUMEAVGVALUE,
3555
3556 geo_INTRINSIC_META_GRID_CLASS,
3557 geo_INTRINSIC_META_GRID_CREATOR,
3558 geo_INTRINSIC_META_IS_LOCAL_SPACE,
3559 geo_INTRINSIC_META_SAVE_HALF_FLOAT,
3560 geo_INTRINSIC_META_VALUE_TYPE,
3561 geo_INTRINSIC_META_VECTOR_TYPE,
3562
3563 geo_NUM_INTRINSICS
3564 };
3565
3566 const UT_FSATable theMetaNames(
3567 geo_INTRINSIC_META_GRID_CLASS, "vdb_class",
3568 geo_INTRINSIC_META_GRID_CREATOR, "vdb_creator",
3569 geo_INTRINSIC_META_IS_LOCAL_SPACE, "vdb_is_local_space",
3570 geo_INTRINSIC_META_SAVE_HALF_FLOAT, "vdb_is_saved_as_half_float",
3571 geo_INTRINSIC_META_VALUE_TYPE, "vdb_value_type",
3572 geo_INTRINSIC_META_VECTOR_TYPE, "vdb_vector_type",
3573 -1, nullptr
3574 );
3575
3576 geo_Size
intrinsicBackgroundTupleSize(const GEO_PrimVDB * p)3577 intrinsicBackgroundTupleSize(const GEO_PrimVDB *p)
3578 {
3579 return UTvdbGetGridTupleSize(p->getStorageType());
3580 }
3581 template <typename GridT> void
intrinsicBackgroundV(const GridT & grid,fpreal64 * v,GA_Size n)3582 intrinsicBackgroundV(const GridT &grid, fpreal64 *v, GA_Size n)
3583 {
3584 typename GridT::ValueType background = grid.background();
3585 for (GA_Size i = 0; i < n; i++)
3586 v[i] = background[i];
3587 }
3588 template <typename GridT> void
intrinsicBackgroundS(const GridT & grid,fpreal64 * v)3589 intrinsicBackgroundS(const GridT &grid, fpreal64 *v)
3590 {
3591 v[0] = (fpreal64)grid.background();
3592 }
3593 geo_Size
intrinsicBackground(const GEO_PrimVDB * p,fpreal64 * v,GA_Size size)3594 intrinsicBackground(const GEO_PrimVDB *p, fpreal64 *v, GA_Size size)
3595 {
3596 UT_VDBType grid_type = p->getStorageType();
3597 GA_Size n = SYSmin(UTvdbGetGridTupleSize(grid_type), size);
3598
3599 UT_ASSERT(n > 0);
3600 UTvdbCallScalarType(grid_type, intrinsicBackgroundS, p->getGrid(), v)
3601 else UTvdbCallVec3Type(grid_type, intrinsicBackgroundV,
3602 p->getGrid(), v, n)
3603 else if (grid_type == UT_VDB_BOOL)
3604 {
3605 intrinsicBackgroundS<openvdb::BoolGrid>(
3606 UTvdbGridCast<openvdb::BoolGrid>(p->getGrid()), v);
3607 }
3608 else
3609 n = 0;
3610
3611 return n;
3612 }
3613
3614 geo_Size
intrinsicVoxelSize(const GEO_PrimVDB * prim,fpreal64 * v,GA_Size size)3615 intrinsicVoxelSize(const GEO_PrimVDB *prim, fpreal64 *v, GA_Size size)
3616 {
3617 openvdb::Vec3d voxel_size = prim->getGrid().voxelSize();
3618 GA_Size n = SYSmin(3, size);
3619 for (GA_Size i = 0; i < n; i++)
3620 v[i] = voxel_size[i];
3621 return n;
3622 }
3623
3624 geo_Size
intrinsicActiveVoxelDim(const GEO_PrimVDB * prim,int64 * v,GA_Size size)3625 intrinsicActiveVoxelDim(const GEO_PrimVDB *prim, int64 *v, GA_Size size)
3626 {
3627 using namespace openvdb;
3628 Coord dim = prim->getGrid().evalActiveVoxelDim();
3629 GA_Size n = SYSmin(3, size);
3630 for (GA_Size i = 0; i < n; i++)
3631 v[i] = dim[i];
3632 return n;
3633 }
3634 int64
intrinsicActiveVoxelCount(const GEO_PrimVDB * prim)3635 intrinsicActiveVoxelCount(const GEO_PrimVDB *prim)
3636 {
3637 return prim->getGrid().activeVoxelCount();
3638 }
3639
3640 geo_Size
intrinsicTransform(const GEO_PrimVDB * prim,fpreal64 * v,GA_Size size)3641 intrinsicTransform(const GEO_PrimVDB *prim, fpreal64 *v, GA_Size size)
3642 {
3643 using namespace openvdb;
3644 const GridBase & grid = prim->getGrid();
3645 const math::Transform & xform = grid.transform();
3646 math::MapBase::ConstPtr bmap = xform.baseMap();
3647 math::AffineMap::Ptr amap = bmap->getAffineMap();
3648 math::Mat4d m4 = amap->getMat4();
3649 const double * data = m4.asPointer();
3650
3651 size = SYSmin(size, 16);
3652 for (int i = 0; i < size; ++i)
3653 v[i] = data[i];
3654 return geo_Size(size);
3655 }
3656 geo_Size
intrinsicSetTransform(GEO_PrimVDB * q,const fpreal64 * v,GA_Size size)3657 intrinsicSetTransform(GEO_PrimVDB *q, const fpreal64 *v, GA_Size size)
3658 {
3659 if (size < 16)
3660 return 0;
3661 UT_DMatrix4 m(v[0], v[1], v[2], v[3],
3662 v[4], v[5], v[6], v[7],
3663 v[8], v[9], v[10], v[11],
3664 v[12], v[13], v[14], v[15]);
3665 q->setTransform4(m);
3666 return 16;
3667 }
3668 const char *
intrinsicVisualMode(const GEO_PrimVDB * p)3669 intrinsicVisualMode(const GEO_PrimVDB *p)
3670 {
3671 return GEOgetVolumeVisToken(p->getVisualization());
3672 }
3673
3674 const char *
intrinsicVisualLod(const GEO_PrimVDB * p)3675 intrinsicVisualLod(const GEO_PrimVDB *p)
3676 {
3677 return GEOgetVolumeVisLodToken(p->getVisLod());
3678 }
3679
3680 openvdb::Metadata::ConstPtr
intrinsicGetMeta(const GEO_PrimVDB * p,geo_Intrinsic id)3681 intrinsicGetMeta(const GEO_PrimVDB *p, geo_Intrinsic id)
3682 {
3683 using namespace openvdb;
3684 return p->getGrid()[theMetaNames.getToken(id) + 4];
3685 }
3686 void
intrinsicSetMeta(GEO_PrimVDB * p,geo_Intrinsic id,const openvdb::Metadata & meta)3687 intrinsicSetMeta(
3688 GEO_PrimVDB *p,
3689 geo_Intrinsic id,
3690 const openvdb::Metadata &meta)
3691 {
3692 using namespace openvdb;
3693
3694 MetaMap &meta_map = p->getMetadata();
3695 const char *name = theMetaNames.getToken(id) + 4;
3696 meta_map.removeMeta(name);
3697 meta_map.insertMeta(name, meta);
3698 }
3699
3700 void
intrinsicGetMetaString(const GEO_PrimVDB * p,geo_Intrinsic id,UT_String & v)3701 intrinsicGetMetaString(
3702 const GEO_PrimVDB *p,
3703 geo_Intrinsic id,
3704 UT_String &v)
3705 {
3706 using namespace openvdb;
3707 Metadata::ConstPtr meta = intrinsicGetMeta(p, id);
3708 if (meta)
3709 v = meta->str();
3710 else
3711 v = "";
3712 }
3713 void
intrinsicSetMetaString(GEO_PrimVDB * p,geo_Intrinsic id,const char * v)3714 intrinsicSetMetaString(
3715 GEO_PrimVDB *p,
3716 geo_Intrinsic id,
3717 const char *v)
3718 {
3719 intrinsicSetMeta(p, id, openvdb::StringMetadata(v));
3720 }
3721
3722 bool
intrinsicGetMetaBool(const GEO_PrimVDB * p,geo_Intrinsic id)3723 intrinsicGetMetaBool(const GEO_PrimVDB *p, geo_Intrinsic id)
3724 {
3725 using namespace openvdb;
3726 Metadata::ConstPtr meta = intrinsicGetMeta(p, id);
3727 if (meta)
3728 return meta->asBool();
3729 else
3730 return false;
3731 }
3732 void
intrinsicSetMetaBool(GEO_PrimVDB * p,geo_Intrinsic id,int64 v)3733 intrinsicSetMetaBool(GEO_PrimVDB *p, geo_Intrinsic id, int64 v)
3734 {
3735 intrinsicSetMeta(p, id, openvdb::BoolMetadata(v != 0));
3736 }
3737
3738 } // namespace anonymous
3739
3740 #define VDB_INTRINSIC_META_STR(CLASS, ID) { \
3741 struct callbacks { \
3742 static geo_Size evalS(const CLASS *o, UT_String &v) \
3743 { intrinsicGetMetaString(o, ID, v); return 1; } \
3744 static geo_Size evalSA(const CLASS *o, UT_StringArray &v) \
3745 { \
3746 UT_String str; \
3747 intrinsicGetMetaString(o, ID, str); \
3748 v.append(str); \
3749 return 1; \
3750 } \
3751 static geo_Size setSS(CLASS *o, const char **v, GA_Size) \
3752 { intrinsicSetMetaString(o, ID, v[0]); return 1; } \
3753 static geo_Size setSA(CLASS *o, const UT_StringArray &a) \
3754 { intrinsicSetMetaString(o, ID, a(0)); return 1; } \
3755 }; \
3756 GA_INTRINSIC_DEF_S(ID, theMetaNames.getToken(ID), 1) \
3757 myEval[ID].myS = callbacks::evalS; \
3758 myEval[ID].mySA = callbacks::evalSA; \
3759 myEval[ID].mySetSS = callbacks::setSS; \
3760 myEval[ID].mySetSA = callbacks::setSA; \
3761 myEval[ID].myReadOnly = false; \
3762 }
3763 #define VDB_INTRINSIC_META_BOOL(CLASS, ID) { \
3764 struct callbacks { \
3765 static geo_Size eval(const CLASS *o, int64 *v, GA_Size) \
3766 { v[0] = intrinsicGetMetaBool(o, ID); return 1; } \
3767 static geo_Size setFunc(CLASS *o, const int64 *v, GA_Size) \
3768 { intrinsicSetMetaBool(o, ID, v[0]); return 1; } \
3769 }; \
3770 GA_INTRINSIC_DEF_I(ID, theMetaNames.getToken(ID), 1) \
3771 myEval[ID].myI = callbacks::eval; \
3772 myEval[ID].mySetI = callbacks::setFunc; \
3773 myEval[ID].myReadOnly = false; \
3774 }
3775
3776 GA_START_INTRINSIC_DEF(GEO_PrimVDB, geo_NUM_INTRINSICS)
3777
3778 GA_INTRINSIC_VARYING_F(GEO_PrimVDB, geo_INTRINSIC_BACKGROUND,
3779 "background", intrinsicBackgroundTupleSize, intrinsicBackground);
3780 GA_INTRINSIC_TUPLE_F(GEO_PrimVDB, geo_INTRINSIC_VOXELSIZE,
3781 "voxelsize", 3, intrinsicVoxelSize);
3782
3783 GA_INTRINSIC_TUPLE_I(GEO_PrimVDB, geo_INTRINSIC_ACTIVEVOXELDIM,
3784 "activevoxeldimensions", 3, intrinsicActiveVoxelDim);
3785 GA_INTRINSIC_I(GEO_PrimVDB, geo_INTRINSIC_ACTIVEVOXELCOUNT,
3786 "activevoxelcount", intrinsicActiveVoxelCount);
3787
3788 GA_INTRINSIC_TUPLE_F(GEO_PrimVDB, geo_INTRINSIC_TRANSFORM,
3789 "transform", 16, intrinsicTransform);
3790 GA_INTRINSIC_SET_TUPLE_F(GEO_PrimVDB, geo_INTRINSIC_TRANSFORM,
3791 intrinsicSetTransform);
3792
3793 GA_INTRINSIC_S(GEO_PrimVDB, geo_INTRINSIC_VOLUMEVISUALMODE,
3794 "volumevisualmode", intrinsicVisualMode)
3795 GA_INTRINSIC_METHOD_F(GEO_PrimVDB, geo_INTRINSIC_VOLUMEVISUALDENSITY,
3796 "volumevisualdensity", getVisDensity)
3797 GA_INTRINSIC_METHOD_F(GEO_PrimVDB, geo_INTRINSIC_VOLUMEVISUALISO,
3798 "volumevisualiso", getVisIso)
3799 GA_INTRINSIC_S(GEO_PrimVDB, geo_INTRINSIC_VOLUMEVISUALLOD,
3800 "volumevisuallod", intrinsicVisualLod)
3801
3802 GA_INTRINSIC_METHOD_F(GEO_PrimVDB, geo_INTRINSIC_VOLUMEMINVALUE,
3803 "volumeminvalue", calcMinimum)
3804 GA_INTRINSIC_METHOD_F(GEO_PrimVDB, geo_INTRINSIC_VOLUMEMAXVALUE,
3805 "volumemaxvalue", calcMaximum)
3806 GA_INTRINSIC_METHOD_F(GEO_PrimVDB, geo_INTRINSIC_VOLUMEAVGVALUE,
3807 "volumeavgvalue", calcAverage)
3808
VDB_INTRINSIC_META_STR(GEO_PrimVDB,geo_INTRINSIC_META_GRID_CLASS)3809 VDB_INTRINSIC_META_STR(GEO_PrimVDB, geo_INTRINSIC_META_GRID_CLASS)
3810 VDB_INTRINSIC_META_STR(GEO_PrimVDB, geo_INTRINSIC_META_GRID_CREATOR)
3811 VDB_INTRINSIC_META_BOOL(GEO_PrimVDB, geo_INTRINSIC_META_IS_LOCAL_SPACE)
3812 VDB_INTRINSIC_META_BOOL(GEO_PrimVDB, geo_INTRINSIC_META_SAVE_HALF_FLOAT)
3813 VDB_INTRINSIC_META_STR(GEO_PrimVDB, geo_INTRINSIC_META_VALUE_TYPE)
3814 VDB_INTRINSIC_META_STR(GEO_PrimVDB, geo_INTRINSIC_META_VECTOR_TYPE)
3815
3816 GA_END_INTRINSIC_DEF(GEO_PrimVDB, GEO_Primitive)
3817
3818 /*static*/ bool
3819 GEO_PrimVDB::isIntrinsicMetadata(const char *name)
3820 {
3821 return theMetaNames.contains(name);
3822 }
3823
3824 #endif // SESI_OPENVDB || SESI_OPENVDB_PRIM
3825