Seamless Penrose Kite and Dart images

[pcallahan]

Note that you don't have to use Penrose tilings. Here is one I made from a patch of a periodic still life (generated from a search).

stilllife.png (316.9 KiB) Viewed 355 times

Maybe I can make custom waffle irons.

[pcallahan]

I realized after a while that I don't want intersecting bubbles. I want the Voronoi cells to have flat tops and fillets on the edges. I was doing that already tediously in FreeCAD and then I realized it is pretty simple to define them in terms of distance from voronoi edges. This will work for any point set not just Penrose.

This is Claude-generated code but there is substantial human input. It took an obviously distorted 3d print for me to realize it was not scaling uniformly (though some of my problems in GIMP and Inkscape made sense in retrospect). I also had to get it to find a reasonably efficient algorithm for larger bitmaps and realized (without its help) that I only need to compare distance with the cells you're actually inside. It figured out how to vectorize numpy with that.

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()

Run with arguments "--clip-center -17.45 1.5388 --clip-radius 4.45 --pts pen7.txt --size 2000 --max-coeff 0.5" Using the attached input.

pen7.txt

(60.64 KiB) Downloaded 4 times

The radial clipping centers on this almost symmetric arrangement.

voronoi_surface.png (767.76 KiB) Viewed 347 times

I used imagetostl to get the surface from it and here's a screenshot. Next step is to print it. It looks like it will be over 3 hours for coaster size but I like this layout and may make a silicone mold.

Screenshot 2026-06-07 104401.png (88.68 KiB) Viewed 347 times

[pcallahan]

The 3d print in PLA.

20260607_151402.jpg (2.75 MiB) Viewed 320 times

The silicone mold and a paraffin cast. But the trick is that I floated it for a few seconds in 60°C water to melt off the 3d print texture. The result is not smooth but mottled.

20260608_193050.jpg (2.08 MiB) Viewed 320 times

Finally after molding and casting the paraffin, here's how it looks in resin.

20260610_075434 (1).jpg (2.97 MiB) Viewed 320 times

Here's an album of my other resin art, mostly Penrose tiles, but you can see some gliders in one picture.

[pcallahan]

I wanted to see how this came out with a glossy lens-like surface. First, here's the silicone mold after coating with more silicone. You can see how the leveling and the meniscus turns all the cells into concave mirrors reflecting the overhead light.

20260613_073732.jpg (1.78 MiB) Viewed 250 times

Cast in clear resin, they turn into lenses, and you can even see the white lawn chair behind it, occluded by the object itself but visible in many cells.

20260614_063024 (1).jpg (3.46 MiB) Viewed 250 times

[Citation needed]

I noticed that a spammer (possibly DeepSeek) has created a copycat of this thread.

[NNlk05]

I noticed that a spammer (possibly DeepSeek) has created a copycat of this thread.

I don’t think they is a LLM/s(p|c)amer.

[pcallahan]

Latest Penrose variation: 3D-printed Wieringa roofs. https://www.printables.com/model/183098 ... ose-tiling

These came out well. I'm surprised I have not gotten more downloads on Printables. There are four pieces consisting of assemblies of golden rhombuses with appropriate tilts, and they snap onto kite and dart jigsaw pieces like this:

20260902_075228 (1).jpg (4.59 MiB) Viewed 99 times

Here is a larger section from two different angles.

20260906_125405.jpg (3.18 MiB) Viewed 99 times

20260906_125328.jpg (3.8 MiB) Viewed 99 times

I'm working on a different approach using angled connectors between flat rhombus faces. The only dihedral angles between faces are 72º, 108º, and 144º.

[pcallahan]

Using printed clips angled 72, 108, 144 degrees and old business cards.

20260912_115157.jpg (2.92 MiB) Viewed 66 times

You can see it roughly casts the expected shadow.

20260912_115659.jpg (2.5 MiB) Viewed 66 times

[pcallahan]

Here's a link to Printables for these clips: https://www.printables.com/model/184338 ... inga-roofs

I did some redesign and added printable panels as well as a guide for cutting business cards.

The pins are 72º, 108º, and 144º, which is the complete set of dihedral angles needed. I printed them in two pieces to avoid having to print a small gap reliably.

20260915_074043.jpg (1.34 MiB) Viewed 22 times

These pins are smaller than before but can be doubled up, positioned closer to corners.

20260914_223040.jpg (2.31 MiB) Viewed 22 times

You can also print panels with notches. These hold with a single pin per edge.

20260913_132418.jpg (2.66 MiB) Viewed 22 times

I wonder if I would be better off printing connectors for each corner. I haven't counted all the cases.