Cone Illumination
In a real camera system, light reaching each pixel comes from the full area of the lens exit pupil, not from a single direction. The ConeIllumination class models this by decomposing the illumination cone into weighted planewaves and integrating the results. This guide is the practical, code-first walkthrough.
Theory background
For the underlying physics — CRA, the F-number-to-cone-angle relationship, pupil weighting, and the angle-induced shift of thin-film filter response — see Pixel Optical Effects and Thin Film Optics.
Cone parameters
Interactive Cone Illumination Viewer
Visualize how Chief Ray Angle (CRA) and cone half-angle affect pixel illumination. The cone is rendered as multiple parallel beam bundles (one per direction); each bundle is refracted by the microlens and the high-index stack via Snell's law, converging to a direction-specific focal point.
A cone is described by three controls: CRA (chief ray angle to the pixel), F-number (which sets the half-cone via
Why cone illumination matters for color filters
A thin-film color filter (or any Fabry-Perot-like cavity in the stack) shifts its transmission peak to shorter wavelengths under oblique incidence -- approximately
Interactive: Fabry-Perot cone-integration
The simulator below reproduces the central result of Goossens et al. (2018) inline. Sweep CRA and F-number to watch the integrated transmittance peak blue-shift and broaden away from the plane-wave curve. The wavelength window and sample count are the same knobs you control in ConeIllumination(cra_deg=..., f_number=..., n_points=...).
Fabry-Perot Cone Illumination Simulator
Reproduces the central result of Goossens et al. (2018), Appl. Opt. 57(26):7539. A single-cavity Fabry-Perot filter on an image sensor is illuminated by a focused cone of light; the simulator integrates the Airy transmittance over the cone defined by the chief-ray angle and the F-number, and compares against the plane-wave result.
Cone samples (top view, sensor frame)
Model details
The Airy transmittance of a single-cavity Fabry-Perot filter is
T(λ, θ) = 1 / [1 + F · sin²(δ(λ, θ) / 2)], with δ = (4π neff d cosθint) / λ, F = 4R / (1 − R)², neff sinθint = sinθ.
The cavity thickness d is fixed by the design condition (first-order peak at λ₀ at normal incidence): d = λ₀ / (2 n_eff). For each cone sample, the local ray direction is rotated by the chief-ray angle around the y-axis to obtain the actual incidence angle θ on the filter; the cone is sampled on a Fibonacci spiral over the half-cone θ_h = arcsin(1 / 2F) and weighted by cosθ (aplanatic pupil).
The cone-integrated transmittance is then
Tcone(λ) = ∑i wi T(λ, θi) / ∑i wi.
This reproduces the two effects emphasized in Goossens et al. (2018): a blue-shift of the centroid (because the average cosθ_int over the cone is less than the chief-ray value) and a broadening of the transmission peak (because each sample peaks at a different wavelength).
Reference: T. Goossens et al., Finite aperture correction for spectral cameras with integrated thin-film Fabry-Perot filters, Appl. Opt. 57(26):7539 (2018). DOI: 10.1364/AO.57.007539.
Top view: footprint on the pixel array
Cone Illumination – Top View
Bird's eye view of cone illumination on a 2×2 Bayer pixel array. Adjust CRA, f-number, and sampling to see how the illumination footprint covers the pixels.
Recommended default: golden-angle, near-uniform coverage without ring artifacts.
The side view above shows the cone geometry in cross-section. The top view provides a complementary perspective: looking down at the pixel array from above, you can see how the illumination cone projects onto the 2x2 Bayer pattern.
Key observations from the top view:
- Footprint diameter: The cone footprint on the focal plane has diameter
, where is the effective propagation height used for the cone spread. A lower F-number produces a wider footprint. - CRA shift: A nonzero CRA shifts the footprint center away from the pixel center. In a real BSI stack this is not the raw air-path value
; refraction in the color-filter, BARL, and silicon layers bends the chief ray toward normal, so the effective shift is smaller. - Sampling coverage: The interactive viewer above compares fibonacci, equal-area rings, Halton and Hammersley low-discrepancy sets, Gauss-Legendre, and legacy polar-grid sampling.
gridis useful as a baseline, but it is usually not the best production choice. - Lens area: The footprint area
where determines how much of the neighboring pixel receives light from the cone, which directly affects crosstalk.
Creating a ConeIllumination instance
from compass.sources.cone_illumination import ConeIllumination
cone = ConeIllumination(
cra_deg=15.0, # Chief Ray Angle in degrees
f_number=2.0, # F-number of the lens
n_points=37, # Number of angular sample points
sampling="fibonacci", # "fibonacci", "rings", "halton", "hammersley", "gauss", or "grid"
weighting="cosine", # "uniform", "cosine", "cos4", or "gaussian"
)
print(f"Half-cone angle: {cone.half_cone_rad * 180 / 3.14159:.1f} degrees")Sampling methods
The get_sampling_points() method returns a list of (theta_deg, phi_deg, weight) tuples. Each tuple represents a planewave direction and its associated integration weight.
Fibonacci sampling
Fibonacci (golden-angle spiral) sampling distributes points quasi-uniformly over the cone area. It provides good coverage with relatively few points.
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=37, sampling="fibonacci"
)
points = cone.get_sampling_points()
print(f"Number of sample points: {len(points)}")
for i, (theta, phi, w) in enumerate(points[:5]):
print(f" Point {i}: theta={theta:.2f} deg, phi={phi:.2f} deg, weight={w:.4f}")Fibonacci sampling is recommended for most cases. Use 19-37 points for quick estimates and 61-91 points for production results.
Equal-area rings
Concentric ring sampling divides the cone cap into equal-area radial bands and allocates more azimuth samples to larger outer rings. It is deterministic and easy to inspect visually, which makes it a good alternative when a regular polar grid looks too structured.
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=37, sampling="rings"
)
points = cone.get_sampling_points()
print(f"Ring sampling: {len(points)} points")Halton sampling
Halton sampling uses a low-discrepancy sequence over the cone cap. It avoids obvious ring or spoke symmetry and is useful for convergence checks where you want a deterministic quasi-random pattern.
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=37, sampling="halton"
)Halton is incremental: adding more samples preserves the earlier sequence. This makes it convenient for progressive convergence checks.
Hammersley sampling
Hammersley sampling is also a low-discrepancy point set, but it assumes the total sample count is known in advance. For fixed-budget simulations, this usually gives a cleaner spread than Halton for the same n_points.
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=37, sampling="hammersley"
)Use Hammersley when you will run one chosen sample count, and use Halton when you plan to grow the sample count progressively.
Gauss-Legendre sampling
Gauss-Legendre sampling uses quadrature nodes in the radial angle and uniform azimuth samples. It is the best choice when angular integration accuracy matters more than having exactly n_points samples.
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=36, sampling="gauss"
)sampling="gaussian_quadrature" is accepted as an alias for compatibility with YAML configs.
Legacy grid sampling
Grid sampling uses a uniform
cone = ConeIllumination(
cra_deg=10.0, f_number=2.8, n_points=36, sampling="grid"
)
points = cone.get_sampling_points()
print(f"Grid sampling: {len(points)} points")Grid sampling produces n_theta x n_phi points where n_theta = sqrt(n_points) and n_phi = n_points / n_theta.
Use grid mainly as a diagnostic baseline. It oversamples the cone center and creates strong radial/azimuthal structure, so it can make convergence look better or worse than it really is.
Weighting functions
The weighting function models the intensity distribution across the pupil:
| Weight | Formula | Physical model |
|---|---|---|
uniform | Flat pupil illumination | |
cosine | Lambertian source / aplanatic lens | |
cos4 | Cos-fourth falloff at image plane | |
gaussian | Apodized pupil |
The default is cosine, which is appropriate for most camera lens systems.
# Compare different weighting functions
for wf in ["uniform", "cosine", "cos4", "gaussian"]:
cone = ConeIllumination(
cra_deg=0.0, f_number=2.0, n_points=37, weighting=wf
)
points = cone.get_sampling_points()
weights = [p[2] for p in points]
print(f"{wf:10s}: max_w={max(weights):.4f}, min_w={min(weights):.4f}")Ray-file cone averaging
ConeIllumination is a compact model when the cone can be described by CRA, F-number, and an analytic pupil weighting. A lens-design workflow often provides something more explicit: a ray file for multiple sensor positions. Each ray usually carries:
| Field | Meaning |
|---|---|
image_x, image_y | Sensor position where the ray bundle lands |
pupil_x, pupil_y | Pupil coordinate, with the chief ray at the pupil center |
theta_deg, phi_deg | Incident angle at the sensor |
intensity | Lens transmission / Fresnel / vignetting factor |
weight | Pupil-area or quadrature weight |
For a ray bundle at one sensor position, cone-averaged QE is:
where
There are two practical ways to get the angular response values:
| Strategy | Use when | Trade-off |
|---|---|---|
| Direct ray simulation | Few positions or few ray samples | Accurate at the exact rays, expensive when repeated across the sensor |
| Angular-grid interpolation | Many positions, many lens rays | Run a structured |
Angular-grid interpolation is usually the scalable choice for camera-level characterization. The angular grid should cover the full CRA/MRA range of the lens and should be refined wherever
def cone_average_from_ray_bundle(qe_lookup, rays, wavelength):
numerator = 0.0
denominator = 0.0
for ray in rays:
qe = qe_lookup.interpolate(
wavelength=wavelength,
theta_deg=ray["theta_deg"],
phi_deg=ray["phi_deg"],
)
ray_weight = ray["intensity"] * ray["weight"]
numerator += qe * ray_weight
denominator += ray_weight
return numerator / max(denominator, 1e-12)INFO
A ray-file workflow is conceptually compatible with Zemax/OpticStudio, custom Python ray tracers, or measured CRA/MRA maps. COMPASS should treat the file as an optical interface: it needs ray angles and weights, not a dependency on any specific lens-design tool.
Integrating with planewave solvers
To compute cone-illuminated QE, run a planewave simulation at each sampled angle and compute the weighted sum:
import numpy as np
from compass.sources.cone_illumination import ConeIllumination
from compass.solvers.base import SolverFactory
# Set up cone
cone = ConeIllumination(cra_deg=15.0, f_number=2.0, n_points=37, weighting="cosine")
points = cone.get_sampling_points()
# Create solver
solver = SolverFactory.create("torcwa", solver_config, device="cuda")
solver.setup_geometry(pixel_stack)
# Run planewave at each sample point
wavelength = 0.55
weighted_qe = {}
for theta_deg, phi_deg, weight in points:
solver.setup_source({
"wavelength": wavelength,
"theta": float(theta_deg),
"phi": float(phi_deg),
"polarization": "unpolarized",
})
result = solver.run()
for pixel_name, qe in result.qe_per_pixel.items():
if pixel_name not in weighted_qe:
weighted_qe[pixel_name] = 0.0
weighted_qe[pixel_name] += weight * float(np.mean(qe))
print("Cone-illuminated QE at 550 nm:")
for pixel_name, qe in weighted_qe.items():
print(f" {pixel_name}: QE = {qe:.3f}")Cone illumination with wavelength sweep
For a full spectral sweep under cone illumination, iterate over wavelengths and angular samples:
wavelengths = np.arange(0.40, 0.701, 0.01)
cone = ConeIllumination(cra_deg=15.0, f_number=2.0, n_points=37)
points = cone.get_sampling_points()
# Initialize QE storage
pixel_names = None
cone_qe = {}
for wl in wavelengths:
wl_qe = {}
for theta_deg, phi_deg, weight in points:
solver.setup_source({
"wavelength": float(wl),
"theta": float(theta_deg),
"phi": float(phi_deg),
"polarization": "unpolarized",
})
result = solver.run()
if pixel_names is None:
pixel_names = list(result.qe_per_pixel.keys())
for pn in pixel_names:
cone_qe[pn] = []
for pn in pixel_names:
if pn not in wl_qe:
wl_qe[pn] = 0.0
wl_qe[pn] += weight * float(np.mean(result.qe_per_pixel[pn]))
for pn in pixel_names:
cone_qe[pn].append(wl_qe[pn])
# Plot results
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 6))
wl_nm = wavelengths * 1000
for pn in pixel_names:
ax.plot(wl_nm, cone_qe[pn], label=pn)
ax.set_xlabel("Wavelength (nm)")
ax.set_ylabel("QE (cone illumination)")
ax.set_title(f"Cone Illumination QE (CRA={cone.cra_deg} deg, F/{cone.f_number})")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()Sampling convergence
Check that the number of sample points is sufficient by comparing results at increasing n_points:
for n in [7, 19, 37, 61, 91]:
cone = ConeIllumination(cra_deg=15.0, f_number=2.0, n_points=n)
points = cone.get_sampling_points()
# ... run weighted sum, record QE
print(f"n_points={n}: avg QE = ...")Typically, 37 points gives results within 1% of the fully converged integral for F/2.0 and below.
Lens area sweep: F-number vs QE and crosstalk
The illumination cone footprint area scales with F-number. Sweeping the F-number reveals how lens speed affects QE and optical crosstalk — a critical trade-off in CIS design.
Footprint area vs F-number
The footprint radius
| F-number | Footprint radius (um) | Area (um²) | |
|---|---|---|---|
| F/1.4 | 20.9 | 1.91 | 11.5 |
| F/2.0 | 14.5 | 1.29 | 5.3 |
| F/2.8 | 10.3 | 0.91 | 2.6 |
| F/4.0 | 7.2 | 0.63 | 1.2 |
| F/5.6 | 5.1 | 0.45 | 0.63 |
(Assuming stack height h = 5.0 um)
Running an F-number sweep
import numpy as np
from compass.sources.cone_illumination import ConeIllumination
from compass.solvers.base import SolverFactory
f_numbers = [1.4, 2.0, 2.8, 4.0, 5.6, 8.0]
wavelength = 0.55
cra_deg = 15.0
solver = SolverFactory.create("torcwa", solver_config, device="cuda")
solver.setup_geometry(pixel_stack)
results = {}
for fn in f_numbers:
cone = ConeIllumination(cra_deg=cra_deg, f_number=fn, n_points=37)
points = cone.get_sampling_points()
weighted_qe = {}
for theta_deg, phi_deg, weight in points:
solver.setup_source({
"wavelength": wavelength,
"theta": float(theta_deg),
"phi": float(phi_deg),
"polarization": "unpolarized",
})
result = solver.run()
for pixel_name, qe in result.qe_per_pixel.items():
if pixel_name not in weighted_qe:
weighted_qe[pixel_name] = 0.0
weighted_qe[pixel_name] += weight * float(np.mean(qe))
results[fn] = weighted_qe
print(f"F/{fn}: {weighted_qe}")Analyzing the trade-off
Faster lenses (lower F-number) collect more light, improving signal. But the wider cone also increases angular spread and crosstalk between adjacent pixels:
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# QE vs F-number for the green pixel
green_qe = [results[fn].get("green_tl", 0) for fn in f_numbers]
ax1.plot(f_numbers, green_qe, "go-", linewidth=2, markersize=8)
ax1.set_xlabel("F-number")
ax1.set_ylabel("QE (green pixel)")
ax1.set_title("QE vs F-number (550 nm)")
ax1.grid(True, alpha=0.3)
ax1.invert_xaxis()
# Crosstalk: ratio of non-target pixel QE to target pixel QE
crosstalk = []
for fn in f_numbers:
green = results[fn].get("green_tl", 1e-9)
red = results[fn].get("red_tr", 0)
xtalk = red / green * 100 # percentage
crosstalk.append(xtalk)
ax2.plot(f_numbers, crosstalk, "rs-", linewidth=2, markersize=8)
ax2.set_xlabel("F-number")
ax2.set_ylabel("Crosstalk (%)")
ax2.set_title("Green→Red crosstalk vs F-number")
ax2.grid(True, alpha=0.3)
ax2.invert_xaxis()
plt.tight_layout()Combined CRA + F-number sweep
For a full lens-area sensitivity study, sweep both CRA and F-number to build a 2D map:
cra_values = [0, 5, 10, 15, 20, 25, 30]
f_numbers = [1.4, 2.0, 2.8, 4.0]
wavelength = 0.55
qe_map = np.zeros((len(cra_values), len(f_numbers)))
for i, cra in enumerate(cra_values):
for j, fn in enumerate(f_numbers):
cone = ConeIllumination(cra_deg=cra, f_number=fn, n_points=37)
points = cone.get_sampling_points()
weighted_qe = 0.0
for theta_deg, phi_deg, weight in points:
solver.setup_source({
"wavelength": wavelength,
"theta": float(theta_deg),
"phi": float(phi_deg),
"polarization": "unpolarized",
})
result = solver.run()
weighted_qe += weight * float(
np.mean(result.qe_per_pixel.get("green_tl", [0]))
)
qe_map[i, j] = weighted_qe
# Plot as heatmap
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(qe_map, aspect="auto", origin="lower",
extent=[f_numbers[0], f_numbers[-1],
cra_values[0], cra_values[-1]])
ax.set_xlabel("F-number")
ax.set_ylabel("CRA (degrees)")
ax.set_title("Green pixel QE: CRA vs F-number")
plt.colorbar(im, ax=ax, label="QE")
plt.tight_layout()This 2D sweep helps identify the CRA and F-number limits for a target QE threshold — essential information for microlens design optimization.
Configuration via YAML
Cone illumination can be specified in the source config:
source:
type: "cone"
cone:
cra_deg: 15.0
f_number: 2.0
n_points: 37
sampling: "fibonacci"
weighting: "cosine"
wavelength:
mode: "sweep"
sweep: {start: 0.40, stop: 0.70, step: 0.01}
polarization: "unpolarized"Next steps
- ROI Sweep -- sweep cone illumination across the sensor
- Microlens & CRA Optimization cookbook -- CRA vs QE study
- First Simulation -- basic planewave setup