Painting → Bas-Relief

Building a depth-derived relief pipeline for Bruegel's Hunters in the Snow
Target: Bambu Lab A1 mini, single filament, 170 × 121 × 7 mm panel.

About this project

The goal. Turn a painting into a bas-relief: a printed panel on which the painting’s composition reads as physical depth — near things standing proud, distant things receding, the drawing carved into the surface. The subject is Bruegel’s Hunters in the Snow (1565), from the Google Art Project scan; the output is a 170 × 121 mm panel from a desktop FDM printer.

The approach, in plain words. A monocular depth model looks at the painting and estimates how far away everything is. That estimate is then conditioned into a printable height map: edges are snapped onto the painted edges, the depth range is compressed so the middle distance survives, the ground is subtracted so objects stand at full height, and the dark linework is cut into the surface as engraving. The height map becomes a watertight mesh, the mesh goes to the slicer, and the printer does the rest. The whole pipeline is six short Python files you type in yourself.

The form. Part guide, part lab diary. §1–11 are the recipe: the environment, the code file by file, and the print settings that produced the first panel — you can follow them start to finish and stop there. §12–13 are the revisions that came after, kept in the order they happened: what was changed, what the printed plastic showed, which ideas failed, and what each failure taught — the failed branches are documented as fully as the working one, because that is where most of the transferable lessons live (above all: the printer is a physical low-pass filter, and a screen preview will lie to you about what it can make). Decisions and their reasons are recorded at the point they were taken, so you can disagree with one and branch off from there.

0. The shape of the thing

painting.jpg
   │
   ├─▶ [A] depth model          Depth Anything V2  →  relative disparity map
   │                            (near = large value)
   ├─▶ [B] tiled refinement     re-run on overlapping crops, affine-align,
   │                            keep only the high frequencies
   ├─▶ [C] conditioning         guided filter → histogram equalise →
   │                            gradient-domain range compression →
   │                            luminance detail → rim fade
   ├─▶ [D] resample             to the print grid (0.25 mm pitch)
   ├─▶ [E] mesh                 heightfield → watertight solid
   └─▶ relief.stl               → Bambu Studio → A1 mini

The one idea that matters. A depth map is not a heightfield. Raw depth is dominated by a handful of huge jumps — the hunters are 20 m away, the mountains are 5 km away. Map that linearly onto 5 mm of travel and the entire village, the skaters, the birds, the whole middle of the painting occupy about 0.2 mm. You get a foreground silhouette glued to a flat wall.

Step [C] is where the print is won. Everything else is plumbing.

1. Why this painting in particular

Three things make Hunters in the Snow a great subject and a hard one:

It has real recession. Bruegel built the composition as a depth ladder — the hunters and dogs in the immediate foreground, the diagonal of trees, the frozen ponds at mid-distance, the village, then those impossible Alpine peaks. A depth model has something genuine to find. Compare this to a portrait, where you get a face-shaped lump and little else.

Luminance lies about depth here, badly. The snow is the brightest thing in the painting and it is simultaneously the nearest thing (foreground bank) and the furthest (distant fields). The trees are the darkest thing and they sit at mid-distance. Any pure brightness-to-height mapping produces nonsense on this image. This is exactly why you want the depth model rather than the lithophane approach.

Its detail is high-frequency and thin. Bare branches, birds, tiny skaters. A depth model at its native 518 px working resolution will smear all of that. Hence step [B], and hence the luminance detail term in step [C].

2. Environment

Hardware

On Apple Silicon, PyTorch uses the mps backend — the machine's GPU cores. No CUDA, no Rosetta, nothing exotic.

2.1 Python

macOS ships a Python you should not use. Get a clean one:

brew install python@3.12

If you don't have Homebrew: /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

2.2 Project folder and virtualenv

The scan's filename has parentheses in it, and your own path may have spaces — quote paths everywhere. The examples use /Users/you/painting2model; put the folder wherever you like and substitute your own path.

cd "/Users/you/painting2model"
/opt/homebrew/bin/python3.12 -m venv .venv
source .venv/bin/activate
python -V          # expect 3.12.x

If that last line reports anything but 3.12. The venv is almost certainly fine — what usually breaks is the name python, not the environment. Check which:

type -a python python3
python3 -V

python3 -V reads the venv through PATH and is hard to fool, so if it says 3.12.x your environment is correct and only the lookup of python is wrong. The two things that shadow it:

• An alias. A line like alias python=/usr/bin/python3 in your ~/.zshrc beats PATH outright, so python silently means Apple's system 3.9 no matter which venv is active. type -a python prints it first when this is the case. Delete the alias, or repoint it at python3 so it resolves through PATH like everything else — grep -n "alias python" ~/.zshrc ~/.zprofile ~/.zshenv finds it. This matters beyond this line: every python in §§2.4, 2.6, 6.5 and 9 would otherwise run the system interpreter, which has none of your packages.
• The wrong activate script. source .venv/bin/activate is bash/zsh syntax. Under fish or csh use activate.fish or activate.csh — both are already in your .venv/bin.

If python3 -V is also wrong, then the venv really was built by the wrong interpreter. cat .venv/pyvenv.cfg names the one that made it; delete the folder and re-run the block above with the full Homebrew path.

2.3 Dependencies

pip install --upgrade pip
pip install torch torchvision
pip install "transformers>=4.45" pillow numpy scipy trimesh

Note what is not here: no opencv, no open3d, no numpy-stl. The guided filter is fifteen lines of scipy.ndimage (§6.1) and trimesh writes STL directly. Fewer moving parts, fewer wheel-compatibility problems on arm64.

2.4 Confirm the GPU is live

python - <<'PY'
import torch
print("torch", torch.__version__)
print("mps available:", torch.backends.mps.is_available())
print("mps built:", torch.backends.mps.is_built())
PY

Both should be True. If mps available is False, the pipeline still runs on CPU — roughly 8–15× slower, which for a handful of inferences means minutes rather than seconds. Not fatal.

Add this to your shell profile or prepend it to runs; a few ops in the DPT head occasionally lack MPS kernels and this lets them silently fall back to CPU instead of crashing:

export PYTORCH_ENABLE_MPS_FALLBACK=1

2.5 Which model size

INFER_LONG_SIDE must be a multiple of 14 (the ViT patch size): 518, 728, 1036, 1428. The model was trained at 518; going higher recovers real detail but with diminishing returns and quadratic memory cost. First run downloads weights (~1.3 GB for Large) into ~/.cache/huggingface. Neither checkpoint is gated — no Hugging Face account or token is needed. Both are released under CC-BY-NC-4.0, so a panel on your own wall is fine and selling prints of one is not.

One coupling to know about: main.py feeds the model a working image at twice the print grid, i.e. 1360 × 968 px. Setting INFER_LONG_SIDE above 1360 would just be upsampling — if you want to go there, raise the grid_w * 2 factor in main.py to * 3 first.

2.6 Record your versions

Two kinds of number appear in this guide. The maths — the Poisson round-trip, the guided filter's behaviour, the compression table — is pure numpy/scipy, and §6.5 lets you re-derive all of it on your own machine in about a second. The model's behaviour is not: it depends on the exact transformers and torch you installed, and the DPT processor's handling of size and do_rescale has moved between releases before. So record what you have, once:

mkdir -p out
python - <<'PY' | tee out/versions.txt
import platform, torch, transformers, numpy, scipy, PIL, trimesh
print("macos       ", platform.mac_ver()[0], platform.machine())
print("python      ", platform.python_version())
print("torch       ", torch.__version__)
print("transformers", transformers.__version__)
print("numpy       ", numpy.__version__)
print("scipy       ", scipy.__version__)
print("pillow      ", PIL.__version__)
print("trimesh     ", trimesh.__version__)
PY

If a reinstall later changes what the model returns, that file tells you what moved, and you can pin from it — pip install torch==<the version in your file> — rather than from memory. Nothing in §§4–8 pins a version, because a pin invented today ages worse than a recorded fact.

2.7 The source image

Nothing so far has put a painting in the folder. Get the highest-resolution scan you can — for this one, the Google Art Project file on Wikimedia Commons, which is public domain (the painting is from 1565) and about 8 MB at 6819 × 4853. Download the full-size original, not a thumbnail: every pixel you skip here is detail the model cannot recover later. Put it in the project folder and check what you actually got:

cd "/Users/you/painting2model"
python - <<'PY'
from PIL import Image
from pathlib import Path
Image.MAX_IMAGE_PIXELS = None
for p in sorted(Path(".").glob("*.jpg")) + sorted(Path(".").glob("*.png")):
    print(f"{p.name}   {Image.open(p).size}")
PY

Write down the filename and the size, because config.py in §4 needs the exact path and the pixel dimensions decide the panel's proportions. The listing in §4 uses 6819 × 4853 — aspect 1.4051 — and every derived number in this guide follows from it. If your file differs, your grid and triangle counts will differ by the same ratio; nothing is wrong, the numbers just won't match line for line.

A different painting works too — see §14 — but then §1's argument about why this one suits a depth model is yours to re-make.

3. Project layout

painting2model/
├── .venv/
├── config.py         all parameters, nothing else
├── heightfield.py    the maths — pure numpy/scipy, no torch
├── depth.py          the model — torch only lives here
├── mesh.py           heightfield → STL
├── main.py           orchestration + checkpoint images
├── out/              generated
└── Pieter_Bruegel_the_Elder_-_Hunters_in_the_Snow_(Winter)_-_Google_Art_Project.jpg

Keeping heightfield.py free of torch is deliberate: it means you can iterate on the conditioning — which is where you'll spend your time — without reloading a 1.3 GB model. That split only pays off if the depth map itself is cached between runs, which is what USE_CACHED_DEPTH does (§8): the first run writes out/depth_cache.npz, and every later run that changes only conditioning parameters reloads it and skips the model entirely. Seconds instead of minutes, which is the difference between tuning §6 and not bothering.

4. config.py

Write this first. Every parameter you'll want to turn is here, and nowhere else.

The file is one continuous listing; it is broken into its five sections here so each group of knobs can be read on its own. Type them in this order into a single config.py.

"""All tunable parameters. Edit here, never in the other modules."""
from pathlib import Path

# ---- input / output -------------------------------------------------------
IMAGE = Path("/Users/you/painting2model/"
             "Pieter_Bruegel_the_Elder_-_Hunters_in_the_Snow_(Winter)_-_Google_Art_Project.jpg")
OUTDIR = IMAGE.parent / "out"

# Crop the scanned frame/margin away, as fractions of width/height.
# (left, top, right, bottom). Check out/00_source.png and adjust.
CROP = (0.0, 0.0, 0.0, 0.0)

Which checkpoint, at what resolution, and how the tiled second pass is shaped.

# ---- depth model ----------------------------------------------------------
MODEL_ID = "depth-anything/Depth-Anything-V2-Large-hf"   # Base-hf on 8 GB Macs
INFER_LONG_SIDE = 1036        # must be a multiple of 14. 518 / 728 / 1036 / 1428
USE_TILED_REFINE = True
TILE = 518                    # multiple of 14
TILE_OVERLAP = 172
TILE_LOWPASS_SIGMA_PX = 24.0  # in working-grid pixels
USE_CACHED_DEPTH = True       # reuse out/depth_cache.npz while tuning the block below

The seven numbers that decide what the relief actually looks like. This is the block you will come back to.

# ---- heightfield conditioning --------------------------------------------
GUIDED_RADIUS = 6             # px, snaps depth edges onto painted edges
GUIDED_EPS = 2e-3
EQUALIZE = 0.55               # 0 = keep model's depth distribution, 1 = flat histogram
COMPRESS_ALPHA = 4.0          # gradient-domain range compression. 0 = off, 2-6 typical
DETAIL_GAIN = 0.12            # brushwork/texture from image luminance. 0 = pure depth
DETAIL_SIGMA_PX = 2.5
INVERT = False                # True if the print comes out with foreground sunken
RIM_MM = 2.0                  # flat recessed border. 0 = relief runs to the edge

Millimetres and sample spacing. Everything downstream is derived from these four.

# ---- physical print -------------------------------------------------------
WIDTH_MM = 170.0              # A1 mini plate is 180x180; leave margin for a brim
BASE_MM = 2.0                 # solid backing plate
RELIEF_MM = 5.0               # peak-to-valley sculpted depth
PITCH_MM = 0.25               # sample spacing; 0.4 mm nozzle, so 0.25 is ~1.6x oversampled

# ---- preview --------------------------------------------------------------
LIGHT_AZIMUTH_DEG = 315.0
LIGHT_ALTITUDE_DEG = 32.0

Where those numbers come from

WIDTH_MM = 170. The A1 mini build volume is 180 × 180 × 180 mm. Going to 178 technically fits but leaves nothing for a brim and puts you inside the plate's exclusion margins. At 170 the panel is 170 × 121.0 mm for the scan in §2.7. Note that that is the scan's aspect, 1.4051, not the painting's — Hunters in the Snow is 162 × 117 cm, aspect 1.385, and the Google Art Project file is a little wider than the stretcher. The script reads the real aspect from your file rather than assuming, so a different scan shifts every number in this section slightly.

PITCH_MM = 0.25. Your nozzle is 0.4 mm, so 0.4 mm is the true XY feature floor. Sampling at 0.25 gives the slicer ~1.6 samples per nozzle width, which is enough to render smooth diagonals without wasting geometry. That's a 680 × 484 grid → 662,886 triangles → a 33.1 MB binary STL. Bambu Studio handles that comfortably. Dropping to 0.15 mm quadruples the file for detail the nozzle physically cannot lay down.

RELIEF_MM = 5.0 and the layer-height trap. This is the parameter people get wrong. Your real Z resolution is not infinite — it is the layer height. At 0.10 mm layers, 5 mm of relief gives you 50 distinct height levels, and that's it. Fifty grey levels, in effect. Shrink the relief to 2 mm and you have 20, and the print will visibly terrace. Go to 8 mm and you get 80 levels but the panel starts to read as a sculpture rather than a relief, and print time climbs. 5 mm at 0.10 mm layers is the sweet spot; if you want more, drop to 0.08 mm layers (62 levels) rather than increasing the depth.

BASE_MM = 2.0. Below about 1.5 mm a 170 mm panel will cup as it cools. 2 mm is 20 layers at 0.10 — stiff enough, and it gives you something to drill or glue a hanger onto.

5. depth.py — running the model

This file runs Depth Anything V2 on the painting. How such a model can read depth out of a single image at all — its architecture, training and metrics — is a subject of its own: How the depth model works, a primer that assumes no machine-learning background.

5.1 Why not pipeline()

The one-liner works:

pipe = pipeline("depth-estimation", model=MODEL_ID, device="mps")
result = pipe(image)

but its result["depth"] is an 8-bit PIL image. Eight bits, quantised, for a field you're about to compress in the gradient domain. You'd be throwing away most of your dynamic range before you started. result["predicted_depth"] is the float tensor, but it comes out at the model's internal resolution and the pipeline gives you no control over what that resolution is.

Use the processor and model directly. Slightly more code, complete control.

5.2 The polarity thing

Depth Anything outputs relative inverse depth — disparity. Larger values are NEARER. That happens to be exactly what you want for a relief (near things stick out), so no inversion is needed. But this is a convention, not a law, and it is the single most likely thing to silently flip on you if you swap models. The INVERT flag in config and checkpoint image 01_depth_global.png exist so you can verify it in two seconds rather than after a twelve-hour print.

Because this is the failure most likely to cost you a print, main.py also checks it for you: after the global pass it compares the mean depth of the bottom quarter of the frame against the top quarter and prints a verdict. That is a heuristic about receding ground planes, not a law — an aerial view, a ceiling fresco or a flat-on close-up will fail it legitimately — so it reports unclear rather than guessing whenever the margin is small. For Hunters in the Snow, with a foreground bank at the bottom and Alpine peaks at the top, it is decisive.

5.3 The file

Four short pieces, in order, make up depth.py.

"""Monocular depth inference on Apple Silicon (MPS)."""
import time
import numpy as np
import torch
import torch.nn.functional as F
from scipy.ndimage import gaussian_filter
from transformers import AutoImageProcessor, AutoModelForDepthEstimation

_CACHE = {}


def pick_device():
    return "mps" if torch.backends.mps.is_available() else "cpu"


def load(model_id, device):
    if model_id not in _CACHE:
        t = time.time()
        print(f"    fetching {model_id}", flush=True)
        proc = AutoImageProcessor.from_pretrained(model_id)
        model = AutoModelForDepthEstimation.from_pretrained(model_id)
        print(f"    weights in RAM after {time.time() - t:.0f}s; moving to {device}", flush=True)
        model = model.to(device).eval()
        _CACHE[model_id] = (proc, model)
        print(f"    model ready in {time.time() - t:.0f}s", flush=True)
    return _CACHE[model_id]


def _round14(n):
    return max(int(round(n / 14)) * 14, 14)

The processor is fed an explicit size so the model runs at a resolution you chose rather than one it picked. _round14 keeps it legal for DPT.

@torch.inference_mode()
def infer(pil_img, model_id, device, long_side):
    """Raw model output as float64 numpy, at the input image's own size.

    Depth Anything outputs relative inverse depth (disparity):
    LARGER VALUES ARE NEARER.
    """
    proc, model = load(model_id, device)
    w, h = pil_img.size
    if w >= h:
        size = {"width": _round14(long_side), "height": _round14(long_side * h / w)}
    else:
        size = {"height": _round14(long_side), "width": _round14(long_side * w / h)}

Output comes back at the model's internal resolution, so it is interpolated back up to the image's own size before leaving the torch world as float64.

    inputs = proc(images=pil_img, return_tensors="pt", size=size, do_rescale=True)
    inputs = {k: v.to(device) for k, v in inputs.items()}
    pred = model(**inputs).predicted_depth              # (1, h', w')
    if pred.dim() == 3:
        pred = pred.unsqueeze(1)
    pred = F.interpolate(pred.float(), size=(h, w), mode="bicubic", align_corners=False)
    return pred[0, 0].detach().to("cpu").numpy().astype(np.float64)

do_rescale=True and the explicit size dict override the processor's defaults; the DPT processor wants multiples of 14, which _round14 guarantees.

Expect roughly 10–40 s for one Large inference at 1036 on an M-series Mac, dominated by the first call (weight loading). A fanless machine will thermal-throttle partway through a long tiled run. That's normal and harmless.

5.4 Tiled refinement (append to depth.py)

One global pass sees the whole painting at ~1036 px and produces coherent but soft depth. Running the model again on overlapping 518 px crops lets it resolve the branches and the skaters — but each crop comes back with its own arbitrary scale and offset, because monocular depth is only relative. Two adjacent crops might report the same tree at 0.3 and at 1.7.

The fix is three moves:

1. Affine-align each crop to the global pass. Least-squares fit a·local + b ≈ global over the crop's footprint. Now every crop lives in the global pass's units.
2. Hann-window blend. Weight each crop by a 2D raised cosine so contributions taper to zero at crop edges. Overlapping crops cross-fade instead of butting.
3. Keep only the high frequencies from the merged result. The global pass is authoritative about large-scale structure; the crops are authoritative about fine structure. result = lowpass(global) + highpass(merged).

Measured on a synthetic scene with known ground truth: high-frequency energy went from 0.0216 (global only) to 0.0475 (refined) against a true value of 0.0530, and correlation with truth rose from 0.973 to 0.9986. Residual seam visibility at tile boundaries was 0.5–1.4× the mean gradient, i.e. lost in the noise.

Three small helpers, then the routine itself split at the point where the tile loop begins. All of it appends to depth.py.

def _hann2d(h, w):
    return np.maximum(np.outer(np.hanning(h + 2)[1:-1], np.hanning(w + 2)[1:-1]), 1e-6)

Fits a·local + b ≈ global, with a mean-offset fallback when the crop is degenerate (a flat sky patch, for instance).

def _fit_affine(local, ref):
    A = np.stack([local.ravel(), np.ones(local.size)], axis=1)
    a, b = np.linalg.lstsq(A, ref.ravel(), rcond=None)[0]
    if not np.isfinite(a) or abs(a) < 1e-8:
        return 1.0, float(ref.mean() - local.mean())
    return float(a), float(b)

Regular stride, plus a final flush-to-the-edge position so the last strip is never dropped.

def _positions(n, tile, overlap):
    step = max(tile - overlap, 1)
    pos = list(range(0, max(n - tile, 0) + 1, step))
    if not pos or pos[-1] != max(n - tile, 0):
        pos.append(max(n - tile, 0))
    return sorted(set(pos))

Accumulators for the weighted sum and the weights, one shared window, and the grid of crop origins.

def tiled_refine(pil_img, global_depth, model_id, device,
                 tile=518, overlap=172, lowpass_sigma=24.0, log=print):
    W, H = pil_img.size
    tile = min(tile, H, W)
    acc = np.zeros((H, W)); wacc = np.zeros((H, W))
    win = _hann2d(tile, tile)
    ys, xs = _positions(H, tile, overlap), _positions(W, tile, overlap)
    total, n = len(ys) * len(xs), 0

Each crop is inferred, affine-aligned, and accumulated under the window. The return line is the whole trick: low frequencies from the global pass, high frequencies from the merged crops.

    for y in ys:
        for x in xs:
            n += 1
            log(f"    tile {n}/{total}  at ({x}, {y})")
            crop = pil_img.crop((x, y, x + tile, y + tile))
            local = infer(crop, model_id, device, long_side=tile)
            a, b = _fit_affine(local, global_depth[y:y + tile, x:x + tile])
            acc[y:y + tile, x:x + tile] += (a * local + b) * win
            wacc[y:y + tile, x:x + tile] += win
    merged = acc / wacc
    return (gaussian_filter(global_depth, lowpass_sigma)
            + merged - gaussian_filter(merged, lowpass_sigma))

At the settings in config this is 12 tiles for Hunters in the Snow — call it 3–10 minutes. Set USE_TILED_REFINE = False while you're tuning the conditioning, turn it on for the final run.

6. heightfield.py — the part that matters

6.1 Guided filter

The depth model works at low resolution and upsamples, so its boundaries are soft and misaligned with the painting's actual edges: a tree trunk's depth edge sits three or four pixels off from the tree trunk. Printed, that reads as a halo.

The guided filter (He et al., 2010) fixes this by locally fitting depth as a linear function of image luminance within a sliding window, then applying that fit. Where the image has a crisp edge, the depth inherits it. Where the image is flat, it averages, killing noise.

Measured on a synthetic test: edge step magnitude went 0.180 → 0.470 (true edge 1.000) while noise in flat regions fell from 0.049 → 0.0017. Sharper and cleaner.

6.2 Range compression — the core of it

Take the depth gradient field, attenuate large gradients logarithmically, leave small ones untouched, then reintegrate by solving a Poisson equation. This is the Weyrich et al. (2007) digital-bas-relief construction and it is the difference between a print that reads as a sculpted scene and a print that reads as a step.

For small |g|, log(1+αg)/α ≈ g — fine detail passes through untouched. For large |g|, it grows logarithmically — a 100× gradient becomes a ~5× gradient. Big steps get crushed, texture survives.

Reintegration is a Poisson solve with Neumann boundaries, which has an exact closed-form solution in the DCT domain: transform the divergence, divide by the Laplacian eigenvalues 2(cos(πi/M) − 1) + 2(cos(πj/N) − 1), zero the DC term, inverse transform. One DCT pair. No iteration, no solver library. Round-trip error measured at 3.9 × 10⁻¹⁴ — exact to floating point.

If that sentence went past too quickly — what a Poisson solve is, why 2D needs one at all, and what “Neumann” buys you — there is a side note that builds the whole thing from arithmetic: The Poisson solve. It is self-contained and assumes no calculus.

Measured effect on a test scene with a hard foreground step plus fine texture:

α = 15 flattens the scene into pure texture — the depth story is gone. α = 2–6 is the usable band. Start at 4.

The same crop — bare branches over the valley — through the whole sweep. The middle three are from a real tuning session; the two ends are rendered from the same depth map to show where the band stops working:

6.3 The file

This is the longest file, so it is broken here into one block per function. They concatenate, in this order, into heightfield.py.

"""Turn a raw relative-depth map into a printable heightfield. No torch here."""
import numpy as np
from scipy.fft import dctn, idctn
from scipy.ndimage import uniform_filter, gaussian_filter, zoom

Normalise to [0,1] with a zero-range guard, and a resize that will not off-by-one you.

# ---------------------------------------------------------------- utilities
def norm01(a):
    a = np.asarray(a, dtype=np.float64)
    lo, hi = float(a.min()), float(a.max())
    return np.zeros_like(a) if hi - lo < 1e-12 else (a - lo) / (hi - lo)


def resize_to(a, shape, order=3):
    zy, zx = shape[0] / a.shape[0], shape[1] / a.shape[1]
    out = zoom(a, (zy, zx), order=order, mode="nearest")
    return out[:shape[0], :shape[1]]

The whole filter is six lines of box means: local variance of the guide, local covariance with the source, then the linear coefficients.

# ------------------------------------------------- edge-aware guided filter
def guided_filter(guide, src, radius=6, eps=2e-3):
    """He et al. 2010. Snaps `src` onto the edges of `guide` and kills noise."""
    box = lambda a: uniform_filter(a, size=2 * radius + 1, mode="nearest")
    I, p = np.asarray(guide, np.float64), np.asarray(src, np.float64)
    mI, mp = box(I), box(p)
    var = box(I * I) - mI * mI
    cov = box(I * p) - mI * mp
    a = cov / (var + eps)
    b = mp - a * mI
    return box(a) * I + box(b)

Divergence of the gradient field, then a DCT, a division by the Laplacian eigenvalues, and a DCT back. Setting lam[0,0] = 1 before dividing and zeroing the DC term afterwards is what pins the arbitrary constant.

# --------------------------------------- gradient-domain range compression
def _laplacian_eigenvalues(shape):
    M, N = shape
    ii = np.arange(M).reshape(-1, 1)
    jj = np.arange(N).reshape(1, -1)
    return 2 * (np.cos(np.pi * ii / M) - 1) + 2 * (np.cos(np.pi * jj / N) - 1)


def poisson_from_gradients(gx, gy):
    """Solve laplacian(h) = div(g), Neumann boundary, via DCT. Exact and fast."""
    div = np.zeros_like(gx)
    div[:, 1:] += gx[:, 1:] - gx[:, :-1]
    div[:, 0] += gx[:, 0]
    div[1:, :] += gy[1:, :] - gy[:-1, :]
    div[0, :] += gy[0, :]

    lam = _laplacian_eigenvalues(div.shape)
    lam[0, 0] = 1.0
    h_hat = dctn(div, type=2, norm="ortho") / lam
    h_hat[0, 0] = 0.0
    return idctn(h_hat, type=2, norm="ortho")

Gradients are normalised by their own mean magnitude first, so alpha means the same thing regardless of the depth map's scale.

def compress_range(h, alpha=4.0):
    """Logarithmically attenuate large gradients, leave small ones alone."""
    if alpha <= 0:
        return np.array(h, dtype=np.float64)
    gx = np.zeros_like(h); gx[:, :-1] = np.diff(h, axis=1)
    gy = np.zeros_like(h); gy[:-1, :] = np.diff(h, axis=0)
    mag = np.hypot(gx, gy)
    nz = mag > 1e-12
    if not np.any(nz):
        return np.array(h, dtype=np.float64)
    s = mag[nz].mean()
    scale = np.ones_like(mag)
    mn = mag[nz] / s
    scale[nz] = np.log1p(alpha * mn) / (alpha * mn)
    return poisson_from_gradients(gx * scale, gy * scale)

A rank transform, blended with the original by strength.

# ------------------------------------------------------------ tonal shaping
def equalize(h, strength=1.0):
    """Push the depth histogram toward flat so every depth band gets relief."""
    if strength <= 0:
        return np.array(h, dtype=np.float64)
    flat = np.asarray(h, np.float64).ravel()
    order = np.argsort(flat)
    ranks = np.empty(flat.size, dtype=np.float64)
    ranks[order] = np.linspace(0.0, 1.0, flat.size)
    return norm01((1 - strength) * norm01(flat) + strength * ranks).reshape(h.shape)

def add_detail(h, luma, gain=0.12, sigma=2.5):
    """Add the painting's high-frequency structure as fine surface texture."""
    if gain <= 0:
        return np.array(h, dtype=np.float64)
    hp = luma - gaussian_filter(luma, sigma)
    denom = np.abs(hp).max()
    if denom < 1e-9:
        return np.array(h, dtype=np.float64)
    return h + gain * (hp / denom)

Distance to the nearest edge, clipped to the rim width, put through a smoothstep so the border fades rather than steps.

def apply_rim(h, rim_mm, pitch_mm):
    """Fade the outermost few mm down to the base plate: a clean recessed border."""
    if rim_mm <= 0:
        return h
    r = max(int(round(rim_mm / pitch_mm)), 1)
    H, W = h.shape
    dy = np.minimum(np.arange(H), np.arange(H)[::-1]).reshape(-1, 1)
    dx = np.minimum(np.arange(W), np.arange(W)[::-1]).reshape(1, -1)
    t = np.clip(np.minimum(dy, dx) / r, 0.0, 1.0)
    return h * (t * t * (3 - 2 * t))          # smoothstep

Lambertian shading of the surface normals. Note the row/world sign flip on dzdy — get it wrong and the preview lights from the opposite side.

# ------------------------------------------------------------------ preview
def hillshade(h01, relief_mm, pitch_mm, azimuth_deg=315.0, altitude_deg=32.0):
    """What the print will look like under a raking lamp."""
    dz_drow, dzdx = np.gradient(h01 * relief_mm, pitch_mm)
    dzdy = -dz_drow                     # image rows run downward, world +y runs up
    nz = np.ones_like(dzdx)
    norm = np.sqrt(dzdx ** 2 + dzdy ** 2 + 1.0)
    az, alt = np.radians(azimuth_deg), np.radians(altitude_deg)
    lx, ly, lz = np.cos(alt) * np.sin(az), np.cos(alt) * np.cos(az), np.sin(alt)
    shade = (-dzdx * lx + -dzdy * ly + nz * lz) / norm
    return np.clip(shade, 0.0, 1.0)

6.4 On equalize and add_detail

EQUALIZE remaps the depth histogram toward flat, so that depth bands containing lots of pixels get proportionally more of your 5 mm. Bruegel's painting has a big population of pixels at mid-distance (the valley floor); without equalisation they share a narrow slab. At 1.0 it's aggressive and the scene starts to look pressure-flattened; 0.5–0.6 is a good landing.

add_detail is the one place luminance re-enters. You chose the depth panel, not the hybrid, so this defaults low — 0.12 puts about 0.6 mm of high-frequency brushwork on top of 5 mm of depth-driven form. It's what makes the bare branches and the roof tiles legible rather than implied. Set it to 0 for a pure depth interpretation, and be aware the result will be noticeably softer.

hillshade is your most valuable checkpoint. It renders the heightfield under a simulated raking light — Lambertian shading of the surface normals. What you see in 06_preview_lit.png is very close to what the physical panel will look like under a lamp. Judge the print from this image, not from the greyscale heightfield, which is misleading (your eye reads it as an image, not a surface).

6.5 Checking your typing

You have just hand-typed a hundred lines of numerical code in which one transposed index or one sign produces a plausible-looking but wrong relief — the kind of bug you would otherwise discover as a disappointing print. Every claim §6.1 and §6.2 make about this code is checkable in a second. Save this as selftest.py beside heightfield.py and run python selftest.py.

"""Checks heightfield.py against the claims in this guide. Run: python selftest.py"""
import numpy as np
from scipy.ndimage import gaussian_filter
import heightfield as hf

N = 256
yy, xx = np.mgrid[0:N, 0:N] / (N - 1.0)
EDGE = 102                                   # the step sits between column 101 and 102

# One hard step (a near slab on the left) with fine texture laid over it.
step = np.where(np.arange(N)[None, :] < EDGE, 0.90, 0.10) * np.ones((N, 1))
scene = step + 0.02 * np.sin(2 * np.pi * 16 * xx) * np.sin(2 * np.pi * 16 * yy)

# What a depth model hands you: the same edge, but soft, and noisy.
rng = np.random.default_rng(0)
guide = np.where(np.arange(N)[None, :] < EDGE, 0.85, 0.15) * np.ones((N, 1))
soft = gaussian_filter(step, 4.0) + 0.05 * rng.standard_normal((N, N))

d_edge = lambda a: float(a[:, EDGE - 1].mean() - a[:, EDGE].mean())
d_noise = lambda a: float(a[N // 4:N // 2, N // 8:N // 4].std())
d_texture = lambda a: float(a[N // 4:3 * N // 4, N // 8:N // 3].std())

print("[1] poisson round-trip")
gx = np.zeros_like(scene); gx[:, :-1] = np.diff(scene, axis=1)
gy = np.zeros_like(scene); gy[:-1, :] = np.diff(scene, axis=0)
back = hf.poisson_from_gradients(gx, gy)
err = np.abs((back - back.mean()) - (scene - scene.mean())).max()
print(f"    max abs error {err:.2e}")
assert err < 1e-9, "reintegration is not inverting the gradient operator"

print("[2] guided filter")
g = hf.guided_filter(guide, soft, radius=6, eps=2e-3)
print(f"    edge  {d_edge(soft):.3f} -> {d_edge(g):.3f}   (true {d_edge(step):.3f})")
print(f"    noise {d_noise(soft):.4f} -> {d_noise(g):.4f}")
assert d_edge(g) > 2.0 * d_edge(soft), "guided filter did not sharpen the edge"
assert d_noise(g) < 0.2 * d_noise(soft), "guided filter did not suppress noise"

print("[3] range compression")
print("    alpha    step   texture   ratio")
rows = []
for a in (0.0, 2.0, 4.0, 6.0, 15.0):
    c = hf.norm01(hf.compress_range(scene, a))
    s, t = d_edge(c), d_texture(c)
    rows.append((s, t))
    print(f"    {a:5.0f}   {s:.3f}   {t:.4f}   {t / s:.3f}")
assert all(rows[i][0] > rows[i + 1][0] for i in range(len(rows) - 1)), "step must shrink"
assert all(rows[i][1] < rows[i + 1][1] for i in range(len(rows) - 1)), "texture must grow"

print("[4] shaping")
e = hf.equalize(hf.norm01(scene), 1.0)
r = hf.apply_rim(np.ones((N, N)), rim_mm=2.0, pitch_mm=0.25)
sh = hf.hillshade(hf.norm01(scene), 5.0, 0.25)
print(f"    equalize(1.0) median {np.median(e):.3f}")
print(f"    rim border {r[0, 0]:.3f}  centre {r[N // 2, N // 2]:.3f}")
print(f"    hillshade in [{sh.min():.3f}, {sh.max():.3f}]")
assert abs(np.median(e) - 0.5) < 0.02, "equalize(1.0) should flatten the histogram"
assert r[0, 0] == 0.0 and r[N // 2, N // 2] == 1.0, "rim should fade edge to zero"
assert np.isfinite(sh).all() and sh.min() >= 0.0 and sh.max() <= 1.0

print("\nall checks passed")

It builds one synthetic scene — a hard step with fine texture over it, plus the softened, noisy version a depth model would hand you — and puts heightfield.py through it. Expected output, about one second:

[1] poisson round-trip
    max abs error 2.69e-13
[2] guided filter
    edge  0.080 -> 0.400   (true 0.800)
    noise 0.0510 -> 0.0019
[3] range compression
    alpha    step   texture   ratio
        0   0.952   0.0119   0.013
        2   0.459   0.1293   0.282
        4   0.383   0.1454   0.380
        6   0.346   0.1524   0.440
       15   0.285   0.1629   0.572
[4] shaping
    equalize(1.0) median 0.500
    rim border 0.000  centre 1.000
    hillshade in [0.000, 0.634]

all checks passed

Read it as four statements. The Poisson solve inverts the gradient operator exactly — 2.7 × 10⁻¹³ is floating-point noise, not approximation error, which is what earns the “no iteration, no solver library” claim in §6.2. The guided filter multiplies the edge step by five (0.080 → 0.400, against a true 0.800) while cutting flat-region noise by twenty-five (0.0510 → 0.0019) — sharper and cleaner, in one pass. The step column falls monotonically as the texture column rises, which is the entire argument of §6.2 reproduced on your machine; note the ratio climbing 0.013 → 0.572 as α goes 0 → 15, and that α = 15 has flattened the step to 0.285 — the depth story is gone, exactly as the table above warns. The shaping helpers do what their names say: a flat histogram, a border at zero, a hillshade that stays in range.

The absolute values belong to this scene and this seed, so treat them as a fingerprint rather than a specification — the asserts carry generous tolerances and will survive a numpy or scipy upgrade. If one fires, the message names the property that broke, and the function to re-read is the one it names.

7. mesh.py

7.1 Building a solid, cheaply

The obvious approach — duplicate the whole grid at z = 0 — doubles vertex count for a face nobody sees. Instead: full grid on top, a boundary ring of vertices at z = 0, quad walls stitching ring to ring, and a triangle fan across the floor from a centre point. That's 331,445 vertices instead of 658,240, and the result is genuinely watertight (verified, not assumed).

The fan needs a centre vertex rather than fanning from a corner. Fanning from a corner produces zero-area triangles for every ring point collinear with it, and degenerate facets make some slicers complain.

7.2 No supports, ever

Worth noticing: a heightfield is by construction a function z = f(x, y). Every point on the top surface is visible from directly above. There are no overhangs anywhere, at any α, no matter how steep the compression leaves the cliffs. This geometry can never need support material. Measured on a real conditioned field: across all fifty relief layers, the number of pixels in a layer not supported by the layer beneath it is zero — guaranteed, because {h ≥ z} only ever shrinks as z rises.

Bambu Studio will nonetheless warn you about floating regions on the first slice, and it is not wrong to. That warning is about toolpaths rather than overhangs — parts of the relief are finer than the nozzle can lay down, so the slicer drops them and then notices material with nothing printed beneath. It is expected, it is harmless, and §10 works through the numbers. The one thing not to do is take it up on the offer of supports: there is nothing for them to hold up, and they would be generated against the relief surface you just spent all this effort on.

7.3 The file

Five blocks, in order, make up mesh.py. The middle three are all one function — heightfield_to_mesh — split at its natural seams.

"""Heightfield -> watertight solid -> binary STL."""
import numpy as np
import trimesh

Grid coordinates in millimetres, one vertex per sample, two triangles per cell. ys counts down so image row 0 lands at the top of the panel.

def heightfield_to_mesh(h01, width_mm, height_mm, base_mm, relief_mm):
    """h01 is (H, W) in [0, 1]; row 0 is the TOP of the image."""
    H, W = h01.shape
    xs = np.linspace(0.0, width_mm, W)
    ys = np.linspace(height_mm, 0.0, H)
    X, Y = np.meshgrid(xs, ys)
    Z = base_mm + relief_mm * np.asarray(h01, np.float64)

    top = np.stack([X.ravel(), Y.ravel(), Z.ravel()], axis=1)
    idx = np.arange(H * W).reshape(H, W)

    # two triangles per cell, wound counter-clockwise seen from +Z
    a = idx[:-1, :-1]; b = idx[:-1, 1:]; c = idx[1:, 1:]; d = idx[1:, :-1]
    top_faces = np.concatenate([
        np.stack([d.ravel(), c.ravel(), b.ravel()], axis=1),
        np.stack([d.ravel(), b.ravel(), a.ravel()], axis=1),
    ])