Geometry Nodes: Poisson disk point distribution node/method

This patch does two things:
* Introduce a Seed to the random distribution method
* Bring in a new distribution method for the point scattering node

Patch Review: https://developer.blender.org/D9787

Note: This commit doesn't not handle doversion. Which means that users
need to manually update their files that were using the Point Distribute
node and reconnect inputs to the "Maximum Density" socket.

Original patch by Sebastian Parborg, with changes to not rely on the cy
libraries and overall cleanup.

Patch review by Jacques Lucke, besides help with the new "heap" system
that was required for this algorithm.

Based on Cem Yuksel. 2015. Sample Elimination for Generating Poisson Disk
Sample. Sets. Computer Graphics Forum 34, 2 (May 2015), 25-32
http://www.cemyuksel.com/research/sampleelimination/

10 files changed, 521 insertions, 9 deletions

diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h
index c34b71a60f9..5e0ea4f2a99 100644
--- a/source/blender/blenlib/BLI_kdopbvh.h
+++ b/source/blender/blenlib/BLI_kdopbvh.h
@@ -178,6 +178,7 @@ int *BLI_bvhtree_intersect_plane(BVHTree *tree, float plane[4], uint *r_intersec
 int BLI_bvhtree_get_len(const BVHTree *tree);
 int BLI_bvhtree_get_tree_type(const BVHTree *tree);
 float BLI_bvhtree_get_epsilon(const BVHTree *tree);
+void BLI_bvhtree_get_bounding_box(BVHTree *tree, float r_bb_min[3], float r_bb_max[3]);
 
 /* find nearest node to the given coordinates
  * (if nearest is given it will only search nodes where
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
index f126c5a977b..0f90ad3a490 100644
--- a/source/blender/blenlib/intern/BLI_kdopbvh.c
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -1076,6 +1076,25 @@ float BLI_bvhtree_get_epsilon(const BVHTree *tree)
   return tree->epsilon;
 }
 
+/**
+ * This function returns the bounding box of the BVH tree.
+ */
+void BLI_bvhtree_get_bounding_box(BVHTree *tree, float r_bb_min[3], float r_bb_max[3])
+{
+  BVHNode *root = tree->nodes[tree->totleaf];
+  if (root != NULL) {
+    const float bb_min[3] = {root->bv[0], root->bv[2], root->bv[4]};
+    const float bb_max[3] = {root->bv[1], root->bv[3], root->bv[5]};
+    copy_v3_v3(r_bb_min, bb_min);
+    copy_v3_v3(r_bb_max, bb_max);
+  }
+  else {
+    BLI_assert(false);
+    zero_v3(r_bb_min);
+    zero_v3(r_bb_max);
+  }
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */
diff --git a/source/blender/editors/space_node/drawnode.c b/source/blender/editors/space_node/drawnode.c
index 421d645d7bd..2857c08cad6 100644
--- a/source/blender/editors/space_node/drawnode.c
+++ b/source/blender/editors/space_node/drawnode.c
@@ -3208,6 +3208,13 @@ static void node_geometry_buts_attribute_mix(uiLayout *layout,
   uiItemR(col, ptr, "input_type_b", DEFAULT_FLAGS, IFACE_("B"), ICON_NONE);
 }
 
+static void node_geometry_buts_attribute_point_distribute(uiLayout *layout,
+                                                          bContext *UNUSED(C),
+                                                          PointerRNA *ptr)
+{
+  uiItemR(layout, ptr, "distribute_method", DEFAULT_FLAGS, "", ICON_NONE);
+}
+
 static void node_geometry_set_butfunc(bNodeType *ntype)
 {
   switch (ntype->type) {
@@ -3235,6 +3242,9 @@ static void node_geometry_set_butfunc(bNodeType *ntype)
     case GEO_NODE_ATTRIBUTE_MIX:
       ntype->draw_buttons = node_geometry_buts_attribute_mix;
       break;
+    case GEO_NODE_POINT_DISTRIBUTE:
+      ntype->draw_buttons = node_geometry_buts_attribute_point_distribute;
+      break;
   }
 }
 
diff --git a/source/blender/makesdna/DNA_node_types.h b/source/blender/makesdna/DNA_node_types.h
index 9cf64743843..34a5372892c 100644
--- a/source/blender/makesdna/DNA_node_types.h
+++ b/source/blender/makesdna/DNA_node_types.h
@@ -1501,6 +1501,11 @@ typedef enum GeometryNodeAttributeInputMode {
   GEO_NODE_ATTRIBUTE_INPUT_COLOR = 3,
 } GeometryNodeAttributeInputMode;
 
+typedef enum GeometryNodePointDistributeMethod {
+  GEO_NODE_POINT_DISTRIBUTE_RANDOM = 0,
+  GEO_NODE_POINT_DISTRIBUTE_POISSON = 1,
+} GeometryNodePointDistributeMethod;
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/source/blender/makesrna/intern/rna_nodetree.c b/source/blender/makesrna/intern/rna_nodetree.c
index 34f4d287d6e..1731a42583b 100644
--- a/source/blender/makesrna/intern/rna_nodetree.c
+++ b/source/blender/makesrna/intern/rna_nodetree.c
@@ -452,6 +452,20 @@ static const EnumPropertyItem rna_node_geometry_attribute_input_type_items[] = {
     {0, NULL, 0, NULL, NULL},
 };
 
+static const EnumPropertyItem rna_node_geometry_point_distribute_method_items[] = {
+    {GEO_NODE_POINT_DISTRIBUTE_RANDOM,
+     "RANDOM",
+     0,
+     "Random",
+     "Distribute points randomly on the surface"},
+    {GEO_NODE_POINT_DISTRIBUTE_POISSON,
+     "POISSON",
+     0,
+     "Poisson Disk",
+     "Project points on the surface evenly with a Poisson disk distribution"},
+    {0, NULL, 0, NULL, NULL},
+};
+
 #endif
 
 #ifdef RNA_RUNTIME
@@ -8481,6 +8495,18 @@ static void def_geo_attribute_mix(StructRNA *srna)
   RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_socket_update");
 }
 
+static void def_geo_point_distribute(StructRNA *srna)
+{
+  PropertyRNA *prop;
+
+  prop = RNA_def_property(srna, "distribute_method", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_sdna(prop, NULL, "custom1");
+  RNA_def_property_enum_items(prop, rna_node_geometry_point_distribute_method_items);
+  RNA_def_property_enum_default(prop, GEO_NODE_POINT_DISTRIBUTE_RANDOM);
+  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 rna_def_shader_node(BlenderRNA *brna)
diff --git a/source/blender/nodes/CMakeLists.txt b/source/blender/nodes/CMakeLists.txt
index 779680a5da7..46536903d55 100644
--- a/source/blender/nodes/CMakeLists.txt
+++ b/source/blender/nodes/CMakeLists.txt
@@ -148,6 +148,7 @@ set(SRC
   geometry/nodes/node_geo_attribute_mix.cc
   geometry/nodes/node_geo_object_info.cc
   geometry/nodes/node_geo_point_distribute.cc
+  geometry/nodes/node_geo_point_distribute_poisson_disk.cc
   geometry/nodes/node_geo_point_instance.cc
   geometry/nodes/node_geo_subdivision_surface.cc
   geometry/nodes/node_geo_transform.cc
diff --git a/source/blender/nodes/NOD_static_types.h b/source/blender/nodes/NOD_static_types.h
index 681286ae5ab..68efb0806fd 100644
--- a/source/blender/nodes/NOD_static_types.h
+++ b/source/blender/nodes/NOD_static_types.h
@@ -271,7 +271,7 @@ DefNode(GeometryNode, GEO_NODE_EDGE_SPLIT, 0, "EDGE_SPLIT", EdgeSplit, "Edge Spl
 DefNode(GeometryNode, GEO_NODE_TRANSFORM, 0, "TRANSFORM", Transform, "Transform", "")
 DefNode(GeometryNode, GEO_NODE_SUBDIVISION_SURFACE, 0, "SUBDIVISION_SURFACE", SubdivisionSurface, "Subdivision Surface", "")
 DefNode(GeometryNode, GEO_NODE_BOOLEAN, def_geo_boolean, "BOOLEAN", Boolean, "Boolean", "")
-DefNode(GeometryNode, GEO_NODE_POINT_DISTRIBUTE, 0, "POINT_DISTRIBUTE", PointDistribute, "Point Distribute", "")
+DefNode(GeometryNode, GEO_NODE_POINT_DISTRIBUTE, def_geo_point_distribute, "POINT_DISTRIBUTE", PointDistribute, "Point Distribute", "")
 DefNode(GeometryNode, GEO_NODE_POINT_INSTANCE, def_geo_point_instance, "POINT_INSTANCE", PointInstance, "Point Instance", "")
 DefNode(GeometryNode, GEO_NODE_OBJECT_INFO, 0, "OBJECT_INFO", ObjectInfo, "Object Info", "")
 DefNode(GeometryNode, GEO_NODE_ATTRIBUTE_RANDOMIZE, def_geo_attribute_randomize, "ATTRIBUTE_RANDOMIZE", AttributeRandomize, "Attribute Randomize", "")
diff --git a/source/blender/nodes/geometry/node_geometry_util.hh b/source/blender/nodes/geometry/node_geometry_util.hh
index ec389961615..c97463cdc22 100644
--- a/source/blender/nodes/geometry/node_geometry_util.hh
+++ b/source/blender/nodes/geometry/node_geometry_util.hh
@@ -42,4 +42,10 @@ namespace blender::nodes {
 void update_attribute_input_socket_availabilities(bNode &node,
                                                   const StringRef name,
                                                   const GeometryNodeAttributeInputMode mode);
-}
+
+void poisson_disk_point_elimination(Vector<float3> const *input_points,
+                                    Vector<float3> *output_points,
+                                    float maximum_distance,
+                                    float3 boundbox);
+
+}  // namespace blender::nodes
diff --git a/source/blender/nodes/geometry/nodes/node_geo_point_distribute.cc b/source/blender/nodes/geometry/nodes/node_geo_point_distribute.cc
index 2f5f7e264bc..8be9636e14d 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_point_distribute.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_point_distribute.cc
@@ -24,6 +24,7 @@
 #include "DNA_meshdata_types.h"
 #include "DNA_pointcloud_types.h"
 
+#include "BKE_bvhutils.h"
 #include "BKE_deform.h"
 #include "BKE_mesh.h"
 #include "BKE_mesh_runtime.h"
@@ -33,8 +34,10 @@
 
 static bNodeSocketTemplate geo_node_point_distribute_in[] = {
     {SOCK_GEOMETRY, N_("Geometry")},
-    {SOCK_FLOAT, N_("Density"), 10.0f, 0.0f, 0.0f, 0.0f, 0.0f, 100000.0f, PROP_NONE},
+    {SOCK_FLOAT, N_("Minimum Distance"), 0.1f, 0.0f, 0.0f, 0.0f, 0.0f, 100000.0f, PROP_NONE},
+    {SOCK_FLOAT, N_("Maximum Density"), 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 100000.0f, PROP_NONE},
     {SOCK_STRING, N_("Density Attribute")},
+    {SOCK_INT, N_("Seed"), 0, 0, 0, 0, -10000, 10000},
     {-1, ""},
 };
 
@@ -43,11 +46,19 @@ static bNodeSocketTemplate geo_node_point_distribute_out[] = {
     {-1, ""},
 };
 
+static void node_point_distribute_update(bNodeTree *UNUSED(ntree), bNode *node)
+{
+  bNodeSocket *sock_min_dist = (bNodeSocket *)BLI_findlink(&node->inputs, 1);
+
+  nodeSetSocketAvailability(sock_min_dist, ELEM(node->custom1, GEO_NODE_POINT_DISTRIBUTE_POISSON));
+}
+
 namespace blender::nodes {
 
-static Vector<float3> scatter_points_from_mesh(const Mesh *mesh,
-                                               const float density,
-                                               const FloatReadAttribute &density_factors)
+static Vector<float3> random_scatter_points_from_mesh(const Mesh *mesh,
+                                                      const float density,
+                                                      const FloatReadAttribute &density_factors,
+                                                      const int seed)
 {
   /* This only updates a cache and can be considered to be logically const. */
   const MLoopTri *looptris = BKE_mesh_runtime_looptri_ensure(const_cast<Mesh *>(mesh));
@@ -71,7 +82,7 @@ static Vector<float3> scatter_points_from_mesh(const Mesh *mesh,
                                          3.0f;
     const float area = area_tri_v3(v0_pos, v1_pos, v2_pos);
 
-    const int looptri_seed = BLI_hash_int(looptri_index);
+    const int looptri_seed = BLI_hash_int(looptri_index + seed);
     RandomNumberGenerator looptri_rng(looptri_seed);
 
     const float points_amount_fl = area * density * looptri_density_factor;
@@ -90,17 +101,158 @@ static Vector<float3> scatter_points_from_mesh(const Mesh *mesh,
   return points;
 }
 
+struct RayCastAll_Data {
+  void *bvhdata;
+
+  BVHTree_RayCastCallback raycast_callback;
+
+  const Mesh *mesh;
+  float base_weight;
+  FloatReadAttribute *density_factors;
+  Vector<float3> *projected_points;
+  float cur_point_weight;
+};
+
+static void project_2d_bvh_callback(void *userdata,
+                                    int index,
+                                    const BVHTreeRay *ray,
+                                    BVHTreeRayHit *hit)
+{
+  struct RayCastAll_Data *data = (RayCastAll_Data *)userdata;
+  data->raycast_callback(data->bvhdata, index, ray, hit);
+  if (hit->index != -1) {
+    /* This only updates a cache and can be considered to be logically const. */
+    const MLoopTri *looptris = BKE_mesh_runtime_looptri_ensure(const_cast<Mesh *>(data->mesh));
+    const MVert *mvert = data->mesh->mvert;
+
+    const MLoopTri &looptri = looptris[index];
+    const FloatReadAttribute &density_factors = data->density_factors[0];
+
+    const int v0_index = data->mesh->mloop[looptri.tri[0]].v;
+    const int v1_index = data->mesh->mloop[looptri.tri[1]].v;
+    const int v2_index = data->mesh->mloop[looptri.tri[2]].v;
+
+    const float v0_density_factor = std::max(0.0f, density_factors[v0_index]);
+    const float v1_density_factor = std::max(0.0f, density_factors[v1_index]);
+    const float v2_density_factor = std::max(0.0f, density_factors[v2_index]);
+
+    /* Calculate barycentric weights for hit point. */
+    float3 weights;
+    interp_weights_tri_v3(
+        weights, mvert[v0_index].co, mvert[v1_index].co, mvert[v2_index].co, hit->co);
+
+    float point_weight = weights[0] * v0_density_factor + weights[1] * v1_density_factor +
+                         weights[2] * v2_density_factor;
+
+    point_weight *= data->base_weight;
+
+    if (point_weight >= FLT_EPSILON && data->cur_point_weight <= point_weight) {
+      data->projected_points->append(hit->co);
+    }
+  }
+}
+
+static Vector<float3> poisson_scatter_points_from_mesh(const Mesh *mesh,
+                                                       const float density,
+                                                       const float minimum_distance,
+                                                       const FloatReadAttribute &density_factors,
+                                                       const int seed)
+{
+  Vector<float3> points;
+
+  if (minimum_distance <= FLT_EPSILON || density <= FLT_EPSILON) {
+    return points;
+  }
+
+  /* Scatter points randomly on the mesh with higher density (5-7) times higher than desired for
+   * good quality possion disk distributions. */
+  int quality = 5;
+  const int output_points_target = 1000;
+  points.resize(output_points_target * quality);
+
+  const float required_area = output_points_target *
+                              (2.0f * sqrtf(3.0f) * minimum_distance * minimum_distance);
+  const float point_scale_multiplier = sqrtf(required_area);
+
+  {
+    const int rnd_seed = BLI_hash_int(seed);
+    RandomNumberGenerator point_rng(rnd_seed);
+
+    for (int i = 0; i < points.size(); i++) {
+      points[i].x = point_rng.get_float() * point_scale_multiplier;
+      points[i].y = point_rng.get_float() * point_scale_multiplier;
+      points[i].z = 0.0f;
+    }
+  }
+
+  /* Eliminate the scattered points until we get a possion distribution. */
+  Vector<float3> output_points(output_points_target);
+
+  const float3 bounds_max = float3(point_scale_multiplier, point_scale_multiplier, 0);
+  poisson_disk_point_elimination(&points, &output_points, 2.0f * minimum_distance, bounds_max);
+  Vector<float3> final_points;
+  final_points.reserve(output_points_target);
+
+  /* Check if we have any points we should remove from the final possion distribition. */
+  BVHTreeFromMesh treedata;
+  BKE_bvhtree_from_mesh_get(&treedata, const_cast<Mesh *>(mesh), BVHTREE_FROM_LOOPTRI, 2);
+
+  float3 bb_min, bb_max;
+  BLI_bvhtree_get_bounding_box(treedata.tree, bb_min, bb_max);
+
+  struct RayCastAll_Data data;
+  data.bvhdata = &treedata;
+  data.raycast_callback = treedata.raycast_callback;
+  data.mesh = mesh;
+  data.projected_points = &final_points;
+  data.density_factors = const_cast<FloatReadAttribute *>(&density_factors);
+  data.base_weight = std::min(
+      1.0f, density / (output_points.size() / (point_scale_multiplier * point_scale_multiplier)));
+
+  const float max_dist = bb_max[2] - bb_min[2] + 2.0f;
+  const float3 dir = float3(0, 0, -1);
+  float3 raystart;
+  raystart.z = bb_max[2] + 1.0f;
+
+  float tile_start_x_coord = bb_min[0];
+  int tile_repeat_x = ceilf((bb_max[0] - bb_min[0]) / point_scale_multiplier);
+
+  float tile_start_y_coord = bb_min[1];
+  int tile_repeat_y = ceilf((bb_max[1] - bb_min[1]) / point_scale_multiplier);
+
+  for (int x = 0; x < tile_repeat_x; x++) {
+    float tile_curr_x_coord = x * point_scale_multiplier + tile_start_x_coord;
+    for (int y = 0; y < tile_repeat_y; y++) {
+      float tile_curr_y_coord = y * point_scale_multiplier + tile_start_y_coord;
+      for (int idx = 0; idx < output_points.size(); idx++) {
+        raystart.x = output_points[idx].x + tile_curr_x_coord;
+        raystart.y = output_points[idx].y + tile_curr_y_coord;
+
+        data.cur_point_weight = (float)idx / (float)output_points.size();
+
+        BLI_bvhtree_ray_cast_all(
+            treedata.tree, raystart, dir, 0.0f, max_dist, project_2d_bvh_callback, &data);
+      }
+    }
+  }
+
+  return final_points;
+}
+
 static void geo_node_point_distribute_exec(GeoNodeExecParams params)
 {
   GeometrySet geometry_set = params.extract_input<GeometrySet>("Geometry");
   GeometrySet geometry_set_out;
 
+  GeometryNodePointDistributeMethod distribute_method =
+      static_cast<GeometryNodePointDistributeMethod>(params.node().custom1);
+
   if (!geometry_set.has_mesh()) {
     params.set_output("Geometry", std::move(geometry_set_out));
     return;
   }
 
-  const float density = params.extract_input<float>("Density");
+  const float density = params.extract_input<float>("Maximum Density");
   const std::string density_attribute = params.extract_input<std::string>("Density Attribute");
 
   if (density <= 0.0f) {
@@ -113,8 +265,19 @@ static void geo_node_point_distribute_exec(GeoNodeExecParams params)
 
   const FloatReadAttribute density_factors = mesh_component.attribute_get_for_read<float>(
       density_attribute, ATTR_DOMAIN_POINT, 1.0f);
+  const int seed = params.get_input<int>("Seed");
 
-  Vector<float3> points = scatter_points_from_mesh(mesh_in, density, density_factors);
+  Vector<float3> points;
+
+  switch (distribute_method) {
+    case GEO_NODE_POINT_DISTRIBUTE_RANDOM:
+      points = random_scatter_points_from_mesh(mesh_in, density, density_factors, seed);
+      break;
+    case GEO_NODE_POINT_DISTRIBUTE_POISSON:
+      const float min_dist = params.extract_input<float>("Minimum Distance");
+      points = poisson_scatter_points_from_mesh(mesh_in, density, min_dist, density_factors, seed);
+      break;
+  }
 
   PointCloud *pointcloud = BKE_pointcloud_new_nomain(points.size());
   memcpy(pointcloud->co, points.data(), sizeof(float3) * points.size());
@@ -135,6 +298,7 @@ void register_node_type_geo_point_distribute()
   geo_node_type_base(
       &ntype, GEO_NODE_POINT_DISTRIBUTE, "Point Distribute", NODE_CLASS_GEOMETRY, 0);
   node_type_socket_templates(&ntype, geo_node_point_distribute_in, geo_node_point_distribute_out);
+  node_type_update(&ntype, node_point_distribute_update);
   ntype.geometry_node_execute = blender::nodes::geo_node_point_distribute_exec;
   nodeRegisterType(&ntype);
 }
diff --git a/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc b/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc
new file mode 100644
index 00000000000..5dfdfce0e68
--- /dev/null
+++ b/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc
@@ -0,0 +1,280 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/*
+ * Based on Cem Yuksel. 2015. Sample Elimination for Generating Poisson Disk Sample
+ * ! Sets. Computer Graphics Forum 34, 2 (May 2015), 25-32.
+ * ! http://www.cemyuksel.com/research/sampleelimination/
+ * Copyright (c) 2016, Cem Yuksel <cem@cemyuksel.com>
+ * All rights reserved.
+ */
+
+#include "BLI_inplace_priority_queue.hh"
+#include "BLI_kdtree.h"
+
+#include "node_geometry_util.hh"
+
+#include <iostream>
+
+namespace blender::nodes {
+
+static void tile_point(Vector<float3> *tiled_points,
+                       Vector<size_t> *indices,
+                       const float maximum_distance,
+                       const float3 boundbox,
+                       float3 const &point,
+                       size_t index,
+                       int dimension = 0)
+{
+  for (int dimension_iter = dimension; dimension_iter < 3; dimension_iter++) {
+    if (boundbox[dimension_iter] - point[dimension_iter] < maximum_distance) {
+      float3 point_tiled = point;
+      point_tiled[dimension_iter] -= boundbox[dimension_iter];
+
+      tiled_points->append(point_tiled);
+      indices->append(index);
+
+      tile_point(tiled_points,
+                 indices,
+                 maximum_distance,
+                 boundbox,
+                 point_tiled,
+                 index,
+                 dimension_iter + 1);
+    }
+
+    if (point[dimension_iter] < maximum_distance) {
+      float3 point_tiled = point;
+      point_tiled[dimension_iter] += boundbox[dimension_iter];
+
+      tiled_points->append(point_tiled);
+      indices->append(index);
+
+      tile_point(tiled_points,
+                 indices,
+                 maximum_distance,
+                 boundbox,
+                 point_tiled,
+                 index,
+                 dimension_iter + 1);
+    }
+  }
+}
+
+/**
+ * Returns the weight the point gets based on the distance to another point.
+ */
+static float point_weight_influence_get(const float maximum_distance,
+                                        const float minimum_distance,
+                                        float distance)
+{
+  const float alpha = 8.0f;
+
+  if (distance < minimum_distance) {
+    distance = minimum_distance;
+  }
+
+  return std::pow(1.0f - distance / maximum_distance, alpha);
+}
+
+/**
+ * Weight each point based on their proximity to its neighbors
+ *
+ * For each index in the weight array add a weight based on the proximity the
+ * corresponding point has with its neighboors.
+ **/
+static void points_distance_weight_calculate(Vector<float> *weights,
+                                             const size_t point_id,
+                                             const float3 *input_points,
+                                             const void *kd_tree,
+                                             const float minimum_distance,
+                                             const float maximum_distance,
+                                             InplacePriorityQueue<float> *heap)
+{
+  KDTreeNearest_3d *nearest_point = nullptr;
+  int neighbors = BLI_kdtree_3d_range_search(
+      (KDTree_3d *)kd_tree, input_points[point_id], &nearest_point, maximum_distance);
+
+  for (int i = 0; i < neighbors; i++) {
+    size_t neighbor_point_id = nearest_point[i].index;
+
+    if (neighbor_point_id >= weights->size()) {
+      continue;
+    }
+
+    /* The point should not influence itself. */
+    if (neighbor_point_id == point_id) {
+      continue;
+    }
+
+    const float weight_influence = point_weight_influence_get(
+        maximum_distance, minimum_distance, nearest_point[i].dist);
+
+    /* In the first pass we just the weights. */
+    if (heap == nullptr) {
+      (*weights)[point_id] += weight_influence;
+    }
+    /* When we run again we need to update the weights and the heap. */
+    else {
+      (*weights)[neighbor_point_id] -= weight_influence;
+      heap->priority_decreased(neighbor_point_id);
+    }
+  }
+
+  if (nearest_point) {
+    MEM_freeN(nearest_point);
+  }
+}
+
+/**
+ * Returns the minimum radius fraction used by the default weight function.
+ */
+static float weight_limit_fraction_get(const size_t input_size, const size_t output_size)
+{
+  const float beta = 0.65f;
+  const float gamma = 1.5f;
+  float ratio = float(output_size) / float(input_size);
+  return (1.0f - std::pow(ratio, gamma)) * beta;
+}
+
+/**
+ * Tile the input points.
+ */
+static void points_tiling(const float3 *input_points,
+                          const size_t input_size,
+                          void **kd_tree,
+                          const float maximum_distance,
+                          const float3 boundbox)
+
+{
+  Vector<float3> tiled_points(input_points, input_points + input_size);
+  Vector<size_t> indices(input_size);
+
+  for (size_t i = 0; i < input_size; i++) {
+    indices[i] = i;
+  }
+
+  /* Tile the tree based on the boundbox. */
+  for (size_t i = 0; i < input_size; i++) {
+    tile_point(&tiled_points, &indices, maximum_distance, boundbox, input_points[i], i);
+  }
+
+  /* Build a new tree with the new indices and tiled points. */
+  *kd_tree = BLI_kdtree_3d_new(tiled_points.size());
+  for (size_t i = 0; i < tiled_points.size(); i++) {
+    BLI_kdtree_3d_insert(*(KDTree_3d **)kd_tree, indices[i], tiled_points[i]);
+  }
+  BLI_kdtree_3d_balance(*(KDTree_3d **)kd_tree);
+}
+
+static void weighted_sample_elimination(const float3 *input_points,
+                                        const size_t input_size,
+                                        float3 *output_points,
+                                        const size_t output_size,
+                                        const float maximum_distance,
+                                        const float3 boundbox,
+                                        const bool do_copy_eliminated)
+{
+  const float minimum_distance = maximum_distance *
+                                 weight_limit_fraction_get(input_size, output_size);
+
+  void *kd_tree = nullptr;
+  points_tiling(input_points, input_size, &kd_tree, maximum_distance, boundbox);
+
+  /* Assign weights to each sample. */
+  Vector<float> weights(input_size, 0.0f);
+  for (size_t point_id = 0; point_id < weights.size(); point_id++) {
+    points_distance_weight_calculate(
+        &weights, point_id, input_points, kd_tree, minimum_distance, maximum_distance, nullptr);
+  }
+
+  /* Remove the points based on their weight. */
+  InplacePriorityQueue<float> heap(weights);
+
+  size_t sample_size = input_size;
+  while (sample_size > output_size) {
+    /* For each sample around it, remove its weight contribution and update the heap. */
+    size_t point_id = heap.pop_index();
+    points_distance_weight_calculate(
+        &weights, point_id, input_points, kd_tree, minimum_distance, maximum_distance, &heap);
+    sample_size--;
+  }
+
+  /* Copy the samples to the output array. */
+  size_t target_size = do_copy_eliminated ? input_size : output_size;
+  for (size_t i = 0; i < target_size; i++) {
+    size_t index = heap.all_indices()[i];
+    output_points[i] = input_points[index];
+  }
+
+  /* Cleanup. */
+  BLI_kdtree_3d_free((KDTree_3d *)kd_tree);
+}
+
+static void progressive_sampling_reorder(Vector<float3> *output_points,
+                                         float maximum_density,
+                                         float3 boundbox)
+{
+  /* Re-order the points for progressive sampling. */
+  Vector<float3> temporary_points(output_points->size());
+  float3 *source_points = output_points->data();
+  float3 *dest_points = temporary_points.data();
+  size_t source_size = output_points->size();
+  size_t dest_size = 0;
+
+  while (source_size >= 3) {
+    dest_size = source_size * 0.5f;
+
+    /* Changes the weight function radius using half of the number of samples.
+     * It is used for progressive sampling. */
+    maximum_density *= std::sqrt(2.0f);
+    weighted_sample_elimination(
+        source_points, source_size, dest_points, dest_size, maximum_density, boundbox, true);
+
+    if (dest_points != output_points->data()) {
+      mempcpy((*output_points)[dest_size],
+              dest_points[dest_size],
+              (source_size - dest_size) * sizeof(float3));
+    }
+
+    /* Swap the arrays around. */
+    float3 *points_iter = source_points;
+    source_points = dest_points;
+    dest_points = points_iter;
+    source_size = dest_size;
+  }
+  if (source_points != output_points->data()) {
+    memcpy(output_points->data(), source_points, dest_size * sizeof(float3));
+  }
+}
+
+void poisson_disk_point_elimination(Vector<float3> const *input_points,
+                                    Vector<float3> *output_points,
+                                    float maximum_density,
+                                    float3 boundbox)
+{
+  weighted_sample_elimination(input_points->data(),
+                              input_points->size(),
+                              output_points->data(),
+                              output_points->size(),
+                              maximum_density,
+                              boundbox,
+                              false);
+
+  progressive_sampling_reorder(output_points, maximum_density, boundbox);
+}
+
+}  // namespace blender::nodes