242 lines
8.5 KiB
Python
242 lines
8.5 KiB
Python
"""Planar geometry helpers for the OSM → asset pipeline.
|
|
|
|
Pure Python: no `bpy`, so this runs and tests outside Blender. All functions
|
|
work in projected metres (see `osmassets.osm.Projector`) unless the name says
|
|
otherwise; `geometry_rings` and `feature_in_bounds` take raw GeoJSON and are the
|
|
two exceptions, operating on lon/lat.
|
|
|
|
Rings are lists of (x, y) tuples. A repeated closing point is tolerated
|
|
everywhere but never required.
|
|
"""
|
|
|
|
import math
|
|
|
|
|
|
def geometry_rings(geometry):
|
|
"""Exterior rings of a GeoJSON Polygon/MultiPolygon; holes are dropped."""
|
|
if not geometry:
|
|
return []
|
|
kind = geometry.get("type")
|
|
coordinates = geometry.get("coordinates", [])
|
|
if kind == "Polygon":
|
|
return coordinates[:1]
|
|
if kind == "MultiPolygon":
|
|
return [polygon[0] for polygon in coordinates if polygon]
|
|
return []
|
|
|
|
|
|
def feature_in_bounds(feature, projector):
|
|
"""True when any coordinate of the feature falls inside the padded bounds."""
|
|
def walk(value):
|
|
if isinstance(value, list) and value and isinstance(value[0], (int, float)):
|
|
return projector.inside(value)
|
|
return any(walk(v) for v in value) if isinstance(value, list) else False
|
|
return walk(feature.get("geometry", {}).get("coordinates", []))
|
|
|
|
|
|
def clip_polygon(ring, xmin, xmax, ymin, ymax):
|
|
"""Sutherland-Hodgman clip of a ring against an axis-aligned box."""
|
|
if len(ring) < 3:
|
|
return []
|
|
|
|
def clip_edge(points, inside, intersection):
|
|
if not points:
|
|
return []
|
|
result = []
|
|
previous = points[-1]
|
|
previous_inside = inside(previous)
|
|
for current in points:
|
|
current_inside = inside(current)
|
|
if current_inside != previous_inside:
|
|
result.append(intersection(previous, current))
|
|
if current_inside:
|
|
result.append(current)
|
|
previous = current
|
|
previous_inside = current_inside
|
|
return result
|
|
|
|
ring = clip_edge(
|
|
ring, lambda p: p[0] >= xmin,
|
|
lambda a, b: (xmin, a[1] + (b[1] - a[1]) * (xmin - a[0]) /
|
|
(b[0] - a[0]) if b[0] != a[0] else a[1]))
|
|
ring = clip_edge(
|
|
ring, lambda p: p[0] <= xmax,
|
|
lambda a, b: (xmax, a[1] + (b[1] - a[1]) * (xmax - a[0]) /
|
|
(b[0] - a[0]) if b[0] != a[0] else a[1]))
|
|
ring = clip_edge(
|
|
ring, lambda p: p[1] >= ymin,
|
|
lambda a, b: (a[0] + (b[0] - a[0]) * (ymin - a[1]) /
|
|
(b[1] - a[1]) if b[1] != a[1] else a[0], ymin))
|
|
ring = clip_edge(
|
|
ring, lambda p: p[1] <= ymax,
|
|
lambda a, b: (a[0] + (b[0] - a[0]) * (ymax - a[1]) /
|
|
(b[1] - a[1]) if b[1] != a[1] else a[0], ymax))
|
|
return ring
|
|
|
|
|
|
def sample_tree_row(points, spacing, height):
|
|
"""Evenly space (x, y, height) samples along a polyline.
|
|
|
|
The trailing point is appended only when the last regular sample stops well
|
|
short of it, so a row does not end in a double-planted tree.
|
|
"""
|
|
if len(points) < 2:
|
|
return []
|
|
samples = [(points[0][0], points[0][1], height)]
|
|
distance_until_next = spacing
|
|
for start, end in zip(points, points[1:]):
|
|
dx = end[0] - start[0]
|
|
dy = end[1] - start[1]
|
|
segment_length = math.hypot(dx, dy)
|
|
if segment_length == 0:
|
|
continue
|
|
while distance_until_next <= segment_length:
|
|
ratio = distance_until_next / segment_length
|
|
samples.append((start[0] + dx * ratio, start[1] + dy * ratio, height))
|
|
distance_until_next += spacing
|
|
distance_until_next -= segment_length
|
|
last = points[-1]
|
|
if math.hypot(samples[-1][0] - last[0], samples[-1][1] - last[1]) > spacing * 0.45:
|
|
samples.append((last[0], last[1], height))
|
|
return samples
|
|
|
|
|
|
def polygon_area(ring):
|
|
"""Unsigned shoelace area; 0.0 for degenerate rings."""
|
|
return abs(signed_polygon_area(ring))
|
|
|
|
|
|
def signed_polygon_area(ring):
|
|
"""Signed shoelace area; positive for counter-clockwise rings."""
|
|
if len(ring) < 3:
|
|
return 0.0
|
|
area = 0.0
|
|
for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
|
|
area += x1 * y2 - x2 * y1
|
|
return area * 0.5
|
|
|
|
|
|
def sample_ring_boundary(ring, spacing, inset=0.0, max_samples=None):
|
|
"""Evenly sample a closed ring's boundary.
|
|
|
|
Returns (x, y, angle, index) samples. ``angle`` follows the local edge
|
|
direction, and ``inset`` moves the sample toward the polygon interior.
|
|
"""
|
|
if len(ring) > 1 and ring[0] == ring[-1]:
|
|
ring = ring[:-1]
|
|
if len(ring) < 3 or spacing <= 0.0:
|
|
return []
|
|
|
|
edges = []
|
|
perimeter = 0.0
|
|
winding = signed_polygon_area(ring)
|
|
for index, (start, end) in enumerate(zip(ring, ring[1:] + ring[:1])):
|
|
dx = end[0] - start[0]
|
|
dy = end[1] - start[1]
|
|
length = math.hypot(dx, dy)
|
|
if length <= 1e-9:
|
|
continue
|
|
ux = dx / length
|
|
uy = dy / length
|
|
# Counter-clockwise rings have their interior on the left side of each
|
|
# edge; clockwise rings have it on the right.
|
|
inward = (-uy, ux) if winding >= 0.0 else (uy, -ux)
|
|
edges.append((perimeter, start, ux, uy, length, inward, index))
|
|
perimeter += length
|
|
if not edges:
|
|
return []
|
|
|
|
count = max(1, int(perimeter / spacing))
|
|
if max_samples:
|
|
count = min(count, max_samples)
|
|
step = perimeter / count
|
|
samples = []
|
|
edge_cursor = 0
|
|
for sample_index in range(count):
|
|
target = (sample_index + 0.5) * step
|
|
while edge_cursor + 1 < len(edges) and (
|
|
edges[edge_cursor][0] + edges[edge_cursor][4] < target
|
|
):
|
|
edge_cursor += 1
|
|
edge_start, start, ux, uy, length, inward, _ = edges[edge_cursor]
|
|
along = max(0.0, min(length, target - edge_start))
|
|
x = start[0] + ux * along
|
|
y = start[1] + uy * along
|
|
sx = x + inward[0] * inset
|
|
sy = y + inward[1] * inset
|
|
if inset > 0.0 and not point_in_polygon((sx, sy), ring):
|
|
sx, sy = x, y
|
|
samples.append((sx, sy, math.atan2(uy, ux), sample_index))
|
|
return samples
|
|
|
|
|
|
def sample_polygon_interior(ring, spacing, edge_clearance=0.0, max_samples=None,
|
|
seed=0):
|
|
"""Jittered interior samples for sparse planting inside a polygon."""
|
|
if len(ring) > 1 and ring[0] == ring[-1]:
|
|
ring = ring[:-1]
|
|
if len(ring) < 3 or spacing <= 0.0 or polygon_area(ring) <= 1e-9:
|
|
return []
|
|
|
|
xmin = min(x for x, _ in ring)
|
|
xmax = max(x for x, _ in ring)
|
|
ymin = min(y for _, y in ring)
|
|
ymax = max(y for _, y in ring)
|
|
cols = max(1, int(math.ceil((xmax - xmin) / spacing)))
|
|
rows = max(1, int(math.ceil((ymax - ymin) / spacing)))
|
|
samples = []
|
|
for col in range(cols):
|
|
for row in range(rows):
|
|
sample_seed = ((col + 1) * 73856093) ^ ((row + 1) * 19349663) ^ seed
|
|
jx = ((sample_seed * 0.61803398875) % 1.0 - 0.5) * spacing * 0.7
|
|
jy = ((sample_seed * 0.41421356237) % 1.0 - 0.5) * spacing * 0.7
|
|
x = xmin + (col + 0.5) * spacing + jx
|
|
y = ymin + (row + 0.5) * spacing + jy
|
|
if not point_in_polygon((x, y), ring):
|
|
continue
|
|
if edge_clearance > 0.0 and distance_to_ring((x, y), ring) < edge_clearance:
|
|
continue
|
|
samples.append((x, y, sample_seed))
|
|
|
|
if max_samples and len(samples) > max_samples:
|
|
samples.sort(key=lambda item: (item[2] * 0.754877666) % 1.0)
|
|
samples = samples[:max_samples]
|
|
return samples
|
|
|
|
|
|
def point_in_polygon(point, ring):
|
|
x, y = point
|
|
inside = False
|
|
j = len(ring) - 1
|
|
for i, (xi, yi) in enumerate(ring):
|
|
xj, yj = ring[j]
|
|
crosses = ((yi > y) != (yj > y))
|
|
if crosses:
|
|
x_at_y = (xj - xi) * (y - yi) / (yj - yi) + xi
|
|
if x < x_at_y:
|
|
inside = not inside
|
|
j = i
|
|
return inside
|
|
|
|
|
|
def distance_to_ring(point, ring):
|
|
"""Shortest distance from a point to the ring's edges (not its interior)."""
|
|
px, py = point
|
|
best = float("inf")
|
|
count = len(ring)
|
|
for index in range(count):
|
|
ax, ay = ring[index]
|
|
bx, by = ring[(index + 1) % count]
|
|
dx = bx - ax
|
|
dy = by - ay
|
|
length_sq = dx * dx + dy * dy
|
|
if length_sq <= 1e-9:
|
|
distance = math.hypot(px - ax, py - ay)
|
|
else:
|
|
t = ((px - ax) * dx + (py - ay) * dy) / length_sq
|
|
t = max(0.0, min(1.0, t))
|
|
distance = math.hypot(px - (ax + t * dx), py - (ay + t * dy))
|
|
if distance < best:
|
|
best = distance
|
|
return best
|