-
-
Save tbbooher/4290b9661cb041c5f245f6fd260a47c2 to your computer and use it in GitHub Desktop.
maketree.py — E57 tree scan + surveyed height -> watertight tree mesh (.blend/.glb): AdTree skeleton, constrained growth, voxel-remesh union
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python | |
| """maketree.py -- E57 scan + surveyed height -> one watertight tree mesh. | |
| .venv/bin/python scripts/maketree.py <scan.e57> <height_ft> [--outdir DIR] | |
| [--stem ID] [--seed 7] [--voxel-fuse 0.03] [--blender EXE] | |
| [--measured-geom grown|adtree-mesh] | |
| The settled 2026-08-30 technique, end to end: | |
| 1 decode (stride 3) + stem inventory e57_read/tree_inventory | |
| 2 isolate best stem, full-res crop, ground strip isolate_tree | |
| 3 2 cm voxel decimate decimate_cloud | |
| 4 AdTree skeleton (macOS binary) reference/AdTree | |
| 5 semantic decimation + repair + connect + cloud-prune structure_* / stub_prune | |
| 6 auto-generate constraints from a species template scaled | |
| to THE TWO INPUTS (height; DBH comes from the inventory), | |
| seeds derived from the measured skeleton | |
| 7 constrained growth + gap gate tree_grow / grown_gap_check | |
| 8 base mesh fit to 10 ft: occupancy grid + marching cubes over ALL points | |
| 9 TreeQSM bridge -> Blender -> voxel-remesh UNION of tubes + base | |
| Output: <outdir>/<stem>-tree.blend (+ .glb), report JSON on stdout per stage. | |
| Assumes: trees/.venv (numpy/scipy/sklearn/skimage), the AdTree arm64 binary, | |
| Blender discoverable (--blender / $BLENDER_BIN / ~/Applications / /Applications | |
| / PATH). Windows note: AdTree is macOS-only; use run_trees.py's TreeAIBox path | |
| there instead. | |
| """ | |
| import argparse | |
| import csv | |
| import json | |
| import math | |
| import os | |
| import shutil | |
| import subprocess | |
| import sys | |
| import time | |
| from pathlib import Path | |
| import numpy as np | |
| SCRIPTS = Path(__file__).resolve().parent | |
| TREES = SCRIPTS.parent | |
| ADTREE = TREES.parent / "reference" / "AdTree" / "build" / "bin" / "AdTree.app" / "Contents" / "MacOS" / "AdTree" | |
| PY = sys.executable | |
| def sh(argv, name, ok_codes=(0,)): | |
| t0 = time.time() | |
| r = subprocess.run([str(a) for a in argv], capture_output=True, text=True) | |
| dt = time.time() - t0 | |
| if r.returncode not in ok_codes: | |
| sys.stderr.write(r.stdout[-3000:] + r.stderr[-3000:]) | |
| sys.exit(f"STAGE FAILED: {name} (exit {r.returncode}, {dt:.0f}s)") | |
| print(f" [{name:<12}] ok {dt:6.1f}s") | |
| return r.stdout | |
| def find_blender(cli): | |
| for c in [cli, os.environ.get("BLENDER_BIN"), | |
| str(Path.home() / "Applications/Blender.app/Contents/MacOS/Blender"), | |
| "/Applications/Blender.app/Contents/MacOS/Blender", | |
| shutil.which("blender")]: | |
| if c and Path(c).is_file(): | |
| return c | |
| sys.exit("no Blender found (--blender / $BLENDER_BIN / ~/Applications / PATH)") | |
| def pick_stem(inventory_csv, wanted): | |
| rows = list(csv.DictReader(open(inventory_csv, encoding="utf-8"))) | |
| if not rows: | |
| sys.exit("inventory found no stems") | |
| if wanted: | |
| row = next((r for r in rows if r["id"] == wanted), None) | |
| if not row: | |
| sys.exit(f"stem {wanted} not in inventory ({[r['id'] for r in rows]})") | |
| return row | |
| # best = quality letter, then largest DBH | |
| return sorted(rows, key=lambda r: (r["quality"], -float(r["dbh_cm"])))[0] | |
| def make_constraints(structure_path, height_m, dbh_cm, out_path): | |
| """Species-template constraints scaled to the two inputs. Template = | |
| the validated cedar-elm S01 set (2026-08-30); envelope scales with height, | |
| seeds derived from the measured skeleton (leaf branches reaching the upper | |
| 40% of the scan; leader if within ~35 deg of vertical).""" | |
| d = json.loads(Path(structure_path).read_text(encoding="utf-8")) | |
| zmax = max(n["p"][2] for b in d["branches"] for n in b["nodes"]) | |
| # Seed every branch whose TIP is unclaimed -- no child attaches within | |
| # 0.3 m of it -- and reaches the upper half of the scan. Having children | |
| # mid-branch must NOT disqualify: the trunk always has children, and a | |
| # trunk without a seed ends in a sliced-off cap (found on the first | |
| # non-cabana run, 2026-08-30 night). | |
| children_bases = {} | |
| for b in d["branches"]: | |
| if b.get("parent"): | |
| children_bases.setdefault(b["parent"], []).append( | |
| np.array(b["nodes"][0]["p"], dtype=float)) | |
| seeds = [] | |
| for b in d["branches"]: | |
| pts = np.array([n["p"] for n in b["nodes"]], dtype=float) | |
| tip = pts[-1] | |
| if tip[2] < 0.5 * zmax: | |
| continue | |
| claimed = any(np.linalg.norm(cb - tip) < 0.3 | |
| for cb in children_bases.get(b["id"], [])) | |
| if claimed: | |
| continue | |
| v = tip - pts[max(len(pts) - 3, 0)] | |
| v = v / (np.linalg.norm(v) + 1e-9) | |
| is_root = not b.get("parent") | |
| seeds.append({"branch": b["id"], | |
| "kind": "leader" if (is_root or v[2] > math.cos(math.radians(35))) | |
| else "scaffold", | |
| "start_r_override": None}) | |
| if not seeds: # degenerate scan: seed the root | |
| root = next(b["id"] for b in d["branches"] if not b.get("parent")) | |
| seeds = [{"branch": root, "kind": "leader", "start_r_override": None}] | |
| s = height_m / 14.0 # template tree was 14.0 m | |
| xy = np.array([n["p"][:2] for b in d["branches"] for n in b["nodes"]]) | |
| centre = [round(float(xy[:, 0].mean()), 2), round(float(xy[:, 1].mean()), 2)] | |
| env = [[3.0, 4.0], [5.0, 5.5], [8.0, 7.0], [11.0, 6.8], [13.0, 5.0], [14.0, 2.5]] | |
| cons = { | |
| "meta": {"tree": Path(structure_path).stem, | |
| "authored": "maketree.py species template (cedar-elm base), auto-scaled", | |
| "hard_measurements": { | |
| "total_height_m": round(height_m, 2), | |
| "height_source": "user input (ft, converted)", | |
| "crown_spread_m": round(0.93 * height_m, 1), | |
| "spread_source": "template ratio [estimate]"}}, | |
| "species": "template: open-grown broadleaf (cedar-elm derived)", | |
| "form": {"habit": "sympodial, co-dominant leaders, vase", | |
| "crown_centre_xy": centre, | |
| "crown_envelope_z_to_radius": [[round(z * s, 2), round(r * s, 2)] for z, r in env], | |
| "clear_trunk_height_m": 0.8, | |
| "branch_angle_secondary_deg": [35, 75], | |
| "apical_dominance": "low", "gravitropism": "moderate"}, | |
| "growth_rules": {"segment_length_m": [0.9, 1.6], | |
| "fork_every_m_leader": [1.6, 2.6], | |
| "fork_every_m_scaffold": [2.0, 3.2], | |
| "fork_children": {"two": 0.75, "three": 0.25}, | |
| "fork_main_share": [0.55, 0.68], | |
| "pipe_model_exponent": 2.49, | |
| "taper_between_forks_r_per_m": 0.008, | |
| "stop_radius_m": 0.03, "max_depth": 4}, | |
| "mechanics": {"green_wood_density_kg_m3": 900, | |
| "allowable_bending_stress_mpa": 8.0, | |
| "max_slenderness_L_over_d": 45}, | |
| "basal": {"ground_radius_m": round(max(0.005 * float(dbh_cm), 0.05) * 0.99, 3), | |
| "note": "= DBH/2 (audit: basal flare ~1.0); the fused base " | |
| "mesh covers the bottom anyway"}, | |
| "continuation_seeds": seeds, | |
| } | |
| Path(out_path).write_text(json.dumps(cons, indent=1), encoding="utf-8") | |
| return len(seeds) | |
| def base_fit(e57_npz, cx, cy, gz, out_ply, z_top=3.05, vox=0.035, | |
| stem_xyz=None): | |
| from scipy import ndimage | |
| from skimage import measure | |
| # Scan-envelope rework (2026-08-30): evidence set = the ISOLATED stem | |
| # cloud when available (isolate_tree already separated the clutter that | |
| # a raw re-crop of cloud.npz re-admits), already in the structure frame. | |
| if stem_xyz is not None and Path(stem_xyz).is_file(): | |
| p = np.fromfile(stem_xyz, sep=" ").reshape(-1, 3).astype(np.float64) | |
| else: | |
| d = np.load(e57_npz) | |
| xyz = d["xyz"] | |
| m = (np.abs(xyz[:, 0] - cx) < 2.0) & (np.abs(xyz[:, 1] - cy) < 2.0) | |
| p = xyz[m].astype(np.float64) | |
| p[:, 0] -= cx; p[:, 1] -= cy; p[:, 2] -= gz | |
| p = p[(p[:, 2] > 0.04) & (p[:, 2] <= z_top + 0.05)] | |
| if len(p) < 20000: | |
| print(f" [base-fit ] skipped: only {len(p):,} points below {z_top} m") | |
| return False | |
| r = np.hypot(p[:, 0] - p[:, 0].mean(), p[:, 1] - p[:, 1].mean()) | |
| p = p[r < 1.6] | |
| lo = p.min(axis=0) - 3 * vox | |
| idx = np.floor((p - lo) / vox).astype(int) | |
| raw = np.zeros(idx.max(axis=0) + 4, dtype=bool) | |
| raw[tuple(idx.T)] = True | |
| # Closing gap limit = 7 cm (2 voxels): the old closing | |
| # (structure(3,2) x 2) bridged 14-24 cm diagonally -- the mechanism that | |
| # fused trunk to off-stem returns and doubled the 0.75-1.50 m band | |
| # footprint. Close conservatively, then the HARD support mask: no voxel | |
| # farther than 2 voxels (city-block) from a scanned voxel survives. | |
| grid = ndimage.binary_closing(raw, ndimage.generate_binary_structure(3, 1), | |
| iterations=1) | |
| grid = ndimage.binary_opening(grid, iterations=1) | |
| support = ndimage.binary_dilation(raw, iterations=2) | |
| unsupported_pre = grid & ~support | |
| grid &= support | |
| lab, _ = ndimage.label(grid) | |
| sizes = np.bincount(lab.ravel()); sizes[0] = 0 | |
| grid = lab == sizes.argmax() | |
| # per-band stats: voxels the closing invented beyond scan support (culled) | |
| zs = lo[2] + (np.arange(grid.shape[2]) + 0.5) * vox | |
| for b0 in np.arange(0.0, z_top + 0.25, 0.25): | |
| zm = (zs >= b0) & (zs < b0 + 0.25) | |
| n_g = int(grid[:, :, zm].sum()) | |
| n_u = int(unsupported_pre[:, :, zm].sum()) | |
| if n_g or n_u: | |
| print(f" [base-fit ] band {b0:4.2f}-{b0 + 0.25:4.2f} m: " | |
| f"{n_g:6d} vox kept, {n_u:5d} unsupported culled") | |
| # smoothing field AFTER the support mask, so the 0.42 iso-surface cannot | |
| # migrate into unsupported space | |
| field = ndimage.gaussian_filter(grid.astype(np.float32), 1.1) | |
| verts, faces, _, _ = measure.marching_cubes(field, level=0.42, spacing=(vox,) * 3) | |
| verts += lo | |
| with open(out_ply, "wb") as fh: | |
| fh.write(f"ply\nformat binary_little_endian 1.0\nelement vertex {len(verts)}\n" | |
| f"property float x\nproperty float y\nproperty float z\n" | |
| f"element face {len(faces)}\nproperty list uchar int vertex_indices\n" | |
| f"end_header\n".encode()) | |
| verts.astype("<f4").tofile(fh) | |
| fc = np.full((len(faces), 1), 3, dtype="u1") | |
| fh.write(np.hstack([fc, faces.astype("<i4").view("u1").reshape(len(faces), 12)]).tobytes()) | |
| print(f" [base-fit ] ok {len(verts):,} verts / {len(faces):,} tris -> {out_ply}") | |
| return True | |
| def main(): | |
| ap = argparse.ArgumentParser() | |
| ap.add_argument("e57") | |
| ap.add_argument("height_ft", type=float) | |
| ap.add_argument("--outdir", default=None) | |
| ap.add_argument("--stem", default=None, help="stem id (default: best quality/largest DBH)") | |
| ap.add_argument("--seed", type=int, default=7) | |
| ap.add_argument("--voxel-fuse", type=float, default=0.03) | |
| ap.add_argument("--blender", default=None) | |
| ap.add_argument("--measured-geom", choices=["grown", "adtree-mesh"], | |
| default="grown", | |
| help="'grown' (default): measured zone from the grown " | |
| "tubes + base fit only. 'adtree-mesh': ALSO union in " | |
| "AdTree's <stem>_branches.obj after a scan-support " | |
| "cull (adtree_mesh_clean.py drops every unwelded " | |
| "component whose median/95th-pct vertex->cloud " | |
| "distance exceeds 0.06/0.15 m -- including the " | |
| "clutter-fat trunk tube and the crumpled crotch " | |
| "cylinders -- and open3d-decimates the rest); the grower " | |
| "still consumes the pruned structure graph for " | |
| "topology/seeds") | |
| a = ap.parse_args() | |
| e57 = Path(a.e57).resolve() | |
| if not e57.is_file(): | |
| sys.exit(f"no such scan: {e57}") | |
| if not (2 < a.height_ft < 300): | |
| sys.exit(f"height {a.height_ft} ft is not a tree height") | |
| height_m = a.height_ft * 0.3048 | |
| if not ADTREE.is_file(): | |
| sys.exit(f"AdTree binary missing at {ADTREE} (macOS-only stage; " | |
| "on Windows use run_trees.py's TreeAIBox path)") | |
| blender = find_blender(a.blender) | |
| out = Path(a.outdir) if a.outdir else TREES / "out" / f"maketree-{e57.stem}".replace(" ", "-") | |
| out.mkdir(parents=True, exist_ok=True) | |
| print(f"maketree: {e57.name} height {a.height_ft} ft = {height_m:.2f} m -> {out}") | |
| npz = out / "cloud.npz" | |
| sh([PY, SCRIPTS / "e57_read.py", e57, "--stride", "3", "--out", npz], "decode") | |
| sh([PY, SCRIPTS / "tree_inventory.py", npz, "--outdir", out], "inventory") | |
| # Stem fallback (found on the two-stem cabana scan 2026-08-30): the | |
| # top-ranked stem's crop can catch a NEIGHBOUR's branches, leaving an | |
| # AdTree skeleton that honestly cannot be repaired to a single root. | |
| # Rather than force parents into existence, fall through to the next stem. | |
| rows = list(csv.DictReader(open(out / "tree-inventory.csv", encoding="utf-8"))) | |
| first = pick_stem(out / "tree-inventory.csv", a.stem) | |
| candidates = [first] + [r for r in sorted(rows, key=lambda r: (r["quality"], -float(r["dbh_cm"]))) | |
| if r["id"] != first["id"]] | |
| stem_base = row = None | |
| for cand in candidates: | |
| stem, low = cand["id"], cand["id"].lower() | |
| print(f" stem {stem}: DBH {cand['dbh_cm']} cm, quality {cand['quality']}") | |
| sh([PY, SCRIPTS / "isolate_tree.py", "--e57", e57, "--inventory", | |
| out / "tree-inventory.csv", "--stem", stem, "--outdir", out / "stems"], "isolate") | |
| small = out / "stems" / f"{low}-2cm.xyz" | |
| sh([PY, SCRIPTS / "decimate_cloud.py", out / "stems" / f"{low}.xyz", | |
| "--voxel", "0.02", "--out", small], "decimate") | |
| ad = out / "adtree" | |
| ad.mkdir(exist_ok=True) | |
| sh([ADTREE, small, ad, "-s"], "adtree", ok_codes=(0, 1)) | |
| skel = ad / f"{low}-2cm_skeleton.ply" | |
| if not skel.is_file(): | |
| skel = next(ad.glob(f"{low}*_skeleton.ply")) | |
| base = str(skel).removesuffix(".ply") | |
| sh([PY, SCRIPTS / "structure_graph.py", skel, "--cloud", small, "--out", base], | |
| "structure") | |
| rep = subprocess.run([PY, str(SCRIPTS / "structure_repair.py"), | |
| base + ".structure.json", | |
| "--out", base + ".structure.repaired.json", | |
| "--report", str(ad / f"{low}-repair-report.txt")], | |
| capture_output=True, text=True) | |
| if rep.returncode == 0: | |
| print(f" [repair ] ok") | |
| stem_base, row = base, cand | |
| break | |
| print(f" [repair ] REFUSED for {stem} (not single-rootable); " | |
| f"trying next stem") | |
| if stem_base is None: | |
| sys.exit("no inventory stem yields a single-rootable skeleton -- " | |
| "inspect the scan (multi-tree overlap?) or pass --stem") | |
| sh([PY, SCRIPTS / "structure_connect.py", stem_base + ".structure.repaired.json", | |
| "--out", stem_base + ".structure.connected.json"], "connect") | |
| sh([PY, SCRIPTS / "stub_prune.py", stem_base + ".structure.connected.json", | |
| "--cloud", small, "--out", stem_base + ".structure.pruned.json"], "prune") | |
| cons = out / f"{low}-constraints.json" | |
| nseeds = make_constraints(stem_base + ".structure.pruned.json", | |
| height_m, row["dbh_cm"], cons) | |
| print(f" [constraints ] ok template scaled to {height_m:.1f} m, {nseeds} auto-seeds") | |
| grown = out / f"{low}-grown" | |
| # full-res stem cloud (not the 2 cm decimation): the scan-envelope clamp | |
| # inside tree_grow needs every return it can get per ring slab | |
| full_cloud = out / "stems" / f"{low}.xyz" | |
| sh([PY, SCRIPTS / "tree_grow.py", stem_base + ".structure.pruned.json", cons, | |
| "--seed", str(a.seed), "--cloud", full_cloud, "--out", grown], "grow") | |
| sh([PY, SCRIPTS / "grown_gap_check.py", str(grown) + ".json"], "gap-gate") | |
| base_ply = out / f"{low}-base-fit.ply" | |
| have_base = base_fit(npz, float(row["x_m"]), float(row["y_m"]), | |
| float(row["ground_z_m"]), base_ply, | |
| stem_xyz=full_cloud) | |
| # Optional AdTree-mesh measured zone (2026-08-30): union in the fitted | |
| # branches.obj cylinders wherever the scan confirms them. The cleaner is | |
| # itself an envelope filter -- every kept component has median | |
| # vertex->cloud distance <= 0.06 m and 95th-pct <= 0.15 m, so the | |
| # clutter-fat trunk tube (median 0.31 m on the big scan), the crumpled | |
| # crotch cylinders, and the below-floor needle spikes are dropped, and | |
| # radii that AdTree over-fit never reach the union. | |
| # Decimation runs only if the interpreter has open3d, so prefer .venv-o3d. | |
| measured_ply = None | |
| if a.measured_geom == "adtree-mesh": | |
| branches_obj = Path(str(skel).removesuffix("_skeleton.ply") + "_branches.obj") | |
| if not branches_obj.is_file(): | |
| branches_obj = next(ad.glob(f"{low}*_branches.obj")) | |
| o3d_py = TREES / ".venv-o3d" / "bin" / "python" | |
| measured_ply = out / f"{low}-measured-mesh.ply" | |
| sh([o3d_py if o3d_py.is_file() else PY, SCRIPTS / "adtree_mesh_clean.py", | |
| branches_obj, "--cloud", small, "--out", measured_ply, | |
| "--decimate-ratio", "0.5", | |
| "--report", out / f"{low}-measured-mesh-report.json"], "mesh-clean") | |
| cyl = out / f"{low}.cylinders.txt" | |
| sh([PY, SCRIPTS / "structure_to_qsm.py", str(grown) + ".json", "--out", cyl], | |
| "export-qsm") | |
| blend = out / f"{low}-tree.blend" | |
| argv = [blender, "-b", "--factory-startup", "--python", | |
| SCRIPTS / "blender_qsm_import.py", "--", "--cylinders", cyl, | |
| "--output", blend, "--fuse", str(a.voxel_fuse), "--glb"] | |
| if have_base: | |
| argv += ["--merge-ply", base_ply] | |
| if measured_ply is not None: | |
| argv += ["--merge-ply", measured_ply] | |
| stdout = sh(argv, "blender-fuse") | |
| # THE LOUD GATE (2026-08-30): wherever the scan has coverage, the fuse | |
| # INPUTS (cylinders + base PLY, which bound the fused mesh to within one | |
| # fuse voxel) must not exceed the scan's occupied hull. Nonzero exit | |
| # kills the run -- a model fatter than the scan never ships. | |
| gate_argv = [PY, SCRIPTS / "scan_envelope_gate.py", "--grown", | |
| str(grown) + ".json", "--cylinders", cyl, | |
| "--cloud", full_cloud, "--fuse-vox", str(a.voxel_fuse), | |
| "--report", out / f"{low}-envelope-gate.json"] | |
| if have_base: | |
| gate_argv += ["--base-ply", base_ply] | |
| sh(gate_argv, "envelope-gate") | |
| line = next((ln for ln in stdout.splitlines() if ln.startswith("QSM_MESH ")), "") | |
| print(f"\nDONE {blend}") | |
| if line: | |
| print(line) | |
| if __name__ == "__main__": | |
| main() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment