Update: Jacques added a great Poisson Disk sampling method that also samples directly onto the mesh that will be incorporated into 2.92 and satisfies most use cases. His implementation is similar to what in the literature they call "dart throwing", but instead of progressively adding samples which would require many KDTree insertions (and balancing), he adds all samples upfront and then removes samples that violate the minimum distance. This is actually a very elegant approach because it circumvents the issue of having to update a growing datastructure for neighborhood lookups.

The method in this patch (and the old implementation) were based on Yuksel's sample elimination. While Yuksel's algorithm is also quite easy to implement, it is harder to understand e.g. setting a good value of r_max is a surprisingly difficult (where good means that the resulting distribution will be nice and that algorithms runs quickly). So for the moment the code in this patch is put in freezer. We might decide to revive Yuksel's method for two main reasons:

1. It could produce Poisson Disks of higher quality = with denser packing for a given Poisson Disk radius than dart throwing (although probably slower). For some use cases this slightly denser packing/slightly smaller gaps might be desirable, so this alone might justify implementing this patch for 2.93. However I find the second reason even more appealling:
2. Yuksel's algorithm can be used to generate different distributions based on the weight function. For example, in his paper he shows how you can easily get magenta noise. By cleverly parametrizing the weight function, we could expose a slider to the user called "Clustering" or something similar. At 0 cluster, you would get a Poisson Disk point set, at slightly higher values you would get magenta noise and at even higher values the points would completely separate into several "islands" or clustered points.

Below this point, the content of this post is outdated.

Parts of this post:

• Problem statement
• Yuksel's algorithm (weighted sample elimination)
• Current implementation: sampling in 2D + projection
• Sampling directly on the mesh
• Patch status

Problem Statement

The goal of Poisson disk sampling a set of random points where no 2 points are closer than a certain minimum distance. The radius of a Poisson disk set is half the distance between the 2 closest points of the set.

The current implementation samples the points in 2D (parallel to the XY plane) and then projects the along the Z axis onto the mesh. This creates the unexpected (and undesirable in most cases?) artefact that "steep" faces have less points.

This patch introduces Poisson disk sampling directly on the input mesh surface. However it is still very slow, which is why the original developers decided to do the sampling in 2D as a first step. Below I will explain why sampling directly on the mesh is slow and what we can do about it. Before that I will also first explain the underlying algorithm and the current 2D implementation in a bit more depth.

Currently implemented, 2D sampling + Z-projection:

This patch, sampling directly on mesh surface:

Yuksel's algorithm (weighted sample elimination)

The both implementations are based on the algorithm described in: Sample Elimination for Generating Poisson Disk Sample Sets (Yuksel 2015). The main idea of Yuksel's algorithm is to first generate 5x too many random points. Then calculate a weight for each point that captures how crowded its neighborhood is. Finally, remove the points with highest weight until you're happy with the result. In pseudocode the algorithm goes something like this:

points = rng.get_random_points(5 * desired_points_amount)
r_max = ...
weights = []

// Initial weight assignment
foreach(point in points) {
  weight = 0.0
  neighbors = find_neighbors_within_radius(point, r_max)
  foreach(neighbor in neighors) {
    distance = distance(point, neighbor)
    weight += calculate_weight(distance, r_max)
  }
}

// Point elimination
heap = Heap(points, weights)
while(heap.size() > desired_points_amount) {
  point = heap.pop()
  neighbors = kdtree.find_neighbors_within_radius(point, r_max)
  foreach(neighbor in neighbors) {
    distance = distance(point, neighbor)
    weight = calculate_weight(distance, r_max)
    heap.decrease_priority(neighbor, weight)
  }
}
remaining_points = heap.items()
return remaining_points
}

Current implementation: sampling in 2D + projection

The current (2D) implementation uses Yuksel's algorithm to generate 5000 points in a square. The size of this square depends only on the minimum distance specified by the user. When points need to be distributed over a larger area than the square, the square is simply tiled in the X and Y directions to cover the larger area. Then these tiled points are projected onto the input mesh.

Ok so how does this method handle the density input? To understand that we need to understand progressive sample reordering. Progressive sample reordering sorts the set of generated Poisson disk points such that when points are added in that order, each subset will look like a Poisson disk set (but with larger radius). After ordering the generated points in that way, all points at the end of the list are filtered to achieve the desired density.

What's nice about this approach is that it's very fast and produces stable results. For a given radius and random seed, it produces the same tiling of the XY-plane. These points are then sorted in the same way each time. When the user changes the density, we just take a smaller/larger part of the list. If the user resizes the mesh to have a larger XY bounding box, simply more tiles are added, without changing the existing tiles.

However for small radius + large mesh XY bounding box, tiling artefacts appear:

I also want to note the following about the current 2D implementation: it does not guarantee the user specified minimum distance for high density. Below is an example of scattering spheres with radius 0.5 using minimum distance 0.5, the spheres overlap in several places.

The idea with these parameters is that you have models/assets that have a certain predetermined radius. So you plug in the the min radius of the assets in the Minimum Distance and tweak the max density to get as sparse/tightly packed distribution as you like <without worrying about overlap between the assets>.

(I added that last part for clarification.)

Sampling directly on the mesh

Generating a set of Poisson disk points directly on the mesh is conceptually simple. First generate 5x too many random points on the mesh (which was already implemented). Then just run Yuksel's algorithm with those points. Yuksel proposes using a KDTree for the neighborhood lookup. I tried it with blender's KDTree implementation, but this was much too slow (e.g. 20k points took 15s). Several optimizations were implemented that together achieved a 40x speedup.

Optimizations

All optimizations in this section revolve around the neighborhood queries. As these were the most expensive part of the algorithm.

Better approximation of r_max for meshes

r_max is used throughout Yuksel's algorithm. It determines the size of neighborhoods and influences the weights assigned to points. Yuksel mentions r_max is "the maximum possible distance for N samples to cover the sampling domain". This maximum distance is achievable with hexagonal packing.

Yuksel gives formulas to calculate r_max for 2D and 3D:

[...] since we are sampling a surface in 3D, the sampling domain should be considered 2D.

I find this a bit weird, because using the 2D formula with the surface area of the mesh will generally underestimate the achievable Poisson disk radius, I thought you ought to be conservative and always use an overestimation. However, lower values are faster (smaller neighborhoods to search) and the underestimation doesn't seem to give any problems so I decided to also use it. I actually even went a step further and used the min(r_max_2D, r_max_3D) where r_max_3D is calculated with the AABB of the mesh.

Grid-based neighborhood search

The KDTree neighborhood queries turned out to be a huge bottleneck, so the KDTree was replaced by a spatial grid with cell size = 2 * r_max. This guarantees that all neighborhoods for a point in grid cell (i, j, k) are in the neighboring cells (e.g. (i+1, j, k), (i+1, j+1, k), ...). So in total 27 cells have to be searched for each point. This might seem inefficient, but for is surprisingly fast: compared to the KDTree this already made the algorithm about 15x faster with better scaling.

Additionally, for each point its neighbors' ids are stored together with the weight for that neighbor. This way, we don't search the neighbors again/recalculate the weights in the point elimination stage of the algorithm.

Parallelizing the neighborhood queries

The algorithm consists of two expensive parts:

1. Initial weight calculation
2. Point elimination

For the initial weight assignment there's no dependencies on other points' weights that loop was fully parallelized.

The second part is more tricky as for each point removed, the points in the vicinity of the removed one need to have their weight values updated. Multi-threading this would require clever sectioning of points (threads should not be able to eliminate points whose weights are being updated by other threads).

Patch Status

In my opinion, the new Poisson Disk sampling on the mesh is certainly fast enough to be usable (generating 100K points takes about 1,8 seconds, but that's already a lot points, about the same as the average amount of hairs on a human head). In some cases I found it even to be faster than the 2D implementation, however with the same level of optimization, the 2D approach will always be faster because you can't beat the tiling :p

Weight painting:

Planned changes before patch is finished:

• Currently I'm calculating the mesh bounding box and surface area just by looping over vertices/triangles. This is actually quite fast and even parallelizable if necessary, so it's not really a problem but I was wondering whether there's a more standard way to do this in Blender. Also for meshes with very high vertex/triangle counts this could become non-negligible.
  • ...

Known limitations:

• The generated points are not stable, i.e. changing the density causes the points to "jump around". Note however that the points are stable during weight painting. Points smoothly pop into and out of existance, without the other points on the mesh moving around.

Future work:

• Parallelize the Point Eliminationpart of the algorithm. This part is still single threaded, so when scattering a lot of points, this becomes the bottleneck. I've talked about this with Sebastian, what he proposes is splitting up space into parts that can be processed in parallel without interference. However combining the results and handling the borders does seem quite difficult so I don't have any concrete plans to implement this.
• Think about how to speed up the algorithm if the user only weight paints small areas. Currently, points are always scattered over the entire mesh, then filtered.

Testing/reviews are very welcome!