From b1fa56ddad12f9499849e9ec4c13a29e5497c64e Mon Sep 17 00:00:00 2001
From: Simon Forman <sforman@hushmail.com>
Date: Thu, 11 Apr 2024 21:05:10 -0700
Subject: [PATCH] Poisson process for 2D random points.

---
 poisson.py | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 poisson.py

diff --git a/poisson.py b/poisson.py
new file mode 100644
index 0000000..3582afe
--- /dev/null
+++ b/poisson.py
@@ -0,0 +1,90 @@
+'''
+Poisson Distribution
+'''
+from math import pi, sqrt, sin, cos, hypot
+from random import random, randint
+
+
+tau = 2 * pi
+
+
+def poisson(w, h, d, test_points=30):
+    '''
+    Return a list of coordinates generated by a Poisson point process.
+
+    Return a list of two-tuples, pairs of ints, representing points in a 2D
+    plane.  The points are bounded by the width and height (w and h) and
+    separated by a minimum distance of d.
+    '''
+    grid, x, y = Grid(d), randint(0, w), randint(0, h)
+    grid.add(x, y)
+    todo = [(x, y)]
+    while todo:
+        x, y = todo.pop(0)
+        todo.extend(
+            (nx, ny)
+            for nx, ny in candidates(x, y, d, test_points)
+            if 0 <= nx <= w and 0 <= ny <= h and grid.add(nx, ny)
+            )
+    return list(grid.grid.values())
+
+
+class Grid(object):
+    '''
+    Keep track of points in a plane and efficiently check that added points
+    are not nearer than a certain distance d to any others.
+    '''
+
+    def __init__(self, d):
+        self.d = d
+        self.r = d / sqrt(2) # A square of side r and diagonal d.
+        self.grid = {}  # Map from grid "bins" to points, 1-to-1.
+
+    def add(self, x, y):
+        '''
+        Add a point pair if it's not too close to any others and return a
+        bool: True if the point was added, False if it was rejected.
+        '''
+        grid_coords = round(x / self.r), round(y / self.r)
+        if any(hypot(x - nx, y - ny) < self.d
+                     for nx, ny in self._neighborhood(grid_coords)):
+            return False
+        self.grid[grid_coords] = x, y
+        return True
+
+    def _neighborhood(self, coordinates):
+        '''Return a list of the points in the neighborhood of (x, y).'''
+        # Note, x and y here are grid coordinates, not plane coordinates.
+        # This method returns plane coordinates.
+        x, y = coordinates
+        return filter(None, map(self.grid.get, [
+                                    (x-1, y-2), (x, y-2), (x+1, y-2),
+            (x-2, y-1), (x-1, y-1), (x, y-1), (x+1, y-1), (x+2, y-1),
+            (x-2, y),   (x-1, y),   (x, y),   (x+1, y),   (x+2, y),
+            (x-2, y+1), (x-1, y+1), (x, y+1), (x+1, y+1), (x+2, y+1),
+                                    (x-1, y+2), (x, y+2), (x+1, y+2),
+            ]))
+
+
+def candidates(x, y, d, test_points):
+    '''Generate points at a distance from (x, y) from d to 2*d.'''
+    for _ in range(test_points):
+        radius = d + random() * d
+        angle = random() * tau
+        yield (
+            int(round(x + radius * cos(angle))),
+            int(round(y + radius * sin(angle))),
+            )
+
+
+if __name__ == '__main__':
+    from random import seed
+
+    seed(23)
+
+    WIDTH, HEIGHT, DISTANCE = 80, 24, 4
+
+    P = set(poisson(WIDTH, HEIGHT, DISTANCE))
+
+    for y in range(HEIGHT + 1):
+        print(''.join(' *'[(x, y) in P] for x in range(WIDTH + 1)))