Surprisingly Mathematica resources have a lack examples of this effective method. In fact, all I found was my imperfect demonstration. Here I will try to fill in this gap and present a sufficiently effective code for generating Poisson disk samples with Mathematica.

Introduction

It is well known that purely random points are unevenly filling the region. However if you use strictly periodic lattice, then the results of such sampling may be incorrect. Let’s take an annulus region (it will be useful for us in the future) and look at results of built-ins RandomPoint and DiscretizeRegion:
In[]:=
testRegion=Annulus[{0,0},{0.5,1.}];​​numSamples=500;​​randomSamples=RandomPoint[testRegion,numSamples];​​meshSamples=MeshCoordinates@DiscretizeRegion[testRegion,​​MaxCellMeasureArea@testRegion/numSamples];​​Row[Graphics[Point@#,ImageSize->150]&/@{randomSamples,meshSamples}," "]
Out[]=
Honestly, I don’t know exactly what algorithm is used in DiscretizeRegion, but it's result looks evenly distributed. However, it is not convenient to cover regions of random samples. The golden ratio algorithm is suitable for the disk and ring. The points are located on spirals with the GoldenRatio and at each generation are rotated at a random angle:
In[]:=
ClearAll[annulusMesh];​​annulusMesh[cnt_,minR_,maxR_,maxRnd_,num_]:=​​ Withq=3-
5
π, ω=RandomReal[{-maxRnd,maxRnd}],​​ α=(maxR^2-minR^2)num,β=minR^2​​ ,Tablecnt+
α j+β
{Cos[q j+ω],Sin[q j+ω]},{j,num};
In[]:=
Row[Graphics[Point@#,ImageSize->150,PlotRange->{{-1,1},{-1,1}}]&/@​​{randomSamples,meshSamples,annulusMesh[{0,0},0.5,1,0.0,numSamples]}," "]
Out[]=
This method is very specific and not universal. But we will need it further.

Fast Poisson Disk Sampling

Approach presented here is based on the algorithm described in this article. It uses background grid (grid variable) for storing samples and accelerating spatial searches. In the code below, pure random test points are replaced by the golden ratio distribution.

Main module

Input:
◼
  • numTries — number of test points, usually 10 < numTries < 30
    • ◼
  • region — the region in which distribution is formed. The more complex the region, the longer the generation.
    • ◼
  • radius — minimal distance between samples
    • ◼
  • randomness — random angle on which the ring of test points turns, 0 < randomness < π
    • Output:
    Progress monitor is displayed during operation.
    ◼
  • output — array of samples with Point[] Head.
    • This is a small hack that allows you to use SparseArray to store the coordinates of points. Output can be used directly in graphics, or apply First to obtain Cartesian coordinates.
    In[]:=
    ClearAll[poissonDisk,insertPoint,isValidPoint];​​poissonDisk[numTries_,region_,radius_,randomness_]:=​​Module​​ κ=
    2.
    ,​​ bounds=RegionBounds@region,​​ numPoints,sizes,​​ origin,widthCells,heightCells,​​ init,grid,active, res,output,​​ position,insertPoint,isValidPoint,​​ regMember=RegionMember@region,​​ ​​ position[p_]:=Reverse@κ(p-origin)radius;​​ ​​ insertPoint[p_]:=​​ (grid〚Sequence@@position@p〛=Point@p);​​ ​​ isValidPoint[p_]:=​​ With[{pos=position@p},​​ regMember@p&&​​ AllTrue[​​ Take[grid,​​ Clip[First@pos+{-2,2},{1,heightCells}],​​ Clip[Last@pos+{-2,2},{1,widthCells}]]["ExplicitValues"],​​ EuclideanDistance[First@#,p]>radius&]​​ ];​​ ​​ origin=First/@bounds;​​ {widthCells, heightCells}=κ Abs/@(sizes=Subtract@@@bounds)radius+1;​​ numPoints = Area@region*widthCells*heightCellsTimes@@sizesπ;​​ init=RandomPoint[region];​​ grid=SparseArray[{},{heightCells,widthCells}];​​ active={init};​​ insertPoint@init;​​ ​​ Monitor​​ Whileactive!={},​​ active=Flatten​​ Reap[Do[ ​​ res=SelectFirst[​​ annulusMesh[elemActive,radius,2 radius,​​ randomness,numTries],​​ isValidPoint,{}]; ​​ If[res!={},​​ Sow@elemActive;Sow@res;​​ insertPoint@res​​ ]​​ ,{elemActive,active}]​​ ]2, 1;​​ (*End main While*)​​ output=grid["ExplicitValues"]​​ ,​​ ProgressIndicator[​​ Length@output,​​ {1,numPoints}];​​ output​​;

    Examples

    In[]:=
    regions=​​ (*Rectangle*)​​ Rectangle[{0.,0.},{1.,1.}],​​ (*Annulus*)​​ Annulus[0.5 {1.,1.},{0.15, 0.5}],​​ (*Star*)​​ BoundaryDiscretizeRegion@WindingPolygon​​ CirclePoints[.5{1.,1.},.5 ,5]{1,3,5,2,4},​​ (*Self-intersected*)​​ BoundaryDiscretizeRegion@ParametricRegion[​​ {.5{1.,1.}+.5 {s,s t},s^2+t^2≤1},{s,t}]​​;

    General usage

    In[]:=
    data1=poissonDisk[20,regions[[1]],0.06,1.5];
    In[]:=
    Graphics[{​​FaceForm[{Yellow,Opacity@0.3}],EdgeForm[Thickness@0.01],​​regions[[1]],​​FaceForm[{LightGray,Opacity@0.5}],​​EdgeForm[{Thickness@0.005,Gray}],​​Disk[First@#,0.03]&/@data1,​​data1},​​PlotRange->All,ImageSize->200]

    Different regions

    Randomness effect

    With small (~0.01—0.1) value of randomness resulted structure looks like polycrystalline, while for bigger values we get really uniform distribution.

    CITE THIS NOTEBOOK

    Poisson disk sampling​
    by Denis Ivanov​
    Wolfram Community, STAFF PICKS, May 26, 2024
    ​https://community.wolfram.com/groups/-/m/t/3182613