Geometry Nodes: Distribute Points in Volume

This commit adds a node to distribute points inside of volume grids.
The "Random" mode usese OpenVDB's "point scatter" implementation, and
there is also a "Grid" mode for uniform distributions. Both methods
operate on all of the float grids in the volume, using every voxel with
a value higher than the threshold. The random method is not stable as
the input volume deforms.

Based on a patch by Angus Stanton (@abstanton), which was based on a
patch by Kenzie (@kenziemac130).

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

9 files changed, 326 insertions, 0 deletions

diff --git a/release/scripts/startup/nodeitems_builtins.py b/release/scripts/startup/nodeitems_builtins.py
index 89b729595db..9c0635d7bd4 100644
--- a/release/scripts/startup/nodeitems_builtins.py
+++ b/release/scripts/startup/nodeitems_builtins.py
@@ -241,6 +241,7 @@ def point_node_items(context):
     space = context.space_data
     if not space:
         return
+    yield NodeItem("GeometryNodeDistributePointsInVolume")
     yield NodeItem("GeometryNodeDistributePointsOnFaces")
     yield NodeItem("GeometryNodePoints")
     yield NodeItem("GeometryNodePointsToVertices")
diff --git a/source/blender/blenkernel/BKE_node.h b/source/blender/blenkernel/BKE_node.h
index 7d019dcc7cb..35d3fbf0b73 100644
--- a/source/blender/blenkernel/BKE_node.h
+++ b/source/blender/blenkernel/BKE_node.h
@@ -1526,6 +1526,7 @@ struct TexResult;
 #define GEO_NODE_EDGE_PATHS_TO_CURVES 1169
 #define GEO_NODE_EDGE_PATHS_TO_SELECTION 1170
 #define GEO_NODE_MESH_FACE_SET_BOUNDARIES 1171
+#define GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME 1172
 
 /** \} */
 
diff --git a/source/blender/blenkernel/intern/node.cc b/source/blender/blenkernel/intern/node.cc
index 1af04d3034c..04f5b131743 100644
--- a/source/blender/blenkernel/intern/node.cc
+++ b/source/blender/blenkernel/intern/node.cc
@@ -4702,6 +4702,7 @@ static void registerGeometryNodes()
   register_node_type_geo_curve_trim();
   register_node_type_geo_deform_curves_on_surface();
   register_node_type_geo_delete_geometry();
+  register_node_type_geo_distribute_points_in_volume();
   register_node_type_geo_distribute_points_on_faces();
   register_node_type_geo_dual_mesh();
   register_node_type_geo_duplicate_elements();
diff --git a/source/blender/makesdna/DNA_node_types.h b/source/blender/makesdna/DNA_node_types.h
index a416e184d2b..7832541e360 100644
--- a/source/blender/makesdna/DNA_node_types.h
+++ b/source/blender/makesdna/DNA_node_types.h
@@ -1595,6 +1595,11 @@ typedef struct NodeGeometryUVUnwrap {
   uint8_t method;
 } NodeGeometryUVUnwrap;
 
+typedef struct NodeGeometryDistributePointsInVolume {
+  /* GeometryNodePointDistributeVolumeMode. */
+  uint8_t mode;
+} NodeGeometryDistributePointsInVolume;
+
 typedef struct NodeFunctionCompare {
   /* NodeCompareOperation */
   int8_t operation;
@@ -2176,6 +2181,11 @@ typedef enum GeometryNodeTriangulateQuads {
   GEO_NODE_TRIANGULATE_QUAD_LONGEDGE = 4,
 } GeometryNodeTriangulateQuads;
 
+typedef enum GeometryNodeDistributePointsInVolumeMode {
+  GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM = 0,
+  GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID = 1,
+} GeometryNodeDistributePointsInVolumeMode;
+
 typedef enum GeometryNodeDistributePointsOnFacesMode {
   GEO_NODE_POINT_DISTRIBUTE_POINTS_ON_FACES_RANDOM = 0,
   GEO_NODE_POINT_DISTRIBUTE_POINTS_ON_FACES_POISSON = 1,
diff --git a/source/blender/makesrna/intern/rna_nodetree.c b/source/blender/makesrna/intern/rna_nodetree.c
index daae1d6f3e4..971043158e1 100644
--- a/source/blender/makesrna/intern/rna_nodetree.c
+++ b/source/blender/makesrna/intern/rna_nodetree.c
@@ -9614,6 +9614,31 @@ static void def_geo_extrude_mesh(StructRNA *srna)
   RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");
 }
 
+static void def_geo_distribute_points_in_volume(StructRNA *srna)
+{
+  PropertyRNA *prop;
+
+  static const EnumPropertyItem mode_items[] = {
+      {GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM,
+       "DENSITY_RANDOM",
+       0,
+       "Random",
+       "Distribute points randomly inside of the volume"},
+      {GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID,
+       "DENSITY_GRID",
+       0,
+       "Grid",
+       "Distribute the points in a grid pattern inside of the volume"},
+      {0, NULL, 0, NULL, NULL},
+  };
+
+  RNA_def_struct_sdna_from(srna, "NodeGeometryDistributePointsInVolume", "storage");
+  prop = RNA_def_property(srna, "mode", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_items(prop, mode_items);
+  RNA_def_property_ui_text(prop, "Distribution Method", "Method to use for scattering points");
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_socket_update");
+}
+
 static void def_geo_distribute_points_on_faces(StructRNA *srna)
 {
   PropertyRNA *prop;
diff --git a/source/blender/nodes/NOD_geometry.h b/source/blender/nodes/NOD_geometry.h
index 63af2c71d45..a5e8ac2b093 100644
--- a/source/blender/nodes/NOD_geometry.h
+++ b/source/blender/nodes/NOD_geometry.h
@@ -49,6 +49,7 @@ void register_node_type_geo_curve_to_points(void);
 void register_node_type_geo_curve_trim(void);
 void register_node_type_geo_deform_curves_on_surface(void);
 void register_node_type_geo_delete_geometry(void);
+void register_node_type_geo_distribute_points_in_volume(void);
 void register_node_type_geo_distribute_points_on_faces(void);
 void register_node_type_geo_dual_mesh(void);
 void register_node_type_geo_duplicate_elements(void);
diff --git a/source/blender/nodes/NOD_static_types.h b/source/blender/nodes/NOD_static_types.h
index 73fb2bf32c8..10eb1721133 100644
--- a/source/blender/nodes/NOD_static_types.h
+++ b/source/blender/nodes/NOD_static_types.h
@@ -305,6 +305,7 @@ DefNode(GeometryNode, GEO_NODE_CURVE_TO_MESH, 0, "CURVE_TO_MESH", CurveToMesh, "
 DefNode(GeometryNode, GEO_NODE_CURVE_TO_POINTS, def_geo_curve_to_points, "CURVE_TO_POINTS", CurveToPoints, "Curve to Points", "Generate a point cloud by sampling positions along curves")
 DefNode(GeometryNode, GEO_NODE_DEFORM_CURVES_ON_SURFACE, 0, "DEFORM_CURVES_ON_SURFACE", DeformCurvesOnSurface, "Deform Curves on Surface", "Translate and rotate curves based on changes between the object's original and evaluated surface mesh")
 DefNode(GeometryNode, GEO_NODE_DELETE_GEOMETRY, def_geo_delete_geometry, "DELETE_GEOMETRY", DeleteGeometry, "Delete Geometry", "Remove selected elements of a geometry")
+DefNode(GeometryNode, GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME, def_geo_distribute_points_in_volume, "DISTRIBUTE_POINTS_IN_VOLUME", DistributePointsInVolume, "Distribute Points In Volume", "Generate points inside a volume")
 DefNode(GeometryNode, GEO_NODE_DISTRIBUTE_POINTS_ON_FACES, def_geo_distribute_points_on_faces, "DISTRIBUTE_POINTS_ON_FACES", DistributePointsOnFaces, "Distribute Points on Faces", "Generate points spread out on the surface of a mesh")
 DefNode(GeometryNode, GEO_NODE_DUAL_MESH, 0, "DUAL_MESH", DualMesh, "Dual Mesh", "Convert Faces into vertices and vertices into faces")
 DefNode(GeometryNode, GEO_NODE_DUPLICATE_ELEMENTS, def_geo_duplicate_elements, "DUPLICATE_ELEMENTS", DuplicateElements, "Duplicate Elements", "Generate an arbitrary number copies of each selected input element")
diff --git a/source/blender/nodes/geometry/CMakeLists.txt b/source/blender/nodes/geometry/CMakeLists.txt
index 8e71bc8bbb8..ce285488109 100644
--- a/source/blender/nodes/geometry/CMakeLists.txt
+++ b/source/blender/nodes/geometry/CMakeLists.txt
@@ -59,6 +59,7 @@ set(SRC
   nodes/node_geo_curve_trim.cc
   nodes/node_geo_deform_curves_on_surface.cc
   nodes/node_geo_delete_geometry.cc
+  nodes/node_geo_distribute_points_in_volume.cc
   nodes/node_geo_distribute_points_on_faces.cc
   nodes/node_geo_dual_mesh.cc
   nodes/node_geo_duplicate_elements.cc
diff --git a/source/blender/nodes/geometry/nodes/node_geo_distribute_points_in_volume.cc b/source/blender/nodes/geometry/nodes/node_geo_distribute_points_in_volume.cc
new file mode 100644
index 00000000000..24f81a81b3e
--- /dev/null
+++ b/source/blender/nodes/geometry/nodes/node_geo_distribute_points_in_volume.cc
@@ -0,0 +1,285 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#ifdef WITH_OPENVDB
+#  include <openvdb/openvdb.h>
+#  include <openvdb/tools/Interpolation.h>
+#  include <openvdb/tools/PointScatter.h>
+#endif
+
+#include "DNA_node_types.h"
+#include "DNA_pointcloud_types.h"
+
+#include "BKE_pointcloud.h"
+#include "BKE_volume.h"
+
+#include "UI_interface.h"
+#include "UI_resources.h"
+
+#include "node_geometry_util.hh"
+
+namespace blender::nodes {
+
+NODE_STORAGE_FUNCS(NodeGeometryDistributePointsInVolume)
+
+static void geo_node_distribute_points_in_volume_declare(NodeDeclarationBuilder &b)
+{
+  b.add_input<decl::Geometry>(N_("Volume")).supported_type(GEO_COMPONENT_TYPE_VOLUME);
+  b.add_input<decl::Float>(N_("Density"))
+      .default_value(1.0f)
+      .min(0.0f)
+      .max(100000.0f)
+      .subtype(PROP_NONE)
+      .description(N_("Number of points to sample per unit volume"));
+  b.add_input<decl::Int>(N_("Seed"))
+      .min(-10000)
+      .max(10000)
+      .description(N_("Seed used by the random number generator to generate random points"));
+  b.add_input<decl::Vector>(N_("Spacing"))
+      .default_value({0.3, 0.3, 0.3})
+      .min(0.0001f)
+      .subtype(PROP_XYZ)
+      .description(N_("Spacing between grid points"));
+  b.add_input<decl::Float>(N_("Threshold"))
+      .default_value(0.1f)
+      .min(0.0f)
+      .max(FLT_MAX)
+      .description(N_("Minimum density of a volume cell to contain a grid point"));
+  b.add_output<decl::Geometry>(N_("Points"));
+}
+
+static void geo_node_distribute_points_in_volume_layout(uiLayout *layout,
+                                                        bContext *UNUSED(C),
+                                                        PointerRNA *ptr)
+{
+  uiItemR(layout, ptr, "mode", 0, "", ICON_NONE);
+}
+
+static void node_distribute_points_in_volume_init(bNodeTree *UNUSED(ntree), bNode *node)
+{
+  NodeGeometryDistributePointsInVolume *data = MEM_cnew<NodeGeometryDistributePointsInVolume>(
+      __func__);
+  data->mode = GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM;
+  node->storage = data;
+}
+
+static void node_distribute_points_in_volume_update(bNodeTree *ntree, bNode *node)
+{
+  const NodeGeometryDistributePointsInVolume &storage = node_storage(*node);
+  GeometryNodeDistributePointsInVolumeMode mode =
+      static_cast<GeometryNodeDistributePointsInVolumeMode>(storage.mode);
+
+  bNodeSocket *sock_density = ((bNodeSocket *)(node->inputs.first))->next;
+  bNodeSocket *sock_seed = sock_density->next;
+  bNodeSocket *sock_spacing = sock_seed->next;
+  bNodeSocket *sock_threshold = sock_spacing->next;
+
+  nodeSetSocketAvailability(
+      ntree, sock_density, mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM);
+  nodeSetSocketAvailability(
+      ntree, sock_seed, mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM);
+  nodeSetSocketAvailability(
+      ntree, sock_spacing, mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID);
+  nodeSetSocketAvailability(
+      ntree, sock_threshold, mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID);
+}
+
+#ifdef WITH_OPENVDB
+/* Implements the interface required by #openvdb::tools::NonUniformPointScatter. */
+class PositionsVDBWrapper {
+ private:
+  float3 offset_fix_;
+  Vector<float3> &vector_;
+
+ public:
+  PositionsVDBWrapper(Vector<float3> &vector, const float3 offset_fix)
+      : offset_fix_(offset_fix), vector_(vector)
+  {
+  }
+  PositionsVDBWrapper(const PositionsVDBWrapper &wrapper) = default;
+
+  void add(const openvdb::Vec3R &pos)
+  {
+    vector_.append((float3((float)pos[0], (float)pos[1], (float)pos[2]) + offset_fix_));
+  }
+};
+
+/* Use #std::mt19937 as a random number generator,
+ * it has a very long period and thus there should be no visible patterns in the generated points.
+ */
+using RNGType = std::mt19937;
+/* Non-uniform scatter allows the amount of points to be scaled with the volume's density. */
+using NonUniformPointScatterVDB =
+    openvdb::tools::NonUniformPointScatter<PositionsVDBWrapper, RNGType>;
+
+static void point_scatter_density_random(const openvdb::FloatGrid &grid,
+                                         const float density,
+                                         const int seed,
+                                         Vector<float3> &r_positions)
+{
+  /* Offset points by half a voxel so that grid points are aligned with world grid points. */
+  const float3 offset_fix = {0.5f * (float)grid.voxelSize().x(),
+                             0.5f * (float)grid.voxelSize().y(),
+                             0.5f * (float)grid.voxelSize().z()};
+  /* Setup and call into OpenVDB's point scatter API. */
+  PositionsVDBWrapper vdb_position_wrapper = PositionsVDBWrapper(r_positions, offset_fix);
+  RNGType random_generator(seed);
+  NonUniformPointScatterVDB point_scatter(vdb_position_wrapper, density, random_generator);
+  point_scatter(grid);
+}
+
+static void point_scatter_density_grid(const openvdb::FloatGrid &grid,
+                                       const float3 spacing,
+                                       const float threshold,
+                                       Vector<float3> &r_positions)
+{
+  const openvdb::Vec3d half_voxel(0.5, 0.5, 0.5);
+  const openvdb::Vec3d voxel_spacing((double)spacing.x / grid.voxelSize().x(),
+                                     (double)spacing.y / grid.voxelSize().y(),
+                                     (double)spacing.z / grid.voxelSize().z());
+
+  /* Abort if spacing is zero. */
+  const double min_spacing = std::min(voxel_spacing.x(),
+                                      std::min(voxel_spacing.y(), voxel_spacing.z()));
+  if (std::abs(min_spacing) < 0.0001) {
+    return;
+  }
+
+  /* Iterate through tiles and voxels on the grid. */
+  for (openvdb::FloatGrid::ValueOnCIter cell = grid.cbeginValueOn(); cell; ++cell) {
+    /* Check if the cell's value meets the minimum threshold. */
+    if (cell.getValue() < threshold) {
+      continue;
+    }
+    /* Compute the bounding box of each tile/voxel. */
+    const openvdb::CoordBBox bbox = cell.getBoundingBox();
+    const openvdb::Vec3d box_min = bbox.min().asVec3d() - half_voxel;
+    const openvdb::Vec3d box_max = bbox.max().asVec3d() + half_voxel;
+
+    /* Pick a starting point rounded up to the nearest possible point. */
+    double abs_spacing_x = std::abs(voxel_spacing.x());
+    double abs_spacing_y = std::abs(voxel_spacing.y());
+    double abs_spacing_z = std::abs(voxel_spacing.z());
+    const openvdb::Vec3d start(ceil(box_min.x() / abs_spacing_x) * abs_spacing_x,
+                               ceil(box_min.y() / abs_spacing_y) * abs_spacing_y,
+                               ceil(box_min.z() / abs_spacing_z) * abs_spacing_z);
+
+    /* Iterate through all possible points in box. */
+    for (double x = start.x(); x < box_max.x(); x += abs_spacing_x) {
+      for (double y = start.y(); y < box_max.y(); y += abs_spacing_y) {
+        for (double z = start.z(); z < box_max.z(); z += abs_spacing_z) {
+          /* Transform with grid matrix and add point. */
+          const openvdb::Vec3d idx_pos(x, y, z);
+          const openvdb::Vec3d local_pos = grid.indexToWorld(idx_pos + half_voxel);
+          r_positions.append({(float)local_pos.x(), (float)local_pos.y(), (float)local_pos.z()});
+        }
+      }
+    }
+  }
+}
+
+#endif /* WITH_OPENVDB */
+
+static void geo_node_distribute_points_in_volume_exec(GeoNodeExecParams params)
+{
+  GeometrySet geometry_set_in = params.extract_input<GeometrySet>("Volume");
+
+  const NodeGeometryDistributePointsInVolume &storage = node_storage(params.node());
+  const GeometryNodeDistributePointsInVolumeMode mode =
+      static_cast<GeometryNodeDistributePointsInVolumeMode>(storage.mode);
+
+#ifdef WITH_OPENVDB
+
+  float density;
+  int seed;
+  float3 spacing{0, 0, 0};
+  float threshold;
+
+  if (mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM) {
+    density = params.extract_input<float>("Density");
+    seed = params.extract_input<int>("Seed");
+  }
+  else if (mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID) {
+    spacing = params.extract_input<float3>("Spacing");
+    threshold = params.extract_input<float>("Threshold");
+  }
+
+  geometry_set_in.modify_geometry_sets([&](GeometrySet &geometry_set) {
+    if (!geometry_set.has_volume()) {
+      geometry_set.keep_only({GEO_COMPONENT_TYPE_POINT_CLOUD, GEO_COMPONENT_TYPE_INSTANCES});
+      return;
+    }
+    const VolumeComponent *component = geometry_set.get_component_for_read<VolumeComponent>();
+    const Volume *volume = component->get_for_read();
+
+    Vector<float3> positions;
+
+    for (const int i : IndexRange(BKE_volume_num_grids(volume))) {
+      const VolumeGrid *volume_grid = BKE_volume_grid_get_for_read(volume, i);
+      if (volume_grid == nullptr) {
+        continue;
+      }
+
+      openvdb::GridBase::ConstPtr base_grid = BKE_volume_grid_openvdb_for_read(volume,
+                                                                               volume_grid);
+      if (!base_grid) {
+        continue;
+      }
+
+      if (!base_grid->isType<openvdb::FloatGrid>()) {
+        continue;
+      }
+
+      const openvdb::FloatGrid::ConstPtr grid = openvdb::gridConstPtrCast<openvdb::FloatGrid>(
+          base_grid);
+
+      if (mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM) {
+        point_scatter_density_random(*grid, density, seed, positions);
+      }
+      else if (mode == GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID) {
+        point_scatter_density_grid(*grid, spacing, threshold, positions);
+      }
+    }
+
+    PointCloud *pointcloud = BKE_pointcloud_new_nomain(positions.size());
+    bke::MutableAttributeAccessor point_attributes = pointcloud->attributes_for_write();
+    bke::SpanAttributeWriter<float3> point_positions =
+        point_attributes.lookup_or_add_for_write_only_span<float3>("position", ATTR_DOMAIN_POINT);
+    bke::SpanAttributeWriter<float> point_radii =
+        point_attributes.lookup_or_add_for_write_only_span<float>("radius", ATTR_DOMAIN_POINT);
+
+    point_positions.span.copy_from(positions);
+    point_radii.span.fill(0.05f);
+    point_positions.finish();
+    point_radii.finish();
+
+    geometry_set.replace_pointcloud(pointcloud);
+    geometry_set.keep_only({GEO_COMPONENT_TYPE_POINT_CLOUD, GEO_COMPONENT_TYPE_INSTANCES});
+  });
+
+  params.set_output("Points", std::move(geometry_set_in));
+
+#else  /* WITH_OPENVDB */
+  params.error_message_add(NodeWarningType::Error, TIP_("Blender is compiled without OpenVDB"));
+#endif /* !WITH_OPENVDB */
+}
+}  // namespace blender::nodes
+
+void register_node_type_geo_distribute_points_in_volume()
+{
+  static bNodeType ntype;
+  geo_node_type_base(&ntype,
+                     GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME,
+                     "Distribute Points in Volume",
+                     NODE_CLASS_GEOMETRY);
+  node_type_storage(&ntype,
+                    "NodeGeometryDistributePointsInVolume",
+                    node_free_standard_storage,
+                    node_copy_standard_storage);
+  node_type_init(&ntype, blender::nodes::node_distribute_points_in_volume_init);
+  node_type_update(&ntype, blender::nodes::node_distribute_points_in_volume_update);
+  node_type_size(&ntype, 170, 100, 320);
+  ntype.declare = blender::nodes::geo_node_distribute_points_in_volume_declare;
+  ntype.geometry_node_execute = blender::nodes::geo_node_distribute_points_in_volume_exec;
+  ntype.draw_buttons = blender::nodes::geo_node_distribute_points_in_volume_layout;
+  nodeRegisterType(&ntype);
+}