Physics-Based Sensor Simulator

Implement camera and lidar sensor models from first principles, adding realistic noise, distortion, and weather artifacts that close the sim-to-real gap for autonomous driving perception.

Track: B -- Synthetic Data & Sensor Sim | Level: Intermediate | Total Time: ~20-25 hours

Overview

Clean, pixel-perfect synthetic images and noiseless lidar point clouds are tempting to work with but they are also dangerously unrealistic. A perception model trained exclusively on pristine renders will fail the moment it encounters the photon shot noise of a real CMOS sensor, the barrel distortion of a wide-angle lens, the range jitter of a spinning lidar, or the visibility degradation caused by rain and fog. Bridging this sensor realism gap is one of the highest-leverage improvements you can make to a synthetic data pipeline, and it requires surprisingly little compute -- just a solid understanding of the underlying physics.

In this project you will build camera and lidar sensor simulators from first principles. You will start with the camera imaging pipeline, tracing the path of light from scene radiance through a lens with geometric distortion, onto a noisy photodetector with shot noise, read noise, and dark current, through an analog-to-digital converter with quantization error, and finally into an 8-bit RGB image. You will then build a lidar model that accounts for time-of-flight measurement noise, beam divergence, intensity falloff with distance, angle-dependent ray dropout, and multipath returns. Finally you will add weather effects -- atmospheric scattering for fog, raindrop rendering for camera, backscatter and attenuation for lidar -- and integrate everything into a configurable sensor simulation pipeline.

The result is a pair of sensor models that can be dropped into any synthetic data pipeline to transform clean renders into data that looks like it came from real hardware. The physics parameters (read noise electrons, lens distortion coefficients, fog extinction coefficient) are all exposed as tunable knobs, so you can match any specific sensor by plugging in its datasheet values.

Every major autonomous driving simulation platform -- Applied Intuition's ADAS toolchain, NVIDIA DRIVE Sim, CARLA, and internal tools at Waymo and Cruise -- includes physics-based sensor models as a core component. The camera and lidar noise models you build here are simplified versions of what those platforms use in production. Understanding these models at a first-principles level will make you a stronger simulation engineer, and the code you write will be directly usable for data augmentation in any perception training pipeline.

Learning Objectives

After completing this project, you will be able to:

Trace the full camera imaging pipeline -- from scene radiance through lens optics, sensor photodetection, and analog-to-digital conversion to a final digital image, understanding where each noise source enters.
Implement photon-level noise models -- generate realistic Poisson shot noise, Gaussian read noise, dark current accumulation, and ADC quantization error, and understand how each depends on exposure time and sensor parameters.
Model geometric lens distortion -- implement the Brown-Conrady radial/tangential distortion model with k1-k3 and p1-p2 coefficients, and build both distort and undistort functions.
Build a physics-based lidar simulator -- model time-of-flight range noise, beam divergence footprints, Lambertian intensity with 1/r-squared falloff, angle-dependent dropout, and multi-return detection.
Simulate weather effects on both modalities -- implement Beer-Lambert atmospheric attenuation for fog, Mie backscatter for lidar, raindrop lens effects for camera, and rain-induced dropouts for lidar.
Parameterize models from real sensor datasheets -- map physical sensor specifications (quantum efficiency, read noise floor, distortion coefficients, beam divergence angle) to model parameters.
Visualize and validate sensor artifacts -- produce side-by-side comparisons of clean vs degraded data, noise power spectra, distortion field plots, and range-dependent error distributions.
Integrate sensor models into a configurable pipeline -- combine all effects into a single parameterized function call with a configuration dictionary that controls every physical parameter.

Prerequisites

Required

Python proficiency -- comfortable with NumPy array operations, function composition, and matplotlib visualization. You will work extensively with 2D arrays (images) and 3D arrays (point clouds).
Basic optics and physics -- understanding of how lenses form images, what photons are, and the inverse-square law for light intensity. You do not need a physics degree, but you should be comfortable with the idea that sensor measurements are noisy samples of a physical quantity.
Signal processing basics -- familiarity with noise as a random process, Gaussian and Poisson distributions, signal-to-noise ratio (SNR), and the concept of a noise floor.

Image processing experience -- prior use of OpenCV or scikit-image for tasks like filtering, warping, or color conversion.
Basic 3D geometry -- understanding of 3D coordinate systems, projection, and the pinhole camera model. Notebook 2 covers this from scratch but prior exposure helps.
Probability and statistics -- comfort with random variables, expected value, variance, and how to sample from common distributions using NumPy.

Deep Dive Reading

Before starting, read the companion deep dive for theoretical background:

Sensor Simulation for Autonomous Driving -- covers the full sensor simulation landscape including physics-based vs neural approaches, lidar models, camera models, radar, and the role of sensor sim in closing the sim-to-real gap.

Key Concepts

The Camera Imaging Pipeline

A real camera does not simply record the scene -- it transforms it through a chain of physical processes, each adding its own signature to the final image. Understanding this chain is the key to building a realistic sensor model.

Scene Radiance → Lens (distortion, vignetting) → Sensor (photon → electron)
    → Amplifier (gain, read noise) → ADC (quantization) → Digital Image

Detailed pipeline:
===============================================================
1. Scene radiance L(x,y,lambda)  [W/m^2/sr/nm]
   - What the world actually emits/reflects
   - Depends on lighting, materials, geometry

2. Lens optics
   - Geometric projection  (pinhole + distortion)
   - Vignetting             (cos^4 falloff at edges)
   - Chromatic aberration    (wavelength-dependent focus)
   - Diffraction             (Airy disk blur)

3. Sensor irradiance E(u,v)  [W/m^2]
   - Image formed on the sensor plane
   - E = L * (pi/4) * (D/f)^2 * cos^4(theta)

4. Photon detection
   - Exposure:    H = E * t_exp            [J/m^2]
   - Photons:     n_photon = H * A_pixel * QE / (h*c/lambda)
   - Signal:      n_electron ~ Poisson(n_photon)

5. Noise addition
   - Shot noise:     already in Poisson sampling
   - Dark current:   n_dark ~ Poisson(D * t_exp)
   - Read noise:     n_read ~ Normal(0, sigma_read)
   - Total:          n_total = n_electron + n_dark + n_read

6. Analog-to-Digital Conversion
   - Gain:   DN = n_total * K_gain
   - Clip:   DN = clip(DN, 0, 2^bit_depth - 1)
   - Quantize: DN = floor(DN)
===============================================================

The key insight is that noise is signal-dependent. Bright regions have high photon counts and therefore high absolute shot noise, but the relative noise (1/sqrt(N)) is small. Dark regions have few photons and are dominated by the read noise floor. This is why real cameras produce cleaner images in daylight than at dusk -- the signal-to-noise ratio is fundamentally limited by photon statistics.

Photon Shot Noise

Shot noise arises from the quantum nature of light. Photons arrive at the sensor as discrete particles following a Poisson process. If the expected number of photons at a pixel is N, the actual count in any single exposure is a random sample from Poisson(N).

The standard deviation of Poisson noise equals sqrt(N), so the signal-to-noise ratio is:

SNR_shot = N / sqrt(N) = sqrt(N)

This means doubling the light (or exposure time) only improves SNR by sqrt(2) = 1.41x. To get 10x better SNR you need 100x more photons. This fundamental relationship governs every camera system ever built.

In code, applying shot noise to a clean image is straightforward:

import numpy as np

def apply_shot_noise(clean_image: np.ndarray, photon_scale: float = 100.0) -> np.ndarray:
    """Apply Poisson shot noise to a clean image.

    Args:
        clean_image: float64 image in [0, 1] range
        photon_scale: mean photon count for a pixel at intensity 1.0
                      Higher = less noise. Typical range: 50-500.
    """
    # Convert normalized intensity to expected photon count
    expected_photons = clean_image * photon_scale

    # Sample actual photon counts from Poisson distribution
    actual_photons = np.random.poisson(expected_photons)

    # Convert back to normalized intensity
    noisy_image = actual_photons / photon_scale

    return noisy_image

Read Noise

Read noise comes from the electronics that convert the accumulated charge in each pixel into a voltage and then a digital number. It is independent of the signal level and adds a constant noise floor to every pixel. Read noise is well-modeled as zero-mean Gaussian with a standard deviation specified in electrons (e-).

def apply_read_noise(electron_count: np.ndarray, sigma_read: float = 5.0) -> np.ndarray:
    """Add Gaussian read noise.

    Args:
        electron_count: signal in electrons
        sigma_read: read noise standard deviation in electrons
                    Typical values: 2-10 e- for modern CMOS sensors
    """
    noise = np.random.normal(0, sigma_read, size=electron_count.shape)
    return electron_count + noise

Modern CMOS sensors achieve read noise as low as 1-2 electrons (e.g., Sony IMX sensors used in many autonomous driving cameras). Older CCD sensors had read noise of 5-15 electrons. The read noise floor determines the minimum detectable signal: a feature must produce more photoelectrons than the read noise to be visible.

Dark Current

Dark current is the trickle of electrons generated by thermal excitation in the silicon, even in complete darkness. It accumulates linearly with exposure time and doubles for every ~6-8 degrees Celsius increase in temperature. For automotive cameras operating in hot environments (engine bay, rooftop in Arizona sun), dark current can be significant.

def apply_dark_current(electron_count: np.ndarray,
                       dark_rate: float = 0.1,
                       exposure_time: float = 0.033) -> np.ndarray:
    """Add dark current noise.

    Args:
        electron_count: signal in electrons
        dark_rate: dark current rate in electrons/pixel/second
                   Typical: 0.01-1.0 e/px/s depending on temperature
        exposure_time: exposure time in seconds
    """
    expected_dark = dark_rate * exposure_time
    dark_electrons = np.random.poisson(expected_dark, size=electron_count.shape)
    return electron_count + dark_electrons

ADC Quantization

The analog-to-digital converter (ADC) maps the continuous voltage from the amplifier into discrete digital numbers. A 10-bit ADC maps to 0-1023, a 12-bit ADC to 0-4095, and a typical camera's 8-bit output to 0-255. Quantization introduces a rounding error uniformly distributed in [-0.5, +0.5] DN.

def apply_adc(electron_count: np.ndarray,
              gain: float = 0.5,
              bit_depth: int = 12,
              black_level: int = 64) -> np.ndarray:
    """Apply analog-to-digital conversion.

    Args:
        electron_count: signal in electrons
        gain: conversion gain in DN/electron. Higher gain = brighter image.
        bit_depth: ADC bit depth (typically 10, 12, or 14)
        black_level: black level offset added to avoid clipping negative noise
    """
    max_dn = 2**bit_depth - 1

    # Convert electrons to digital numbers
    dn = electron_count * gain + black_level

    # Clip and quantize
    dn = np.clip(dn, 0, max_dn)
    dn = np.floor(dn).astype(np.uint16)

    return dn

Lens Distortion Models

Real lenses do not produce perfect pinhole projections. The most common deviation is radial distortion, where straight lines in the world appear curved in the image. Wide-angle lenses used in autonomous driving exhibit significant barrel distortion (lines bow outward from center).

The standard Brown-Conrady model parameterizes distortion using radial coefficients (k1, k2, k3) and tangential coefficients (p1, p2):

Radial distortion:
  r^2 = x_n^2 + y_n^2
  x_radial = x_n * (1 + k1*r^2 + k2*r^4 + k3*r^6)
  y_radial = y_n * (1 + k1*r^2 + k2*r^4 + k3*r^6)

Tangential distortion:
  x_tangential = 2*p1*x_n*y_n + p2*(r^2 + 2*x_n^2)
  y_tangential = p1*(r^2 + 2*y_n^2) + 2*p2*x_n*y_n

Combined:
  x_distorted = x_radial + x_tangential
  y_distorted = y_radial + y_tangential

Here (x_n, y_n) are normalized image coordinates (after dividing by focal length and subtracting the principal point). Negative k1 produces barrel distortion; positive k1 produces pincushion distortion.

def distort_points(points_normalized: np.ndarray,
                   k1: float, k2: float, k3: float,
                   p1: float, p2: float) -> np.ndarray:
    """Apply Brown-Conrady distortion to normalized image points.

    Args:
        points_normalized: (N, 2) array of [x, y] normalized coordinates
        k1, k2, k3: radial distortion coefficients
        p1, p2: tangential distortion coefficients

    Returns:
        (N, 2) array of distorted normalized coordinates
    """
    x = points_normalized[:, 0]
    y = points_normalized[:, 1]

    r2 = x**2 + y**2
    r4 = r2**2
    r6 = r2**3

    # Radial distortion factor
    radial = 1 + k1 * r2 + k2 * r4 + k3 * r6

    # Tangential distortion
    dx_tang = 2 * p1 * x * y + p2 * (r2 + 2 * x**2)
    dy_tang = p1 * (r2 + 2 * y**2) + 2 * p2 * x * y

    x_dist = x * radial + dx_tang
    y_dist = y * radial + dy_tang

    return np.column_stack([x_dist, y_dist])

Typical distortion coefficients for automotive cameras:

Lens Type	k1	k2	k3	p1	p2
Narrow (50mm)	-0.05	0.01	0	0.001	-0.001
Standard (30mm)	-0.15	0.05	-0.01	0.002	-0.002
Wide-angle (12mm)	-0.35	0.15	-0.03	0.003	-0.003
Fisheye (6mm)	-0.50	0.30	-0.10	0.005	-0.005

Lidar Time-of-Flight Physics

A lidar measures distance by timing how long it takes a laser pulse to travel to a target and back. The range is:

R = c * t_round_trip / 2

where c is the speed of light (2.998 x 10^8 m/s). The precision of the time measurement determines the range noise. Modern lidar sensors achieve timing precision of 0.1-1 nanosecond, corresponding to range noise of 1.5-15 cm.

Range noise increases with distance for two reasons: (1) the return signal weakens as 1/R^2, reducing the signal-to-noise ratio of the timing circuit, and (2) beam divergence increases the footprint on the target, averaging over a larger area and adding geometric uncertainty.

def lidar_range_noise(ranges: np.ndarray,
                      base_sigma: float = 0.02,
                      distance_factor: float = 0.001) -> np.ndarray:
    """Generate distance-dependent range noise for lidar measurements.

    Args:
        ranges: (N,) array of true ranges in meters
        base_sigma: baseline range noise std in meters (at zero distance)
        distance_factor: additional noise per meter of distance
                        Total sigma = base_sigma + distance_factor * range

    Returns:
        (N,) array of noisy ranges
    """
    sigma = base_sigma + distance_factor * ranges
    noise = np.random.normal(0, sigma)
    noisy_ranges = ranges + noise
    return np.maximum(noisy_ranges, 0)  # range cannot be negative

Lidar Beam Divergence

A lidar beam is not infinitely thin -- it has a divergence angle (typically 0.1-0.5 milliradians for automotive lidars). At distance R, the beam illuminates a circular footprint of diameter:

d_footprint = 2 * R * tan(theta_div / 2) ≈ R * theta_div

For a 0.3 mrad beam at 100m, the footprint is 3 cm -- roughly the size of a tennis ball. At short range the footprint is smaller than most objects, but at long range it can span object boundaries, producing mixed returns and edge effects.

def beam_footprint_diameter(distance: float,
                            divergence_mrad: float = 0.3) -> float:
    """Calculate beam footprint diameter at a given distance.

    Args:
        distance: distance to target in meters
        divergence_mrad: full-angle beam divergence in milliradians

    Returns:
        footprint diameter in meters
    """
    divergence_rad = divergence_mrad * 1e-3
    return distance * divergence_rad

Lidar Intensity Model

The intensity of the return signal depends on the target's reflectance, the incidence angle, and the distance. For a Lambertian (diffuse) reflector:

I_return = I_0 * rho * cos(alpha) / R^2

where rho is the surface reflectance (0 to 1), alpha is the angle between the beam and the surface normal, and R is the distance. In practice lidar sensors report intensity as an 8-bit value (0-255) after automatic gain control.

def lidar_intensity(reflectance: np.ndarray,
                    incidence_angle: np.ndarray,
                    distance: np.ndarray,
                    reference_distance: float = 10.0) -> np.ndarray:
    """Calculate lidar return intensity using Lambertian model.

    Args:
        reflectance: (N,) surface reflectance values in [0, 1]
        incidence_angle: (N,) angle between beam and surface normal in radians
        distance: (N,) range to target in meters
        reference_distance: distance at which reflectance=1 gives intensity=255

    Returns:
        (N,) intensity values in [0, 255]
    """
    cos_alpha = np.clip(np.cos(incidence_angle), 0, 1)
    raw_intensity = reflectance * cos_alpha * (reference_distance / distance)**2
    intensity = np.clip(raw_intensity * 255, 0, 255).astype(np.uint8)
    return intensity

Lidar Ray Dropout

Not every lidar pulse produces a valid return. Rays can miss entirely (hitting the sky), be absorbed by dark surfaces, or return too weakly to detect. The dropout probability depends on:

Distance -- return signal drops as 1/R^2, so distant targets are more likely to drop out
Incidence angle -- highly oblique angles reflect very little energy back to the sensor
Surface reflectance -- dark or absorptive surfaces (matte black vehicles, wet asphalt) cause more dropouts

def compute_dropout_probability(distance: np.ndarray,
                                incidence_angle: np.ndarray,
                                reflectance: np.ndarray,
                                max_range: float = 120.0,
                                base_dropout: float = 0.01) -> np.ndarray:
    """Compute per-point dropout probability.

    Args:
        distance: (N,) range in meters
        incidence_angle: (N,) angle in radians
        reflectance: (N,) reflectance in [0, 1]
        max_range: maximum sensor range in meters
        base_dropout: baseline dropout rate for ideal conditions

    Returns:
        (N,) dropout probability in [0, 1]
    """
    # Distance factor: increases quadratically near max range
    dist_factor = (distance / max_range)**2

    # Angle factor: dropout increases at grazing angles
    cos_angle = np.clip(np.cos(incidence_angle), 0.01, 1.0)
    angle_factor = 1.0 - cos_angle

    # Reflectance factor: low reflectance increases dropout
    refl_factor = 1.0 - reflectance

    # Combined probability (clamped to [0, 1])
    prob = base_dropout + 0.3 * dist_factor + 0.3 * angle_factor + 0.2 * refl_factor
    return np.clip(prob, 0, 1)

Weather Effects on Cameras

Fog and haze reduce visibility by scattering light in the atmosphere. The physics is governed by the Beer-Lambert law and the Koschmieder atmospheric model:

I_observed(x) = I_scene(x) * T(d) + A * (1 - T(d))

where:
  T(d) = exp(-beta * d)    -- transmittance at distance d
  beta                      -- extinction coefficient (m^-1)
  A                         -- atmospheric light (airlight color, typically white-gray)
  d                         -- depth (distance from camera to scene point)

Objects far away have low transmittance T(d) and appear washed out, converging to the airlight color A. Objects close to the camera are barely affected. The extinction coefficient beta controls the density of the fog:

Condition	Visibility (m)	beta (m^-1)
Clear	> 10,000	< 0.0003
Light haze	2,000-5,000	0.0006-0.0015
Moderate fog	500-1,000	0.003-0.006
Dense fog	50-200	0.015-0.06
Very dense fog	< 50	> 0.06

Rain on camera creates two effects: (1) raindrops on the lens surface act as tiny lenses themselves, creating bright blurry spots, and (2) rain streaks in the air appear as bright diagonal lines due to motion blur during exposure.

Weather Effects on Lidar

Fog and rain affect lidar differently than cameras:

Attenuation -- the laser pulse loses energy as it passes through the medium, reducing the maximum effective range. The attenuation follows Beer-Lambert: the return signal is reduced by exp(-2 * beta * R) (factor of 2 because the pulse travels to the target and back).
Backscatter / false returns -- fog droplets and rain can scatter the laser pulse back to the sensor before it reaches the target, creating false return points at short range (typically 0-10m). These appear as a "noise cloud" around the sensor.
Random dropout increase -- reduced signal strength means more points fall below the detection threshold, increasing the dropout rate across the entire scan.

def fog_lidar_attenuation(ranges: np.ndarray,
                          intensities: np.ndarray,
                          beta: float = 0.01) -> tuple:
    """Apply fog attenuation to lidar returns.

    Args:
        ranges: (N,) range measurements in meters
        intensities: (N,) intensity values
        beta: extinction coefficient in m^-1

    Returns:
        Tuple of (attenuated_intensities, additional_dropout_mask)
    """
    # Two-way attenuation
    attenuation = np.exp(-2 * beta * ranges)

    # Attenuate intensities
    attenuated = intensities * attenuation

    # Points with very low intensity are dropped
    min_detectable = 5.0  # minimum detectable intensity
    dropout_mask = attenuated < min_detectable

    return attenuated, dropout_mask

Step-by-Step Implementation Guide

Step 1: Environment Setup

Create a virtual environment and install the required dependencies.

# Navigate to the project directory
cd projects/track-b-synthetic-data/physics-based-sensor-simulator

# Create and activate virtual environment
python -m venv .venv
source .venv/bin/activate

# Install dependencies
pip install -r requirements.txt

# Register Jupyter kernel
python -m ipykernel install --user --name sensor-sim --display-name "Physics Sensor Sim"

# Launch notebooks
jupyter notebook notebooks/

The core dependencies are:

Package	Purpose
numpy	Array operations, noise sampling
scipy	Signal processing, interpolation
matplotlib	Visualization
opencv-python	Image warping, calibration
Pillow	Image I/O
scikit-image	Image quality metrics, test images
open3d	3D point cloud visualization
tqdm	Progress bars

No GPU is required for this project. All computations run on CPU with NumPy.

Step 2: Camera Noise Model (Notebook 01)

Duration: ~60 minutes

In this notebook you implement the full camera imaging noise pipeline from photon arrival to digital number output. The key concept is that sensor noise is a physical process with well-understood statistics, not arbitrary random perturbation.

2.1 Setting Up the Clean Input

Start with a clean synthetic image. You can use any test image, but for autonomous driving realism, a scene with both bright sky and dark shadows works best because it exercises the full dynamic range of the noise model.

import numpy as np
from skimage import data, img_as_float64
import matplotlib.pyplot as plt

# Load a test image and convert to float64 [0, 1]
clean = img_as_float64(data.chelsea())  # Cat image, 300x451x3
print(f"Image shape: {clean.shape}, dtype: {clean.dtype}")
print(f"Intensity range: [{clean.min():.3f}, {clean.max():.3f}]")

# For driving scenes, you could also create a synthetic gradient image
# that simulates bright sky at the top and dark road at the bottom
H, W = 480, 640
gradient = np.zeros((H, W, 3), dtype=np.float64)
for i in range(H):
    intensity = 1.0 - (i / H) * 0.9  # bright at top, dark at bottom
    gradient[i, :] = intensity
gradient[:, :, 0] *= 0.6  # slight color tint
gradient[:, :, 2] *= 0.8

2.2 Photon Shot Noise

The most fundamental noise source. Every pixel independently samples from a Poisson distribution whose mean is the expected photon count.

def photon_to_electron(image_normalized: np.ndarray,
                       quantum_efficiency: float = 0.7,
                       well_capacity: int = 10000,
                       exposure_factor: float = 1.0) -> np.ndarray:
    """Convert normalized image to expected electron count.

    The well capacity represents the maximum number of electrons a pixel
    can hold before saturating (full well capacity). A pixel at intensity
    1.0 maps to well_capacity * exposure_factor electrons.
    """
    expected_electrons = image_normalized * well_capacity * exposure_factor * quantum_efficiency
    return expected_electrons


def apply_shot_noise(expected_electrons: np.ndarray) -> np.ndarray:
    """Sample actual electron counts from Poisson distribution."""
    # For very high counts, Poisson ≈ Gaussian (faster)
    # For low counts, use true Poisson sampling
    actual = np.where(
        expected_electrons > 1000,
        np.random.normal(expected_electrons, np.sqrt(np.maximum(expected_electrons, 1))),
        np.random.poisson(np.maximum(expected_electrons, 0).astype(np.int64))
    )
    return np.maximum(actual, 0).astype(np.float64)

The signal-to-noise ratio demonstrates the sqrt(N) relationship:

# Demonstrate SNR vs photon count
photon_counts = np.logspace(1, 5, 50)
measured_snr = []
for n in photon_counts:
    samples = np.random.poisson(n, size=10000)
    snr = np.mean(samples) / np.std(samples)
    measured_snr.append(snr)

plt.loglog(photon_counts, measured_snr, 'b.', label='Measured')
plt.loglog(photon_counts, np.sqrt(photon_counts), 'r-', label='sqrt(N) theory')
plt.xlabel('Expected photon count')
plt.ylabel('SNR')
plt.legend()
plt.title('Shot Noise SNR: sqrt(N) law')
plt.grid(True, alpha=0.3)
plt.show()

2.3 Read Noise and Dark Current

After adding photon shot noise, layer on the electronics noise:

def apply_read_noise(electrons: np.ndarray, sigma_read: float = 5.0) -> np.ndarray:
    """Add Gaussian read noise to electron count."""
    read_noise = np.random.normal(0, sigma_read, size=electrons.shape)
    return electrons + read_noise


def apply_dark_current(electrons: np.ndarray,
                       dark_rate: float = 0.5,
                       exposure_time: float = 0.033) -> np.ndarray:
    """Add Poisson dark current noise."""
    expected_dark = dark_rate * exposure_time
    dark_electrons = np.random.poisson(
        max(expected_dark, 0),
        size=electrons.shape
    ).astype(np.float64)
    return electrons + dark_electrons

2.4 Full Pipeline

Combine all stages into a single configurable function:

from dataclasses import dataclass

@dataclass
class CameraSensorParams:
    """Physical parameters of a camera sensor."""
    quantum_efficiency: float = 0.7    # Fraction of photons detected
    well_capacity: int = 10000         # Max electrons per pixel
    read_noise: float = 5.0            # Read noise std in electrons
    dark_current_rate: float = 0.5     # Dark current in e/pixel/second
    exposure_time: float = 0.033       # Exposure time in seconds (30fps)
    adc_bit_depth: int = 12            # ADC resolution
    adc_gain: float = 1.0              # DN per electron
    black_level: int = 64              # Black level offset

    @property
    def max_dn(self) -> int:
        return 2**self.adc_bit_depth - 1


def camera_noise_pipeline(clean_image: np.ndarray,
                          params: CameraSensorParams,
                          exposure_factor: float = 1.0) -> np.ndarray:
    """Apply the full camera noise pipeline to a clean image.

    Args:
        clean_image: float64 image in [0, 1] range, shape (H, W, 3)
        params: camera sensor parameters
        exposure_factor: multiplier on exposure (1.0 = normal, 0.5 = dark)

    Returns:
        uint16 image with realistic sensor noise
    """
    # Step 1: Convert to expected electron count
    expected_e = photon_to_electron(
        clean_image, params.quantum_efficiency,
        params.well_capacity, exposure_factor
    )

    # Step 2: Apply shot noise (Poisson)
    electrons = apply_shot_noise(expected_e)

    # Step 3: Add dark current
    electrons = apply_dark_current(
        electrons, params.dark_current_rate, params.exposure_time
    )

    # Step 4: Add read noise
    electrons = apply_read_noise(electrons, params.read_noise)

    # Step 5: ADC conversion
    dn = electrons * params.adc_gain + params.black_level
    dn = np.clip(dn, 0, params.max_dn)
    dn = np.floor(dn).astype(np.uint16)

    return dn

2.5 Comparison at Different Exposure Levels

The most revealing test is to render the same scene at multiple exposure levels:

params = CameraSensorParams()
exposure_levels = [0.1, 0.3, 1.0, 3.0]

fig, axes = plt.subplots(1, len(exposure_levels), figsize=(16, 4))
for ax, exp in zip(axes, exposure_levels):
    noisy = camera_noise_pipeline(clean, params, exposure_factor=exp)
    # Normalize to [0, 1] for display
    display = (noisy.astype(np.float64) - params.black_level) / (params.max_dn - params.black_level)
    display = np.clip(display, 0, 1)
    ax.imshow(display)
    ax.set_title(f'Exposure: {exp}x')
    ax.axis('off')
plt.suptitle('Camera Noise vs Exposure Level')
plt.tight_layout()
plt.show()

At low exposure (0.1x) the image should be visibly noisy, especially in dark regions. At high exposure (3.0x) it should be cleaner but with clipped highlights. This is exactly what happens with real cameras.

2.6 Exercise

Exercise 1: Implement a compute_snr_map function that takes a clean image and sensor parameters, runs the noise pipeline 100 times, and computes the per-pixel SNR (mean / std across trials). Visualize the SNR map as a heatmap. Verify that bright regions have higher SNR than dark regions.

Exercise 2: Compare noise characteristics for three sensor configurations: (a) a cheap dashcam (read noise = 15 e, well capacity = 5000, 8-bit ADC), (b) a standard automotive camera (read noise = 5 e, well capacity = 10000, 12-bit ADC), and (c) a premium sensor (read noise = 1.5 e, well capacity = 30000, 14-bit ADC). Generate the same scene with all three and plot the noise-vs-intensity curves.

Step 3: Lens Distortion (Notebook 02)

Duration: ~60 minutes

In this notebook you implement the Brown-Conrady lens distortion model, build both forward (distort) and inverse (undistort) functions, and visualize distortion fields for different lens types.

3.1 Pinhole Camera Model Review

The pinhole model projects a 3D point in camera coordinates to a 2D pixel coordinate:

def pinhole_project(points_3d: np.ndarray, K: np.ndarray) -> np.ndarray:
    """Project 3D points to 2D using pinhole model.

    Args:
        points_3d: (N, 3) points in camera frame [X, Y, Z]
        K: (3, 3) camera intrinsic matrix

    Returns:
        (N, 2) pixel coordinates [u, v]
    """
    # Perspective division
    x_norm = points_3d[:, 0] / points_3d[:, 2]
    y_norm = points_3d[:, 1] / points_3d[:, 2]

    # Apply intrinsics
    u = K[0, 0] * x_norm + K[0, 2]
    v = K[1, 1] * y_norm + K[1, 2]

    return np.column_stack([u, v])


# Example: construct intrinsic matrix for a 1920x1200 camera
fx, fy = 1000.0, 1000.0  # focal length in pixels
cx, cy = 960.0, 600.0    # principal point (image center)
K = np.array([[fx,  0, cx],
              [ 0, fy, cy],
              [ 0,  0,  1]])

3.2 Forward Distortion

The forward distortion function takes ideal (undistorted) normalized coordinates and produces distorted coordinates:

def apply_distortion(x_norm: np.ndarray, y_norm: np.ndarray,
                     k1: float, k2: float, k3: float,
                     p1: float, p2: float) -> tuple:
    """Apply Brown-Conrady distortion model.

    Args:
        x_norm, y_norm: normalized image coordinates (after removing K)
        k1, k2, k3: radial distortion coefficients
        p1, p2: tangential distortion coefficients

    Returns:
        (x_dist, y_dist) distorted normalized coordinates
    """
    r2 = x_norm**2 + y_norm**2
    r4 = r2**2
    r6 = r2**3

    # Radial distortion
    radial = 1 + k1 * r2 + k2 * r4 + k3 * r6
    x_rad = x_norm * radial
    y_rad = y_norm * radial

    # Tangential distortion
    x_tan = 2 * p1 * x_norm * y_norm + p2 * (r2 + 2 * x_norm**2)
    y_tan = p1 * (r2 + 2 * y_norm**2) + 2 * p2 * x_norm * y_norm

    x_dist = x_rad + x_tan
    y_dist = y_rad + y_tan

    return x_dist, y_dist

3.3 Distortion Field Visualization

A distortion field shows how each pixel is displaced by the lens distortion. This is the most intuitive way to understand what a set of distortion coefficients actually does:

def visualize_distortion_field(k1, k2, k3, p1, p2,
                               image_size=(640, 480),
                               K=None, grid_spacing=20):
    """Draw the distortion field as a vector plot."""
    W, H = image_size
    if K is None:
        K = np.array([[W/2, 0, W/2], [0, H/2, H/2], [0, 0, 1]])

    # Create grid of pixel coordinates
    u = np.arange(0, W, grid_spacing)
    v = np.arange(0, H, grid_spacing)
    uu, vv = np.meshgrid(u, v)

    # Convert to normalized coordinates
    x_norm = (uu - K[0, 2]) / K[0, 0]
    y_norm = (vv - K[1, 2]) / K[1, 1]

    # Apply distortion
    x_dist, y_dist = apply_distortion(x_norm, y_norm, k1, k2, k3, p1, p2)

    # Convert back to pixels
    u_dist = x_dist * K[0, 0] + K[0, 2]
    v_dist = y_dist * K[1, 1] + K[1, 2]

    # Displacement vectors
    du = u_dist - uu
    dv = v_dist - vv

    fig, ax = plt.subplots(figsize=(10, 7))
    magnitude = np.sqrt(du**2 + dv**2)
    ax.quiver(uu, vv, du, dv, magnitude, cmap='hot', scale=1, scale_units='xy')
    ax.set_xlim(0, W)
    ax.set_ylim(H, 0)  # flip y for image convention
    ax.set_aspect('equal')
    ax.set_title(f'Distortion Field (k1={k1}, k2={k2})')
    plt.colorbar(ax.collections[0], label='Displacement (pixels)')
    plt.tight_layout()
    plt.show()

3.4 Image Warping with Distortion

To apply distortion to an actual image, you need to compute a remap -- for each pixel in the output (distorted) image, find where it maps to in the input (undistorted) image. This is an inverse warp:

import cv2

def distort_image(image: np.ndarray, K: np.ndarray,
                  dist_coeffs: np.ndarray) -> np.ndarray:
    """Apply lens distortion to a clean image.

    Uses OpenCV's initUndistortRectifyMap in reverse:
    we want to go from undistorted -> distorted.
    """
    H, W = image.shape[:2]

    # Create pixel coordinate grids
    u = np.arange(W, dtype=np.float32)
    v = np.arange(H, dtype=np.float32)
    uu, vv = np.meshgrid(u, v)

    # Normalize coordinates
    x_norm = (uu - K[0, 2]) / K[0, 0]
    y_norm = (vv - K[1, 2]) / K[1, 1]

    k1, k2, p1, p2, k3 = dist_coeffs

    # Apply distortion
    x_dist, y_dist = apply_distortion(
        x_norm, y_norm, k1, k2, k3, p1, p2
    )

    # Convert back to pixel coordinates
    map_x = (x_dist * K[0, 0] + K[0, 2]).astype(np.float32)
    map_y = (y_dist * K[1, 1] + K[1, 2]).astype(np.float32)

    # Remap
    distorted = cv2.remap(image, map_x, map_y, cv2.INTER_LINEAR,
                          borderMode=cv2.BORDER_CONSTANT, borderValue=0)
    return distorted

3.5 Undistortion

Undistortion is the inverse operation -- given a distorted image from a real camera, produce a rectified pinhole image. The challenge is that the distortion function is not easily invertible analytically, so we use iterative refinement:

def undistort_points_iterative(x_dist: np.ndarray, y_dist: np.ndarray,
                                k1: float, k2: float, k3: float,
                                p1: float, p2: float,
                                n_iterations: int = 10) -> tuple:
    """Iteratively undistort normalized coordinates.

    Uses fixed-point iteration: start with distorted coords as initial
    guess for undistorted coords, then refine.
    """
    # Initialize with distorted coordinates
    x = x_dist.copy()
    y = y_dist.copy()

    for _ in range(n_iterations):
        # Apply distortion to current estimate
        x_d, y_d = apply_distortion(x, y, k1, k2, k3, p1, p2)

        # Update: correct by the error
        x = x + (x_dist - x_d)
        y = y + (y_dist - y_d)

    return x, y

OpenCV provides cv2.undistort() for image undistortion and cv2.undistortPoints() for point undistortion. In practice you should use these optimized implementations, but understanding the iterative algorithm helps debug calibration issues.

3.6 Calibration Basics

Camera calibration estimates both intrinsics (K) and distortion coefficients from images of a known pattern (checkerboard):

def calibrate_camera_from_checkerboard(images: list,
                                       pattern_size: tuple = (9, 6),
                                       square_size: float = 0.025):
    """Calibrate camera from checkerboard images.

    Args:
        images: list of grayscale images containing checkerboard
        pattern_size: (cols, rows) inner corners of checkerboard
        square_size: physical size of each square in meters

    Returns:
        K, dist_coeffs, rvecs, tvecs
    """
    # Prepare object points (3D coordinates of checkerboard corners)
    objp = np.zeros((pattern_size[0] * pattern_size[1], 3), np.float32)
    objp[:, :2] = np.mgrid[0:pattern_size[0],
                            0:pattern_size[1]].T.reshape(-1, 2) * square_size

    obj_points = []  # 3D points in world
    img_points = []  # 2D points in image

    for img in images:
        ret, corners = cv2.findChessboardCorners(img, pattern_size, None)
        if ret:
            # Refine corner positions to subpixel accuracy
            criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER,
                       30, 0.001)
            corners_refined = cv2.cornerSubPix(img, corners, (11, 11),
                                               (-1, -1), criteria)
            obj_points.append(objp)
            img_points.append(corners_refined)

    # Calibrate
    ret, K, dist_coeffs, rvecs, tvecs = cv2.calibrateCamera(
        obj_points, img_points, img.shape[::-1], None, None
    )

    print(f"Calibration RMS error: {ret:.4f} pixels")
    print(f"Intrinsic matrix K:\n{K}")
    print(f"Distortion coefficients: {dist_coeffs.ravel()}")

    return K, dist_coeffs, rvecs, tvecs

3.7 Exercise

Exercise 1: Generate a synthetic checkerboard image (no distortion), apply your distort_image function with k1=-0.3, k2=0.1, then use cv2.undistort to recover the original. Compute the mean pixel error between the original and recovered images. How does the error change with more extreme distortion values?

Exercise 2: Implement a vignetting model that reduces brightness near the edges of the image following the cos^4(theta) law, where theta is the angle from the optical axis. Combine vignetting with barrel distortion and camera noise to produce a triple-degraded image. Compare it to the original.

Step 4: Lidar Physics Model (Notebook 03)

Duration: ~60 minutes

In this notebook you build a physics-based lidar simulator. You will create a simple 3D scene using geometric primitives (planes, boxes, spheres), cast rays through a configurable scan pattern, compute intersections, and then apply realistic noise models for range, intensity, and dropout.

4.1 Lidar Scan Pattern

A spinning lidar produces a characteristic scan pattern. Each beam fires at a specific elevation angle, and the sensor rotates horizontally to sweep 360 degrees:

from dataclasses import dataclass
from typing import Optional
import numpy as np

@dataclass
class LidarConfig:
    """Configuration for a spinning lidar sensor."""
    n_beams: int = 32                    # Number of vertical beams
    fov_up: float = 10.0                 # Upper vertical FOV (degrees)
    fov_down: float = -30.0              # Lower vertical FOV (degrees)
    horizontal_resolution: float = 0.2   # Degrees between horizontal firings
    max_range: float = 120.0             # Maximum range (meters)
    min_range: float = 0.5              # Minimum range (meters)
    beam_divergence_mrad: float = 0.3    # Beam divergence (milliradians)

    @property
    def n_horizontal(self) -> int:
        return int(360.0 / self.horizontal_resolution)

    @property
    def elevation_angles(self) -> np.ndarray:
        """Elevation angles for each beam in radians."""
        return np.linspace(
            np.radians(self.fov_down),
            np.radians(self.fov_up),
            self.n_beams
        )

    @property
    def azimuth_angles(self) -> np.ndarray:
        """Azimuth angles for full rotation in radians."""
        return np.linspace(0, 2 * np.pi, self.n_horizontal, endpoint=False)


def generate_ray_directions(config: LidarConfig) -> np.ndarray:
    """Generate unit direction vectors for all lidar rays.

    Returns:
        (N, 3) array where N = n_beams * n_horizontal
    """
    elevations = config.elevation_angles
    azimuths = config.azimuth_angles

    # Create meshgrid of all angle combinations
    az, el = np.meshgrid(azimuths, elevations)
    az = az.ravel()
    el = el.ravel()

    # Convert to unit direction vectors
    x = np.cos(el) * np.cos(az)
    y = np.cos(el) * np.sin(az)
    z = np.sin(el)

    return np.column_stack([x, y, z])


config = LidarConfig(n_beams=32)
rays = generate_ray_directions(config)
print(f"Total rays per scan: {len(rays)}")
print(f"Horizontal steps: {config.n_horizontal}")
print(f"Points expected: {config.n_beams} x {config.n_horizontal} = {config.n_beams * config.n_horizontal}")

4.2 Simple Scene with Geometric Primitives

For testing the lidar model without a full rendering engine, we use analytic ray-primitive intersections:

def ray_plane_intersection(origins: np.ndarray, directions: np.ndarray,
                           plane_point: np.ndarray, plane_normal: np.ndarray,
                           reflectance: float = 0.5) -> dict:
    """Intersect rays with an infinite plane.

    Returns dict with 'distances', 'normals', 'reflectances', 'hit_mask'.
    """
    N = len(origins)
    denom = directions @ plane_normal
    hit_mask = np.abs(denom) > 1e-8

    distances = np.full(N, np.inf)
    normals = np.tile(plane_normal, (N, 1))
    reflectances = np.full(N, reflectance)

    valid = hit_mask
    t = ((plane_point - origins[valid]) @ plane_normal) / denom[valid]
    distances[valid] = t
    hit_mask[valid] &= (t > 0)

    return {
        'distances': distances,
        'normals': normals,
        'reflectances': reflectances,
        'hit_mask': hit_mask
    }


def ray_box_intersection(origins: np.ndarray, directions: np.ndarray,
                         box_min: np.ndarray, box_max: np.ndarray,
                         reflectance: float = 0.5) -> dict:
    """Intersect rays with an axis-aligned bounding box (AABB).

    Uses the slab method for efficient ray-AABB intersection.
    """
    N = len(origins)
    distances = np.full(N, np.inf)
    normals = np.zeros((N, 3))
    hit_mask = np.zeros(N, dtype=bool)

    # Slab method
    inv_dir = np.where(np.abs(directions) > 1e-10,
                       1.0 / directions,
                       np.sign(directions) * 1e10)

    t_min_vec = (box_min - origins) * inv_dir
    t_max_vec = (box_max - origins) * inv_dir

    t_enter = np.minimum(t_min_vec, t_max_vec)
    t_exit = np.maximum(t_min_vec, t_max_vec)

    t_near = np.max(t_enter, axis=1)
    t_far = np.min(t_exit, axis=1)

    valid = (t_near <= t_far) & (t_far > 0)
    t_hit = np.where(t_near > 0, t_near, t_far)

    hit_mask[valid] = True
    distances[valid] = t_hit[valid]

    # Compute normals at hit points
    for i in np.where(valid)[0]:
        hit_point = origins[i] + t_hit[i] * directions[i]
        # Find which face was hit
        for axis in range(3):
            if abs(hit_point[axis] - box_min[axis]) < 0.01:
                normals[i, axis] = -1
                break
            elif abs(hit_point[axis] - box_max[axis]) < 0.01:
                normals[i, axis] = 1
                break

    return {
        'distances': distances,
        'normals': normals,
        'reflectances': np.full(N, reflectance),
        'hit_mask': hit_mask
    }

4.3 Complete Scene and Raycasting

Build a simple driving scene and cast all lidar rays:

def build_test_scene():
    """Create a simple test scene with ground, walls, and boxes."""
    scene = []

    # Ground plane at z=0
    scene.append({
        'type': 'plane',
        'point': np.array([0, 0, 0]),
        'normal': np.array([0, 0, 1]),
        'reflectance': 0.3,  # asphalt
        'label': 'ground'
    })

    # Vehicle-like boxes
    vehicles = [
        {'center': [15, 2, 0.8], 'size': [4.5, 2.0, 1.6], 'refl': 0.6},
        {'center': [30, -1, 0.8], 'size': [4.5, 2.0, 1.6], 'refl': 0.4},
        {'center': [20, -8, 1.0], 'size': [6.5, 2.5, 2.0], 'refl': 0.5},
    ]
    for v in vehicles:
        c = np.array(v['center'])
        s = np.array(v['size']) / 2
        scene.append({
            'type': 'box',
            'min': c - s,
            'max': c + s,
            'reflectance': v['refl'],
            'label': 'vehicle'
        })

    # Building wall
    scene.append({
        'type': 'box',
        'min': np.array([5, 12, 0]),
        'max': np.array([40, 12.5, 5]),
        'reflectance': 0.7,
        'label': 'building'
    })

    return scene


def cast_rays(origins: np.ndarray, directions: np.ndarray,
              scene: list, max_range: float = 120.0) -> dict:
    """Cast rays against all scene primitives, return closest hits."""
    N = len(origins)
    closest_dist = np.full(N, np.inf)
    closest_normal = np.zeros((N, 3))
    closest_refl = np.zeros(N)
    any_hit = np.zeros(N, dtype=bool)

    for obj in scene:
        if obj['type'] == 'plane':
            result = ray_plane_intersection(
                origins, directions,
                obj['point'], obj['normal'], obj['reflectance']
            )
        elif obj['type'] == 'box':
            result = ray_box_intersection(
                origins, directions,
                obj['min'], obj['max'], obj['reflectance']
            )
        else:
            continue

        # Keep closest intersection
        closer = result['hit_mask'] & (result['distances'] < closest_dist)
        closer &= (result['distances'] > 0) & (result['distances'] < max_range)
        closest_dist[closer] = result['distances'][closer]
        closest_normal[closer] = result['normals'][closer]
        closest_refl[closer] = result['reflectances'][closer]
        any_hit[closer] = True

    return {
        'distances': closest_dist,
        'normals': closest_normal,
        'reflectances': closest_refl,
        'hit_mask': any_hit
    }

4.4 Applying Lidar Noise

Layer all noise sources onto the clean raycasting results:

@dataclass
class LidarNoiseParams:
    """Physical noise parameters for a lidar sensor."""
    range_noise_base: float = 0.02       # Base range noise std (meters)
    range_noise_distance_factor: float = 0.001  # Additional noise per meter
    intensity_noise_std: float = 5.0     # Intensity noise std (DN)
    base_dropout_rate: float = 0.02      # Baseline dropout probability
    max_range: float = 120.0             # Maximum detectable range
    min_detectable_intensity: float = 3.0  # Minimum intensity for detection


def apply_lidar_noise(clean_distances: np.ndarray,
                      clean_normals: np.ndarray,
                      clean_reflectances: np.ndarray,
                      hit_mask: np.ndarray,
                      ray_directions: np.ndarray,
                      params: LidarNoiseParams) -> dict:
    """Apply realistic noise to clean lidar measurements."""
    N = len(clean_distances)
    noisy_dist = clean_distances.copy()
    noisy_int = np.zeros(N)
    valid = hit_mask.copy()

    # --- Range noise (distance-dependent Gaussian) ---
    sigma = params.range_noise_base + params.range_noise_distance_factor * clean_distances
    range_noise = np.random.normal(0, sigma)
    noisy_dist[valid] += range_noise[valid]
    noisy_dist = np.clip(noisy_dist, 0, params.max_range)

    # --- Intensity model ---
    incidence_angles = np.zeros(N)
    for i in np.where(valid)[0]:
        cos_angle = abs(np.dot(-ray_directions[i], clean_normals[i]))
        incidence_angles[i] = np.arccos(np.clip(cos_angle, 0, 1))

    raw_intensity = (clean_reflectances * np.cos(incidence_angles) *
                     (10.0 / np.maximum(clean_distances, 0.1))**2 * 255)
    noisy_int = raw_intensity + np.random.normal(0, params.intensity_noise_std, N)
    noisy_int = np.clip(noisy_int, 0, 255)

    # --- Dropout model ---
    dist_factor = (clean_distances / params.max_range)**2
    cos_angles = np.cos(incidence_angles)
    angle_factor = 1.0 - np.clip(cos_angles, 0, 1)
    refl_factor = 1.0 - clean_reflectances

    dropout_prob = (params.base_dropout_rate +
                    0.3 * dist_factor +
                    0.3 * angle_factor +
                    0.2 * refl_factor)
    dropout_prob = np.clip(dropout_prob, 0, 1)

    dropout = np.random.random(N) < dropout_prob
    valid &= ~dropout

    # Also drop points with very low intensity
    valid &= noisy_int > params.min_detectable_intensity

    return {
        'distances': noisy_dist,
        'intensities': noisy_int,
        'valid_mask': valid
    }

4.5 Visualization

Convert polar measurements to Cartesian points and visualize:

def lidar_to_point_cloud(distances: np.ndarray,
                         ray_directions: np.ndarray,
                         sensor_origin: np.ndarray = np.zeros(3)) -> np.ndarray:
    """Convert lidar measurements to 3D point cloud."""
    points = sensor_origin + distances[:, np.newaxis] * ray_directions
    return points


def visualize_lidar_bev(points: np.ndarray, intensities: np.ndarray,
                        xlim=(-5, 50), ylim=(-20, 20),
                        title="Lidar BEV"):
    """Bird's-eye view of lidar point cloud."""
    fig, ax = plt.subplots(figsize=(12, 6))
    scatter = ax.scatter(
        points[:, 0], points[:, 1],
        c=intensities, cmap='viridis',
        s=0.5, alpha=0.8, vmin=0, vmax=255
    )
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    ax.set_aspect('equal')
    ax.set_xlabel('X (forward) [m]')
    ax.set_ylabel('Y (left) [m]')
    ax.set_title(title)
    ax.set_facecolor('#1a1a2e')
    plt.colorbar(scatter, label='Intensity')
    plt.tight_layout()
    plt.show()

4.6 Side-by-Side Comparison

Generate both clean and noisy scans and compare:

# Generate clean scan
config = LidarConfig(n_beams=32, horizontal_resolution=0.4)
ray_dirs = generate_ray_directions(config)
origins = np.tile(np.array([0, 0, 1.8]), (len(ray_dirs), 1))  # sensor at 1.8m height

scene = build_test_scene()
clean_result = cast_rays(origins, ray_dirs, scene, config.max_range)

# Clean point cloud
clean_mask = clean_result['hit_mask']
clean_pts = lidar_to_point_cloud(
    clean_result['distances'][clean_mask],
    ray_dirs[clean_mask],
    np.array([0, 0, 1.8])
)

# Noisy scan
noise_params = LidarNoiseParams()
noisy_result = apply_lidar_noise(
    clean_result['distances'],
    clean_result['normals'],
    clean_result['reflectances'],
    clean_result['hit_mask'],
    ray_dirs,
    noise_params
)

noisy_mask = noisy_result['valid_mask']
noisy_pts = lidar_to_point_cloud(
    noisy_result['distances'][noisy_mask],
    ray_dirs[noisy_mask],
    np.array([0, 0, 1.8])
)

print(f"Clean points: {clean_mask.sum()}")
print(f"Noisy points: {noisy_mask.sum()}")
print(f"Dropout rate: {1 - noisy_mask.sum()/clean_mask.sum():.2%}")

4.7 Exercise

Exercise 1: Add a multi-return model to the lidar. When a beam hits an object edge where the footprint spans two surfaces at different depths, the sensor can report two returns (first and last). Implement a function that takes the beam footprint size at the target distance and, if the footprint exceeds the object width at the hit point, generates a second return at a farther distance.

Exercise 2: Implement a range-dependent noise histogram analysis. Bin all points by distance (0-20m, 20-40m, 40-60m, 60-80m, 80-120m) and for each bin, compute the mean and standard deviation of range error. Plot these against distance and verify they follow the linear model sigma = base + factor * distance.

Step 5: Weather Effects and Integration (Notebook 04)

Duration: ~60 minutes

In this notebook you implement weather effects on both camera and lidar, then integrate all sensor models into a unified pipeline with a configuration system.

5.1 Fog on Camera

The Koschmieder model adds depth-dependent haze to camera images:

def apply_fog_camera(image: np.ndarray, depth_map: np.ndarray,
                     beta: float = 0.01,
                     airlight: np.ndarray = np.array([0.8, 0.8, 0.85])) -> np.ndarray:
    """Apply fog to a camera image using the Koschmieder atmospheric model.

    Args:
        image: (H, W, 3) float64 image in [0, 1]
        depth_map: (H, W) depth in meters
        beta: extinction coefficient (higher = denser fog)
        airlight: (3,) color of atmospheric light

    Returns:
        (H, W, 3) foggy image
    """
    # Transmittance: fraction of scene radiance that reaches the camera
    transmittance = np.exp(-beta * depth_map)
    transmittance = transmittance[:, :, np.newaxis]  # broadcast to 3 channels

    # Koschmieder equation
    foggy = image * transmittance + airlight * (1 - transmittance)

    return np.clip(foggy, 0, 1)


# Visibility distance (where transmittance drops to 5%)
def beta_from_visibility(visibility_m: float) -> float:
    """Convert meteorological visibility to extinction coefficient."""
    return -np.log(0.05) / visibility_m  # ≈ 3.0 / visibility

5.2 Rain on Camera

Rain produces two effects: raindrops on the lens and streaks in the air.

def add_rain_streaks(image: np.ndarray,
                     n_streaks: int = 200,
                     streak_length: tuple = (15, 40),
                     streak_width: float = 1,
                     intensity: float = 0.3,
                     angle_range: tuple = (70, 110)) -> np.ndarray:
    """Add rain streak artifacts to a camera image.

    Args:
        image: (H, W, 3) float64 image in [0, 1]
        n_streaks: number of rain streaks to render
        streak_length: (min, max) length of streaks in pixels
        streak_width: width of each streak in pixels
        intensity: brightness of streaks (0 = invisible, 1 = white)
        angle_range: (min, max) angle of streaks in degrees (90 = vertical)
    """
    H, W = image.shape[:2]
    rain_layer = np.zeros((H, W), dtype=np.float64)

    for _ in range(n_streaks):
        # Random start position
        x0 = np.random.randint(0, W)
        y0 = np.random.randint(0, H)

        # Random length and angle
        length = np.random.randint(*streak_length)
        angle = np.radians(np.random.uniform(*angle_range))

        # End position
        x1 = int(x0 + length * np.cos(angle))
        y1 = int(y0 + length * np.sin(angle))

        # Draw line with anti-aliasing using Bresenham
        cv2.line(rain_layer, (x0, y0), (x1, y1),
                 intensity, thickness=int(streak_width),
                 lineType=cv2.LINE_AA)

    # Apply slight Gaussian blur to soften streaks
    rain_layer = cv2.GaussianBlur(rain_layer, (3, 3), 0.5)

    # Add rain layer to image
    rainy = image + rain_layer[:, :, np.newaxis]
    return np.clip(rainy, 0, 1)


def add_lens_raindrops(image: np.ndarray,
                       n_drops: int = 30,
                       radius_range: tuple = (5, 20),
                       blur_strength: float = 5.0) -> np.ndarray:
    """Simulate raindrops on the camera lens.

    Each raindrop acts as a small lens, creating a blurred/distorted patch.
    """
    H, W = image.shape[:2]
    result = image.copy()

    for _ in range(n_drops):
        cx = np.random.randint(0, W)
        cy = np.random.randint(0, H)
        radius = np.random.randint(*radius_range)

        # Create circular mask
        y, x = np.ogrid[-cy:H-cy, -cx:W-cx]
        mask = x**2 + y**2 <= radius**2

        # Extract patch and apply heavy blur (lens effect)
        if mask.any():
            blurred = cv2.GaussianBlur(
                image, (0, 0), sigmaX=blur_strength
            )
            # Blend blurred patch with slight brightness increase
            result[mask] = blurred[mask] * 0.7 + 0.3  # brightened refraction

    return np.clip(result, 0, 1)

5.3 Fog on Lidar

Fog affects lidar through two mechanisms: attenuation (weakened returns) and backscatter (false close-range returns):

def apply_fog_lidar(distances: np.ndarray,
                    intensities: np.ndarray,
                    valid_mask: np.ndarray,
                    beta: float = 0.01,
                    backscatter_rate: float = 0.05,
                    backscatter_range: tuple = (0.5, 8.0)) -> dict:
    """Apply fog effects to lidar measurements.

    Args:
        distances: (N,) range measurements
        intensities: (N,) intensity values
        valid_mask: (N,) boolean mask of valid returns
        beta: extinction coefficient (same as camera fog)
        backscatter_rate: fraction of rays that get a false backscatter return
        backscatter_range: (min, max) distance range for backscatter returns

    Returns:
        dict with modified distances, intensities, valid_mask, and
        additional backscatter points
    """
    N = len(distances)

    # --- Attenuation ---
    # Two-way path through fog
    attenuation = np.exp(-2 * beta * distances)
    attenuated_intensity = intensities * attenuation

    # Drop returns that fall below detection threshold
    min_detectable = 3.0
    fog_dropout = valid_mask & (attenuated_intensity < min_detectable)
    new_valid = valid_mask & ~fog_dropout

    # --- Backscatter (false returns from fog particles) ---
    n_backscatter = int(N * backscatter_rate)
    bs_indices = np.random.choice(N, size=n_backscatter, replace=False)
    bs_distances = np.random.uniform(*backscatter_range, size=n_backscatter)
    bs_intensities = np.random.uniform(3, 30, size=n_backscatter)

    return {
        'distances': distances,
        'intensities': attenuated_intensity,
        'valid_mask': new_valid,
        'backscatter_distances': bs_distances,
        'backscatter_intensities': bs_intensities,
        'backscatter_indices': bs_indices,
        'fog_dropout_count': fog_dropout.sum(),
    }

5.4 Rain on Lidar

Rain causes random dropouts and intensity reduction:

def apply_rain_lidar(distances: np.ndarray,
                     intensities: np.ndarray,
                     valid_mask: np.ndarray,
                     rain_rate: float = 10.0) -> dict:
    """Apply rain effects to lidar measurements.

    Args:
        distances: (N,) range measurements
        intensities: (N,) intensity values
        valid_mask: (N,) boolean mask
        rain_rate: rainfall rate in mm/hr. Light=2, Moderate=10, Heavy=25

    Returns:
        dict with modified measurements
    """
    N = len(distances)

    # Rain causes a general extinction (less severe than fog)
    # Empirical model: beta_rain ≈ 0.01 * rain_rate^0.6
    beta_rain = 0.01 * rain_rate**0.6

    # Attenuation
    attenuation = np.exp(-2 * beta_rain * distances)
    attenuated_intensity = intensities * attenuation

    # Random dropouts (raindrops intercepting individual beams)
    # Probability increases with rain rate
    dropout_prob = 0.005 * rain_rate**0.5
    rain_dropout = np.random.random(N) < dropout_prob

    new_valid = valid_mask & ~rain_dropout
    new_valid &= attenuated_intensity > 3.0  # detection threshold

    return {
        'distances': distances,
        'intensities': attenuated_intensity,
        'valid_mask': new_valid,
        'rain_dropout_count': rain_dropout.sum()
    }

5.5 Combined Weather Pipeline

Integrate all weather effects into a single configurable pipeline:

from dataclasses import dataclass, field
from enum import Enum

class WeatherCondition(Enum):
    CLEAR = "clear"
    LIGHT_FOG = "light_fog"
    DENSE_FOG = "dense_fog"
    LIGHT_RAIN = "light_rain"
    HEAVY_RAIN = "heavy_rain"
    FOG_AND_RAIN = "fog_and_rain"


WEATHER_PRESETS = {
    WeatherCondition.CLEAR: {
        'fog_beta': 0.0,
        'rain_rate': 0.0,
        'n_rain_streaks': 0,
        'n_lens_drops': 0,
    },
    WeatherCondition.LIGHT_FOG: {
        'fog_beta': 0.005,
        'rain_rate': 0.0,
        'n_rain_streaks': 0,
        'n_lens_drops': 0,
    },
    WeatherCondition.DENSE_FOG: {
        'fog_beta': 0.03,
        'rain_rate': 0.0,
        'n_rain_streaks': 0,
        'n_lens_drops': 0,
    },
    WeatherCondition.LIGHT_RAIN: {
        'fog_beta': 0.001,
        'rain_rate': 5.0,
        'n_rain_streaks': 100,
        'n_lens_drops': 10,
    },
    WeatherCondition.HEAVY_RAIN: {
        'fog_beta': 0.003,
        'rain_rate': 25.0,
        'n_rain_streaks': 400,
        'n_lens_drops': 40,
    },
    WeatherCondition.FOG_AND_RAIN: {
        'fog_beta': 0.015,
        'rain_rate': 10.0,
        'n_rain_streaks': 200,
        'n_lens_drops': 20,
    },
}


def apply_weather_camera(image: np.ndarray, depth_map: np.ndarray,
                         weather: WeatherCondition) -> np.ndarray:
    """Apply weather effects to a camera image."""
    params = WEATHER_PRESETS[weather]
    result = image.copy()

    # Apply fog
    if params['fog_beta'] > 0:
        result = apply_fog_camera(result, depth_map, beta=params['fog_beta'])

    # Apply rain streaks
    if params['n_rain_streaks'] > 0:
        result = add_rain_streaks(result, n_streaks=params['n_rain_streaks'])

    # Apply lens raindrops
    if params['n_lens_drops'] > 0:
        result = add_lens_raindrops(result, n_drops=params['n_lens_drops'])

    return result


def apply_weather_lidar(distances: np.ndarray,
                        intensities: np.ndarray,
                        valid_mask: np.ndarray,
                        weather: WeatherCondition) -> dict:
    """Apply weather effects to lidar measurements."""
    params = WEATHER_PRESETS[weather]

    result = {
        'distances': distances.copy(),
        'intensities': intensities.copy(),
        'valid_mask': valid_mask.copy()
    }

    # Apply fog effects
    if params['fog_beta'] > 0:
        fog_result = apply_fog_lidar(
            result['distances'], result['intensities'],
            result['valid_mask'], beta=params['fog_beta']
        )
        result.update(fog_result)

    # Apply rain effects
    if params['rain_rate'] > 0:
        rain_result = apply_rain_lidar(
            result['distances'], result['intensities'],
            result['valid_mask'], rain_rate=params['rain_rate']
        )
        result['intensities'] = rain_result['intensities']
        result['valid_mask'] &= rain_result['valid_mask']

    return result

5.6 Full Sensor Pipeline

Bring camera noise, distortion, and weather together into a single call:

@dataclass
class SensorPipelineConfig:
    """Complete configuration for the sensor simulation pipeline."""
    # Camera noise
    camera: CameraSensorParams = field(default_factory=CameraSensorParams)
    exposure_factor: float = 1.0

    # Lens distortion
    k1: float = -0.15
    k2: float = 0.05
    k3: float = 0.0
    p1: float = 0.001
    p2: float = -0.001

    # Lidar
    lidar: LidarConfig = field(default_factory=LidarConfig)
    lidar_noise: LidarNoiseParams = field(default_factory=LidarNoiseParams)

    # Weather
    weather: WeatherCondition = WeatherCondition.CLEAR


def run_camera_pipeline(clean_image: np.ndarray,
                        depth_map: np.ndarray,
                        K: np.ndarray,
                        config: SensorPipelineConfig) -> np.ndarray:
    """Run the full camera simulation pipeline.

    Order: weather -> noise -> distortion
    """
    # 1. Weather effects (operate on scene radiance)
    image = apply_weather_camera(clean_image, depth_map, config.weather)

    # 2. Sensor noise
    noisy = camera_noise_pipeline(image, config.camera, config.exposure_factor)

    # 3. Normalize back to float for distortion
    noisy_float = (noisy.astype(np.float64) - config.camera.black_level)
    noisy_float /= (config.camera.max_dn - config.camera.black_level)
    noisy_float = np.clip(noisy_float, 0, 1)

    # 4. Lens distortion
    dist_coeffs = np.array([config.k1, config.k2, config.p1, config.p2, config.k3])
    distorted = distort_image((noisy_float * 255).astype(np.uint8),
                              K, dist_coeffs)

    return distorted

5.7 Before/After Gallery

The most compelling way to demonstrate the sensor models is a side-by-side gallery across weather conditions:

conditions = [
    WeatherCondition.CLEAR,
    WeatherCondition.LIGHT_FOG,
    WeatherCondition.DENSE_FOG,
    WeatherCondition.HEAVY_RAIN,
]

fig, axes = plt.subplots(2, len(conditions), figsize=(20, 8))
for i, weather in enumerate(conditions):
    config = SensorPipelineConfig(weather=weather)

    # Camera
    cam_result = run_camera_pipeline(clean_image, depth_map, K, config)
    axes[0, i].imshow(cam_result)
    axes[0, i].set_title(f'Camera: {weather.value}')
    axes[0, i].axis('off')

    # Lidar
    lidar_result = apply_weather_lidar(
        clean_distances, clean_intensities, clean_valid, weather
    )
    valid = lidar_result['valid_mask']
    pts = lidar_to_point_cloud(
        lidar_result['distances'][valid],
        ray_dirs[valid]
    )
    axes[1, i].scatter(pts[:, 0], pts[:, 1], c='cyan', s=0.3, alpha=0.5)
    axes[1, i].set_title(f'Lidar: {weather.value} ({valid.sum()} pts)')
    axes[1, i].set_facecolor('#1a1a2e')
    axes[1, i].set_xlim(-5, 50)
    axes[1, i].set_ylim(-20, 20)
    axes[1, i].set_aspect('equal')

plt.suptitle('Sensor Degradation Across Weather Conditions', fontsize=14)
plt.tight_layout()
plt.show()

5.8 Exercise

Exercise 1: Implement a snow model for both camera and lidar. Snow on camera should include snowflakes (bright dots scattered across the image) and accumulated snow on surfaces (brighten the upper portions of objects). Snow on lidar should produce backscatter similar to fog but at higher altitudes and with lower intensity. Define a SNOW weather preset with appropriate parameters.

Exercise 2: Build a sensor degradation sweep that generates the same scene with fog beta ranging from 0 to 0.05 in 20 steps. For each step, record: (a) the number of surviving lidar points, (b) the mean camera image SSIM compared to the clean image, and (c) the mean lidar range error. Plot all three metrics against beta on a single figure with three y-axes. At what fog density does the camera become more degraded than the lidar, and vice versa?

Notebook Overview

#	Notebook	Focus	Duration
1	`01_camera_noise_model.ipynb`	Shot noise, read noise, dark current, ADC pipeline	60 min
2	`02_lens_distortion.ipynb`	Brown-Conrady model, distort/undistort, calibration	60 min
3	`03_lidar_physics.ipynb`	Range noise, beam divergence, intensity, dropout, raycasting	60 min
4	`04_weather_effects.ipynb`	Fog, rain on camera and lidar, integrated pipeline	60 min

What You Will Build

By the end of this project you will have:

A camera noise simulator that converts clean renders into physically realistic images with signal-dependent shot noise, read noise, dark current, and ADC quantization -- parameterized by real sensor datasheet values.
A lens distortion engine that applies and inverts the Brown-Conrady model, letting you simulate any automotive lens from narrow telephoto to extreme wide-angle.
A physics-based lidar simulator that generates realistic point clouds from 3D scenes, complete with distance-dependent range noise, Lambertian intensity modeling, beam divergence effects, and angle/distance/reflectance-dependent dropout.
A weather effects pipeline covering fog (Beer-Lambert attenuation), rain (streaks, lens drops), and their distinct impacts on both camera and lidar modalities.
A unified sensor configuration system that bundles all parameters into a single config object, enabling one-line sensor degradation for any synthetic scene.

The code is pure NumPy/SciPy with no heavy dependencies, making it easy to integrate into any synthetic data pipeline or use as a standalone data augmentation tool for perception model training.

Tips and Troubleshooting

Poisson sampling is slow for large arrays: Use the Gaussian approximation (mean=N, std=sqrt(N)) when the expected count exceeds ~1000. The error is negligible.
Distortion remap uses inverse mapping: To produce a distorted image, you compute where each output pixel maps in the input (not the other way around). Getting this backwards is the most common bug.
Lidar noise parameters are sensor-specific: A Velodyne VLP-16 and a Hesai AT128 have very different noise characteristics. Always check the datasheet.
Fog beta vs visibility: Remember that visibility = 3/beta (approximately). A beta of 0.01 corresponds to 300m visibility (moderate fog), not "slightly hazy."
Combine effects in the right order: Weather affects the scene radiance (apply first), then sensor noise acts on the received signal, then lens distortion warps the geometry. Reversing the order produces unrealistic results.

Physics-Based Sensor Simulator

Overview

Learning Objectives

Prerequisites

Required

Recommended

Deep Dive Reading

Key Concepts

The Camera Imaging Pipeline

Photon Shot Noise

Read Noise

Dark Current

ADC Quantization

Lens Distortion Models

Lidar Time-of-Flight Physics

Lidar Beam Divergence

Lidar Intensity Model

Lidar Ray Dropout

Weather Effects on Cameras

Weather Effects on Lidar

Step-by-Step Implementation Guide

Step 1: Environment Setup

Step 2: Camera Noise Model (Notebook 01)

2.1 Setting Up the Clean Input

2.2 Photon Shot Noise

2.3 Read Noise and Dark Current

2.4 Full Pipeline

2.5 Comparison at Different Exposure Levels

2.6 Exercise

Step 3: Lens Distortion (Notebook 02)

3.1 Pinhole Camera Model Review

3.2 Forward Distortion

3.3 Distortion Field Visualization

3.4 Image Warping with Distortion

3.5 Undistortion

3.6 Calibration Basics

3.7 Exercise

Step 4: Lidar Physics Model (Notebook 03)

4.1 Lidar Scan Pattern

4.2 Simple Scene with Geometric Primitives

4.3 Complete Scene and Raycasting

4.4 Applying Lidar Noise

4.5 Visualization

4.6 Side-by-Side Comparison

4.7 Exercise

Step 5: Weather Effects and Integration (Notebook 04)

5.1 Fog on Camera

5.2 Rain on Camera

5.3 Fog on Lidar

5.4 Rain on Lidar

5.5 Combined Weather Pipeline

5.6 Full Sensor Pipeline

5.7 Before/After Gallery

5.8 Exercise

Notebook Overview

What You Will Build

Further Reading

Tips and Troubleshooting