"""Poisson disc sampling layout using the algorithm of Bridson (2007).
Reference: https://doi.org/10.1145/1278780.1278807
"""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
import shapely
from shapely.affinity import scale as shp_scale
from shapely.affinity import translate as shp_translate
from .._random import farthest_point_sample
from .._validation import _validate_integer
from .pattern import _PatternLayout
if TYPE_CHECKING: # pragma: no cover
from .fitters import RigidPatternFitter
__all__ = ["PoissonDisc"]
#: Bridson's method at minimum separation ``r`` places about
#: ``_AREA_DENSITY * area / r**2 + _PERIMETER_DENSITY * perimeter / r`` points,
#: a least-squares fit over circles, rectangles, annuli and organic shapes.
_AREA_DENSITY = 0.63
_PERIMETER_DENSITY = 0.135
#: Points sampled per turbine kept. Few spare points leave the closest pair
#: of the farthest-point selection wider and the spacing more even.
_OVERSAMPLING = 1.05
#: A pattern of at most ``_FEW`` turbines is selected from about
#: ``_FEW_OVERSAMPLING`` times as many points, which widens its closest pair.
_FEW = 4
_FEW_OVERSAMPLING = 3.0
#: Largest and smallest factor by which one retry shrinks the separation.
_MAX_SHRINK = 0.98
_MIN_SHRINK = 0.5
_MINIMUM_DISTANCE_FLOOR = 1e-3
#: Grid cells of padding around the bounding box. A candidate lies less than
#: three cells from an accepted point, which lies inside the box, and its
#: 5x5 neighborhood reaches two cells further.
_PADDING = 5
#: Candidate offsets drawn per refill of the random buffer.
_BUFFER_ROWS = 4096
_MAX_K = 10_000
#: Most sampling-grid cells a fill allocates, which bounds its memory to
#: under a hundred megabytes and its count to about two million turbines.
_MAX_FILL_CELLS = 4_000_000
def _normalized_geometry(boundary):
"""``boundary`` geometry moved and scaled into the unit square."""
minx, miny, maxx, maxy = boundary.bounds
factor = 1.0 / max(maxx - minx, maxy - miny)
geometry = shp_translate(boundary.geometry, -minx, -miny)
return shp_scale(geometry, factor, factor, origin=(0, 0))
def _grid_shape(bounds, min_dist: float) -> tuple[int, int]:
"""Rows and columns of the padded sampling grid over ``bounds``."""
minx, miny, maxx, maxy = bounds
cell = min_dist / np.sqrt(2)
rows = int(np.ceil((maxy - miny) / cell)) + 1 + 2 * _PADDING
cols = int(np.ceil((maxx - minx) / cell)) + 1 + 2 * _PADDING
return rows, cols
def _expected_count(area: float, perimeter: float, min_dist: float) -> float:
"""Points Bridson's method places at ``min_dist``, from the fitted densities."""
return (
_AREA_DENSITY * area / min_dist**2 + _PERIMETER_DENSITY * perimeter / min_dist
)
def _distance_for_count(area: float, perimeter: float, count: float) -> float:
"""The ``min_dist`` at which :func:`_expected_count` equals ``count``."""
# Positive root of count r**2 - P_d perimeter r - A_d area = 0, in the
# form without cancellation.
a = _AREA_DENSITY * area
b = _PERIMETER_DENSITY * perimeter
return 2.0 * a / (np.sqrt(b * b + 4.0 * a * count) - b)
def _cell_classes(geom, origin, cell: float, shape: tuple[int, int]) -> np.ndarray:
"""Per grid cell 1 when inside ``geom``, 0 when outside, -1 when undecided.
A cell the boundary may cross is undecided and needs a point test. Every
other cell lies wholly on one side, the side of its center.
"""
rows, cols = shape
xs = origin[0] + (np.arange(cols) + 0.5) * cell
ys = origin[1] + (np.arange(rows) + 0.5) * cell
classes = shapely.contains_xy(geom, xs[None, :], ys[:, None]).astype(np.int8)
# Vertices a quarter cell apart come within an eighth of a cell of every
# boundary point, so each cell the boundary crosses holds a vertex or
# borders a cell that does.
ring = shapely.segmentize(geom.boundary, cell / 4.0)
ix = np.floor((shapely.get_coordinates(ring) - origin) / cell).astype(np.intp)
crossed = np.zeros((rows + 2, cols + 2), dtype=bool)
crossed[ix[:, 1] + 1, ix[:, 0] + 1] = True
near = np.zeros(shape, dtype=bool)
for dy in range(3):
for dx in range(3):
near |= crossed[dy : dy + rows, dx : dx + cols]
classes[near] = -1
return classes
def _bridson(geom, min_dist: float, k: int, rng) -> np.ndarray:
"""Bridson Poisson-disc sampling inside ``geom``, batched ``k`` at a time.
Each step draws ``k`` candidates around a random active point and accepts
the first that is inside ``geom`` and ``min_dist`` from every point, as
Bridson's one-at-a-time loop does.
"""
minx, miny, _, _ = geom.bounds
cell = min_dist / np.sqrt(2)
origin = np.array([minx, miny]) - _PADDING * cell
shape = _grid_shape(geom.bounds, min_dist)
width = shape[1]
classes = _cell_classes(geom, origin, cell, shape).ravel()
# At most one point fits in each cell. An empty cell holds an infinite
# position, which no candidate comes near.
cell_x = np.full(shape[0] * width, np.inf)
cell_y = np.full(shape[0] * width, np.inf)
placed: list[int] = []
active: list[int] = []
neighbors = np.array(
[dy * width + dx for dy in range(-2, 3) for dx in range(-2, 3)],
dtype=np.intp,
)
strides = np.array([1, width], dtype=np.intp)
min_dist_sq = min_dist**2
def _cells_of(points: np.ndarray) -> np.ndarray:
# The padding keeps every point right of and above the origin, so
# truncation rounds down.
return ((points - origin) / cell).astype(np.intp) @ strides
def _nearest_sq(points: np.ndarray, flat: np.ndarray) -> np.ndarray:
"""Squared distance from each point to the nearest placed point."""
around = flat[:, None] + neighbors
dx = cell_x[around]
dx -= points[:, :1]
dy = cell_y[around]
dy -= points[:, 1:]
dx *= dx
dy *= dy
dx += dy
return dx.min(axis=1)
def _place(point: np.ndarray, flat: int) -> None:
cell_x[flat] = point[0]
cell_y[flat] = point[1]
placed.append(flat)
active.append(flat)
# Seed from a guaranteed interior point of every polygon, so each part
# fills. A single polygon gets one seed, as Bridson's method has.
for part in getattr(geom, "geoms", (geom,)):
interior = part.representative_point()
point = np.array([[interior.x, interior.y]], dtype=float)
flat = _cells_of(point)
if _nearest_sq(point, flat)[0] >= min_dist_sq:
_place(point[0], int(flat[0]))
# A refill holds at least one batch.
rows = max(_BUFFER_ROWS, k)
def _offsets() -> np.ndarray:
draws = rng.random((rows, 2))
radii = min_dist * (1.0 + draws[:, 0])
angles = (2.0 * np.pi) * draws[:, 1]
return np.column_stack([radii * np.cos(angles), radii * np.sin(angles)])
offsets = _offsets()
picks = rng.random(_BUFFER_ROWS)
used = 0
picked = 0
while active:
if picked == len(picks):
picks = rng.random(_BUFFER_ROWS)
picked = 0
ai = int(picks[picked] * len(active))
picked += 1
if used + k > len(offsets):
offsets = np.concatenate([offsets[used:], _offsets()])
used = 0
flat = active[ai]
candidates = offsets[used : used + k] + (cell_x[flat], cell_y[flat])
cells = _cells_of(candidates)
kind = classes[cells]
open_ = (_nearest_sq(candidates, cells) >= min_dist_sq) & (kind != 0)
accepted = int(open_.argmax()) if open_.any() else -1
if accepted >= 0 and kind[accepted] < 0:
# The first open candidate in a cell wholly inside is accepted
# unless an earlier open one in an undecided cell passes the
# point test.
inside = open_ & (kind > 0)
stop = int(inside.argmax()) if inside.any() else k
undecided = np.flatnonzero(open_[:stop])
hits = np.flatnonzero(
shapely.contains_xy(
geom, candidates[undecided, 0], candidates[undecided, 1]
)
)
accepted = (
int(undecided[hits[0]]) if len(hits) else (stop if stop < k else -1)
)
if accepted >= 0:
# Only the draws up to the accepted candidate are consumed, as
# Bridson's one-at-a-time loop would.
used += accepted + 1
_place(candidates[accepted], int(cells[accepted]))
else:
used += k
active[ai] = active[-1]
active.pop()
return np.column_stack([cell_x[placed], cell_y[placed]])
[docs]
class PoissonDisc(_PatternLayout):
"""Place turbines by Poisson-disc sampling, a blue-noise pattern.
Bridson's algorithm, with ``k`` candidates per active point, fills a copy
of the boundary scaled into the unit square. Farthest-point sampling keeps
exactly ``n_turbines`` points, and the pattern is then scaled and translated
to fill the boundary. With ``minimum_spacing``, Bridson's algorithm fills the
boundary itself with that minimum distance, starting in every polygon, and
the count is however many points fit.
The minimum separation starts where Bridson's method is expected to place
5% more points than ``n_turbines``, or three times as many for at most
four turbines, from the area and perimeter of the normalized boundary.
While too few fit, each retry rescales that estimate by the count the
last run placed, shrinking by a factor of 0.5 to 0.98. Few spare points
keep the closest pair of the selection wide. The search stops at 1e-3 for
a square bounding box and lower for an elongated one, by the square root
of the short side over the long side, so the sampling grid keeps the same
number of cells. The value used is reported as
``metadata["fit"]["minimum_distance"]``. Raises ``ValueError`` if too few
points fit even then.
"""
fit_method = "poisson_disc"
_fills_directly = True
def __init__(
self,
n_turbines: int | None = None,
*,
density: dict[str, float] | None = None,
k: int = 30,
fitter: RigidPatternFitter | None = None,
seed: int | np.random.Generator | None = None,
):
self.k = _validate_integer("k", k)
if self.k > _MAX_K:
raise ValueError(f"k must be a positive integer <= {_MAX_K}")
super().__init__(
n_turbines,
density=density,
fitter=fitter,
seed=seed,
)
def _fill_positions(self, boundary, spacing, rng):
geometry = boundary.geometry
rows, cols = _grid_shape(geometry.bounds, spacing)
if rows * cols > _MAX_FILL_CELLS:
raise ValueError(
f"a spacing of {spacing:.4g} m is too small for this boundary, "
f"which would hold about {int(2 * geometry.area / spacing**2):,} "
"turbines"
)
positions = _bridson(geometry, spacing, self.k, rng)
return positions, {
"method": self.fit_method,
"minimum_distance": spacing,
"k": self.k,
}
def _pattern(self, boundary, rng):
n = self._required_n
# Sample a normalized boundary, then scale and translate the result.
geometry = _normalized_geometry(boundary)
# An elongated farm lowers the floor, keeping the grid's cell count.
minx, miny, maxx, maxy = geometry.bounds
floor = _MINIMUM_DISTANCE_FLOOR * np.sqrt(
min(maxx - minx, maxy - miny) / max(maxx - minx, maxy - miny)
)
area = geometry.area
perimeter = geometry.length
target = (_FEW_OVERSAMPLING if n <= _FEW else _OVERSAMPLING) * n
min_dist = max(_distance_for_count(area, perimeter, target), floor)
pts = _bridson(geometry, min_dist, self.k, rng)
while len(pts) < n and min_dist > floor:
# The fitted count, scaled to what this boundary placed, gives
# the separation at which the target fits.
placed = len(pts) / _expected_count(area, perimeter, min_dist)
shrink = _distance_for_count(area, perimeter, target / placed) / min_dist
min_dist *= min(max(shrink, _MIN_SHRINK), _MAX_SHRINK)
pts = _bridson(geometry, min_dist, self.k, rng)
if len(pts) < n:
raise ValueError(
f"Poisson-disc sampling fitted only {len(pts)} of {n} turbines, "
"a larger k may fit more"
)
positions = farthest_point_sample(pts, n, rng)
# Report inverse normalization as a fit hint for rigid placement.
minx, miny, maxx, maxy = boundary.bounds
hint_scale = max(maxx - minx, maxy - miny)
hint_translation = positions.mean(axis=0) * hint_scale + np.array([minx, miny])
return (
positions,
{
"minimum_distance": float(min_dist),
"k": self.k,
"fit_hint": (hint_scale, hint_translation),
},
None,
)
