Particles: spawn particles on mesh surface

8 files changed, 403 insertions, 46 deletions

diff --git a/source/blender/blenkernel/BKE_persistent_data_handle.hh b/source/blender/blenkernel/BKE_persistent_data_handle.hh
index 721132560e3..884e4c00766 100644
--- a/source/blender/blenkernel/BKE_persistent_data_handle.hh
+++ b/source/blender/blenkernel/BKE_persistent_data_handle.hh
@@ -85,45 +85,45 @@ class PersistentObjectHandle : public PersistentIDHandle {
 
 class PersistentDataHandleMap {
  private:
-  Map<int32_t, const ID *> id_by_handle_;
-  Map<const ID *, int32_t> handle_by_id_;
+  Map<int32_t, ID *> id_by_handle_;
+  Map<ID *, int32_t> handle_by_id_;
 
  public:
-  void add(int32_t handle, const ID &id)
+  void add(int32_t handle, ID &id)
   {
     BLI_assert(handle >= 0);
     handle_by_id_.add(&id, handle);
     id_by_handle_.add(handle, &id);
   }
 
-  PersistentIDHandle lookup(const ID *id) const
+  PersistentIDHandle lookup(ID *id) const
   {
     const int handle = handle_by_id_.lookup_default(id, -1);
     return PersistentIDHandle(handle);
   }
 
-  PersistentObjectHandle lookup(const Object *object) const
+  PersistentObjectHandle lookup(Object *object) const
   {
-    const int handle = handle_by_id_.lookup_default((const ID *)object, -1);
+    const int handle = handle_by_id_.lookup_default((ID *)object, -1);
     return PersistentObjectHandle(handle);
   }
 
-  const ID *lookup(const PersistentIDHandle &handle) const
+  ID *lookup(const PersistentIDHandle &handle) const
   {
-    const ID *id = id_by_handle_.lookup_default(handle.handle_, nullptr);
+    ID *id = id_by_handle_.lookup_default(handle.handle_, nullptr);
     return id;
   }
 
-  const Object *lookup(const PersistentObjectHandle &handle) const
+  Object *lookup(const PersistentObjectHandle &handle) const
   {
-    const ID *id = this->lookup((const PersistentIDHandle &)handle);
+    ID *id = this->lookup((const PersistentIDHandle &)handle);
     if (id == nullptr) {
       return nullptr;
     }
     if (GS(id->name) != ID_OB) {
       return nullptr;
     }
-    return (const Object *)id;
+    return (Object *)id;
   }
 };
 
diff --git a/source/blender/blenlib/BLI_float3.hh b/source/blender/blenlib/BLI_float3.hh
index b2633985ac7..13ec563578b 100644
--- a/source/blender/blenlib/BLI_float3.hh
+++ b/source/blender/blenlib/BLI_float3.hh
@@ -143,6 +143,17 @@ struct float3 {
     return normalize_v3(*this);
   }
 
+  /**
+   * Normalizes the vector inplace.
+   */
+  void normalize()
+  {
+    normalize_v3(*this);
+  }
+
+  /**
+   * Returns a normalized vector. The original vector is not changed.
+   */
   float3 normalized() const
   {
     float3 result;
diff --git a/source/blender/blenlib/BLI_float4x4.hh b/source/blender/blenlib/BLI_float4x4.hh
index 185cffd13ac..5698d06f05d 100644
--- a/source/blender/blenlib/BLI_float4x4.hh
+++ b/source/blender/blenlib/BLI_float4x4.hh
@@ -71,8 +71,8 @@ struct float4x4 {
 
   float4x4 inverted() const
   {
-    float result[4][4];
-    invert_m4_m4(result, values);
+    float4x4 result;
+    invert_m4_m4(result.values, values);
     return result;
   }
 
@@ -86,6 +86,13 @@ struct float4x4 {
     return this->inverted();
   }
 
+  float4x4 transposed() const
+  {
+    float4x4 result;
+    transpose_m4_m4(result.values, values);
+    return result;
+  }
+
   struct float3x3_ref {
     const float4x4 &data;
 
diff --git a/source/blender/blenlib/BLI_rand.hh b/source/blender/blenlib/BLI_rand.hh
index 612ac0bbe19..7a98ee0f2bb 100644
--- a/source/blender/blenlib/BLI_rand.hh
+++ b/source/blender/blenlib/BLI_rand.hh
@@ -63,6 +63,15 @@ class RandomNumberGenerator {
   }
 
   /**
+   * \return Random value (0..N), but never N.
+   */
+  int32_t get_int32(int32_t max_exclusive)
+  {
+    BLI_assert(max_exclusive > 0);
+    return this->get_int32() % max_exclusive;
+  }
+
+  /**
    * \return Random value (0..1), but never 1.0.
    */
   double get_double()
@@ -78,6 +87,35 @@ class RandomNumberGenerator {
     return (float)this->get_int32() / 0x80000000;
   }
 
+  template<typename T> void shuffle(MutableSpan<T> values)
+  {
+    /* Cannot shuffle arrays of this size yet. */
+    BLI_assert(values.size() <= INT32_MAX);
+
+    for (int i = values.size() - 1; i >= 2; i--) {
+      int j = this->get_int32(i);
+      if (i != j) {
+        std::swap(values[i], values[j]);
+      }
+    }
+  }
+
+  /**
+   * Compute uniformly distributed barycentric coordinates.
+   */
+  float3 get_barycentric_coordinates()
+  {
+    float rand1 = this->get_float();
+    float rand2 = this->get_float();
+
+    if (rand1 + rand2 > 1.0f) {
+      rand1 = 1.0f - rand1;
+      rand2 = 1.0f - rand2;
+    }
+
+    return float3(rand1, rand2, 1.0f - rand1 - rand2);
+  }
+
   float2 get_unit_float2();
   float3 get_unit_float3();
   float2 get_triangle_sample(float2 v1, float2 v2, float2 v3);
diff --git a/source/blender/blenlib/BLI_vector_adaptor.hh b/source/blender/blenlib/BLI_vector_adaptor.hh
new file mode 100644
index 00000000000..cadffc0b445
--- /dev/null
+++ b/source/blender/blenlib/BLI_vector_adaptor.hh
@@ -0,0 +1,105 @@
+/*
+ * 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.
+ */
+
+#ifndef __BLI_VECTOR_ADAPTOR_HH__
+#define __BLI_VECTOR_ADAPTOR_HH__
+
+/** \file
+ * \ingroup bli
+ *
+ * A `blender::VectorAdaptor` is a container with a fixed maximum size and does not own the
+ * underlying memory. When an adaptor is constructed, you have to provide it with an uninitialized
+ * array that will be filled when elements are added to the vector. The vector adaptor is not able
+ * to grow. Therefore, it is undefined behavior to add more elements than fit into the provided
+ * buffer.
+ */
+
+#include "BLI_span.hh"
+
+namespace blender {
+
+template<typename T> class VectorAdaptor {
+ private:
+  T *begin_;
+  T *end_;
+  T *capacity_end_;
+
+ public:
+  VectorAdaptor() : begin_(nullptr), end_(nullptr), capacity_end_(nullptr)
+  {
+  }
+
+  VectorAdaptor(T *data, int64_t capacity, int64_t size = 0)
+      : begin_(data), end_(data + size), capacity_end_(data + capacity)
+  {
+  }
+
+  VectorAdaptor(MutableSpan<T> span) : VectorAdaptor(span.data(), span.size(), 0)
+  {
+  }
+
+  void append(const T &value)
+  {
+    BLI_assert(end_ < capacity_end_);
+    new (end_) T(value);
+    end_++;
+  }
+
+  void append(T &&value)
+  {
+    BLI_assert(end_ < capacity_end_);
+    new (end_) T(std::move(value));
+    end_++;
+  }
+
+  void append_n_times(const T &value, int64_t n)
+  {
+    BLI_assert(end_ + n <= capacity_end_);
+    uninitialized_fill_n(end_, n, value);
+    end_ += n;
+  }
+
+  void extend(Span<T> values)
+  {
+    BLI_assert(end_ + values.size() <= capacity_end_);
+    uninitialized_copy_n(values.data(), values.size(), end_);
+    end_ += values.size();
+  }
+
+  int64_t capacity() const
+  {
+    return capacity_end_ - begin_;
+  }
+
+  int64_t size() const
+  {
+    return end_ - begin_;
+  }
+
+  bool is_empty() const
+  {
+    return begin_ == end_;
+  }
+
+  bool is_full() const
+  {
+    return end_ == capacity_end_;
+  }
+};
+
+}  // namespace blender
+
+#endif /* __BLI_VECTOR_ADAPTOR_HH__ */
diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt
index d7d1d701bc8..db44ebe8e55 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -267,6 +267,7 @@ set(SRC
   BLI_utility_mixins.hh
   BLI_uvproject.h
   BLI_vector.hh
+  BLI_vector_adaptor.hh
   BLI_vector_set.hh
   BLI_vector_set_slots.hh
   BLI_vfontdata.h
diff --git a/source/blender/nodes/intern/node_socket.cc b/source/blender/nodes/intern/node_socket.cc
index 49868505d68..57c8d51975e 100644
--- a/source/blender/nodes/intern/node_socket.cc
+++ b/source/blender/nodes/intern/node_socket.cc
@@ -586,10 +586,10 @@ static bNodeSocketType *make_socket_type_string()
 
 class ObjectSocketMultiFunction : public blender::fn::MultiFunction {
  private:
-  const Object *object_;
+  Object *object_;
 
  public:
-  ObjectSocketMultiFunction(const Object *object) : object_(object)
+  ObjectSocketMultiFunction(Object *object) : object_(object)
   {
     blender::fn::MFSignatureBuilder signature = this->get_builder("Object Socket");
     signature.depends_on_context();
@@ -631,7 +631,7 @@ static bNodeSocketType *make_socket_type_object()
     return blender::fn::MFDataType::ForSingle<blender::bke::PersistentObjectHandle>();
   };
   socktype->expand_in_mf_network = [](blender::nodes::SocketMFNetworkBuilder &builder) {
-    const Object *object = builder.socket_default_value<bNodeSocketValueObject>()->value;
+    Object *object = builder.socket_default_value<bNodeSocketValueObject>()->value;
     builder.construct_generator_fn<ObjectSocketMultiFunction>(object);
   };
   return socktype;
diff --git a/source/blender/simulation/intern/particle_mesh_emitter.cc b/source/blender/simulation/intern/particle_mesh_emitter.cc
index ab2ee0c81ed..486cd1142c3 100644
--- a/source/blender/simulation/intern/particle_mesh_emitter.cc
+++ b/source/blender/simulation/intern/particle_mesh_emitter.cc
@@ -18,6 +18,9 @@
 
 #include "BLI_float4x4.hh"
 #include "BLI_rand.hh"
+#include "BLI_vector_adaptor.hh"
+
+#include "BKE_mesh_runtime.h"
 
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
@@ -26,14 +29,14 @@
 namespace blender::sim {
 
 struct EmitterSettings {
-  const Object *object;
+  Object *object;
   float rate;
 };
 
-static void compute_birth_times(float rate,
-                                TimeInterval emit_interval,
-                                ParticleMeshEmitterSimulationState &state,
-                                Vector<float> &r_birth_times)
+static BLI_NOINLINE void compute_birth_times(float rate,
+                                             TimeInterval emit_interval,
+                                             ParticleMeshEmitterSimulationState &state,
+                                             Vector<float> &r_birth_times)
 {
   const float time_between_particles = 1.0f / rate;
   int counter = 0;
@@ -51,51 +54,237 @@ static void compute_birth_times(float rate,
   }
 }
 
-static void compute_new_particle_attributes(EmitterSettings &settings,
-                                            TimeInterval emit_interval,
-                                            ParticleMeshEmitterSimulationState &state,
-                                            Vector<float3> &r_positions,
-                                            Vector<float3> &r_velocities,
-                                            Vector<float> &r_birth_times)
+static BLI_NOINLINE Span<MLoopTri> get_mesh_triangles(Mesh &mesh)
 {
-  if (settings.object == nullptr) {
+  const MLoopTri *triangles = BKE_mesh_runtime_looptri_ensure(&mesh);
+  int amount = BKE_mesh_runtime_looptri_len(&mesh);
+  return Span(triangles, amount);
+}
+
+static BLI_NOINLINE void compute_triangle_areas(Mesh &mesh,
+                                                Span<MLoopTri> triangles,
+                                                MutableSpan<float> r_areas)
+{
+  assert_same_size(triangles, r_areas);
+
+  for (int i : triangles.index_range()) {
+    const MLoopTri &tri = triangles[i];
+
+    const float3 v1 = mesh.mvert[mesh.mloop[tri.tri[0]].v].co;
+    const float3 v2 = mesh.mvert[mesh.mloop[tri.tri[1]].v].co;
+    const float3 v3 = mesh.mvert[mesh.mloop[tri.tri[2]].v].co;
+
+    const float area = area_tri_v3(v1, v2, v3);
+    r_areas[i] = area;
+  }
+}
+
+static BLI_NOINLINE void compute_triangle_weights(Mesh &mesh,
+                                                  Span<MLoopTri> triangles,
+                                                  MutableSpan<float> r_weights)
+{
+  assert_same_size(triangles, r_weights);
+  compute_triangle_areas(mesh, triangles, r_weights);
+}
+
+static BLI_NOINLINE void compute_cumulative_distribution(Span<float> weights,
+                                                         MutableSpan<float> r_cumulative_weights)
+{
+  BLI_assert(weights.size() + 1 == r_cumulative_weights.size());
+
+  r_cumulative_weights[0] = 0;
+  for (int i : weights.index_range()) {
+    r_cumulative_weights[i + 1] = r_cumulative_weights[i] + weights[i];
+  }
+}
+
+static void sample_cumulative_distribution_recursive(RandomNumberGenerator &rng,
+                                                     int amount,
+                                                     int start,
+                                                     int one_after_end,
+                                                     Span<float> cumulative_weights,
+                                                     VectorAdaptor<int> &r_sampled_indices)
+{
+  BLI_assert(start <= one_after_end);
+  const int size = one_after_end - start;
+  if (size == 0) {
+    BLI_assert(amount == 0);
+  }
+  else if (amount == 0) {
     return;
   }
+  else if (size == 1) {
+    r_sampled_indices.append_n_times(start, amount);
+  }
+  else {
+    const int middle = start + size / 2;
+    const float left_weight = cumulative_weights[middle] - cumulative_weights[start];
+    const float right_weight = cumulative_weights[one_after_end] - cumulative_weights[middle];
+    BLI_assert(left_weight >= 0.0f && right_weight >= 0.0f);
+    const float weight_sum = left_weight + right_weight;
+    BLI_assert(weight_sum > 0.0f);
+
+    const float left_factor = left_weight / weight_sum;
+    const float right_factor = right_weight / weight_sum;
+
+    int left_amount = amount * left_factor;
+    int right_amount = amount * right_factor;
+
+    if (left_amount + right_amount < amount) {
+      BLI_assert(left_amount + right_amount + 1 == amount);
+      const float weight_per_item = weight_sum / amount;
+      const float total_remaining_weight = weight_sum -
+                                           (left_amount + right_amount) * weight_per_item;
+      const float left_remaining_weight = left_weight - left_amount * weight_per_item;
+      const float left_remaining_factor = left_remaining_weight / total_remaining_weight;
+      if (rng.get_float() < left_remaining_factor) {
+        left_amount++;
+      }
+      else {
+        right_amount++;
+      }
+    }
+
+    sample_cumulative_distribution_recursive(
+        rng, left_amount, start, middle, cumulative_weights, r_sampled_indices);
+    sample_cumulative_distribution_recursive(
+        rng, right_amount, middle, one_after_end, cumulative_weights, r_sampled_indices);
+  }
+}
+
+static BLI_NOINLINE void sample_cumulative_distribution(RandomNumberGenerator &rng,
+                                                        Span<float> cumulative_weights,
+                                                        MutableSpan<int> r_samples)
+{
+  VectorAdaptor<int> sampled_indices(r_samples);
+  sample_cumulative_distribution_recursive(rng,
+                                           r_samples.size(),
+                                           0,
+                                           cumulative_weights.size() - 1,
+                                           cumulative_weights,
+                                           sampled_indices);
+  BLI_assert(sampled_indices.is_full());
+}
+
+static BLI_NOINLINE bool sample_weighted_buckets(RandomNumberGenerator &rng,
+                                                 Span<float> weights,
+                                                 MutableSpan<int> r_samples)
+{
+  Array<float> cumulative_weights(weights.size() + 1);
+  compute_cumulative_distribution(weights, cumulative_weights);
+
+  if (r_samples.size() > 0 && cumulative_weights.as_span().last() == 0.0f) {
+    /* All weights are zero. */
+    return false;
+  }
+
+  sample_cumulative_distribution(rng, cumulative_weights, r_samples);
+  return true;
+}
+
+static BLI_NOINLINE void sample_looptris(RandomNumberGenerator &rng,
+                                         Mesh &mesh,
+                                         Span<MLoopTri> triangles,
+                                         Span<int> triangles_to_sample,
+                                         MutableSpan<float3> r_sampled_positions,
+                                         MutableSpan<float3> r_sampled_normals)
+{
+  assert_same_size(triangles_to_sample, r_sampled_positions, r_sampled_normals);
+
+  MLoop *loops = mesh.mloop;
+  MVert *verts = mesh.mvert;
+
+  for (uint i : triangles_to_sample.index_range()) {
+    const uint triangle_index = triangles_to_sample[i];
+    const MLoopTri &triangle = triangles[triangle_index];
+
+    const float3 v1 = verts[loops[triangle.tri[0]].v].co;
+    const float3 v2 = verts[loops[triangle.tri[1]].v].co;
+    const float3 v3 = verts[loops[triangle.tri[2]].v].co;
+
+    const float3 bary_coords = rng.get_barycentric_coordinates();
+
+    float3 position;
+    interp_v3_v3v3v3(position, v1, v2, v3, bary_coords);
+
+    float3 normal;
+    normal_tri_v3(normal, v1, v2, v3);
+
+    r_sampled_positions[i] = position;
+    r_sampled_normals[i] = normal;
+  }
+}
+
+static BLI_NOINLINE bool compute_new_particle_attributes(EmitterSettings &settings,
+                                                         TimeInterval emit_interval,
+                                                         ParticleMeshEmitterSimulationState &state,
+                                                         Vector<float3> &r_positions,
+                                                         Vector<float3> &r_velocities,
+                                                         Vector<float> &r_birth_times)
+{
+  if (settings.object == nullptr) {
+    return false;
+  }
   if (settings.rate <= 0.000001f) {
-    return;
+    return false;
   }
   if (settings.object->type != OB_MESH) {
-    return;
+    return false;
   }
-  const Mesh &mesh = *(Mesh *)settings.object->data;
+  Mesh &mesh = *(Mesh *)settings.object->data;
   if (mesh.totvert == 0) {
-    return;
+    return false;
   }
 
-  float4x4 local_to_world = settings.object->obmat;
+  const float4x4 local_to_world = settings.object->obmat;
+  float4x4 local_to_world_normal = local_to_world.inverted_affine().transposed();
 
   const float start_time = emit_interval.start();
   const uint32_t seed = DefaultHash<StringRef>{}(state.head.name);
   RandomNumberGenerator rng{(*(uint32_t *)&start_time) ^ seed};
 
   compute_birth_times(settings.rate, emit_interval, state, r_birth_times);
-  if (r_birth_times.is_empty()) {
-    return;
+  const int particle_amount = r_birth_times.size();
+  if (particle_amount == 0) {
+    return false;
   }
 
-  state.last_birth_time = r_birth_times.last();
+  rng.shuffle(r_birth_times.as_mutable_span());
+
+  Span<MLoopTri> triangles = get_mesh_triangles(mesh);
+  if (triangles.is_empty()) {
+    return false;
+  }
+
+  Array<float> triangle_weights(triangles.size());
+  compute_triangle_weights(mesh, triangles, triangle_weights);
 
-  for (int i : r_birth_times.index_range()) {
-    UNUSED_VARS(i);
-    const int vertex_index = rng.get_int32() % mesh.totvert;
-    float3 vertex_position = mesh.mvert[vertex_index].co;
-    r_positions.append(local_to_world * vertex_position);
-    r_velocities.append(rng.get_unit_float3());
+  Array<int> triangles_to_sample(particle_amount);
+  if (!sample_weighted_buckets(rng, triangle_weights, triangles_to_sample)) {
+    return false;
   }
+
+  Array<float3> local_positions(particle_amount);
+  Array<float3> local_normals(particle_amount);
+  sample_looptris(rng, mesh, triangles, triangles_to_sample, local_positions, local_normals);
+
+  r_positions.reserve(particle_amount);
+  r_velocities.reserve(particle_amount);
+  for (int i : IndexRange(particle_amount)) {
+    float3 position = local_to_world * local_positions[i];
+    float3 normal = local_to_world_normal.ref_3x3() * local_normals[i];
+    normal.normalize();
+
+    r_positions.append(position);
+    r_velocities.append(normal);
+  }
+
+  state.last_birth_time = r_birth_times.last();
+  return true;
 }
 
-static EmitterSettings compute_settings(const fn::MultiFunction &inputs_fn,
-                                        ParticleEmitterContext &context)
+static BLI_NOINLINE EmitterSettings compute_settings(const fn::MultiFunction &inputs_fn,
+                                                     ParticleEmitterContext &context)
 {
   EmitterSettings parameters;
 
@@ -126,8 +315,14 @@ void ParticleMeshEmitter::emit(ParticleEmitterContext &context) const
   Vector<float3> new_velocities;
   Vector<float> new_birth_times;
 
-  compute_new_particle_attributes(
-      settings, context.emit_interval(), *state, new_positions, new_velocities, new_birth_times);
+  if (!compute_new_particle_attributes(settings,
+                                       context.emit_interval(),
+                                       *state,
+                                       new_positions,
+                                       new_velocities,
+                                       new_birth_times)) {
+    return;
+  }
 
   for (StringRef name : particle_names_) {
     ParticleAllocator *allocator = context.try_get_particle_allocator(name);