959 lines
36 KiB
Python
959 lines
36 KiB
Python
# coding=utf-8
|
|
"""
|
|
Noyau de calcul de l'extension « Voronoi Fill », sans dependance a inkex.
|
|
|
|
Remplit une forme avec un pavage de Voronoi dont les cellules sont separees
|
|
par un filet de largeur donnee. Les cellules sont des surfaces pleines ; le
|
|
filet est le vide qui les separe (ou, en sortie, la forme privee des cellules).
|
|
|
|
Testable avec pytest seul et reutilisable hors Inkscape (scripts, schema
|
|
des parametres).
|
|
|
|
Conventions :
|
|
- un point est un tuple (x, y) ; l'axe y est oriente vers le bas (repere SVG) ;
|
|
- un anneau (`ring`) est une liste de points, ferme implicitement ;
|
|
- une region (`rings`) est une liste d'anneaux lue selon la regle pair-impair
|
|
(un anneau interieur a un autre est un trou), quelle que soit l'orientation ;
|
|
- une cellule produite est une region (liste d'anneaux) : en general un seul
|
|
anneau convexe, plusieurs quand la forme la decoupe ou la troue.
|
|
"""
|
|
|
|
import math
|
|
import random
|
|
|
|
# Distance minimale entre germes de Poisson, rapportee a la taille de cellule :
|
|
# calibree pour que le nombre de cellules soit celui d'un pavage hexagonal de
|
|
# meme pas (voir test_distributions_have_similar_density).
|
|
POISSON_RATIO = 0.75
|
|
|
|
DISTRIBUTIONS = ("random", "poisson", "hexagonal")
|
|
|
|
|
|
class FillError(ValueError):
|
|
"""Erreur de parametrage ; `code` permet a la couche inkex de traduire."""
|
|
|
|
def __init__(self, code, count=0):
|
|
ValueError.__init__(self, code)
|
|
self.code = code
|
|
self.count = count
|
|
|
|
|
|
class _Degenerate(Exception):
|
|
"""Contact non transversal (sommet sur une arete, aretes confondues)."""
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Utilitaires geometriques
|
|
# --------------------------------------------------------------------------
|
|
|
|
def parse_color(value, default=("#b3b3b3", 1.0)):
|
|
"""Couleur Inkscape (entier RGBA decimal ou 0x..., ou #rrggbb[aa]).
|
|
|
|
Renvoie (couleur CSS #rrggbb, opacite entre 0 et 1). Le parametre
|
|
« color » d'Inkscape arrive sous forme d'entier RGBA ; on le decode
|
|
nous-memes pour ne pas dependre de l'API couleur d'inkex, qui a change
|
|
entre les versions 1.x.
|
|
"""
|
|
text = str(value).strip()
|
|
try:
|
|
if text.startswith("#"):
|
|
digits = text[1:]
|
|
if len(digits) == 3:
|
|
digits = "".join(c * 2 for c in digits)
|
|
if len(digits) == 6:
|
|
digits += "ff"
|
|
if len(digits) != 8:
|
|
return default
|
|
number = int(digits, 16)
|
|
else:
|
|
number = int(text, 0)
|
|
except ValueError:
|
|
return default
|
|
number &= 0xFFFFFFFF
|
|
red, green, blue = (number >> 24) & 255, (number >> 16) & 255, (number >> 8) & 255
|
|
alpha = (number & 255) / 255.0
|
|
return "#{:02x}{:02x}{:02x}".format(red, green, blue), round(alpha, 4)
|
|
|
|
|
|
def polylines_to_d(polylines, precision=4):
|
|
"""Donnees `d` d'un chemin SVG : une polyligne par sous-chemin."""
|
|
fmt = "{:.%df},{:.%df}" % (precision, precision)
|
|
parts = []
|
|
for polyline in polylines:
|
|
if len(polyline) < 2:
|
|
continue
|
|
parts.append("M " + fmt.format(*polyline[0]))
|
|
parts.extend("L " + fmt.format(*point) for point in polyline[1:])
|
|
return " ".join(parts)
|
|
|
|
|
|
def rings_to_d(rings, precision=4):
|
|
"""Donnees `d` d'un chemin SVG : un sous-chemin ferme par anneau."""
|
|
parts = [polylines_to_d([ring], precision) + " Z" for ring in rings if len(ring) >= 3]
|
|
return " ".join(parts)
|
|
|
|
|
|
def ring_area(ring):
|
|
"""Aire signee (formule du lacet) ; le signe donne l'orientation."""
|
|
total = 0.0
|
|
n = len(ring)
|
|
for k in range(n):
|
|
x1, y1 = ring[k]
|
|
x2, y2 = ring[(k + 1) % n]
|
|
total += x1 * y2 - x2 * y1
|
|
return total / 2.0
|
|
|
|
|
|
def region_area(rings):
|
|
"""Aire d'une region pair-impair (profondeur d'imbrication de chaque anneau)."""
|
|
total = 0.0
|
|
for k, ring in enumerate(rings):
|
|
others = [r for m, r in enumerate(rings) if m != k]
|
|
depth = 1 if point_in_rings(_inner_point(ring), others) else 0
|
|
total += abs(ring_area(ring)) * (-1 if depth else 1)
|
|
return total
|
|
|
|
|
|
def _inner_point(ring):
|
|
"""Point proche du premier sommet, du cote interieur de l'anneau."""
|
|
a, b, c = ring[-1], ring[0], ring[1]
|
|
sign = 1.0 if ring_area(ring) > 0 else -1.0
|
|
# Bissectrice interieure au sommet b, petit pas relatif a la taille locale.
|
|
ux, uy = _unit_vec(a[0] - b[0], a[1] - b[1])
|
|
vx, vy = _unit_vec(c[0] - b[0], c[1] - b[1])
|
|
step = 1e-6 * max(math.hypot(a[0] - b[0], a[1] - b[1]),
|
|
math.hypot(c[0] - b[0], c[1] - b[1]))
|
|
wx, wy = ux + vx, uy + vy
|
|
if math.hypot(wx, wy) < 1e-12:
|
|
wx, wy = -vy * sign, vx * sign
|
|
wx, wy = _unit_vec(wx, wy)
|
|
candidate = (b[0] + step * wx, b[1] + step * wy)
|
|
if point_in_rings(candidate, [ring]):
|
|
return candidate
|
|
return (b[0] - step * wx, b[1] - step * wy)
|
|
|
|
|
|
def _unit_vec(x, y):
|
|
length = math.hypot(x, y)
|
|
if length == 0:
|
|
return 0.0, 0.0
|
|
return x / length, y / length
|
|
|
|
|
|
def point_in_rings(point, rings):
|
|
"""Regle pair-impair : demi-droite horizontale vers les x croissants."""
|
|
x, y = point
|
|
inside = False
|
|
for ring in rings:
|
|
n = len(ring)
|
|
for k in range(n):
|
|
x1, y1 = ring[k]
|
|
x2, y2 = ring[(k + 1) % n]
|
|
if (y1 > y) != (y2 > y):
|
|
if x < x1 + (y - y1) * (x2 - x1) / (y2 - y1):
|
|
inside = not inside
|
|
return inside
|
|
|
|
|
|
def bbox_of(rings):
|
|
"""(xmin, ymin, xmax, ymax) de tous les points, None si vide."""
|
|
xs = [p[0] for ring in rings for p in ring]
|
|
ys = [p[1] for ring in rings for p in ring]
|
|
if not xs:
|
|
return None
|
|
return min(xs), min(ys), max(xs), max(ys)
|
|
|
|
|
|
def _boxes_overlap(a, b):
|
|
return a[0] <= b[2] and b[0] <= a[2] and a[1] <= b[3] and b[1] <= a[3]
|
|
|
|
|
|
def clean_rings(rings, eps=1e-9):
|
|
"""Retire les points doubles et alignes, et les anneaux degeneres.
|
|
|
|
Des sommets alignes produiraient, pour le filet sur le contour, des bords
|
|
de gelules confondus : on les supprime des le depart.
|
|
"""
|
|
result = []
|
|
for ring in rings:
|
|
points = []
|
|
for p in ring:
|
|
p = (float(p[0]), float(p[1]))
|
|
if not points or math.hypot(p[0] - points[-1][0], p[1] - points[-1][1]) > eps:
|
|
points.append(p)
|
|
while len(points) > 1 and math.hypot(points[0][0] - points[-1][0],
|
|
points[0][1] - points[-1][1]) <= eps:
|
|
points.pop()
|
|
changed = True
|
|
while changed and len(points) >= 3:
|
|
changed = False
|
|
n = len(points)
|
|
for k in range(n):
|
|
a, b, c = points[k - 1], points[k], points[(k + 1) % n]
|
|
cross = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0])
|
|
scale = math.hypot(b[0] - a[0], b[1] - a[1]) * math.hypot(c[0] - b[0], c[1] - b[1])
|
|
if abs(cross) <= 1e-12 * max(scale, 1e-300):
|
|
del points[k]
|
|
changed = True
|
|
break
|
|
if len(points) >= 3 and abs(ring_area(points)) > eps * eps:
|
|
result.append(points)
|
|
return result
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Germes
|
|
# --------------------------------------------------------------------------
|
|
|
|
def hex_density(size):
|
|
"""Nombre de germes par unite d'aire d'un pavage hexagonal de pas `size`."""
|
|
return 2.0 / (math.sqrt(3.0) * size * size)
|
|
|
|
|
|
def random_points(bbox, size, rng):
|
|
"""Germes uniformes, en nombre egal a celui d'un pavage hexagonal de pas `size`."""
|
|
x0, y0, x1, y1 = bbox
|
|
count = max(1, int(round((x1 - x0) * (y1 - y0) * hex_density(size))))
|
|
return [(rng.uniform(x0, x1), rng.uniform(y0, y1)) for _ in range(count)]
|
|
|
|
|
|
def hex_points(bbox, size, irregularity, rng):
|
|
"""Grille hexagonale de pas `size`, chaque germe deplace dans un disque.
|
|
|
|
`irregularity` en % : rayon du disque rapporte a la moitie du pas ;
|
|
0 donne un nid d'abeille parfait.
|
|
"""
|
|
x0, y0, x1, y1 = bbox
|
|
step_y = size * math.sqrt(3.0) / 2.0
|
|
radius = max(0.0, irregularity) / 100.0 * size / 2.0
|
|
points = []
|
|
row = 0
|
|
y = y0
|
|
while y <= y1 + 1e-9:
|
|
x = x0 + (size / 2.0 if row % 2 else 0.0)
|
|
while x <= x1 + 1e-9:
|
|
if radius > 0:
|
|
angle = rng.uniform(0.0, 2.0 * math.pi)
|
|
r = radius * math.sqrt(rng.random())
|
|
points.append((x + r * math.cos(angle), y + r * math.sin(angle)))
|
|
else:
|
|
points.append((x, y))
|
|
x += size
|
|
y += step_y
|
|
row += 1
|
|
return points
|
|
|
|
|
|
def poisson_points(bbox, size, rng, attempts=30):
|
|
"""Echantillonnage de Poisson (Bridson) : germes au moins a POISSON_RATIO*size."""
|
|
x0, y0, x1, y1 = bbox
|
|
radius = POISSON_RATIO * size
|
|
cell = radius / math.sqrt(2.0)
|
|
cols = max(1, int(math.ceil((x1 - x0) / cell)))
|
|
rows = max(1, int(math.ceil((y1 - y0) / cell)))
|
|
grid = [[None] * cols for _ in range(rows)]
|
|
points = []
|
|
|
|
def grid_pos(p):
|
|
return (min(rows - 1, int((p[1] - y0) / cell)),
|
|
min(cols - 1, int((p[0] - x0) / cell)))
|
|
|
|
def fits(p):
|
|
if not (x0 <= p[0] <= x1 and y0 <= p[1] <= y1):
|
|
return False
|
|
gi, gj = grid_pos(p)
|
|
for i in range(max(0, gi - 2), min(rows, gi + 3)):
|
|
for j in range(max(0, gj - 2), min(cols, gj + 3)):
|
|
k = grid[i][j]
|
|
if k is not None:
|
|
q = points[k]
|
|
if (q[0] - p[0]) ** 2 + (q[1] - p[1]) ** 2 < radius * radius:
|
|
return False
|
|
return True
|
|
|
|
def add(p):
|
|
gi, gj = grid_pos(p)
|
|
grid[gi][gj] = len(points)
|
|
points.append(p)
|
|
|
|
add((rng.uniform(x0, x1), rng.uniform(y0, y1)))
|
|
active = [0]
|
|
while active:
|
|
index = rng.randrange(len(active))
|
|
base = points[active[index]]
|
|
for _ in range(attempts):
|
|
angle = rng.uniform(0.0, 2.0 * math.pi)
|
|
r = radius * (1.0 + rng.random())
|
|
candidate = (base[0] + r * math.cos(angle), base[1] + r * math.sin(angle))
|
|
if fits(candidate):
|
|
add(candidate)
|
|
active.append(len(points) - 1)
|
|
break
|
|
else:
|
|
active[index] = active[-1]
|
|
active.pop()
|
|
return points
|
|
|
|
|
|
def make_points(distribution, bbox, size, irregularity=30.0, seed=1):
|
|
"""Germes selon la distribution choisie (graine reproductible)."""
|
|
rng = random.Random(seed)
|
|
if distribution == "random":
|
|
return random_points(bbox, size, rng)
|
|
if distribution == "hexagonal":
|
|
return hex_points(bbox, size, irregularity, rng)
|
|
if distribution == "poisson":
|
|
return poisson_points(bbox, size, rng)
|
|
raise ValueError("distribution inconnue : {}".format(distribution))
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Cellules de Voronoi (retrecies)
|
|
# --------------------------------------------------------------------------
|
|
|
|
def clip_half_plane(poly, nx, ny, c):
|
|
"""Sutherland-Hodgman : garde les points x tels que n.x <= c."""
|
|
result = []
|
|
n = len(poly)
|
|
for k in range(n):
|
|
p = poly[k]
|
|
q = poly[(k + 1) % n]
|
|
dp = nx * p[0] + ny * p[1] - c
|
|
dq = nx * q[0] + ny * q[1] - c
|
|
if dp <= 0:
|
|
result.append(p)
|
|
if (dp < 0 < dq) or (dq < 0 < dp):
|
|
t = dp / (dp - dq)
|
|
result.append((p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1])))
|
|
return result if len(result) >= 3 else []
|
|
|
|
|
|
def voronoi_cells(points, bbox, gap=0.0):
|
|
"""Cellule de chaque germe, retrecie de gap/2 de chaque cote de ses aretes.
|
|
|
|
Chaque cellule est le rectangle `bbox` decoupe par les mediatrices avec
|
|
les germes voisins, decalees de gap/2 vers le germe : deux cellules voisines
|
|
sont ainsi separees exactement de `gap`. Les voisins sont parcourus par
|
|
couronnes d'une grille ; on s'arrete quand aucun germe plus lointain ne
|
|
peut plus couper la cellule (rayon de surete).
|
|
|
|
Renvoie une liste alignee sur `points` : polygone convexe ou None.
|
|
"""
|
|
if not points:
|
|
return []
|
|
x0, y0, x1, y1 = bbox
|
|
area = max((x1 - x0) * (y1 - y0), 1e-12)
|
|
step = math.sqrt(area / len(points))
|
|
cols = max(1, int(math.ceil((x1 - x0) / step)))
|
|
rows = max(1, int(math.ceil((y1 - y0) / step)))
|
|
grid = {}
|
|
positions = []
|
|
for k, p in enumerate(points):
|
|
gi = min(rows - 1, max(0, int((p[1] - y0) / step)))
|
|
gj = min(cols - 1, max(0, int((p[0] - x0) / step)))
|
|
grid.setdefault((gi, gj), []).append(k)
|
|
positions.append((gi, gj))
|
|
max_ring = max(rows, cols)
|
|
half_gap = gap / 2.0
|
|
|
|
cells = []
|
|
for k, p in enumerate(points):
|
|
poly = [(x0, y0), (x1, y0), (x1, y1), (x0, y1)]
|
|
gi, gj = positions[k]
|
|
ring = 0
|
|
while poly and ring <= max_ring:
|
|
for i in range(gi - ring, gi + ring + 1):
|
|
for j in range(gj - ring, gj + ring + 1):
|
|
if max(abs(i - gi), abs(j - gj)) != ring:
|
|
continue
|
|
for m in grid.get((i, j), ()):
|
|
if m == k:
|
|
continue
|
|
q = points[m]
|
|
dx, dy = q[0] - p[0], q[1] - p[1]
|
|
dist = math.hypot(dx, dy)
|
|
if dist == 0:
|
|
continue
|
|
nx, ny = dx / dist, dy / dist
|
|
c = nx * (p[0] + q[0]) / 2.0 + ny * (p[1] + q[1]) / 2.0 - half_gap
|
|
poly = clip_half_plane(poly, nx, ny, c)
|
|
if not poly:
|
|
break
|
|
if not poly:
|
|
break
|
|
if not poly:
|
|
break
|
|
if not poly:
|
|
break
|
|
# Germes non visites : a plus de ring*step de p. Ils ne coupent la
|
|
# cellule que si dist/2 - gap/2 < rayon de la cellule.
|
|
reach = max(math.hypot(v[0] - p[0], v[1] - p[1]) for v in poly)
|
|
if ring * step >= 2.0 * reach + gap:
|
|
break
|
|
ring += 1
|
|
cells.append(poly if poly else None)
|
|
return cells
|
|
|
|
|
|
def arc_steps(radius, angle, tolerance):
|
|
"""Nombre de segments pour qu'une corde s'ecarte de l'arc d'au plus `tolerance`."""
|
|
if radius <= tolerance:
|
|
return max(1, int(math.ceil(abs(angle) / (math.pi / 2))))
|
|
step = 2.0 * math.acos(max(-1.0, 1.0 - tolerance / radius))
|
|
return max(1, int(math.ceil(abs(angle) / step)))
|
|
|
|
|
|
def erode_convex(poly, radius):
|
|
"""Retrecit un polygone convexe de `radius` (aretes decalees vers l'interieur).
|
|
|
|
Exact pour un convexe : intersection des demi-plans de ses aretes
|
|
reculees de `radius`. Renvoie [] si le polygone disparait.
|
|
"""
|
|
if radius <= 0 or not poly:
|
|
return poly
|
|
if ring_area(poly) < 0:
|
|
poly = poly[::-1]
|
|
result = poly
|
|
n = len(poly)
|
|
for k in range(n):
|
|
a, b = poly[k], poly[(k + 1) % n]
|
|
ux, uy = _unit_vec(b[0] - a[0], b[1] - a[1])
|
|
nx, ny = uy, -ux # normale exterieure (aire positive)
|
|
if nx == 0 and ny == 0:
|
|
continue
|
|
result = clip_half_plane(result, nx, ny, nx * a[0] + ny * a[1] - radius)
|
|
if not result:
|
|
return []
|
|
return result
|
|
|
|
|
|
def inradius(poly, iterations=18):
|
|
"""Rayon du plus grand disque inscrit dans un polygone convexe (dichotomie)."""
|
|
if not poly or len(poly) < 3:
|
|
return 0.0
|
|
cx = sum(p[0] for p in poly) / len(poly)
|
|
cy = sum(p[1] for p in poly) / len(poly)
|
|
low, high = 0.0, max(math.hypot(p[0] - cx, p[1] - cy) for p in poly)
|
|
for _ in range(iterations):
|
|
middle = (low + high) / 2.0
|
|
if erode_convex(poly, middle):
|
|
low = middle
|
|
else:
|
|
high = middle
|
|
return low
|
|
|
|
|
|
def round_cell(poly, corner_radius, roundness, tolerance):
|
|
"""Arrondit une cellule convexe sans la faire deborder.
|
|
|
|
Rayon applique : le plus grand de `corner_radius` (longueur fixe) et de
|
|
`roundness` (0 a 1) fois le rayon inscrit de la cellule, ce qui arrondit
|
|
chaque cellule en proportion de sa taille ; a 1, la cellule devient la
|
|
forme la plus ronde contenue dans la cellule d'origine (un disque pour un
|
|
polygone regulier). Erosion puis dilatation du meme rayon : le resultat
|
|
reste inclus dans la cellule, l'ecart du filet est donc respecte.
|
|
Un rayon fixe plus grand que la cellule la fait disparaitre (None).
|
|
"""
|
|
radius = corner_radius
|
|
if roundness > 0:
|
|
rho = inradius(poly)
|
|
# Un peu en dessous du rayon inscrit : le polygone erode garde une
|
|
# taille non nulle, et donc des normales d'aretes bien definies.
|
|
radius = max(radius, min(roundness, 0.97) * rho)
|
|
if radius <= 0:
|
|
return poly
|
|
core = erode_convex(poly, radius)
|
|
if not core:
|
|
return None
|
|
scale = max(abs(p[0]) + abs(p[1]) for p in poly)
|
|
core = clean_rings([core], eps=1e-9 * max(scale, 1.0))
|
|
if not core:
|
|
return None
|
|
return round_convex(core[0], radius, tolerance)
|
|
|
|
|
|
def round_convex(poly, radius, tolerance):
|
|
"""Dilate un polygone convexe de `radius` : aretes decalees, coins en arcs.
|
|
|
|
Applique a une cellule deja retrecie de `radius` en plus, on obtient la
|
|
cellule attendue avec des coins arrondis de rayon `radius`.
|
|
"""
|
|
if radius <= 0 or not poly:
|
|
return poly
|
|
if ring_area(poly) < 0:
|
|
poly = poly[::-1]
|
|
n = len(poly)
|
|
normals = []
|
|
for k in range(n):
|
|
a, b = poly[k], poly[(k + 1) % n]
|
|
ux, uy = _unit_vec(b[0] - a[0], b[1] - a[1])
|
|
# Aire positive : l'exterieur est a droite de chaque arete.
|
|
normals.append((uy, -ux))
|
|
result = []
|
|
for k in range(n):
|
|
v = poly[k]
|
|
n_in = normals[k - 1]
|
|
n_out = normals[k]
|
|
a0 = math.atan2(n_in[1], n_in[0])
|
|
a1 = math.atan2(n_out[1], n_out[0])
|
|
while a1 < a0:
|
|
a1 += 2.0 * math.pi
|
|
steps = arc_steps(radius, a1 - a0, tolerance)
|
|
for s in range(steps + 1):
|
|
a = a0 + (a1 - a0) * s / steps
|
|
result.append((v[0] + radius * math.cos(a), v[1] + radius * math.sin(a)))
|
|
return result
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Decoupe d'une region pair-impair par un polygone convexe
|
|
# --------------------------------------------------------------------------
|
|
|
|
class Region(object):
|
|
"""Region pair-impair indexee : aretes par grille, test d'appartenance par bandes.
|
|
|
|
L'index evite de parcourir tous les sommets de la forme pour chaque cellule
|
|
(texte vectorise, courbes finement aplaties).
|
|
"""
|
|
|
|
def __init__(self, rings, cell=None):
|
|
self.rings = rings
|
|
self.bbox = bbox_of(rings)
|
|
self.edges = []
|
|
for r, ring in enumerate(rings):
|
|
n = len(ring)
|
|
for i in range(n):
|
|
self.edges.append((r, i, ring[i], ring[(i + 1) % n]))
|
|
self.indexed = len(self.edges) > 48 and self.bbox is not None
|
|
if not self.indexed:
|
|
return
|
|
x0, y0, x1, y1 = self.bbox
|
|
if cell is None:
|
|
cell = max(x1 - x0, y1 - y0) / math.sqrt(len(self.edges)) * 2.0
|
|
self.cell = max(cell, 1e-9)
|
|
self.grid = {}
|
|
self.bands = {}
|
|
for e, (_r, _i, a, b) in enumerate(self.edges):
|
|
j0, j1 = self._col(min(a[0], b[0])), self._col(max(a[0], b[0]))
|
|
i0, i1 = self._row(min(a[1], b[1])), self._row(max(a[1], b[1]))
|
|
for i in range(i0, i1 + 1):
|
|
self.bands.setdefault(i, []).append(e)
|
|
for j in range(j0, j1 + 1):
|
|
self.grid.setdefault((i, j), []).append(e)
|
|
|
|
def _col(self, x):
|
|
return int(math.floor((x - self.bbox[0]) / self.cell))
|
|
|
|
def _row(self, y):
|
|
return int(math.floor((y - self.bbox[1]) / self.cell))
|
|
|
|
def edges_in(self, box):
|
|
"""Indices des aretes dont la boite rencontre `box`."""
|
|
if self.bbox is None or not _boxes_overlap(box, self.bbox):
|
|
return []
|
|
if not self.indexed:
|
|
found = range(len(self.edges))
|
|
else:
|
|
found = set()
|
|
for i in range(self._row(box[1]), self._row(box[3]) + 1):
|
|
for j in range(self._col(box[0]), self._col(box[2]) + 1):
|
|
found.update(self.grid.get((i, j), ()))
|
|
result = []
|
|
for e in found:
|
|
a, b = self.edges[e][2], self.edges[e][3]
|
|
if _boxes_overlap(box, (min(a[0], b[0]), min(a[1], b[1]),
|
|
max(a[0], b[0]), max(a[1], b[1]))):
|
|
result.append(e)
|
|
return sorted(result)
|
|
|
|
def contains(self, point):
|
|
"""Appartenance pair-impair."""
|
|
if self.bbox is None:
|
|
return False
|
|
x, y = point
|
|
if not self.indexed:
|
|
candidates = range(len(self.edges))
|
|
else:
|
|
candidates = self.bands.get(self._row(y), ())
|
|
inside = False
|
|
for e in candidates:
|
|
(x1, y1), (x2, y2) = self.edges[e][2], self.edges[e][3]
|
|
if (y1 > y) != (y2 > y):
|
|
if x < x1 + (y - y1) * (x2 - x1) / (y2 - y1):
|
|
inside = not inside
|
|
return inside
|
|
|
|
|
|
def _walk(start_pos, end_pos, n, forward):
|
|
"""Indices des sommets rencontres entre deux positions (indice + fraction).
|
|
|
|
Le sommet k est a la position k ; l'arete k va du sommet k au sommet k+1.
|
|
"""
|
|
i, t = int(start_pos), start_pos - int(start_pos)
|
|
j, s = int(end_pos), end_pos - int(end_pos)
|
|
if forward:
|
|
count = (j - i) % n
|
|
if i == j and s < t:
|
|
count = n
|
|
return [(i + 1 + m) % n for m in range(count)]
|
|
count = (i - j) % n
|
|
if i == j and s > t:
|
|
count = n
|
|
return [(i - m) % n for m in range(count)]
|
|
|
|
|
|
def clip_convex(region, convex, inside=True, eps=1e-9):
|
|
"""Intersection (inside=True) ou difference (False) d'une region et d'un convexe.
|
|
|
|
`region` : Region ou liste d'anneaux (pair-impair). Renvoie une liste
|
|
d'anneaux, a lire aussi en pair-impair.
|
|
|
|
Parcours de type Greiner-Hormann specialise : les croisements du bord de
|
|
la region avec le bord du convexe alternent l'appartenance le long de
|
|
chacun des deux bords. Le resultat est borde par les morceaux de la region
|
|
interieurs (ou exterieurs) au convexe et par les arcs du convexe
|
|
interieurs a la region. Un contact non transversal leve _Degenerate ;
|
|
on reessaie alors avec un convexe tres legerement reduit et decale.
|
|
Si rien n'y fait, on renvoie une region vide : la cellule concernee
|
|
disparait (le filet la recouvre) plutot que de produire un contour faux.
|
|
"""
|
|
if not isinstance(region, Region):
|
|
region = Region(region)
|
|
if len(convex) < 3:
|
|
return [] if inside else [list(r) for r in region.rings]
|
|
if ring_area(convex) < 0:
|
|
convex = convex[::-1]
|
|
cx = sum(p[0] for p in convex) / len(convex)
|
|
cy = sum(p[1] for p in convex) / len(convex)
|
|
scale = max(max(math.hypot(p[0] - cx, p[1] - cy) for p in convex), 1e-12)
|
|
for attempt in range(8):
|
|
if attempt:
|
|
# Perturbation relative infime, donc invisible : reduction (le
|
|
# convexe reste dans la cellule) et decalage selon une direction
|
|
# irrationnelle, car une reduction seule garde les symetries
|
|
# (sommet d'hexagone pile sur un bord vertical, par exemple).
|
|
delta = 1e-9 * 7.3 ** attempt
|
|
f = 1.0 - delta
|
|
dx, dy = 0.618034 * delta * scale, 0.414214 * delta * scale
|
|
shape = [(cx + (p[0] - cx) * f + dx, cy + (p[1] - cy) * f + dy) for p in convex]
|
|
else:
|
|
shape = convex
|
|
try:
|
|
return _clip_convex(region, shape, inside, eps * scale)
|
|
except _Degenerate:
|
|
continue
|
|
return []
|
|
|
|
|
|
def _clip_convex(region, convex, inside, eps):
|
|
m = len(convex)
|
|
box = bbox_of([convex])
|
|
edge_ids = region.edges_in(box)
|
|
|
|
def in_convex(p):
|
|
for j in range(m):
|
|
a, b = convex[j], convex[(j + 1) % m]
|
|
if (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0]) < 0:
|
|
return False
|
|
return True
|
|
|
|
# --- Croisements -------------------------------------------------------
|
|
crossings = [] # [point, ring, subject_pos, convex_pos, entering]
|
|
touched = set()
|
|
for e in edge_ids:
|
|
r, i, p, q = region.edges[e]
|
|
touched.add(r)
|
|
dx, dy = q[0] - p[0], q[1] - p[1]
|
|
seg_len = math.hypot(dx, dy)
|
|
for j in range(m):
|
|
a, b = convex[j], convex[(j + 1) % m]
|
|
ex, ey = b[0] - a[0], b[1] - a[1]
|
|
denom = dx * ey - dy * ex
|
|
wx, wy = a[0] - p[0], a[1] - p[1]
|
|
edge_len = math.hypot(ex, ey)
|
|
if abs(denom) <= 1e-12 * seg_len * edge_len:
|
|
# Paralleles : degenere seulement si confondus et superposes.
|
|
if abs(wx * dy - wy * dx) <= eps * seg_len:
|
|
t0 = (wx * dx + wy * dy) / (seg_len * seg_len)
|
|
t1 = ((b[0] - p[0]) * dx + (b[1] - p[1]) * dy) / (seg_len * seg_len)
|
|
if max(t0, t1) >= 0 and min(t0, t1) <= 1:
|
|
raise _Degenerate()
|
|
continue
|
|
t = (wx * ey - wy * ex) / denom
|
|
s = (wx * dy - wy * dx) / denom
|
|
tol_t = eps / seg_len
|
|
tol_s = eps / edge_len
|
|
if -tol_t <= t <= 1 + tol_t and -tol_s <= s <= 1 + tol_s:
|
|
if t <= tol_t or t >= 1 - tol_t or s <= tol_s or s >= 1 - tol_s:
|
|
raise _Degenerate()
|
|
point = (p[0] + t * dx, p[1] + t * dy)
|
|
entering = ex * dy - ey * dx > 0
|
|
crossings.append([point, r, i + t, j + s, entering])
|
|
|
|
rings = region.rings
|
|
result = []
|
|
|
|
# --- Aucun croisement : anneaux entiers et convexe entier ---------------
|
|
if not crossings:
|
|
if inside:
|
|
for r in sorted(touched):
|
|
if in_convex(rings[r][0]):
|
|
result.append(list(rings[r]))
|
|
# Bord du convexe sans croisement : entierement dans la region
|
|
# ou entierement dehors ; son centre, lui, peut tomber dans un trou.
|
|
if region.contains(convex[0]):
|
|
result.append(list(convex))
|
|
else:
|
|
for r, ring in enumerate(rings):
|
|
if r not in touched or not in_convex(ring[0]):
|
|
result.append(list(ring))
|
|
# Bord du convexe sans croisement : entierement dans la region
|
|
# ou entierement dehors ; son centre, lui, peut tomber dans un trou.
|
|
if region.contains(convex[0]):
|
|
result.append(list(convex))
|
|
return result
|
|
|
|
# --- Chainages le long de chaque anneau et le long du convexe ----------
|
|
by_ring = {}
|
|
for c, x in enumerate(crossings):
|
|
by_ring.setdefault(x[1], []).append(c)
|
|
ring_next, ring_prev = {}, {}
|
|
for r, ids in by_ring.items():
|
|
ids.sort(key=lambda c: crossings[c][2])
|
|
for k, c in enumerate(ids):
|
|
ring_next[c] = ids[(k + 1) % len(ids)]
|
|
ring_prev[c] = ids[k - 1]
|
|
order = sorted(range(len(crossings)), key=lambda c: crossings[c][3])
|
|
conv_next, conv_prev = {}, {}
|
|
for k, c in enumerate(order):
|
|
conv_next[c] = order[(k + 1) % len(order)]
|
|
conv_prev[c] = order[k - 1]
|
|
|
|
# Appartenance a la region de l'arc du convexe qui suit chaque croisement :
|
|
# un seul test, puis alternance (chaque croisement franchit le bord).
|
|
first, second = order[0], conv_next[order[0]]
|
|
arc_points = _walk(crossings[first][3], crossings[second][3], m, True)
|
|
if arc_points:
|
|
probe = convex[arc_points[0]]
|
|
else:
|
|
pa, pb = crossings[first][0], crossings[second][0]
|
|
probe = ((pa[0] + pb[0]) / 2.0, (pa[1] + pb[1]) / 2.0)
|
|
state = region.contains(probe)
|
|
arc_in = {}
|
|
for c in order:
|
|
arc_in[c] = state
|
|
state = not state
|
|
|
|
# --- Parcours ----------------------------------------------------------
|
|
visited = set()
|
|
for start in range(len(crossings)):
|
|
if start in visited:
|
|
continue
|
|
loop = []
|
|
c = start
|
|
guard = 0
|
|
while True:
|
|
guard += 1
|
|
if guard > 4 * len(crossings) + 4:
|
|
raise _Degenerate()
|
|
visited.add(c)
|
|
point, r, pos, _cpos, entering = crossings[c]
|
|
ring = rings[r]
|
|
# Morceau de la region du bon cote du convexe.
|
|
forward = entering if inside else not entering
|
|
nxt = ring_next[c] if forward else ring_prev[c]
|
|
loop.append(point)
|
|
loop.extend(ring[k] for k in _walk(pos, crossings[nxt][2], len(ring), forward))
|
|
visited.add(nxt)
|
|
# Arc du convexe interieur a la region.
|
|
c2 = nxt
|
|
loop.append(crossings[c2][0])
|
|
if arc_in[c2]:
|
|
end = conv_next[c2]
|
|
loop.extend(convex[k] for k in _walk(crossings[c2][3], crossings[end][3], m, True))
|
|
else:
|
|
end = conv_prev[c2]
|
|
loop.extend(convex[k] for k in _walk(crossings[c2][3], crossings[end][3], m, False))
|
|
c = end
|
|
if c == start:
|
|
break
|
|
if c in visited:
|
|
raise _Degenerate()
|
|
if len(loop) >= 3:
|
|
result.append(loop)
|
|
|
|
# Anneaux sans croisement : entiers, selon leur position.
|
|
for r, ring in enumerate(rings):
|
|
if r in by_ring:
|
|
continue
|
|
if inside:
|
|
if r in touched and in_convex(ring[0]):
|
|
result.append(list(ring))
|
|
elif r not in touched or not in_convex(ring[0]):
|
|
result.append(list(ring))
|
|
return result
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Filet sur le contour : erosion locale de la forme
|
|
# --------------------------------------------------------------------------
|
|
|
|
def _edge_capsule(a, b, radius):
|
|
"""Rectangle de demi-largeur `radius` autour de [a, b], raccourci d'un rien.
|
|
|
|
Le leger raccourcissement evite que ses petits cotes passent exactement
|
|
par les sommets de la forme ; le disque du sommet couvre ce vide.
|
|
"""
|
|
ux, uy = _unit_vec(b[0] - a[0], b[1] - a[1])
|
|
length = math.hypot(b[0] - a[0], b[1] - a[1])
|
|
shrink = min(1e-6 * radius, length * 1e-3)
|
|
a = (a[0] + shrink * ux, a[1] + shrink * uy)
|
|
b = (b[0] - shrink * ux, b[1] - shrink * uy)
|
|
nx, ny = -uy * radius, ux * radius
|
|
return [(a[0] + nx, a[1] + ny), (b[0] + nx, b[1] + ny),
|
|
(b[0] - nx, b[1] - ny), (a[0] - nx, a[1] - ny)]
|
|
|
|
|
|
def _vertex_disk(center, radius, tolerance):
|
|
"""Polygone circonscrit au disque (contient le disque vrai)."""
|
|
count = max(8, arc_steps(radius, 2.0 * math.pi, tolerance))
|
|
outer = radius / math.cos(math.pi / count) * (1.0 + 1e-7)
|
|
# Dephasage irrationnel : evite l'alignement avec les rectangles voisins.
|
|
phase = 0.318309886
|
|
return [(center[0] + outer * math.cos(phase + 2.0 * math.pi * k / count),
|
|
center[1] + outer * math.sin(phase + 2.0 * math.pi * k / count))
|
|
for k in range(count)]
|
|
|
|
|
|
def erode_near_boundary(piece, region, radius, tolerance):
|
|
"""Retire de `piece` tout ce qui est a moins de `radius` du bord de `region`.
|
|
|
|
`piece` est deja incluse dans la region ; on lui soustrait, pour chaque
|
|
arete du bord assez proche, un rectangle et un disque a chaque extremite
|
|
(union = voisinage du bord), ce qui revient a l'intersecter avec la forme
|
|
erodee sans jamais calculer cette derniere en entier.
|
|
"""
|
|
box = bbox_of(piece)
|
|
if box is None:
|
|
return []
|
|
near = (box[0] - radius, box[1] - radius, box[2] + radius, box[3] + radius)
|
|
vertices = set()
|
|
for e in region.edges_in(near):
|
|
r, i, a, b = region.edges[e]
|
|
vertices.add((r, i))
|
|
vertices.add((r, (i + 1) % len(region.rings[r])))
|
|
piece = _subtract(piece, _edge_capsule(a, b, radius))
|
|
if not piece:
|
|
return []
|
|
for r, i in sorted(vertices):
|
|
piece = _subtract(piece, _vertex_disk(region.rings[r][i], radius, tolerance))
|
|
if not piece:
|
|
return []
|
|
return piece
|
|
|
|
|
|
def _subtract(piece, convex):
|
|
"""Difference piece - convexe, ignoree si les boites sont disjointes."""
|
|
if not _boxes_overlap(bbox_of(piece), bbox_of([convex])):
|
|
return piece
|
|
return [ring for ring in clip_convex(piece, convex, inside=False) if len(ring) >= 3]
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Remplissage complet
|
|
# --------------------------------------------------------------------------
|
|
|
|
def fill_shape(rings, cell_size, net_width, distribution="poisson", irregularity=30.0,
|
|
seed=1, border=True, corner_radius=0.0, roundness=0.0, tolerance=None,
|
|
max_cells=20000, min_area=None, min_thickness=None):
|
|
"""Cellules de Voronoi retrecies qui remplissent la region `rings`.
|
|
|
|
Renvoie une liste de cellules, chacune une liste d'anneaux (pair-impair).
|
|
Deux cellules voisines sont separees de `net_width` ; avec `border`, les
|
|
cellules restent aussi a net_width du contour de la forme : le filet
|
|
borde la forme d'un cadre de meme largeur que ses brins.
|
|
|
|
`corner_radius` (longueur) et `roundness` (0 a 1, en proportion du rayon
|
|
inscrit de chaque cellule) arrondissent les cellules, le plus grand des
|
|
deux l'emportant (voir round_cell) ; les coins issus du contour de la
|
|
forme restent vifs.
|
|
|
|
Les eclats trop petits (`min_area`) ou trop fins (`min_thickness`,
|
|
epaisseur moyenne 2*aire/perimetre) laisses par le contour sont ecartes :
|
|
le filet y est simplement un peu plus large.
|
|
|
|
Leve FillError("too_many_cells") si le nombre de germes depasse `max_cells`.
|
|
"""
|
|
if cell_size <= 0:
|
|
raise FillError("cell_size")
|
|
if net_width < 0 or corner_radius < 0 or roundness < 0:
|
|
raise FillError("negative")
|
|
rings = clean_rings(rings)
|
|
if not rings:
|
|
return []
|
|
shape_box = bbox_of(rings)
|
|
if tolerance is None:
|
|
tolerance = cell_size / 200.0
|
|
if min_area is None:
|
|
min_area = 0.002 * cell_size * cell_size
|
|
if min_thickness is None:
|
|
min_thickness = 0.5 * net_width
|
|
|
|
# Germes sur la boite de la forme elargie d'une cellule : les cellules du
|
|
# bord ont ainsi des voisins exterieurs et une taille normale.
|
|
margin = cell_size
|
|
seed_box = (shape_box[0] - margin, shape_box[1] - margin,
|
|
shape_box[2] + margin, shape_box[3] + margin)
|
|
estimate = (seed_box[2] - seed_box[0]) * (seed_box[3] - seed_box[1]) * hex_density(cell_size)
|
|
if estimate > max_cells:
|
|
raise FillError("too_many_cells", int(estimate))
|
|
points = make_points(distribution, seed_box, cell_size, irregularity, seed)
|
|
if len(points) > max_cells:
|
|
raise FillError("too_many_cells", len(points))
|
|
|
|
clip_box = (seed_box[0] - 2 * margin, seed_box[1] - 2 * margin,
|
|
seed_box[2] + 2 * margin, seed_box[3] + 2 * margin)
|
|
cells = voronoi_cells(points, clip_box, net_width)
|
|
|
|
region = Region(rings, cell=cell_size)
|
|
result = []
|
|
for poly in cells:
|
|
if not poly:
|
|
continue
|
|
if corner_radius > 0 or roundness > 0:
|
|
poly = round_cell(poly, corner_radius, roundness, tolerance)
|
|
if not poly:
|
|
continue
|
|
if not _boxes_overlap(bbox_of([poly]), shape_box):
|
|
continue
|
|
piece = clip_convex(region, poly, inside=True)
|
|
if border and net_width > 0 and piece:
|
|
piece = erode_near_boundary(piece, region, net_width, tolerance)
|
|
piece = [ring for ring in piece if _substantial(ring, min_area, min_thickness)]
|
|
if piece:
|
|
result.append(piece)
|
|
return result
|
|
|
|
|
|
def _substantial(ring, min_area, min_thickness):
|
|
"""Anneau assez grand et assez epais pour etre garde."""
|
|
area = abs(ring_area(ring))
|
|
if area < min_area:
|
|
return False
|
|
n = len(ring)
|
|
perimeter = sum(math.hypot(ring[(k + 1) % n][0] - ring[k][0],
|
|
ring[(k + 1) % n][1] - ring[k][1]) for k in range(n))
|
|
return 2.0 * area / perimeter >= min_thickness
|
|
|
|
|
|
def net_rings(shape_rings, cells):
|
|
"""Filet sous forme de region pair-impair : forme + anneaux de toutes les cellules.
|
|
|
|
Les cellules sont disjointes et interieures a la forme : chaque anneau de
|
|
cellule y perce un trou, sans aucune operation booleenne.
|
|
"""
|
|
result = [list(r) for r in clean_rings(shape_rings)]
|
|
for cell in cells:
|
|
result.extend(cell)
|
|
return result
|