inkscape.circlePacking/packing_core.py

#!/usr/bin/env python3
# coding=utf-8
"""
Noyau geometrique du circle packing.

Ce module ne depend pas d'inkex : il peut etre importe et teste sans Inkscape.
Un contour est represente par une liste d'anneaux (`rings`), chaque anneau
etant une liste de points (x, y) decrivant un polygone ferme implicitement
(le dernier point est relie au premier).
"""

import math
import random


def rings_bbox(rings):
    """Boite englobante (xmin, ymin, xmax, ymax), ou 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 bounding_circle(rings):
    """Centre de la boite englobante et rayon circonscrit, ou None si vide.

    Sert a inscrire un motif quelconque dans l'emplacement circulaire calcule
    par `pack_circles` : une copie mise a l'echelle `r / radius` autour de
    (cx, cy) tient dans le disque de rayon r, quelle que soit sa rotation.
    """
    bbox = rings_bbox(rings)
    if bbox is None:
        return None
    xmin, ymin, xmax, ymax = bbox
    cx = 0.5 * (xmin + xmax)
    cy = 0.5 * (ymin + ymax)
    radius = max(math.hypot(p[0] - cx, p[1] - cy)
                 for ring in rings for p in ring)
    return cx, cy, radius


def point_in_rings(x, y, rings):
    """Test d'appartenance par crossing number, regle pair-impair.

    La parite du nombre de traversees sur l'ensemble des anneaux donne
    directement la gestion des trous : un point situe dans un trou traverse un
    nombre pair de bords depuis l'infini.
    """
    inside = False
    for ring in rings:
        n = len(ring)
        for i in range(n):
            x1, y1 = ring[i]
            x2, y2 = ring[(i + 1) % n]
            if (y1 > y) != (y2 > y):
                # abscisse de l'intersection du segment avec la droite Y = y
                xint = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
                if x < xint:
                    inside = not inside
    return inside


def distance_point_segment(px, py, x1, y1, x2, y2):
    """Distance du point (px, py) au segment [(x1, y1), (x2, y2)]."""
    dx = x2 - x1
    dy = y2 - y1
    seg_len2 = dx * dx + dy * dy
    if seg_len2 == 0.0:
        return math.hypot(px - x1, py - y1)
    t = ((px - x1) * dx + (py - y1) * dy) / seg_len2
    t = max(0.0, min(1.0, t))
    return math.hypot(px - (x1 + t * dx), py - (y1 + t * dy))


def distance_to_rings(x, y, rings):
    """Distance minimale du point au contour (bord exterieur ou trous)."""
    best = float('inf')
    for ring in rings:
        n = len(ring)
        for i in range(n):
            x1, y1 = ring[i]
            x2, y2 = ring[(i + 1) % n]
            d = distance_point_segment(x, y, x1, y1, x2, y2)
            if d < best:
                best = d
    return best


class SpatialGrid:
    """Grille de hachage spatial pour retrouver rapidement les cercles voisins.

    Sans elle, chaque tentative comparerait le point candidat a tous les cercles
    deja poses, ce qui rend les dizaines de milliers de tentatives inutilisables.
    """

    def __init__(self, cell_size):
        self.cell_size = max(cell_size, 1e-9)
        self.cells = {}
        self.circles = []

    def _key(self, x, y):
        return (int(math.floor(x / self.cell_size)),
                int(math.floor(y / self.cell_size)))

    def add(self, x, y, r):
        index = len(self.circles)
        self.circles.append((x, y, r))
        self.cells.setdefault(self._key(x, y), []).append(index)

    def neighbors(self, x, y, reach):
        """Cercles dont le centre tombe dans les cellules couvrant `reach`."""
        span = int(math.ceil(reach / self.cell_size))
        ci, cj = self._key(x, y)
        for i in range(ci - span, ci + span + 1):
            for j in range(cj - span, cj + span + 1):
                for index in self.cells.get((i, j), ()):
                    yield self.circles[index]

    def __len__(self):
        return len(self.circles)


def order_by_proximity(circles, cell_size=None):
    """Reordonne les cercles en chaine du plus proche voisin.

    Deux cercles consecutifs dans la liste sont physiquement proches, ce qui
    minimise les deplacements a vide d'une tete de decoupe laser suivant
    l'ordre du document. Le depart est le cercle de (y, x) minimal, soit le
    coin haut-gauche : le resultat est deterministe.

    La recherche du voisin le plus proche balaie des anneaux de cellules de
    plus en plus larges autour du point courant, et s'arrete des que l'anneau
    suivant ne peut plus contenir mieux. Les cercles consommes sont retires de
    leur cellule, donc la grille se vide au fil du parcours.
    """
    if len(circles) < 2:
        return list(circles)

    if cell_size is None:
        xs = [c[0] for c in circles]
        ys = [c[1] for c in circles]
        width = max(xs) - min(xs)
        height = max(ys) - min(ys)
        area = max(width * height, 1e-9)
        # vise environ un cercle par cellule
        cell_size = math.sqrt(area / len(circles))
    cell_size = max(cell_size, 1e-9)

    def key(x, y):
        return (int(math.floor(x / cell_size)), int(math.floor(y / cell_size)))

    cells = {}
    for index, (x, y, _r) in enumerate(circles):
        cells.setdefault(key(x, y), []).append(index)

    def take(index):
        """Retire un cercle de sa cellule."""
        x, y, _r = circles[index]
        bucket = cells[key(x, y)]
        bucket.remove(index)
        if not bucket:
            del cells[key(x, y)]

    # depart : coin haut-gauche
    current = min(range(len(circles)),
                  key=lambda i: (circles[i][1], circles[i][0]))
    take(current)
    order = [current]

    imin = min(c[0] for c in cells)
    imax = max(c[0] for c in cells)
    jmin = min(c[1] for c in cells)
    jmax = max(c[1] for c in cells)

    remaining = len(circles) - 1
    while remaining:
        cx, cy, _r = circles[current]
        ci, cj = key(cx, cy)
        # au-dela, les anneaux sortent completement de la grille
        k_max = max(ci - imin, imax - ci, cj - jmin, jmax - cj)
        best_index = None
        best_dist = float('inf')
        k = 0
        while k <= k_max:
            # borne inferieure de la distance aux cellules de l'anneau k
            if best_index is not None and (k - 1) * cell_size >= best_dist:
                break
            for i in range(ci - k, ci + k + 1):
                for j in range(cj - k, cj + k + 1):
                    # anneau seul : l'interieur a deja ete balaye
                    if k and abs(i - ci) != k and abs(j - cj) != k:
                        continue
                    for index in cells.get((i, j), ()):
                        x, y, _r2 = circles[index]
                        d = math.hypot(x - cx, y - cy)
                        if d < best_dist:
                            best_dist = d
                            best_index = index
            k += 1

        current = best_index
        take(current)
        order.append(current)
        remaining -= 1

    return [circles[i] for i in order]


def pack_circles(rings, min_radius, max_radius, gap=0.0, margin=0.0,
                 attempts=20000, max_circles=0, rng=None):
    """Greedy random circle packing. Renvoie la liste des (x, y, r) places.

    A chaque tentative, un point est tire au hasard dans la boite englobante ;
    s'il est interieur, on y place le plus grand cercle qui respecte `margin`
    vis-a-vis du contour et `gap` vis-a-vis des cercles deja poses.

    La liste renvoyee est ordonnee par `order_by_proximity`.
    """
    rng = rng or random
    bbox = rings_bbox(rings)
    if bbox is None:
        return []
    xmin, ymin, xmax, ymax = bbox
    if xmax - xmin <= 0 or ymax - ymin <= 0:
        return []

    grid = SpatialGrid(max_radius + gap)
    # Un voisin ne peut contraindre le rayon que si son centre est a moins de
    # son rayon (<= max_radius) + max_radius + gap.
    reach = 2.0 * max_radius + gap

    for _ in range(attempts):
        if max_circles and len(grid) >= max_circles:
            break
        x = rng.uniform(xmin, xmax)
        y = rng.uniform(ymin, ymax)
        if not point_in_rings(x, y, rings):
            continue

        radius = min(max_radius, distance_to_rings(x, y, rings) - margin)
        if radius < min_radius:
            continue

        for cx, cy, cr in grid.neighbors(x, y, reach):
            radius = min(radius, math.hypot(x - cx, y - cy) - cr - gap)
            if radius < min_radius:
                break
        if radius >= min_radius:
            grid.add(x, y, radius)

    # ordre de sortie pense pour la decoupe laser, pas pour le rendu
    return order_by_proximity(grid.circles)