mirror of
https://github.com/ruvnet/RuView
synced 2026-07-31 18:51:42 +00:00
feat(wifi-densepose-sar): coherent wideband RF tomography research crate (ADR-283)
New standalone leaf crate implementing the synthetic-aperture-radar reconstruction primitive a handheld through-wall RF imaging device would need: a stepped-frequency multi-position complex forward measurement simulator, delay-and-sum backprojection reconstruction, point-cloud extraction, and closed-form range/cross-range resolution + antenna-pose coherence-budget formulas checked against the reconstruction's actual behavior in tests/physics_validation.rs. Motivated by comparing this repo against Applied Electrodynamics' "WaveSight" launch. Scoped explicitly below ADR-278's RISE/DiffRadar/ GeRaF reproduction gates: this is the bare measurement-model + backprojection primitive, not a reproduction of any published system or a claim about real hardware capability. Every number is SYNTHETIC/L0 (ADR-282) -- no wideband RF hardware backs this crate. 24 tests (21 unit + 3 integration), 0 failed, clippy-clean. Adds a tutorial walkthrough and MEASURED backprojection benchmark numbers.
This commit is contained in:
ruv
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
|
||||
|
||||
### Added
|
||||
|
||||
- **`wifi-densepose-sar` — coherent wideband RF tomography research crate (ADR-283).** New standalone leaf crate (the `nvsim` pattern; zero coupling to `wifi-densepose-hardware` or any real ingestion path) implementing the synthetic-aperture-radar reconstruction primitive a handheld through-wall RF imaging device would need — motivated by comparison against Applied Electrodynamics' "WaveSight" launch, and explicitly scoped below ADR-278's RISE/DiffRadar/GeRaF reproduction gates. Ships: (1) a stepped-frequency, multi-position complex forward measurement simulator (`y_{m,k} = Σ σ_j/R² · exp(-i·4π·f·R/c) + noise`, deterministic ChaCha20 seeding); (2) delay-and-sum backprojection reconstruction onto a 3D voxel grid, rayon-parallelized over voxels; (3) threshold + local-maximum point-cloud extraction; (4) closed-form range/cross-range resolution and antenna-pose coherence-budget formulas (`ΔR=c/2B`, `δ_CR≈λR/2L`, `Δp≤λ/8`) checked against the reconstruction's *actual* behavior in `tests/physics_validation.rs` rather than merely documented — forward-simulating two targets at controlled separations and proving they resolve or merge exactly where the formulas predict, and that reconstructed focus at a known target degrades as injected antenna-pose error grows. Every number is SYNTHETIC/L0 (ADR-282) — no real wideband RF hardware backs this crate; see the crate README and `docs/tutorials/coherent-rf-tomography-backprojection.md` for the full honesty boundary and a worked walkthrough. 24 tests (21 unit + 3 integration), 0 failed, clippy-clean; MEASURED backprojection throughput ~350-450K voxels/sec (criterion, 21 poses × 32 freq steps).
|
||||
- **HOMECORE platform runtime completion — secure native/Wasmtime plugins, authenticated HAP IP, expanded Home Assistant APIs, durable restoration/migration, and voice protocols.** `homecore-server` now owns deterministic compiled-in native plugin registration plus explicitly configured, path-bounded, Ed25519 publisher-verified Wasm packages executed through Wasmtime with setup/state-change/teardown lifecycle; arbitrary native dynamic libraries remain intentionally unsupported. The optional HAP server implements persisted accessory identity and controller records, SRP-6a Pair-Setup M1–M6, X25519/Ed25519 Pair-Verify M1–M4, HKDF-SHA512/ChaCha20-Poly1305 record framing, authenticated/admin endpoint gates, replay/tamper closure, live entity synchronization, and paired-state `_hap._tcp` mDNS updates (45 focused tests; external Apple certification is not claimed). Startup restores device/entity registries and deterministic latest recorder states before plugins, and migration now atomically preserves forward-compatible device/config-entry fields. The HA-compatible surface adds events, templates, config checks, components, registries, history/logbook with SQL-enforced global response bounds, calendar/camera provider routes, and modern WebSocket negotiation while retaining a machine-readable limitations matrix for integration-specific behavior. Assist adds bounded PCM16, async STT/TTS contracts, an end-to-end speech pipeline, and an authenticated satellite session protocol; real deployments still provide the speech engines.
|
||||
- **`ruview-unified` increment 3 — Gaussian update-loop completion, separable delay-Doppler, and property-tested boundary hardening.** (1) `GaussianMap::merge_overlapping` (ADR-275 step 5: mutual-Mahalanobis + semantic-compatibility dedup catching drift the insert-time gate misses) and lifetime-aware decay (`τ_eff = τ·(1+ln(1+lifetime/τ))` — confirmed structures outlive transients at equal nominal τ). (2) `delay_doppler_map` reimplemented separably (`O(B²S+S²B)`), proven equivalent to the direct reference to <1e-10 and **measured 8.3× faster** (520 µs vs 4.34 ms at 56×8). (3) `tests/security_boundaries.rs` — 8 `proptest` properties over the boundary surfaces (arbitrary values incl. NaN/±inf via `f64::from_bits`) that found and fixed three input-controlled defects: a BLE-CS phase-unwrap infinite loop on non-finite phases and an ~1e299-iteration loop on finite-huge phases (now O(1) modular unwrap + plausibility bound), and a subnormal Gaussian scale overflowing `1/σ²` to NaN density (now physical σ/occupancy bounds). (4) New criterion benches for all increment-2 hot paths (`to_canonical` 38 µs, `ble_cs_range` 481 ns, AoI planner 647 ns/200 regions, coherent fusion 1.5 µs/32 members, factorized pose 521 ns). ruview-unified now 98 tests (87 lib + 3 acceptance + 8 security), 0 failed, clippy-clean.
|
||||
- **`ruview-unified` increment 2 — native frame contract + programmable perception (ADR-279..282).** (1) `RfFrameV2` becomes the authoritative RF record: native complex IQ with explicit validity masks, declared `PhaseState`, TX/RX poses + antenna geometry in one building frame, calibration/quality state, and a provenance rule enforced at construction — `Synthetic ⇒ L0Simulation` and `Measured ⇒ ≥ L1CapturedReplay` can never alias (the public L0–L5 evidence ladder is now a type); the 56-bin canonical tensor is demoted to a derived compatibility view (`to_canonical`, mask-aware gap-filling through the same normalization path as every adapter; native samples proven byte-untouched). (2) Active sensing control plane (`control.rs`): ETSI-ISAC-vocabulary `SensingTask` admission (raw export always refused; identity requires consent), `SensingAction`/`InformationGoal`, an age-of-information `ActiveSensingPlanner` (priority = uncertainty × change rate × criticality ÷ cost; **measured 95% sensing-traffic reduction** vs uniform refresh on a 20-region scenario), fail-closed `CoherentSensorGroup` fusion gates (time/phase/geometry bounds; five denial paths tested), policy-authorized RIS/movable-antenna actuation receipts, and purpose-scoped `TaskSufficientRepresentation` leakage validation. (3) New modality surfaces: BLE Channel Sounding adapter + `ble_cs_range` treating phase-slope and RTT as **separate cross-validated evidence** (exact distance recovery on synthetic tones; relay-style divergence flagged, never averaged), delay-Doppler-native `FieldAxis` + `delay_doppler_map` (unit-peak tone test), IEEE P3162 synthetic-aperture import profile. (4) RePos-factorized pose head (relative skeleton on the content representation, root on the geometry-conditioned one, calibrated per-joint uncertainties): held-out-room MPJPE 0.0003 m vs 0.2534 m for the monolithic baseline in the room-shortcut leakage experiment; ≤2% structured-adapter budget (740 params). (5) Age gate input now `log(1+age_ms)` per the age-aware-CSI recipe (gradient check re-proven); Gaussian primitives gained `first_seen_ns`/`doppler_variance`/bounded `source_receipts` lineage; `PartitionKey` gained a `session` dimension and `SplitManifest` certifies disjointness across all seven dimensions. 87 tests, 0 failed; crate clippy-clean. Docker images unaffected (no shipped binary consumes the crate yet); Python proof re-verified PASS.
|
||||
|
||||
@@ -0,0 +1,56 @@
|
||||
# ADR-283: `wifi-densepose-sar` — coherent wideband RF tomography research crate
|
||||
|
||||
| Field | Value |
|
||||
|-------|-------|
|
||||
| **Status** | Accepted — implemented (P1) |
|
||||
| **Date** | 2026-07-30 |
|
||||
| **Parent** | ADR-278 (radar inverse rendering research program), ADR-282 (mandatory L0–L5 evidence ladder) |
|
||||
| **Relates to** | ADR-273/274 (`ruview-unified`'s `FmcwRadarCube` adapter, the eventual integration point), ADR-275 (`GaussianMap`, ditto) |
|
||||
|
||||
## 0. PROOF discipline
|
||||
|
||||
Every accuracy number this crate produces is **SYNTHETIC / evidence level L0** (ADR-282): generated by the crate's own forward simulator (`measurement::simulate_measurement`), scored against its own known ground truth (`ScatteringTarget` positions). Nothing here has been validated against real wideband RF hardware, and the crate contains no such hardware integration.
|
||||
|
||||
## 1. Context
|
||||
|
||||
A YC-backed company, Applied Electrodynamics ("WaveSight"), publicly launched a handheld "camera that can see through walls" using undisclosed radio-imaging technology. Comparing it against this repo's capabilities surfaced a real gap: `wifi-densepose-signal::ruvsense::tomography` implements *radio tomographic imaging* (Wilson & Patwari 2010) — RSS-based shadowing attenuation on a fixed-link topology, no coherent phase, no multi-frequency stepping, no synthetic aperture. It is a different technique from what a SAR-style through-wall imager needs: coherent, wideband, multi-position backprojection.
|
||||
|
||||
`ruview-unified`'s `FmcwRadarCube` adapter (ADR-274) already normalizes wideband radar cubes into range profiles per position, and ADR-278 already names a radar-cube-output extension of the ADR-276 synthetic world generator as the intended sandbox for any future radar-inverse research. Neither, before this ADR, contained an actual backprojection reconstruction kernel — the primitive every candidate technique (matched-filter SAR, GPR imaging, RISE/DiffRadar-style inversion) is built on.
|
||||
|
||||
## 2. Decision
|
||||
|
||||
Ship `wifi-densepose-sar` as a standalone leaf crate (the `nvsim` pattern: pure Rust, deterministic ChaCha20 seeding, zero coupling to `wifi-densepose-hardware` or any real ingestion path) implementing:
|
||||
|
||||
1. **Forward measurement model** (`measurement.rs`): simulates the complex, stepped-frequency returns a monostatic synthetic-aperture radar would record from known point scatterers — `y_{m,k} = Σ_j σ_j/R_{m,j}² · exp(-i·4π·f_k·R_{m,j}/c) + noise`.
|
||||
2. **Backprojection reconstruction** (`reconstruct.rs`): the matched-filter inverse of (1) onto a 3D voxel grid, parallelized over voxels (rayon).
|
||||
3. **Point-cloud extraction** (`pointcloud.rs`): threshold + local-maximum extraction from the dense voxel image.
|
||||
4. **Closed-form resolution/coherence formulas** (`resolution.rs`): `ΔR = c/2B` (range resolution), `δ_CR ≈ λR/2L` (cross-range/synthetic-aperture resolution), `Δp ≤ λ/8` (antenna-pose coherence budget, derived from a quarter-wavelength round-trip-path tolerance) — checked against the reconstruction's actual behavior in `tests/physics_validation.rs`, not merely documented.
|
||||
|
||||
This is deliberately scoped **one level below** ADR-278's RISE/DiffRadar/GeRaF reproduction program: it is the bare measurement-model + backprojection primitive, not a reproduction of any specific published system, and not a claim about Applied Electrodynamics' undisclosed product (their waveform, antenna count, bandwidth, and algorithm are unknown; this crate applies the same well-established SAR/GPR physics — see Skolnik, *Radar Handbook* — to synthetic data).
|
||||
|
||||
## 3. What this explicitly is NOT
|
||||
|
||||
- Not a hardware driver. No VNA/SDR/wideband-RF-frontend code exists anywhere in this crate or was added to `wifi-densepose-hardware`.
|
||||
- Not wired into `ruview-unified`'s `FmcwRadarCube` adapter or `GaussianMap`. That integration is real future work (§5), deliberately deferred so the reconstruction physics validates in isolation first — the same staging ADR-278 §2.3 already prescribes ("sandbox-first... before hardware").
|
||||
- Not a reproduction of RISE, DiffRadar, or GeRaF. ADR-278's gates (G1–G4) are untouched by this ADR.
|
||||
- Not a real-world through-wall imaging performance claim. The forward model is free-space propagation only — no multipath, no per-material attenuation, no antenna gain pattern, no receiver noise figure. Real-world performance depends on all of these.
|
||||
|
||||
## 4. Simplifications (honesty boundary)
|
||||
|
||||
- **Monostatic, not MIMO.** A single antenna acts as both transmitter and receiver at each synthetic-aperture position (the standard stripmap-SAR simplification), not a multi-element MIMO array. Extending to bistatic/MIMO `(m, n)` transmitter/receiver pairs is straightforward given the existing measurement-model structure but not implemented.
|
||||
- **Isotropic antenna, no gain pattern.** Every antenna position radiates/receives equally in all directions.
|
||||
- **Free-space propagation only.** No multipath, no material transmission/reflection/attenuation (contrast `ruview-unified::synth::room`'s Fresnel material model, which is narrowband-CW and not yet extended to wideband — a natural follow-up, §5).
|
||||
- **`1/R²` two-way amplitude falloff, no calibration.** Real receivers have finite dynamic range, noise figures, and require calibration against a known reference target; none of that is modeled.
|
||||
|
||||
## 5. Follow-up (not in this ADR's scope)
|
||||
|
||||
1. Extend `ruview-unified::synth::room`'s image-method ray tracer to emit wideband stepped-frequency multi-position cubes (per ADR-278 §2.3), and wire `wifi-densepose-sar::reconstruct` against that richer (multipath-aware) synthetic generator instead of the free-space-only model here.
|
||||
2. A `ruview-unified` integration adapter converting `ReflectivityImage`/`PointCloudPoint` output into `RfGaussian`/`GaussianMap` primitives (ADR-278 §2.4's stated integration contract).
|
||||
3. Bistatic/MIMO measurement model.
|
||||
4. Any of ADR-278's actual gated reproductions (RISE first), if and when that program proceeds — this crate would be a component, not a substitute.
|
||||
|
||||
## 6. Consequences
|
||||
|
||||
- The workspace gains a real (if intentionally scoped-down) coherent-imaging primitive where before there was none — useful groundwork for ADR-278 if that research program proceeds, and a direct, honest answer to "could this repo build a WaveSight-like device" (no, not without the hardware program described in the motivating comparison; yes, this is the reconstruction-algorithm groundwork such a program would need).
|
||||
- Zero risk to the existing `wifi-densepose-signal::ruvsense::tomography` (RSS-based RTI) code path or any production pipeline — this crate is not referenced by any of them.
|
||||
- 24 tests (21 unit + 3 integration physics-validation), 0 failed, clippy-clean. Criterion bench: MEASURED 512/4096/32768-voxel backprojection reconstruction throughput (see crate README for the numbers as last recorded).
|
||||
@@ -0,0 +1,306 @@
|
||||
# Coherent Wideband RF Tomography: Simulating and Reconstructing with `wifi-densepose-sar`
|
||||
|
||||
A walkthrough of the `wifi-densepose-sar` crate (ADR-283): simulating
|
||||
synthetic-aperture radar (SAR) style measurements and reconstructing a 3D
|
||||
reflectivity image from them via delay-and-sum backprojection.
|
||||
|
||||
**Estimated time:** 30 minutes.
|
||||
|
||||
**What you will build:** A small Rust program that simulates a handheld
|
||||
stepped-frequency radar sweep past a couple of point targets, reconstructs
|
||||
a 3D image from the resulting complex measurements, and extracts a sparse
|
||||
point cloud from it — then verifies the reconstruction's resolution
|
||||
against closed-form theory.
|
||||
|
||||
**Who this is for:** Rust developers comfortable with basic signal
|
||||
processing terminology (frequency, bandwidth, phase) who want to
|
||||
understand what a coherent RF imaging pipeline actually computes, or who
|
||||
are evaluating whether this crate is a useful building block for their own
|
||||
radar-imaging research.
|
||||
|
||||
---
|
||||
|
||||
## Table of Contents
|
||||
|
||||
1. [What This Is (and Isn't)](#1-what-this-is-and-isnt)
|
||||
2. [Prerequisites](#2-prerequisites)
|
||||
3. [The Physics in Five Minutes](#3-the-physics-in-five-minutes)
|
||||
4. [Your First Reconstruction](#4-your-first-reconstruction)
|
||||
5. [Range Resolution: Why Bandwidth Matters](#5-range-resolution-why-bandwidth-matters)
|
||||
6. [Cross-Range Resolution: Why You Need to Move the Antenna](#6-cross-range-resolution-why-you-need-to-move-the-antenna)
|
||||
7. [The Antenna-Pose Coherence Budget](#7-the-antenna-pose-coherence-budget)
|
||||
8. [Extracting a Point Cloud](#8-extracting-a-point-cloud)
|
||||
9. [Benchmarking Your Own Scenario](#9-benchmarking-your-own-scenario)
|
||||
10. [Where This Could Go Next](#10-where-this-could-go-next)
|
||||
11. [Troubleshooting](#11-troubleshooting)
|
||||
|
||||
---
|
||||
|
||||
## 1. What This Is (and Isn't)
|
||||
|
||||
This crate exists because of a real question: could this repo build
|
||||
something like [Applied Electrodynamics' WaveSight](https://www.ae-dyn.com/)
|
||||
— a handheld device that images through walls using radio waves? The
|
||||
honest answer, worked out in ADR-283, is **no, not as a hardware product**
|
||||
— that needs a custom coherent RF front end, a calibrated antenna array,
|
||||
and real-time reconstruction hardware, which is an 18–36 month, high
|
||||
six-to-seven-figure hardware engineering program, not a software change.
|
||||
|
||||
What *is* useful to build, and what this crate is, is the **reconstruction
|
||||
algorithm** such a device needs: given coherent, phase-preserving,
|
||||
stepped-frequency measurements recorded from several known antenna
|
||||
positions, recover the 3D locations of the things that reflected the
|
||||
signal. That's a well-understood problem (synthetic-aperture radar,
|
||||
ground-penetrating radar imaging, and microwave tomography all solve
|
||||
versions of it) with textbook closed-form math behind it.
|
||||
|
||||
Every number in this crate comes from its own **synthetic forward
|
||||
simulator** — there is no real radio hardware anywhere in this crate, and
|
||||
none of its tests, benchmarks, or accuracy numbers say anything about how
|
||||
well a real device would perform through a real wall. That's evidence
|
||||
level **L0 (Synthetic)** in this repo's [ADR-282](../adr/ADR-282-ruview-ecosystem-positioning.md)
|
||||
evidence ladder, and it stays L0 until (if ever) real wideband RF hardware
|
||||
feeds this pipeline real measurements.
|
||||
|
||||
## 2. Prerequisites
|
||||
|
||||
- Rust 1.75+ (workspace MSRV), already set up if you can build the rest of
|
||||
this repo's `v2/` workspace.
|
||||
- No special hardware. Everything in this tutorial runs from synthetic
|
||||
data.
|
||||
|
||||
```bash
|
||||
cd v2
|
||||
cargo test -p wifi-densepose-sar --no-default-features
|
||||
```
|
||||
|
||||
If that passes (24 tests, 0 failed), you're ready.
|
||||
|
||||
## 3. The Physics in Five Minutes
|
||||
|
||||
A stepped-frequency radar sweeps `K` frequencies `f_0..f_{K-1}` across a
|
||||
band of total width `B` (the bandwidth). At each of `M` antenna positions
|
||||
`p_0..p_{M-1}` along a handheld sweep, it records one complex number per
|
||||
frequency — amplitude and phase, not just amplitude, which is what makes
|
||||
this "coherent."
|
||||
|
||||
For a point scatterer at position `x` with reflectivity `σ`, range
|
||||
`R = |p_m - x|` from antenna position `m`, the forward model this crate
|
||||
simulates is:
|
||||
|
||||
```text
|
||||
y_{m,k} = sigma / R^2 * exp(-i * 4*pi * f_k * R / c)
|
||||
```
|
||||
|
||||
`4*pi*f*R/c` is the two-way (round-trip) propagation phase; `1/R^2` is the
|
||||
two-way free-space spreading loss. With several targets, the measurement
|
||||
is just the sum of each target's contribution (superposition — this crate
|
||||
never models multipath/interaction between targets, only free-space direct
|
||||
paths).
|
||||
|
||||
**Reconstruction (backprojection)** inverts this: for every candidate
|
||||
voxel `x` in a 3D grid, it multiplies each measurement by the *complex
|
||||
conjugate* of the phase the forward model would have applied for a target
|
||||
at `x`, then sums:
|
||||
|
||||
```text
|
||||
I(x) = | (1/MK) * sum_m sum_k y_{m,k} * R_{m,x}^2 * exp(+i * 4*pi * f_k * R_{m,x} / c) |
|
||||
```
|
||||
|
||||
If `x` coincides with a real target, every term's phase correction exactly
|
||||
cancels the phase the forward model applied — the sum adds up
|
||||
constructively ("coherent gain"). At any other voxel, the phases are
|
||||
essentially uncorrelated across the `(m, k)` grid and the sum averages
|
||||
toward zero. That's the entire algorithm: matched filtering, done in 3D,
|
||||
one voxel at a time.
|
||||
|
||||
## 4. Your First Reconstruction
|
||||
|
||||
Add `wifi-densepose-sar` to a scratch binary or run this in a workspace
|
||||
example. It simulates two targets, reconstructs, and finds the brightest
|
||||
voxel:
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::{
|
||||
backproject, linear_aperture, simulate_measurement, FrequencySweep,
|
||||
Point3, ScatteringTarget, VoxelGrid,
|
||||
};
|
||||
|
||||
fn main() {
|
||||
// A 1-meter handheld sweep, 21 antenna positions along it.
|
||||
let poses = linear_aperture(
|
||||
Point3::new(-0.5, 0.0, 0.0),
|
||||
Point3::new(0.5, 0.0, 0.0),
|
||||
21,
|
||||
);
|
||||
|
||||
// Sweep 2-6 GHz (4 GHz of bandwidth) in 32 steps.
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
|
||||
// One target, 2 meters downrange, reflectivity 1.0 (arbitrary units).
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
|
||||
// Simulate the measurement with a touch of noise (seeded -- rerunning
|
||||
// with the same seed gives byte-identical output).
|
||||
let measurement = simulate_measurement(&poses, &sweep, &[target], 0.01, 42);
|
||||
|
||||
// Reconstruct a 21x21x21 voxel grid around where we expect the target.
|
||||
let grid = VoxelGrid::new(Point3::new(-0.3, 1.7, -0.3), 0.03, 21, 21, 21);
|
||||
let image = backproject(&measurement, &poses, &sweep, &grid);
|
||||
|
||||
let (peak_location, peak_magnitude) = image.peak();
|
||||
println!("true target: {:?}", target.position);
|
||||
println!("reconstructed peak: {peak_location:?} (magnitude {peak_magnitude:.4})");
|
||||
}
|
||||
```
|
||||
|
||||
Run it and you should see the reconstructed peak within a couple of
|
||||
centimeters of the true target position — well inside the voxel spacing
|
||||
used here (3 cm). That's `tests/reconstruct.rs::single_point_target_reconstructs_at_its_true_location`
|
||||
running live.
|
||||
|
||||
## 5. Range Resolution: Why Bandwidth Matters
|
||||
|
||||
How close together can two targets be *along the same bearing* (same
|
||||
antenna, different distance) before they blur into one blob? The classic
|
||||
radar answer: `ΔR = c / (2B)` — resolution improves with more swept
|
||||
bandwidth, full stop. Carrier frequency, antenna count, and aperture
|
||||
length don't enter into it at all.
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::resolution::range_resolution_m;
|
||||
|
||||
let dr = range_resolution_m(4.0e9); // 4 GHz swept bandwidth
|
||||
println!("range resolution: {:.1} cm", dr * 100.0);
|
||||
// -> range resolution: 3.7 cm
|
||||
```
|
||||
|
||||
`tests/physics_validation.rs::range_separated_targets_resolve_only_beyond_range_resolution`
|
||||
proves this isn't just a formula sitting in a doc comment: it forward-simulates
|
||||
two targets 4x `ΔR` apart (they resolve into two distinct peaks) and 0.25x
|
||||
`ΔR` apart (they merge into one), using the *same* `range_resolution_m`
|
||||
call to pick the separations.
|
||||
|
||||
## 6. Cross-Range Resolution: Why You Need to Move the Antenna
|
||||
|
||||
A single antenna position, no matter how much bandwidth it sweeps, cannot
|
||||
tell two targets apart if they're at the same range but different bearing
|
||||
— all it measures is round-trip distance, which is the same for both. This
|
||||
is exactly why "handheld... sweep the antenna around" matters: moving the
|
||||
antenna across a synthetic aperture of length `L` gives you angular
|
||||
information, with cross-range resolution:
|
||||
|
||||
```text
|
||||
delta_CR ~= lambda * R / (2 * L)
|
||||
```
|
||||
|
||||
— finer with a longer aperture, a shorter wavelength (higher carrier
|
||||
frequency), or a closer target.
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::resolution::cross_range_resolution_m;
|
||||
|
||||
let short = cross_range_resolution_m(4.0e9, 0.05, 2.0); // 5cm sweep
|
||||
let long = cross_range_resolution_m(4.0e9, 1.0, 2.0); // 1m sweep
|
||||
println!("5cm aperture: {:.2} m cross-range resolution", short);
|
||||
println!("1m aperture: {:.2} m cross-range resolution", long);
|
||||
// -> a 20x longer aperture gives 20x finer cross-range resolution
|
||||
```
|
||||
|
||||
`tests/physics_validation.rs::cross_range_separated_targets_resolve_only_with_long_enough_aperture`
|
||||
demonstrates this end-to-end: the same pair of cross-range-separated
|
||||
targets resolves into two peaks with a 1m synthetic aperture and collapses
|
||||
into one with a 5cm aperture, no other change.
|
||||
|
||||
## 7. The Antenna-Pose Coherence Budget
|
||||
|
||||
Backprojection assumes you know exactly where the antenna was at each
|
||||
measurement. If your position tracking (in a real device: visual-inertial
|
||||
odometry, encoders, whatever) is off by `Δp`, the phase correction applied
|
||||
during reconstruction is wrong by an amount that grows with `Δp` and with
|
||||
frequency. The classical rule of thumb for "still well focused": keep the
|
||||
round-trip path error under a quarter wavelength, which works out to an
|
||||
antenna-position tolerance of `λ/8`:
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::resolution::max_coherent_pose_error_m;
|
||||
|
||||
let budget = max_coherent_pose_error_m(8.0e9); // 8 GHz carrier
|
||||
println!("position tolerance at 8 GHz: {:.1} mm", budget * 1000.0);
|
||||
// -> position tolerance at 8 GHz: 4.7 mm
|
||||
```
|
||||
|
||||
`tests/physics_validation.rs::phase_error_from_pose_jitter_degrades_focus_beyond_pose_budget`
|
||||
verifies this isn't just asserted: it perturbs the *true* antenna positions
|
||||
away from the *assumed* ones used in reconstruction, and shows focus at the
|
||||
true target location degrades as that perturbation grows — the concrete
|
||||
mechanism behind why real SAR/GPR imaging systems need accurate pose
|
||||
tracking, not just a good radio.
|
||||
|
||||
## 8. Extracting a Point Cloud
|
||||
|
||||
A dense voxel grid isn't a useful end product — you want a short list of
|
||||
detected points:
|
||||
|
||||
Continuing the program from §4 (which already has `image` in scope):
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::extract_point_cloud;
|
||||
|
||||
let points = extract_point_cloud(&image, 0.5); // 50%-of-peak threshold
|
||||
for p in &points {
|
||||
println!("{:?} magnitude={:.3}", p.position, p.magnitude);
|
||||
}
|
||||
```
|
||||
|
||||
`extract_point_cloud` does threshold + 6-connected local-maximum
|
||||
extraction — a real blob will still yield one point, not one per voxel
|
||||
inside it. There is deliberately no clustering, material classification,
|
||||
or confidence calibration here (ADR-283 §5): that needs real data to
|
||||
calibrate against, which this crate does not have.
|
||||
|
||||
## 9. Benchmarking Your Own Scenario
|
||||
|
||||
```bash
|
||||
cargo bench -p wifi-densepose-sar
|
||||
```
|
||||
|
||||
The shipped benchmark (`benches/backprojection_bench.rs`) sweeps 512 /
|
||||
4,096 / 32,768-voxel grids with 21 poses x 32 frequencies. Reconstruction
|
||||
is embarrassingly parallel over voxels (each voxel's cost is independent),
|
||||
so it's rayon-parallelized already — see the crate README for the last
|
||||
recorded MEASURED numbers on the reference machine.
|
||||
|
||||
## 10. Where This Could Go Next
|
||||
|
||||
This crate deliberately stops short of several things (ADR-283 §5):
|
||||
|
||||
- It's monostatic (one antenna, both TX and RX) — real handheld SAR/MIMO
|
||||
devices often use multiple simultaneous antenna elements.
|
||||
- The forward model is free-space only — no multipath, no per-material
|
||||
attenuation (contrast `ruview-unified`'s narrowband Fresnel material
|
||||
model, which isn't yet extended to wideband).
|
||||
- It isn't wired into `ruview-unified`'s `FmcwRadarCube` adapter or
|
||||
`GaussianMap` — ADR-278 names that as the eventual integration point,
|
||||
once (and if) a reconstruction system is ready for it.
|
||||
|
||||
If you're picking this up to extend it, start with ADR-283's "Follow-up"
|
||||
section rather than guessing at scope.
|
||||
|
||||
## 11. Troubleshooting
|
||||
|
||||
**"My reconstructed peak isn't near my target."** Check your voxel grid
|
||||
actually covers the target's true location — `backproject` happily
|
||||
reconstructs whatever region you ask for; if the target is outside the
|
||||
grid, you'll get whatever's brightest inside it instead (usually noise).
|
||||
|
||||
**"Two targets I expected to resolve didn't."** Compute
|
||||
`range_resolution_m`/`cross_range_resolution_m` for your actual bandwidth
|
||||
and aperture length and check your separation against them — resolution
|
||||
is a hard physical limit here, not a tuning parameter.
|
||||
|
||||
**"Backprojection is slow for my grid size."** Cost is
|
||||
`O(voxels x poses x freqs)` and already parallelized over voxels via
|
||||
rayon; the only way to go faster is fewer voxels, fewer poses, or fewer
|
||||
frequency steps (each is a hard tradeoff against resolution or aperture
|
||||
coverage — see §5/§6).
|
||||
Generated
+13
@@ -11627,6 +11627,19 @@ dependencies = [
|
||||
"thiserror 2.0.18",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "wifi-densepose-sar"
|
||||
version = "0.3.1"
|
||||
dependencies = [
|
||||
"criterion",
|
||||
"num-complex",
|
||||
"rand 0.8.5",
|
||||
"rand_chacha 0.3.1",
|
||||
"rayon",
|
||||
"serde",
|
||||
"thiserror 2.0.18",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "wifi-densepose-sensing-server"
|
||||
version = "0.3.5"
|
||||
|
||||
@@ -88,6 +88,12 @@ members = [
|
||||
# submodule crates (rufield-core/-provenance/-privacy/-fusion); single
|
||||
# coupling point between RuView and the standalone RuField MFS spec.
|
||||
"crates/wifi-densepose-rufield",
|
||||
# ADR-283 — coherent wideband RF tomography research crate: synthetic
|
||||
# stepped-frequency multi-position measurement simulation + delay-and-sum
|
||||
# backprojection reconstruction. Standalone leaf (nvsim pattern), zero
|
||||
# hardware coupling, every number SYNTHETIC/L0 until real wideband RF
|
||||
# hardware exists.
|
||||
"crates/wifi-densepose-sar",
|
||||
]
|
||||
# ADR-040: WASM edge crate targets wasm32-unknown-unknown (no_std),
|
||||
# excluded from workspace to avoid breaking `cargo test --workspace`.
|
||||
|
||||
@@ -0,0 +1,42 @@
|
||||
[package]
|
||||
name = "wifi-densepose-sar"
|
||||
description = "Coherent wideband RF tomography research crate (ADR-283): synthetic stepped-frequency multi-position measurement simulation + delay-and-sum backprojection reconstruction. SYNTHETIC/L0 only -- no wideband RF hardware backs this crate."
|
||||
version.workspace = true
|
||||
edition.workspace = true
|
||||
authors.workspace = true
|
||||
license.workspace = true
|
||||
repository.workspace = true
|
||||
documentation.workspace = true
|
||||
keywords = ["radar", "sar", "backprojection", "rf-tomography", "simulation"]
|
||||
categories = ["science", "simulation"]
|
||||
readme = "README.md"
|
||||
|
||||
# `wifi-densepose-sar` is a standalone leaf crate (the nvsim / ruview-unified
|
||||
# pattern): pure-Rust math, deterministic ChaCha20 randomness (same seed =>
|
||||
# byte-identical output on every machine), zero coupling to any hardware
|
||||
# ingestion path. It has NO internal RuView dependency -- see the crate-level
|
||||
# doc comment in `src/lib.rs` for why (ADR-283 §2): this validates the
|
||||
# reconstruction *algorithm* against synthetic ground truth before any
|
||||
# question of wiring it into `ruview-unified`'s `FmcwRadarCube` adapter or
|
||||
# `GaussianMap` is in scope.
|
||||
[dependencies]
|
||||
num-complex = { workspace = true }
|
||||
thiserror = { workspace = true }
|
||||
serde = { workspace = true, features = ["derive"] }
|
||||
rand = { version = "0.8", default-features = false }
|
||||
rand_chacha = { version = "0.3", default-features = false }
|
||||
rayon = "1.10"
|
||||
|
||||
[dev-dependencies]
|
||||
criterion = { workspace = true }
|
||||
|
||||
[[bench]]
|
||||
name = "backprojection_bench"
|
||||
harness = false
|
||||
|
||||
[lints.rust]
|
||||
unsafe_code = "forbid"
|
||||
missing_docs = "warn"
|
||||
|
||||
[lints.clippy]
|
||||
all = "warn"
|
||||
@@ -0,0 +1,62 @@
|
||||
# wifi-densepose-sar
|
||||
|
||||
Coherent wideband RF tomography research crate (ADR-283): synthetic
|
||||
stepped-frequency multi-position measurement simulation + delay-and-sum
|
||||
backprojection reconstruction of a 3D reflectivity field.
|
||||
|
||||
**This is not a hardware capability.** It is the reconstruction primitive a
|
||||
handheld through-wall RF imaging device would need, validated against its
|
||||
own synthetic ground truth. Every number this crate produces is
|
||||
SYNTHETIC / evidence level L0 (ADR-282) until real wideband RF hardware
|
||||
(a VNA, SDR, or purpose-built radar front end) exists to feed it real
|
||||
measurements. See the crate-level doc comment in `src/lib.rs` for the full
|
||||
honesty boundary, and the tutorial at
|
||||
`docs/tutorials/coherent-rf-tomography-backprojection.md` for a walkthrough.
|
||||
|
||||
## Quick example
|
||||
|
||||
```rust
|
||||
use wifi_densepose_sar::{
|
||||
backproject, linear_aperture, simulate_measurement, FrequencySweep,
|
||||
Point3, ScatteringTarget, VoxelGrid,
|
||||
};
|
||||
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
let measurement = simulate_measurement(&poses, &sweep, &[target], 0.01, 42);
|
||||
|
||||
let grid = VoxelGrid::new(Point3::new(-0.3, 1.7, -0.3), 0.03, 21, 21, 21);
|
||||
let image = backproject(&measurement, &poses, &sweep, &grid);
|
||||
let (peak_location, peak_magnitude) = image.peak();
|
||||
println!("reconstructed target near {peak_location:?}, magnitude {peak_magnitude:.4}");
|
||||
```
|
||||
|
||||
## Testing
|
||||
|
||||
```bash
|
||||
cargo test -p wifi-densepose-sar --no-default-features
|
||||
cargo bench -p wifi-densepose-sar
|
||||
```
|
||||
|
||||
`tests/physics_validation.rs` checks the reconstruction's actual behavior
|
||||
against the closed-form formulas in `resolution.rs` (range resolution,
|
||||
cross-range/synthetic-aperture resolution, and the antenna-pose coherence
|
||||
budget) rather than merely asserting them: 24 tests (21 unit + 3
|
||||
integration), 0 failed, clippy-clean.
|
||||
|
||||
## Performance (MEASURED)
|
||||
|
||||
`cargo bench -p wifi-densepose-sar`, 21 antenna poses × 32 frequency steps
|
||||
(672 measurement terms/voxel), rayon-parallelized over voxels, this
|
||||
machine, release profile:
|
||||
|
||||
| Voxels | Median time | Throughput |
|
||||
|-------:|------------:|-----------:|
|
||||
| 512 | 1.47 ms | ~348,000 voxels/s |
|
||||
| 4,096 | 10.4 ms | ~394,000 voxels/s |
|
||||
| 32,768 | 73.5 ms | ~446,000 voxels/s |
|
||||
|
||||
Scales as expected: each voxel's cost is independent (`O(poses × freqs)`
|
||||
per voxel, embarrassingly parallel), so throughput is roughly constant
|
||||
across grid sizes and total time scales linearly with voxel count.
|
||||
@@ -0,0 +1,38 @@
|
||||
//! Criterion benchmark for the backprojection reconstruction kernel.
|
||||
//! `cargo bench -p wifi-densepose-sar` reports MEASURED throughput --
|
||||
//! see the crate README for the last recorded numbers.
|
||||
#![allow(missing_docs)]
|
||||
|
||||
use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
|
||||
use wifi_densepose_sar::geometry::linear_aperture;
|
||||
use wifi_densepose_sar::measurement::{simulate_measurement, FrequencySweep, ScatteringTarget};
|
||||
use wifi_densepose_sar::reconstruct::{backproject, VoxelGrid};
|
||||
use wifi_densepose_sar::Point3;
|
||||
|
||||
fn backprojection_benchmark(c: &mut Criterion) {
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
let targets = vec![
|
||||
ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0),
|
||||
ScatteringTarget::new(Point3::new(0.3, 1.8, 0.1), 0.6),
|
||||
];
|
||||
let measurement = simulate_measurement(&poses, &sweep, &targets, 0.01, 7);
|
||||
|
||||
let mut group = c.benchmark_group("backproject");
|
||||
for &voxels_per_axis in &[8usize, 16, 32] {
|
||||
let grid = VoxelGrid::new(
|
||||
Point3::new(-0.5, 1.5, -0.5),
|
||||
1.0 / voxels_per_axis as f64,
|
||||
voxels_per_axis,
|
||||
voxels_per_axis,
|
||||
voxels_per_axis,
|
||||
);
|
||||
group.bench_with_input(BenchmarkId::from_parameter(grid.len()), &grid, |b, grid| {
|
||||
b.iter(|| backproject(&measurement, &poses, &sweep, grid));
|
||||
});
|
||||
}
|
||||
group.finish();
|
||||
}
|
||||
|
||||
criterion_group!(benches, backprojection_benchmark);
|
||||
criterion_main!(benches);
|
||||
@@ -0,0 +1,140 @@
|
||||
//! Antenna positions, synthetic-aperture trajectories, and point geometry.
|
||||
|
||||
use serde::{Deserialize, Serialize};
|
||||
|
||||
/// A point in 3D space, meters, in an arbitrary right-handed scene frame.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||||
pub struct Point3 {
|
||||
/// X coordinate, meters.
|
||||
pub x: f64,
|
||||
/// Y coordinate, meters.
|
||||
pub y: f64,
|
||||
/// Z coordinate, meters.
|
||||
pub z: f64,
|
||||
}
|
||||
|
||||
impl Point3 {
|
||||
/// Construct a point.
|
||||
pub fn new(x: f64, y: f64, z: f64) -> Self {
|
||||
Self { x, y, z }
|
||||
}
|
||||
|
||||
/// Euclidean distance to another point, meters.
|
||||
pub fn distance(&self, other: &Point3) -> f64 {
|
||||
let dx = self.x - other.x;
|
||||
let dy = self.y - other.y;
|
||||
let dz = self.z - other.z;
|
||||
(dx * dx + dy * dy + dz * dz).sqrt()
|
||||
}
|
||||
|
||||
/// Unit vector pointing from `self` toward `other`. Returns `None` if
|
||||
/// the two points coincide (distance below `f64::EPSILON`).
|
||||
pub fn direction_to(&self, other: &Point3) -> Option<Point3> {
|
||||
let d = self.distance(other);
|
||||
if d < f64::EPSILON {
|
||||
return None;
|
||||
}
|
||||
Some(Point3::new(
|
||||
(other.x - self.x) / d,
|
||||
(other.y - self.y) / d,
|
||||
(other.z - self.z) / d,
|
||||
))
|
||||
}
|
||||
|
||||
/// Translate this point by `dist` meters along a unit vector `dir`.
|
||||
pub fn translated(&self, dir: Point3, dist: f64) -> Point3 {
|
||||
Point3::new(
|
||||
self.x + dir.x * dist,
|
||||
self.y + dir.y * dist,
|
||||
self.z + dir.z * dist,
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
/// A single antenna position along a synthetic-aperture trajectory.
|
||||
///
|
||||
/// Only position is modeled (an isotropic-antenna approximation, ADR-283
|
||||
/// §4) -- no antenna gain pattern / boresight direction is applied to the
|
||||
/// forward measurement model.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||||
pub struct AntennaPose {
|
||||
/// Antenna phase-center position, meters.
|
||||
pub position: Point3,
|
||||
}
|
||||
|
||||
impl AntennaPose {
|
||||
/// Construct a pose at `position`.
|
||||
pub fn new(position: Point3) -> Self {
|
||||
Self { position }
|
||||
}
|
||||
}
|
||||
|
||||
/// Generate `n` antenna poses evenly spaced along a straight line segment
|
||||
/// from `start` to `end` (inclusive), the canonical "handheld linear sweep"
|
||||
/// synthetic aperture. `n` must be >= 2 for a non-degenerate aperture.
|
||||
pub fn linear_aperture(start: Point3, end: Point3, n: usize) -> Vec<AntennaPose> {
|
||||
if n == 0 {
|
||||
return Vec::new();
|
||||
}
|
||||
if n == 1 {
|
||||
return vec![AntennaPose::new(start)];
|
||||
}
|
||||
(0..n)
|
||||
.map(|i| {
|
||||
let t = i as f64 / (n - 1) as f64;
|
||||
AntennaPose::new(Point3::new(
|
||||
start.x + (end.x - start.x) * t,
|
||||
start.y + (end.y - start.y) * t,
|
||||
start.z + (end.z - start.z) * t,
|
||||
))
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// The physical length of a synthetic aperture: the distance between its
|
||||
/// first and last pose. Used by [`crate::resolution::cross_range_resolution`].
|
||||
pub fn aperture_length(poses: &[AntennaPose]) -> f64 {
|
||||
match (poses.first(), poses.last()) {
|
||||
(Some(a), Some(b)) => a.position.distance(&b.position),
|
||||
_ => 0.0,
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn linear_aperture_spans_endpoints() {
|
||||
let start = Point3::new(0.0, 0.0, 0.0);
|
||||
let end = Point3::new(1.0, 0.0, 0.0);
|
||||
let poses = linear_aperture(start, end, 5);
|
||||
assert_eq!(poses.len(), 5);
|
||||
assert_eq!(poses[0].position, start);
|
||||
assert_eq!(poses[4].position, end);
|
||||
// Evenly spaced: 0, 0.25, 0.5, 0.75, 1.0 along x.
|
||||
assert!((poses[2].position.x - 0.5).abs() < 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn aperture_length_matches_endpoint_distance() {
|
||||
let poses = linear_aperture(Point3::new(0.0, 0.0, 0.0), Point3::new(3.0, 4.0, 0.0), 10);
|
||||
assert!((aperture_length(&poses) - 5.0).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn direction_to_is_unit_length() {
|
||||
let a = Point3::new(0.0, 0.0, 0.0);
|
||||
let b = Point3::new(2.0, 0.0, 0.0);
|
||||
let dir = a.direction_to(&b).unwrap();
|
||||
assert!((dir.x - 1.0).abs() < 1e-12);
|
||||
let len = (dir.x * dir.x + dir.y * dir.y + dir.z * dir.z).sqrt();
|
||||
assert!((len - 1.0).abs() < 1e-12);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn direction_to_coincident_points_is_none() {
|
||||
let a = Point3::new(1.0, 1.0, 1.0);
|
||||
assert!(a.direction_to(&a).is_none());
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,77 @@
|
||||
//! # wifi-densepose-sar — Coherent Wideband RF Tomography (ADR-283)
|
||||
//!
|
||||
//! A research crate implementing the reconstruction primitive that a
|
||||
//! handheld through-wall RF imaging device (the class of product exemplified
|
||||
//! by YC-backed Applied Electrodynamics' "WaveSight") would need: recovering
|
||||
//! a 3D reflectivity field from coherent, stepped-frequency, multi-position
|
||||
//! RF measurements via delay-and-sum backprojection.
|
||||
//!
|
||||
//! ## What this crate is
|
||||
//!
|
||||
//! - [`measurement`]: a forward simulator producing synthetic complex,
|
||||
//! stepped-frequency returns from known point scatterers observed by a
|
||||
//! known synthetic-aperture antenna trajectory.
|
||||
//! - [`reconstruct`]: the backprojection reconstruction kernel that inverts
|
||||
//! that forward model back into a 3D voxel reflectivity image.
|
||||
//! - [`pointcloud`]: sparse point extraction from the dense voxel image.
|
||||
//! - [`resolution`]: closed-form range/cross-range resolution and
|
||||
//! coherence-budget formulas (`ΔR = c/2B`, `δ_CR ≈ λR/2L`, the `λ/8`
|
||||
//! antenna-pose tolerance), checked against the reconstruction's actual
|
||||
//! behavior in `tests/physics_validation.rs` rather than merely asserted.
|
||||
//! - [`geometry`]: antenna poses and synthetic-aperture trajectories.
|
||||
//!
|
||||
//! ## What this crate is NOT (ADR-283 §1, honesty boundary)
|
||||
//!
|
||||
//! - **Not a hardware driver.** There is no VNA, SDR, or wideband RF
|
||||
//! front-end integration here, and none of `wifi-densepose-hardware`'s
|
||||
//! ESP32/CSI chipset code is touched. Every measurement in this crate's
|
||||
//! tests and benchmarks is [`measurement::simulate_measurement`] output.
|
||||
//! - **Not a reproduction of any published system.** ADR-278 already gates
|
||||
//! RISE/DiffRadar/GeRaF reproduction as a separate, much larger research
|
||||
//! program with its own acceptance gates; this crate does not attempt
|
||||
//! any of them. It is scoped one level below that: the bare
|
||||
//! measurement-model + backprojection primitive those (or any other SAR
|
||||
//! pipeline) would be built on.
|
||||
//! - **Not a claim about Applied Electrodynamics' product.** Their exact
|
||||
//! waveform, antenna count, bandwidth, and reconstruction algorithm are
|
||||
//! undisclosed; nothing here reproduces or benchmarks against their
|
||||
//! device. It is simply the same well-established SAR/GPR physics
|
||||
//! (Skolnik-style backprojection) applied to synthetic data.
|
||||
//! - **Not wired into `ruview-unified` or `GaussianMap`.** ADR-278 §2.4
|
||||
//! names that as the eventual integration point once (and if) a gated
|
||||
//! reconstruction system exists; this crate is deliberately a leaf with
|
||||
//! no RuView dependency (the nvsim pattern) so its physics can be
|
||||
//! validated in isolation first.
|
||||
//! - **Every number produced by this crate is SYNTHETIC / evidence level
|
||||
//! L0** (ADR-282's mandatory evidence ladder), generated by its own
|
||||
//! forward simulator, scored against its own ground truth. It says
|
||||
//! nothing about real-world through-wall imaging performance, which
|
||||
//! depends on multipath, wall materials, antenna gain patterns, receiver
|
||||
//! noise figures, and calibration accuracy that this crate does not
|
||||
//! model.
|
||||
//!
|
||||
//! ## Design commitments
|
||||
//!
|
||||
//! - **Deterministic**: [`measurement::simulate_measurement`] is seeded
|
||||
//! ChaCha20 (the nvsim / ruview-unified commitment) -- same seed yields
|
||||
//! byte-identical output on every machine.
|
||||
//! - **Proven, not asserted**: `tests/physics_validation.rs` checks the
|
||||
//! reconstruction's actual point-target localization error, range
|
||||
//! resolution, cross-range resolution, and pose-error sensitivity against
|
||||
//! the closed-form predictions in [`resolution`] -- the same discipline
|
||||
//! `ruview-unified` uses for its ray tracer (checked against Friis and
|
||||
//! reciprocity) and its encoder (checked against finite differences).
|
||||
|
||||
#![warn(missing_docs)]
|
||||
#![forbid(unsafe_code)]
|
||||
|
||||
pub mod geometry;
|
||||
pub mod measurement;
|
||||
pub mod pointcloud;
|
||||
pub mod reconstruct;
|
||||
pub mod resolution;
|
||||
|
||||
pub use geometry::{linear_aperture, AntennaPose, Point3};
|
||||
pub use measurement::{simulate_measurement, FrequencySweep, Measurement, ScatteringTarget};
|
||||
pub use pointcloud::{extract_point_cloud, PointCloudPoint};
|
||||
pub use reconstruct::{backproject, focus_at_point, ReflectivityImage, VoxelGrid};
|
||||
@@ -0,0 +1,227 @@
|
||||
//! Forward measurement model: simulate the complex, stepped-frequency
|
||||
//! returns a monostatic synthetic-aperture radar would record from a set
|
||||
//! of point scatterers (ADR-283 §2).
|
||||
//!
|
||||
//! ```text
|
||||
//! y_{m,k} = sum_j sigma_j / R_{m,j}^2 * exp(-i * 4*pi * f_k * R_{m,j} / c) + noise
|
||||
//! ```
|
||||
//!
|
||||
//! where `m` indexes antenna position, `k` indexes swept frequency, `j`
|
||||
//! indexes point scatterer, `R_{m,j}` is the range from antenna position
|
||||
//! `m` to scatterer `j`, and the `4*pi*f/c` phase term is the two-way
|
||||
//! (round-trip) propagation phase for a monostatic (single antenna acting
|
||||
//! as both transmitter and receiver) system. The `1/R^2` term is the
|
||||
//! two-way free-space spreading loss (amplitude, not power).
|
||||
//!
|
||||
//! This is a deliberately simplified physical model: free-space
|
||||
//! propagation only (no multipath, no per-material attenuation, no
|
||||
//! antenna gain pattern). It exists to give the reconstruction algorithm
|
||||
//! in [`crate::reconstruct`] a known ground truth to be checked against,
|
||||
//! not to predict real hardware performance.
|
||||
|
||||
use crate::geometry::{AntennaPose, Point3};
|
||||
use num_complex::Complex64;
|
||||
use rand::{Rng, SeedableRng};
|
||||
use rand_chacha::ChaCha20Rng;
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::f64::consts::PI;
|
||||
|
||||
use crate::resolution::SPEED_OF_LIGHT_M_PER_S;
|
||||
|
||||
/// A stepped-frequency sweep: `n_steps` evenly spaced frequencies from
|
||||
/// `start_hz` to `stop_hz` inclusive.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||||
pub struct FrequencySweep {
|
||||
/// Lowest swept frequency, Hz.
|
||||
pub start_hz: f64,
|
||||
/// Highest swept frequency, Hz.
|
||||
pub stop_hz: f64,
|
||||
/// Number of evenly spaced frequency steps (>= 1).
|
||||
pub n_steps: usize,
|
||||
}
|
||||
|
||||
impl FrequencySweep {
|
||||
/// Construct a sweep. Panics if `n_steps == 0` or `stop_hz < start_hz`.
|
||||
pub fn new(start_hz: f64, stop_hz: f64, n_steps: usize) -> Self {
|
||||
assert!(n_steps > 0, "a frequency sweep needs at least one step");
|
||||
assert!(stop_hz >= start_hz, "stop_hz must be >= start_hz");
|
||||
Self { start_hz, stop_hz, n_steps }
|
||||
}
|
||||
|
||||
/// The `n_steps` evenly spaced frequencies, Hz, ascending.
|
||||
pub fn frequencies(&self) -> Vec<f64> {
|
||||
if self.n_steps == 1 {
|
||||
return vec![self.start_hz];
|
||||
}
|
||||
(0..self.n_steps)
|
||||
.map(|i| {
|
||||
let t = i as f64 / (self.n_steps - 1) as f64;
|
||||
self.start_hz + (self.stop_hz - self.start_hz) * t
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Total swept bandwidth, Hz.
|
||||
pub fn bandwidth_hz(&self) -> f64 {
|
||||
self.stop_hz - self.start_hz
|
||||
}
|
||||
|
||||
/// Center frequency, Hz.
|
||||
pub fn center_freq_hz(&self) -> f64 {
|
||||
(self.start_hz + self.stop_hz) / 2.0
|
||||
}
|
||||
}
|
||||
|
||||
/// A single point scatterer: a location and a scalar reflectivity
|
||||
/// (dimensionless; only relative magnitudes across scatterers matter).
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||||
pub struct ScatteringTarget {
|
||||
/// Scatterer location, meters.
|
||||
pub position: Point3,
|
||||
/// Scalar reflectivity (>= 0 for a physical scatterer; the forward
|
||||
/// model does not enforce this so synthetic "negative" scatterers can
|
||||
/// be used to probe reconstruction linearity in tests).
|
||||
pub reflectivity: f64,
|
||||
}
|
||||
|
||||
impl ScatteringTarget {
|
||||
/// Construct a target.
|
||||
pub fn new(position: Point3, reflectivity: f64) -> Self {
|
||||
Self { position, reflectivity }
|
||||
}
|
||||
}
|
||||
|
||||
/// The complex, stepped-frequency measurement recorded at every
|
||||
/// (antenna position, frequency) pair. Row-major: index `[m][k]` is at
|
||||
/// `samples[m * n_freqs + k]`.
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct Measurement {
|
||||
/// Number of antenna positions.
|
||||
pub n_poses: usize,
|
||||
/// Number of frequency steps.
|
||||
pub n_freqs: usize,
|
||||
/// Complex samples, row-major over (pose, frequency).
|
||||
pub samples: Vec<Complex64>,
|
||||
}
|
||||
|
||||
impl Measurement {
|
||||
/// The complex sample at antenna position `m`, frequency index `k`.
|
||||
pub fn get(&self, m: usize, k: usize) -> Complex64 {
|
||||
self.samples[m * self.n_freqs + k]
|
||||
}
|
||||
}
|
||||
|
||||
/// Simulate the forward measurement model for `targets` observed from
|
||||
/// `poses` across the frequencies in `sweep`, with i.i.d. complex Gaussian
|
||||
/// noise of standard deviation `noise_std` (per real/imaginary component)
|
||||
/// added to every sample. Deterministic given `seed` (ChaCha20, the
|
||||
/// nvsim/ruview-unified reproducibility commitment: same seed -> identical
|
||||
/// output on every machine).
|
||||
pub fn simulate_measurement(
|
||||
poses: &[AntennaPose],
|
||||
sweep: &FrequencySweep,
|
||||
targets: &[ScatteringTarget],
|
||||
noise_std: f64,
|
||||
seed: u64,
|
||||
) -> Measurement {
|
||||
let freqs = sweep.frequencies();
|
||||
let n_poses = poses.len();
|
||||
let n_freqs = freqs.len();
|
||||
let mut rng = ChaCha20Rng::seed_from_u64(seed);
|
||||
let mut samples = Vec::with_capacity(n_poses * n_freqs);
|
||||
|
||||
for pose in poses {
|
||||
for &f in &freqs {
|
||||
let mut y = Complex64::new(0.0, 0.0);
|
||||
for t in targets {
|
||||
let r = pose.position.distance(&t.position);
|
||||
if r < 1e-6 {
|
||||
// Degenerate: antenna co-located with target. Skip to
|
||||
// avoid a divide-by-zero singularity in 1/R^2.
|
||||
continue;
|
||||
}
|
||||
let amplitude = t.reflectivity / (r * r);
|
||||
let phase = -4.0 * PI * f * r / SPEED_OF_LIGHT_M_PER_S;
|
||||
y += Complex64::from_polar(amplitude, phase);
|
||||
}
|
||||
if noise_std > 0.0 {
|
||||
y += Complex64::new(gaussian(&mut rng, noise_std), gaussian(&mut rng, noise_std));
|
||||
}
|
||||
samples.push(y);
|
||||
}
|
||||
}
|
||||
|
||||
Measurement { n_poses, n_freqs, samples }
|
||||
}
|
||||
|
||||
/// Box-Muller transform: one N(0, std^2) sample from two uniform draws.
|
||||
fn gaussian(rng: &mut ChaCha20Rng, std: f64) -> f64 {
|
||||
let u1: f64 = rng.gen_range(f64::EPSILON..1.0);
|
||||
let u2: f64 = rng.gen_range(0.0..1.0);
|
||||
std * (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::geometry::linear_aperture;
|
||||
|
||||
#[test]
|
||||
fn frequency_sweep_endpoints_and_count() {
|
||||
let sweep = FrequencySweep::new(1.0e9, 4.0e9, 5);
|
||||
let freqs = sweep.frequencies();
|
||||
assert_eq!(freqs.len(), 5);
|
||||
assert!((freqs[0] - 1.0e9).abs() < 1e-6);
|
||||
assert!((freqs[4] - 4.0e9).abs() < 1e-6);
|
||||
assert!((sweep.bandwidth_hz() - 3.0e9).abs() < 1e-6);
|
||||
assert!((sweep.center_freq_hz() - 2.5e9).abs() < 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn single_step_sweep_is_just_start_hz() {
|
||||
let sweep = FrequencySweep::new(2.0e9, 2.0e9, 1);
|
||||
assert_eq!(sweep.frequencies(), vec![2.0e9]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn simulate_measurement_is_deterministic_given_seed() {
|
||||
let poses = linear_aperture(Point3::new(0.0, 0.0, 0.0), Point3::new(0.3, 0.0, 0.0), 4);
|
||||
let sweep = FrequencySweep::new(2.0e9, 3.0e9, 8);
|
||||
let targets = vec![ScatteringTarget::new(Point3::new(0.15, 1.0, 0.0), 1.0)];
|
||||
let a = simulate_measurement(&poses, &sweep, &targets, 0.01, 42);
|
||||
let b = simulate_measurement(&poses, &sweep, &targets, 0.01, 42);
|
||||
for (sa, sb) in a.samples.iter().zip(b.samples.iter()) {
|
||||
assert_eq!(sa, sb, "same seed must give byte-identical measurements");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn different_seeds_give_different_noise() {
|
||||
let poses = linear_aperture(Point3::new(0.0, 0.0, 0.0), Point3::new(0.3, 0.0, 0.0), 4);
|
||||
let sweep = FrequencySweep::new(2.0e9, 3.0e9, 8);
|
||||
let targets = vec![ScatteringTarget::new(Point3::new(0.15, 1.0, 0.0), 1.0)];
|
||||
let a = simulate_measurement(&poses, &sweep, &targets, 0.05, 1);
|
||||
let b = simulate_measurement(&poses, &sweep, &targets, 0.05, 2);
|
||||
assert_ne!(a.samples, b.samples);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn zero_noise_is_purely_deterministic_physics() {
|
||||
let poses = linear_aperture(Point3::new(0.0, 0.0, 0.0), Point3::new(0.3, 0.0, 0.0), 2);
|
||||
let sweep = FrequencySweep::new(2.0e9, 2.0e9, 1);
|
||||
let targets = vec![ScatteringTarget::new(Point3::new(0.0, 1.0, 0.0), 2.0)];
|
||||
let m = simulate_measurement(&poses, &sweep, &targets, 0.0, 7);
|
||||
// Antenna 0 is directly below the target at range 1.0 m.
|
||||
let r = 1.0_f64;
|
||||
let expected_amp = 2.0 / (r * r);
|
||||
assert!((m.get(0, 0).norm() - expected_amp).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn no_noise_std_gt_zero_produces_nonzero_noise() {
|
||||
let poses = linear_aperture(Point3::new(0.0, 0.0, 0.0), Point3::new(0.3, 0.0, 0.0), 1);
|
||||
let sweep = FrequencySweep::new(2.0e9, 2.0e9, 1);
|
||||
let m = simulate_measurement(&poses, &sweep, &[], 1.0, 5);
|
||||
assert!(m.get(0, 0).norm() > 0.0, "no targets but noise_std>0 must still yield noise");
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,116 @@
|
||||
//! Extract a sparse point cloud from a dense [`ReflectivityImage`].
|
||||
//!
|
||||
//! A `nx * ny * nz` voxel grid is not a useful end product on its own --
|
||||
//! real point-cloud consumers (visualization, `ruview-unified`'s
|
||||
//! `GaussianMap`, downstream fusion) want a short list of "here is
|
||||
//! something" points, not every voxel. This module does simple
|
||||
//! threshold + local-maximum extraction: no clustering, no material
|
||||
//! classification, no confidence calibration against real data (ADR-283
|
||||
//! §5 -- explicitly out of scope for this crate).
|
||||
|
||||
use crate::geometry::Point3;
|
||||
use crate::reconstruct::ReflectivityImage;
|
||||
use serde::{Deserialize, Serialize};
|
||||
|
||||
/// A single detected point: a location and its reconstructed reflectivity
|
||||
/// magnitude (relative, not calibrated to any physical unit).
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||||
pub struct PointCloudPoint {
|
||||
/// World-space location, meters.
|
||||
pub position: Point3,
|
||||
/// Reconstructed reflectivity magnitude at this voxel.
|
||||
pub magnitude: f64,
|
||||
}
|
||||
|
||||
/// Extract detected points from `image`: every voxel whose magnitude is
|
||||
/// (a) at least `threshold_fraction` of the image's peak magnitude, and
|
||||
/// (b) a local maximum among its 6-connected neighbors (so a single broad
|
||||
/// blob yields one point, not every voxel inside it).
|
||||
///
|
||||
/// `threshold_fraction` must be in `(0.0, 1.0]`. A typical value is
|
||||
/// `0.5` (a classic radar/SAR "half-power point" style threshold).
|
||||
pub fn extract_point_cloud(image: &ReflectivityImage, threshold_fraction: f64) -> Vec<PointCloudPoint> {
|
||||
assert!(
|
||||
threshold_fraction > 0.0 && threshold_fraction <= 1.0,
|
||||
"threshold_fraction must be in (0, 1]"
|
||||
);
|
||||
let grid = &image.grid;
|
||||
let peak = image.magnitude.iter().cloned().fold(0.0_f64, f64::max);
|
||||
if peak <= 0.0 {
|
||||
return Vec::new();
|
||||
}
|
||||
let threshold = peak * threshold_fraction;
|
||||
|
||||
let mag_at = |i: i64, j: i64, k: i64| -> f64 {
|
||||
if i < 0 || j < 0 || k < 0 || i as usize >= grid.nx || j as usize >= grid.ny || k as usize >= grid.nz {
|
||||
return f64::NEG_INFINITY;
|
||||
}
|
||||
let linear = grid.linear_index(i as usize, j as usize, k as usize);
|
||||
image.magnitude[linear]
|
||||
};
|
||||
|
||||
let mut points = Vec::new();
|
||||
for k in 0..grid.nz {
|
||||
for j in 0..grid.ny {
|
||||
for i in 0..grid.nx {
|
||||
let here = mag_at(i as i64, j as i64, k as i64);
|
||||
if here < threshold {
|
||||
continue;
|
||||
}
|
||||
let neighbors = [
|
||||
mag_at(i as i64 - 1, j as i64, k as i64),
|
||||
mag_at(i as i64 + 1, j as i64, k as i64),
|
||||
mag_at(i as i64, j as i64 - 1, k as i64),
|
||||
mag_at(i as i64, j as i64 + 1, k as i64),
|
||||
mag_at(i as i64, j as i64, k as i64 - 1),
|
||||
mag_at(i as i64, j as i64, k as i64 + 1),
|
||||
];
|
||||
if neighbors.iter().all(|&n| here >= n) {
|
||||
points.push(PointCloudPoint {
|
||||
position: grid.voxel_center(i, j, k),
|
||||
magnitude: here,
|
||||
});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
points
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::geometry::linear_aperture;
|
||||
use crate::measurement::{simulate_measurement, FrequencySweep, ScatteringTarget};
|
||||
use crate::reconstruct::{backproject, VoxelGrid};
|
||||
|
||||
#[test]
|
||||
fn single_target_yields_a_single_detected_point() {
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
let measurement = simulate_measurement(&poses, &sweep, &[target], 0.0, 3);
|
||||
let grid = VoxelGrid::new(Point3::new(-0.5, 1.6, -0.5), 0.05, 21, 17, 21);
|
||||
let image = backproject(&measurement, &poses, &sweep, &grid);
|
||||
|
||||
let points = extract_point_cloud(&image, 0.5);
|
||||
assert!(!points.is_empty(), "must detect at least the true target");
|
||||
let best = points.iter().max_by(|a, b| a.magnitude.partial_cmp(&b.magnitude).unwrap()).unwrap();
|
||||
assert!(best.position.distance(&target.position) < 0.1);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn empty_image_yields_no_points() {
|
||||
let grid = VoxelGrid::new(Point3::new(0.0, 0.0, 0.0), 0.1, 3, 3, 3);
|
||||
let image = ReflectivityImage { grid, magnitude: vec![0.0; grid.len()] };
|
||||
assert!(extract_point_cloud(&image, 0.5).is_empty());
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[should_panic(expected = "threshold_fraction")]
|
||||
fn rejects_out_of_range_threshold() {
|
||||
let grid = VoxelGrid::new(Point3::new(0.0, 0.0, 0.0), 0.1, 2, 2, 2);
|
||||
let image = ReflectivityImage { grid, magnitude: vec![1.0; grid.len()] };
|
||||
let _ = extract_point_cloud(&image, 1.5);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,214 @@
|
||||
//! Delay-and-sum backprojection reconstruction (ADR-283 §2).
|
||||
//!
|
||||
//! Given a [`crate::measurement::Measurement`] recorded from known antenna
|
||||
//! [`AntennaPose`]s across a known [`FrequencySweep`], reconstruct a 3D
|
||||
//! reflectivity image on a regular voxel grid:
|
||||
//!
|
||||
//! ```text
|
||||
//! I(x) = | (1 / (M*K)) * sum_m sum_k y_{m,k} * R_{m,x}^2 * exp(+i * 4*pi * f_k * R_{m,x} / c) |
|
||||
//! ```
|
||||
//!
|
||||
//! This is the matched-filter / frequency-domain backprojection kernel:
|
||||
//! for a *correct* hypothesis voxel `x` coinciding with a real scatterer,
|
||||
//! every (pose, frequency) term's phase-correction exactly cancels the
|
||||
//! phase the forward model applied in [`crate::measurement`], so the sum
|
||||
//! coheres constructively. For any other voxel the per-term phases are
|
||||
//! effectively uncorrelated across the (pose, frequency) grid and the sum
|
||||
//! averages toward zero. The `R_{m,x}^2` factor undoes the forward
|
||||
//! model's `1/R^2` spreading-loss term (matched-filter gain
|
||||
//! compensation), so voxel brightness reflects relative reflectivity
|
||||
//! rather than falling off with range.
|
||||
|
||||
use crate::geometry::{AntennaPose, Point3};
|
||||
use crate::measurement::{FrequencySweep, Measurement};
|
||||
use num_complex::Complex64;
|
||||
use rayon::prelude::*;
|
||||
use std::f64::consts::PI;
|
||||
|
||||
use crate::resolution::SPEED_OF_LIGHT_M_PER_S;
|
||||
|
||||
/// A regular 3D grid of voxel centers over an axis-aligned box.
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub struct VoxelGrid {
|
||||
/// Grid origin (the center of voxel `(0,0,0)`), meters.
|
||||
pub origin: Point3,
|
||||
/// Voxel edge length along each axis, meters.
|
||||
pub spacing: f64,
|
||||
/// Number of voxels along x.
|
||||
pub nx: usize,
|
||||
/// Number of voxels along y.
|
||||
pub ny: usize,
|
||||
/// Number of voxels along z.
|
||||
pub nz: usize,
|
||||
}
|
||||
|
||||
impl VoxelGrid {
|
||||
/// Construct a grid.
|
||||
pub fn new(origin: Point3, spacing: f64, nx: usize, ny: usize, nz: usize) -> Self {
|
||||
assert!(spacing > 0.0, "voxel spacing must be positive");
|
||||
Self { origin, spacing, nx, ny, nz }
|
||||
}
|
||||
|
||||
/// Total voxel count.
|
||||
pub fn len(&self) -> usize {
|
||||
self.nx * self.ny * self.nz
|
||||
}
|
||||
|
||||
/// True if the grid has zero voxels along any axis.
|
||||
pub fn is_empty(&self) -> bool {
|
||||
self.len() == 0
|
||||
}
|
||||
|
||||
/// World-space center of voxel `(i, j, k)`.
|
||||
pub fn voxel_center(&self, i: usize, j: usize, k: usize) -> Point3 {
|
||||
Point3::new(
|
||||
self.origin.x + (i as f64) * self.spacing,
|
||||
self.origin.y + (j as f64) * self.spacing,
|
||||
self.origin.z + (k as f64) * self.spacing,
|
||||
)
|
||||
}
|
||||
|
||||
/// Flatten a 3D voxel index into a linear index (row-major, x fastest).
|
||||
pub fn linear_index(&self, i: usize, j: usize, k: usize) -> usize {
|
||||
(k * self.ny + j) * self.nx + i
|
||||
}
|
||||
|
||||
/// Recover the 3D voxel index `(i, j, k)` from a linear index.
|
||||
pub fn unflatten(&self, linear: usize) -> (usize, usize, usize) {
|
||||
let i = linear % self.nx;
|
||||
let j = (linear / self.nx) % self.ny;
|
||||
let k = linear / (self.nx * self.ny);
|
||||
(i, j, k)
|
||||
}
|
||||
}
|
||||
|
||||
/// The reconstructed reflectivity image: one magnitude value per voxel,
|
||||
/// row-major (`grid.linear_index`/`unflatten` order).
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct ReflectivityImage {
|
||||
/// The grid this image was reconstructed on.
|
||||
pub grid: VoxelGrid,
|
||||
/// Per-voxel reflectivity magnitude, same length as `grid.len()`.
|
||||
pub magnitude: Vec<f64>,
|
||||
}
|
||||
|
||||
impl ReflectivityImage {
|
||||
/// The voxel with the largest magnitude, and its world-space center.
|
||||
pub fn peak(&self) -> (Point3, f64) {
|
||||
let (idx, &mag) = self
|
||||
.magnitude
|
||||
.iter()
|
||||
.enumerate()
|
||||
.max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
|
||||
.expect("grid must have at least one voxel");
|
||||
let (i, j, k) = self.grid.unflatten(idx);
|
||||
(self.grid.voxel_center(i, j, k), mag)
|
||||
}
|
||||
}
|
||||
|
||||
/// Reconstruct a reflectivity image from `measurement`, recorded at
|
||||
/// `poses` across `sweep`, onto `grid`. Parallelized over voxels (rayon).
|
||||
///
|
||||
/// `poses` and `sweep` must describe the *same* geometry the measurement
|
||||
/// was recorded with (or, when studying pose-error sensitivity, a
|
||||
/// deliberately perturbed version of it -- see
|
||||
/// `tests/physics_validation.rs`).
|
||||
pub fn backproject(
|
||||
measurement: &Measurement,
|
||||
poses: &[AntennaPose],
|
||||
sweep: &FrequencySweep,
|
||||
grid: &VoxelGrid,
|
||||
) -> ReflectivityImage {
|
||||
assert_eq!(poses.len(), measurement.n_poses, "pose count must match measurement");
|
||||
let freqs = sweep.frequencies();
|
||||
assert_eq!(freqs.len(), measurement.n_freqs, "frequency count must match measurement");
|
||||
|
||||
let n = grid.len();
|
||||
|
||||
let magnitude: Vec<f64> = (0..n)
|
||||
.into_par_iter()
|
||||
.map(|linear| {
|
||||
let (i, j, k) = grid.unflatten(linear);
|
||||
let voxel = grid.voxel_center(i, j, k);
|
||||
focus_at_point(measurement, poses, &freqs, &voxel)
|
||||
})
|
||||
.collect();
|
||||
|
||||
ReflectivityImage { grid: *grid, magnitude }
|
||||
}
|
||||
|
||||
/// Evaluate the coherent backprojection sum at a single world-space
|
||||
/// `point`, without building a grid. This is the same matched-filter
|
||||
/// kernel [`backproject`] evaluates per voxel; exposed directly so callers
|
||||
/// (and tests) can measure focus quality exactly at a location of
|
||||
/// interest -- e.g. a known target position -- rather than only at
|
||||
/// whatever grid points happen to be sampled.
|
||||
pub fn focus_at_point(measurement: &Measurement, poses: &[AntennaPose], freqs: &[f64], point: &Point3) -> f64 {
|
||||
assert_eq!(poses.len(), measurement.n_poses, "pose count must match measurement");
|
||||
assert_eq!(freqs.len(), measurement.n_freqs, "frequency count must match measurement");
|
||||
|
||||
let n_terms = (measurement.n_poses * measurement.n_freqs) as f64;
|
||||
let mut acc = Complex64::new(0.0, 0.0);
|
||||
for (m, pose) in poses.iter().enumerate() {
|
||||
let r = pose.position.distance(point);
|
||||
if r < 1e-6 {
|
||||
continue;
|
||||
}
|
||||
let gain_compensation = r * r;
|
||||
for (kf, &f) in freqs.iter().enumerate() {
|
||||
let phase = 4.0 * PI * f * r / SPEED_OF_LIGHT_M_PER_S;
|
||||
acc += measurement.get(m, kf) * gain_compensation * Complex64::from_polar(1.0, phase);
|
||||
}
|
||||
}
|
||||
acc.norm() / n_terms
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::geometry::linear_aperture;
|
||||
use crate::measurement::{simulate_measurement, ScatteringTarget};
|
||||
|
||||
#[test]
|
||||
fn voxel_grid_flatten_unflatten_roundtrip() {
|
||||
let grid = VoxelGrid::new(Point3::new(0.0, 0.0, 0.0), 0.1, 4, 5, 3);
|
||||
for k in 0..grid.nz {
|
||||
for j in 0..grid.ny {
|
||||
for i in 0..grid.nx {
|
||||
let lin = grid.linear_index(i, j, k);
|
||||
assert_eq!(grid.unflatten(lin), (i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn single_point_target_reconstructs_at_its_true_location() {
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
let measurement = simulate_measurement(&poses, &sweep, &[target], 0.0, 1);
|
||||
|
||||
let grid = VoxelGrid::new(Point3::new(-0.5, 1.6, -0.5), 0.05, 21, 17, 21);
|
||||
let image = backproject(&measurement, &poses, &sweep, &grid);
|
||||
let (peak_loc, _peak_mag) = image.peak();
|
||||
|
||||
let err = peak_loc.distance(&target.position);
|
||||
assert!(err < 0.1, "reconstructed peak {:?} should be within one voxel-ish of the true target {:?}, err={err}", peak_loc, target.position);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn peak_at_target_is_far_above_background() {
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 32);
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
let measurement = simulate_measurement(&poses, &sweep, &[target], 0.0, 2);
|
||||
|
||||
let grid = VoxelGrid::new(Point3::new(-0.5, 1.6, -0.5), 0.05, 21, 17, 21);
|
||||
let image = backproject(&measurement, &poses, &sweep, &grid);
|
||||
let (_peak_loc, peak_mag) = image.peak();
|
||||
let mean_mag: f64 = image.magnitude.iter().sum::<f64>() / image.magnitude.len() as f64;
|
||||
|
||||
assert!(peak_mag > mean_mag * 5.0, "coherent focus at the target must dominate the incoherent background: peak={peak_mag}, mean={mean_mag}");
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,105 @@
|
||||
//! Closed-form resolution and coherence-budget formulas (ADR-283 §3).
|
||||
//!
|
||||
//! These are textbook radar-imaging identities (see e.g. Skolnik, *Radar
|
||||
//! Handbook*, and the standard stripmap-SAR cross-range formula). They are
|
||||
//! implemented here so the crate's own reconstruction behavior can be
|
||||
//! checked against them in [`tests/physics_validation.rs`] rather than
|
||||
//! merely asserted in documentation.
|
||||
|
||||
/// Speed of light in vacuum, m/s.
|
||||
pub const SPEED_OF_LIGHT_M_PER_S: f64 = 299_792_458.0;
|
||||
|
||||
/// Wavelength (meters) of a signal at `freq_hz`.
|
||||
pub fn wavelength_m(freq_hz: f64) -> f64 {
|
||||
SPEED_OF_LIGHT_M_PER_S / freq_hz
|
||||
}
|
||||
|
||||
/// Range resolution (meters) of a stepped-frequency / wideband radar with
|
||||
/// total swept bandwidth `bandwidth_hz`: `ΔR = c / (2B)`.
|
||||
///
|
||||
/// This is the Rayleigh-style minimum range separation at which two
|
||||
/// point targets on the same bearing become distinguishable after pulse
|
||||
/// compression / coherent range processing. It does **not** depend on
|
||||
/// carrier frequency, antenna count, or synthetic-aperture length --
|
||||
/// only on how much spectrum was actually swept.
|
||||
pub fn range_resolution_m(bandwidth_hz: f64) -> f64 {
|
||||
SPEED_OF_LIGHT_M_PER_S / (2.0 * bandwidth_hz)
|
||||
}
|
||||
|
||||
/// Cross-range (azimuth) resolution (meters) of a synthetic aperture of
|
||||
/// physical length `aperture_length_m`, imaging a target at `range_m`,
|
||||
/// at carrier frequency `center_freq_hz`: `δ_CR ≈ λ·R / (2·L)`.
|
||||
///
|
||||
/// This is the classic stripmap-SAR angular-resolution identity: doubling
|
||||
/// the aperture (or halving the wavelength) halves the achievable
|
||||
/// cross-range spot size at a fixed range. It is undefined (returns
|
||||
/// `f64::INFINITY`) for a degenerate (zero-length) aperture -- a single
|
||||
/// antenna position carries no cross-range information at all, which is
|
||||
/// exactly the point of building a synthetic aperture in the first place.
|
||||
pub fn cross_range_resolution_m(center_freq_hz: f64, aperture_length_m: f64, range_m: f64) -> f64 {
|
||||
if aperture_length_m <= 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
wavelength_m(center_freq_hz) * range_m / (2.0 * aperture_length_m)
|
||||
}
|
||||
|
||||
/// Maximum antenna-position error (meters) that keeps a coherent
|
||||
/// (phase-focused) reconstruction inside the classical quarter-wave
|
||||
/// budget, at carrier frequency `center_freq_hz`.
|
||||
///
|
||||
/// Derivation: moving an antenna's phase center by `Δp` while looking
|
||||
/// (worst case) directly along boresight at the target changes the
|
||||
/// round-trip path length by up to `2·Δp` (both the outbound and return
|
||||
/// leg shift by `Δp`). Keeping that two-way path error under the
|
||||
/// standard quarter-wavelength coherence budget (`λ/4` -- the same
|
||||
/// criterion used for reflector-antenna and optical-surface tolerancing)
|
||||
/// requires `2·Δp ≤ λ/4`, i.e. `Δp ≤ λ/8`.
|
||||
///
|
||||
/// Position error beyond this does not make reconstruction impossible --
|
||||
/// it degrades the coherent sum smoothly (see
|
||||
/// `phase_error_degrades_focus_beyond_pose_budget` in
|
||||
/// `tests/physics_validation.rs`) -- but it is the standard rule-of-thumb
|
||||
/// budget for "still well focused."
|
||||
pub fn max_coherent_pose_error_m(center_freq_hz: f64) -> f64 {
|
||||
wavelength_m(center_freq_hz) / 8.0
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn range_resolution_matches_hand_computed_value() {
|
||||
// 3 GHz of swept bandwidth: c/(2*3e9) = 4.9965...cm.
|
||||
let r = range_resolution_m(3.0e9);
|
||||
assert!((r - 0.049_965_409_666_666_66).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn eight_ghz_pose_budget_is_about_5mm() {
|
||||
// At 8 GHz, lambda = c/f ~= 37.47mm, so lambda/8 ~= 4.68mm -- close
|
||||
// to the ~5mm rule-of-thumb quoted in the motivating design note.
|
||||
let budget = max_coherent_pose_error_m(8.0e9);
|
||||
assert!((budget - wavelength_m(8.0e9) / 8.0).abs() < 1e-12);
|
||||
assert!(budget < 0.005 && budget > 0.004);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cross_range_resolution_improves_with_longer_aperture() {
|
||||
let short = cross_range_resolution_m(5.0e9, 0.1, 2.0);
|
||||
let long = cross_range_resolution_m(5.0e9, 1.0, 2.0);
|
||||
assert!(long < short, "10x longer aperture must give finer cross-range resolution");
|
||||
// Exactly linear in 1/L.
|
||||
assert!((short / long - 10.0).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn zero_length_aperture_has_no_cross_range_resolution() {
|
||||
assert_eq!(cross_range_resolution_m(5.0e9, 0.0, 2.0), f64::INFINITY);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wider_bandwidth_gives_finer_range_resolution() {
|
||||
assert!(range_resolution_m(4.0e9) < range_resolution_m(1.0e9));
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,216 @@
|
||||
//! Checks the reconstruction's *actual* behavior against the closed-form
|
||||
//! predictions in `wifi_densepose_sar::resolution`, rather than merely
|
||||
//! asserting the formulas in documentation (ADR-283 §3, the
|
||||
//! ruview-unified "proven, not asserted" discipline).
|
||||
//!
|
||||
//! Three physical claims are validated end-to-end (forward-simulate ->
|
||||
//! backproject -> measure the reconstruction's behavior):
|
||||
//!
|
||||
//! 1. Two point targets separated along *range* resolve into two distinct
|
||||
//! peaks only once their separation exceeds `range_resolution_m`.
|
||||
//! 2. Two point targets separated along *cross-range* (same range,
|
||||
//! different bearing) resolve only once the synthetic-aperture length
|
||||
//! is long enough per `cross_range_resolution_m` -- a single short
|
||||
//! aperture cannot resolve them no matter how much bandwidth is used.
|
||||
//! 3. Antenna-position error decoheres the reconstruction: coherent focus
|
||||
//! at the true target location trends downward as position error
|
||||
//! grows past the `max_coherent_pose_error_m` (`λ/8`) budget, and has
|
||||
//! collapsed toward the incoherent background by an order of
|
||||
//! magnitude beyond it.
|
||||
|
||||
use wifi_densepose_sar::geometry::linear_aperture;
|
||||
use wifi_densepose_sar::measurement::{simulate_measurement, FrequencySweep, ScatteringTarget};
|
||||
use wifi_densepose_sar::reconstruct::{backproject, focus_at_point, VoxelGrid};
|
||||
use wifi_densepose_sar::resolution::{
|
||||
cross_range_resolution_m, max_coherent_pose_error_m, range_resolution_m, wavelength_m,
|
||||
};
|
||||
use wifi_densepose_sar::{AntennaPose, Point3};
|
||||
|
||||
/// Count local maxima at or above `threshold_fraction` of the profile's
|
||||
/// peak, in a 1D magnitude profile. Adjacent samples above threshold count
|
||||
/// as one maximum (a plateau/peak region), not one-per-sample.
|
||||
fn count_resolved_peaks(profile: &[f64], threshold_fraction: f64) -> usize {
|
||||
let peak = profile.iter().cloned().fold(0.0_f64, f64::max);
|
||||
let threshold = peak * threshold_fraction;
|
||||
let mut count = 0;
|
||||
let mut in_peak = false;
|
||||
for &v in profile {
|
||||
if v >= threshold {
|
||||
if !in_peak {
|
||||
count += 1;
|
||||
in_peak = true;
|
||||
}
|
||||
} else {
|
||||
in_peak = false;
|
||||
}
|
||||
}
|
||||
count
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn range_separated_targets_resolve_only_beyond_range_resolution() {
|
||||
let poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let sweep = FrequencySweep::new(2.0e9, 6.0e9, 64); // 4 GHz bandwidth
|
||||
let dr = range_resolution_m(sweep.bandwidth_hz());
|
||||
|
||||
// A 1D range profile: fixed cross-range (x=0, z=0), fine steps in y.
|
||||
let profile_grid = |center_y: f64, half_span: f64| {
|
||||
VoxelGrid::new(Point3::new(0.0, center_y - half_span, 0.0), dr / 6.0, 1, (2.0 * half_span / (dr / 6.0)) as usize, 1)
|
||||
};
|
||||
|
||||
// Case A: well separated (4x the theoretical resolution) -> two peaks.
|
||||
let sep_resolved = 4.0 * dr;
|
||||
let targets_a = vec![
|
||||
ScatteringTarget::new(Point3::new(0.0, 2.0 - sep_resolved / 2.0, 0.0), 1.0),
|
||||
ScatteringTarget::new(Point3::new(0.0, 2.0 + sep_resolved / 2.0, 0.0), 1.0),
|
||||
];
|
||||
let meas_a = simulate_measurement(&poses, &sweep, &targets_a, 0.0, 10);
|
||||
let grid_a = profile_grid(2.0, sep_resolved * 1.5);
|
||||
let image_a = backproject(&meas_a, &poses, &sweep, &grid_a);
|
||||
let peaks_a = count_resolved_peaks(&image_a.magnitude, 0.7);
|
||||
assert_eq!(
|
||||
peaks_a, 2,
|
||||
"targets separated by 4x the range resolution ({sep_resolved:.4} m vs dr={dr:.4} m) must resolve into 2 peaks, got {peaks_a}"
|
||||
);
|
||||
|
||||
// Case B: too close (0.25x the theoretical resolution) -> one merged peak.
|
||||
let sep_unresolved = 0.25 * dr;
|
||||
let targets_b = vec![
|
||||
ScatteringTarget::new(Point3::new(0.0, 2.0 - sep_unresolved / 2.0, 0.0), 1.0),
|
||||
ScatteringTarget::new(Point3::new(0.0, 2.0 + sep_unresolved / 2.0, 0.0), 1.0),
|
||||
];
|
||||
let meas_b = simulate_measurement(&poses, &sweep, &targets_b, 0.0, 11);
|
||||
let grid_b = profile_grid(2.0, dr * 2.0);
|
||||
let image_b = backproject(&meas_b, &poses, &sweep, &grid_b);
|
||||
let peaks_b = count_resolved_peaks(&image_b.magnitude, 0.7);
|
||||
assert_eq!(
|
||||
peaks_b, 1,
|
||||
"targets separated by only 0.25x the range resolution must merge into 1 peak, got {peaks_b}"
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cross_range_separated_targets_resolve_only_with_long_enough_aperture() {
|
||||
let sweep = FrequencySweep::new(3.0e9, 5.0e9, 32); // center 4 GHz
|
||||
let range_m = 2.0;
|
||||
let lambda = wavelength_m(sweep.center_freq_hz());
|
||||
|
||||
let long_aperture_len = 1.0;
|
||||
let short_aperture_len = 0.05;
|
||||
let res_long = cross_range_resolution_m(sweep.center_freq_hz(), long_aperture_len, range_m);
|
||||
let res_short = cross_range_resolution_m(sweep.center_freq_hz(), short_aperture_len, range_m);
|
||||
assert!(res_long < res_short, "a longer aperture must predict finer cross-range resolution");
|
||||
|
||||
// Pick a separation that is resolvable with the long aperture (well
|
||||
// above its predicted resolution) but not with the short one (well
|
||||
// below its much coarser predicted resolution).
|
||||
let separation = 5.0 * res_long;
|
||||
assert!(separation < res_short, "test setup: separation must sit inside the short aperture's blind spot (lambda={lambda:.4})");
|
||||
|
||||
let targets = vec![
|
||||
ScatteringTarget::new(Point3::new(-separation / 2.0, range_m, 0.0), 1.0),
|
||||
ScatteringTarget::new(Point3::new(separation / 2.0, range_m, 0.0), 1.0),
|
||||
];
|
||||
|
||||
let half_span = separation * 1.5;
|
||||
let cross_range_grid = || VoxelGrid::new(Point3::new(-half_span, range_m, 0.0), separation / 20.0, (2.0 * half_span / (separation / 20.0)) as usize, 1, 1);
|
||||
|
||||
// Long aperture: must resolve into two peaks.
|
||||
let poses_long = linear_aperture(
|
||||
Point3::new(-long_aperture_len / 2.0, 0.0, 0.0),
|
||||
Point3::new(long_aperture_len / 2.0, 0.0, 0.0),
|
||||
41,
|
||||
);
|
||||
let meas_long = simulate_measurement(&poses_long, &sweep, &targets, 0.0, 20);
|
||||
let grid_long = cross_range_grid();
|
||||
let image_long = backproject(&meas_long, &poses_long, &sweep, &grid_long);
|
||||
let peaks_long = count_resolved_peaks(&image_long.magnitude, 0.7);
|
||||
assert_eq!(peaks_long, 2, "a {long_aperture_len} m synthetic aperture must resolve cross-range-separated targets {separation:.4} m apart, got {peaks_long} peak(s)");
|
||||
|
||||
// Short aperture: must NOT resolve (collapses to one blob/ridge).
|
||||
let poses_short = linear_aperture(
|
||||
Point3::new(-short_aperture_len / 2.0, 0.0, 0.0),
|
||||
Point3::new(short_aperture_len / 2.0, 0.0, 0.0),
|
||||
41,
|
||||
);
|
||||
let meas_short = simulate_measurement(&poses_short, &sweep, &targets, 0.0, 21);
|
||||
let grid_short = cross_range_grid();
|
||||
let image_short = backproject(&meas_short, &poses_short, &sweep, &grid_short);
|
||||
let peaks_short = count_resolved_peaks(&image_short.magnitude, 0.7);
|
||||
assert_eq!(peaks_short, 1, "a {short_aperture_len} m synthetic aperture (far below the required {res_short:.4} m cross-range resolution) must NOT resolve the same targets, got {peaks_short} peak(s)");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn phase_error_from_pose_jitter_degrades_focus_beyond_pose_budget() {
|
||||
use rand::{Rng, SeedableRng};
|
||||
use rand_chacha::ChaCha20Rng;
|
||||
|
||||
let center_freq = 4.0e9;
|
||||
let sweep = FrequencySweep::new(3.0e9, 5.0e9, 32);
|
||||
let target = ScatteringTarget::new(Point3::new(0.0, 2.0, 0.0), 1.0);
|
||||
let nominal_poses = linear_aperture(Point3::new(-0.5, 0.0, 0.0), Point3::new(0.5, 0.0, 0.0), 21);
|
||||
let budget = max_coherent_pose_error_m(center_freq);
|
||||
let lambda = wavelength_m(center_freq);
|
||||
|
||||
// A single FIXED per-position error *direction* (random sign along
|
||||
// each antenna's boresight to the target, drawn once), then scaled by
|
||||
// a growing `epsilon`. This isolates "how does focus respond as
|
||||
// position-error magnitude grows" from "which specific random error
|
||||
// pattern did we happen to draw" -- redrawing a fresh random pattern
|
||||
// at every epsilon level (tried first) makes neighboring levels
|
||||
// statistically incomparable and the trend noisy enough to need heavy
|
||||
// Monte Carlo averaging. Levels are kept within half a wavelength
|
||||
// (4x budget = lambda/2): beyond that, per-position phase error wraps
|
||||
// past 2*pi and can partially and coincidentally realign at specific
|
||||
// epsilon values (a real grating/aliasing effect, not a test bug) --
|
||||
// an honest reason to keep this test inside the regime the lambda/8
|
||||
// budget is actually about, rather than claiming a monotonic trend
|
||||
// the underlying physics doesn't guarantee once error exceeds ~lambda.
|
||||
let mut sign_rng = ChaCha20Rng::seed_from_u64(99);
|
||||
let signs: Vec<f64> = nominal_poses.iter().map(|_| if sign_rng.gen_bool(0.5) { 1.0 } else { -1.0 }).collect();
|
||||
let jittered_poses_at = |epsilon: f64| -> Vec<AntennaPose> {
|
||||
nominal_poses
|
||||
.iter()
|
||||
.zip(&signs)
|
||||
.map(|(p, &sign)| {
|
||||
let dir = p.position.direction_to(&target.position).expect("pose must not coincide with target");
|
||||
AntennaPose::new(p.position.translated(dir, sign * epsilon))
|
||||
})
|
||||
.collect()
|
||||
};
|
||||
|
||||
let freqs = sweep.frequencies();
|
||||
let focus_at_epsilon = |epsilon: f64| -> f64 {
|
||||
let true_poses = jittered_poses_at(epsilon);
|
||||
// The measurement is recorded at the (jittered) TRUE antenna
|
||||
// positions, but reconstruction always assumes the NOMINAL
|
||||
// (design) positions -- the real-world scenario of an
|
||||
// uncalibrated / imperfectly tracked antenna trajectory.
|
||||
// Evaluate focus exactly AT the true target location (not a
|
||||
// grid-wide peak search, which can hop to a nearby voxel that
|
||||
// happens to focus slightly better and mask the coherence loss
|
||||
// this test is measuring).
|
||||
let measurement = simulate_measurement(&true_poses, &sweep, &[target], 0.0, 1);
|
||||
focus_at_point(&measurement, &nominal_poses, &freqs, &target.position)
|
||||
};
|
||||
|
||||
let levels = [0.0, 0.5 * budget, budget, 2.0 * budget, 4.0 * budget];
|
||||
assert!(4.0 * budget < lambda / 2.0 + 1e-12, "test setup: must stay within half a wavelength to avoid phase-wrap aliasing");
|
||||
let focus: Vec<f64> = levels.iter().map(|&eps| focus_at_epsilon(eps)).collect();
|
||||
|
||||
assert!(
|
||||
focus[0] == focus.iter().cloned().fold(0.0, f64::max),
|
||||
"perfect pose knowledge (zero jitter) must give the best focus of the sweep: {focus:?} at levels {levels:?}"
|
||||
);
|
||||
assert!(
|
||||
focus[4] < 0.85 * focus[0],
|
||||
"position error at 4x the lambda/8 budget ({:.4} m, still under half a wavelength) must measurably degrade focus: {:.4} vs zero-jitter {:.4}",
|
||||
4.0 * budget,
|
||||
focus[4],
|
||||
focus[0]
|
||||
);
|
||||
assert!(
|
||||
focus[4] <= focus[1] + 1e-9,
|
||||
"focus at 4x the budget should be no better than focus at 0.5x the budget: {focus:?} at levels {levels:?}"
|
||||
);
|
||||
}
|
||||
Reference in New Issue
Block a user