CMPlanning/scripts/build_skeleton.py

#!/usr/bin/env python3
"""Task 7: axis + standard-section sweep skeleton of the main line (본선), 4+745..5+480.

Inputs (all confirmed by the user, 2026-10-21):
  - centreline nodes  data/alignment_plan_001_centre.csv   (midpoint of side ticks)
  - profile           data/profile_gl_002_draft.csv        (A = 지반고 GL, B = 계획고 = 노면)
  - section table     DOC/03_설계/설계_표준단면치수표_2026-10-21_v0.1.00.md
      R1 6,982 mm (lining inner), lining 300 mm -> outer R 7,282, excavation R 7,402
      crown A = 6.615 m above road (F = 0.000 = road reference)
  - cross slope: 4+745~5+105 -2 %, 5+105~5+165 linear -2 -> +2 %, 5+165~5+480 +2 %
Model (skeleton, not as-built):
  - Section = arch (concentric centre c = 6.615 - 6.982) + vertical walls down to the road.
    Walls and invert are APPROXIMATIONS and are flagged in the metadata. No invert.
  - Superelevation: z += (S/100) * x across the section (sign convention ASSUMED, flagged).
  - Coordinates: plan XY in m (station spacing 100 m matches plan units), origin at 4+745.
Outputs: model/tunnel_skeleton.glb, model/tunnel_skeleton_meta.json, data/verify_skeleton.json
"""
import csv
import json
import math
import sys

import numpy as np
import trimesh

ROOT = sys.argv[1]  # Code/3d-drawing
C_CSV = f"{ROOT}/data/alignment_plan_001_centre.csv"
P_CSV = f"{ROOT}/data/profile_gl_002_draft.csv"
R1, R_OUT, R_EX = 6.982, 7.282, 7.402           # m (mm / 1000)
CROWN_A = 6.615                                  # m above road
CENTRE_C = CROWN_A - R1                          # arch centre height above road (m)
S0, S1 = 4745, 5480                              # modelled range (fully determined)
DS = 5.0                                         # station step for rings (m)
ANG = np.radians(np.linspace(40.0, 140.0, 41))   # arch arc, right springline -> left

def read_centre():
    rows = [r for r in csv.DictReader(open(C_CSV, encoding="utf-8"))]
    return np.array([[int(r["station_m"]), float(r["centre_x"]), float(r["centre_y"])] for r in rows])

def read_profile():
    gl, dg = {}, {}
    for r in csv.DictReader(open(P_CSV, encoding="utf-8")):
        sta = int(r["station_anchor"])
        (gl if float(r["across_m"]) < -25 else dg)[sta] = float(r["value"])
    return gl, dg

def slope_pct(s):
    if s < 5105: return -2.0
    if s < 5165: return -2.0 + 4.0 * (s - 5105) / 60.0
    return 2.0

cn = read_centre()
gl_d, dg_d = read_profile()
gl_s = np.array(sorted(gl_d)); gl_v = np.array([gl_d[k] for k in gl_s])
dg_s = np.array(sorted(dg_d)); dg_v = np.array([dg_d[k] for k in dg_s])
cn_s = cn[:, 0]

def axis_xy(s):
    x = np.interp(s, cn_s, cn[:, 1]); y = np.interp(s, cn_s, cn[:, 2])
    i = min(max(np.searchsorted(cn_s, s) - 1, 0), len(cn_s) - 2)
    dx, dy = cn[i + 1, 1] - cn[i, 1], cn[i + 1, 2] - cn[i, 2]
    L = math.hypot(dx, dy)
    return x, y, dx / L, dy / L

ORIGIN = np.array(axis_xy(S0)[:2])
stations = np.arange(S0, S1 + 0.001, DS)

def section(R, s, slope):
    """ring (x lateral m, z above road m) for radius R, including superelevation."""
    xs = np.concatenate([[R * math.cos(ANG[0])], R * np.cos(ANG), [R * math.cos(ANG[-1])]])
    zs = np.concatenate([[0.0], CENTRE_C + R * np.sin(ANG), [0.0]])
    zs = zs + slope * xs
    return xs, zs

def ring_world(R, s):
    x, y, tx, ty = axis_xy(s)
    nx, ny = -ty, tx                           # left normal (horizontal)
    zr = np.interp(s, dg_s, dg_v)              # road (계획고 = 노면)
    xs, zs = section(R, s, slope_pct(s) / 100.0)
    P = np.column_stack([x - ORIGIN[0] + xs * nx, y - ORIGIN[1] + xs * ny, zr + zs])
    return P

def make_band(R, color, name):
    rings = [ring_world(R, s) for s in stations]
    n = len(rings[0]); m = len(rings)
    V = np.vstack(rings)
    F = []
    for i in range(m - 1):
        for j in range(n - 1):
            a = i * n + j; b = a + 1; c = a + n; d = c + 1
            F += [[a, b, c], [b, d, c]]
    mesh = trimesh.Trimesh(vertices=V, faces=np.array(F), process=False)
    mesh.visual.face_colors = np.tile(np.array(color, dtype=np.uint8), (len(F), 1))
    mesh.metadata["name"] = name
    return mesh

def make_ribbon(z_fn, half_w, color, name):
    V, F = [], []
    for k, s in enumerate(stations):
        x, y, tx, ty = axis_xy(s)
        nx, ny = -ty, tx
        zr = z_fn(s)
        V.append([x - ORIGIN[0] - half_w * nx, y - ORIGIN[1] - half_w * ny, zr])
        V.append([x - ORIGIN[0] + half_w * nx, y - ORIGIN[1] + half_w * ny, zr])
    V = np.array(V)
    for k in range(len(stations) - 1):
        a = 2 * k; b = a + 1; c = a + 2; d = a + 3
        F += [[a, c, b], [b, c, d]]
    mesh = trimesh.Trimesh(vertices=V, faces=np.array(F), process=False)
    mesh.visual.face_colors = np.tile(np.array(color, dtype=np.uint8), (len(F), 1))
    mesh.metadata["name"] = name
    return mesh

parts = [
    make_band(R_OUT, (190, 190, 186, 255), "lining_outer_R7282"),
    make_band(R1,    (235, 235, 230, 255), "lining_inner_R6982"),
    make_band(R_EX,  (150, 120, 90, 110),  "excavation_R7402"),
    make_ribbon(lambda s: float(np.interp(s, dg_s, dg_v)), 5.0, (70, 70, 75, 255), "road_surface_계획고"),
    make_ribbon(lambda s: float(np.interp(s, gl_s, gl_v)), 20.0, (110, 150, 90, 200), "ground_GL_지반고"),
]
scene = trimesh.Scene()
for p in parts:
    scene.add_geometry(p, geom_name=p.metadata["name"])
glb = f"{ROOT}/model/tunnel_skeleton.glb"
scene.export(glb)

# ---- verification (assertions recorded, not forced) ----
checks = {}
# radius check in section plane: distance from arch centre (lateral=0, z=centre_c)
xs, zs = section(R1, S0, 0.0)
rad = np.hypot(xs[1:-1], zs[1:-1] - CENTRE_C)
checks["inner_radius_max_dev_m"] = float(np.max(np.abs(rad - R1)))
checks["lining_thickness_crown_m"] = round(R_OUT - R1, 3)
crown_z = ring_world(R1, 5000)[:, 2].max() - np.interp(5000, dg_s, dg_v)
checks["crown_above_road_m (expect 6.615 at S=0 station 5000)"] = round(float(crown_z), 4)
checks["road_z_4745_m"] = round(float(np.interp(S0, dg_s, dg_v)), 3)
checks["road_z_5480_m"] = round(float(np.interp(S1, dg_s, dg_v)), 3)
L_chord = float(np.sum([math.dist(cn[i, 1:], cn[i + 1, 1:]) for i in range(len(cn) - 1)]))
checks["plan_chord_length_total_m (4+700..5+500)"] = round(L_chord, 3)
checks["station_range_modelled"] = f"{S0} .. {S1} (ring step {DS} m, {len(stations)} rings)"
checks["glb_bytes"] = int(__import__("os").path.getsize(glb))
reload = trimesh.load(glb)
checks["glb_reload_geometries"] = len(reload.geometry)
json.dump(checks, open(f"{ROOT}/data/verify_skeleton.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2)
meta = {
  "model": "tunnel_skeleton.glb", "scope": f"{S0}~{S1} (4+745~5+480)",
  "units": "m", "origin_plan_xy": [round(ORIGIN[0], 3), round(ORIGIN[1], 3)],
  "sources": {"centreline": "C0060302-001 CA-ALGN-STSM (side-tick midpoint)",
              "road_grade_B": "C0060302-002 CR-TABL-TEX2 (계획고=노면)",
              "GL_A": "C0060302-002 CR-TABL-TEX2 (지반고)",
              "section": "C0060303-001 CS-TABL-TEX1/2, CS-CONC-MAJR, CS-EXCV"},
  "approximations": ["walls vertical down to road; invert not modelled",
                     "lining band thickness at walls is radial-consistent but not as-built",
                     "centreline and grade linearly interpolated between nodes",
                     "superelevation sign convention assumed (z += S/100 * x, x = left of tangent)"],
  "layers": [p.metadata["name"] for p in parts],
  "open": ["4+700~4+745 and 5+480~5+500 excluded ([미확정])", "pattern zones (03-003~008) pending Task 5",
           "lining details (04-004~009) pending Task 6"]
}
json.dump(meta, open(f"{ROOT}/model/tunnel_skeleton_meta.json", "w", encoding="utf-8"), ensure_ascii=False, indent=2)
print(json.dumps(checks, ensure_ascii=False, indent=1))