Rework radar system
This commit is contained in:
Ayzen
@@ -1,18 +0,0 @@
|
||||
"""Extract production CUDA kernels and host helpers for the test-only adapter."""
|
||||
|
||||
from pathlib import Path
|
||||
import re
|
||||
import sys
|
||||
|
||||
|
||||
source = (Path(__file__).resolve().parents[1] / "processing/processors/src/sar3d_gpu.cpp").read_text()
|
||||
kernel = re.search(r'kKernelSource\s*=\s*R"CUDA\((.*?)\)CUDA";', source, re.DOTALL)
|
||||
if kernel is None:
|
||||
raise RuntimeError("SAR CUDA kernel literal was not found")
|
||||
first = source.index("constexpr std::size_t kMaxIndexedItems")
|
||||
last = source.index("} // namespace", first)
|
||||
destination = Path(sys.argv[1])
|
||||
destination.mkdir(parents=True, exist_ok=True)
|
||||
(destination / "sar3d_kernels.inc").write_text(kernel.group(1))
|
||||
helpers = re.sub(r"^#(?:ifdef SAR3D_WITH_CUDA|endif)$", "", source[first:last], flags=re.MULTILINE)
|
||||
(destination / "sar3d_helpers.inc").write_text(helpers)
|
||||
@@ -0,0 +1,204 @@
|
||||
"""Generate small SAR3D references by executing notebook cells 0 and 4.
|
||||
|
||||
Requires the original, untracked Multi_antenna_3D_clean_v1.ipynb; a missing notebook
|
||||
is an error. NumPy/SciPy execute its function/class ASTs unchanged. Only parameter
|
||||
values are replaced; plotting imports and the interactive top-level runs are
|
||||
excluded. No production SAR implementation or hand-transcribed oracle is used.
|
||||
|
||||
Run from the repository root:
|
||||
.venv/bin/python data_acq_and_processing/tests/generate_sar3d_reference.py
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import ast
|
||||
import hashlib
|
||||
import json
|
||||
from pathlib import Path
|
||||
import tempfile
|
||||
|
||||
import numpy as np
|
||||
import scipy
|
||||
|
||||
ROOT = Path(__file__).resolve().parents[2]
|
||||
REFERENCE = Path(__file__).with_name("sar3d_reference.npz")
|
||||
|
||||
|
||||
def notebook_functions(notebook: Path, parameters: dict) -> tuple[dict, dict]:
|
||||
cells = json.loads(notebook.read_text(encoding="utf-8"))["cells"]
|
||||
namespace = {"__name__": __name__}
|
||||
hashes = {}
|
||||
for index in (0, 4):
|
||||
source = "".join(cells[index]["source"])
|
||||
hashes[str(index)] = hashlib.sha256(source.encode()).hexdigest()
|
||||
tree = ast.parse(source, filename=f"{notebook.name}:cell{index}")
|
||||
definitions = []
|
||||
for node in tree.body:
|
||||
if isinstance(node, (ast.FunctionDef, ast.ClassDef, ast.ImportFrom)):
|
||||
definitions.append(node)
|
||||
elif isinstance(node, ast.Import):
|
||||
if not any(alias.name.startswith(("matplotlib", "plotly")) for alias in node.names):
|
||||
definitions.append(node)
|
||||
elif isinstance(node, ast.Assign) and len(node.targets) == 1 and isinstance(node.targets[0], ast.Name):
|
||||
try:
|
||||
namespace[node.targets[0].id] = ast.literal_eval(node.value)
|
||||
except (ValueError, TypeError):
|
||||
pass # Dataset loading and the interactive run are supplied below.
|
||||
namespace.update(parameters)
|
||||
exec(compile(ast.Module(definitions, type_ignores=[]), f"{notebook.name}:cell{index}",
|
||||
"exec", dont_inherit=True), namespace)
|
||||
return namespace, hashes
|
||||
|
||||
|
||||
def scene() -> dict:
|
||||
"""Deterministic two-echo input with motion, vertical baselines and four pairs."""
|
||||
request = dict(
|
||||
tx=[[0, -.4, 0, 0], [1, .45, .2, .05]],
|
||||
rx=[[0, -.2, .1, 0], [1, .3, .35, .02]], f_start=1e9, f_stop=4e9,
|
||||
min_depth=1., max_depth=3., speed=.4, pitch=13., ground=.3,
|
||||
voxel=.4, radius=1, power=1.5, compensation=True, normalize=True,
|
||||
)
|
||||
# Match the float32 settings/geometry accepted by the native fixture.
|
||||
for key in ("tx", "rx"):
|
||||
request[key] = np.asarray(request[key], dtype=np.float32).astype(float).tolist()
|
||||
frequency = np.linspace(1e9, 4e9, 201, dtype=np.float32)
|
||||
members = []
|
||||
for frame, offset in enumerate((0, 180_000_000, 510_000_000)):
|
||||
traces = []
|
||||
start = 10_000_000_000 + offset
|
||||
for tx in range(2):
|
||||
for rx in range(2):
|
||||
samples = ((1 + .1 * rx) * np.exp(-2j * np.pi * frequency.astype(float)
|
||||
* (4.8 + .08 * tx - .03 * frame) / 3e8)
|
||||
+ .35 * np.exp(-2j * np.pi * frequency.astype(float) * (5.6 - .05 * rx) / 3e8)).astype(np.complex64)
|
||||
traces.append(dict(input=rx, output=tx, start=start, end=start + 165_000_000,
|
||||
frequency=frequency.tolist(), s21=np.column_stack((samples.real, samples.imag)).tolist()))
|
||||
members.append(dict(id=frame + 1, timestamp=start + 165_000_000, traces=traces))
|
||||
request.update(previous=members[:-1], latest=members[-1])
|
||||
return request
|
||||
|
||||
|
||||
def write_dataset(directory: Path, request: dict) -> list[Path]:
|
||||
"""Supply the notebook's original config_profile/meta/NPY input contract."""
|
||||
profile = {"gpr": {}, "run": {"combos": []}}
|
||||
for kind, position in (("tx", "output_pos"), ("rx", "input_pos")):
|
||||
profile["gpr"][f"{kind}_geometry"] = [
|
||||
{position: int(index), "x_m": x, "y_m": y, "z_m": z}
|
||||
for index, x, y, z in request[kind]
|
||||
]
|
||||
profile["run"]["combos"] = [
|
||||
{"input": trace["input"], "output": trace["output"]} for trace in request["latest"]["traces"]
|
||||
]
|
||||
(directory / "config_profile.json").write_text(json.dumps(profile), encoding="utf-8")
|
||||
snapshots = []
|
||||
for index, member in enumerate([*request["previous"], request["latest"]]):
|
||||
snapshot = directory / f'{index:04d}_id{member["id"]}_ns{member["timestamp"]}'
|
||||
snapshot.mkdir()
|
||||
metadata = {"traces": []}
|
||||
for trace in member["traces"]:
|
||||
stem = f'i{trace["input"]}_o{trace["output"]}'
|
||||
values = np.asarray(trace["s21"], dtype=np.float32)
|
||||
np.save(snapshot / f"{stem}_s21.npy", values[:, 0] + 1j * values[:, 1])
|
||||
np.save(snapshot / f"{stem}_freq.npy", np.asarray(trace["frequency"], dtype=np.float32))
|
||||
metadata["traces"].append(dict(input=trace["input"], output=trace["output"],
|
||||
capture_start_ns=trace["start"], capture_end_ns=trace["end"]))
|
||||
(snapshot / "meta.json").write_text(json.dumps(metadata), encoding="utf-8")
|
||||
snapshots.append(snapshot)
|
||||
return snapshots
|
||||
|
||||
|
||||
def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dict]:
|
||||
snapshots = write_dataset(directory, request)
|
||||
center = snapshots[max(len(snapshots) - 1 - request["radius"], 0)]
|
||||
parameters = dict(
|
||||
INPUT_IDX=[int(row[0]) for row in request["rx"]],
|
||||
OUTPUT_IDX=[int(row[0]) for row in request["tx"]],
|
||||
F_START=request["f_start"], F_STOP=request["f_stop"],
|
||||
MIN_DEPTH=request["min_depth"], MAX_DEPTH=request["max_depth"],
|
||||
PAIR_NORM_DEPTH_MIN=request["min_depth"], PAIR_NORM_DEPTH_MAX=request["max_depth"],
|
||||
SPEED_M_S=request["speed"], LOOK_ANGLE_DEG=request["pitch"],
|
||||
BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"],
|
||||
SAR_VOXEL_M=(request["voxel"],) * 3, SAR_FRAME_RADIUS=request["radius"],
|
||||
SAR_USE_COMPENSATION=request["compensation"], SAR_NORMALIZE_EACH_FRAME=request["normalize"],
|
||||
SAR_SCAN_ROOT=directory, DATA_PATH=center, BG_PATH=directory, v=3e8,
|
||||
)
|
||||
ns, hashes = notebook_functions(notebook, parameters)
|
||||
tx, rx, order, tx_order, pair_count = ns["load_antenna_config_profile"](
|
||||
directory / "config_profile.json", ns["INPUT_IDX"], ns["OUTPUT_IDX"])
|
||||
ns["N_tx"], ns["N_rx"] = len(tx), len(rx)
|
||||
for kind, positions in (("tx", tx), ("rx", rx)):
|
||||
for axis, values in zip("xyz", ns["positions_to_xyz"](positions).T):
|
||||
ns[f"{axis}_{kind}"] = values
|
||||
antenna_x = np.concatenate((ns["x_tx"], ns["x_rx"]))
|
||||
ns["x_grid_bp"] = np.linspace(antenna_x.min() - 2, antenna_x.max() + 2, ns["NX_BP"])
|
||||
sweep_time, _, _ = ns["estimate_sweep_time_from_meta"](directory, n_pairs=pair_count)
|
||||
config = ns["MotionConfig"](
|
||||
speed_m_s=ns["SPEED_M_S"], look_angle_deg=ns["LOOK_ANGLE_DEG"],
|
||||
tx_sweep_time_s=sweep_time, frame_acq_time_s=sweep_time * pair_count,
|
||||
pair_order_phys=order, tx_order_full=tuple(tx_order),
|
||||
n_freq_hint=len(request["latest"]["traces"][0]["frequency"]),
|
||||
apply_pair_motion_correction=ns["APPLY_PAIR_MOTION_CORRECTION"],
|
||||
apply_freq_phase_correction=ns["APPLY_FREQ_PHASE_CORRECTION"],
|
||||
pair_motion_correction_mode=ns["PAIR_MOTION_CORRECTION_MODE"],
|
||||
motion_phase_sign=ns["MOTION_PHASE_SIGN"],
|
||||
)
|
||||
timestamps, _ = ns["compute_pair_timestamps_for_snapshot"](center, config)
|
||||
result = ns["sar_compute_volume"](center, None, config, timestamps, verbose=False)
|
||||
volume, grid = result["volume"], result["grid"]
|
||||
iy, iz, ix = ns["sar_global_max_index"](volume, grid)
|
||||
expected = {axis: grid[axis] for axis in "xyz"}
|
||||
expected.update(front=volume.max(axis=1), side=volume.max(axis=2), top=volume.max(axis=0),
|
||||
peak=np.array([result["level"], grid["x"][ix], grid["y"][iy], grid["z"][iz]]))
|
||||
return expected, hashes
|
||||
|
||||
|
||||
def generate(notebook: Path, output: Path) -> None:
|
||||
# Read up front so a missing notebook cannot produce an incomplete artifact.
|
||||
notebook.read_text(encoding="utf-8")
|
||||
cases = [
|
||||
dict(compensation=True, normalize=True),
|
||||
dict(compensation=True, normalize=False, speed=-.2, pitch=-7., ground=.5),
|
||||
dict(compensation=False, normalize=True, f_start=1.3e9, f_stop=3.7e9),
|
||||
dict(compensation=False, normalize=False, speed=0., pitch=0., radius=0),
|
||||
]
|
||||
arrays, configurations = {}, []
|
||||
for index, settings in enumerate(cases):
|
||||
request = scene()
|
||||
request.update(settings)
|
||||
for key in ("min_depth", "max_depth", "speed", "pitch", "ground", "voxel", "power"):
|
||||
request[key] = float(np.float32(request[key]))
|
||||
if index == 1:
|
||||
# A separate, unequal capture window per output instead of interleaving.
|
||||
for member in [*request["previous"], request["latest"]]:
|
||||
for trace in member["traces"]:
|
||||
if trace["output"] == 0:
|
||||
trace["end"] = trace["start"] + 80_000_000
|
||||
else:
|
||||
trace["start"] += 95_000_000
|
||||
with tempfile.TemporaryDirectory(prefix="sar3d_notebook_") as directory:
|
||||
expected, hashes = reference(notebook, request, Path(directory))
|
||||
if index == 0:
|
||||
arrays["request"] = np.array(json.dumps(request, separators=(",", ":")))
|
||||
configurations.append(dict(
|
||||
settings={key: request[key] for key in settings},
|
||||
windows=[[[trace["start"], trace["end"]] for trace in member["traces"]]
|
||||
for member in [*request["previous"], request["latest"]]],
|
||||
))
|
||||
arrays.update({f"case_{index}_{key}": value for key, value in expected.items()})
|
||||
arrays["cases"] = np.array(json.dumps(configurations, separators=(",", ":")))
|
||||
arrays["provenance"] = np.array(json.dumps(dict(
|
||||
notebook=notebook.name, cell_source_sha256=hashes, cells=[0, 4], cases=len(cases),
|
||||
numpy=np.__version__, scipy=scipy.__version__,
|
||||
method="Unchanged notebook AST functions; max projections of sar_compute_volume output",
|
||||
), sort_keys=True))
|
||||
output.parent.mkdir(parents=True, exist_ok=True)
|
||||
np.savez_compressed(output, **arrays)
|
||||
print(f"Wrote {output}: {output.stat().st_size} bytes, {len(cases)} notebook references")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
parser = argparse.ArgumentParser(description=__doc__)
|
||||
parser.add_argument("--notebook", type=Path, default=ROOT / "Multi_antenna_3D_clean_v1.ipynb")
|
||||
parser.add_argument("--output", type=Path, default=REFERENCE)
|
||||
args = parser.parse_args()
|
||||
generate(args.notebook, args.output)
|
||||
@@ -2,10 +2,10 @@
|
||||
#include "channel_bundle_support.hpp"
|
||||
#include "data_preprocessor.hpp"
|
||||
#include "locator/tcp_server.hpp"
|
||||
#include "object_approach_filter.hpp"
|
||||
#include "processing_live_config.hpp"
|
||||
#include "runtime_wait.hpp"
|
||||
#include "shm_ring.hpp"
|
||||
#include "sweep_orchestrator.hpp"
|
||||
#include "tty_protocol_writer.h"
|
||||
|
||||
#include <sys/mman.h>
|
||||
@@ -13,6 +13,7 @@
|
||||
#include <sys/wait.h>
|
||||
#include <unistd.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <atomic>
|
||||
#include <chrono>
|
||||
@@ -65,9 +66,9 @@ struct RingName {
|
||||
std::string name;
|
||||
explicit RingName(const std::string& suffix = {})
|
||||
: name("/radar_native_test_" + std::to_string(::getpid()) + suffix) {
|
||||
ipc::ShmRing::unlink_ring(name);
|
||||
::shm_unlink(name.c_str());
|
||||
}
|
||||
~RingName() { ipc::ShmRing::unlink_ring(name); }
|
||||
~RingName() { ::shm_unlink(name.c_str()); }
|
||||
};
|
||||
|
||||
auto save_bundle(const Workspace& workspace, const char* name, const ipc::RawSweepCollection& collection)
|
||||
@@ -84,13 +85,15 @@ void write_u32(std::vector<std::uint8_t>& bytes, std::size_t offset, std::uint32
|
||||
std::memcpy(bytes.data() + offset, &value, sizeof(value));
|
||||
}
|
||||
|
||||
void test_wire_decoders() {
|
||||
void test_wire_format() {
|
||||
ipc::RawSweepCollection raw{.collection_id = 7, .monotonic_ns = 42,
|
||||
.capture_start_ns = 21, .capture_end_ns = 40,
|
||||
.traces = {{.combo = {1, 2}, .frequency_hz = {1, 2},
|
||||
.s11 = {{1, 2}, {3, 4}}, .s21 = {{5, 6}, {7, 8}},
|
||||
.s21 = {{5, 6}, {7, 8}},
|
||||
.capture_start_ns = 25, .capture_end_ns = 39}}};
|
||||
ipc::ResultCollection results{.collection_id = 7, .processing_duration_ns = 17,
|
||||
ipc::ResultCollection results{.collection_id = 7, .monotonic_ns = 0x0102030405060708ULL,
|
||||
.processing_duration_ns = 17, .settings_revision = 0xfedcba9876543210ULL,
|
||||
.capture_start_ns = 21, .capture_end_ns = 40, .speed_m_s = -1.5F,
|
||||
.collection_payloads = {
|
||||
{.kind = ipc::ResultKind::TraceComplex, .frequency_hz = {1}, .trace = {{2, 3}}},
|
||||
{.kind = ipc::ResultKind::ScalarF32, .scalar_value = 4},
|
||||
@@ -103,32 +106,39 @@ void test_wire_decoders() {
|
||||
const auto raw_bytes = ipc::serialize_raw_collection(raw);
|
||||
const auto pre_bytes = ipc::serialize_preprocessed_collection(raw);
|
||||
const auto result_bytes = ipc::serialize_result_collection(results);
|
||||
const std::array<std::uint8_t, 64> result_header{
|
||||
'R', 'S', 'L', '3', 7, 0, 0, 0, 0, 0, 0, 0,
|
||||
8, 7, 6, 5, 4, 3, 2, 1, 17, 0, 0, 0, 0, 0, 0, 0,
|
||||
0x10, 0x32, 0x54, 0x76, 0x98, 0xba, 0xdc, 0xfe,
|
||||
21, 0, 0, 0, 0, 0, 0, 0, 40, 0, 0, 0, 0, 0, 0, 0,
|
||||
0, 0, 0xc0, 0xbf, 6, 0, 0, 0, 1, 0, 0, 0,
|
||||
};
|
||||
require(std::equal(result_header.begin(), result_header.end(), result_bytes.begin()),
|
||||
"RSL3 must preserve exact acquisition identity and settings in its 64-byte header");
|
||||
require(std::memcmp(raw_bytes.data(), "RAW3", 4) == 0 &&
|
||||
std::memcmp(pre_bytes.data(), "PRP3", 4) == 0,
|
||||
"trace collections must use the S21-only version 3 format");
|
||||
// Collection header + one trace header + two 12-byte points + capture trailer.
|
||||
require(raw_bytes.size() == 24 + 12 + 2 * 12 + 16 + 4 + 16 && pre_bytes.size() == raw_bytes.size(),
|
||||
"trace wire format must contain only frequency and S21 samples");
|
||||
auto old_raw = raw_bytes;
|
||||
auto old_pre = pre_bytes;
|
||||
write_u32(old_raw, 0, 0x32574152U);
|
||||
write_u32(old_pre, 0, 0x32525050U);
|
||||
expect_error([&] { (void)ipc::deserialize_raw_collection(old_raw); });
|
||||
expect_error([&] { (void)ipc::deserialize_preprocessed_collection(old_pre); });
|
||||
require(ipc::serialize_raw_collection(ipc::deserialize_raw_collection(raw_bytes)) == raw_bytes,
|
||||
"raw roundtrip changed data");
|
||||
require(ipc::serialize_preprocessed_collection(ipc::deserialize_preprocessed_collection(pre_bytes)) == pre_bytes,
|
||||
"preprocessed roundtrip changed data");
|
||||
require(ipc::serialize_result_collection(ipc::deserialize_result_collection(result_bytes)) == result_bytes,
|
||||
"result roundtrip changed data or empty shape");
|
||||
|
||||
for (std::size_t size = 0; size < result_bytes.size(); ++size) {
|
||||
expect_error([&] { (void)ipc::deserialize_result_collection(std::span(result_bytes).first(size)); });
|
||||
}
|
||||
for (std::size_t count_offset : {28U, 32U, 39U}) {
|
||||
auto corrupt = result_bytes;
|
||||
write_u32(corrupt, count_offset, std::numeric_limits<std::uint32_t>::max());
|
||||
expect_error([&] { (void)ipc::deserialize_result_collection(corrupt); });
|
||||
}
|
||||
for (auto kind : {ipc::ResultKind::ImageF32, ipc::ResultKind::TableF32}) {
|
||||
ipc::ResultCollection empty{.collection_payloads = {{.kind = kind}}};
|
||||
auto corrupt = ipc::serialize_result_collection(empty);
|
||||
write_u32(corrupt, 39, std::numeric_limits<std::uint32_t>::max());
|
||||
write_u32(corrupt, 43, std::numeric_limits<std::uint32_t>::max());
|
||||
expect_error([&] { (void)ipc::deserialize_result_collection(corrupt); });
|
||||
for (const std::size_t tail : {36U, 20U}) {
|
||||
expect_error([&] { (void)ipc::deserialize_raw_collection(std::span(raw_bytes).first(raw_bytes.size() - tail)); });
|
||||
expect_error([&] { (void)ipc::deserialize_preprocessed_collection(std::span(pre_bytes).first(pre_bytes.size() - tail)); });
|
||||
}
|
||||
// Deterministic random mutations include wire-compatible valid frames as seeds.
|
||||
std::mt19937 random(0x52414452U);
|
||||
for (unsigned iteration = 0; iteration < 10'000U; ++iteration) {
|
||||
auto bytes = iteration % 3U == 0U ? raw_bytes : iteration % 3U == 1U ? pre_bytes : result_bytes;
|
||||
auto bytes = iteration % 2U == 0U ? raw_bytes : pre_bytes;
|
||||
for (unsigned edit = 0; edit < 1U + random() % 5U; ++edit) {
|
||||
bytes[random() % bytes.size()] = static_cast<std::uint8_t>(random());
|
||||
}
|
||||
@@ -136,9 +146,8 @@ void test_wire_decoders() {
|
||||
bytes.resize(random() % bytes.size());
|
||||
}
|
||||
try {
|
||||
if (iteration % 3U == 0U) { (void)ipc::deserialize_raw_collection(bytes); }
|
||||
else if (iteration % 3U == 1U) { (void)ipc::deserialize_preprocessed_collection(bytes); }
|
||||
else { (void)ipc::deserialize_result_collection(bytes); }
|
||||
if (iteration % 2U == 0U) { (void)ipc::deserialize_raw_collection(bytes); }
|
||||
else { (void)ipc::deserialize_preprocessed_collection(bytes); }
|
||||
} catch (const std::bad_alloc&) { throw; }
|
||||
catch (const std::runtime_error&) {}
|
||||
}
|
||||
@@ -151,12 +160,35 @@ void test_configuration(const std::string& seed_path, const Workspace& workspace
|
||||
const auto save = [&](const Json& config) { std::ofstream(config_path) << config; };
|
||||
save(seed);
|
||||
const auto config = config::load_run_config(config_path.string());
|
||||
auto no_taps = seed;
|
||||
no_taps["rings"]["preprocessed_tap"]["name"] = "";
|
||||
save(no_taps);
|
||||
const auto disabled = config::load_run_config(config_path.string()).rings;
|
||||
require(disabled.preprocessed_tap.name.empty(),
|
||||
"explicit empty tap names must disable the optional rings");
|
||||
for (const auto* key : {"raw", "preprocessed", "results"}) {
|
||||
auto invalid = no_taps;
|
||||
invalid["rings"][key]["name"] = "";
|
||||
save(invalid);
|
||||
expect_error([&] { (void)config::load_run_config(config_path.string()); });
|
||||
}
|
||||
for (const auto* port : {"port1", "port2"}) {
|
||||
auto with_driver = seed;
|
||||
with_driver["switches"][port]["driver"] = "h7992";
|
||||
save(with_driver);
|
||||
(void)config::load_run_config(config_path.string());
|
||||
for (const auto& driver : {Json("hmc349a"), Json("unknown"), Json(7)}) {
|
||||
with_driver["switches"][port]["driver"] = driver;
|
||||
save(with_driver);
|
||||
expect_error([&] { (void)config::load_run_config(config_path.string()); });
|
||||
}
|
||||
}
|
||||
auto with_object_calibration = seed;
|
||||
with_object_calibration["preprocess"]["s21"]["object_calibration"] =
|
||||
Json{{"set_name", "reflector"}, {"bundle_path", "corrections.bin"}};
|
||||
save(with_object_calibration);
|
||||
const auto selected = config::load_run_config(config_path.string()).preprocess.s21.object_calibration;
|
||||
require(selected.set_name == "reflector" && selected.bundle_path == "corrections.bin",
|
||||
require(selected.bundle_path == "corrections.bin",
|
||||
"object calibration asset did not load");
|
||||
with_object_calibration["preprocess"]["s21"].erase("object_calibration");
|
||||
save(with_object_calibration);
|
||||
@@ -167,10 +199,8 @@ void test_configuration(const std::string& seed_path, const Workspace& workspace
|
||||
save(with_object_calibration);
|
||||
expect_error([&] { (void)config::load_run_config(config_path.string()); });
|
||||
}
|
||||
for (const auto& [pointer, value] : std::array<std::pair<const char*, double>, 8>{{
|
||||
{"/radar/sweep/start_hz", -1}, {"/radar/sweep/stop_hz", 1e100},
|
||||
{"/radar/sweep/if_bandwidth_hz", 0}, {"/radar/sweep/points", 0},
|
||||
{"/radar/sweep/points", 4294967296.0}, {"/run/combos/0/input", 0.0000001},
|
||||
for (const auto& [pointer, value] : std::array<std::pair<const char*, double>, 3>{{
|
||||
{"/run/combos/0/input", 0.0000001},
|
||||
{"/rings/raw/slot_size_bytes", 65}, {"/gpr/rx_geometry/0/x_m", 1e100}}}) {
|
||||
auto invalid = seed;
|
||||
invalid[Json::json_pointer(pointer)] = value;
|
||||
@@ -190,11 +220,6 @@ void test_configuration(const std::string& seed_path, const Workspace& workspace
|
||||
require(current.bscan_gain == 2 && loader.revision() == revision,
|
||||
"invalid live config replaced the last valid settings");
|
||||
}
|
||||
const std::array<float, 1> nan_axis{std::numeric_limits<float>::quiet_NaN()};
|
||||
require(!preprocessing::frequency_axes_match(nan_axis, nan_axis), "NaN frequency axes matched");
|
||||
preprocessing::ChannelTrace invalid_trace{.frequency_hz = {1},
|
||||
.samples = {{std::numeric_limits<float>::infinity(), 0}}};
|
||||
expect_error([&] { preprocessing::validate_channel_trace_layout(invalid_trace, "test", {}); });
|
||||
}
|
||||
|
||||
void test_ring() {
|
||||
@@ -245,38 +270,16 @@ void test_ring() {
|
||||
}
|
||||
}
|
||||
|
||||
void test_partial_driver_open() {
|
||||
struct Radar final : drivers::RadarDriver {
|
||||
unsigned closes = 0;
|
||||
void open() override {}
|
||||
void close() override { ++closes; }
|
||||
auto acquire_sweep() -> drivers::SweepTrace override { return {}; }
|
||||
} radar;
|
||||
struct Switch final : drivers::SwitchDriver {
|
||||
std::atomic<bool>& stop;
|
||||
bool fail_open = false;
|
||||
bool fail_close = false;
|
||||
unsigned closes = 0;
|
||||
explicit Switch(std::atomic<bool>& stop_flag) : stop(stop_flag) {}
|
||||
void open() override {
|
||||
if (fail_open) { stop.store(true); throw std::runtime_error("open failure"); }
|
||||
}
|
||||
void close() override { ++closes; if (fail_close) { throw std::runtime_error("close failure"); } }
|
||||
auto position_count() const -> std::uint32_t override { return 1; }
|
||||
void switch_to(std::uint32_t) override {}
|
||||
};
|
||||
std::atomic<bool> stop{false};
|
||||
Switch input(stop), output(stop);
|
||||
input.fail_open = true;
|
||||
output.fail_close = true;
|
||||
config::RunConfig config;
|
||||
config.radar.sweep.points = 1;
|
||||
config.run_combos = {{0, 0}};
|
||||
RingName name;
|
||||
auto ring = ipc::ShmRing::open_or_create(name.name, 2, 256);
|
||||
acq::SweepOrchestrator(config, radar, input, output, ring).run(stop);
|
||||
require(radar.closes == 1 && input.closes == 1 && output.closes == 1,
|
||||
"partial driver startup did not close each driver exactly once");
|
||||
void test_acquisition_identity() {
|
||||
processing::ObjectApproachFilter filter;
|
||||
const std::vector<processing::ObjectApproachFilter::Point> points{{0, 2}};
|
||||
require(!filter.confirm(points, 1, 1'000'000'000ULL, 0, 0, 2).front(), "first frame confirmed a track");
|
||||
require(filter.confirm(points, 1, 2'000'000'000ULL, 0, 0, 2).front(),
|
||||
"reused collection id replaced a distinct acquisition");
|
||||
require(filter.confirm(points, 1, 2'000'000'000ULL, 0, 0, 2).front(),
|
||||
"reprocessing duplicated a frame or lost its predecessor");
|
||||
require(!filter.confirm(points, 7, 1'000'000'000ULL, 0, 0, 2).front(),
|
||||
"an earlier capture time failed to reset replay history");
|
||||
}
|
||||
|
||||
void test_interruptible_wait() {
|
||||
@@ -291,6 +294,19 @@ void test_interruptible_wait() {
|
||||
"large wait delayed shutdown");
|
||||
}
|
||||
|
||||
void test_cosine_notch_taper() {
|
||||
// Band 2.0-2.1 GHz with 0.1 GHz shoulders: 1 at the outer edges, 0.5 halfway,
|
||||
// 0 from the band edges inwards.
|
||||
const std::vector<float> frequency_hz{1.85e9F, 1.9e9F, 1.95e9F, 2.0e9F, 2.05e9F,
|
||||
2.1e9F, 2.15e9F, 2.2e9F, 2.25e9F};
|
||||
const std::vector<float> expected{1.0F, 1.0F, 0.5F, 0.0F, 0.0F, 0.0F, 0.5F, 1.0F, 1.0F};
|
||||
const auto mask = preprocessing::build_notch_mask(frequency_hz,
|
||||
{.enabled = true, .bands_hz = {{2.0e9F, 2.1e9F}}, .taper_width_hz = 1e8F, .taper_type = "cosine"});
|
||||
for (std::size_t index = 0U; index < expected.size(); ++index) {
|
||||
require(std::abs(mask[index] - expected[index]) < 1e-3F, "cosine notch shoulders are not 1 -> 0 towards the band");
|
||||
}
|
||||
}
|
||||
|
||||
void test_csv_failures(const Workspace& workspace) {
|
||||
const auto csv = workspace.path / "capture.csv";
|
||||
const auto spool = workspace.path / "capture.capture_spool.bin";
|
||||
@@ -328,37 +344,32 @@ void test_csv_failures(const Workspace& workspace) {
|
||||
void test_asset_reload_preserves_data(const Workspace& workspace) {
|
||||
const auto save = [&](const char* name, float sample) {
|
||||
return save_bundle(workspace, name, {.traces = {
|
||||
{.frequency_hz = {1}, .s11 = {{sample, 0}}, .s21 = {{sample, 0}}}}});
|
||||
{.frequency_hz = {1}, .s21 = {{sample, 0}}}}});
|
||||
};
|
||||
const auto open = save("open.bin", 1);
|
||||
const auto short_standard = save("short.bin", -1);
|
||||
const auto load = save("load.bin", 0);
|
||||
const auto nan = save("nan.bin", std::numeric_limits<float>::quiet_NaN());
|
||||
preprocessing::S21ReferenceBundle s21_reference;
|
||||
s21_reference.load(open);
|
||||
expect_error([&] { s21_reference.load(nan); });
|
||||
require(s21_reference.resolve({}).samples.front().re == 1, "failed S21 reload lost previous reference");
|
||||
preprocessing::S21CalibrationBundle s21_calibration;
|
||||
s21_calibration.load(open);
|
||||
s21_calibration.load(save("through.bin", 2));
|
||||
expect_error([&] { s21_calibration.load(nan); });
|
||||
require(s21_calibration.is_enabled(), "failed S21 reload disabled previous calibration");
|
||||
preprocessing::S11ReferenceBundle s11_reference;
|
||||
s11_reference.load(open);
|
||||
expect_error([&] { s11_reference.load(nan); });
|
||||
require(s11_reference.resolve({}).samples.front().re == 1, "failed S11 reload lost previous reference");
|
||||
preprocessing::S11CalibrationBundle s11_calibration;
|
||||
s11_calibration.load(open, short_standard, load);
|
||||
expect_error([&] { s11_calibration.load(open, open, load); });
|
||||
const std::array<float, 1> axis{1};
|
||||
require(s11_calibration.apply({}, axis, {{0.5, 0}}).front().re == 0.5F,
|
||||
"failed OSL reload lost previous calibration");
|
||||
const std::array<float, 1> frequency{1};
|
||||
std::array<ipc::Complex32, 1> samples{{{6, 8}}};
|
||||
s21_calibration.apply({}, frequency, samples);
|
||||
require(samples[0].re == 3 && samples[0].im == 4, "failed S21 reload changed previous calibration");
|
||||
const std::array<float, 1> nan_axis{std::numeric_limits<float>::quiet_NaN()};
|
||||
expect_error([&] { s21_calibration.apply({}, nan_axis, samples); });
|
||||
preprocessing::S21ReferenceBundle reference;
|
||||
reference.load(open, s21_calibration, {{}});
|
||||
expect_error([&] { reference.load(nan, s21_calibration, {{}}); });
|
||||
ipc::SweepTraceBlock trace{.frequency_hz = {1}, .s21 = {{2, 0}}};
|
||||
reference.apply(trace);
|
||||
require(trace.s21[0].re == 1.5F, "failed reference reload changed the calibrated reference");
|
||||
}
|
||||
|
||||
void test_object_calibration(const Workspace& workspace) {
|
||||
const ipc::SweepTraceBlock first{.combo = {7, 3}, .frequency_hz = {1e9F, 2e9F},
|
||||
.s11 = {{1, 2}, {3, 4}}, .s21 = {{2, 3}, {4, 5}}, .capture_start_ns = 5, .capture_end_ns = 15};
|
||||
.s21 = {{2, 3}, {4, 5}}, .capture_start_ns = 5, .capture_end_ns = 15};
|
||||
const ipc::SweepTraceBlock second{.combo = {1, 9}, .frequency_hz = {3e9F, 4e9F, 5e9F},
|
||||
.s11 = {{}, {}, {}}, .s21 = {{6, 7}, {8, 9}, {10, 11}}};
|
||||
.s21 = {{6, 7}, {8, 9}, {10, 11}}};
|
||||
auto factors_first = first;
|
||||
factors_first.s21 = {{0, 2}, {2, 0}};
|
||||
auto factors_second = second;
|
||||
@@ -390,7 +401,7 @@ void test_object_calibration(const Workspace& workspace) {
|
||||
trace = first;
|
||||
correction.apply(trace);
|
||||
require(ipc::serialize_raw_collection({.traces = {trace}}) == ipc::serialize_raw_collection({.traces = {expected}}),
|
||||
"object calibration changed S11, timing, pair identity, or multiplication after a failed reload");
|
||||
"object calibration changed timing, pair identity, or multiplication after a failed reload");
|
||||
trace = second;
|
||||
correction.apply(trace);
|
||||
require(trace.s21[0].re == 13 && trace.s21[0].im == 1 && trace.s21[1].re == -8 && trace.s21[1].im == -9 &&
|
||||
@@ -406,19 +417,18 @@ void test_object_calibration(const Workspace& workspace) {
|
||||
|
||||
void test_object_calibration_pipeline(const Workspace& workspace) {
|
||||
const ipc::SweepTraceBlock trace{.combo = {7, 3}, .frequency_hz = {1e9F, 2e9F},
|
||||
.s11 = {{7, 3}, {8, 4}}, .s21 = {{2, 4}, {6, 8}}, .capture_start_ns = 5, .capture_end_ns = 15};
|
||||
.s21 = {{2, 4}, {6, 8}}, .capture_start_ns = 5, .capture_end_ns = 15};
|
||||
auto standard = trace;
|
||||
standard.s21 = {{2, 0}, {2, 0}};
|
||||
auto reference = trace;
|
||||
reference.s21 = {{1, 2}, {2, 2}};
|
||||
auto factors = trace;
|
||||
factors.s21 = {{0, 2}, {5, 1}};
|
||||
preprocessing::CalibrationMaster calibration(preprocessing::make_s21_through_calibrator());
|
||||
calibration.load_s21_calibration_bundle(save_bundle(workspace, "pipeline_through.bin", {.traces = {standard}}));
|
||||
preprocessing::ReferenceMaster references;
|
||||
references.load_s21_reference_bundle(save_bundle(workspace, "pipeline_reference.bin", {.traces = {reference}}));
|
||||
references.prepare_calibrated(calibration, {trace.combo});
|
||||
config::RunConfig config;
|
||||
config.preprocess.s21.calibration.bundle_path =
|
||||
save_bundle(workspace, "pipeline_through.bin", {.traces = {standard}});
|
||||
config.preprocess.s21.reference.bundle_path =
|
||||
save_bundle(workspace, "pipeline_reference.bin", {.traces = {reference}});
|
||||
config.run_combos = {trace.combo};
|
||||
config.preprocess.s21.object_calibration.bundle_path =
|
||||
save_bundle(workspace, "pipeline_object.bin", {.traces = {factors}});
|
||||
@@ -426,7 +436,7 @@ void test_object_calibration_pipeline(const Workspace& workspace) {
|
||||
RingName raw_name("_object_raw"), pre_name("_object_pre");
|
||||
auto raw_ring = ipc::ShmRing::open_or_create(raw_name.name, 2, 4096);
|
||||
auto pre_ring = ipc::ShmRing::open_or_create(pre_name.name, 2, 4096);
|
||||
preprocessing::DataPreprocessor preprocessor(config, calibration, references, raw_ring, pre_ring);
|
||||
preprocessing::DataPreprocessor preprocessor(config, raw_ring, pre_ring);
|
||||
require(raw_ring.push(ipc::serialize_raw_collection({.collection_id = 41, .capture_start_ns = 5,
|
||||
.capture_end_ns = 15, .traces = {trace}})), "test raw ring could not hold a sweep");
|
||||
std::atomic<bool> stop{false};
|
||||
@@ -449,8 +459,7 @@ void test_object_calibration_pipeline(const Workspace& workspace) {
|
||||
"object calibration lost capture identity or timing");
|
||||
require(output.s21[0].re == -2 && output.s21[0].im == 1 && output.s21[1].re == 0 && output.s21[1].im == 0,
|
||||
"object calibration ran before existing preprocessing or was applied more than once");
|
||||
require(output.s11[0].re == 7 && output.s11[0].im == 3 && output.s11[1].re == 0 && output.s11[1].im == 0,
|
||||
"object calibration changed S11 or bypassed the notch filter");
|
||||
|
||||
}
|
||||
|
||||
void test_locator_packet_errors() {
|
||||
@@ -506,14 +515,15 @@ int main(int argc, char** argv) {
|
||||
try {
|
||||
require(argc == 2, "Pass the simulator config fixture path");
|
||||
Workspace workspace;
|
||||
test_wire_decoders();
|
||||
test_wire_format();
|
||||
test_configuration(argv[1], workspace);
|
||||
test_asset_reload_preserves_data(workspace);
|
||||
test_object_calibration(workspace);
|
||||
test_object_calibration_pipeline(workspace);
|
||||
test_ring();
|
||||
test_partial_driver_open();
|
||||
test_acquisition_identity();
|
||||
test_interruptible_wait();
|
||||
test_cosine_notch_taper();
|
||||
test_csv_failures(workspace);
|
||||
test_locator_packet_errors();
|
||||
test_unread_tty_shutdown(workspace);
|
||||
|
||||
@@ -1,7 +1,9 @@
|
||||
#include "sar3d_processor.hpp"
|
||||
#include "../processing/processors/src/sar3d_gpu.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <complex>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <nlohmann/json.hpp>
|
||||
|
||||
using Json = nlohmann::json;
|
||||
@@ -29,60 +31,57 @@ auto collection(const Json& source) -> ipc::PreprocessedCollection {
|
||||
for (const auto& sample : complex_samples(entry.at("s21"))) {
|
||||
trace.s21.push_back({sample.real(), sample.imag()});
|
||||
}
|
||||
if (result.traces.empty() || trace.capture_start_ns < result.capture_start_ns) {
|
||||
result.capture_start_ns = trace.capture_start_ns;
|
||||
}
|
||||
result.capture_end_ns = std::max(result.capture_end_ns, trace.capture_end_ns);
|
||||
result.traces.push_back(std::move(trace));
|
||||
}
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
// Runs the sar3d processor once: the JSON request on stdin holds the antenna
|
||||
// geometry, the live settings and the latest and previous collections; the
|
||||
// published payloads go to stdout as JSON.
|
||||
int main() {
|
||||
try {
|
||||
const Json input = Json::parse(std::cin);
|
||||
Json output;
|
||||
if (input.at("operation") == "processor") {
|
||||
config::RunConfig config;
|
||||
for (const auto& row : input.at("tx")) {
|
||||
config.gpr.tx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)});
|
||||
}
|
||||
for (const auto& row : input.at("rx")) {
|
||||
config.gpr.rx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)});
|
||||
}
|
||||
processing::ProcessingLiveConfig live;
|
||||
live.gpr_start_freq_mhz = input.at("f_start").get<float>() / 1e6F;
|
||||
live.gpr_stop_freq_mhz = input.at("f_stop").get<float>() / 1e6F;
|
||||
live.gpr_min_depth_m = input.at("min_depth");
|
||||
live.gpr_max_depth_m = input.at("max_depth");
|
||||
live.gpr_speed_m_s = input.at("speed");
|
||||
live.gpr_look_angle_deg = input.at("pitch");
|
||||
live.sar3d_voxel_m = input.at("voxel");
|
||||
live.sar3d_bp_ground_depth_m = input.at("ground");
|
||||
live.sar3d_frame_radius = input.at("radius");
|
||||
live.sar3d_comp_range_power = input.at("power");
|
||||
live.sar3d_use_compensation = input.at("compensation");
|
||||
live.sar3d_normalize_each_frame = input.at("normalize");
|
||||
std::vector<ipc::PreprocessedCollection> previous;
|
||||
for (const auto& member : input.at("previous")) {
|
||||
previous.push_back(collection(member));
|
||||
}
|
||||
processing::Sar3dProcessor processor;
|
||||
const auto result = processor.process_collection(config, collection(input.at("latest")), previous, live);
|
||||
output["payloads"] = Json::array();
|
||||
for (const auto& payload : result.collection_payloads) {
|
||||
output["payloads"].push_back({{"name", payload.processing_name},
|
||||
{"x", payload.image_x_axis}, {"y", payload.image_y_axis},
|
||||
{"image", payload.image_values}, {"table", payload.table_values},
|
||||
{"columns", payload.table_columns}});
|
||||
}
|
||||
} else {
|
||||
processing::sar3d::Engine engine(1.5);
|
||||
const auto samples = engine.inverse_dft(complex_samples(input.at("bands")),
|
||||
input.at("points"), input.at("first_bin"), input.at("size"));
|
||||
output["ascans"] = Json::array();
|
||||
for (const auto& sample : samples) {
|
||||
output["ascans"].push_back({sample.real(), sample.imag()});
|
||||
}
|
||||
config::RunConfig config;
|
||||
config.rings.results.slot_size_bytes = input.value("slot", config.rings.results.slot_size_bytes);
|
||||
for (const auto& row : input.at("tx")) {
|
||||
config.gpr.tx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)});
|
||||
}
|
||||
std::cout << output << '\n';
|
||||
for (const auto& row : input.at("rx")) {
|
||||
config.gpr.rx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)});
|
||||
}
|
||||
processing::ProcessingLiveConfig live;
|
||||
live.gpr_start_freq_mhz = input.at("f_start").get<float>() / 1e6F;
|
||||
live.gpr_stop_freq_mhz = input.at("f_stop").get<float>() / 1e6F;
|
||||
live.gpr_min_depth_m = input.at("min_depth");
|
||||
live.gpr_max_depth_m = input.at("max_depth");
|
||||
live.gpr_speed_m_s = input.at("speed");
|
||||
live.gpr_look_angle_deg = input.at("pitch");
|
||||
live.sar3d_voxel_m = input.at("voxel");
|
||||
live.sar3d_bp_ground_depth_m = input.at("ground");
|
||||
live.sar3d_frame_radius = input.at("radius");
|
||||
live.sar3d_comp_range_power = input.at("power");
|
||||
live.sar3d_use_compensation = input.at("compensation");
|
||||
live.sar3d_normalize_each_frame = input.at("normalize");
|
||||
std::vector<ipc::PreprocessedCollection> previous;
|
||||
for (const auto& member : input.at("previous")) {
|
||||
previous.push_back(collection(member));
|
||||
}
|
||||
processing::Sar3dProcessor processor;
|
||||
const auto result = processor.process_collection(config, collection(input.at("latest")), previous, live);
|
||||
Json payloads = Json::array();
|
||||
for (const auto& payload : result.collection_payloads) {
|
||||
payloads.push_back({{"name", payload.processing_name},
|
||||
{"x", payload.image_x_axis}, {"y", payload.image_y_axis},
|
||||
{"image", payload.image_values}, {"table", payload.table_values},
|
||||
{"columns", payload.table_columns}});
|
||||
}
|
||||
std::cout << Json{{"payloads", payloads}} << '\n';
|
||||
} catch (const std::exception& error) {
|
||||
std::cout << Json{{"error", error.what()}} << '\n';
|
||||
return 1;
|
||||
|
||||
@@ -1,136 +0,0 @@
|
||||
// Test-only execution of the exact CUDA kernel source, one logical thread at a time.
|
||||
// Linked solely into sar3d_math_fixture; never used as a production backend.
|
||||
#include "../processing/processors/src/sar3d_gpu.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <bit>
|
||||
#include <climits>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
|
||||
#define __global__
|
||||
using std::isfinite;
|
||||
using std::min;
|
||||
struct Index { unsigned x; };
|
||||
Index blockIdx{0}, blockDim{1}, threadIdx{0};
|
||||
struct float2 { float x, y; };
|
||||
auto make_float2(float x, float y) -> float2 { return {x, y}; }
|
||||
void sincospif(float x, float* sine, float* cosine) {
|
||||
*sine = std::sin(x * 3.14159265358979323846F);
|
||||
*cosine = std::cos(x * 3.14159265358979323846F);
|
||||
}
|
||||
auto __powf(float x, float y) -> float { return std::pow(x, y); }
|
||||
auto __float_as_uint(float x) -> unsigned int { return std::bit_cast<unsigned int>(x); }
|
||||
auto __double2int_rn(double x) -> int { return static_cast<int>(std::nearbyint(x)); }
|
||||
void atomicMax(unsigned int* destination, unsigned int value) { *destination = std::max(*destination, value); }
|
||||
#include "sar3d_kernels.inc"
|
||||
|
||||
auto probe_ascan_upper_bound() -> float {
|
||||
blockIdx.x = 0;
|
||||
const double x[]{1073741824.0}, zero[]{0.0};
|
||||
const ::Contribution contribution{{}, {}, 1.0F, 1.0F};
|
||||
float volume = -1.0F;
|
||||
::backproject(nullptr, INT_MAX, &contribution, 1, x, zero, zero, 1, 1, 1,
|
||||
1.0, 0.0, 0.0, 0.0, 1.0, 1.0F, 0, 0.0F, 0.0F, &volume);
|
||||
return volume;
|
||||
}
|
||||
|
||||
namespace {
|
||||
template <class Operation>
|
||||
void each(int count, Operation operation) {
|
||||
for (blockIdx.x = 0; blockIdx.x < static_cast<unsigned>(count); ++blockIdx.x) { operation(); }
|
||||
}
|
||||
|
||||
auto float_pairs(std::span<const std::complex<float>> values) -> std::vector<float2> {
|
||||
std::vector<float2> pairs;
|
||||
pairs.reserve(values.size());
|
||||
for (const auto value : values) { pairs.push_back({value.real(), value.imag()}); }
|
||||
return pairs;
|
||||
}
|
||||
}
|
||||
|
||||
namespace radar::processing::sar3d {
|
||||
namespace {
|
||||
#include "sar3d_helpers.inc"
|
||||
}
|
||||
|
||||
struct Engine::Impl { std::vector<float> taps; };
|
||||
Engine::Engine(double sigma) : impl_(std::make_unique<Impl>()) { impl_->taps = gaussian_taps(sigma); }
|
||||
Engine::~Engine() = default;
|
||||
|
||||
auto Engine::inverse_dft(std::span<const std::complex<float>> bands, std::size_t points,
|
||||
std::size_t first_bin, std::size_t size) -> std::vector<std::complex<float>> {
|
||||
if (points == 0U || bands.empty() || bands.size() % points != 0U || !std::has_single_bit(size) ||
|
||||
size > kMaxIndexedItems || first_bin >= size || points > size - first_bin) {
|
||||
throw std::runtime_error("sar3d invalid inverse FFT dimensions");
|
||||
}
|
||||
const auto count = checked_product(bands.size() / points, size, kMaxIndexedItems, "A-scan stack");
|
||||
const auto input = float_pairs(bands);
|
||||
std::vector<float2> output(count);
|
||||
each(static_cast<int>(count), [&] {
|
||||
::inverse_dft(input.data(), points, first_bin, size, count, output.data());
|
||||
});
|
||||
std::vector<std::complex<float>> result;
|
||||
result.reserve(count);
|
||||
for (const auto value : output) { result.emplace_back(value.x, value.y); }
|
||||
return result;
|
||||
}
|
||||
|
||||
auto Engine::project(const VolumeRequest& request) -> Projections {
|
||||
const int nx = request.x.size(), ny = request.y.size(), nz = request.z.size();
|
||||
const int plane = checked_product(nx, nz, kMaxIndexedItems, "grid plane");
|
||||
const int voxels = checked_product(plane, ny, kMaxIndexedItems, "grid");
|
||||
std::vector<::Contribution> contributions;
|
||||
for (const auto& term : request.contributions) {
|
||||
contributions.push_back(std::bit_cast<::Contribution>(term));
|
||||
}
|
||||
const auto ascans = float_pairs(request.ascans);
|
||||
std::vector<float> volume(voxels), scratch(voxels), levels(ny);
|
||||
std::vector<unsigned int> indices(ny);
|
||||
each(voxels, [&] {
|
||||
::backproject(ascans.data(), request.samples_per_ascan, contributions.data(), contributions.size(),
|
||||
request.x.data(), request.y.data(), request.z.data(), nx, nz, voxels,
|
||||
std::cos(request.pitch_rad), std::sin(request.pitch_rad), request.ground_depth_m,
|
||||
request.min_depth_m, request.max_depth_m, request.samples_per_metre,
|
||||
request.use_compensation, request.comp_range_power, request.comp_weight_max, volume.data());
|
||||
});
|
||||
each(ny, [&] { ::peak_planes(volume.data(), plane, ny, levels.data(), indices.data()); });
|
||||
Projections raw;
|
||||
choose_peak(raw, request, levels);
|
||||
const int radius = impl_->taps.size() / 2;
|
||||
each(voxels, [&] {
|
||||
::smooth(volume.data(), scratch.data(), impl_->taps.data(), radius, nx, nz, voxels, 1);
|
||||
});
|
||||
each(voxels, [&] {
|
||||
::smooth(scratch.data(), volume.data(), impl_->taps.data(), radius, nx, nz, voxels, 0);
|
||||
});
|
||||
each(ny, [&] { ::peak_planes(volume.data(), plane, ny, levels.data(), indices.data()); });
|
||||
Projections result;
|
||||
result.raw_level = raw.level;
|
||||
choose_peak(result, request, levels, indices);
|
||||
result.front.resize(ny * nx);
|
||||
result.side.resize(ny * nz);
|
||||
result.top.resize(plane);
|
||||
const auto& view = request.oblique;
|
||||
std::vector<unsigned int> oblique(checked_product(view.nu, view.nv, kMaxIndexedItems, "oblique image"));
|
||||
each(result.front.size(), [&] {
|
||||
::max_project(volume.data(), result.front.data(), result.front.size(), nx, plane, 1, nz, nx);
|
||||
});
|
||||
each(result.side.size(), [&] {
|
||||
::max_project(volume.data(), result.side.data(), result.side.size(), nz, plane, nx, nx, 1);
|
||||
});
|
||||
each(result.top.size(), [&] {
|
||||
::max_project(volume.data(), result.top.data(), result.top.size(), nx, nx, 1, ny, plane);
|
||||
});
|
||||
each(voxels, [&] {
|
||||
::project_oblique(volume.data(), request.x.data(), request.y.data(), request.z.data(), nx, nz, voxels,
|
||||
view.u[0], view.u[1], view.u[2], view.v[0], view.v[1], view.v[2], view.u0, view.v0,
|
||||
1 / view.pixel_m, view.nu, view.nv, oblique.data());
|
||||
});
|
||||
result.oblique.reserve(oblique.size());
|
||||
for (const auto value : oblique) { result.oblique.push_back(std::bit_cast<float>(value)); }
|
||||
return result;
|
||||
}
|
||||
}
|
||||
@@ -1,85 +0,0 @@
|
||||
// Kernel arithmetic checks linked only with the host emulation fixture.
|
||||
#include "../processing/processors/src/sar3d_gpu.hpp"
|
||||
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
#include <string_view>
|
||||
|
||||
extern "C" void peak_planes(const float*, int, int, float*, unsigned int*);
|
||||
struct Index { unsigned x; };
|
||||
extern Index blockIdx;
|
||||
auto probe_ascan_upper_bound() -> float;
|
||||
|
||||
namespace {
|
||||
void rejects_unsigned_padding() {
|
||||
blockIdx.x = 0x80000000U;
|
||||
peak_planes(nullptr, 1, 1, nullptr, nullptr);
|
||||
}
|
||||
|
||||
void rejects_rounded_sample_bound() {
|
||||
if (probe_ascan_upper_bound() != 0.0F) {
|
||||
throw std::runtime_error("out-of-range float A-scan position was accepted");
|
||||
}
|
||||
}
|
||||
|
||||
void rejects_nan_peak() {
|
||||
blockIdx.x = 0;
|
||||
const float volume[]{std::numeric_limits<float>::quiet_NaN(), 0.0F};
|
||||
float level = -1.0F;
|
||||
unsigned int index = 0U;
|
||||
peak_planes(volume, 2, 1, &level, &index);
|
||||
if (std::isfinite(level)) { throw std::runtime_error("peak kernel silently ignored NaN"); }
|
||||
}
|
||||
|
||||
void keeps_first_joint_peak() {
|
||||
blockIdx.x = 0;
|
||||
const float volume[]{0, 1, 0, 1, 0, 0};
|
||||
float level = 0;
|
||||
unsigned int index = 0U;
|
||||
peak_planes(volume, 6, 1, &level, &index);
|
||||
if (level != 1.0F || index != 1U) { throw std::runtime_error("joint peak tie chose a nonexistent voxel"); }
|
||||
}
|
||||
|
||||
void rejects_finite_coherent_overflow() {
|
||||
using namespace radar::processing::sar3d;
|
||||
Engine engine(1.5);
|
||||
const std::vector<double> x{-1, 1}, y{-.5, .5}, z{2, 3};
|
||||
std::vector<Contribution> terms(3);
|
||||
for (auto& term : terms) { term.weight = 1; term.inv_reference_ranges = 1; }
|
||||
std::vector<std::complex<float>> scans(3 * 64);
|
||||
for (int i = 0; i < 64; ++i) {
|
||||
scans[i] = {1e38F, 0};
|
||||
scans[64 + i] = {-1e38F, 0};
|
||||
scans[128 + i] = {1, 0};
|
||||
}
|
||||
const VolumeRequest request{
|
||||
.x = x, .y = y, .z = z, .min_depth_m = 1, .max_depth_m = 4,
|
||||
.contributions = terms, .ascans = scans, .samples_per_ascan = 64, .samples_per_metre = 1,
|
||||
.use_compensation = true, .comp_range_power = 1, .comp_weight_max = 5, .norm_edge_margin_m = .3,
|
||||
.oblique = {.u = {1, 0, 0}, .v = {0, 0, 1}, .u0 = -1, .v0 = 2, .pixel_m = 1, .nu = 3, .nv = 2},
|
||||
};
|
||||
try {
|
||||
(void)engine.project(request);
|
||||
} catch (const std::runtime_error& error) {
|
||||
if (std::string_view(error.what()).find("non-finite") == std::string_view::npos) { throw; }
|
||||
return;
|
||||
}
|
||||
throw std::runtime_error("finite coherent overflow was silently accepted as no signal");
|
||||
}
|
||||
}
|
||||
|
||||
int main() {
|
||||
try {
|
||||
rejects_unsigned_padding();
|
||||
rejects_rounded_sample_bound();
|
||||
rejects_nan_peak();
|
||||
keeps_first_joint_peak();
|
||||
rejects_finite_coherent_overflow();
|
||||
std::cout << "SAR CUDA kernel arithmetic regressions passed\n";
|
||||
} catch (const std::exception& error) {
|
||||
std::cerr << error.what() << '\n';
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
Binary file not shown.
@@ -0,0 +1,139 @@
|
||||
"""The sar3d processor on the real CUDA backend: notebook reference and contracts.
|
||||
|
||||
Every case runs build/native_tests/sar3d_fixture (SAR3D_FIXTURE overrides it).
|
||||
Without a usable GPU the processor publishes sar3d_status 3: then only that
|
||||
warning is checked and the remaining tests are skipped. The notebook is needed
|
||||
only to regenerate sar3d_reference.npz with generate_sar3d_reference.py.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import copy
|
||||
import json
|
||||
import os
|
||||
from pathlib import Path
|
||||
import subprocess
|
||||
import unittest
|
||||
|
||||
import numpy as np
|
||||
|
||||
ROOT = Path(__file__).resolve().parents[2]
|
||||
FIXTURE = Path(os.environ.get("SAR3D_FIXTURE", ROOT / "build/native_tests/sar3d_fixture"))
|
||||
REFERENCE = Path(__file__).with_name("sar3d_reference.npz")
|
||||
VIEWS = ("sar3d_front", "sar3d_side", "sar3d_top", "sar3d_oblique")
|
||||
|
||||
|
||||
def run(request: dict) -> dict[str, dict]:
|
||||
output = subprocess.run([str(FIXTURE)], input=json.dumps(request), text=True,
|
||||
capture_output=True, timeout=30, check=True)
|
||||
return {item["name"]: item for item in json.loads(output.stdout)["payloads"]}
|
||||
|
||||
|
||||
def status(payloads: dict[str, dict]) -> tuple[int, str]:
|
||||
table = payloads["sar3d_status"]["table"]
|
||||
return int(table[0]), bytes(map(int, table[1:])).decode("utf-8")
|
||||
|
||||
|
||||
def reference_cases(archive) -> list[dict]:
|
||||
"""Requests of the reference archive, one per notebook parameter case."""
|
||||
cases = []
|
||||
for case in json.loads(archive["cases"].item()):
|
||||
request = json.loads(archive["request"].item())
|
||||
request.update(case["settings"])
|
||||
for member, windows in zip([*request["previous"], request["latest"]], case["windows"], strict=True):
|
||||
for trace, (start, end) in zip(member["traces"], windows, strict=True):
|
||||
trace.update(start=start, end=end)
|
||||
cases.append(request)
|
||||
return cases
|
||||
|
||||
|
||||
class Sar3dTest(unittest.TestCase):
|
||||
@classmethod
|
||||
def setUpClass(cls) -> None:
|
||||
with np.load(REFERENCE, allow_pickle=False) as archive:
|
||||
cls.reference = {key: archive[key] for key in archive.files}
|
||||
cls.cases = reference_cases(cls.reference)
|
||||
cls.scene_payloads = run(cls.cases[0])
|
||||
code, message = status(cls.scene_payloads)
|
||||
cls.unavailable = message if code == 3 else None
|
||||
|
||||
def scene(self, **settings) -> dict:
|
||||
if self.unavailable:
|
||||
self.skipTest(self.unavailable)
|
||||
request = copy.deepcopy(self.cases[0])
|
||||
request.update(settings)
|
||||
return request
|
||||
|
||||
def test_unavailable_gpu_is_a_visible_warning(self) -> None:
|
||||
code, message = status(self.scene_payloads)
|
||||
if self.unavailable:
|
||||
self.assertIn("CUDA", message)
|
||||
self.assertFalse(any(item["image"] for item in self.scene_payloads.values()))
|
||||
else:
|
||||
self.assertEqual(code, 0, message)
|
||||
|
||||
def test_cuda_matches_notebook(self) -> None:
|
||||
if self.unavailable:
|
||||
self.skipTest(self.unavailable)
|
||||
for index, request in enumerate(self.cases):
|
||||
with self.subTest(case=index):
|
||||
payloads = run(request)
|
||||
self.assertEqual(status(payloads)[0], 0, status(payloads))
|
||||
for view, horizontal, vertical in (("front", "x", "y"), ("side", "z", "y"), ("top", "x", "z")):
|
||||
payload = payloads[f"sar3d_{view}"]
|
||||
expected = self.reference[f"case_{index}_{view}"]
|
||||
np.testing.assert_allclose(payload["x"], self.reference[f"case_{index}_{horizontal}"],
|
||||
rtol=2e-6, atol=1e-6)
|
||||
np.testing.assert_allclose(payload["y"], self.reference[f"case_{index}_{vertical}"],
|
||||
rtol=2e-6, atol=1e-6)
|
||||
np.testing.assert_allclose(np.asarray(payload["image"]).reshape(expected.shape),
|
||||
expected, rtol=2e-3, atol=3e-4)
|
||||
peak = payloads["sar3d_peak"]["table"]
|
||||
expected_peak = self.reference[f"case_{index}_peak"]
|
||||
np.testing.assert_allclose(peak[1:], expected_peak[1:], rtol=0, atol=2e-6)
|
||||
np.testing.assert_allclose(peak[0], expected_peak[0], rtol=2e-3, atol=0)
|
||||
|
||||
def test_invalid_input_is_an_error_with_a_reason(self) -> None:
|
||||
def no_traces(request): request["latest"]["traces"] = []
|
||||
def empty_window(request): request["latest"]["traces"][0]["end"] = request["latest"]["traces"][0]["start"]
|
||||
def falling_frequency(request): request["latest"]["traces"][0]["frequency"].reverse()
|
||||
for name, settings, damage in (
|
||||
("no traces", {}, no_traces),
|
||||
("empty capture window", {}, empty_window),
|
||||
("falling frequency", {}, falling_frequency),
|
||||
("grid beyond 32-bit indexing", {"voxel": 1e-8}, None),
|
||||
("vertical look angle", {"pitch": 90}, None),
|
||||
("result larger than the ring slot", {"slot": 1024}, None),
|
||||
):
|
||||
with self.subTest(name):
|
||||
request = self.scene(**settings)
|
||||
if damage:
|
||||
damage(request)
|
||||
code, message = status(run(request))
|
||||
self.assertEqual(code, 2, message)
|
||||
self.assertTrue(message)
|
||||
|
||||
def test_empty_signal_is_reported_without_a_peak(self) -> None:
|
||||
request = self.scene(previous=[], radius=0)
|
||||
for trace in request["latest"]["traces"]:
|
||||
trace["s21"] = [[0, 0]] * len(trace["frequency"])
|
||||
payloads = run(request)
|
||||
self.assertEqual(status(payloads)[0], 1)
|
||||
self.assertNotIn("sar3d_peak", payloads)
|
||||
for name in VIEWS:
|
||||
self.assertFalse(np.any(payloads[name]["image"]))
|
||||
|
||||
def test_weak_signal_stays_finite(self) -> None:
|
||||
request = self.scene(previous=[], radius=0, normalize=False)
|
||||
for trace in request["latest"]["traces"]:
|
||||
trace["s21"] = (np.asarray(trace["s21"]) * 1e-38).tolist()
|
||||
payloads = run(request)
|
||||
self.assertEqual(status(payloads)[0], 0, status(payloads))
|
||||
for name in VIEWS:
|
||||
values = np.asarray(payloads[name]["image"])
|
||||
self.assertTrue(np.isfinite(values).all())
|
||||
self.assertGreater(values.max(), 0)
|
||||
self.assertLess(values.max(), 1e-10)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
@@ -1,297 +0,0 @@
|
||||
"""Compare the native SAR port with notebook cells 0/4, without dataset or GPU fixtures."""
|
||||
from __future__ import annotations
|
||||
|
||||
import json
|
||||
from pathlib import Path
|
||||
import re
|
||||
import subprocess
|
||||
import tempfile
|
||||
import unittest
|
||||
|
||||
import numpy as np
|
||||
|
||||
ROOT = Path(__file__).resolve().parents[2]
|
||||
FIXTURE = ROOT / "build/native_tests/sar3d_math_fixture"
|
||||
CUDA_FIXTURE = ROOT / "build/native_tests/sar3d_fixture"
|
||||
C = 3e8
|
||||
|
||||
|
||||
def smooth(volume):
|
||||
# scipy gaussian_filter(sigma=(0,1.5,1.5), mode='reflect', truncate=4).
|
||||
taps = np.exp(-.5 * (np.arange(-6, 7) / 1.5) ** 2)
|
||||
taps /= taps.sum()
|
||||
for axis in (1, 2):
|
||||
volume = np.apply_along_axis(lambda row: np.convolve(
|
||||
np.pad(row, 6, mode="symmetric"), taps, mode="valid"), axis, volume)
|
||||
return volume
|
||||
|
||||
|
||||
def native(request, executable=None):
|
||||
run = subprocess.run([str(executable or FIXTURE)], input=json.dumps(request), text=True,
|
||||
capture_output=True, timeout=30, check=False)
|
||||
if not run.stdout:
|
||||
raise AssertionError(run.stderr)
|
||||
return json.loads(run.stdout)
|
||||
|
||||
|
||||
def fixture(base=2**64 - 3_000_000_000):
|
||||
request = dict(operation="processor", tx=[[0, -.4, 0, 0], [1, .1, .2, 0], [2, .45, .1, .05]],
|
||||
rx=[[0, -.2, .1, 0], [1, .3, .35, .02]], f_start=1e9, f_stop=4e9,
|
||||
min_depth=1., max_depth=3., speed=.4, pitch=13., ground=.3,
|
||||
voxel=.4, radius=1, power=1.5, compensation=True, normalize=True)
|
||||
members = []
|
||||
frequency = np.linspace(1e9, 4e9, 64, dtype=np.float32)
|
||||
for frame, offset in enumerate([0, 180_000_000, 510_000_000]):
|
||||
traces = []
|
||||
for tx in range(3):
|
||||
for rx in range(2):
|
||||
# Two transmitters interleave in one real capture window; the
|
||||
# third has a longer independent segment, matching meta timing.
|
||||
start = base + offset + (0 if tx < 2 else 70_000_000)
|
||||
end = base + offset + (60_000_000 if tx < 2 else 165_000_001)
|
||||
samples = ((1 + .1 * rx) * np.exp(-2j * np.pi * frequency *
|
||||
(4.8 + .08 * tx - .03 * frame) / C) +
|
||||
.35 * np.exp(-2j * np.pi * frequency * (5.6 - .05 * rx) / C)).astype(np.complex64)
|
||||
traces.append(dict(input=rx, output=tx, start=start, end=end,
|
||||
frequency=frequency.tolist(), s21=np.column_stack((samples.real, samples.imag)).tolist()))
|
||||
members.append(dict(id=frame + 1, timestamp=base + offset, traces=traces))
|
||||
request.update(previous=members[:-1], latest=members[-1])
|
||||
return request
|
||||
|
||||
|
||||
def oracle(request):
|
||||
"""Notebook ground/frequency-domain mode, with actual per-trace capture windows."""
|
||||
pitch = np.deg2rad(float(np.float32(request["pitch"])))
|
||||
cos_p, sin_p = np.cos(pitch), np.sin(pitch)
|
||||
speed, voxel, ground = (float(np.float32(request[key])) for key in ("speed", "voxel", "ground"))
|
||||
axis = lambda lo, hi: np.linspace(lo, hi, max(round((hi - lo) / voxel) + 1, 2))
|
||||
x, y = axis(-2, 2), axis(-.5, 2)
|
||||
limits = [(depth + (height - ground) * sin_p) / cos_p
|
||||
for height in (y[0], y[-1]) for depth in (request["min_depth"], request["max_depth"])]
|
||||
z = axis(min(limits), max(limits))
|
||||
xx, zz = np.meshgrid(x, z)
|
||||
frames, centres = [], []
|
||||
for member in request["previous"] + [request["latest"]]:
|
||||
traces = member["traces"]
|
||||
first, last = min(t["start"] for t in traces), max(t["end"] for t in traces)
|
||||
centres.append((first + last) / 2 if first < 2**53 else (first, last))
|
||||
ascans = []
|
||||
for trace in traces:
|
||||
frequency = np.asarray(trace["frequency"], dtype=np.float32)
|
||||
samples = np.asarray(trace["s21"])
|
||||
samples = samples[:, 0] + 1j * samples[:, 1]
|
||||
mask = (frequency >= request["f_start"]) & (frequency <= request["f_stop"])
|
||||
points = np.flatnonzero(mask)
|
||||
stages = sorted({t["output"] for t in traces
|
||||
if (t["start"], t["end"]) == (trace["start"], trace["end"])})
|
||||
# Relative integer subtraction retains precision at uint64 limits.
|
||||
t = ((trace["start"] - first) - .5 * (last - first)) * 1e-9
|
||||
t += (len(stages) * points + stages.index(trace["output"])) * (
|
||||
(trace["end"] - trace["start"]) * 1e-9 / len(stages) / len(frequency))
|
||||
band = samples[mask] * np.exp(-4j * np.pi * speed * cos_p * frequency[mask] * t / C)
|
||||
df = float(np.median(np.diff(frequency[mask].astype(float))))
|
||||
k0, n = round(float(frequency[mask][0]) / df), len(band)
|
||||
size = 2 ** (2 * (k0 + n - 1) - 1).bit_length() * 8
|
||||
spectrum = np.zeros(size, dtype=complex)
|
||||
spectrum[k0:k0+n] = band * np.hanning(n)
|
||||
h = np.fft.ifft(spectrum)
|
||||
dt = 1 / (size * df)
|
||||
if request["normalize"]:
|
||||
depth = np.arange(size) * dt * C / 2
|
||||
gate = (depth >= request["min_depth"]) & (depth <= request["max_depth"])
|
||||
scale = np.percentile(np.abs(h)[gate], 50)
|
||||
h /= (scale if np.isfinite(scale) and scale > 1e-15 else 1) + 1e-15
|
||||
ascans.append(h)
|
||||
frames.append(ascans)
|
||||
# All members use identical grids. Native centre trails latest by radius.
|
||||
centre_index = max(len(frames) - 1 - request["radius"], 0)
|
||||
def relative(a, b):
|
||||
if isinstance(a, tuple):
|
||||
return .5e-9 * ((a[0] - b[0]) + (a[1] - b[1]))
|
||||
return 1e-9 * (a - b)
|
||||
offsets = np.array([relative(t, centres[centre_index]) for t in centres])
|
||||
farthest = np.max(np.abs(offsets))
|
||||
weights = np.ones(len(frames)) if farthest <= 1e-12 else .8 + .1 * (1 + np.cos(np.pi * abs(offsets) / farthest))
|
||||
acc = np.zeros((len(y), len(z), len(x)), dtype=complex)
|
||||
count = np.zeros(acc.shape)
|
||||
tx = np.asarray(request["tx"], dtype=np.float32)[:, 1:].astype(float)
|
||||
rx = np.asarray(request["rx"], dtype=np.float32)[:, 1:].astype(float)
|
||||
for scans, offset, weight in zip(frames, offsets, weights):
|
||||
shift = np.array([0, speed * offset * sin_p, speed * offset * cos_p])
|
||||
for iy, height in enumerate(y):
|
||||
yr, zr = (height - ground) * cos_p + zz * sin_p, -(height - ground) * sin_p + zz * cos_p
|
||||
depth_ok = (zr >= request["min_depth"]) & (zr <= request["max_depth"])
|
||||
for it, antenna_t in enumerate(tx):
|
||||
for ir, antenna_r in enumerate(rx):
|
||||
def distance(antenna):
|
||||
return np.sqrt((xx-antenna[0])**2 + (yr-antenna[1])**2 + (zr-antenna[2])**2)
|
||||
rt, rr = distance(antenna_t + shift), distance(antenna_r + shift)
|
||||
position = (rt + rr) / (C * dt)
|
||||
h = scans[it * len(rx) + ir]
|
||||
valid = depth_ok & (position <= len(h) - 1)
|
||||
sample = np.interp(position, np.arange(len(h)), h)
|
||||
if request["compensation"]:
|
||||
point = np.array([.5 * (antenna_t[0] + antenna_r[0]),
|
||||
(-ground + 5 * sin_p) / cos_p, 5])
|
||||
ref = 1 / (np.linalg.norm(point-antenna_t) * np.linalg.norm(point-antenna_r) + 1e-12) + 1e-30
|
||||
ratio = (1 / (rt * rr + 1e-12) + 1e-30) / ref
|
||||
sample *= np.clip(1 / (ratio ** request["power"] + 1e-12), 0, 5)
|
||||
acc[iy] += weight * np.where(valid, sample, 0)
|
||||
count[iy] += weight * valid
|
||||
eligible = (y >= y[0]+.3) & (y <= y[-1]-.3)
|
||||
raw_volume = np.abs(acc / (count + 1e-12))
|
||||
raw_level = raw_volume[eligible].max() + 1e-12
|
||||
volume = smooth(raw_volume / raw_level)
|
||||
smooth_level = volume[eligible].max() + 1e-12
|
||||
level = raw_level * smooth_level
|
||||
peak = np.unravel_index(np.argmax(np.where(eligible[:, None, None], volume, -1)), volume.shape)
|
||||
volume /= smooth_level
|
||||
iy, iz, ix = peak
|
||||
eye = np.array([1.7, -2.2, .85])
|
||||
forward = -eye / np.linalg.norm(eye)
|
||||
right = np.cross(forward, [0, 0, 1])
|
||||
right /= np.linalg.norm(right)
|
||||
up = np.cross(right, forward)
|
||||
u, v = right[[0, 2, 1]], up[[0, 2, 1]]
|
||||
points = np.array(np.meshgrid(x, y, z, indexing="ij")).reshape(3, -1).T
|
||||
us, vs = points @ u, points @ v
|
||||
u0, v0 = us.min(), vs.min()
|
||||
nu, nv = round((us.max() - u0) / voxel) + 1, round((vs.max() - v0) / voxel) + 1
|
||||
oblique = np.zeros((nv, nu))
|
||||
pu, pv = np.rint((us-u0)/voxel).astype(int), np.rint((vs-v0)/voxel).astype(int)
|
||||
valid = (pu >= 0) & (pu < nu) & (pv >= 0) & (pv < nv)
|
||||
np.maximum.at(oblique, (pv[valid], pu[valid]), volume.transpose(2, 0, 1).ravel()[valid])
|
||||
return dict(x=x, y=y, z=z, volume=volume, level=level, peak=[x[ix], y[iy], z[iz]],
|
||||
peak_index=peak, oblique=oblique, u=u0 + np.arange(nu)*voxel, v=v0 + np.arange(nv)*voxel)
|
||||
|
||||
|
||||
class Sar3dOracleTest(unittest.TestCase):
|
||||
def test_cuda_branch_cpp_syntax(self):
|
||||
# Hardware-free name/type structure check. Variadic API placeholders do
|
||||
# not validate the NVIDIA ABI or replace testing NVRTC on an actual GPU.
|
||||
source_path = ROOT / "data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp"
|
||||
source = source_path.read_text()
|
||||
with tempfile.TemporaryDirectory(prefix="radar_cuda_syntax_") as directory:
|
||||
path = Path(directory)
|
||||
functions = set(re.findall(r"\b(cu[A-Z]\w*)\s*\(", source))
|
||||
constants = set(re.findall(r"\bCU_DEVICE_ATTRIBUTE_\w+", source))
|
||||
(path / "cuda.h").write_text(
|
||||
"#pragma once\n#include <cstdint>\nenum CUresult {CUDA_SUCCESS=0};\n"
|
||||
"using CUdevice=int;using CUdeviceptr=std::uintptr_t;using CUcontext=struct Context*;"
|
||||
"using CUmodule=struct Module*;using CUfunction=struct Function*;\n" +
|
||||
"".join(f"constexpr int {name}=0;\n" for name in constants) +
|
||||
"".join(f"inline CUresult {name}(...) {{return CUDA_SUCCESS;}}\n" for name in functions))
|
||||
functions = set(re.findall(r"\b(nvrtc[A-Z]\w*)\s*\(", source)) | {
|
||||
"nvrtcGetCUBINSize", "nvrtcGetPTXSize", "nvrtcGetCUBIN", "nvrtcGetPTX"}
|
||||
(path / "nvrtc.h").write_text(
|
||||
"#pragma once\nenum nvrtcResult {NVRTC_SUCCESS=0};using nvrtcProgram=struct Program*;\n"
|
||||
'inline const char* nvrtcGetErrorString(...) {return "stub";}\n' +
|
||||
"".join(f"inline nvrtcResult {name}(...) {{return NVRTC_SUCCESS;}}\n"
|
||||
for name in functions if name != "nvrtcGetErrorString"))
|
||||
result = subprocess.run(["g++", "-std=c++20", "-pthread", "-fsyntax-only", "-Wall", "-Wextra",
|
||||
"-Wpedantic", "-Werror", "-DSAR3D_WITH_CUDA", f"-I{path}", str(source_path)],
|
||||
text=True, capture_output=True, timeout=30, check=False)
|
||||
self.assertEqual(result.returncode, 0, result.stderr)
|
||||
|
||||
def ready(self, result):
|
||||
payloads = {p["name"]: p for p in result["payloads"]}
|
||||
status = payloads["sar3d_status"]["table"]
|
||||
self.assertEqual(payloads["sar3d_status"]["columns"], len(status))
|
||||
self.assertEqual(status[0], 0)
|
||||
return payloads
|
||||
|
||||
def test_inverse_fft_matches_numpy(self):
|
||||
random = np.random.default_rng(209)
|
||||
bands = random.normal(size=(2, 48)) + 1j * random.normal(size=(2, 48))
|
||||
result = native(dict(operation="ifft", bands=np.column_stack((bands.ravel().real, bands.ravel().imag)).tolist(),
|
||||
points=48, first_bin=17, size=512))
|
||||
actual = np.asarray(result["ascans"])
|
||||
spectrum = np.zeros((2, 512), dtype=complex)
|
||||
spectrum[:, 17:65] = bands.astype(np.complex64)
|
||||
expected = np.fft.ifft(spectrum).ravel()
|
||||
np.testing.assert_allclose(actual[:, 0] + 1j * actual[:, 1], expected, rtol=3e-4, atol=2e-7)
|
||||
|
||||
def test_production_backend_reports_cuda_warning(self):
|
||||
result = native(fixture(), CUDA_FIXTURE)
|
||||
payloads = {p["name"]: p for p in result["payloads"]}
|
||||
status = payloads["sar3d_status"]["table"]
|
||||
if status[0] == 3:
|
||||
self.assertIn("CUDA", bytes(map(int, status[1:])).decode("utf-8"))
|
||||
self.assertFalse(any(item["image"] for item in payloads.values()))
|
||||
else:
|
||||
# On a CUDA-equipped test host, also run the real compiled backend.
|
||||
self.ready(result)
|
||||
|
||||
def test_motion_volume_matches_notebook(self):
|
||||
for compensation, normalize in ((True, True), (True, False), (False, True), (False, False)):
|
||||
with self.subTest(compensation=compensation, normalize=normalize):
|
||||
request = fixture()
|
||||
request.update(compensation=compensation, normalize=normalize)
|
||||
payloads = self.ready(native(request))
|
||||
expected = oracle(request)
|
||||
for name, image, x, y in (
|
||||
("sar3d_front", expected["volume"].max(axis=1), expected["x"], expected["y"]),
|
||||
("sar3d_side", expected["volume"].max(axis=2), expected["z"], expected["y"]),
|
||||
("sar3d_top", expected["volume"].max(axis=0), expected["x"], expected["z"]),
|
||||
("sar3d_oblique", expected["oblique"], expected["u"], expected["v"])):
|
||||
np.testing.assert_allclose(payloads[name]["x"], x, rtol=2e-6, atol=1e-6)
|
||||
np.testing.assert_allclose(payloads[name]["y"], y, rtol=2e-6, atol=1e-6)
|
||||
np.testing.assert_allclose(np.asarray(payloads[name]["image"]).reshape(image.shape), image,
|
||||
rtol=2e-3, atol=3e-4)
|
||||
np.testing.assert_allclose(payloads["sar3d_peak"]["table"][1:], expected["peak"], atol=2e-6)
|
||||
np.testing.assert_allclose(payloads["sar3d_peak"]["table"][0], expected["level"], rtol=2e-3)
|
||||
|
||||
def test_errors_are_published(self):
|
||||
for malformed in ("empty", "capture", "frequency", "axis", "angle"):
|
||||
with self.subTest(malformed=malformed):
|
||||
request = fixture()
|
||||
if malformed == "empty":
|
||||
request["latest"]["traces"] = []
|
||||
elif malformed == "capture":
|
||||
request["latest"]["traces"][0]["end"] = request["latest"]["traces"][0]["start"]
|
||||
elif malformed == "frequency":
|
||||
request["latest"]["traces"][0]["frequency"].reverse()
|
||||
elif malformed == "axis":
|
||||
request["voxel"] = 1e-8
|
||||
else:
|
||||
request["pitch"] = 90
|
||||
result = native(request)
|
||||
payloads = {p["name"]: p for p in result["payloads"]}
|
||||
status = payloads["sar3d_status"]["table"]
|
||||
self.assertEqual(payloads["sar3d_status"]["columns"], len(status))
|
||||
self.assertEqual(status[0], 2)
|
||||
self.assertTrue(bytes(map(int, status[1:])).decode())
|
||||
|
||||
def test_empty_signal_is_reported(self):
|
||||
request = fixture()
|
||||
request.update(previous=[], radius=0)
|
||||
for trace in request["latest"]["traces"]:
|
||||
trace["s21"] = [[0, 0]] * len(trace["frequency"])
|
||||
payloads = {p["name"]: p for p in native(request)["payloads"]}
|
||||
self.assertEqual(payloads["sar3d_status"]["table"][0], 1)
|
||||
self.assertNotIn("sar3d_peak", payloads)
|
||||
for name in ("sar3d_front", "sar3d_side", "sar3d_top", "sar3d_oblique"):
|
||||
self.assertTrue(np.all(np.asarray(payloads[name]["image"]) == 0))
|
||||
|
||||
def test_invalid_fft_dimensions_are_rejected(self):
|
||||
for points, first_bin, size in ((0, 0, 8), (1, 0, 7), (1, 8, 8), (2, 7, 8), (1, 0, 2**31)):
|
||||
with self.subTest(points=points, first_bin=first_bin, size=size):
|
||||
result = native(dict(operation="ifft", bands=[[1, 0]],
|
||||
points=points, first_bin=first_bin, size=size))
|
||||
self.assertIn("error", result)
|
||||
|
||||
def test_weak_signal_stays_finite(self):
|
||||
request = fixture()
|
||||
request.update(previous=[], radius=0, normalize=False)
|
||||
for trace in request["latest"]["traces"]:
|
||||
trace["s21"] = (np.asarray(trace["s21"]) * 1e-38).tolist()
|
||||
payloads = self.ready(native(request))
|
||||
for name in ("sar3d_front", "sar3d_side", "sar3d_top", "sar3d_oblique"):
|
||||
values = np.asarray(payloads[name]["image"])
|
||||
self.assertTrue(np.isfinite(values).all())
|
||||
self.assertGreater(values.max(), 0)
|
||||
self.assertLess(values.max(), 1e-10)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user