Algorithms
This section describes the computational algorithms implemented in L-SURF.
Ray-Surface Intersection
Overview
Ray-surface intersection is the most critical operation in raytracing. L-SURF implements specialized algorithms for each surface type.
Ray Parametrization
A ray is defined by:
where:
\(\vec{r}_0\) = ray origin
\(\vec{d}\) = ray direction (normalized)
\(t\) = distance parameter (\(t \geq 0\))
Planar Surface Intersection
For a plane defined by point \(\vec{p}\) and normal \(\vec{n}\):
Algorithm:
Compute ray-plane angle:
\[\cos\alpha = \vec{d} \cdot \vec{n}\]If \(|\cos\alpha| < \epsilon\), ray is parallel (no intersection)
Compute intersection distance:
\[t = \frac{(\vec{p} - \vec{r}_0) \cdot \vec{n}}{\vec{d} \cdot \vec{n}}\]If \(t < 0\), intersection is behind ray origin
Intersection point: \(\vec{r}_{\text{hit}} = \vec{r}_0 + t \vec{d}\)
Complexity: \(O(1)\) - constant time
Numerical Stability:
Use double precision for \(t\) calculation
Threshold \(\epsilon \approx 10^{-10}\) for parallelism test
Offset intersection point by \(\epsilon \vec{n}\) to avoid self-intersection
Spherical Surface Intersection
For a sphere with center \(\vec{c}\) and radius \(R\):
Implicit equation:
Algorithm:
Transform to sphere-local coordinates: \(\vec{o} = \vec{r}_0 - \vec{c}\)
Solve quadratic equation \(at^2 + bt + c = 0\):
\[\begin{split}a &= \vec{d} \cdot \vec{d} = 1 \\ b &= 2 \vec{o} \cdot \vec{d} \\ c &= \vec{o} \cdot \vec{o} - R^2\end{split}\]Compute discriminant: \(\Delta = b^2 - 4ac\)
If \(\Delta < 0\), no intersection
Solutions:
\[t_{\pm} = \frac{-b \pm \sqrt{\Delta}}{2a}\]Choose appropriate root:
For exterior ray: \(t_1 = \min(t_+, t_-)\) (first intersection)
For interior ray: \(t_2 = \max(t_+, t_-)\) (exit point)
Complexity: \(O(1)\) - constant time
Numerical Stability:
Use numerically stable quadratic solver:
\[ \begin{align}\begin{aligned}q = -\frac{1}{2}(b + \text{sign}(b)\sqrt{\Delta})\\t_1 = q/a, \quad t_2 = c/q\end{aligned}\end{align} \]
Wave Surface Intersection
For parametric wave surfaces (e.g., Gerstner waves), no closed-form solution exists. We use iterative root-finding.
Algorithm (Newton-Raphson):
Define implicit function:
\[f(t) = z(\vec{r}_0 + t\vec{d}) - (z_0 + t d_z)\]where \(z(\vec{r})\) is the wave surface height.
Initial guess from planar approximation:
\[t_0 = \frac{-z_0}{d_z}\]Iterate until convergence:
\[t_{n+1} = t_n - \frac{f(t_n)}{f'(t_n)}\]Derivative:
\[f'(t) = \nabla z \cdot \vec{d} - d_z\]Convergence criterion: \(|f(t_n)| < \epsilon\) (typically \(\epsilon = 10^{-6}\) m)
Maximum iterations: 20 (typically converges in 3-5 iterations)
Complexity: \(O(N_{\text{iter}} \cdot N_{\text{waves}})\) where \(N_{\text{waves}}\) is number of wave components
Robustness:
Fallback to bisection if Newton fails to converge
Check for multiple intersections in rough seas
Use adaptive step size for steep waves
GPU Implementation
For GPU execution, wave surface intersection is parallelized:
@cuda.jit
def intersect_wave_surface_kernel(ray_positions, ray_directions,
wave_params, intersection_points):
idx = cuda.grid(1)
if idx < ray_positions.shape[0]:
# Newton-Raphson iteration in parallel
t = initial_guess(ray_positions[idx], ray_directions[idx])
for iteration in range(MAX_ITERATIONS):
pos = ray_positions[idx] + t * ray_directions[idx]
z_surf = evaluate_wave_height(pos, wave_params)
f = z_surf - pos[2]
if abs(f) < EPSILON:
break
df_dt = evaluate_wave_gradient(pos, wave_params).dot(
ray_directions[idx]
) - ray_directions[idx][2]
t -= f / df_dt
intersection_points[idx] = ray_positions[idx] + t * ray_directions[idx]
Reflection Vector Calculation
Law of Reflection
Given incident direction \(\vec{d}_i\) and surface normal \(\vec{n}\):
Algorithm:
Ensure normal points toward incident medium:
if dot(incident_dir, normal) > 0: normal = -normal
Compute projection: \(p = \vec{d}_i \cdot \vec{n}\)
Reflect: \(\vec{d}_r = \vec{d}_i - 2p\vec{n}\)
Normalize (if needed): \(\vec{d}_r = \vec{d}_r / |\vec{d}_r|\)
Verification:
Check angle equality:
Refraction Vector Calculation
Snell’s Law in Vector Form
Given incident direction \(\vec{d}_i\), surface normal \(\vec{n}\), and refractive indices \(n_1\), \(n_2\):
where:
Algorithm:
Compute refractive index ratio: \(\eta = n_1/n_2\)
Compute incident angle: \(\cos\theta_i = -\vec{d}_i \cdot \vec{n}\)
Compute discriminant: \(k = 1 - \eta^2(1 - \cos^2\theta_i)\)
Check for total internal reflection:
if k < 0: return None # Total internal reflection
Compute refracted direction:
cos_theta_t = sqrt(k) d_t = eta * d_i + (eta * cos_theta_i - cos_theta_t) * normal
Numerical Stability:
Use \(k = 1 - \eta^2 + \eta^2\cos^2\theta_i\) form to avoid cancellation
Threshold for TIR: \(k < -\epsilon\) with \(\epsilon = 10^{-10}\)
Intensity and Amplitude Tracking
Ray Intensity Update
At each surface interaction:
Reflection:
Refraction:
where \(R\) and \(T\) are Fresnel coefficients.
Propagation through absorbing medium:
where \(\alpha\) is the absorption coefficient and \(d\) is distance traveled.
Polarization State Update
For polarized rays, track s and p components separately:
Algorithm:
Compute local s and p reference frames at intersection:
# Plane of incidence plane_normal = cross(incident_dir, surface_normal) plane_normal /= norm(plane_normal) # s-direction (perpendicular to plane of incidence) s_dir = plane_normal # p-direction (in plane of incidence) p_dir = cross(s_dir, incident_dir)
Compute Fresnel coefficients: \(R_s, R_p, T_s, T_p\)
Update intensities:
I_s_refl = I_s_inc * R_s I_p_refl = I_p_inc * R_p I_s_trans = I_s_inc * T_s I_p_trans = I_p_inc * T_p
Rotate polarization frame for refracted ray (if needed)
Multi-Bounce Ray Tracing
Algorithm Overview
Input: Initial rays, list of surfaces, max bounces
Output: Final ray states, bounce history
Algorithm:
function trace_rays_multi_bounce(rays, surfaces, max_bounces):
active_rays = rays
all_bounces = []
for bounce in range(max_bounces):
if no active_rays:
break
# Find next intersection for each ray
intersections = []
for ray in active_rays:
nearest_t = infinity
nearest_surface = None
for surface in surfaces:
t, hit = surface.intersect(ray)
if hit and t < nearest_t:
nearest_t = t
nearest_surface = surface
intersections.append((nearest_t, nearest_surface))
# Update ray states
new_rays = []
for ray, (t, surface) in zip(active_rays, intersections):
if surface is None:
continue # Ray escaped
# Compute interaction point
hit_point = ray.position + t * ray.direction
normal = surface.normal_at(hit_point)
# Compute Fresnel coefficients
n1 = surface.material_front.n(ray.wavelength, hit_point)
n2 = surface.material_back.n(ray.wavelength, hit_point)
R, T = fresnel_coefficients(ray.direction, normal, n1, n2)
# Generate reflected ray
if ray.intensity * R > threshold:
refl_dir = reflect(ray.direction, normal)
refl_ray = Ray(
position=hit_point + epsilon * normal,
direction=refl_dir,
intensity=ray.intensity * R,
time=ray.time + t * n1 / c
)
new_rays.append(refl_ray)
# Generate refracted ray
if ray.intensity * T > threshold:
refr_dir = refract(ray.direction, normal, n1, n2)
if refr_dir is not None: # Check TIR
refr_ray = Ray(
position=hit_point - epsilon * normal,
direction=refr_dir,
intensity=ray.intensity * T,
time=ray.time + t * n1 / c
)
new_rays.append(refr_ray)
all_bounces.append(BounceRecord(ray, surface, hit_point))
active_rays = new_rays
return active_rays, all_bounces
Computational Complexity
For \(N\) rays, \(S\) surfaces, \(B\) bounces:
Naïve intersection: \(O(N \cdot S \cdot B)\)
With spatial acceleration: \(O(N \cdot \log S \cdot B)\)
GPU parallelization reduces to:
Per-ray cost: \(O(S \cdot B / N_{\text{threads}})\)
Accuracy and Error Analysis
Floating Point Precision
L-SURF uses double precision (float64) for:
Intersection distance calculations
Time-of-flight accumulation
Fresnel coefficient computation
Single precision (float32) for:
Ray positions and directions (after intersection)
Intensity values
Surface coordinates (GPU arrays)
Intersection Epsilon
To avoid self-intersection, rays are offset from surfaces by \(\epsilon\):
Reflected rays: offset in +normal direction
Refracted rays: offset in -normal direction
Typical value: \(\epsilon = 10^{-6}\) m
Convergence Criteria
Wave surface intersection:
Position error: \(|f(t)| < 10^{-6}\) m
Angle error: \(< 10^{-4}\) rad
Fresnel coefficients:
Energy conservation: \(|R + T - 1| < 10^{-10}\)
Intensity tracking:
Termination threshold: \(I < 10^{-6} I_0\)