"""Strategies for maximally fitting layout geometry into a boundary."""
from __future__ import annotations
import weakref
from abc import ABC, abstractmethod
from functools import cache
from itertools import combinations
from typing import TYPE_CHECKING, Any
import numpy as np
import shapely
from .._random import farthest_point_sample
from .._repr import _ArgumentRepr
from .._validation import _validate_integer, _validate_positive
if TYPE_CHECKING: # pragma: no cover
from ..boundaries.base import Boundary
from .base import LayoutAlgorithm
from .candidate import _CandidateLayout
__all__ = [
"Fitter",
"RigidPatternFitter",
"SelectedCandidateFitter",
"SelectedCircularFitter",
"SelectedContourFitter",
"SelectedGridFitter",
]
#: Each scale the rigid fit scans is this fraction of the one before it.
_SCAN_RATIO = 0.9
#: Enough scan steps to reach about a thousandth of the upper bound.
_SCAN_STEPS = 64
#: Points of a rigid pattern ordered to spread out before the rest, so the
#: first few cut the region of fitting translations the most.
_SPREAD_POINTS = 32
#: The rigid fit tries the scale its region's edges predict this fraction of
#: the tolerance either side, so a right prediction closes the bracket.
_MEETING_MARGIN = 1e-3
#: Most edge lines of a region the rigid fit predicts a scale from.
_MEETING_LINES = 24
[docs]
class Fitter(_ArgumentRepr, ABC):
"""Base strategy for fitting one layout family to one polygonal farm."""
[docs]
@abstractmethod
def fit(
self,
layout: Any,
boundary: Boundary,
pattern: np.ndarray | None,
rng: np.random.Generator,
hint: tuple[float, np.ndarray] | None = None,
) -> tuple[np.ndarray, dict]:
"""Return fitted positions and JSON-compatible fit metadata.
Pattern layouts pass their ``(N, 2)`` point ``pattern`` and an optional
verified ``(scale, translation)`` ``hint``. Candidate layouts build
their own candidate field, so they pass ``pattern=None`` and no hint.
"""
[docs]
class RigidPatternFitter(Fitter):
"""Maximally scale and translate an exact finite point pattern.
The scale is searched for at most ``max_iterations`` steps, stopping when
the bracket is within ``tolerance`` relative to the upper bound. A step
tries the scale at which the lines of the fitting region's edges meet,
just either side of it, or bisects when that does not halve the bracket.
"""
def __init__(self, *, max_iterations: int = 50, tolerance: float = 1e-6):
self.max_iterations = _validate_integer("max_iterations", max_iterations)
self.tolerance = _validate_positive("tolerance", tolerance)
[docs]
def fit(
self,
layout: LayoutAlgorithm,
boundary: Boundary,
pattern: np.ndarray | None,
rng: np.random.Generator,
hint: tuple[float, np.ndarray] | None = None,
) -> tuple[np.ndarray, dict]:
"""Return the largest scaled and translated ``pattern`` inside ``boundary``.
``pattern`` has shape ``(N, 2)``. Positions come back in meters with fit
metadata. An optional ``(scale, translation)`` ``hint`` that is verified
feasible becomes a lower bound on the scale. A missing or empty
``pattern`` raises ``ValueError``.
"""
if pattern is None:
raise ValueError("RigidPatternFitter requires a point pattern")
pattern = np.asarray(pattern, dtype=float).reshape(-1, 2)
if not len(pattern):
raise ValueError("cannot fit an empty pattern")
geom = boundary.geometry
minx, miny, maxx, maxy = boundary.bounds
width, height = maxx - minx, maxy - miny
pattern_width = float(np.ptp(pattern[:, 0]))
pattern_height = float(np.ptp(pattern[:, 1]))
limits = []
if pattern_width > 0:
limits.append(width / pattern_width)
if pattern_height > 0:
limits.append(height / pattern_height)
upper = min(limits) if limits else max(width, height)
lower = 0.0
best_scale = 0.0
# Use an interior point for the defensive zero-scale fallback.
fallback_point = geom.representative_point()
best_translation = np.array([fallback_point.x, fallback_point.y])
# The region of fitting translations at ``best_scale`` once built, or
# ``None`` where rounding left it empty, the points whose shifted
# boundaries were cut to build it, and the scale it was built at.
best_region = best_used = None
region_scale = None
def fits(scale, translation):
return bool(boundary._contains(pattern * scale + translation).all())
# The translations that fit are the boundary shifted back by each
# scaled point, all intersected. Points far apart cut that region the
# most, so they go first, behind those that cut it at the last scale.
shapely.prepare(geom)
spread = _spread_order(pattern)
# Inside a convex boundary the pattern fits exactly where the corners
# of its hull do, so only they cut the region, which a farm of 380
# turbines otherwise cut once per turbine.
corners = _hull_corners(pattern) if _is_convex(geom) else None
if corners is not None:
spread = spread[np.isin(spread, corners)]
cutting = spread[:0]
edges = _edges(geom)
def feasible_region(scale):
"""The translations fitting the scaled pattern inside, or ``None``."""
nonlocal cutting
# Only the turbines must fit, not the pattern's centroid, which may
# sit over a hole or between the parts of a multipart boundary.
order = np.concatenate([cutting, spread[~np.isin(spread, cutting)]])
shifts = scale * pattern
region = _shifted(geom, -shifts[order[0]])
used = [order[0]]
pending = order[1:]
# Points whose shifted boundary holds the region's box cannot cut
# it, and the box only shrinks, so once held they are dropped for
# good. Checked after every second cut, a pattern of 1,238 along
# the edges of a holed boundary checked a million boxes per fit,
# so the checks follow cuts 1, 3, 7, 15 and on.
check = 1
while len(pending):
index, pending = pending[0], pending[1:]
region = region.intersection(_shifted(geom, -shifts[index]))
used.append(index)
if region.is_empty:
cutting = np.array(used)
return None
if len(used) - 1 == check and len(pending):
check = 2 * check + 1
pending = pending[~_holds(geom, region.bounds, shifts[pending])]
cutting = np.array(used)
return region
def placement(scale):
"""A fitting translation and the region it came from, if built."""
# Keep the last verified placement while it still fits, so a
# marginally larger scale does not relocate the whole pattern.
if best_scale > 0 and fits(scale, best_translation):
return best_translation, None
region = feasible_region(scale)
if region is None:
return None, None
# representative_point() stays inside non-convex regions, and the
# direct check rejects points that rounding left a hair outside.
point = region.representative_point()
translation = np.array([point.x, point.y])
return (translation, region) if fits(scale, translation) else (None, None)
def probe(scale):
"""Whether the scaled pattern fits, narrowing the bracket either way."""
nonlocal lower, upper, best_scale, best_translation
translation, region = placement(scale)
if translation is None:
upper = scale
return False
lower = best_scale = scale
best_translation = translation
if region is not None:
built(region)
return True
def built(region):
"""Keep ``region`` as the one at ``best_scale``."""
nonlocal best_region, best_used, region_scale
best_region, best_used, region_scale = region, cutting, best_scale
def best():
"""The region at ``best_scale``, built if it was not."""
if region_scale != best_scale:
built(feasible_region(best_scale))
return best_region
# A hint that passes the direct containment check bounds the scale.
if hint is not None:
hint_scale, hint_translation = hint
hint_scale = float(hint_scale)
hint_translation = np.asarray(hint_translation, dtype=float)
if 0 < hint_scale < upper:
if fits(hint_scale, hint_translation):
lower = best_scale = hint_scale
best_translation = hint_translation
# Feasibility is not monotone in scale on concave, holed or multipart
# boundaries, so bisecting from zero can miss the largest fit. Scan
# down from the upper bound, then search the step that first fits.
scale = upper
for _ in range(_SCAN_STEPS):
if scale <= lower or probe(scale):
break
scale *= _SCAN_RATIO
# Near the largest scale the region is a small polygon, each edge on a
# boundary edge shifted back by one scaled point. Those lines move
# with the scale, so where they meet predicts the largest scale.
predict = True
for _ in range(self.max_iterations):
if upper - lower <= self.tolerance * max(upper, 1.0):
break
meeting = None
if predict and best_scale > 0 and best() is not None:
meeting = _meeting_scale(
best_region, best_scale, pattern, best_used, edges
)
if meeting is None or not lower < meeting[0] < upper:
probe(0.5 * (lower + upper))
predict = True
continue
width = upper - lower
scale, cutting = meeting
# The points whose lines meet cut the region first either side.
margin = _MEETING_MARGIN * self.tolerance * max(scale, 1.0)
if (scale - margin <= lower or probe(scale - margin)) and (
scale + margin < upper
):
probe(scale + margin)
# A prediction that does not halve the bracket is followed by a
# bisection.
predict = upper - lower <= 0.5 * width
# The search keeps the first placement that fits, which leaves the
# pattern against one edge wherever its region has slack. The centroid
# of the region at the final scale shares that slack out.
if best_scale > 0:
region = best()
if region is not None:
center = np.asarray(region.centroid.coords[0])
if fits(best_scale, center):
best_translation = center
positions = pattern * best_scale + best_translation
return positions, {
"method": "rigid_pattern",
"scale": float(best_scale),
"translation": best_translation.tolist(),
}
[docs]
class SelectedCandidateFitter(Fitter):
"""Fit a scalable candidate field and occupy exactly ``N`` positions.
The field is tried at ``n_phases`` by ``n_phases`` offsets. The largest
feasible scale is found on ``search_intervals`` geometric steps, down to a
bottom on a ladder of quarter octaves so that fits at nearby counts share
their steps, then on ``max_iterations`` linear steps inside the best one.
A step at which the layout bounds every field below ``N`` candidates is
passed over without building its fields. The fields of the last few
boundaries are kept, so fitting a nearby count again, as a spacing search
does, builds few new ones, and its result is that of a first fit. Fields
above ``candidate_cap`` points are skipped. ``max_evaluations`` caps the
field evaluations, kept ones included, and defaults to the most the search
can need. The selected positions are then enlarged by ``pattern_fitter``
when the layout allows.
"""
fit_method = "selected_candidates"
def __init__(
self,
*,
n_phases: int = 6,
max_iterations: int = 40,
pattern_fitter: RigidPatternFitter | None = None,
candidate_cap: int = 100_000,
max_evaluations: int | None = None,
search_intervals: int = 128,
):
n_phases = _validate_integer("n_phases", n_phases)
max_iterations = _validate_integer("max_iterations", max_iterations)
candidate_cap = _validate_integer("candidate_cap", candidate_cap)
if max_evaluations is not None:
max_evaluations = _validate_integer("max_evaluations", max_evaluations)
search_intervals = _validate_integer(
"search_intervals", search_intervals, minimum=2
)
# Each phase is probed at most once per scale, at the initial estimate,
# 32 shrinking floors, the scan's search_intervals + 1 scales and the
# max_iterations refinements.
required = n_phases**2 * (search_intervals + max_iterations + 34)
if max_evaluations is None:
max_evaluations = required
elif required > max_evaluations:
raise ValueError(
f"max_evaluations of {max_evaluations} is too small for "
f"{n_phases}x{n_phases} phases, which need {required}"
)
self.n_phases = n_phases
self.max_iterations = max_iterations
self.pattern_fitter = pattern_fitter or RigidPatternFitter()
self.candidate_cap = candidate_cap
self.max_evaluations = max_evaluations
self.search_intervals = search_intervals
def _phases(self):
fractions = (np.arange(self.n_phases) + 0.5) / self.n_phases
return tuple((float(x), float(y)) for x in fractions for y in fractions)
@staticmethod
def _contained(points, boundary) -> bool:
"""Whether every one of ``points`` lies inside ``boundary``."""
if not len(points):
return True
candidates = np.asarray(points, dtype=float).reshape(-1, 2)
return bool(boundary._contains(candidates).all())
[docs]
def fit(
self,
layout: _CandidateLayout,
boundary: Boundary,
pattern: np.ndarray | None,
rng: np.random.Generator,
hint: tuple[float, np.ndarray] | None = None,
) -> tuple[np.ndarray, dict]:
"""Return exactly ``N`` positions from a bounded candidate-scale search.
Candidate count is not monotone in scale on concave polygons, so the
whole scale interval is scanned for the largest feasible scale before
refining it, passing over the scales no field can fill. Candidate
allocation and field evaluations are both capped.
"""
# The layout supplies the candidate field, so there is no pattern or hint.
del pattern, hint
n_turbines = layout._required_n
if n_turbines > self.candidate_cap:
raise ValueError(f"n_turbines exceeds candidate_cap={self.candidate_cap}")
phases = self._phases()
required = layout._required_candidates(boundary)
evaluations = 0
cache: dict[float, list[np.ndarray | None]] = {}
memo = _field_memo(self, layout, boundary, phases)
def candidates(scale, hopeful=None):
"""Return phase fields. ``None`` means the allocation cap was hit.
A phase ``hopeful`` rules out holds fewer than ``N`` candidates, so
its field is left empty without being built.
"""
nonlocal evaluations
scale = float(scale)
if scale in cache:
return cache[scale]
built = range(len(phases)) if hopeful is None else np.flatnonzero(hopeful)
if evaluations + len(built) > self.max_evaluations:
raise ValueError(
f"candidate search exceeded max_evaluations={self.max_evaluations}"
)
# A kept field counts as evaluated, so the search and its record
# are those of a first fit.
evaluations += len(built)
known = memo.get(scale)
if known is not None and known[1] <= n_turbines:
cache[scale] = known[0]
return known[0]
if len(built) == len(phases):
fields = layout._candidate_fields(
boundary, scale, phases, self.candidate_cap, n_turbines
)
else:
fields: list[np.ndarray | None] = [np.empty((0, 2))] * len(phases)
subset = tuple(phases[index] for index in built)
for index, field in zip(
built,
layout._candidate_fields(
boundary, scale, subset, self.candidate_cap, n_turbines
),
strict=True,
):
fields[index] = field
# Every field built in full serves any count. One left empty, or
# every one cut short below N, shows no more than that N do not fit.
whole = hopeful is None and any(
field is None or len(field) for field in fields
)
memo[scale] = (fields, 0 if whole else n_turbines)
cache[scale] = fields
return fields
def feasible(scale, hopeful=None):
return any(
points is None or len(points) >= n_turbines
for points in candidates(scale, hopeful)
)
initial = max(float(layout.estimate_spacing(boundary)), 1e-12)
if not np.isfinite(initial):
raise ValueError("estimated candidate spacing must be finite")
if len(required) == n_turbines and self._contained(required, boundary):
origin = layout._candidate_origin(boundary)
metadata = {
"method": self.fit_method,
"scale": initial,
"translation": [0.0, 0.0],
"phase": list(phases[0]),
"candidate_count": int(n_turbines),
"candidate_evaluations": 0,
"candidate_cap": self.candidate_cap,
}
if origin is not None:
metadata["origin"] = np.asarray(origin, dtype=float).tolist()
return np.asarray(required, dtype=float).reshape(-1, 2), metadata
floor = initial
if not feasible(floor):
for _ in range(32):
floor /= 1.5
if feasible(floor):
break
else:
raise ValueError(f"could not generate {n_turbines} candidates")
minx, miny, maxx, maxy = boundary.bounds
diagonal = float(np.hypot(maxx - minx, maxy - miny))
upper = max(initial * 8.0, diagonal * 2.0, floor)
# The grid's bottom lies on a ladder of quarter octaves below its top,
# so fits at nearby counts share their scales and so their fields.
# The feasible floor joins the scan.
bottom = upper * 2.0 ** (-np.ceil(4.0 * np.log2(upper / floor)) / 4.0)
scales = np.unique(
np.append(np.geomspace(bottom, upper, self.search_intervals + 1), floor)
)
# Scales at which no phase can hold N candidates are passed over
# without building their fields.
bound = layout._candidate_bounds(boundary, floor, phases, self.candidate_cap)
def possible(scales):
"""Indices of ``scales`` where a phase can hold ``N`` candidates.
Each comes with those phases, or ``None`` where every phase can.
"""
if bound is None:
return [(index, None) for index in range(len(scales))]
able = bound(scales) >= n_turbines
some, every = able.any(axis=1), able.all(axis=1)
return [
(index, None if every[index] else able[index])
for index in np.flatnonzero(some)
]
# Only the largest feasible scale matters, so scan from the top and stop
# at the first one that fits.
best_index = next(
(
index
for index, hopeful in possible(scales)[::-1]
if feasible(scales[index], hopeful)
),
None,
)
if best_index is None:
raise ValueError(f"could not generate {n_turbines} candidates")
best_scale = float(scales[best_index])
if best_index < len(scales) - 1:
local_scales = np.linspace(
best_scale,
float(scales[best_index + 1]),
self.max_iterations + 2,
)[1:-1]
for index, hopeful in possible(local_scales)[::-1]:
if feasible(local_scales[index], hopeful):
best_scale = float(local_scales[index])
break
fields = candidates(best_scale)
options = [
(index, points)
for index, points in enumerate(fields)
if points is not None and len(points) >= n_turbines
]
if not options:
# Fields that span the boundary's bounding box, as a spiral's or
# a ring's do, need far more points than turbines on a thin boundary.
box = (maxx - minx) * (maxy - miny)
covered = float(boundary.area) / box if box > 0 else 0.0
raise ValueError(
f"every field that fits {n_turbines} candidates exceeds "
f"candidate_cap={self.candidate_cap} points, and the boundary "
f"covers {covered:.3g} of its bounding box. "
f"Raise the cap with {type(layout).__name__}(fitter="
f"{type(self).__name__}(candidate_cap={2 * self.candidate_cap}))"
)
minimum_surplus = min(len(points) - n_turbines for _, points in options)
options = [
option
for option in options
if len(option[1]) - n_turbines == minimum_surplus
]
phase_index, points = options[int(rng.integers(len(options)))]
selected = _select_candidates(points, n_turbines, rng, required)
selected_center = selected.mean(axis=0)
selected_scale = best_scale
if layout.refine_candidates:
normalized = (selected - selected_center) / best_scale
# The selected candidates fit at ``best_scale``, so refinement
# starts there and can only enlarge them.
positions, rigid = self.pattern_fitter.fit(
layout, boundary, normalized, rng, hint=(best_scale, selected_center)
)
# Rounding can fail the hint's containment check, so keep the
# clipped candidates if refinement still shrank them.
if float(rigid["scale"]) < selected_scale * (1.0 - 1e-9):
positions = selected
rigid = {"scale": selected_scale, "translation": [0.0, 0.0]}
else:
positions = selected
rigid = {
"scale": best_scale,
"translation": [0.0, 0.0],
}
metadata = {
"method": self.fit_method,
"scale": rigid["scale"],
"translation": rigid["translation"],
"phase": list(phases[phase_index]),
"candidate_count": int(len(points)),
"candidate_evaluations": evaluations,
"candidate_cap": self.candidate_cap,
}
origin = layout._candidate_origin(boundary)
if origin is not None:
origin = np.asarray(origin, dtype=float)
if layout.refine_candidates:
origin = (origin - selected_center) / best_scale
origin = origin * rigid["scale"] + rigid["translation"]
metadata["origin"] = origin.tolist()
return positions, metadata
[docs]
class SelectedGridFitter(SelectedCandidateFitter):
"""Select exact, well-distributed intersections from an aligned grid.
Arguments are as for :class:`SelectedCandidateFitter`.
"""
fit_method = "selected_grid"
[docs]
class SelectedCircularFitter(SelectedCandidateFitter):
"""Select exact positions from boundary-clipped concentric rings.
Arguments are as for :class:`SelectedCandidateFitter`.
"""
fit_method = "selected_circular"
[docs]
class SelectedContourFitter(SelectedCandidateFitter):
"""Select exact positions from inward offsets of a farm boundary.
Arguments are as for :class:`SelectedCandidateFitter`, with three phases
by default.
"""
fit_method = "selected_contours"
def __init__(self, *, n_phases: int = 3, **kwargs: Any) -> None:
super().__init__(n_phases=n_phases, **kwargs)
#: Fields each fitter evaluated, per layout and boundary, see :func:`_field_memo`.
_FIELD_MEMOS: weakref.WeakKeyDictionary = weakref.WeakKeyDictionary()
#: Boundaries whose fields one fitter keeps, the most recently fitted.
_MEMO_BOUNDARIES = 8
def _field_memo(fitter, layout, boundary, phases) -> dict:
"""The fields ``fitter`` evaluated before for ``layout`` on ``boundary``.
A spacing fits one layout at many counts on one boundary, and the fields
do not depend on the count, so they are kept by scale, each with the
least count it serves. The memo lives beside the fitter, not on it, so a
pickled sampler does not carry it.
"""
memos = _FIELD_MEMOS.setdefault(fitter, [])
key = (layout._field_key(), phases, fitter.candidate_cap)
geometry = boundary.geometry
for index, (entry_key, entry_geometry, scales) in enumerate(memos):
if entry_key == key and (
entry_geometry is geometry or entry_geometry == geometry
):
memos.append(memos.pop(index))
return scales
scales: dict = {}
memos.append((key, geometry, scales))
del memos[:-_MEMO_BOUNDARIES]
return scales
#: Hull area beyond a polygon's own, as a share of it, below which the
#: polygon is taken as convex. Containment is checked directly anyway.
_CONVEX_SLACK = 1e-12
def _is_convex(geometry) -> bool:
"""Whether ``geometry`` is one polygon without holes equal to its hull."""
if geometry.geom_type != "Polygon" or geometry.interiors:
return False
return geometry.convex_hull.area <= geometry.area * (1.0 + _CONVEX_SLACK)
def _hull_corners(points) -> np.ndarray | None:
"""Indices of the corners of the hull of ``points``, or ``None`` if flat."""
from scipy.spatial import ConvexHull, QhullError
if len(points) < 4:
return None
try:
return ConvexHull(points).vertices
except (QhullError, ValueError):
return None
def _select_candidates(points, count, rng, required):
"""Farthest-point selection while retaining required candidate points."""
points = np.asarray(points, dtype=float).reshape(-1, 2)
required = np.asarray(required, dtype=float).reshape(-1, 2)
if not len(required):
return farthest_point_sample(points, count, rng)
if len(required) > count:
raise ValueError("required candidates cannot exceed n_turbines")
available = np.ones(len(points), dtype=bool)
selected = [point.copy() for point in required]
minimum_distances = np.full(len(points), np.inf)
for required_point in required:
deltas = np.abs(points - required_point)
tolerance = np.maximum(
np.spacing(np.abs(points)), np.spacing(np.abs(required_point))
)
matches = np.all(deltas <= tolerance, axis=1)
if not np.any(matches):
raise ValueError("a required point is not in the candidate field")
available &= ~matches
distances = np.linalg.norm(points - required_point, axis=1)
minimum_distances = np.minimum(minimum_distances, distances)
minimum_distances[~available] = -np.inf
# Picks are marked rather than masked again, unavailable candidates are
# already below any maximum, and the distances are the norm's own
# products, so the picks are those of the plain loop, several times faster.
x = np.ascontiguousarray(points[:, 0])
y = np.ascontiguousarray(points[:, 1])
dx = np.empty_like(x)
dy = np.empty_like(y)
while len(selected) < count:
maximum = float(np.max(minimum_distances))
if not np.isfinite(maximum):
raise ValueError("candidate field does not contain enough distinct points")
lower_ulp = np.nextafter(maximum, -np.inf)
choices = np.flatnonzero(minimum_distances >= lower_ulp)
choice = int(rng.choice(choices))
selected.append(points[choice].copy())
np.subtract(x, x[choice], out=dx)
np.subtract(y, y[choice], out=dy)
np.multiply(dx, dx, out=dx)
np.multiply(dy, dy, out=dy)
np.add(dx, dy, out=dx)
np.sqrt(dx, out=dx)
np.minimum(minimum_distances, dx, out=minimum_distances)
minimum_distances[choice] = -np.inf
return np.asarray(selected)
def _spread_order(points):
"""Indices of ``points`` with the first few in farthest-point order.
The order starts at the point farthest from the origin, and the points
after the first :data:`_SPREAD_POINTS` keep their order.
"""
count = len(points)
head = min(count, _SPREAD_POINTS)
order = np.empty(head, dtype=int)
order[0] = int(np.argmax(np.einsum("ij,ij->i", points, points)))
distances = np.full(count, np.inf)
# A chosen point is never chosen again, even among repeated points.
distances[order[0]] = -np.inf
for position in range(1, head):
delta = points - points[order[position - 1]]
distances = np.minimum(distances, np.einsum("ij,ij->i", delta, delta))
order[position] = int(np.argmax(distances))
distances[order[position]] = -np.inf
rest = np.ones(count, dtype=bool)
rest[order] = False
return np.concatenate([order, np.flatnonzero(rest)])
def _shifted(geometry, offset):
"""``geometry`` translated by the ``(dx, dy)`` ``offset``."""
return shapely.transform(geometry, lambda coordinates: coordinates + offset)
def _holds(geometry, bounds, shifts):
"""Whether ``geometry`` holds the box ``bounds`` moved by each of ``shifts``.
The box must lie strictly inside, off the edge.
"""
xmin, ymin, xmax, ymax = bounds
x, y = np.asarray(shifts).T
boxes = shapely.box(xmin + x, ymin + y, xmax + x, ymax + y)
return shapely.contains_properly(geometry, boxes)
def _edges(geometry):
"""Start, direction and squared length of each edge of ``geometry``.
The interior lies to the left of every edge.
"""
rings = shapely.get_rings(shapely.get_parts(shapely.orient_polygons(geometry)))
coordinates, ring = shapely.get_coordinates(rings, return_index=True)
joined = ring[1:] == ring[:-1]
starts = coordinates[:-1][joined]
deltas = np.diff(coordinates, axis=0)[joined]
return starts, deltas, np.einsum("ij,ij->i", deltas, deltas)
def _meeting_scale(region, scale, pattern, used, edges):
"""The largest scale at which the lines of ``region``'s edges still meet.
``region`` holds the translations that fit ``pattern`` at ``scale``, and
each of its edges lies on one of the boundary's ``edges`` shifted back by
one scaled point among ``used``. Each such line moves with the scale.
Return that scale with the points whose lines meet there, or ``None`` when
no three lines meet or ``region`` has more than :data:`_MEETING_LINES`
edges.
"""
if region.geom_type != "Polygon":
return None
ring = shapely.get_coordinates(region.exterior)
if not 3 < len(ring) <= _MEETING_LINES + 1:
return None
starts, deltas, squared = edges
used = np.asarray(used)
# The boundary edge and the point that each region edge's middle lies on.
middles = 0.5 * (ring[:-1] + ring[1:])
touching = middles[:, None, :] + scale * pattern[used]
x = touching[..., 0].reshape(-1, 1) - starts[:, 0]
y = touching[..., 1].reshape(-1, 1) - starts[:, 1]
along = np.clip((x * deltas[:, 0] + y * deltas[:, 1]) / squared, 0.0, 1.0)
x -= along * deltas[:, 0]
y -= along * deltas[:, 1]
distance = (x * x + y * y).reshape(len(middles), -1)
point, edge = np.divmod(np.unique(distance.argmin(axis=1)), len(starts))
if len(point) < 3:
return None
point = used[point]
# Each line keeps the translations t with n . (t + s * p) <= n . start.
normals = np.column_stack([deltas[edge, 1], -deltas[edge, 0]])
normals /= np.sqrt(squared[edge])[:, None]
lines = np.column_stack([normals, np.einsum("ij,ij->i", normals, pattern[point])])
bounds = np.einsum("ij,ij->i", normals, starts[edge] - middles.mean(axis=0))
# Cramer's rule solves every three lines for the translation and scale.
triples = _triples(len(lines))
rows = lines[triples]
cofactors = np.cross(rows[:, [1, 2, 0]], rows[:, [2, 0, 1]])
determinant = np.einsum("ij,ij->i", rows[:, 0], cofactors[:, 0])
size = np.sqrt(np.einsum("tij,tij->ti", rows, rows)).prod(axis=1)
solvable = np.flatnonzero(np.abs(determinant) > 1e-9 * size)
vertices = (
np.einsum("ti,tij->tj", bounds[triples[solvable]], cofactors[solvable])
/ determinant[solvable, None]
)
tolerance = 1e-9 * max(np.abs(bounds).max(), scale * np.abs(lines[:, 2]).max())
inside = np.flatnonzero((vertices @ lines.T <= bounds + tolerance).all(axis=1))
if not len(inside):
return None
best = inside[vertices[inside, 2].argmax()]
return float(vertices[best, 2]), np.unique(point[triples[solvable[best]]])
@cache
def _triples(count):
"""Every three distinct indices below ``count``, one row each."""
return np.array(list(combinations(range(count), 3)), dtype=int).reshape(-1, 3)
