code_public/lidar_rendu
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Accélérer le ray-tracing (SVF, openness) ~10x via accumulation des tangentes

Browse Source
This commit is contained in:
Antoine Jacquin
2026-09-16 20:12:00 +02:00
parent 264ccb028f
commit 0887d240f7
1 changed files with 47 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

81
lidar_pipeline/visualizations.py
View File

@ -298,10 +298,18 @@ def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=Non
recording the max upward angle (positive openness) and max downward angle
(negative openness) reached at each radius checkpoint.
Padding is done on CPU (numpy) to avoid GPU memory pressure and
pre-compiled kernel mismatches (CUDA_ERROR_NO_BINARY_FOR_GPU on sm_89).
The padded array is transferred to GPU once, then each direction is
processed and results are streamed back to CPU.
Optimisation : on accumule la TANGENTE de l'angle (dz/dist) au lieu de
l'angle lui-même — atan étant strictement croissante, max(angles) =
atan(max(tangentes)). L'arctan (coûteuse, pleine image) n'est donc plus
appliquée qu'aux checkpoints de rayon, pas à chaque pas de rayon.
Un seul couple de max cumulés est maintenu, snapshoté à chaque checkpoint
(les rayons étant emboîtés, chaque checkpoint réutilisait avant le même
calcul 3 fois). fmax ignore les NaN du padding : plus de nan_to_num/where.
Padding on CPU (numpy) to avoid GPU memory pressure and pre-compiled
kernel mismatches (CUDA_ERROR_NO_BINARY_FOR_GPU on sm_89). The padded
array is transferred to GPU once, then each direction is processed and
results are streamed back to CPU.
Args:
dem: CPU numpy array — filled DEM (no NaN), shape (rows, cols).
@ -338,6 +346,12 @@ def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=Non
# Free the CPU copy — we don't need it anymore
del padded_np
# Checkpoints triés par pas : (step, r_idx). Les snapshots sont pris quand
# le pas courant atteint le pas du checkpoint — les rayons ne dépassent
# donc pas le plus grand checkpoint demandé (équivalent au break d'avant).
checkpoints = sorted((radii_steps[r_idx], r_idx) for r_idx in range(n_radii))
last_step = checkpoints[-1][0]
# Process one direction at a time to limit GPU memory.
# Store results as flat CPU arrays — transfer back to GPU at the end.
pos_results = [None] * n_dirs
@ -346,9 +360,9 @@ def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=Non
for d_idx in range(n_dirs):
ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]
# Pre-compute valid steps for this direction
# Pre-compute valid steps for this direction (jusqu'au dernier checkpoint)
valid_steps = []
for step in range(1, max_dist + 1):
for step in range(1, last_step + 1):
px = int(round(ddx * step))
py = int(round(ddy * step))
dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
@ -356,11 +370,11 @@ def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=Non
continue
valid_steps.append((step, px, py, dist_m))
# Running max angles per radius
running_pos = xp.zeros((n_radii, rows, cols))
running_neg = xp.zeros((n_radii, rows, cols))
# Track which radius checkpoints have been passed
radii_remaining = set(range(n_radii))
# Max cumulé des tangentes (float32 : moitié de VRAM vs float64)
running_pos = xp.zeros((rows, cols), dtype=np.float32)
running_neg = xp.zeros((rows, cols), dtype=np.float32)
snapshots = {}
cp_queue = list(checkpoints)
for step, px, py, dist_m in valid_steps:
# Slice from padded array, subtract original dem
@ -369,34 +383,33 @@ def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=Non
elev_diff = view - dem
del view # free slice reference
# Positive: angle to terrain above viewer
pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m)
# Negative: angle to terrain below viewer
neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m)
# Tangentes des angles (positive : terrain au-dessus, négative : en
# dessous). fmax propage le non-NaN : le bord de padding ne compte
# pas, comme avec l'ancien where(isnan) — en une seule opération.
running_pos = xp.fmax(running_pos,
xp.maximum(elev_diff, 0) / dist_m)
running_neg = xp.fmax(running_neg,
xp.maximum(-elev_diff, 0) / dist_m)
del elev_diff # free intermediate
# Update running max for all radius checkpoints still active
for r_idx in radii_remaining:
pos_angle_safe = xp.nan_to_num(pos_angle, nan=0)
neg_angle_safe = xp.nan_to_num(neg_angle, nan=0)
running_pos[r_idx] = xp.where(xp.isnan(pos_angle), running_pos[r_idx],
xp.maximum(running_pos[r_idx], pos_angle_safe))
running_neg[r_idx] = xp.where(xp.isnan(neg_angle), running_neg[r_idx],
xp.maximum(running_neg[r_idx], neg_angle_safe))
# Check which radii have been passed
new_remaining = set()
for r_idx in radii_remaining:
if step < radii_steps[r_idx]:
new_remaining.add(r_idx)
radii_remaining = new_remaining
if not radii_remaining:
# Snapshot du checkpoint atteint : conversion en angle UNE fois
while cp_queue and step >= cp_queue[0][0]:
_, r_idx = cp_queue.pop(0)
snapshots[r_idx] = (xp.arctan(running_pos),
xp.arctan(running_neg))
if not cp_queue:
break
# Checkpoints jamais atteints (steps invalides) : état final du balayage
while cp_queue:
_, r_idx = cp_queue.pop(0)
snapshots[r_idx] = (xp.arctan(running_pos),
xp.arctan(running_neg))
# Store results on CPU, free GPU memory before next direction
pos_results[d_idx] = to_cpu(running_pos)
neg_results[d_idx] = to_cpu(running_neg)
del running_pos, running_neg
pos_results[d_idx] = to_cpu(xp.stack([snapshots[r][0] for r in range(n_radii)]))
neg_results[d_idx] = to_cpu(xp.stack([snapshots[r][1] for r in range(n_radii)]))
del running_pos, running_neg, snapshots
gpu_cleanup()
# Free the large padded array