**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:
{F9548872}
This patch, sampling directly on mesh surface:
{F9548873}
Yuksel's algorithm (weighted sample elimination)
=============================================
The both implementations are based on the algorithm described in: [[ http://www.cemyuksel.com/research/sampleelimination/sampleelimination.pdf | 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:
{F9557037}
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.
{F9557046}
I mention this because Sebastian stated this as the reason why having a `Minimum Distance` input is desirable for users:
> 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:
{F9557084}
Yuksel gives formulas to calculate `r_max` for 2D and 3D:
{F9557086}
But this left me puzzled: which formula do I use for a 2D mesh surface in 3D space? On his [[ http://www.cemyuksel.com/cyCodeBase/soln/poisson_disk_sampling.html | website ]] (section 7.) he states:
> [...] 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:
# Initial weight calculation
# 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
{F9558558}
In the dropdown of Point Distribute node, I've used the name "Poisson Disk" for the method from this patch, and added an option "Poisson Disk 2D" for the other method:
{F9558544}
Note that the new method currently only has a `Density` slider and no `Minimum Distance` like the 2D version. The reason is just that it was simpler to implementation wise, also I some ways I find it more user friendly. When using the 2D method I often find myself going back and forth between the `Minimum Distance` slider and the `Density` slider and being frustrated so I really think this design should be reconsidered (sorry for bringing this up again though). Maybe a solution is to 1) also use Density as the main driver for the 2D node 2) have a checkbox `Ensure minimum distance` for when you're scatter things with a known fixed radius that really should not overlap?
Weight painting:
{F9558893}
Planned changes before patch is finished:
- There's no integration yet with vertex painting yet, but this should not be that hard. One solution that should be quite fast is just run the algorithm and then afterwards filter points per triangle based on the average vertex weight. An important optimization could be to eliminate triangles with zero weight before the algorithm runs. This could be very impactful when a user only want particles on a small part of the mesh, e.g. on top of the head of a person.
- 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 Elimination`part 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!