New pdPoissonDisc class: a fast(er) mechanism for rapid supersampling

Thank you to @miorsoft for first recommending Poisson Disc sampling to me:

miorsoft/Site#4 (comment)

The original paper on Bridson's algorithm is just one page long, and easily understandable:

https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

I wish more algorithms were this simple.

I hope to use this class to accelerate some effects and adjustments that currently experience egregious slowdowns when large radii are used.  Poisson disc sampling is cool for a lot of reasons, but one of the ones I'm most excited about is that you can easily mask the sample area to restrict the list of sampled points to an arbitrary geometric shape.  PD could use this for something like a fast bokeh effect (hypothetically... let's see how it pans out in practice).

Note that supersampling like this doesn't necessarily free a function from perf constraints tied to radius; instead, it greatly reduces the slope of the penalty associated with large r values.  Specific effects (e.g. gaussian blur(s)) will still benefit from application-specific implementations.

Classes/pdPoissonDisc.cls

@@ -0,0 +1,217 @@
+VERSION 1.0 CLASS
+BEGIN
+  MultiUse = -1  'True
+  Persistable = 0  'NotPersistable
+  DataBindingBehavior = 0  'vbNone
+  DataSourceBehavior  = 0  'vbNone
+  MTSTransactionMode  = 0  'NotAnMTSObject
+END
+Attribute VB_Name = "pdPoissonDisc"
+Attribute VB_GlobalNameSpace = False
+Attribute VB_Creatable = True
+Attribute VB_PredeclaredId = False
+Attribute VB_Exposed = False
+'***************************************************************************
+'PhotoDemon Poisson Disc Sampler
+'Copyright 2019-2019 by Tanner Helland
+'Created: 14/November/19
+'Last updated: 14/November/19
+'Last update: initial build
+'Dependencies: pdRandomize (for random number generation)
+'
+'Poisson disc sampling (PDS) is a fast way to supersample a region of 2D space:
+' https://en.wikipedia.org/wiki/Supersampling#Poisson_disc
+'
+'This class uses a fast generative technique known as Bridson's algorithm:
+' https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
+'
+'What's the relevance of this to a photo editor?  On some adjustments and effects,
+' PD uses supersampling to accelerate the pixel sampling process.  Some classes of
+' filters work just fine with a representative (but not comprehensive) sample of pixels,
+' and this class produces a supersample mapping guaranteed to produce a representative
+' sample for any arbitrary radius (r).
+'
+'An easy example of is the bilateral filter, a la the 2012 paper "A Low-Memory,
+' Straightforward and Fast Bilateral Filter Through Subsampling in Spatial Domain":
+' http://vcg.isti.cnr.it/Publications/2012/BCCS12/j.1467-8659.2011.02078.x.pdf
+'
+'Those authors compared PNSR between true bilateral results and supersampled results,
+' and PDS was found to be the highest-quality sampler - hence why I use it here.
+' Specifically, they said: "Note that the use of Poisson-disk sampling produces the
+' closest approximation to the full bilateral filter. The other sampling strategies create
+' visual artefacts and pattern-like artefacts, particularly when a regular sampling pattern
+' is used."
+'
+'Distances are currently calculated using standard cartesian distance.  In the future,
+' it may be interesting to add manhattan and chebyshev to the mix.
+'
+'All source code in this file is licensed under a modified BSD license.  This means you may use the code in your own
+' projects IF you provide attribution.  For more information, please visit https://photodemon.org/license/
+'
+'***************************************************************************
+
+Option Explicit
+
+'Use our own internal randomizer for point distribution
+Private m_Randomize As pdRandomize
+
+'Produce a list of points (indexed using the passed grid), evenly sampled at min radius r,
+' for a source surface of width WxH.
+'
+'Returns: FALSE if the requested radius is too small to produce a usable grid
+Friend Function GetDisc(ByRef dstPoints() As PointFloat, ByRef dstNumPoints As Long, ByRef outGrid() As Long, ByRef outGridWidth As Long, ByRef outGridHeight As Long, ByVal radius As Double, ByVal srcWidth As Long, ByVal srcHeight As Long) As Boolean
+
+    GetDisc = True
+
+    'Initialize the destination point list to a reasonable default
+    dstNumPoints = 64
+    ReDim dstPoints(0 To dstNumPoints - 1) As PointFloat
+
+    'We also need an internal set of "active" points (points for which we are still trying to assign
+    ' new neighbors).  This is sorta like a stack.
+    Dim activePoints() As PointFloat
+    ReDim activePoints(0 To dstNumPoints - 1) As PointFloat
+    Dim numActivePoints As Long
+
+    'Next, we need to set up the grid.  Per Bridson's algorithm, the size of each cell should be
+    ' r / sqrt(2).
+    Dim cellSize As Long
+    cellSize = Int(radius / Sqr(2#))
+
+    'If cells are too small, the caller needs to brute-force the algorithm instead
+    If (cellSize < 1) Then
+        GetDisc = False
+        Exit Function
+    End If
+
+    Dim rSquared As Double
+    rSquared = radius * radius
+
+    'Use the cell size and the source width/height to calculate width/height for the lookup grid.
+    ' (Note that we allocate an extra cell in either direction to improve boundary behavior;
+    ' in the future, something similar could be done in the left/top direction as well, though this
+    ' would add a minor perf hit for array lookups.)
+    outGridWidth = Int(srcWidth / cellSize + 0.9999999) + 1
+    outGridHeight = Int(srcHeight / cellSize + 0.9999999) + 1
+    ReDim outGrid(0 To outGridWidth - 1, 0 To outGridHeight - 1) As Long
+
+    'Initialize the grid to some kind of "null" flag.  (Cells can be empty.)
+    Dim x As Long, y As Long
+    For y = 0 To outGridHeight - 1
+    For x = 0 To outGridWidth - 1
+        outGrid(x, y) = -1
+    Next x
+    Next y
+
+    'Add an initial random point to both the the active collection.
+    ' (Note that we deliberately initialize it to be somewhere in the middle quadrant of
+    ' the grid; this improves performance over an initial point added to the periphery,
+    ' where points are more likely to fail multiple boundary checks on each addition.)
+    Dim tmpPoint As PointFloat
+    tmpPoint.x = (outGridWidth \ 4) + (m_Randomize.GetRandomFloat_WH() * (outGridWidth \ 2))
+    tmpPoint.y = (outGridHeight \ 4) + (m_Randomize.GetRandomFloat_WH() * (outGridHeight \ 2))
+
+    dstPoints(0) = tmpPoint
+    dstNumPoints = 1
+
+    activePoints(0) = tmpPoint
+    numActivePoints = 1
+
+    Dim apIndex As Long, ptAdded As Boolean
+    Dim rndTheta As Double, rndRadius As Double
+    Dim xMin As Long, xMax As Long, yMin As Long, yMax As Long, xGrd As Long, yGrd As Long
+
+    'Now we loop endlessly; as long as there are active points to process, we'll try and
+    ' add more points to the collection
+    Do While (numActivePoints > 0)
+
+        'Select a new active point
+        apIndex = m_Randomize.GetRandomIntRange_WH(0, numActivePoints - 1)
+
+        'We are now going to attempt to add new points around the current one.  If we succeed,
+        ' this point gets to stay in the active list, as do any point(s) we add.
+
+        'In his original paper (https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf)
+        ' Bridson suggests - without evidence, for better or worse - 30 sampling attempts per pixel.
+        ' You can obviously try fewer and risk a poorer distribution, at some boost to perf.
+        Const MAX_SAMPLING_ATTEMPTS As Long = 30
+        ptAdded = False
+
+        Dim k As Long
+        For k = 1 To MAX_SAMPLING_ATTEMPTS
+
+            'Generate a random theta and radius for this test point
+            rndTheta = m_Randomize.GetRandomFloat_WH() * PI_DOUBLE
+            rndRadius = radius + m_Randomize.GetRandomFloat_WH() * radius
+
+            'Construct a matching cartesian position
+            tmpPoint.x = activePoints(apIndex).x + rndRadius * Cos(rndTheta)
+            tmpPoint.y = activePoints(apIndex).y + rndRadius * Sin(rndTheta)
+
+            'Test this point; if it's good, add it and exit the loop
+
+            'First, test it against source image dimensions
+            If (tmpPoint.x < 0) Or (tmpPoint.y < 0) Then GoTo TryAnotherPoint
+            If (tmpPoint.x > srcWidth) Or (tmpPoint.y > srcHeight) Then GoTo TryAnotherPoint
+
+            'Next, figure out where this points (hypothetical) grid index lies
+            xGrd = Int(tmpPoint.x / cellSize)
+            yGrd = Int(tmpPoint.y / cellSize)
+
+            'Calculate corresponding loop intervals
+            xMin = xGrd - 1
+            If (xMin < 0) Then xMin = 0
+            yMin = yGrd - 1
+            If (yMin < 0) Then yMin = 0
+            xMax = xMin + 2
+            If (xMax >= outGridWidth) Then xMax = outGridWidth - 1
+            yMax = yMin + 2
+            If (yMax >= outGridHeight) Then yMax = outGridHeight - 1
+
+            For y = yMin To yMax
+            For x = xMin To xMax
+
+                'Only test non-null points
+                If (outGrid(x, y) >= 0) Then
+                    If (PDMath.DistanceTwoPointsShortcut(tmpPoint.x, tmpPoint.y, dstPoints(outGrid(x, y)).x, dstPoints(outGrid(x, y)).y) < rSquared) Then GoTo TryAnotherPoint
+                End If
+
+            Next x
+            Next y
+
+            'If we're still here, this point is valid!  Add it to both lists and mark its
+            ' corresponding grid index.
+            outGrid(xGrd, yGrd) = dstNumPoints
+
+            If (dstNumPoints > UBound(dstPoints)) Then ReDim Preserve dstPoints(0 To dstNumPoints * 2 - 1) As PointFloat
+            dstPoints(dstNumPoints) = tmpPoint
+            dstNumPoints = dstNumPoints + 1
+
+            If (numActivePoints > UBound(activePoints)) Then ReDim Preserve activePoints(0 To numActivePoints * 2 - 1) As PointFloat
+            activePoints(numActivePoints) = tmpPoint
+            numActivePoints = numActivePoints + 1
+
+            ptAdded = True
+            Exit For
+
+TryAnotherPoint:
+        Next k
+
+        'If a point was *not* added, this point should no longer be active.  Remove it from the list.
+        If (Not ptAdded) Then
+
+            'Fast remove; swap with the trailing point
+            activePoints(apIndex) = activePoints(numActivePoints - 1)
+            numActivePoints = numActivePoints - 1
+
+        End If
+
+    Loop
+
+End Function
+
+'The point of this class is to return
+Private Sub Class_Initialize()
+    Set m_Randomize = New pdRandomize
+    m_Randomize.SetSeed_AutomaticAndRandom
+End Sub

Classes/pdRandomize.cls

@@ -154,6 +154,10 @@ Friend Function GetRandomInt_VB() As Long
     GetRandomInt_VB = Int((m_HighBound - m_LowBound + 1) * Rnd + m_LowBound)
 End Function
 
+Friend Function GetRandomIntRange_VB(ByVal lowBound As Long, ByVal highBound As Long) As Long
+    GetRandomIntRange_VB = Int((highBound - lowBound + 1) * Rnd + lowBound)
+End Function
+
 'Return a random float using VB's internal randomize engine.  Bounds are ignored.  This is kind of a stupid function, as it would
 ' be faster to just use Rnd yourself, but it's included here for completeness.
 Friend Function GetRandomFloat_VB() As Double
@@ -162,7 +166,11 @@ End Function
 
 'Return a random integer using the Wichmann-Hill PRNG.  If supplied earlier, bounds are used.
 Friend Function GetRandomInt_WH() As Long
-    GetRandomInt_WH = Int((m_HighBound - m_LowBound + 1) * GetRandomFloat_WH + m_LowBound)
+    GetRandomInt_WH = Int((m_HighBound - m_LowBound + 1) * GetRandomFloat_WH() + m_LowBound)
+End Function
+
+Friend Function GetRandomIntRange_WH(ByVal lowBound As Long, ByVal highBound As Long) As Long
+    GetRandomIntRange_WH = Int((highBound - lowBound + 1) * GetRandomFloat_WH() + lowBound)
 End Function
 
 'Return a random float using the Wichmann-Hill PRNG.  Pretty fast, good distribution too.

Forms/MainWindow.frm

@@ -1891,6 +1891,31 @@ Private Sub MnuTest_Click()
 
     On Error GoTo StopTestImmediately
 
+    Dim cPoisson As pdPoissonDisc
+    Set cPoisson = New pdPoissonDisc
+
+    Dim listPoints() As PointFloat, numPoints As Long
+    Dim ptGrid() As Long, gridWidth As Long, gridHeight As Long
+    cPoisson.GetDisc listPoints, numPoints, ptGrid, gridWidth, gridHeight, 20#, PDImages.GetActiveImage.GetActiveDIB.GetDIBWidth, PDImages.GetActiveImage.GetActiveDIB.GetDIBHeight
+
+    'Draw all points?
+    Dim cSurface As pd2DSurface, cBrush As pd2DBrush
+    Set cSurface = New pd2DSurface
+    cSurface.WrapSurfaceAroundPDDIB PDImages.GetActiveImage.GetActiveDIB
+
+    Set cBrush = New pd2DBrush
+    cBrush.SetBrushColor RGB(255, 0, 0)
+
+    Dim i As Long
+    For i = 0 To numPoints - 1
+        PD2D.FillCircleF cSurface, cBrush, listPoints(i).x, listPoints(i).y, 3!
+    Next i
+
+    Set cSurface = Nothing
+    PDImages.GetActiveImage.NotifyImageChanged UNDO_Layer, PDImages.GetActiveImage.GetActiveLayerIndex
+
+    ViewportEngine.Stage2_CompositeAllLayers PDImages.GetActiveImage, FormMain.MainCanvas(0)
+
     'Filters_Scientific.InternalFFTTest
 
     'Want to test a new dialog?  Call it here, using a line like the following: