Source code: Tiling a large 3D-printed terrain model: jigsaw joints and dowel pins in Python

"""Tile a large 3D-printed terrain model: the central Lake District as six interlocking jigsaw tiles.

A 45 x 30 km Copernicus GLO-30 DEM is scaled once for the whole model (so seams match), the plan is
cut into a 3 x 2 grid of 150 mm tiles with jigsaw knobs and a printing clearance, and every tile is
proved watertight and non-overlapping. A second set of tiles uses straight edges with holes for
1.75 mm filament dowels instead.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/tile-large-3d-print-terrain/

Data: Copernicus DEM GLO-30 (© DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018,
provided under COPERNICUS by the European Union and ESA).
"""
import sys
from importlib.metadata import version
from pathlib import Path

import numpy as np
import rasterio
from rasterio.transform import from_origin
from rasterio.warp import reproject, Resampling
from shapely.geometry import Point, box
from shapely.geometry.polygon import orient
from shapely.ops import unary_union
import manifold3d as mf
import pyvista as pv
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9,
                     "figure.facecolor": "white", "axes.facecolor": "white"})
pv.OFF_SCREEN = True
OUT = Path("output"); OUT.mkdir(exist_ok=True)

DEM_URLS = [f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N54_00_W00{w}_00_DEM/"
            f"Copernicus_DSM_COG_10_N54_00_W00{w}_00_DEM.tif" for w in (4, 3)]
AOI = (310_000, 495_000, 355_000, 525_000)      # 45 x 30 km of the central Lake District, BNG
RES = 40.0
NX, NY = 3, 2                                    # tiles across and up
TILE_MM = 150.0                                  # printed size of each square tile
Z_EXAGGERATION = 2.0
BASE_MM = 5.0
CLEARANCE_MM = 0.2                               # total gap between mating edges
KNOB_NECK_MM, KNOB_R_MM, KNOB_CENTRE_MM = 6.0, 4.0, 6.0
PIN_HOLE_D_MM, PIN_DEPTH_MM = 1.9, 10.0          # for 1.75 mm filament used as dowels

X0, Y0, X1, Y1 = AOI
COLS, ROWS = int((X1 - X0) / RES) + 1, int((Y1 - Y0) / RES) + 1   # vertices on cell centres
DST = from_origin(X0 - RES / 2, Y1 + RES / 2, RES, RES)


def section(title):
    print(f"\n=== {title} ===")


# ------------------------------------------------------------------ mesh helpers (as in Part 1)
def grid_triangles(r, c):
    idx = np.arange(r * c).reshape(r, c)
    a, b = idx[:-1, :-1].ravel(), idx[:-1, 1:].ravel()
    d, e = idx[1:, :-1].ravel(), idx[1:, 1:].ravel()
    return np.concatenate([np.column_stack([a, d, b]), np.column_stack([b, d, e])])


def heightfield_solid(Z, cell_mm, x0=0.0, y0=0.0):
    rows, cols = Z.shape
    X, Y = np.meshgrid(x0 + np.arange(cols) * cell_mm, y0 + (rows - 1 - np.arange(rows)) * cell_mm)
    top = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    bot = top.copy()
    bot[:, 2] = 0.0
    n = rows * cols
    tri_top = grid_triangles(rows, cols)
    idx = np.arange(n).reshape(rows, cols)
    ring = np.concatenate([idx[0, :], idx[1:, -1], idx[-1, -2::-1], idx[-2:0:-1, 0]])
    p, q = ring, np.roll(ring, -1)
    walls = np.concatenate([np.column_stack([p, q, p + n]), np.column_stack([q, q + n, p + n])])
    return np.vstack([top, bot]), np.vstack([tri_top, tri_top[:, ::-1] + n, walls])


def to_manifold(verts, tris):
    m = mf.Manifold(mf.Mesh(vert_properties=np.ascontiguousarray(verts, dtype=np.float32),
                            tri_verts=np.ascontiguousarray(tris, dtype=np.uint32)))
    assert m.status() == mf.Error.NoError, m.status()
    return m


def from_manifold(m):
    mesh = m.to_mesh()
    return np.asarray(mesh.vert_properties)[:, :3].astype("float64"), np.asarray(mesh.tri_verts).astype("int64")


def check_mesh(verts, tris):
    e = np.concatenate([tris[:, [0, 1]], tris[:, [1, 2]], tris[:, [2, 0]]])
    n = len(verts)
    und = np.unique(np.sort(e, axis=1) @ np.array([n, 1]), return_counts=True)[1]
    drc = np.unique(e @ np.array([n, 1]), return_counts=True)[1]
    t = verts[tris]
    vol = np.einsum("ij,ij->i", t[:, 0], np.cross(t[:, 1], t[:, 2])).sum() / 6
    return bool((und == 2).all()), bool((drc == 1).all()), vol


def write_stl(path, verts, tris, header):
    t = verts[tris]
    nrm = np.cross(t[:, 1] - t[:, 0], t[:, 2] - t[:, 0])
    ln = np.linalg.norm(nrm, axis=1, keepdims=True)
    rec = np.zeros(len(tris), dtype=[("normal", "<f4", 3), ("v", "<f4", (3, 3)), ("attr", "<u2")])
    rec["normal"] = np.divide(nrm, ln, out=np.zeros_like(nrm), where=ln > 0)
    rec["v"] = t
    with open(path, "wb") as f:
        f.write(header.encode("ascii")[:80].ljust(80, b" "))
        f.write(np.uint32(len(rec)).tobytes())
        f.write(rec.tobytes())


def prism(poly, height):
    polys = getattr(poly, "geoms", [poly])
    contours = []
    for p in polys:
        p = orient(p, 1.0)
        contours += [np.array(p.exterior.coords)[:-1]] + [np.array(h.coords)[:-1] for h in p.interiors]
    return mf.Manifold.extrude(mf.CrossSection(contours), height + 10).translate((0, 0, -5))


section("versions")
print(f"python {sys.version.split()[0]} | numpy {np.__version__} | rasterio {rasterio.__version__} | "
      f"shapely {version('shapely')} | manifold3d {version('manifold3d')} | pyvista {pv.__version__}")

# ------------------------------------------------------------------ Step 1: one DEM from two source tiles
section("Step 1: mosaic two GLO-30 tiles onto a 40 m grid")
dem = np.full((ROWS, COLS), np.nan, dtype="float32")
with rasterio.Env(GDAL_HTTP_MAX_RETRY="4", GDAL_HTTP_RETRY_DELAY="5"):
    for url in DEM_URLS:
        with rasterio.open(url) as src:
            reproject(rasterio.band(src, 1), dem, dst_transform=DST, dst_crs="EPSG:27700",
                      resampling=Resampling.bilinear, dst_nodata=np.nan, init_dest_nodata=False)
        print(f"{url.rsplit('/', 1)[1]:<44s} NaNs left: {int(np.isnan(dem).sum()):,}")
print(f"grid {COLS} x {ROWS} vertices at {RES:.0f} m | {np.nanmin(dem):.1f} - {np.nanmax(dem):.1f} m")

# ------------------------------------------------------------------ Step 2: one scale for every tile
section("Step 2: one horizontal and vertical scale for the whole model")
scale = NX * TILE_MM / (X1 - X0)                 # mm per metre
cell_mm = RES * scale
Z = BASE_MM + (dem - dem.min()) * scale * Z_EXAGGERATION
W_mm, D_mm = (COLS - 1) * cell_mm, (ROWS - 1) * cell_mm
print(f"scale 1:{1000 / scale:,.0f} | cell {cell_mm:.2f} mm | whole model {W_mm:.0f} x {D_mm:.0f} x {Z.max():.1f} mm")
print(f"{NX} x {NY} tiles of {TILE_MM:.0f} x {TILE_MM:.0f} mm = {TILE_MM / scale / 1000:.0f} x "
      f"{TILE_MM / scale / 1000:.0f} km each")
assert abs(D_mm - NY * TILE_MM) < 1e-6

# ------------------------------------------------------------------ Step 3: jigsaw outlines
section("Step 3: interlocking outlines with a printing clearance")


def knob(x, y, direction):
    """A jigsaw knob rooted on a tile edge at (x, y), pointing 'E' or 'N'."""
    dx, dy = (1, 0) if direction == "E" else (0, 1)
    cx, cy = x + dx * KNOB_CENTRE_MM, y + dy * KNOB_CENTRE_MM
    half = KNOB_NECK_MM / 2
    if direction == "E":
        neck = box(x - 0.01, y - half, cx, y + half)
    else:
        neck = box(x - half, y - 0.01, x + half, cy)
    return neck.union(Point(cx, cy).buffer(KNOB_R_MM, quad_segs=16))


squares = {(i, j): box(i * TILE_MM, j * TILE_MM, (i + 1) * TILE_MM, (j + 1) * TILE_MM)
           for i in range(NX) for j in range(NY)}
exact = dict(squares)
for (i, j) in squares:
    for frac in (1 / 3, 2 / 3):
        if i < NX - 1:                                     # knobs from this tile into its east neighbour
            k = knob((i + 1) * TILE_MM, (j + frac) * TILE_MM, "E")
            exact[(i, j)] = exact[(i, j)].union(k)
            exact[(i + 1, j)] = exact[(i + 1, j)].difference(k)
        if j < NY - 1:                                     # and into its north neighbour
            k = knob((i + frac) * TILE_MM, (j + 1) * TILE_MM, "N")
            exact[(i, j)] = exact[(i, j)].union(k)
            exact[(i, j + 1)] = exact[(i, j + 1)].difference(k)
whole = box(0, 0, W_mm, D_mm)
pair_overlap = sum(exact[a].intersection(exact[b]).area for a in exact for b in exact if a < b)
print(f"exact partition: sum of tile areas {sum(p.area for p in exact.values()):,.3f} mm2 | "
      f"whole model {whole.area:,.3f} mm2 | overlap between tiles {pair_overlap:.6f} mm2")
# shrink every outline by half the clearance, so that each joint has the full clearance between parts
outlines = {k: p.buffer(-CLEARANCE_MM / 2, join_style="mitre") for k, p in exact.items()}
gaps = [outlines[a].distance(outlines[b]) for a in outlines for b in outlines
        if a < b and exact[a].touches(exact[b])]
print(f"clearance between neighbouring tiles: min {min(gaps):.3f} mm, max {max(gaps):.3f} mm "
      f"over {len(gaps)} joints")
name = {(i, j): f"{'ABC'[i]}{j + 1}" for i in range(NX) for j in range(NY)}

# ------------------------------------------------------------------ Step 4: cut the tiles
section("Step 4: cut six tiles and prove they fit")


def tile_solid(poly):
    """Heightfield over the outline's bounding box (plus a cell), cut to the outline."""
    minx, miny, maxx, maxy = poly.bounds
    c0, c1 = max(int(minx / cell_mm) - 1, 0), min(int(np.ceil(maxx / cell_mm)) + 2, COLS)
    r0 = max(ROWS - 1 - int(np.ceil(maxy / cell_mm)) - 1, 0)
    r1 = min(ROWS - int(miny / cell_mm) + 1, ROWS)
    hf = to_manifold(*heightfield_solid(Z[r0:r1, c0:c1], cell_mm, c0 * cell_mm, (ROWS - r1) * cell_mm))
    return hf ^ prism(poly, Z.max())


tiles = {k: tile_solid(p) for k, p in outlines.items()}
whole_terrain = to_manifold(*heightfield_solid(Z, cell_mm))
for k, m in tiles.items():
    v, t = from_manifold(m)
    wt, ori, vol = check_mesh(v, t)
    assert wt and ori and vol > 0
    lo, hi = v.min(axis=0), v.max(axis=0)
    print(f"tile {name[k]}: {len(t):>7,} triangles | watertight {wt} | oriented {ori} | {vol / 1000:6.1f} cm3 | "
          f"footprint {hi[0] - lo[0]:.1f} x {hi[1] - lo[1]:.1f} mm, height {hi[2]:.1f} mm")
    write_stl(OUT / f"lake_district_jigsaw_{name[k]}.stl", v, t,
              f"TechnicalGIS Lake District tile {name[k]} - Copernicus GLO-30")
keys = list(tiles)
overlap = sum((tiles[a] ^ tiles[b]).volume() for i, a in enumerate(keys) for b in keys[i + 1:])
lost = whole_terrain.volume() - sum(m.volume() for m in tiles.values())
print(f"pairwise overlap between tiles: {overlap:.4f} mm3")
print(f"whole model {whole_terrain.volume() / 1000:.2f} cm3 | tiles {sum(m.volume() for m in tiles.values()) / 1000:.2f} cm3 | "
      f"removed by the clearance {lost / 1000:.2f} cm3")

# ------------------------------------------------------------------ Step 5: the dowel-pin alternative
section("Step 5: straight edges with holes for filament dowels")


def pin_hole(x, y, axis):
    """A horizontal hole, PIN_DEPTH_MM either side of a tile edge, half way up the base."""
    cyl = mf.Manifold.cylinder(2 * PIN_DEPTH_MM, PIN_HOLE_D_MM / 2, circular_segments=24, center=True)
    cyl = cyl.rotate((0, 90, 0)) if axis == "x" else cyl.rotate((90, 0, 0))
    return cyl.translate((x, y, BASE_MM / 2))


holes, centres = [], []
for (i, j) in squares:
    for frac in (1 / 3, 2 / 3):
        if i < NX - 1:
            holes.append(pin_hole((i + 1) * TILE_MM, (j + frac) * TILE_MM, "x"))
            centres.append(((i + 1) * TILE_MM, (j + frac) * TILE_MM))
        if j < NY - 1:
            holes.append(pin_hole((i + frac) * TILE_MM, (j + 1) * TILE_MM, "y"))
            centres.append(((i + frac) * TILE_MM, (j + 1) * TILE_MM))
all_holes = mf.Manifold.batch_boolean(holes, mf.OpType.Add)
pinned = {k: tile_solid(sq.buffer(-CLEARANCE_MM / 2, join_style="mitre")) - all_holes for k, sq in squares.items()}
for k, m in pinned.items():
    v, t = from_manifold(m)
    wt, ori, vol = check_mesh(v, t)
    assert wt and ori and vol > 0
    write_stl(OUT / f"lake_district_pinned_{name[k]}.stl", v, t,
              f"TechnicalGIS Lake District pinned tile {name[k]} - Copernicus GLO-30")
print(f"{len(holes)} dowels, each {2 * PIN_DEPTH_MM - CLEARANCE_MM:.1f} mm long; every pinned tile watertight "
      f"and oriented: True")
print("pin holes per tile: " + ", ".join(
    f"{name[k]} {sum(sq.boundary.distance(Point(c)) < 1e-6 for c in centres)}" for k, sq in squares.items()))
wall = (BASE_MM - PIN_HOLE_D_MM) / 2
print(f"plastic above and below each hole: {wall:.2f} mm")

# ------------------------------------------------------------------ Figures
section("figures")
km = 1 / scale / 1000


def to_km(x, y):
    return X0 / 1000 + np.asarray(x) * km, Y0 / 1000 + np.asarray(y) * km


ext = (X0 / 1000 - RES / 2000, X1 / 1000 + RES / 2000, Y0 / 1000 - RES / 2000, Y1 / 1000 + RES / 2000)
hs = LightSource(azdeg=315, altdeg=45).hillshade(dem, vert_exag=1.5, dx=RES, dy=RES)
fig, ax = plt.subplots(figsize=(10, 6.8))
im = ax.imshow(dem, cmap="Spectral_r", extent=ext)
ax.imshow(hs, cmap="gray", extent=ext, alpha=0.35)
for k, p in exact.items():
    xs, ys = to_km(*p.exterior.xy)
    ax.plot(xs, ys, color="black", lw=1.0)
    cx, cy = to_km(*np.array(squares[k].centroid.coords[0]))
    ax.text(cx, cy, name[k], ha="center", va="center", fontsize=16, weight="bold", color="white",
            bbox=dict(boxstyle="round", fc="black", alpha=0.45, lw=0))
for label, e, n in [("Scafell Pike", 321_540, 507_210), ("Helvellyn", 334_210, 515_140)]:
    ax.plot(e / 1000, n / 1000, "k^", ms=6)
    ax.annotate(label, (e / 1000, n / 1000), xytext=(5, 5), textcoords="offset points", fontsize=8)
plt.colorbar(im, ax=ax, shrink=0.7, label="Elevation (m)")
ax.set_xlabel("Easting (km, British National Grid)")
ax.set_ylabel("Northing (km)")
ax.set_title(f"Copernicus GLO-30, 40 m: {NX} x {NY} tiles of {TILE_MM / scale / 1000:.0f} km "
             f"({TILE_MM:.0f} mm printed)")
fig.savefig(OUT / "01-tile-layout.png"); plt.close(fig)

# a close-up of one joint, showing the clearance
fig, axes = plt.subplots(1, 2, figsize=(10, 4.6))
jx, jy = TILE_MM, TILE_MM / 3
for ax, (half_w, title) in zip(axes, [(14, "One jigsaw knob (tile A1 into B1)"), (1.2, "Zoomed to the 0.2 mm clearance")]):
    for k, col in [((0, 0), "#3288bd"), ((1, 0), "#f46d43")]:
        p = outlines[k]
        for g in getattr(p, "geoms", [p]):
            ax.fill(*g.exterior.xy, color=col, alpha=0.35, lw=0)
            ax.plot(*g.exterior.xy, color=col, lw=1.2)
        ax.plot(*exact[k].exterior.xy, color="black", lw=0.6, ls=":")
    cx = jx + (KNOB_CENTRE_MM + KNOB_R_MM if half_w < 2 else KNOB_CENTRE_MM)
    ax.set_xlim(cx - half_w, cx + half_w)
    ax.set_ylim(jy - half_w * 0.8, jy + half_w * 0.8)
    ax.set_aspect("equal")
    ax.set_title(title)
    ax.set_xlabel("x on the print (mm)")
axes[0].set_ylabel("y on the print (mm)")
axes[0].legend(handles=[plt.Line2D([], [], color="#3288bd", lw=2, label="tile A1"),
                        plt.Line2D([], [], color="#f46d43", lw=2, label="tile B1"),
                        plt.Line2D([], [], color="black", lw=0.8, ls=":", label="exact cut line")],
               fontsize=8, loc="lower left")
fig.savefig(OUT / "02-joint-clearance.png"); plt.close(fig)


def save_cropped(p, png, pad=20):
    img = p.screenshot(return_img=True)
    ink = np.where((img < 250).any(axis=2))
    r0, r1 = max(ink[0].min() - pad, 0), min(ink[0].max() + pad, img.shape[0])
    c0, c1 = max(ink[1].min() - pad, 0), min(ink[1].max() + pad, img.shape[1])
    plt.imsave(OUT / png, img[r0:r1, c0:c1])
    p.close()


def look(p, centre, dist, elev_deg, azim_deg):
    el, az = np.radians(elev_deg), np.radians(azim_deg)
    c = np.asarray(centre, dtype=float)
    p.camera_position = [tuple(c + dist * np.array([np.cos(el) * np.cos(az), np.cos(el) * np.sin(az), np.sin(el)])),
                         tuple(c), (0, 0, 1)]
    p.camera.view_angle = 30


clim = (float(Z.min()), float(Z.max()))
p = pv.Plotter(off_screen=True, window_size=(2000, 1300))
p.set_background("white")
for k in tiles:
    mesh = pv.read(OUT / f"lake_district_jigsaw_{name[k]}.stl").translate((k[0] * 8.0, k[1] * 8.0, 0))
    p.add_mesh(mesh, scalars=mesh.points[:, 2], cmap="Spectral_r", clim=clim, show_scalar_bar=False,
               ambient=0.2, diffuse=0.8)
p.enable_lightkit()
look(p, (W_mm / 2 + 8, D_mm / 2 + 4, 5), 1.75 * W_mm, 38, -70)
save_cropped(p, "03-stl-preview-tiles.png")

fig_imgs = []
p = pv.Plotter(off_screen=True, window_size=(1800, 900), shape=(1, 2), border=False)
p.subplot(0, 0)
p.set_background("white")
for k, off, colour in [((0, 0), -4.0, "#d9c8a9"), ((1, 0), 4.0, "#9ecae1")]:
    mesh = pv.read(OUT / f"lake_district_jigsaw_{name[k]}.stl")
    p.add_mesh(mesh.clip_box((TILE_MM - 30, TILE_MM + 30, 25, 75, -1, 100), invert=False).translate((off, 0, 0)),
               color=colour, ambient=0.25, diffuse=0.75)
p.add_text("Jigsaw knob, tiles A1 and B1 pulled 8 mm apart", font_size=12, color="#222222")
p.enable_lightkit()
look(p, (TILE_MM, 50, 6), 125, 58, -95)
p.subplot(0, 1)
p.set_background("white")
mesh = pv.read(OUT / f"lake_district_pinned_{name[(1, 0)]}.stl")
p.add_mesh(mesh.clip_box((TILE_MM - 1, TILE_MM + 20, 82, 118, -1, 100), invert=False),
           color="#9ecae1", ambient=0.25, diffuse=0.75)
p.add_text("Dowel hole, west edge of pinned tile B1", font_size=12, color="#222222")
p.enable_lightkit()
look(p, (TILE_MM, 100, 8), 80, 16, 205)
save_cropped(p, "04-joint-closeups.png")
print("wrote 01-tile-layout, 02-joint-clearance, 03-stl-preview-tiles, 04-joint-closeups")