310 lines
17 KiB
Python
310 lines
17 KiB
Python
"""Generate SAR3D/GPR references from notebook cells 1/14 and Horns coherent BP.
|
|
|
|
Requires the original, untracked Multi_antenna_demonstrate_depthnorm.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_backprojection_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 reference_functions(source_path: Path, parameters: dict, cells: tuple[int, ...] = (1, 14)) -> tuple[dict, dict]:
|
|
contents = source_path.read_text(encoding="utf-8")
|
|
sources = [(str(i), "".join(json.loads(contents)["cells"][i]["source"])) for i in cells] if cells else [("script", contents)]
|
|
namespace = {"__name__": __name__}
|
|
hashes = {}
|
|
for index, source in sources:
|
|
hashes[str(index)] = hashlib.sha256(source.encode()).hexdigest()
|
|
tree = ast.parse(source, filename=f"{source_path.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"{source_path.name}:cell{index}",
|
|
"exec", dont_inherit=True), namespace)
|
|
return namespace, hashes
|
|
|
|
|
|
def scene(f_start=1e9, f_stop=4e9, points=201) -> 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=f_start, f_stop=f_stop,
|
|
min_depth=1., max_depth=3., speed=.4, pitch=13., ground=.3,
|
|
voxel=.4, radius=2, power=1.5, compensation=True, normalize=True,
|
|
equalize=True, subtract_mean=False, mean_count=1, motion="int_minus",
|
|
)
|
|
# 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(f_start, f_stop, points, 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 = []
|
|
members = [*request["previous"], request["latest"]]
|
|
for index, member in enumerate(members):
|
|
background = previous_mean({**request, "previous": members[:index]})
|
|
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"])
|
|
samples = values[:, 0] + 1j * values[:, 1] - background.get((trace["output"], trace["input"]), 0.)
|
|
np.save(snapshot / f"{stem}_s21.npy", samples)
|
|
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[-1]
|
|
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"],
|
|
PAIR_DEPTH_EQUALIZE=request["equalize"],
|
|
SAR_SCAN_ROOT=directory, DATA_PATH=center, BG_PATH=directory, v=3e8,
|
|
)
|
|
ns, hashes = reference_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_ant"] = antenna_x
|
|
ns["x_min"], ns["x_max"] = antenna_x.min() - 2, antenna_x.max() + 2
|
|
ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], 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)
|
|
grid = ns["sar_build_full_region_grid"](config, y_range_m=(request.get("y_min", -.5), request.get("y_max", 2.)))
|
|
for axis in ("x", "z"):
|
|
size = request.get(f"{axis}_size", 0.)
|
|
if size:
|
|
center_axis = .5 * (grid[axis][0] + grid[axis][-1])
|
|
grid[axis] = np.linspace(center_axis - size/2, center_axis + size/2, max(round(size/request["voxel"])+1, 2))
|
|
grid["XX"], grid["ZZ"] = np.meshgrid(grid["x"], grid["z"])
|
|
ns["x_min"], ns["x_max"] = grid["x"][0], grid["x"][-1]
|
|
result = ns["sar_compute_volume"](center, None, config, timestamps, grid=grid, verbose=False)
|
|
volume, grid = result["volume"], result["grid"]
|
|
iy, iz, ix = ns["sar_global_max_index"](volume, grid)
|
|
from scipy.signal import peak_widths
|
|
widths = []
|
|
for axis, peak_index, profile in (("x", ix, volume[iy, iz, :]), ("y", iy, volume[:, iz, ix]), ("z", iz, volume[iy, :, ix])):
|
|
profile = np.asarray(profile, dtype=np.float64)
|
|
width, *_ = peak_widths(profile, [peak_index], prominence_data=(np.array([profile[peak_index]]),
|
|
np.array([0], dtype=np.intp), np.array([len(profile)-1], dtype=np.intp)))
|
|
widths.append(max(float(width[0]), 1.) * (grid[axis][1]-grid[axis][0]))
|
|
expected = {axis: grid[axis] for axis in "xyz"}
|
|
expected["target_width"] = np.asarray(widths)
|
|
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 previous_mean(request: dict) -> dict:
|
|
samples = {}
|
|
if request["subtract_mean"] and request["mean_count"]:
|
|
for member in request["previous"][-request["mean_count"]:]:
|
|
for trace in member["traces"]:
|
|
values = np.asarray(trace["s21"])
|
|
samples.setdefault((trace["output"], trace["input"]), []).append(values[:, 0] + 1j * values[:, 1])
|
|
return {key: np.mean(values, axis=0) for key, values in samples.items()}
|
|
|
|
|
|
def gpr_reference(notebook: Path, script: Path, request: dict) -> tuple[np.ndarray, dict]:
|
|
"""Execute Horns coherent BP unchanged, adding the notebook's equalizer before its scalar norm."""
|
|
parameters = dict(F_START=request["f_start"], F_STOP=request["f_stop"], MIN_DEPTH=request["min_depth"],
|
|
MAX_DEPTH=request["max_depth"], LOOK_ANGLE_DEG=request["pitch"],
|
|
BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"], v=3e8)
|
|
ns, hashes = reference_functions(script, parameters, ())
|
|
ns["N_tx"], ns["N_rx"] = len(request["tx"]), len(request["rx"])
|
|
ns["OUTPUT_IDX"], ns["INPUT_IDX"] = [int(row[0]) for row in request["tx"]], [int(row[0]) for row in request["rx"]]
|
|
for kind in ("tx", "rx"):
|
|
for axis, values in zip("xyz", np.asarray(request[kind])[:, 1:].T):
|
|
ns[f"{axis}_{kind}"] = values
|
|
xs = np.concatenate((ns["x_tx"], ns["x_rx"]))
|
|
ns["x_ant"] = xs
|
|
ns["x_min"], ns["x_max"] = xs.min()-2, xs.max()+2
|
|
ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], ns["NX_BP"])
|
|
ns["z_grid_bp"] = np.linspace(.1, request["max_depth"], ns["NZ_BP"])
|
|
ns["XX_bp"], ns["ZZ_bp"] = np.meshgrid(ns["x_grid_bp"], ns["z_grid_bp"])
|
|
ns["depth_gate"] = (ns["ZZ_bp"] >= request["min_depth"]) & (ns["ZZ_bp"] <= request["max_depth"])
|
|
member = request["latest"]
|
|
duration = (max(t["end"] for t in member["traces"]) - min(t["start"] for t in member["traces"])) * 1e-9
|
|
config = ns["MotionConfig"](speed_m_s=request["speed"], look_angle_deg=request["pitch"],
|
|
tx_per_segment=len(request["tx"]), tx_order_full=tuple(ns["OUTPUT_IDX"]),
|
|
frame_acq_time_s=duration, n_freq_hint=len(member["traces"][0]["frequency"]),
|
|
intra_phase_focus_only=request["motion"] == "int_focus", motion_phase_sign=-1.,
|
|
pair_order_phys=[(t["output"], t["input"]) for t in member["traces"]])
|
|
clocks, _ = ns["compute_pair_timestamps"](config, ns["INPUT_IDX"], ns["OUTPUT_IDX"])
|
|
mean = previous_mean(request)
|
|
H, T = {}, {}
|
|
for trace in member["traces"]:
|
|
key = (trace["output"], trace["input"])
|
|
values = np.asarray(trace["s21"], dtype=np.float32)
|
|
samples = values[:, 0].astype(float) + 1j * values[:, 1].astype(float)
|
|
samples -= mean.get(key, 0.)
|
|
frequency = np.asarray(trace["frequency"], dtype=np.float32)
|
|
samples, _ = ns["apply_intra_sweep_phase_correction"](samples, frequency, clocks[key], config,
|
|
ns["v"], ns["F_START"], ns["F_STOP"])
|
|
t, _, h, _, _ = ns["compute_ascan_bp"](samples, frequency, ns["F_START"], ns["F_STOP"])
|
|
H[key], T[key] = h, t
|
|
if request["equalize"]:
|
|
eq, notebook_hashes = reference_functions(notebook, {**parameters, "PAIR_DEPTH_EQUALIZE": True})
|
|
hashes.update({f"notebook_cell_{cell}": value for cell, value in notebook_hashes.items()})
|
|
eq.update(x_ant=xs, x_min=ns["x_min"], x_max=ns["x_max"])
|
|
H = eq["equalize_pairs_by_depth"](H, T)
|
|
for key, h in H.items():
|
|
depth = T[key] * ns["v"] / 2
|
|
mask = (depth >= request["min_depth"]) & (depth <= request["max_depth"])
|
|
scale = np.percentile(np.abs(h)[mask] if mask.any() else np.abs(h), 50.)
|
|
H[key] = h / ((scale if np.isfinite(scale) and scale > 1e-15 else 1.) + 1e-15)
|
|
raw, _ = ns["backproject_coherent"](H, T)
|
|
return ns["normalize_bp_map"](raw).astype(np.float32), hashes
|
|
|
|
|
|
def generate_gpr(notebook: Path, script: Path, output: Path) -> None:
|
|
cases = [dict(equalize=False), dict(equalize=True),
|
|
dict(equalize=True, subtract_mean=True, mean_count=1, speed=-.2, pitch=-7., motion="int_focus", power=2.)]
|
|
arrays = {}
|
|
for i, settings in enumerate(cases):
|
|
request = scene(1.01e9, 4.01e9, 151) if i == 1 else scene()
|
|
request.update(settings, mode="gpr")
|
|
arrays[f"case_{i}_image"], hashes = gpr_reference(notebook, script, request)
|
|
arrays[f"case_{i}_request"] = np.array(json.dumps(request))
|
|
arrays["provenance"] = np.array(json.dumps(dict(script=script.name, source_sha256=hashes,
|
|
notebook=notebook.name, numpy=np.__version__, scipy=scipy.__version__)))
|
|
np.savez_compressed(output, **arrays)
|
|
print(f"Wrote {output}: {output.stat().st_size} bytes, {len(cases)} GPR references")
|
|
|
|
|
|
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, equalize=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),
|
|
dict(compensation=True, normalize=True, subtract_mean=True, mean_count=1),
|
|
dict(compensation=True, normalize=True, x_size=3., z_size=2., y_min=-.2, y_max=1.4),
|
|
]
|
|
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=[1, 14], cases=len(cases),
|
|
numpy=np.__version__, scipy=scipy.__version__,
|
|
method="Unchanged notebook AST functions; per-frame NumPy previous-mean subtraction; 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_demonstrate_depthnorm.ipynb")
|
|
parser.add_argument("--output", type=Path, default=REFERENCE)
|
|
parser.add_argument("--gpr-script", type=Path, default=ROOT / "Horns_motion_4x4_ground.py")
|
|
parser.add_argument("--gpr-output", type=Path, default=REFERENCE.with_name("gpr_reference.npz"))
|
|
args = parser.parse_args()
|
|
generate(args.notebook, args.output)
|
|
generate_gpr(args.notebook, args.gpr_script, args.gpr_output)
|