Cycles: add support for volume motion blur

This adds support for rendering motion blur for volumes, using their
velocity field. This works for fluid simulations and imported VDB
volumes. For the latter, the name of the velocity field can be set per
volume object, with automatic detection of velocity fields that are
split into 3 scalar grids.

A new parameter is also added to scale velocity for more artistic control.

Like for Alembic and USD caches, a parameter to set the unit of time in
which the velocity vectors are expressed is also added. For Blender gas
simulations, the velocity unit should always be in seconds, so this is
only exposed for volume objects which may come from external OpenVDB
files.

These parameters are available under the `Render` panels for the fluid
domain and the volume object data properties respectively.

Credits: kernel advection code from Tangent Animation's Blackbird based
on earlier work by Geraldine Chua

Differential Revision: https://developer.blender.org/D14629

12 files changed, 179 insertions, 21 deletions

diff --git a/intern/cycles/scene/attribute.cpp b/intern/cycles/scene/attribute.cpp
index 0ca602362bc..df01189a54b 100644
--- a/intern/cycles/scene/attribute.cpp
+++ b/intern/cycles/scene/attribute.cpp
@@ -360,6 +360,12 @@ const char *Attribute::standard_name(AttributeStandard std)
       return "temperature";
     case ATTR_STD_VOLUME_VELOCITY:
       return "velocity";
+    case ATTR_STD_VOLUME_VELOCITY_X:
+      return "velocity_x";
+    case ATTR_STD_VOLUME_VELOCITY_Y:
+      return "velocity_y";
+    case ATTR_STD_VOLUME_VELOCITY_Z:
+      return "velocity_z";
     case ATTR_STD_POINTINESS:
       return "pointiness";
     case ATTR_STD_RANDOM_PER_ISLAND:
@@ -587,6 +593,9 @@ Attribute *AttributeSet::add(AttributeStandard std, ustring name)
       case ATTR_STD_VOLUME_FLAME:
       case ATTR_STD_VOLUME_HEAT:
       case ATTR_STD_VOLUME_TEMPERATURE:
+      case ATTR_STD_VOLUME_VELOCITY_X:
+      case ATTR_STD_VOLUME_VELOCITY_Y:
+      case ATTR_STD_VOLUME_VELOCITY_Z:
         attr = add(name, TypeDesc::TypeFloat, ATTR_ELEMENT_VOXEL);
         break;
       case ATTR_STD_VOLUME_COLOR:
diff --git a/intern/cycles/scene/camera.cpp b/intern/cycles/scene/camera.cpp
index 6aca2fcbb81..710f1c5ee90 100644
--- a/intern/cycles/scene/camera.cpp
+++ b/intern/cycles/scene/camera.cpp
@@ -397,6 +397,7 @@ void Camera::update(Scene *scene)
 
   /* motion blur */
   kcam->shuttertime = (need_motion == Scene::MOTION_BLUR) ? shuttertime : -1.0f;
+  kcam->motion_position = motion_position;
 
   /* type */
   kcam->type = camera_type;
diff --git a/intern/cycles/scene/camera.h b/intern/cycles/scene/camera.h
index 97bee430588..c150405acc2 100644
--- a/intern/cycles/scene/camera.h
+++ b/intern/cycles/scene/camera.h
@@ -30,18 +30,6 @@ class Camera : public Node {
  public:
   NODE_DECLARE
 
-  /* Specifies an offset for the shutter's time interval. */
-  enum MotionPosition {
-    /* Shutter opens at the current frame. */
-    MOTION_POSITION_START = 0,
-    /* Shutter is fully open at the current frame. */
-    MOTION_POSITION_CENTER = 1,
-    /* Shutter closes at the current frame. */
-    MOTION_POSITION_END = 2,
-
-    MOTION_NUM_POSITIONS,
-  };
-
   /* Specifies rolling shutter effect. */
   enum RollingShutterType {
     /* No rolling shutter effect. */
diff --git a/intern/cycles/scene/geometry.cpp b/intern/cycles/scene/geometry.cpp
index 351ec4f09ae..349d8ad39c7 100644
--- a/intern/cycles/scene/geometry.cpp
+++ b/intern/cycles/scene/geometry.cpp
@@ -1541,7 +1541,7 @@ void GeometryManager::device_update_preprocess(Device *device, Scene *scene, Pro
       }
 
       Volume *volume = static_cast<Volume *>(geom);
-      create_volume_mesh(volume, progress);
+      create_volume_mesh(scene, volume, progress);
 
       /* always reallocate when we have a volume, as we need to rebuild the BVH */
       device_update_flags |= DEVICE_MESH_DATA_NEEDS_REALLOC;
diff --git a/intern/cycles/scene/geometry.h b/intern/cycles/scene/geometry.h
index 0c2e70d483d..6210a64509a 100644
--- a/intern/cycles/scene/geometry.h
+++ b/intern/cycles/scene/geometry.h
@@ -216,7 +216,7 @@ class GeometryManager {
  protected:
   bool displace(Device *device, Scene *scene, Mesh *mesh, Progress &progress);
 
-  void create_volume_mesh(Volume *volume, Progress &progress);
+  void create_volume_mesh(const Scene *scene, Volume *volume, Progress &progress);
 
   /* Attributes */
   void update_osl_attributes(Device *device,
diff --git a/intern/cycles/scene/image_vdb.cpp b/intern/cycles/scene/image_vdb.cpp
index 9906606a959..8f30433c229 100644
--- a/intern/cycles/scene/image_vdb.cpp
+++ b/intern/cycles/scene/image_vdb.cpp
@@ -66,6 +66,11 @@ struct ToNanoOp {
 #  endif
 #endif
 
+VDBImageLoader::VDBImageLoader(openvdb::GridBase::ConstPtr grid_, const string &grid_name)
+    : grid_name(grid_name), grid(grid_)
+{
+}
+
 VDBImageLoader::VDBImageLoader(const string &grid_name) : grid_name(grid_name)
 {
 }
diff --git a/intern/cycles/scene/image_vdb.h b/intern/cycles/scene/image_vdb.h
index c851ef6250d..4c1aba80af8 100644
--- a/intern/cycles/scene/image_vdb.h
+++ b/intern/cycles/scene/image_vdb.h
@@ -17,6 +17,7 @@ CCL_NAMESPACE_BEGIN
 
 class VDBImageLoader : public ImageLoader {
  public:
+  VDBImageLoader(openvdb::GridBase::ConstPtr grid_, const string &grid_name);
   VDBImageLoader(const string &grid_name);
   ~VDBImageLoader();
 
diff --git a/intern/cycles/scene/object.cpp b/intern/cycles/scene/object.cpp
index 55d89fc3673..c57d26464ed 100644
--- a/intern/cycles/scene/object.cpp
+++ b/intern/cycles/scene/object.cpp
@@ -439,6 +439,14 @@ void ObjectManager::device_update_object_transform(UpdateObjectTransformState *s
       flag |= SD_OBJECT_HAS_VERTEX_MOTION;
     }
   }
+  else if (geom->is_volume()) {
+    Volume *volume = static_cast<Volume *>(geom);
+    if (volume->attributes.find(ATTR_STD_VOLUME_VELOCITY) &&
+        volume->get_velocity_scale() != 0.0f) {
+      flag |= SD_OBJECT_HAS_VOLUME_MOTION;
+      kobject.velocity_scale = volume->get_velocity_scale();
+    }
+  }
 
   if (state->need_motion == Scene::MOTION_PASS) {
     /* Clear motion array if there is no actual motion. */
@@ -488,7 +496,7 @@ void ObjectManager::device_update_object_transform(UpdateObjectTransformState *s
   kobject.dupli_generated[2] = ob->dupli_generated[2];
   kobject.numkeys = (geom->geometry_type == Geometry::HAIR) ?
                         static_cast<Hair *>(geom)->get_curve_keys().size() :
-                    (geom->geometry_type == Geometry::POINTCLOUD) ?
+                        (geom->geometry_type == Geometry::POINTCLOUD) ?
                         static_cast<PointCloud *>(geom)->num_points() :
                         0;
   kobject.dupli_uv[0] = ob->dupli_uv[0];
diff --git a/intern/cycles/scene/scene.cpp b/intern/cycles/scene/scene.cpp
index b6b53004816..b35242139ea 100644
--- a/intern/cycles/scene/scene.cpp
+++ b/intern/cycles/scene/scene.cpp
@@ -381,7 +381,7 @@ void Scene::device_update(Device *device_, Progress &progress)
   }
 }
 
-Scene::MotionType Scene::need_motion()
+Scene::MotionType Scene::need_motion() const
 {
   if (integrator->get_motion_blur())
     return MOTION_BLUR;
@@ -407,6 +407,10 @@ bool Scene::need_global_attribute(AttributeStandard std)
     return need_motion() != MOTION_NONE;
   else if (std == ATTR_STD_MOTION_VERTEX_NORMAL)
     return need_motion() == MOTION_BLUR;
+  else if (std == ATTR_STD_VOLUME_VELOCITY || std == ATTR_STD_VOLUME_VELOCITY_X ||
+           std == ATTR_STD_VOLUME_VELOCITY_Y || std == ATTR_STD_VOLUME_VELOCITY_Z) {
+    return need_motion() != MOTION_NONE;
+  }
 
   return false;
 }
diff --git a/intern/cycles/scene/scene.h b/intern/cycles/scene/scene.h
index 54c7d9d93ea..a0d2f4a6c06 100644
--- a/intern/cycles/scene/scene.h
+++ b/intern/cycles/scene/scene.h
@@ -257,7 +257,7 @@ class Scene : public NodeOwner {
   void need_global_attributes(AttributeRequestSet &attributes);
 
   enum MotionType { MOTION_NONE = 0, MOTION_PASS, MOTION_BLUR };
-  MotionType need_motion();
+  MotionType need_motion() const;
   float motion_shutter_time();
 
   bool need_update();
diff --git a/intern/cycles/scene/volume.cpp b/intern/cycles/scene/volume.cpp
index 0555fdd2fad..1357d5ed762 100644
--- a/intern/cycles/scene/volume.cpp
+++ b/intern/cycles/scene/volume.cpp
@@ -10,9 +10,9 @@
 #  include <openvdb/tools/Dense.h>
 #  include <openvdb/tools/GridTransformer.h>
 #  include <openvdb/tools/Morphology.h>
+#  include <openvdb/tools/Statistics.h>
 #endif
 
-#include "util/foreach.h"
 #include "util/hash.h"
 #include "util/log.h"
 #include "util/openvdb.h"
@@ -28,6 +28,7 @@ NODE_DEFINE(Volume)
   SOCKET_FLOAT(clipping, "Clipping", 0.001f);
   SOCKET_FLOAT(step_size, "Step Size", 0.0f);
   SOCKET_BOOLEAN(object_space, "Object Space", false);
+  SOCKET_FLOAT(velocity_scale, "Velocity Scale", 1.0f);
 
   return type;
 }
@@ -482,11 +483,141 @@ static openvdb::GridBase::ConstPtr openvdb_grid_from_device_texture(device_textu
 
   return sparse;
 }
+
+static int estimate_required_velocity_padding(openvdb::GridBase::ConstPtr grid,
+                                              float velocity_scale)
+{
+  /* TODO: we may need to also find outliers and clamp them to avoid adding too much padding. */
+  openvdb::math::Extrema extrema;
+  openvdb::Vec3d voxel_size;
+
+  /* External .vdb files have a vec3 type for velocity, but the Blender exporter creates a vec4. */
+  if (grid->isType<openvdb::Vec3fGrid>()) {
+    openvdb::Vec3fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast<openvdb::Vec3fGrid>(grid);
+    extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn());
+    voxel_size = vel_grid->voxelSize();
+  }
+  else if (grid->isType<openvdb::Vec4fGrid>()) {
+    openvdb::Vec4fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast<openvdb::Vec4fGrid>(grid);
+    extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn());
+    voxel_size = vel_grid->voxelSize();
+  }
+  else {
+    assert(0);
+    return 0;
+  }
+
+  /* We should only have uniform grids, so x = y = z, but we never know. */
+  const double max_voxel_size = openvdb::math::Max(voxel_size.x(), voxel_size.y(), voxel_size.z());
+  if (max_voxel_size == 0.0) {
+    return 0;
+  }
+
+  const double estimated_padding = extrema.max() * static_cast<double>(velocity_scale) /
+                                   max_voxel_size;
+
+  return static_cast<int>(std::ceil(estimated_padding));
+}
+
+static openvdb::FloatGrid::ConstPtr get_vdb_for_attribute(Volume *volume, AttributeStandard std)
+{
+  Attribute *attr = volume->attributes.find(std);
+  if (!attr) {
+    return nullptr;
+  }
+
+  ImageHandle &handle = attr->data_voxel();
+  VDBImageLoader *vdb_loader = handle.vdb_loader();
+  if (!vdb_loader) {
+    return nullptr;
+  }
+
+  openvdb::GridBase::ConstPtr grid = vdb_loader->get_grid();
+  if (!grid) {
+    return nullptr;
+  }
+
+  if (!grid->isType<openvdb::FloatGrid>()) {
+    return nullptr;
+  }
+
+  return openvdb::gridConstPtrCast<openvdb::FloatGrid>(grid);
+}
+
+class MergeScalarGrids {
+  typedef openvdb::FloatTree ScalarTree;
+
+  openvdb::tree::ValueAccessor<const ScalarTree> m_acc_x, m_acc_y, m_acc_z;
+
+ public:
+  MergeScalarGrids(const ScalarTree *x_tree, const ScalarTree *y_tree, const ScalarTree *z_tree)
+      : m_acc_x(*x_tree), m_acc_y(*y_tree), m_acc_z(*z_tree)
+  {
+  }
+
+  MergeScalarGrids(const MergeScalarGrids &other)
+      : m_acc_x(other.m_acc_x), m_acc_y(other.m_acc_y), m_acc_z(other.m_acc_z)
+  {
+  }
+
+  void operator()(const openvdb::Vec3STree::ValueOnIter &it) const
+  {
+    using namespace openvdb;
+
+    const math::Coord xyz = it.getCoord();
+    float x = m_acc_x.getValue(xyz);
+    float y = m_acc_y.getValue(xyz);
+    float z = m_acc_z.getValue(xyz);
+
+    it.setValue(math::Vec3s(x, y, z));
+  }
+};
+
+static void merge_scalar_grids_for_velocity(const Scene *scene, Volume *volume)
+{
+  if (volume->attributes.find(ATTR_STD_VOLUME_VELOCITY)) {
+    /* A vector grid for velocity is already available. */
+    return;
+  }
+
+  openvdb::FloatGrid::ConstPtr vel_x_grid = get_vdb_for_attribute(volume,
+                                                                  ATTR_STD_VOLUME_VELOCITY_X);
+  openvdb::FloatGrid::ConstPtr vel_y_grid = get_vdb_for_attribute(volume,
+                                                                  ATTR_STD_VOLUME_VELOCITY_Y);
+  openvdb::FloatGrid::ConstPtr vel_z_grid = get_vdb_for_attribute(volume,
+                                                                  ATTR_STD_VOLUME_VELOCITY_Z);
+
+  if (!(vel_x_grid && vel_y_grid && vel_z_grid)) {
+    return;
+  }
+
+  openvdb::Vec3fGrid::Ptr vecgrid = openvdb::Vec3SGrid::create(openvdb::Vec3s(0.0f));
+
+  /* Activate voxels in the vector grid based on the scalar grids to ensure thread safety during
+   * the merge. */
+  vecgrid->tree().topologyUnion(vel_x_grid->tree());
+  vecgrid->tree().topologyUnion(vel_y_grid->tree());
+  vecgrid->tree().topologyUnion(vel_z_grid->tree());
+
+  MergeScalarGrids op(&vel_x_grid->tree(), &vel_y_grid->tree(), &vel_z_grid->tree());
+  openvdb::tools::foreach (vecgrid->beginValueOn(), op, true, false);
+
+  /* Assume all grids have the same transformation. */
+  openvdb::math::Transform::Ptr transform = openvdb::ConstPtrCast<openvdb::math::Transform>(
+      vel_x_grid->transformPtr());
+  vecgrid->setTransform(transform);
+
+  /* Make an attribute for it. */
+  Attribute *attr = volume->attributes.add(ATTR_STD_VOLUME_VELOCITY);
+  ImageLoader *loader = new VDBImageLoader(vecgrid, "merged_velocity");
+  ImageParams params;
+  attr->data_voxel() = scene->image_manager->add_image(loader, params);
+}
 #endif
 
 /* ************************************************************************** */
 
-void GeometryManager::create_volume_mesh(Volume *volume, Progress &progress)
+void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Progress &progress)
 {
   string msg = string_printf("Computing Volume Mesh %s", volume->name.c_str());
   progress.set_status("Updating Mesh", msg);
@@ -495,7 +626,7 @@ void GeometryManager::create_volume_mesh(Volume *volume, Progress &progress)
   Shader *volume_shader = NULL;
   int pad_size = 0;
 
-  foreach (Node *node, volume->get_used_shaders()) {
+  for (Node *node : volume->get_used_shaders()) {
     Shader *shader = static_cast<Shader *>(node);
 
     if (!shader->has_volume) {
@@ -529,7 +660,9 @@ void GeometryManager::create_volume_mesh(Volume *volume, Progress &progress)
   VolumeMeshBuilder builder;
 
 #ifdef WITH_OPENVDB
-  foreach (Attribute &attr, volume->attributes.attributes) {
+  merge_scalar_grids_for_velocity(scene, volume);
+
+  for (Attribute &attr : volume->attributes.attributes) {
     if (attr.element != ATTR_ELEMENT_VOXEL) {
       continue;
     }
@@ -567,9 +700,17 @@ void GeometryManager::create_volume_mesh(Volume *volume, Progress &progress)
     }
 
     if (grid) {
+      /* Add padding based on the maximum velocity vector. */
+      if (attr.std == ATTR_STD_VOLUME_VELOCITY && scene->need_motion() != Scene::MOTION_NONE) {
+        pad_size = max(pad_size,
+                       estimate_required_velocity_padding(grid, volume->get_velocity_scale()));
+      }
+
       builder.add_grid(grid, do_clipping, volume->get_clipping());
     }
   }
+#else
+  (void)scene;
 #endif
 
   /* If nothing to build, early out. */
diff --git a/intern/cycles/scene/volume.h b/intern/cycles/scene/volume.h
index cae0f5f5bce..2b94aaf5253 100644
--- a/intern/cycles/scene/volume.h
+++ b/intern/cycles/scene/volume.h
@@ -18,6 +18,7 @@ class Volume : public Mesh {
   NODE_SOCKET_API(float, clipping)
   NODE_SOCKET_API(float, step_size)
   NODE_SOCKET_API(bool, object_space)
+  NODE_SOCKET_API(float, velocity_scale)
 
   virtual void clear(bool preserve_shaders = false) override;
 };