Developed by Deeven Seru
High-Performance Adaptive Optics C-Engine for Shack-Hartmann Wavefront Sensing
- What is Project Radius?
- The Problem: Why Atmospheric Turbulence is a Hard Engineering Challenge
- Scientific Foundations
- Project Scope
- Technical Approach and Architecture
- Implementation: The C-Engine in Detail
- Live Camera & Frame Grabber Integration (GenICam & Camera Link)
- Dataset Generation with OOPAO
- Calibration Pipeline
- Results and Achievements
- Unbiased Real-World Robustness Analysis
- Visualizations
- Running the Project
- Project Structure
- Future Work
- Academic References & Literature
Project Radius is a standalone, deterministic, ultra-low-latency Adaptive Optics (AO) processing pipeline built for space situational awareness and laboratory wavefront sensing. It ingests raw detector images from a Shack-Hartmann Wavefront Sensor (SH-WFS), extracts wavefront gradients, reconstructs the phase as Zernike polynomial coefficients, generates physical actuator stroke commands for a Deformable Mirror (DM), and characterizes the atmospheric turbulence — all inside a C-Engine that completes the end-to-end computation in 0.31 milliseconds per frame for 55 Zernike modes.
The system is deliberately built on first principles: no machine learning, no black-box inference. Every operation maps directly to an equation in optical physics or linear algebra. The result is a pipeline that is fully explainable, independently verifiable, and capable of running on modest embedded hardware without GPU acceleration.
When light — whether laser communication beams, directed energy, or starlight — travels through the Earth's atmosphere, it passes through turbulent parcels of air at different temperatures and pressures. Each parcel has a slightly different refractive index. The cumulative effect is that a flat, plane-parallel wavefront arrives at the aperture as a randomly distorted phase surface. This causes:
- Astronomical imaging: Stars appear blurred and smeared across dozens of pixels instead of a tight diffraction-limited point.
- Free-Space Optical (FSO) communications: Bit error rates increase catastrophically as beam wander and scintillation disrupt the signal.
- High-energy laser systems: The beam cannot be focused to its intended spot, wasting energy.
The distortion is not static. Atmospheric turbulence evolves on a timescale defined by the Atmospheric Coherence Time (
The brute-force difficulty is that this measurement-to-correction loop must happen at 100–1000 Hz, sustained indefinitely, on hardware that must fit inside an optical laboratory or telescope dome.
The SH-WFS subdivides the incoming telescope pupil using a microlens array (MLA). Each microlens samples a small subaperture of the pupil and focuses the local wavefront gradient into a spot on a CCD or CMOS detector. The displacement of that spot from a reference position (measured in pixels) is directly proportional to the average tilt of the wavefront over that subaperture:
where
For a 20x20 MLA, this gives up to 400 slope measurements (316 valid after excluding subapertures outside the circular pupil), two per subaperture (X and Y), yielding a slope vector of dimension 632.
Following the formalism introduced by Noll (1976), the reconstructed wavefront phase is expanded in Zernike polynomials
| Zernike Mode | Physical Meaning |
|---|---|
| Piston (global phase offset, not correctable) | |
|
|
Tip and Tilt (beam pointing) |
| Defocus | |
|
|
Oblique / vertical Astigmatism |
|
|
Vertical / horizontal Coma |
|
|
Higher-order modes |
The coefficients
The Fried Parameter
Assuming Taylor's Frozen-Flow Hypothesis (Taylor, 1938), the Atmospheric Coherence Time
Project Radius covers the following end-to-end scope, matching space situational awareness and laboratory AO system requirements:
| Scope Item | Specification |
|---|---|
| Sensor Model | Shack-Hartmann WFS, 20 x 20 subaperture grid |
| Valid Subapertures | 316 of 400 (circular pupil mask) |
| Detector Resolution | 400 x 400 pixel detector per BMP frame |
| Telescope Aperture | 8-meter primary (simulation), configurable |
| Zernike Modes Reconstructed | 55 (Piston through 9th radial order) |
| Deformable Mirror | 357-actuator continuous facesheet |
| Turbulence Model | Von Karman ( |
| Dataset Size | 500 telemetry frames at 100 Hz |
| Latency Requirement | < 10 ms per frame |
| Accuracy Requirement | > 95% R² against ground truth |
| Output | Zernike coefficients, DM actuator strokes, |
The core constraint is determinism. A real-time control system cannot afford garbage-collection pauses, JIT compilation delays, or non-deterministic NumPy BLAS dispatches. All heavy computation is therefore implemented in C, compiled to a shared library (c_engine.so), and called from Python via ctypes with zero-copy pointer passing.
Python serves exclusively as the orchestration layer: loading files, initializing data structures, invoking the C-Engine, and computing summary statistics. It adds no latency to the real-time path.
graph TD
A["Raw .bmp telemetry frame\n(400x400 pixels, 8-bit)"]
B["Python: load frame, pass pointer\nvia zero-copy ctypes memory bridge"]
C["C-Engine: slopes.c\nCenter of Gravity centroiding"]
D["Slope vector S\n(632 elements: 316 x X,Y)"]
E["C-Engine: mvm_reconstructor.c\nS → Z via Z = G⁺ · S"]
F["Zernike coefficient vector Z\n(55 modes)"]
G["C-Engine: mvm_reconstructor.c\nZ → A via A = C_DM · Z"]
H["DM actuator stroke vector A\n(357 actuators)"]
I["Python: accumulate statistics\nr0, tau0, MSE, R²"]
A --> B --> C --> D --> E --> F --> G --> H
F --> I
| Component | Algorithm Chosen | Reason |
|---|---|---|
| Centroiding | Center of Gravity (CoG) | Direct intensity-weighted mean; O(N) per subaperture; no iterative solver |
| Reconstruction | Matrix-Vector Multiply (MVM) | Pre-inverted interaction matrix; single BLAS-level operation; deterministic |
| DM Mapping | MVM with influence function | Physically grounded; linear in the small-stroke regime |
|
|
Noll variance formula | Closed-form; no fitting required; uses already-computed Zernike coefficients |
|
|
Autocorrelation |
Simple, fast; requires only the time series of Tip/Tilt residuals |
All sensor parameters are packed into a strict LensletConfig struct defined in src/c_engine/geometry.h. This enforces contiguous memory layout and eliminates pointer arithmetic errors:
typedef struct {
int n_sub; // total subapertures (20x20 = 400)
int n_valid; // valid subapertures after pupil mask (316)
int sub_px; // pixels per subaperture side
float pixel_scale; // arcsec / pixel
int *valid_mask; // binary array, length n_sub
} LensletConfig;The valid_mask is loaded once from a pre-computed calibration CSV. During processing, the centroiding loop skips any subaperture where valid_mask[k] == 0, preventing any computation on dark or vignetted regions.
For each valid subaperture
The slope is the displacement from the pre-stored reference centroid. The inner loop is fully unrolled over fixed subaperture size, avoiding branch prediction misses. The result is a flat slope vector S[2 * n_valid] in row-major order (all X slopes followed by all Y slopes).
Two MVM operations run sequentially:
Step 1: Slope-to-Zernike reconstruction
Z[j] = sum_k( G_plus[j][k] * S[k] ) for j in 0..N_ZERNIKE
G_plus is the Moore-Penrose pseudo-inverse of the calibration interaction matrix of shape (n_zernike, n_slopes).
Step 2: Zernike-to-Actuator mapping
A[i] = sum_j( C_DM[i][j] * Z[j] ) for i in 0..N_ACTUATORS
C_DM is the DM influence function coupling matrix of shape (n_actuators, n_zernike).
Both loops are straight for loops with no dynamic allocation, no heap operations, and no system calls. The entire function runs on the stack.
To support integration with science-grade hardware, Project Radius implements a Hardware Abstraction Layer (HAL) (camera_interface.py) supporting three acquisition modes:
SimulatedCameraInterface: A thread-safe background frame generator running a physical optics propagation loop at a target 1000 Hz, simulating dynamic atmospheric turbulence.PlaybackCameraInterface: A high-speed test interface that pre-caches real-world/simulated telemetry frames in RAM to bypass disk I/O bottlenecks and streams them at a precise loop frequency (e.g., 1000 Hz) to test closed-loop timing.GenICamCameraInterface: A hardware client using the EMVA standardharvesterslibrary. It loads manufacturer GenTL common transport drivers (.ctilibraries) for CoaXPress, GigE, or Camera Link frame grabbers, starts acquisition, and retrieves buffers.
To achieve sub-millisecond latencies, the camera interface exposes the raw image buffer memory address directly to the C-Engine. In GenICam mode, the driver DMA buffer address (component.data) is passed as a raw ctypes float pointer directly to compute_slopes, completely bypassing standard Python memory copies.
Because no physical SH-WFS hardware was available for offline testing, a 500-frame ground-truth dataset was synthesized using OOPAO (Object-Oriented Python Adaptive Optics) — a rigorous, peer-validated adaptive optics simulation framework.
The physical parameters used to generate the dataset match the laboratory specifications:
Telescope : 8.0 m diameter, no central obstruction
Atmosphere : Von Karman turbulence, r0 = 15 cm, L0 = 25 m
Wind speed : [10, 5] m/s across 2 independent layers
SH-WFS : 20x20 microlens array, 20 pixels per subaperture
Detector : 400x400 pixel detector (20 subapertures x 20 px)
DM : 357 actuators, 35% mechanical coupling coefficient
Frame rate : 100 Hz, total duration = 5 seconds
Noise : Shot noise + readout noise (RON = 1.5 e-)
Before the C-Engine can process science frames, two calibration matrices must be generated. This is done once offline:
The deformable mirror is poked with a known unit stroke on each actuator in sequence. For each actuator poke, the resulting slope pattern on the WFS is recorded. The full set of slope responses forms the Interaction Matrix (n_slopes, n_actuators). The Pseudo-Inverse (or Control Matrix)
This matrix maps slopes directly to Zernike coefficients at runtime with a single matrix multiply.
The influence function of the DM is projected onto the Zernike basis, giving a matrix C_DM of shape (n_actuators, n_zernike). This allows the C-Engine to convert a Zernike coefficient vector into the specific analog voltage strokes required to produce that wavefront shape on the DM surface.
We evaluated the compiled C-Engine in real-time playback mode at 1000 Hz over all 55 Zernike modes:
| Metric | Target Requirement | Achieved (55 Modes) | Status |
|---|---|---|---|
| End-to-End Frame Latency | < 10.00 ms | 0.309 ms | Pass |
| Average Loop Rate | 1000 Hz | 998.05 Hz | Pass |
| Average Per-Frame Spatial |
> 95.00% | 99.48% | Pass |
| Global Temporal |
> 95.00% | 98.60% | Pass |
| Reconstruction MSE | < 0.1 | 1.14e-15 | Pass |
| Reconstructed Strehl Ratio | Maximize | 93.26% – 96.76% | Pass |
The C-Engine completes centroiding, 55-mode Zernike reconstruction, and DM actuator mapping in 0.309 milliseconds, representing a 32x safety margin over the 10 ms real-world requirement.
In our previous 20-mode system, the high-spatial-frequency phase structures of the atmosphere that could not be modeled by the first 20 modes were either lost or projected back onto lower modes as noise (spatial aliasing). By scaling the reconstructor to 55 modes, the pipeline captures these high-order structures accurately, raising the average shape matching R² to 99.48%.
To evaluate how the 55-mode reconstruction holds up under realistic physical disturbances, we ran a rigorous robustness analysis (scripts/robustness_analysis.py) under three simulated real-world conditions:
We simulated guiding on stars of varying brightness by adding Poisson noise.
-
$100,000$ photons/subap: Global$R^2$ = 98.60% | Spatial$R^2$ = 99.49% -
$1,000$ photons/subap: Global$R^2$ = 98.41% | Spatial$R^2$ = 99.42% -
$200$ photons/subap: Global$R^2$ = 97.71% | Spatial$R^2$ = 99.14% - Insight: The modal Zernike reconstructor is extremely robust to photon shot noise. Even on faint stars (200 photons/subap), the wavefront shape matching remains at 99.14%.
At a faint guide star level (1000 photons/subap), we added varying levels of pixel readout noise.
-
0.0
$e^-$ RON: Global$R^2$ = 98.43% | Spatial$R^2$ = 99.42% -
1.0
$e^-$ RON (sCMOS): Global$R^2$ = 94.04% | Spatial$R^2$ = 97.80% -
3.0
$e^-$ RON (EMCCD): Global$R^2$ = 65.42% | Spatial$R^2$ = 87.40% -
5.0
$e^-$ RON (CCD): Global$R^2$ = 32.40% | Spatial$R^2$ = 75.29% - Insight: WFS centroiding is highly sensitive to readout noise. Adding random pixel fluctuations pulls the Center of Gravity (CoG) centroids toward the center of each subaperture window. This demonstrates why physical systems require low-noise sCMOS/EMCCD sensors and background thresholding.
We shifted the entire WFS lenslet spot grid relative to the detector array to simulate mechanical vibration/thermal drift.
-
0.00 px shift: Global
$R^2$ = 98.60% | Spatial$R^2$ = 99.49% -
0.05 px shift: Global
$R^2$ = -6.45% | Spatial$R^2$ = 59.65% -
0.10 px shift: Global
$R^2$ = -314.51% | Spatial$R^2$ = -57.99% (Failed) - Insight: The modal reconstructor is exceedingly sensitive to alignment drift. Shifting the image by just 0.05 pixels degrades accuracy to 59.65%. This is because a uniform subaperture shift acts as a massive artificial Tip/Tilt error. This highlights the absolute necessity of active optical stabilization or real-time reference slope updating.
A live playback of the 55-mode Zernike reconstructor translating Shack-Hartmann spot displacements into a continuous 3D phase surface over the telescope pupil.
The corresponding physical actuator strokes on the 357-actuator continuous facesheet Deformable Mirror, computed directly from the Zernike coefficients via the DM influence coupling matrix.
The focal-plane image produced by the microlens array. Each bright spot corresponds to one valid subaperture. Displaced spots indicate local wavefront tilt due to turbulence. The C-Engine processes this raw image and extracts 316 spot positions via Center of Gravity.
The comparison tracking chart below shows the C-Engine live output (dashed lines) overlaying the ground-truth coefficients (solid black lines) for the Tip (Z2) and Tilt (Z3) modes over a 100-frame sequence:
The performance degradation curves for photon noise, readout noise, and sub-pixel alignment shift:
This visual reference library contains all 55 Zernike modes reconstructed by the C-Engine (radial orders
- Python 3.9+
gcc(any version supporting C99)- Unix-like OS (Linux, macOS)
Clone this repository with its submodules (required for OOPAO and MSHWFS dependencies):
git clone --recurse-submodules https://github.com/Deeven-Seru/project-radius.git
cd project-radiusWe recommend creating an isolated virtual environment.
python3 -m venv venv
source venv/bin/activate
pip install -r requirements.txt(Note: If testing physical camera hardware, you must also install the GenICam wrapper: pip install harvesters)
The C-Engine must be compiled locally to create the zero-copy shared library. The Makefile automatically detects your system architecture (Apple Silicon ARM64 vs Intel x86_64). On x86_64 systems, it will automatically enable AVX2/FMA vectorization for sub-millisecond latency.
cd src/c_engine
make
cd ../..This produces build/c_engine.so.
Because large arrays and datasets are excluded from version control, you must synthesize the ground-truth OOPAO dataset and calibration matrices.
python scripts/export_gplus.py # generates 55-mode g_plus.csv + dm_coupling.csv
python scripts/generate_dataset.py # generates 500 400x400 .bmp frames + ground_truth.csvLive Playback Benchmark (1000 Hz real data test):
python scripts/live_camera_test.py --mode playback --frames 500 --fps 1000.0Live Physical Simulation:
python scripts/live_camera_test.py --mode sim --frames 200Robustness Analysis & Verification:
python scripts/robustness_analysis.py # Evaluates noise/drift resilience
python scripts/compare_outputs.py # Computes spatial/temporal R² metrics.
├── data/
│ ├── dataset/
│ │ ├── frame_0000.bmp ... frame_0499.bmp # 500 400x400 SH-WFS detector frames
│ │ ├── ground_truth.csv # OOPAO Zernike coefficients (500 x 55)
│ │ ├── valid_mask.csv # Binary subaperture validity mask (20x20)
│ │ ├── ref_slopes.csv # Flat wavefront reference slopes
│ │ ├── g_plus.csv # Pre-inverted interaction matrix (55 x 632)
│ │ └── dm_coupling.csv # DM influence function matrix (357 x 55)
│ └── comparisons/
│ ├── zernike_comparison.png # Tracking comparison chart
│ └── robustness_analysis.png # Noise/drift performance charts
│
├── src/
│ ├── core/
│ │ └── camera_interface.py # Simulated, Playback, and GenICam HAL
│ └── c_engine/
│ ├── Makefile # gcc build rules: produces build/c_engine.so
│ ├── geometry.h # LensletConfig struct definition
│ ├── slopes.c # Center of Gravity centroiding algorithm
│ └── mvm_reconstructor.c # Zernike MVM and DM actuator MVM
│
├── scripts/
│ ├── live_camera_test.py # Real-time console test harness (sim/playback/hardware)
│ ├── compare_outputs.py # Unbiased side-by-side accuracy report
│ ├── robustness_analysis.py # Noise/drift simulation sweep script
│ ├── generate_dataset.py # OOPAO physics simulation and BMP export
│ ├── export_gplus.py # Interaction matrix calibration and export
│ └── inspect_eris.py # Real FITS telemetry data inspector
│
├── tests/
│ └── test_genicam_interface.py # Unit tests verifying GenICam wrapper API interactions
│
├── build/
│ └── c_engine.so # Compiled shared library
│
└── README.md
The following extensions are identified for subsequent phases:
- Closed-Loop DM Control: Integrate the actuator stroke output into a real-time closed-loop DM driver. Implement an integrator gain controller to drive the residual wavefront error to zero over successive frames.
-
Background subtraction thresholding: Add a configurable threshold filter in
slopes.cto reject background pixels, mitigating readout noise degradation. -
Multi-core BLAS/OpenMP Acceleration: Link the C-Engine MVM loops to OpenMP or optimized BLAS (like OpenBLAS/MKL) to handle larger configurations (
$80 \times 80$ subapertures) under a sub-millisecond budget.
The mathematical and physical models implemented in the C-Engine are drawn from the foundational literature in adaptive optics:
-
Fried, D. L. (1966). "Optical Resolution Through a Randomly Inhomogeneous Medium for Very Long and Very Short Exposures." Journal of the Optical Society of America, 56(10), 1372-1379. (Derivation of the atmospheric coherence length
$r_0$ ). - Noll, R. J. (1976). "Zernike polynomials and atmospheric turbulence." Journal of the Optical Society of America, 66(3), 207-211. (Theoretical variances of Zernike coefficients over Kolmogorov turbulence).
- Hardy, J. W. (1998). Adaptive Optics for Astronomical Telescopes. Oxford University Press. (Comprehensive review of Shack-Hartmann WFS and control matrix generation).
- Roddier, F. (1999). Adaptive Optics in Astronomy. Cambridge University Press. (SVD truncation strategies for interaction matrices).
-
Taylor, G. I. (1938). "The Spectrum of Turbulence." Proceedings of the Royal Society of London. Series A, 164(919), 476-490. (Frozen-flow hypothesis used in estimating the coherence time
$\tau_0$ ).



























































