This commit is contained in:
Ayzen
2026-06-13 12:07:23 +03:00
parent f0d095de80
commit e7f2d25585
20 changed files with 903 additions and 148 deletions
@@ -319,6 +319,38 @@ void normalize_in_place(std::vector<double>& values) {
}
}
// Run `body(row_begin, row_end)` over a partition of [0, row_count) across the
// available hardware threads. Each call owns a disjoint, contiguous row range, so
// a body that writes only its own rows needs no synchronization. The calling
// thread runs the first chunk while spawned workers handle the rest. Falls back to
// a single serial call when there is one row or no concurrency is reported.
template <typename Body>
void parallel_for_rows(std::size_t row_count, const Body& body) {
if (row_count == 0U) {
return;
}
const unsigned int detected = std::thread::hardware_concurrency();
const std::size_t worker_count = std::clamp<std::size_t>(
detected == 0U ? 1U : static_cast<std::size_t>(detected), 1U, row_count
);
if (worker_count == 1U) {
body(0U, row_count);
return;
}
const std::size_t chunk = (row_count + worker_count - 1U) / worker_count;
std::vector<std::thread> workers;
workers.reserve(worker_count - 1U);
for (std::size_t begin = chunk; begin < row_count; begin += chunk) {
workers.emplace_back(body, begin, std::min(begin + chunk, row_count));
}
body(0U, std::min(chunk, row_count));
for (auto& worker : workers) {
worker.join();
}
}
[[nodiscard]] auto build_gaussian_kernel(double sigma) -> std::vector<double> {
if (!(sigma > 0.0)) {
return {1.0};
@@ -355,27 +387,36 @@ void normalize_in_place(std::vector<double>& values) {
std::vector<double> temp(values.size(), 0.0);
std::vector<double> output(values.size(), 0.0);
for (std::size_t row = 0U; row < height; ++row) {
for (std::size_t col = 0U; col < width; ++col) {
double sum = 0.0;
for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
const auto sample_col = reflect_index(static_cast<std::ptrdiff_t>(col) + offset, width);
sum += values[(row * width) + sample_col] * kernel[static_cast<std::size_t>(offset + radius)];
// Separable convolution: every output row depends only on its own row (pass 1)
// or only on already-finished `temp` (pass 2), so each pass parallelizes over
// disjoint rows. parallel_for_rows joins between passes — that join is the
// barrier guaranteeing `temp` is complete before pass 2 reads it. Per-cell
// arithmetic is unchanged, so the result is identical to a serial sweep.
parallel_for_rows(height, [&](std::size_t row_begin, std::size_t row_end) {
for (std::size_t row = row_begin; row < row_end; ++row) {
for (std::size_t col = 0U; col < width; ++col) {
double sum = 0.0;
for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
const auto sample_col = reflect_index(static_cast<std::ptrdiff_t>(col) + offset, width);
sum += values[(row * width) + sample_col] * kernel[static_cast<std::size_t>(offset + radius)];
}
temp[(row * width) + col] = sum;
}
temp[(row * width) + col] = sum;
}
}
});
for (std::size_t row = 0U; row < height; ++row) {
for (std::size_t col = 0U; col < width; ++col) {
double sum = 0.0;
for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
const auto sample_row = reflect_index(static_cast<std::ptrdiff_t>(row) + offset, height);
sum += temp[(sample_row * width) + col] * kernel[static_cast<std::size_t>(offset + radius)];
parallel_for_rows(height, [&](std::size_t row_begin, std::size_t row_end) {
for (std::size_t row = row_begin; row < row_end; ++row) {
for (std::size_t col = 0U; col < width; ++col) {
double sum = 0.0;
for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
const auto sample_row = reflect_index(static_cast<std::ptrdiff_t>(row) + offset, height);
sum += temp[(sample_row * width) + col] * kernel[static_cast<std::size_t>(offset + radius)];
}
output[(row * width) + col] = sum;
}
output[(row * width) + col] = sum;
}
}
});
return output;
}
@@ -901,6 +942,32 @@ void normalize_pair_ascans(
return grid;
}
// Compact signature that fully determines the imaging grid: the antenna layout
// plus the depth window and imaging plane. Two collections with an equal signature
// produce an identical grid, so the (expensive) grid build can be memoized.
[[nodiscard]] auto grid_signature(
const GeometrySelection& selection,
double max_depth_m,
double min_z_m,
double imaging_plane_y_m
) -> std::vector<double> {
std::vector<double> signature;
signature.reserve((selection.x_tx.size() + selection.x_rx.size()) * 3U + 3U);
const auto append = [&](const std::vector<double>& axis) {
signature.insert(signature.end(), axis.begin(), axis.end());
};
append(selection.x_tx);
append(selection.y_tx);
append(selection.z_tx);
append(selection.x_rx);
append(selection.y_rx);
append(selection.z_rx);
signature.push_back(max_depth_m);
signature.push_back(min_z_m);
signature.push_back(imaging_plane_y_m);
return signature;
}
[[nodiscard]] auto interpolate_complex(const AscanResult& ascan, double tau_s) -> std::complex<double> {
if (ascan.samples.empty() || !(ascan.dt_s > 0.0) || tau_s < 0.0) {
return {0.0, 0.0};
@@ -982,38 +1049,6 @@ void normalize_pair_ascans(
return std::clamp(range_weight * angle_weight, 0.0, kTotalWeightMax);
}
// Run `body(row_begin, row_end)` over a partition of [0, row_count) across the
// available hardware threads. Each call owns a disjoint, contiguous row range, so
// a body that writes only its own rows needs no synchronization. The calling
// thread runs the first chunk while spawned workers handle the rest. Falls back to
// a single serial call when there is one row or no concurrency is reported.
template <typename Body>
void parallel_for_rows(std::size_t row_count, const Body& body) {
if (row_count == 0U) {
return;
}
const unsigned int detected = std::thread::hardware_concurrency();
const std::size_t worker_count = std::clamp<std::size_t>(
detected == 0U ? 1U : static_cast<std::size_t>(detected), 1U, row_count
);
if (worker_count == 1U) {
body(0U, row_count);
return;
}
const std::size_t chunk = (row_count + worker_count - 1U) / worker_count;
std::vector<std::thread> workers;
workers.reserve(worker_count - 1U);
for (std::size_t begin = chunk; begin < row_count; begin += chunk) {
workers.emplace_back(body, begin, std::min(begin + chunk, row_count));
}
body(0U, std::min(chunk, row_count));
for (auto& worker : workers) {
worker.join();
}
}
[[nodiscard]] auto backproject_coherent(
const std::vector<SelectedTrace>& selected_traces,
const std::unordered_map<PairKey, AscanResult>& ascans_by_pair,
@@ -1857,7 +1892,20 @@ void add_bp_score_metrics(
normalize_pair_ascans(ascans_by_pair, velocity_mps, min_depth_m, max_depth_m);
const double imaging_plane_y_m = static_cast<double>(live_config.gpr_imaging_plane_y_m);
const auto grid = build_grid(selection, max_depth_m, kGridZMinM, imaging_plane_y_m);
// The imaging grid (90k cells x 6 distance fields) only depends on antenna
// geometry and the depth window, which are constant across a run, so it is
// memoized between collections and rebuilt only when that signature changes.
// process_backprojection_gpr is entered by a single processing thread (the
// backprojection workers join before it returns), so a local static is safe.
static std::vector<double> cached_grid_signature;
static GridDefinition cached_grid;
const auto signature = grid_signature(selection, max_depth_m, kGridZMinM, imaging_plane_y_m);
if (cached_grid.x_grid.empty() || signature != cached_grid_signature) {
cached_grid = build_grid(selection, max_depth_m, kGridZMinM, imaging_plane_y_m);
cached_grid_signature = signature;
}
const GridDefinition& grid = cached_grid;
if (grid.x_grid.empty() || grid.z_grid.empty()) {
return results;
}