ConwayLife.com

Code: Select all

import numpy as np
from scipy.spatial import Voronoi, cKDTree
from PIL import Image, ImageDraw
import argparse
import time
from collections import deque

# Default Penrose rhomb tiling vertices (same as before)
DEFAULT_POINTS = [
    (0.000000, 0.000000), (1.000000, 0.000000), (0.500000, -0.363271),
    (0.309017, -0.951057), (-0.500000, -0.363271), (-0.690983, -0.951057),
    (-0.500000, -1.538842), (0.000000, -1.902113), (0.500000, -1.538842),
    (1.000000, -1.902113), (0.190983, -2.489898), (0.690983, -2.853170),
    (1.309017, -2.853170), (1.809017, -2.489898), (1.618034, -1.902113),
    (2.118034, -1.538842), (2.427051, -2.489898), (2.927051, -2.126627),
    (3.118034, -1.538842), (2.927051, -0.951057), (2.309017, -0.951057),
    (2.118034, -0.363271), (3.118034, -0.363271), (2.927051, 0.224514),
    (2.427051, 0.587785), (1.809017, 0.587785), (1.618034, -0.000000),
    (1.309017, 0.951057), (0.690983, 0.951057), (0.190983, 0.587785),
    (-0.809017, 0.587785), (-0.000000, 1.175571), (0.500000, 1.538842),
    (1.118034, 1.538842), (2.118034, 1.538842), (2.927051, 0.951057),
    (3.427051, 0.587785), (3.618034, -0.000000), (3.927051, -0.951057),
    (3.618034, -1.902113), (3.427051, -2.489898), (2.927051, -2.853170),
    (2.118034, -3.440955), (1.118034, -3.440955), (0.500000, -3.440955),
    (0.000000, -3.077684), (-0.809017, -2.489898), (-1.118034, -1.538842),
    (-1.309017, -0.951057), (-1.118034, -0.363271), (1.309017, -0.951057),
    (-0.190983, -0.587785), (-0.190983, -1.314328), (0.500000, -2.265384),
    (1.190983, -2.489898), (2.309017, -2.126627), (2.736068, -1.538842),
    (2.736068, -0.363271), (2.309017, 0.224514), (1.190983, 0.587785),
    (0.500000, 0.363271), (-0.190983, 0.587785), (1.618034, 1.175571),
    (2.309017, 0.951057), (3.427051, -0.587785), (3.427051, -1.314328),
    (2.309017, -2.853170), (1.618034, -3.077684), (-0.190983, -2.489898),
    (-0.618034, -1.902113), (-0.618034, -0.000000), (1.118034, -0.363271),
    (0.690983, -0.951057), (1.118034, -1.538842), (1.809017, -1.314328),
    (1.809017, -0.587785),
]

def read_points(filename):
    points = []
    with open(filename, 'r') as f:
        for line in f:
            line = line.strip()
            if line:
                x, y = map(float, line.split())
                points.append([x, y])
    return np.array(points)

def segment_distance(px, py, ax, ay, bx, by):
    """Distance from point (px,py) to line segment (ax,ay)-(bx,by)."""
    dx, dy = bx - ax, by - ay
    len2 = dx*dx + dy*dy
    if len2 == 0:
        return np.sqrt((px-ax)**2 + (py-ay)**2)
    t = np.clip(((px-ax)*dx + (py-ay)*dy) / len2, 0.0, 1.0)
    cx, cy = ax + t*dx, ay + t*dy
    return np.sqrt((px-cx)**2 + (py-cy)**2)



def main():
    parser = argparse.ArgumentParser(
        description="Voronoi surface heightmap with quarter-circle fillets")
    parser.add_argument("--pts", metavar="FILE", default=None,
                        help="input point file (x y per line); uses default Penrose set if omitted")
    parser.add_argument("--size", type=int, default=800,
                        help="image size in pixels (default: 800)")
    parser.add_argument("--max-coeff", type=float, default=0.4,
                        help="fillet radius as fraction of median Voronoi edge length (default: 0.2)")
    parser.add_argument("--edges", action="store_true",
                        help="draw finite Voronoi edges in red")
    parser.add_argument("--dots", action="store_true",
                        help="draw input points in blue")
    parser.add_argument("--output", default="voronoi_surface.png",
                        help="output image filename (default: voronoi_surface.png)")
    parser.add_argument("--clip-center", nargs=2, type=float, metavar=("X", "Y"),
                        default=None, help="center of circular clip region")
    parser.add_argument("--clip-radius", type=float, default=None,
                        help="radius: keep only complete cells inside circle")
    args = parser.parse_args()

    # Load points
    if args.pts:
        points = read_points(args.pts)
        print(f"Loaded {len(points)} points from {args.pts}")
    else:
        points = np.array(DEFAULT_POINTS)
        print(f"Using default Penrose point set ({len(points)} points)")

    # Compute Voronoi diagram
    vor = Voronoi(points)

    # Bounding box with margin for rendering — uniform scaling enforced
    margin = 0.5
    raw_x_lo = points[:, 0].min() - margin
    raw_x_hi = points[:, 0].max() + margin
    raw_y_lo = points[:, 1].min() - margin
    raw_y_hi = points[:, 1].max() + margin
    # Make square by expanding the shorter axis to match the longer
    x_span = raw_x_hi - raw_x_lo
    y_span = raw_y_hi - raw_y_lo
    span = max(x_span, y_span)
    x_lo = (raw_x_lo + raw_x_hi) / 2 - span / 2
    x_hi = (raw_x_lo + raw_x_hi) / 2 + span / 2
    y_lo = (raw_y_lo + raw_y_hi) / 2 - span / 2
    y_hi = (raw_y_lo + raw_y_hi) / 2 + span / 2

    # Print point set geometry info
    pt_cx = (points[:, 0].min() + points[:, 0].max()) / 2
    pt_cy = (points[:, 1].min() + points[:, 1].max()) / 2
    pt_w  = points[:, 0].max() - points[:, 0].min()
    pt_h  = points[:, 1].max() - points[:, 1].min()
    print(f"Point set bbox: ({points[:,0].min():.4f},{points[:,1].min():.4f}) to ({points[:,0].max():.4f},{points[:,1].max():.4f})")
    print(f"Point set centroid: ({pt_cx:.4f}, {pt_cy:.4f})")
    print(f"Point set size: {pt_w:.4f} x {pt_h:.4f}")

    # Clip box: input bounding box expanded by factor 2
    # Voronoi vertices outside this are from degenerate or boundary cells
    # and define the clipped Voronoi diagram region of interest
    cx = (points[:, 0].min() + points[:, 0].max()) / 2
    cy = (points[:, 1].min() + points[:, 1].max()) / 2
    hw = (points[:, 0].max() - points[:, 0].min()) / 2 * 2
    hh = (points[:, 1].max() - points[:, 1].min()) / 2 * 2
    def in_clip(v):
        return abs(v[0] - cx) <= hw and abs(v[1] - cy) <= hh

    # Keep only edges where both vertices are valid and inside the clip box
    finite_edges = []
    clipped_point_indices = set()
    for (p0, p1), ridge in zip(vor.ridge_points, vor.ridge_vertices):
        if ridge[0] >= 0 and ridge[1] >= 0:
            v0 = vor.vertices[ridge[0]]
            v1 = vor.vertices[ridge[1]]
            if in_clip(v0) and in_clip(v1):
                finite_edges.append((v0, v1))
            else:
                # Either input point adjacent to an out-of-clip edge is boundary
                clipped_point_indices.add(p0)
                clipped_point_indices.add(p1)
        else:
            # Infinite ridge — mark both adjacent points as boundary
            clipped_point_indices.add(p0)
            clipped_point_indices.add(p1)
    print(f"Finite Voronoi edges (after clipping): {len(finite_edges)}")
    print(f"Boundary/clipped input points: {len(clipped_point_indices)}")

    # Optional circular clip: keep only cells all of whose vertices are inside
    if args.clip_radius is not None:
        ccx = args.clip_center[0] if args.clip_center else pt_cx
        ccy = args.clip_center[1] if args.clip_center else pt_cy
        cr  = args.clip_radius
        print(f"Circular clip: center=({ccx:.4f},{ccy:.4f}) radius={cr:.4f}")
        def in_circle(v):
            return (v[0]-ccx)**2 + (v[1]-ccy)**2 <= cr**2
        # Find cells with any vertex outside circle -> exclude
        outside_cells = set()
        touched_cells = set()
        for (p0, p1), ridge in zip(vor.ridge_points, vor.ridge_vertices):
            if ridge[0] >= 0 and ridge[1] >= 0:
                v0 = vor.vertices[ridge[0]]
                v1 = vor.vertices[ridge[1]]
                if in_clip(v0) and in_clip(v1):
                    touched_cells.add(p0)
                    touched_cells.add(p1)
                    if not in_circle(v0) or not in_circle(v1):
                        outside_cells.add(p0)
                        outside_cells.add(p1)
        complete_cells = touched_cells - outside_cells
        clipped_point_indices = clipped_point_indices | (
            set(range(len(points))) - complete_cells)
        finite_edges = []
        for (p0, p1), ridge in zip(vor.ridge_points, vor.ridge_vertices):
            if ridge[0] >= 0 and ridge[1] >= 0:
                v0 = vor.vertices[ridge[0]]
                v1 = vor.vertices[ridge[1]]
                if (p0 in complete_cells and p1 in complete_cells
                        and in_clip(v0) and in_clip(v1)):
                    finite_edges.append((v0, v1))
        print(f"Complete cells inside circle: {len(complete_cells)}")
        print(f"Edges after circular clip: {len(finite_edges)}")

        # Recompute bounding box from clipped cell vertices — uniform scaling
        if finite_edges:
            all_verts = np.array([v for e in finite_edges for v in e])
            raw_x_lo = all_verts[:,0].min() - margin
            raw_x_hi = all_verts[:,0].max() + margin
            raw_y_lo = all_verts[:,1].min() - margin
            raw_y_hi = all_verts[:,1].max() + margin
            span = max(raw_x_hi - raw_x_lo, raw_y_hi - raw_y_lo)
            x_lo = (raw_x_lo + raw_x_hi) / 2 - span / 2
            x_hi = (raw_x_lo + raw_x_hi) / 2 + span / 2
            y_lo = (raw_y_lo + raw_y_hi) / 2 - span / 2
            y_hi = (raw_y_lo + raw_y_hi) / 2 + span / 2
            print(f"Recentered rendered area: {x_hi-x_lo:.4f} x {y_hi-y_lo:.4f} point units")

    # Statistical breakdown of edge lengths
    el = np.array([np.linalg.norm(v1 - v0) for v0, v1 in finite_edges])
    print(f"Edge lengths: min={el.min():.4f} max={el.max():.4f} mean={el.mean():.4f} std={el.std():.4f}")
    # Round to 4 decimal places to group near-identical lengths
    from collections import Counter
    length_counts = Counter(round(x, 4) for x in el)
    print("Edge length distribution (decreasing count):")
    for length, count in sorted(length_counts.items(), key=lambda x: -x[1]):
        print(f"  {length:.4f}: {count}")

    # Compute edge lengths and derive fillet max distance
    edge_lengths = [np.linalg.norm(v1 - v0) for v0, v1 in finite_edges]
    median_edge_len = np.median(edge_lengths)
    max_dist = args.max_coeff * median_edge_len
    print(f"Median edge length: {median_edge_len:.4f}")
    print(f"Fillet max distance (max): {max_dist:.4f}  (coeff={args.max_coeff})")
    print(f"Rendered area: {x_hi-x_lo:.4f} x {y_hi-y_lo:.4f} point units")
    print(f"Pixels per point unit: {args.size/(x_hi-x_lo):.2f}")
    # Tight bounding box of actual cell vertices (true geometry extent)
    if finite_edges:
        all_verts = np.array([v for e in finite_edges for v in e])
        gx_lo, gx_hi = all_verts[:,0].min(), all_verts[:,0].max()
        gy_lo, gy_hi = all_verts[:,1].min(), all_verts[:,1].max()
        g_w, g_h = gx_hi - gx_lo, gy_hi - gy_lo
        print(f"Geometry bounding box (true cell extent, point units):")
        print(f"  origin: ({gx_lo:.4f}, {gy_lo:.4f})")
        print(f"  width:  {g_w:.4f}")
        print(f"  height: {g_h:.4f}")
        print(f"  max_dist (fillet radius): {max_dist:.4f}")
        print(f"  fillet as fraction of width: {max_dist/g_w:.4f}")
        print(f"  aspect ratio: {g_w/g_h:.4f}")

    size = args.size
    xs = np.linspace(x_lo, x_hi, size)
    ys = np.linspace(y_lo, y_hi, size)
    X, Y = np.meshgrid(xs, ys)

    # For each pixel find its nearest input point (= its Voronoi cell owner)
    tree = cKDTree(points)
    pixel_coords = np.column_stack([X.ravel(), Y.ravel()])
    _, owner = tree.query(pixel_coords)
    owner = owner.reshape(size, size)
    is_background = np.isin(owner, list(clipped_point_indices))

    # Build map from input point index to its bounding Voronoi edges
    # Each ridge separates two input points — add edge to both owners
    point_edges = {i: [] for i in range(len(points))}
    for (p0, p1), ridge in zip(vor.ridge_points, vor.ridge_vertices):
        if ridge[0] >= 0 and ridge[1] >= 0:
            v0 = vor.vertices[ridge[0]]
            v1 = vor.vertices[ridge[1]]
            if in_clip(v0) and in_clip(v1):
                point_edges[p0].append((v0, v1))
                point_edges[p1].append((v0, v1))

    # Distance field: precompute pixel index lists per cell (one pass),
    # then process each cell directly without any masking scan
    print("Computing distance field...")
    t0 = time.time()
    D = np.full(size * size, np.inf)
    owner_flat = owner.ravel()
    # Precompute pixel coordinates
    col_idx = np.arange(size * size) % size
    row_idx = np.arange(size * size) // size
    PX = xs[col_idx]
    PY = ys[row_idx]
    # Build cell -> pixel index list in one pass (no per-cell scan)
    cell_pixels = {i: [] for i in range(len(points))}
    for flat_idx, cell_idx in enumerate(owner_flat):
        cell_pixels[cell_idx].append(flat_idx)
    # Process each cell: direct index into pixel arrays, no masking
    for cell_idx, edges in point_edges.items():
        if not edges:
            continue
        idxs = np.array(cell_pixels[cell_idx])
        if len(idxs) == 0:
            continue
        px = PX[idxs]
        py = PY[idxs]
        best = np.full(len(idxs), np.inf)
        for (v0, v1) in edges:
            dx, dy = v1[0]-v0[0], v1[1]-v0[1]
            len2 = dx*dx + dy*dy
            if len2 == 0:
                d2 = (px-v0[0])**2 + (py-v0[1])**2
            else:
                t = np.clip(((px-v0[0])*dx + (py-v0[1])*dy) / len2, 0.0, 1.0)
                cx = v0[0] + t*dx
                cy = v0[1] + t*dy
                d2 = (px-cx)**2 + (py-cy)**2
            np.minimum(best, d2, out=best)
        D[idxs] = np.sqrt(best)
    D = D.reshape(size, size)
    print(f"Distance field computed in {time.time() - t0:.2f}s")



    # Apply quarter-circle fillet profile:
    # For d < max: height = sqrt(d * (2*max - d))
    #   (derived from Pythagoras: leg of right triangle with hypotenuse=max, other leg=max-d)
    #   rises from 0 at the Voronoi edge to max at distance max, with vertical tangent at edge
    #   and horizontal tangent joining the plateau — a true quarter-circle cross-section
    # For d >= max: height = max  (flat plateau)
    H = np.where(D < max_dist,
                 np.sqrt(np.clip(D * (2*max_dist - D), 0, None)),
                 max_dist)

    # Background (infinite cells) set to zero height
    H[is_background] = 0.0

    # Normalise to 0-1 for greyscale (0=edge/black, 1=plateau/white)
    H_norm = H / max_dist

    # Render greyscale image
    grey = (H_norm * 255).astype(np.uint8)
    img = Image.fromarray(grey, 'L').convert('RGB')
    draw = ImageDraw.Draw(img)

    def to_pixel(px, py):
        ix = int((px - x_lo) / (x_hi - x_lo) * size)
        iy = int((py - y_lo) / (y_hi - y_lo) * size)
        return ix, iy

    # Draw Voronoi edges in red
    if args.edges:
        for (v0, v1) in finite_edges:
            x1, y1 = to_pixel(v0[0], v0[1])
            x2, y2 = to_pixel(v1[0], v1[1])
            draw.line([x1, y1, x2, y2], fill=(220, 30, 30), width=1)

    # Draw input points in blue
    if args.dots:
        r_dot = max(2, size // 300)
        for (px, py) in points:
            cx, cy = to_pixel(px, py)
            draw.ellipse([cx-r_dot, cy-r_dot, cx+r_dot, cy+r_dot],
                         fill=(30, 30, 200), outline=(0, 0, 0))

    img.save(args.output)
    print(f"Saved to {args.output}")

if __name__ == "__main__":
    main()