Het artikel bespreekt Poisson Disk Sampling, een techniek uit de computergraphics om punten willekeurig in een ruimte te plaatsen terwijl een minimale afstand ($r$) tussen de punten wordt gewaarborgd. Dit voorkomt ongewenste clustering die optreedt bij standaard random sampling.
De kern van het artikel richt zich op het algoritme van Bridson, dat gebruikmaakt van een raster voor efficiënte botsingsdetectie en een 'active list' om nieuwe punten te samplen uit een annulus (ringvorm) rond bestaande punten. De auteur beschrijft twee optimalisaties voor dit algoritme:
- Parentale optimalisatie (2D): Het uitsluiten van hoeken die te dicht bij het ouderpunt liggen.
- Afstandsoptimalisatie: Het aanpassen van de radius-distributie via een constante $c$ om een balans te vinden tussen dichtheid en willekeur.
Daarnaast wordt de toepassing van stippling genoemd, waarbij de afstand $r$ wordt gekoppeld aan de helderheid van pixels in een afbeelding. Tot slot wordt de methode van Scott A. Mitchell (2022) geïntroduceerd, die superieur is aan die van Bridson omdat deze deterministisch is en maximale uniformiteit garandeert zonder gebruik te maken van rejection sampling.
Over Poisson Disk Sampling
In 2024 publiceerde een team van negen wiskundigen een monsterlijk, bijna 1.000 pagina's tellend bewijs van het geometrische Langlands-vermoeden. Het is een kroonprestatie in de zuivere wiskunde, en ik heb geaccepteerd dat ik nooit zelfs maar de stellingen die zij bewezen zal begrijpen, laat staan het bewijs zelf.
Aan het tegenovergestelde uiteinde van het spectrum publiceerde Robert Bridson in 2007 een artikel van één pagina dat bijna 1.000 citaties heeft en minder dan 10 minuten kost om volledig te begrijpen. Het presenteert een eenvoudige oplossing voor een probleem dat veelvuldig voorkomt in computergraphics en simulaties: dingen willekeurig plaatsen, maar niet te dicht bij elkaar.
Stel dat je procedureel een bos probeert te genereren en een manier nodig hebt om de bomen te plaatsen. Het probleem met eenvoudige willekeurige sampling (random sampling) is duidelijk: sommige bomen zouden bovenop elkaar staan. Wat we nodig hebben is de mogelijkheid om een minimale afstand tussen twee willekeurige bomen in te stellen. Een distributie van bomen die aan deze regel voldoet, wordt een Poisson disk distribution genoemd.
We zouden een naïve rejection sampling-aanpak kunnen proberen, waarbij we willekeurige punten "werpen" en elk punt verwerpen dat binnen die minimale afstand van een ander punt valt. Maar zonder een slimme datastructuur kost het controleren van botsingen voor elk sample lineaire tijd, en het verwerpingpercentage nadert snel de honderd procent. Het algoritme van Bridson biedt ons een efficiënte manier om dit te doen.
Het algoritme van Bridson
Stel dat de gewenste minimale afstand tussen punten $r$ is en dat we werken in een $d$-dimensionale ruimte. Het algoritme van Bridson werkt als volgt:
- Verdeel de ruimte in een raster (grid) met een zijdelengte van $\frac{r}{\sqrt{d}}$. Dit garandeert dat elke rastercel maximaal één punt kan bevatten.
- Initialiseer een lijst
active met één willekeurig punt dat uniform uit de ruimte is gekozen.
- Zolang
active niet leeg is:
- Selecteer een element $p$ uit
active op basis van een uniforme willekeurige keuze.
- Sample uniform de annulus (ringvorm) gecentreerd rond $p$ met een binnenradius $r$ en een buitenradius $2r$, maximaal $k$ keer.
- Als er een geldig Poisson disk-sample wordt gevonden (met gebruik van het raster voor efficiënte botsingsdetectie), voeg dit dan toe aan
active en kies een nieuwe $p$.
- Als er binnen $k$ pogingen geen geldig punt wordt gevonden, verwijder dan $p$ uit
active. Bridson adviseert om $k=30$ in te stellen.
De eenvoudigste manier om de annulus uniform te samplen is door een willekeurige eenheidsvector $\vec{v} \in \mathbb{R}^d$ te genereren en een getal $x$ uniform te kiezen uit het interval $[1/2^d, 1)$, waarbij het uiteindelijke sample $2rx^{1/d} \cdot \vec{v}$ is. In twee dimensies is het kiezen van een eenheidsvector gelijk aan het kiezen van een hoek $\theta \in [0, 2\pi)$. In hogere dimensies kun je een vector normaliseren waarbij elke component wordt gesampled uit een normale distributie.
Verbeteringen
Er zijn twee eenvoudige verbeteringen aan het algoritme van Bridson die het aantal benodigde iteraties om hetzelfde aantal punten te genereren drastisch verminderen. De eerste werkt in twee dimensies, de tweede werkt ook in hogere dimensies.
Parentale optimalisatie (2D)
Bij de tweedimensionale verbetering kijken we naar het moment waarop het algoritme een punt $p$ plaatst en vervolgens zijn annulus sampled om een nieuw punt $q$ te krijgen. We noemen $p$ de ouder van $q$. Er is waardevolle informatie opgeslagen in de relatie tussen deze punten. Wanneer we onvermijdelijk de annulus rond $q$ samplen, is er een heel bereik aan hoeken die we niet hoeven te overwegen, omdat de punten daarbinnen te dicht bij $p$ zouden liggen.
Hoewel de visuele intuïtie eenvoudig is, is het vertalen hiervan naar een formule een tedious trigonometrie-oefening. De conus gevormd door de uitsluitingszone is gecentreerd op hoek $\alpha$ en heeft een breedte van $2\beta$, waarbij:
$$\alpha = \operatorname{atan2}(py - qy, px - qx)$$ $$\beta = \min\left(\arccos\frac{|p-q|^2+3r^2}{4r \cdot |p-q|}, \arccos\frac{|p-q|}{2r}\right)$$
Het interessante deel van deze formule is het minimum in de vergelijking voor $\beta$. Dit houdt rekening met het feit dat ofwel de binnenste of de buitenste cirkel van de annulus de conus kan begrenzen, afhankelijk van de afstand tussen $p$ en $q$. We moeten het minimum kiezen om te garanderen dat de intersectie van de conus en de annulus volledig in de cirkel is bevat. De grenspunten van de conus springen van de buitenste cirkel naar de binnenste wanneer de afstand tussen de punten $3 \cdot r$ overschrijdt.
De implementatie hiervan vereist enkel het opslaan van de ouder van elk punt. Vervolgens kun je de hoeken van de conus berekenen en de hoek $\theta$ voor het volgende sample kiezen in het bereik buiten de conus.
Deze verbetering zou waarschijnlijk kunnen worden gegeneraliseerd naar hogere dimensies, maar dat zou meer ruimte vereisen om de contactvectoren van de annuli op te slaan. Bovendien zouden de voordelen waarschijnlijk afnemen, omdat het volume van de intersectie van een annulus met een sfeer proportioneel insignificant wordt in hogere dimensies.
Optimalisatie van de afstand
De tweede verbetering verandert hoe we de afstand tot het volgende sample kiezen, in plaats van de hoek. Beschouw de distributie van de afstanden van elk punt in de annulus tot het centrum. De cumulatieve distributiefunctie (CDF) is proportioneel aan $x^d$ op het interval $[r, 2r]$.
Wat gebeurt er als we de exponent veranderen naar een constante $c$ anders dan $d$? Dan kunnen we de punten dichter bij of verder weg van het centrum plaatsen. Voor $c \neq 0$ is de exacte CDF:
$$F_c(x) = \begin{cases} 0 & x \leq r, \\ \frac{x^c - r^c}{(2^c - 1)r^c} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$
Voor $c = 0$ is $F0(x)$ ongedefinieerd vanwege een deling door nul. Om dit op te lossen definiëren we $F0(x) = \lim{c \to 0} Fc(x)$, wat resulteert in:
$$F0(x) = \begin{cases} 0 & x \leq r, \\ \log2{x} - \log_2{r} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$
Om de radius te samplen met een willekeurige waarde van $c$, kunnen we inverse transform sampling toepassen. Wanneer $c = 0$ is, moet de radius $r \cdot 2^x$ zijn, waarbij $x$ een uniforme willekeurige variabele is op het interval $[0, 1)$. Anders gebruiken we $2ry^{1/c}$, waarbij $y$ een uniforme willekeurige variabele is tussen $1$ en $1/2^c$.
Het lijkt erop dat we $c$ op een zeer negatief getal moeten instellen, of zelfs de limiet moeten nemen waarbij $c$ naar negatief oneindig gaat, waardoor elk punt precies op afstand $r$ van zijn ouder komt te staan. Hoewel dit het aantal gegenereerde punten maximaliseert, gaat dit ten koste van het "willekeurige" gevoel van de distributie. In het extreme geval van $c = -\infty$ ontstaan er artefacten zoals lange reeksen punten en gaten waar de beperkte afstand niet kan reiken.
Er moet dus een balans worden gevonden tussen het maximaliseren van de dichtheid en het behouden van willekeur. Voor een $k$ tussen 15 en 40 is het optimaal om $c = -1.4 - \frac{17}{\sqrt{k}}$ in te stellen bij gebruik van de parentale optimalisatie. Deze vergelijking is empirisch afgeleid om het aantal gegenereerde punten ongeveer te laten overeenkomen met de verwachte output van een uniforme en maximale Poisson disk-sampler. Specifiek is gezocht naar een waarde van $c$ waarbij cirkels met radius $r/2$ rond elk punt samen 54,7% van het totale oppervlak beslaan (de verzadigde dekking van cirkelvormige schijven onder het random sequential adsorption-model). Dit geldt enkel voor twee dimensies.
Stippling
Tot nu toe hebben we $r$ constant gehouden, maar dat is niet noodzakelijk. We kunnen de minimale afstand tussen punten dynamisch instellen volgens een functie $r: \mathbb{R}^d \to \mathbb{R}$. Na het plaatsen van een punt $p$ samplen we een annulus met binnenradius $r(p)$.
Een leuke toepassing hiervan is om $r$ te definiëren als de helderheid van elke pixel in een afbeelding om een stippling-effect te produceren.
Het algoritme van Bridson is inherent sequentieel, maar er zijn andere algoritmen die parallel kunnen worden uitgevoerd voor enorme prestatieverbeteringen. Een voorbeeld is PixelPie, dat volledig op de GPU draait.
Maximaliteit en Uniformiteit
In 2022 publiceerde Scott A. Mitchell een methode voor het genereren van Poisson disk-samples in twee dimensies. Zijn algoritme is beter dan dat van Bridson op drie punten:
- Maximaliteit: Na afronding is het gegarandeerd onmogelijk om nog een punt toe te voegen zonder de Poisson disk-eigenschap te schenden.
- Uniformiteit: Het algoritme sampled een uniforme distributie over alle maximale sets van Poisson disk-samples.
- Determinisme: Het algoritme vertrouwt niet op rejection sampling, wat betekent dat er geen mislukte pogingen zijn om punten te plaatsen.
Hoewel maximaliteit en uniformiteit ook door andere algoritmen zijn bereikt (zoals hierarchical dart throwing), is Mitchells methode de eerste die dit doet zonder enige vorm van rejection sampling. De implementatie is ongeveer net zo snel en genereert ongeveer evenveel punten als een geoptimaliseerde versie van Bridson (met $k=20$ en $c=-5.2$), maar is aanzienlijk complexer.
Een sterk vereenvoudigde samenvatting van de methode:
- Verdeel de ruimte in een raster met zijdelengte $\frac{r}{\sqrt{d}}$.
- Zolang er ruimte is om een nieuw punt te plaatsen:
- Selecteer willekeurig een cel $c$ uit het raster, gewogen naar de resterende oppervlakten.
- Ontbind $c$ in disjuncte driehoeken en chocks [^1].
- Selecteer een willekeurige driehoek of chock $t$, gewogen naar oppervlakte.
- Sample uniform een punt $p$ uit $t$ en voeg dit toe aan de uiteindelijke set.
- Snijd een cirkel met radius $r$ gecentreerd op $p$ uit het raster.
***
[^1]: Een chock is een driedelige vorm begrensd door een cirkel, een radiale straal en een raaklijn.
Over Poisson Disk Sampling
In 2024 publiceerde een team van negen wiskundigen een monsterlijk, bijna 1.000 pagina's tellend bewijs van het geometrische Langlands-vermoeden. Het is een kroonprestatie in de zuivere wiskunde, en ik heb geaccepteerd dat ik nooit zelfs maar de stellingen die zij bewezen zal begrijpen, laat staan het bewijs zelf.
Aan het tegenovergestelde uiteinde van het spectrum publiceerde Robert Bridson in 2007 een artikel van één pagina dat bijna 1.000 citaties heeft en minder dan 10 minuten kost om volledig te begrijpen. Het presenteert een eenvoudige oplossing voor een probleem dat veelvuldig voorkomt in computergraphics en simulaties: dingen willekeurig plaatsen, maar niet te dicht bij elkaar.
Stel dat je procedureel een bos probeert te genereren en een manier nodig hebt om de bomen te plaatsen. Het probleem met eenvoudige willekeurige sampling (random sampling) is duidelijk: sommige bomen zouden bovenop elkaar staan. Wat we nodig hebben is de mogelijkheid om een minimale afstand tussen twee willekeurige bomen in te stellen. Een distributie van bomen die aan deze regel voldoet, wordt een Poisson disk distribution genoemd.
We zouden een naïve rejection sampling-aanpak kunnen proberen, waarbij we willekeurige punten "werpen" en elk punt verwerpen dat binnen die minimale afstand van een ander punt valt. Maar zonder een slimme datastructuur kost het controleren van botsingen voor elk sample lineaire tijd, en het verwerpingpercentage nadert snel de honderd procent. Het algoritme van Bridson biedt ons een efficiënte manier om dit te doen.
Het algoritme van Bridson
Stel dat de gewenste minimale afstand tussen punten $r$ is en dat we werken in een $d$-dimensionale ruimte. Het algoritme van Bridson werkt als volgt:
- Verdeel de ruimte in een raster (grid) met een zijdelengte van $\frac{r}{\sqrt{d}}$. Dit garandeert dat elke rastercel maximaal één punt kan bevatten.
- Initialiseer een lijst
active met één willekeurig punt dat uniform uit de ruimte is gekozen.
- Zolang
active niet leeg is:
- Selecteer een element $p$ uit
active op basis van een uniforme willekeurige keuze.
- Sample uniform de annulus (ringvorm) gecentreerd rond $p$ met een binnenradius $r$ en een buitenradius $2r$, maximaal $k$ keer.
- Als er een geldig Poisson disk-sample wordt gevonden (met gebruik van het raster voor efficiënte botsingsdetectie), voeg dit dan toe aan
active en kies een nieuwe $p$.
- Als er binnen $k$ pogingen geen geldig punt wordt gevonden, verwijder dan $p$ uit
active. Bridson adviseert om $k=30$ in te stellen.
De eenvoudigste manier om de annulus uniform te samplen is door een willekeurige eenheidsvector $\vec{v} \in \mathbb{R}^d$ te genereren en een getal $x$ uniform te kiezen uit het interval $[1/2^d, 1)$, waarbij het uiteindelijke sample $2rx^{1/d} \cdot \vec{v}$ is. In twee dimensies is het kiezen van een eenheidsvector gelijk aan het kiezen van een hoek $\theta \in [0, 2\pi)$. In hogere dimensies kun je een vector normaliseren waarbij elke component wordt gesampled uit een normale distributie.
Verbeteringen
Er zijn twee eenvoudige verbeteringen aan het algoritme van Bridson die het aantal benodigde iteraties om hetzelfde aantal punten te genereren drastisch verminderen. De eerste werkt in twee dimensies, de tweede werkt ook in hogere dimensies.
Parentale optimalisatie (2D)
Bij de tweedimensionale verbetering kijken we naar het moment waarop het algoritme een punt $p$ plaatst en vervolgens zijn annulus sampled om een nieuw punt $q$ te krijgen. We noemen $p$ de ouder van $q$. Er is waardevolle informatie opgeslagen in de relatie tussen deze punten. Wanneer we onvermijdelijk de annulus rond $q$ samplen, is er een heel bereik aan hoeken die we niet hoeven te overwegen, omdat de punten daarbinnen te dicht bij $p$ zouden liggen.
Hoewel de visuele intuïtie eenvoudig is, is het vertalen hiervan naar een formule een tedious trigonometrie-oefening. De conus gevormd door de uitsluitingszone is gecentreerd op hoek $\alpha$ en heeft een breedte van $2\beta$, waarbij:
$$\alpha = \operatorname{atan2}(py - qy, px - qx)$$ $$\beta = \min\left(\arccos\frac{|p-q|^2+3r^2}{4r \cdot |p-q|}, \arccos\frac{|p-q|}{2r}\right)$$
Het interessante deel van deze formule is het minimum in de vergelijking voor $\beta$. Dit houdt rekening met het feit dat ofwel de binnenste of de buitenste cirkel van de annulus de conus kan begrenzen, afhankelijk van de afstand tussen $p$ en $q$. We moeten het minimum kiezen om te garanderen dat de intersectie van de conus en de annulus volledig in de cirkel is bevat. De grenspunten van de conus springen van de buitenste cirkel naar de binnenste wanneer de afstand tussen de punten $3 \cdot r$ overschrijdt.
De implementatie hiervan vereist enkel het opslaan van de ouder van elk punt. Vervolgens kun je de hoeken van de conus berekenen en de hoek $\theta$ voor het volgende sample kiezen in het bereik buiten de conus.
Deze verbetering zou waarschijnlijk kunnen worden gegeneraliseerd naar hogere dimensies, maar dat zou meer ruimte vereisen om de contactvectoren van de annuli op te slaan. Bovendien zouden de voordelen waarschijnlijk afnemen, omdat het volume van de intersectie van een annulus met een sfeer proportioneel insignificant wordt in hogere dimensies.
Optimalisatie van de afstand
De tweede verbetering verandert hoe we de afstand tot het volgende sample kiezen, in plaats van de hoek. Beschouw de distributie van de afstanden van elk punt in de annulus tot het centrum. De cumulatieve distributiefunctie (CDF) is proportioneel aan $x^d$ op het interval $[r, 2r]$.
Wat gebeurt er als we de exponent veranderen naar een constante $c$ anders dan $d$? Dan kunnen we de punten dichter bij of verder weg van het centrum plaatsen. Voor $c \neq 0$ is de exacte CDF:
$$F_c(x) = \begin{cases} 0 & x \leq r, \\ \frac{x^c - r^c}{(2^c - 1)r^c} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$
Voor $c = 0$ is $F0(x)$ ongedefinieerd vanwege een deling door nul. Om dit op te lossen definiëren we $F0(x) = \lim{c \to 0} Fc(x)$, wat resulteert in:
$$F0(x) = \begin{cases} 0 & x \leq r, \\ \log2{x} - \log_2{r} & r < x \leq 2r, \\ 1 & x > 2r. \end{cases}$$
Om de radius te samplen met een willekeurige waarde van $c$, kunnen we inverse transform sampling toepassen. Wanneer $c = 0$ is, moet de radius $r \cdot 2^x$ zijn, waarbij $x$ een uniforme willekeurige variabele is op het interval $[0, 1)$. Anders gebruiken we $2ry^{1/c}$, waarbij $y$ een uniforme willekeurige variabele is tussen $1$ en $1/2^c$.
Het lijkt erop dat we $c$ op een zeer negatief getal moeten instellen, of zelfs de limiet moeten nemen waarbij $c$ naar negatief oneindig gaat, waardoor elk punt precies op afstand $r$ van zijn ouder komt te staan. Hoewel dit het aantal gegenereerde punten maximaliseert, gaat dit ten koste van het "willekeurige" gevoel van de distributie. In het extreme geval van $c = -\infty$ ontstaan er artefacten zoals lange reeksen punten en gaten waar de beperkte afstand niet kan reiken.
Er moet dus een balans worden gevonden tussen het maximaliseren van de dichtheid en het behouden van willekeur. Voor een $k$ tussen 15 en 40 is het optimaal om $c = -1.4 - \frac{17}{\sqrt{k}}$ in te stellen bij gebruik van de parentale optimalisatie. Deze vergelijking is empirisch afgeleid om het aantal gegenereerde punten ongeveer te laten overeenkomen met de verwachte output van een uniforme en maximale Poisson disk-sampler. Specifiek is gezocht naar een waarde van $c$ waarbij cirkels met radius $r/2$ rond elk punt samen 54,7% van het totale oppervlak beslaan (de verzadigde dekking van cirkelvormige schijven onder het random sequential adsorption-model). Dit geldt enkel voor twee dimensies.
Stippling
Tot nu toe hebben we $r$ constant gehouden, maar dat is niet noodzakelijk. We kunnen de minimale afstand tussen punten dynamisch instellen volgens een functie $r: \mathbb{R}^d \to \mathbb{R}$. Na het plaatsen van een punt $p$ samplen we een annulus met binnenradius $r(p)$.
Een leuke toepassing hiervan is om $r$ te definiëren als de helderheid van elke pixel in een afbeelding om een stippling-effect te produceren.
Het algoritme van Bridson is inherent sequentieel, maar er zijn andere algoritmen die parallel kunnen worden uitgevoerd voor enorme prestatieverbeteringen. Een voorbeeld is PixelPie, dat volledig op de GPU draait.
Maximaliteit en Uniformiteit
In 2022 publiceerde Scott A. Mitchell een methode voor het genereren van Poisson disk-samples in twee dimensies. Zijn algoritme is beter dan dat van Bridson op drie punten:
- Maximaliteit: Na afronding is het gegarandeerd onmogelijk om nog een punt toe te voegen zonder de Poisson disk-eigenschap te schenden.
- Uniformiteit: Het algoritme sampled een uniforme distributie over alle maximale sets van Poisson disk-samples.
- Determinisme: Het algoritme vertrouwt niet op rejection sampling, wat betekent dat er geen mislukte pogingen zijn om punten te plaatsen.
Hoewel maximaliteit en uniformiteit ook door andere algoritmen zijn bereikt (zoals hierarchical dart throwing), is Mitchells methode de eerste die dit doet zonder enige vorm van rejection sampling. De implementatie is ongeveer net zo snel en genereert ongeveer evenveel punten als een geoptimaliseerde versie van Bridson (met $k=20$ en $c=-5.2$), maar is aanzienlijk complexer.
Een sterk vereenvoudigde samenvatting van de methode:
- Verdeel de ruimte in een raster met zijdelengte $\frac{r}{\sqrt{d}}$.
- Zolang er ruimte is om een nieuw punt te plaatsen:
- Selecteer willekeurig een cel $c$ uit het raster, gewogen naar de resterende oppervlakten.
- Ontbind $c$ in disjuncte driehoeken en chocks [^1].
- Selecteer een willekeurige driehoek of chock $t$, gewogen naar oppervlakte.
- Sample uniform een punt $p$ uit $t$ en voeg dit toe aan de uiteindelijke set.
- Snijd een cirkel met radius $r$ gecentreerd op $p$ uit het raster.
***
[^1]: Een chock is een driedelige vorm begrensd door een cirkel, een radiale straal en een raaklijn.
Groetjes,