"""Independent geometric and continuous-beam LiDAR references; public LESS runs.""" from pathlib import Path import argparse import json import less import numpy as np from scipy.special import ndtr, roots_legendre DATA = Path(__file__).resolve().parents[1] / "data/lidar" C = 299792458.0 BACKENDS = ("optix", "vulkan", "embree") def edges(config): count = int((config["range_max_m"]-config["range_min_m"]) / config["range_resolution_m"]) + 1 return config["range_min_m"] + np.arange(count+1)*config["range_resolution_m"] def sigma_range(config): return config["pulse_fwhm_ns"]*1e-9*C/(2*np.sqrt(8*np.log(2))) def pulse_bins(config, ranges): return np.diff(ndtr((edges(config)[None, :]-np.atleast_1d(ranges)[:, None]) / sigma_range(config)), axis=1) def width_ns(config, wave): if not wave.sum(): return 0. centres = (edges(config)[1:]+edges(config)[:-1])/2 mean = np.sum(centres*wave)/wave.sum() variance = np.sum((centres-mean)**2*wave)/wave.sum() return float(np.sqrt(variance)*np.sqrt(8*np.log(2))*2/C*1e9) def compact(wave): nonzero = np.flatnonzero(wave > 0) if not len(nonzero): return dict(first_bin=0, values=[]) first, last = max(0, int(nonzero[0])-1), min(len(wave), int(nonzero[-1])+2) return dict(first_bin=first, values=wave[first:last].tolist()) def expand(record, count): result = np.zeros(count) first = record["first_bin"] result[first:first+len(record["values"])] = record["values"] return result def direction(view): theta, phi = np.deg2rad([view["zenith_degrees"], view["azimuth_degrees"]]) return np.array([np.sin(theta)*np.sin(phi), -np.sin(theta)*np.cos(phi), np.cos(theta)]) def geometric_reference(config, view): """Closed ray/plane intersections; no mesh parsing or native query reuse.""" origin, ray = np.array(view["position_m"], float), direction(view) planes = [("ground", [0, 0, 1], 0, [0, config["size_m"][0], 0, config["size_m"][1]]), ("roof", [0, 0, 1], config["roof_height_m"], config["roof_bounds_xy_m"]), ("slope", [-config["slope_dz_dx"], 0, 1], config["slope_z_intercept_m"], config["slope_bounds_xy_m"])] hits = [] for name, normal, offset, bounds in planes: normal = np.array(normal, float) denominator = normal @ ray if abs(denominator) < 1e-12: continue distance = (offset-normal @ origin)/denominator point = origin+distance*ray if distance <= 0 or not (bounds[0] <= point[0] <= bounds[1] and bounds[2] <= point[1] <= bounds[3]): continue cosine = abs(denominator)/np.linalg.norm(normal) energy = (config["pulse_energy_J"]*config["reflectance"][name]*cosine / np.pi * config["receiver_area_m2"]/distance**2) hits.append(dict(surface=name, range_m=float(distance), point=point.tolist(), cosine=float(cosine), energy=float(energy), delay_ns=2*distance/C*1e9)) return min(hits, key=lambda h: h["range_m"]) if hits else None def finite_reference(config, offset, order): """Exact azimuth fractions and independent radial Gauss-Legendre quadrature.""" nodes, weights = roots_legendre(order) cut = min(1.0, abs(offset)) intervals = [(0, cut), (cut, 1)] if 0 < cut < 1 else [(0, 1)] wave = np.zeros(len(edges(config))-1) energies = np.zeros(2) a = np.tan(config["finite_half_angle_mrad"]*1e-3) k = -np.log(config["beam_edge_fraction"]) norm = 2*k/(1-np.exp(-k)) for low, high in intervals: radial = (nodes+1)*(high-low)/2+low quadrature = weights*(high-low)/2 * norm*radial*np.exp(-k*radial**2) roof_fraction = np.arccos(np.clip(-offset/radial, -1, 1))/np.pi for index, (surface, fraction) in enumerate(( ("roof", roof_fraction), ("ground", 1-roof_fraction))): height = config["sensor_height_m"] - ( config["roof_height_m"] if surface == "roof" else 0) cosine = 1/np.sqrt(1+(a*radial)**2) distance = height/cosine energy = (config["pulse_energy_J"]*config["receiver_area_m2"] * config["reflectance"][surface]/np.pi/height**2 * cosine**3) weighted = quadrature*fraction*energy energies[index] += weighted.sum() wave += weighted @ pulse_bins(config, distance) return wave, energies def scene_for(config, backend): scene = less.Scene(backend=backend) scene.size = tuple(config["size_m"]) scene.repetitive = False scene.terrain = less.Terrain(property=less.Lambertian( reflectance=config["reflectance"]["ground"])) scene.illumination = less.Illumination(source=less.Sun(irradiance=0), atmosphere=less.NoAtmosphere()) for name in ("roof", "slope"): obj = less.Object(name, mesh=str(DATA/f"{name}.obj")) obj.set_property(less.Lambertian(reflectance=config["reflectance"][name])) scene.add(obj, positions=[(0, 0, 0)]) return scene def occlusion_scene(config, backend): settings = config["occlusion"] scene = less.Scene(backend=backend) scene.size = tuple(config["size_m"]) scene.repetitive = False scene.terrain = less.Terrain(property=less.Lambertian( reflectance=settings["reflectance"]["ground"])) scene.illumination = less.Illumination(source=less.Sun(irradiance=0), atmosphere=less.NoAtmosphere()) for mesh, surface in (("aperture", "upper"), ("lower", "lower")): obj = less.Object(surface, mesh=str(DATA/f"{mesh}.obj")) obj.set_property(less.Lambertian(reflectance=settings["reflectance"][surface])) scene.add(obj, positions=[(0, 0, 0)]) return scene def occlusion_geometry(config, view): """Nearest admissible plane intersection, including the upper slit.""" settings = config["occlusion"] origin, ray = np.array(view["position_m"]), direction(view) for surface, height, bounds in ( ("upper", settings["upper_height_m"], settings["upper_bounds_xy_m"]), ("lower", settings["lower_height_m"], settings["lower_bounds_xy_m"]), ("ground", 0, [0, config["size_m"][0], 0, config["size_m"][1]])): distance = (height-origin[2])/ray[2] point = origin+distance*ray if not (bounds[0] <= point[0] <= bounds[1] and bounds[2] <= point[1] <= bounds[3]): continue if surface == "upper" and settings["opening_x_m"][0] < point[0] < settings["opening_x_m"][1]: continue energy = (config["pulse_energy_J"]*config["receiver_area_m2"] *settings["reflectance"][surface]*abs(ray[2])/np.pi/distance**2) return dict(surface=surface, range_m=float(distance), point=point.tolist(), energy=float(energy), delay_ns=float(2*distance/C*1e9)) raise ValueError("The validation view must intersect one of the bounded planes") def occlusion_reference(config, view, order): """Continuous disk integration split at exact plane/slit angular crossings. This reference integrates polar coordinates, independently of the simulator's Cartesian beam grid. All views lie in the XZ plane; the footprint stays away from the finite panels' Y edges. """ settings = config["occlusion"] origin, axis = np.array(view["position_m"], float), direction(view) assert abs(axis[1]) < 1e-12 and axis[2] < 0 sine, cosine = axis[0], -axis[2] tangent = np.tan(config["finite_half_angle_mrad"]*1e-3) surfaces = ("upper", "lower", "ground") heights = np.array([settings["upper_height_m"], settings["lower_height_m"], 0]) gaps = origin[2]-heights cuts = [] for index, boundaries in ((0, settings["opening_x_m"]), (1, [settings["lower_bounds_xy_m"][0]])): for boundary in boundaries: delta = boundary-origin[0] cut = (delta*cosine-gaps[index]*sine)/(tangent*(gaps[index]*cosine+delta*sine)) if -1 < cut < 1: cuts.append(float(cut)) radial_edges = sorted(set([0., 1.]+[abs(cut) for cut in cuts])) radial_nodes, radial_weights = roots_legendre(order) angular_nodes, angular_weights = roots_legendre(max(16, order//4)) k = -np.log(config["beam_edge_fraction"]) waves = np.zeros((3, len(edges(config))-1)) energies = np.zeros(3) for low, high in zip(radial_edges, radial_edges[1:]): radii = low+(radial_nodes+1)*(high-low)/2 weights = radial_weights*(high-low)/2*2*k*radii*np.exp(-k*radii**2)/(1-np.exp(-k)) for radius, weight in zip(radii, weights): angles = sorted([0., np.pi]+[np.arccos(cut/radius) for cut in cuts if abs(cut)settings["upper_bounds_xy_m"][0]) & (x_uppersettings["upper_bounds_xy_m"][2]) & (y_uppersettings["lower_bounds_xy_m"][2]) & (y_lowersettings["opening_x_m"][1]) lower = x_lower>=settings["lower_bounds_xy_m"][0] component = np.where(upper, 0, np.where(lower, 1, 2)) for index in np.unique(component): selected = component==index ranges = gaps[index]/z_cosine[selected] energy = (config["pulse_energy_J"]*config["receiver_area_m2"] *settings["reflectance"][surfaces[index]]/np.pi*z_cosine[selected]/ranges**2) # The second half of the disk has identical distances. weighted = weight*angular_weights[selected]*(end-start)/(2*np.pi)*energy energies[index] += weighted.sum() waves[index] += weighted @ pulse_bins(config, ranges) return waves, energies def run_occlusion(config): settings, limits = config["occlusion"], config["thresholds"] references, records, convergence = [], [], [] count = len(edges(config))-1 centres = (edges(config)[:-1]+edges(config)[1:])/2 for view in settings["finite_views"]: components, energy = occlusion_reference(config, view, settings["reference_orders"][-1]) coarse, _ = occlusion_reference(config, view, settings["reference_orders"][-2]) change = float(abs(components-coarse).sum()/energy.sum()) assert change <= limits["reference_quadrature_relative_change"] references.append(dict(view=view["name"], waveform=compact(components.sum(axis=0)), components=[compact(wave) for wave in components], energy=energy.tolist(), quadrature_change=change)) for backend in BACKENDS: scene = occlusion_scene(config, backend) for view in settings["zero_views"]: exact = occlusion_geometry(config, view) target = exact["energy"]*pulse_bins(config, exact["range_m"])[0] direct = scene.simulate(survey_for(config, view, 0, 1), product="fast_point_cloud") direct_identity = identity(direct, backend, "lidar_point_cloud") assert direct.num_points == 1 position_error = float(np.linalg.norm(direct.points[0]-exact["point"])) energy_error = abs(float(direct.intensity[0])/exact["energy"]-1) timeline = direct.return_timeline delay_error = abs((timeline["min_arrival_seconds"]- timeline["min_emission_seconds"])*1e9-exact["delay_ns"]) for samples in config["zero_beam_samples"]: product = scene.simulate(survey_for(config, view, 0, samples), product="waveform") execution = identity(product, backend, "lidar_waveform") wave = np.asarray(product.waveform["data"], float)[0, :, 0] error = float(abs(wave-target).sum()/target.sum()) points = point_records(product, np.array(view["position_m"])) range_error = abs(points[0]["range_m"]-exact["range_m"]) if len(points)==1 else None width_error = abs(width_ns(config, wave)-width_ns(config, target)) passed = (position_error<=limits["direct_position_error_m"] and energy_error<=limits["direct_energy_relative_error"] and delay_error<=limits["direct_delay_error_ns"] and error<=limits["zero_waveform_l1_relative_error"] and width_error<=limits["pulse_width_error_ns"] and range_error is not None and range_error<=config["range_resolution_m"]*.5+3e-5) if samples==config["zero_beam_samples"][-1]: records.append(dict(backend=backend, mode="zero", view=view["name"], reference=exact, waveform=compact(wave), points=points, direct_points=direct.points.tolist(), direct_execution=direct_identity, execution=execution, position_error_m=position_error, energy_relative_error=energy_error, delay_error_ns=delay_error, waveform_relative_error=error, return_range_error_m=range_error, width_error_ns=width_error, passed=bool(passed))) convergence.append(dict(backend=backend, mode="zero", view=view["name"], samples=samples, waveform_relative_error=error)) for view, reference in zip(settings["finite_views"], references): target = expand(reference["waveform"], count) expected_energy = np.array(reference["energy"]) components = np.array([expand(wave, count) for wave in reference["components"]]) nominal = (config["sensor_height_m"]-np.array([ settings["upper_height_m"], settings["lower_height_m"], 0]))/abs(direction(view)[2]) splits = (nominal[:-1]+nominal[1:])/2 labels = np.searchsorted(splits, centres) masks = [labels==index for index in range(3)] active = expected_energy>0 expected_ranges = [float(centres[np.argmax(wave)]) for wave, present in zip(components, active) if present] previous = None for samples in settings["finite_beam_samples"]: product = scene.simulate(survey_for(config, view, config["finite_half_angle_mrad"], samples), product="waveform") execution = identity(product, backend, "lidar_waveform") wave = np.asarray(product.waveform["data"], float)[0, :, 0] energy = np.array([wave[mask].sum() for mask in masks]) error = float(abs(wave-target).sum()/target.sum()) energy_error = float(max(abs(energy[active]/expected_energy[active]-1))) leakage = float(energy[~active].sum()/target.sum()) change = None if previous is None else float(abs(wave-previous).sum()/target.sum()) points = point_records(product, np.array(view["position_m"])) actual_ranges = [point["range_m"] for point in points] range_error = (max(abs(a-b) for a,b in zip(actual_ranges, expected_ranges)) if len(points)==len(expected_ranges) else None) separation_error = (float(max(abs(np.diff(actual_ranges)-np.diff(expected_ranges)))) if len(points)==len(expected_ranges)>1 else 0.) width_error = max(abs(width_ns(config, wave*mask)-width_ns(config, ref)) for mask, ref, present in zip(masks, components, active) if present) numbering = ([p["return_number"] for p in points]==list(range(1,len(points)+1)) and all(p["num_returns"]==len(points) for p in points)) passed = (error<=limits["finite_waveform_l1_relative_error"] and energy_error<=limits["finite_component_energy_relative_error"] and leakage<=limits["blocked_energy_fraction"] and width_error<=limits["pulse_width_error_ns"] and numbering and range_error is not None and range_error<=limits["waveform_return_range_error_m"] and separation_error<=limits["echo_separation_error_m"] and change is not None and change<=limits["final_budget_waveform_relative_change"]) convergence.append(dict(backend=backend, mode="finite", view=view["name"], samples=samples, waveform_relative_error=error, component_energy_relative_error=energy_error, budget_change=change)) if samples==settings["finite_beam_samples"][-1]: records.append(dict(backend=backend, mode="finite", view=view["name"], waveform=compact(wave), energy=energy.tolist(), points=points, execution=execution, waveform_relative_error=error, component_energy_relative_error=energy_error, blocked_energy_fraction=leakage, width_error_ns=width_error, return_range_error_m=range_error, echo_separation_error_m=separation_error, budget_change=change, passed=bool(passed))) previous = wave print(backend, view["name"], error, energy_error, change, passed, flush=True) return dict(references=references, records=records, convergence=convergence, passed=all(row["passed"] for row in records)) def survey_for(config, view, divergence, samples): return less.LiDARSurvey(scanner=less.Scanner( beam_divergence=divergence, beam_samples=samples, angular_resolution=.1, wavelengths=[config["wavelength_nm"]], pulse_energy=config["pulse_energy_J"], sensor_area=config["receiver_area_m2"], pulse_fwhm=config["pulse_fwhm_ns"]), platform=less.Tripod(height=config["sensor_height_m"]), legs=[less.ScanPosition(*view["position_m"], azimuth_range=(view["azimuth_degrees"], view["azimuth_degrees"]+.01), zenith_range=(view["zenith_degrees"], view["zenith_degrees"]+.01))], range_min=config["range_min_m"], range_max=config["range_max_m"], range_resolution=config["range_resolution_m"]) def identity(product, backend, workflow): value = product.metadata["execution"] assert value["backend"] == backend and value["native"] is True assert value["workflow"] == workflow and value["executor"] == f"{backend}.slang.lidar" return {key: value[key] for key in ("backend", "workflow", "executor", "native")} def point_records(product, origin): if product.num_points == 0: return [] return [dict(position=p.tolist(), range_m=float(np.linalg.norm(p-origin)), return_number=int(number), num_returns=int(count)) for p, number, count in zip(product.points, product.return_number, product.num_returns)] def run(config): thresholds = config["thresholds"] references, records = [], [] for offset in config["edge_offsets_in_roof_radius"]: wave, energy = finite_reference(config, offset, config["reference_orders"][-1]) coarse, _ = finite_reference(config, offset, config["reference_orders"][-2]) change = float(abs(wave-coarse).sum()/wave.sum()) assert change <= thresholds["reference_quadrature_relative_change"] references.append(dict(offset=offset, waveform=compact(wave), energy=energy.tolist(), quadrature_change=change)) for backend in BACKENDS: scene = scene_for(config, backend) for view in config["zero_views"]: exact = geometric_reference(config, view) survey = survey_for(config, view, 0, 1) direct = scene.simulate(survey, product="fast_point_cloud") execution = identity(direct, backend, "lidar_point_cloud") position_error = energy_error = delay_error = 0. if exact: assert direct.num_points == 1 position_error = float(np.linalg.norm(direct.points[0]-exact["point"])) energy_error = abs(float(direct.intensity[0])/exact["energy"]-1) timeline = direct.return_timeline actual_delay = (timeline["min_arrival_seconds"] - timeline["min_emission_seconds"])*1e9 delay_error = abs(actual_delay-exact["delay_ns"]) else: assert direct.num_points == 0 target = (exact["energy"]*pulse_bins(config, exact["range_m"])[0] if exact else np.zeros(len(edges(config))-1)) for samples in config["zero_beam_samples"]: survey.scanner.beam_samples = samples product = scene.simulate(survey, product="waveform") wave = np.asarray(product.waveform["data"], float)[0, :, 0] wave_execution = identity(product, backend, "lidar_waveform") error = float(abs(wave-target).sum()/target.sum()) if exact else float(abs(wave).sum()) points = point_records(product, np.array(view["position_m"])) range_error = (max(abs(p["range_m"]-exact["range_m"]) for p in points) if exact and points else 0.) width_error = abs(width_ns(config, wave)-width_ns(config, target)) passed = (position_error <= thresholds["direct_position_error_m"] and energy_error <= thresholds["direct_energy_relative_error"] and delay_error <= thresholds["direct_delay_error_ns"] and width_error <= thresholds["pulse_width_error_ns"] and error <= thresholds["zero_waveform_l1_relative_error"] and len(points) == (1 if exact else 0) and range_error <= config["range_resolution_m"]*.5+3e-5) row = dict(backend=backend, mode="zero", view=view["name"], samples=samples, reference=exact, direct_points=direct.points.tolist(), position_error_m=position_error, energy_relative_error=energy_error, waveform_relative_error=error, delay_error_ns=delay_error, width_error_ns=width_error, return_range_error_m=range_error, points=points, direct_execution=execution, execution=wave_execution, passed=bool(passed)) if samples == config["zero_beam_samples"][-1]: row["waveform"] = compact(wave) else: for key in ("passed", "reference", "direct_points", "direct_execution", "points"): row.pop(key) records.append(row) print(backend, view["name"], position_error, energy_error, error, passed, flush=True) radius = ((config["sensor_height_m"]-config["roof_height_m"]) * np.tan(config["finite_half_angle_mrad"]*1e-3)) for reference in references: offset = reference["offset"] x = config["roof_bounds_xy_m"][0]+offset*radius view = dict(position_m=[x, 2.5, config["sensor_height_m"]], zenith_degrees=180, azimuth_degrees=0) target = expand(reference["waveform"], len(edges(config))-1) centres = (edges(config)[:-1]+edges(config)[1:])/2 expected_energies = np.array(reference["energy"]) split = (2*config["sensor_height_m"]-config["roof_height_m"])/2 masks = [centres < split, centres >= split] expected_ranges = [float(centres[np.flatnonzero(mask)[np.argmax(target[mask])]]) for mask, energy in zip(masks, expected_energies) if energy > 0] previous = None for samples in config["finite_beam_samples"]: survey = survey_for(config, view, config["finite_half_angle_mrad"], samples) product = scene.simulate(survey, product="waveform") execution = identity(product, backend, "lidar_waveform") wave = np.asarray(product.waveform["data"], float)[0, :, 0] error = float(abs(wave-target).sum()/target.sum()) energies = np.array([wave[mask].sum() for mask in masks]) nonzero = expected_energies > 0 energy_error = float(np.max(abs(energies[nonzero]/expected_energies[nonzero]-1))) width_error = max(abs(width_ns(config, wave*mask)-width_ns(config, target*mask)) for mask, active in zip(masks, nonzero) if active) change = None if previous is None else float(abs(wave-previous).sum()/target.sum()) points = point_records(product, np.array(view["position_m"])) range_error = (max(abs(p["range_m"]-r) for p,r in zip(points, expected_ranges)) if len(points) == len(expected_ranges) else None) passed = (error <= thresholds["finite_waveform_l1_relative_error"] and energy_error <= thresholds["finite_component_energy_relative_error"] and width_error <= thresholds["pulse_width_error_ns"] and range_error is not None and range_error <= thresholds["waveform_return_range_error_m"] and (samples != config["finite_beam_samples"][-1] or change <= thresholds["final_budget_waveform_relative_change"])) row = dict(backend=backend, mode="finite", offset=offset, samples=samples, waveform_relative_error=error, energy=energies.tolist(), component_energy_relative_error=energy_error, budget_change=change, width_error_ns=width_error, return_range_error_m=range_error, points=points, origin=view["position_m"], execution=execution, passed=bool(passed)) if samples == config["finite_beam_samples"][-1]: row["waveform"] = compact(wave) else: row.pop("passed") row.pop("points") records.append(row) previous = wave print(backend, "finite", offset, error, energy_error, change, passed, flush=True) final = [r for r in records if r["samples"] == config[ "zero_beam_samples" if r["mode"] == "zero" else "finite_beam_samples"][-1]] return dict(input=config, finite_reference=references, records=records, passed=all(r["passed"] for r in final)) def plot_results(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.size": 10, "svg.fonttype": "none"}) config = result["input"] grid = edges(config) centres = (grid[1:]+grid[:-1])/2 colors = dict(optix="#0072b2", vulkan="#d55e00", embree="#009e73") markers = dict(optix="o", vulkan="x", embree="+") final = [r for r in result["records"] if "waveform" in r] fig, axes = plt.subplots(2, 2, figsize=(11, 7), layout="constrained") for ax, name in zip(axes[0], ("ground_oblique", "slope_nadir")): rows = [r for r in final if r["mode"] == "zero" and r["view"] == name] reference = rows[0]["reference"] target = reference["energy"]*pulse_bins(config, reference["range_m"])[0] ax.plot(centres, target*1e6, "k-", label="Analytic") for row in rows: backend = row["backend"] ax.plot(centres, expand(row["waveform"], len(centres))*1e6, color=colors[backend], marker=markers[backend], ls="none", ms=4, label=backend) ax.set(xlim=(reference["range_m"]-.3, reference["range_m"]+.3), title=name.replace("_", " ").capitalize(), xlabel="Range (m)", ylabel="Energy per range bin (µJ)") ax.grid(alpha=.2) ax = axes[1, 0] reference = next(r for r in result["finite_reference"] if r["offset"] == -.35) ax.plot(centres, expand(reference["waveform"], len(centres))*1e6, "k-", label="Continuous-beam reference") for row in final: if row["mode"] == "finite" and row["offset"] == -.35: b = row["backend"] ax.plot(centres, expand(row["waveform"], len(centres))*1e6, color=colors[b], marker=markers[b], ls="none", ms=3, label=b) ax.set(xlim=(5.6, 10.4), title="Finite beam at board edge: two returns", xlabel="Range (m)", ylabel="Energy per range bin (µJ)") ax = axes[1, 1] refs = result["finite_reference"] offsets = [r["offset"] for r in refs] fractions = [r["energy"][0]/sum(r["energy"]) for r in refs] ax.plot(offsets, fractions, "k-", label="Continuous-beam reference") for b in BACKENDS: rows = [r for r in final if r["mode"] == "finite" and r["backend"] == b] ax.plot(offsets, [r["energy"][0]/sum(r["energy"]) for r in rows], color=colors[b], marker=markers[b], ls="none", label=b) ax.set(title="Occlusion and return-energy partition", xlabel="Beam-centre offset / footprint radius at board", ylabel="Board fraction of received energy", ylim=(-.04, 1.04)) for ax in axes.flat: ax.grid(alpha=.2) ax.legend(frameon=False, fontsize=8) fig.savefig(output/"waveforms.svg") fig.savefig(output/"waveforms.png", dpi=140) plt.close(fig) fig, axes = plt.subplots(1, 2, figsize=(11, 3.8), layout="constrained") for b in BACKENDS: rows = [r for r in result["records"] if r["mode"] == "finite" and r["backend"] == b] for ax, metric in zip(axes, ("waveform_relative_error", "component_energy_relative_error")): errors = [max(r[metric] for r in rows if r["samples"] == n)*100 for n in config["finite_beam_samples"]] ax.plot(config["finite_beam_samples"], errors, color=colors[b], marker=markers[b], label=b) for ax, title in zip(axes, ("Waveform: worst of five footprint positions", "Echo energy: worst of both surfaces")): ax.axhline(.5, color="k", ls="--", label="0.5% threshold") ax.set(xscale="log", yscale="log", xlabel="Beam grid intervals per dimension", ylabel="Relative error (%)", title=title) ax.grid(alpha=.2) ax.legend(frameon=False, fontsize=8) fig.savefig(output/"convergence.svg") fig.savefig(output/"convergence.png", dpi=140) plt.close(fig) def plot_occlusion(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.size": 10, "svg.fonttype": "none"}) config, data = result["input"], result["occlusion"] centres = (edges(config)[:-1]+edges(config)[1:])/2 colors = dict(optix="#0072b2", vulkan="#d55e00", embree="#009e73") markers = dict(optix="o", vulkan="x", embree="+") fig, axes = plt.subplots(2, 2, figsize=(11, 7), layout="constrained") for ax, name in zip(axes.flat[:3], ("three_returns", "oblique_lower", "oblique_ground")): reference = next(row for row in data["references"] if row["view"]==name) target = expand(reference["waveform"], len(centres)) ax.plot(centres, target*1e6, "k-", label="Continuous-beam reference") for backend in BACKENDS: row = next(row for row in data["records"] if row["backend"]==backend and row["mode"]=="finite" and row["view"]==name) ax.plot(centres, expand(row["waveform"], len(centres))*1e6, color=colors[backend], marker=markers[backend], ls="none", ms=3, label=backend) ax.set(xlim=(3.5, 11.5), title=name.replace("_", " ").capitalize(), xlabel="Range (m)", ylabel="Energy per range bin (µJ)") ax.legend(frameon=False, fontsize=8) ax = axes.flat[3] for backend in BACKENDS: levels = [row for row in data["convergence"] if row["backend"]==backend and row["mode"]=="finite"] budgets = config["occlusion"]["finite_beam_samples"] ax.plot(budgets, [max(row["waveform_relative_error"] for row in levels if row["samples"]==budget)*100 for budget in budgets], color=colors[backend], marker=markers[backend], label=backend) ax.axhline(.5, color="k", ls="--", label="0.5% threshold") ax.set(xscale="log", yscale="log", xlabel="Beam grid intervals per dimension", ylabel="Maximum waveform error (%)", title="Five occlusion views") ax.legend(frameon=False, fontsize=8) for ax in axes.flat: ax.grid(alpha=.2) fig.savefig(output/"occlusion.svg") fig.savefig(output/"occlusion.png", dpi=140) plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, default=Path("build/tutorial-lidar")) parser.add_argument("--plot-only", type=Path) parser.add_argument("--occlusion-only", action="store_true") args = parser.parse_args() args.output.mkdir(parents=True, exist_ok=True) if args.plot_only: result = json.loads(args.plot_only.read_text()) plot_results(result, args.output) if "occlusion" in result: plot_occlusion(result, args.output) return config = json.loads((DATA/"input.json").read_text()) if args.occlusion_only: result = json.loads((DATA/"results.json").read_text()) assert {k:v for k,v in result["input"].items() if k not in ("occlusion", "thresholds")} == { k:v for k,v in config.items() if k not in ("occlusion", "thresholds")} assert all(config["thresholds"][key]==value for key,value in result["input"]["thresholds"].items()) else: result = run(config) original_passed = result.get("basic_passed", result["passed"]) result["input"] = config result["occlusion"] = run_occlusion(config) result["basic_passed"] = original_passed result["passed"] = bool(original_passed and result["occlusion"]["passed"]) (args.output/"results.json").write_text(json.dumps(result, indent=2, allow_nan=False)) raise SystemExit(0 if result["passed"] else 1) if __name__ == "__main__": main()