# 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