RANSAC deck plane visualization¶

Visualize the RANSAC-fitted deck plane for random source LAZ samples. Uses the same fitting logic as the weak-supervision pipeline. Configure number of random HUCs (M) and random bridges per HUC (N) below. Outputs 3D figures (points colored by inlier/outlier + plane surface) to notebooks/outputs/.

Config¶

Paths, sampling counts, display limit, and output settings.

In [1]:

import os
from pathlib import Path

SOURCE_LAZ_DIR = Path(os.environ.get("BRIDGE_SOURCE_LAZ_DIR", "data/ml-data/silver_training"))

N_RANDOM_HUCS = 3
N_RANDOM_BRIDGES_PER_HUC = 3
SEED = 27

MAX_POINTS_DISPLAY = 50_000
POINT_SIZE = 6
EDGE_LINEWIDTH = 0.0

# 3D viewpoint
VIEW_ELEV = 30
VIEW_AZIM = 225

try:
    _nb_dir = os.path.dirname(os.path.abspath(__file__))
except NameError:
    _nb_dir = os.getcwd()
if os.path.basename(_nb_dir) != "notebooks":
    _nb_dir = os.path.join(_nb_dir, "notebooks")
OUTPUT_DIR = Path(os.environ.get("BRIDGE_VIZ_OUTPUT_DIR", os.path.join(_nb_dir, "outputs")))
OUTPUT_PREFIX = "ransac_deck_plane"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

import os from pathlib import Path SOURCE_LAZ_DIR = Path(os.environ.get("BRIDGE_SOURCE_LAZ_DIR", "data/ml-data/silver_training")) N_RANDOM_HUCS = 3 N_RANDOM_BRIDGES_PER_HUC = 3 SEED = 27 MAX_POINTS_DISPLAY = 50_000 POINT_SIZE = 6 EDGE_LINEWIDTH = 0.0 # 3D viewpoint VIEW_ELEV = 30 VIEW_AZIM = 225 try: _nb_dir = os.path.dirname(os.path.abspath(__file__)) except NameError: _nb_dir = os.getcwd() if os.path.basename(_nb_dir) != "notebooks": _nb_dir = os.path.join(_nb_dir, "notebooks") OUTPUT_DIR = Path(os.environ.get("BRIDGE_VIZ_OUTPUT_DIR", os.path.join(_nb_dir, "outputs"))) OUTPUT_PREFIX = "ransac_deck_plane" OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

Imports¶

Repo root and pipeline module (filename has hyphen; use importlib).

In [ ]:

import random
import sys
import json
import numpy as np
import matplotlib.pyplot as plt
import pdal

_repo_root = Path.cwd() if (Path.cwd() / "src").is_dir() else Path.cwd().parent
if str(_repo_root) not in sys.path:
    sys.path.insert(0, str(_repo_root))

from src.weak_supervision import BridgeProcessingConfig, fit_ransac_from_arrays

import random import sys import json import numpy as np import matplotlib.pyplot as plt import pdal _repo_root = Path.cwd() if (Path.cwd() / "src").is_dir() else Path.cwd().parent if str(_repo_root) not in sys.path: sys.path.insert(0, str(_repo_root)) from src.weak_supervision import BridgeProcessingConfig, fit_ransac_from_arrays

Discovery¶

List (huc_id, stem) from source LAZ dir; sample M HUCs, then N bridges per HUC.

In [3]:

def discover_source_stems(source_root: Path):
    pairs = []
    for child in sorted(source_root.iterdir()):
        if not child.is_dir():
            continue
        huc_id = child.name
        for p in child.glob("*.laz"):
            pairs.append((huc_id, p.stem))
    return pairs


def sample_pairs(pairs: list, n_hucs: int, n_per_huc: int, seed: int):
    rng = random.Random(seed)
    by_huc = {}
    for huc_id, stem in pairs:
        by_huc.setdefault(huc_id, []).append(stem)
    huc_ids = list(by_huc.keys())
    if not huc_ids:
        return []
    n_hucs = min(n_hucs, len(huc_ids))
    chosen_hucs = rng.sample(huc_ids, n_hucs)
    result = []
    for huc_id in chosen_hucs:
        stems = by_huc[huc_id]
        k = min(n_per_huc, len(stems))
        for stem in rng.sample(stems, k):
            result.append((huc_id, stem))
    return result


all_pairs = discover_source_stems(Path(_repo_root / SOURCE_LAZ_DIR))
pairs_to_plot = sample_pairs(all_pairs, N_RANDOM_HUCS, N_RANDOM_BRIDGES_PER_HUC, SEED)

if not pairs_to_plot:
    print("No source LAZ files found. Check SOURCE_LAZ_DIR.")
else:
    hucs_used = sorted(set(p[0] for p in pairs_to_plot))
    print(f"Chosen HUCs: {hucs_used}")
    print(f"Pairs to plot: {len(pairs_to_plot)} — {pairs_to_plot}")

def discover_source_stems(source_root: Path): pairs = [] for child in sorted(source_root.iterdir()): if not child.is_dir(): continue huc_id = child.name for p in child.glob("*.laz"): pairs.append((huc_id, p.stem)) return pairs def sample_pairs(pairs: list, n_hucs: int, n_per_huc: int, seed: int): rng = random.Random(seed) by_huc = {} for huc_id, stem in pairs: by_huc.setdefault(huc_id, []).append(stem) huc_ids = list(by_huc.keys()) if not huc_ids: return [] n_hucs = min(n_hucs, len(huc_ids)) chosen_hucs = rng.sample(huc_ids, n_hucs) result = [] for huc_id in chosen_hucs: stems = by_huc[huc_id] k = min(n_per_huc, len(stems)) for stem in rng.sample(stems, k): result.append((huc_id, stem)) return result all_pairs = discover_source_stems(Path(_repo_root / SOURCE_LAZ_DIR)) pairs_to_plot = sample_pairs(all_pairs, N_RANDOM_HUCS, N_RANDOM_BRIDGES_PER_HUC, SEED) if not pairs_to_plot: print("No source LAZ files found. Check SOURCE_LAZ_DIR.") else: hucs_used = sorted(set(p[0] for p in pairs_to_plot)) print(f"Chosen HUCs: {hucs_used}") print(f"Pairs to plot: {len(pairs_to_plot)} — {pairs_to_plot}")

Chosen HUCs: ['01080102', '01080206', '01090002']
Pairs to plot: 9 — [('01080206', 'bridge_96744504_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018'), ('01080206', 'bridge_839739315_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018'), ('01080206', 'bridge_250569218_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018'), ('01080102', 'bridge_607138210_USGS_LPC_VT_Connecticut_River_2016_LAS_2018'), ('01080102', 'bridge_265909373_USGS_LPC_CT_River_Lot6_Winnipesaulee_2015_LAS_2017'), ('01080102', 'bridge_458466374_USGS_LPC_VT_Connecticut_River_2016_LAS_2018'), ('01090002', 'bridge_102809372_ARRA-LFTNE_Massachusetts_2011'), ('01090002', 'bridge_305358293_ARRA-LFTNE_Massachusetts_2011'), ('01090002', 'bridge_880364612_MA_CentralEastern_2_2021')]

Visualize¶

Load each source LAZ, run RANSAC via fit_ransac_from_arrays, then plot points (inlier/outlier) and fitted plane surface.

In [4]:

def load_laz_arrays(path: Path):
    pipeline_json = {"pipeline": [{"type": "readers.las", "filename": str(path)}]}
    pipeline = pdal.Pipeline(json.dumps(pipeline_json))
    pipeline.execute()
    return pipeline.arrays[0]


def set_equal_aspect_3d(ax):
    xlim, ylim, zlim = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    xr, yr, zr = xlim[1] - xlim[0], ylim[1] - ylim[0], zlim[1] - zlim[0]
    max_range = max(xr, yr, zr)
    xc = (xlim[0] + xlim[1]) / 2
    yc = (ylim[0] + ylim[1]) / 2
    zc = (zlim[0] + zlim[1]) / 2
    ax.set_xlim([xc - max_range / 2, xc + max_range / 2])
    ax.set_ylim([yc - max_range / 2, yc + max_range / 2])
    ax.set_zlim([zc - max_range / 2, zc + max_range / 2])


def style_3d_ax(ax):
    """Apply clean styling to a 3D axes."""
    ax.xaxis.pane.fill = True
    ax.yaxis.pane.fill = True
    ax.zaxis.pane.fill = True
    ax.xaxis.pane.set_facecolor((0.95, 0.95, 0.95, 1.0))
    ax.yaxis.pane.set_facecolor((0.92, 0.92, 0.92, 1.0))
    ax.zaxis.pane.set_facecolor((0.90, 0.90, 0.90, 1.0))
    ax.xaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5})
    ax.yaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5})
    ax.zaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5})
    ax.tick_params(axis="both", labelsize=7, pad=1)
    ax.set_xlabel("X (m)", fontsize=9, labelpad=4)
    ax.set_ylabel("Y (m)", fontsize=9, labelpad=4)
    ax.set_zlabel("Z (m)", fontsize=9, labelpad=4)
    ax.view_init(elev=VIEW_ELEV, azim=VIEW_AZIM)


def short_bridge_name(stem: str, max_len: int = 40) -> str:
    parts = stem.split("_")
    if len(parts) >= 2 and parts[0] == "bridge":
        bridge_id = parts[1]
        rest = "_".join(parts[2:])
        if len(rest) > max_len - len(bridge_id) - 8:
            rest = rest[: max_len - len(bridge_id) - 11] + "..."
        return f"bridge_{bridge_id}_{rest}"
    return stem[:max_len] + ("..." if len(stem) > max_len else "")


# Colors
INLIER_COLOR = "#27ae60"   # green
OUTLIER_COLOR = "#bdc3c7"  # light gray
PLANE_COLOR = "#3498db"    # blue

config = BridgeProcessingConfig()

if not pairs_to_plot:
    print("Skipping visualization (no pairs).")
else:
    for huc_id, stem in pairs_to_plot:
        path = _repo_root / SOURCE_LAZ_DIR / huc_id / f"{stem}.laz"
        try:
            arrays = load_laz_arrays(path)
        except Exception as e:
            print(f"Failed to load {huc_id}/{stem}: {e}")
            continue
        result = fit_ransac_from_arrays(arrays, config)
        if not result.get("success"):
            print(f"RANSAC failed {huc_id}/{stem}: {result.get('error', 'unknown')}")
            continue

        X_local = result["X_local"]
        Y_local = result["Y_local"]
        Z = result["Z"]
        inlier = result["inlier_mask_full"]
        coef_x = result["coef_x"]
        coef_y = result["coef_y"]
        intercept = result["intercept"]

        n_pts = len(X_local)
        if MAX_POINTS_DISPLAY and n_pts > MAX_POINTS_DISPLAY:
            rng = np.random.RandomState(SEED)
            idx = rng.choice(n_pts, size=MAX_POINTS_DISPLAY, replace=False)
            X_local, Y_local, Z = X_local[idx], Y_local[idx], Z[idx]
            inlier = inlier[idx]

        n_inlier = int(inlier.sum())
        n_outlier = int((~inlier).sum())

        fig = plt.figure(figsize=(13, 9), facecolor="white")
        ax = fig.add_subplot(111, projection="3d")

        # Plot outliers first (behind) then inliers on top
        ax.scatter(
            X_local[~inlier], Y_local[~inlier], Z[~inlier],
            c=OUTLIER_COLOR, s=POINT_SIZE, alpha=0.35,
            linewidths=EDGE_LINEWIDTH, rasterized=True,
            label=f"Outlier ({n_outlier:,})",
        )
        ax.scatter(
            X_local[inlier], Y_local[inlier], Z[inlier],
            c=INLIER_COLOR, s=POINT_SIZE, alpha=0.7,
            linewidths=EDGE_LINEWIDTH, rasterized=True,
            label=f"Inlier ({n_inlier:,})",
        )

        # Fitted plane surface
        xmin, xmax = X_local.min(), X_local.max()
        ymin, ymax = Y_local.min(), Y_local.max()
        margin = 0.05 * max(xmax - xmin, ymax - ymin, 1.0)
        xg = np.linspace(xmin - margin, xmax + margin, 50)
        yg = np.linspace(ymin - margin, ymax + margin, 50)
        Xg, Yg = np.meshgrid(xg, yg)
        Zg = coef_x * Xg + coef_y * Yg + intercept
        ax.plot_surface(
            Xg, Yg, Zg, alpha=0.25, color=PLANE_COLOR,
            edgecolor=PLANE_COLOR, linewidth=0.1,
        )

        display_name = short_bridge_name(stem)
        ax.set_title(
            f"RANSAC Deck Plane\nHUC {huc_id}  |  {display_name}",
            fontsize=12, fontweight="bold", pad=14, fontfamily="monospace",
        )

        leg = ax.legend(
            loc="upper left", fontsize=9, frameon=True,
            fancybox=True, framealpha=0.9, edgecolor="#cccccc",
            markerscale=5, handletextpad=0.5,
            borderpad=0.6, labelspacing=0.4,
        )
        leg.get_frame().set_facecolor("white")

        style_3d_ax(ax)
        set_equal_aspect_3d(ax)

        plt.tight_layout()
        out_path = OUTPUT_DIR / f"{OUTPUT_PREFIX}_{huc_id}_{stem}.png"
        fig.savefig(out_path, dpi=180, bbox_inches="tight", facecolor="white")
        plt.show()
        print(f"Saved: {out_path}")

def load_laz_arrays(path: Path): pipeline_json = {"pipeline": [{"type": "readers.las", "filename": str(path)}]} pipeline = pdal.Pipeline(json.dumps(pipeline_json)) pipeline.execute() return pipeline.arrays[0] def set_equal_aspect_3d(ax): xlim, ylim, zlim = ax.get_xlim(), ax.get_ylim(), ax.get_zlim() xr, yr, zr = xlim[1] - xlim[0], ylim[1] - ylim[0], zlim[1] - zlim[0] max_range = max(xr, yr, zr) xc = (xlim[0] + xlim[1]) / 2 yc = (ylim[0] + ylim[1]) / 2 zc = (zlim[0] + zlim[1]) / 2 ax.set_xlim([xc - max_range / 2, xc + max_range / 2]) ax.set_ylim([yc - max_range / 2, yc + max_range / 2]) ax.set_zlim([zc - max_range / 2, zc + max_range / 2]) def style_3d_ax(ax): """Apply clean styling to a 3D axes.""" ax.xaxis.pane.fill = True ax.yaxis.pane.fill = True ax.zaxis.pane.fill = True ax.xaxis.pane.set_facecolor((0.95, 0.95, 0.95, 1.0)) ax.yaxis.pane.set_facecolor((0.92, 0.92, 0.92, 1.0)) ax.zaxis.pane.set_facecolor((0.90, 0.90, 0.90, 1.0)) ax.xaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5}) ax.yaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5}) ax.zaxis._axinfo["grid"].update({"color": (0.8, 0.8, 0.8, 0.5), "linewidth": 0.5}) ax.tick_params(axis="both", labelsize=7, pad=1) ax.set_xlabel("X (m)", fontsize=9, labelpad=4) ax.set_ylabel("Y (m)", fontsize=9, labelpad=4) ax.set_zlabel("Z (m)", fontsize=9, labelpad=4) ax.view_init(elev=VIEW_ELEV, azim=VIEW_AZIM) def short_bridge_name(stem: str, max_len: int = 40) -> str: parts = stem.split("_") if len(parts) >= 2 and parts[0] == "bridge": bridge_id = parts[1] rest = "_".join(parts[2:]) if len(rest) > max_len - len(bridge_id) - 8: rest = rest[: max_len - len(bridge_id) - 11] + "..." return f"bridge_{bridge_id}_{rest}" return stem[:max_len] + ("..." if len(stem) > max_len else "") # Colors INLIER_COLOR = "#27ae60" # green OUTLIER_COLOR = "#bdc3c7" # light gray PLANE_COLOR = "#3498db" # blue config = BridgeProcessingConfig() if not pairs_to_plot: print("Skipping visualization (no pairs).") else: for huc_id, stem in pairs_to_plot: path = _repo_root / SOURCE_LAZ_DIR / huc_id / f"{stem}.laz" try: arrays = load_laz_arrays(path) except Exception as e: print(f"Failed to load {huc_id}/{stem}: {e}") continue result = fit_ransac_from_arrays(arrays, config) if not result.get("success"): print(f"RANSAC failed {huc_id}/{stem}: {result.get('error', 'unknown')}") continue X_local = result["X_local"] Y_local = result["Y_local"] Z = result["Z"] inlier = result["inlier_mask_full"] coef_x = result["coef_x"] coef_y = result["coef_y"] intercept = result["intercept"] n_pts = len(X_local) if MAX_POINTS_DISPLAY and n_pts > MAX_POINTS_DISPLAY: rng = np.random.RandomState(SEED) idx = rng.choice(n_pts, size=MAX_POINTS_DISPLAY, replace=False) X_local, Y_local, Z = X_local[idx], Y_local[idx], Z[idx] inlier = inlier[idx] n_inlier = int(inlier.sum()) n_outlier = int((~inlier).sum()) fig = plt.figure(figsize=(13, 9), facecolor="white") ax = fig.add_subplot(111, projection="3d") # Plot outliers first (behind) then inliers on top ax.scatter( X_local[~inlier], Y_local[~inlier], Z[~inlier], c=OUTLIER_COLOR, s=POINT_SIZE, alpha=0.35, linewidths=EDGE_LINEWIDTH, rasterized=True, label=f"Outlier ({n_outlier:,})", ) ax.scatter( X_local[inlier], Y_local[inlier], Z[inlier], c=INLIER_COLOR, s=POINT_SIZE, alpha=0.7, linewidths=EDGE_LINEWIDTH, rasterized=True, label=f"Inlier ({n_inlier:,})", ) # Fitted plane surface xmin, xmax = X_local.min(), X_local.max() ymin, ymax = Y_local.min(), Y_local.max() margin = 0.05 * max(xmax - xmin, ymax - ymin, 1.0) xg = np.linspace(xmin - margin, xmax + margin, 50) yg = np.linspace(ymin - margin, ymax + margin, 50) Xg, Yg = np.meshgrid(xg, yg) Zg = coef_x * Xg + coef_y * Yg + intercept ax.plot_surface( Xg, Yg, Zg, alpha=0.25, color=PLANE_COLOR, edgecolor=PLANE_COLOR, linewidth=0.1, ) display_name = short_bridge_name(stem) ax.set_title( f"RANSAC Deck Plane\nHUC {huc_id} | {display_name}", fontsize=12, fontweight="bold", pad=14, fontfamily="monospace", ) leg = ax.legend( loc="upper left", fontsize=9, frameon=True, fancybox=True, framealpha=0.9, edgecolor="#cccccc", markerscale=5, handletextpad=0.5, borderpad=0.6, labelspacing=0.4, ) leg.get_frame().set_facecolor("white") style_3d_ax(ax) set_equal_aspect_3d(ax) plt.tight_layout() out_path = OUTPUT_DIR / f"{OUTPUT_PREFIX}_{huc_id}_{stem}.png" fig.savefig(out_path, dpi=180, bbox_inches="tight", facecolor="white") plt.show() print(f"Saved: {out_path}")

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080206_bridge_96744504_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080206_bridge_839739315_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080206_bridge_250569218_USGS_LPC_MA_ME_MA_QL2_UTM18_L1_2015_LAS_2018.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080102_bridge_607138210_USGS_LPC_VT_Connecticut_River_2016_LAS_2018.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080102_bridge_265909373_USGS_LPC_CT_River_Lot6_Winnipesaulee_2015_LAS_2017.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01080102_bridge_458466374_USGS_LPC_VT_Connecticut_River_2016_LAS_2018.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01090002_bridge_102809372_ARRA-LFTNE_Massachusetts_2011.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01090002_bridge_305358293_ARRA-LFTNE_Massachusetts_2011.png

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/ransac_deck_plane_01090002_bridge_880364612_MA_CentralEastern_2_2021.png