Skip to content

subset indexing deterministic across openMP threads. no early termination of multiwindow neighbour list generation. #372

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JoelPhys merged 4 commits into from
Aug 28, 2026
Merged

subset indexing deterministic across openMP threads. no early termination of multiwindow neighbour list generation. #372

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
18bd9b2
small bugfix. inactive subsets are skipped for MULTIWINDOW method.
JoelPhys Jul 28, 2026
c548074
cleared pointers for partial multiwindow.
JoelPhys Aug 28, 2026
d36356e
small var name changes in dicmultiwindow_util to be consistent with o…
JoelPhys Aug 28, 2026
c20976b
Make subset indexing deterministic across OpenMP threads.
JoelPhys Aug 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions src/pyvale/dic/cpp/dicmain.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,9 +191,9 @@ void engine(const py::array_t<bool>& img_roi_arr,
results_ref_l = results_def_l;
interp_ref_l = make_interp(conf.interp_routine, conf.fullpaths[img_num_ref_l]);

bool* roi_updated = propagate_roi(img_roi, results_def_l, conf, ss_grid_l);
std::unique_ptr<bool[]> roi_updated(propagate_roi(img_roi, results_def_l, conf, ss_grid_l));
multiwindow_l.clear();
multiwindow_init_partial(multiwindow_l, roi_updated, conf, mwconf, saveconf,
multiwindow_init_partial(multiwindow_l, roi_updated.get(), conf, mwconf, saveconf,
mwconf.overlap.size() - 1);

WindowLevel last_level;
Expand Down Expand Up @@ -354,9 +354,9 @@ void engine(const py::array_t<bool>& img_roi_arr,
results_ref_r = results_def_r;


bool* roi_updated = propagate_roi(img_roi, results_def_l, conf, ss_grid_l);
std::unique_ptr<bool[]> roi_updated(propagate_roi(img_roi, results_def_l, conf, ss_grid_l));
multiwindow_l.clear();
multiwindow_init_partial(multiwindow_l, roi_updated, conf, mwconf, saveconf,
multiwindow_init_partial(multiwindow_l, roi_updated.get(), conf, mwconf, saveconf,
mwconf.overlap.size() - 1);

WindowLevel last_level;
Expand Down
8 changes: 8 additions & 0 deletions src/pyvale/dic/cpp/dicmultiwindow_only.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,14 @@ void multiwindow_only(const Interpolator &interp_ref,

// append fourier results to master result vectors
OptResult res(conf.num_params);

// Keep inactive subsets at their reset/default result state, matching
// multiwindow_rg and calc_rigid_displacements.
if (!ss_grid.active_ss[ss]) {
results_def.append(res, ss);
continue;
}

res.u = multiwindow.back().u[ss];
res.p[0] = multiwindow.back().u[ss];
res.v = multiwindow.back().v[ss];
Expand Down
43 changes: 23 additions & 20 deletions src/pyvale/dic/cpp/dicmultiwindow_util.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
#include <algorithm>
#include <omp.h>
#include <csignal>
#include <cstdio>
#include <stdexcept>

// Common Header files
#include "../../common_cpp/progressbar.hpp"
Expand Down Expand Up @@ -63,7 +65,7 @@ void multiwindow_init_partial(std::vector<WindowLevel> &level,

for (size_t lvl = 0; lvl < num_levels; lvl++) {

const bool is_last = (lvl == num_levels - 1);
const bool is_last = (lvl == mwconf.overlap.size() - 1);
const subset::Grid *prev = (lvl > 0) ? &level[lvl-1].layout : nullptr;

level.emplace_back(img_roi,
Expand All @@ -87,19 +89,19 @@ void WindowLevel::gen_neighlist(const subset::Grid &layout_prev) {

// a list containing the number of neighbours from the previous
// window size for each subset in the current window size
num_neigh_list.resize(layout.num);
num_neigh_list.assign(layout.num, 0);

// we know the neigh_list is going to be a max size of
// max_neigh*num_ss. we can resize this later once populated
neigh_list.resize(max_num_neigh*layout.num);
neigh_list.assign(max_num_neigh*layout.num, -1);

// shared error state
std::atomic<bool> failed(false);

struct ErrorInfo {
int ss = -1;
int ss_x = 0;
int ss_y = 0;
double cx = 0.0;
double cy = 0.0;
size_t num_found = 0;
};

Expand All @@ -115,16 +117,16 @@ void WindowLevel::gen_neighlist(const subset::Grid &layout_prev) {

if (!layout.active_ss[ss]) continue;

// corner of subset
const int ss_x = layout.coords[2*ss];
const int ss_y = layout.coords[2*ss+1];
// centre of subset
const double cx = layout.coords[2*ss];
const double cy = layout.coords[2*ss+1];

// Vector to store pairs of (distance, index)
std::vector<std::pair<double, int>> dist_index_list;

// loop over a 10x10 section from the previous window
int idx_x = (ss_x / prev_step);
int idx_y = (ss_y / prev_step);
const int idx_x = static_cast<int>(std::floor(cx / prev_step));
const int idx_y = static_cast<int>(std::floor(cy / prev_step));

// range of neighbour search
int min_x = std::max(0,idx_x-5);
Expand All @@ -139,11 +141,11 @@ void WindowLevel::gen_neighlist(const subset::Grid &layout_prev) {
int nss_idx = layout_prev.mask[y*layout_prev.num_ss_x+x];
if (nss_idx == -1) continue;

int nss_x = layout_prev.coords[2*nss_idx];
int nss_y = layout_prev.coords[2*nss_idx+1];
const double nss_x = layout_prev.coords[2*nss_idx];
const double nss_y = layout_prev.coords[2*nss_idx+1];

double dx = (nss_x) - ss_x;
double dy = (nss_y) - ss_y;
double dx = (nss_x) - cx;
double dy = (nss_y) - cy;
double dist_sq = dx*dx + dy*dy;

dist_index_list.emplace_back(dist_sq, nss_idx);
Expand All @@ -159,8 +161,8 @@ void WindowLevel::gen_neighlist(const subset::Grid &layout_prev) {
// only first thread records error
if (failed.compare_exchange_strong(expected, true)) {
error.ss = ss;
error.ss_x = ss_x;
error.ss_y = ss_y;
error.cx = cx;
error.cy = cy;
error.num_found = dist_index_list.size();
}
continue;
Expand Down Expand Up @@ -189,10 +191,10 @@ void WindowLevel::gen_neighlist(const subset::Grid &layout_prev) {
// snprintf(msg,
// sizeof(msg),
// "Could not find any neighbours from the previous FFT "
// "window size for subset (%d, %d). "
// "window size for subset (%.3f, %.3f). "
// "Found %zu neighbours.",
// error.ss_x,
// error.ss_y,
// error.cx,
// error.cy,
// error.num_found);
//
// throw std::runtime_error(msg);
Expand Down Expand Up @@ -275,7 +277,7 @@ void WindowLevel::calc_rigid_displacements(const WindowLevel &prev,
#ifdef _MSC_VER
#pragma omp parallel
#else
#pragma omp parallel shared(stop_request, level, prev, interp_def, img_num_ref, search_area, u, v, max_val)
#pragma omp parallel shared(stop_request, level, prev, interp_ref, interp_def, img_num_ref, search_area, u, v, max_val)
#endif
{
if (fft_precision == util::FFTPrecision::FLOAT32) {
Expand Down Expand Up @@ -322,6 +324,7 @@ void WindowLevel::calc_rigid_displacements(const WindowLevel &prev,
fout << "\n";

for (int ss = 0; ss < layout.num; ss++){
if (!layout.active_ss[ss]) continue;
fout << layout.coords[2*ss] << saveconf.delimiter;
fout << layout.coords[2*ss+1] << saveconf.delimiter;
fout << u[ss] << saveconf.delimiter;
Expand Down
88 changes: 21 additions & 67 deletions src/pyvale/dic/cpp/dicsubset.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,10 +214,11 @@ namespace subset {
ss_grid.coords.resize(2*ss_grid.num_in_mask, -1);


// temp array for storing subset coords for each thread
std::vector<int> thread_counts(omp_get_max_threads(), 0);
// Store validity by grid location so subset indices are independent of
// OpenMP scheduling and thread count.
std::vector<unsigned char> valid_grid(ss_grid.num_in_mask, 0);

// First pass: count valid subsets per thread
// First pass: determine which grid locations contain valid subsets.
#pragma omp parallel for collapse(2)
for (int j = 0; j < num_ss_y; j++) {
for (int i = 0; i < num_ss_x; i++) {
Expand Down Expand Up @@ -258,83 +259,36 @@ namespace subset {
valid = (valid_count >= (ss_size_x * ss_size_y) * 0.70);
}

if (valid) {
int tid = omp_get_thread_num();
thread_counts[tid]++;
}
valid_grid[j * num_ss_x + i] = static_cast<unsigned char>(valid);
}
}

// Compute prefix sum to get offsets
std::vector<int> thread_offsets(omp_get_max_threads(), 0);
for (int t = 1; t < thread_offsets.size(); t++)
thread_offsets[t] = thread_offsets[t-1] + thread_counts[t-1];
// Compute deterministic row-major subset indices.
int total_valid = 0;
for (int grid_idx = 0; grid_idx < ss_grid.num_in_mask; ++grid_idx) {
if (valid_grid[grid_idx]) {
ss_grid.mask[grid_idx] = total_valid++;
}
}

int total_valid = thread_offsets.back() + thread_counts.back();
ss_grid.coords.resize(2 * total_valid);
ss_grid.num = total_valid;
ss_grid.active_ss.resize(total_valid, true);
ss_grid.active_total = total_valid;

// Reset thread counts to use as writing indices
std::fill(thread_counts.begin(), thread_counts.end(), 0);

// Populate coordinates using the deterministic indices in mask.
#pragma omp parallel for collapse(2)
for (int j = 0; j < num_ss_y; j++) {
for (int i = 0; i < num_ss_x; i++) {

// calculate the coordinates of the subset
const int ss_x = i * ss_step;
const int ss_y = j * ss_step;

// pixel range of subset
const int xmin = ss_x;
const int ymin = ss_y;
const int xmax = ss_x + ss_size_x-1;
const int ymax = ss_y + ss_size_y-1;

// check if subset is within image and ROI.
bool valid = true;
int valid_count = 0;

for (int px_y = ymin; px_y <= ymax && valid; px_y++) {
for (int px_x = xmin; px_x <= xmax && valid; px_x++) {

// When no partial subset filling all px must be within roi
if (!partial) {
if (!px_in_img_dims(px_x, px_y, px_hori, px_vert) ||
!px_in_roi(px_x, px_y, px_hori, px_vert, img_roi)) {
valid = false;
break;
}
}

// When partial count num of px in roi. if its outside
// the image its still not valid
else {
if (!px_in_img_dims(px_x, px_y, px_hori, px_vert)) {
valid = false;
break;
}
if (px_in_roi(px_x, px_y, px_hori, px_vert, img_roi)) valid_count++;
}
}

if (!valid && !partial) break;
}

if (partial && valid) {
valid = (valid_count >= (ss_size_x * ss_size_y) * 0.70);
}

// if its a valid subset. add it to a list of coordinates
if (valid) {
const int tid = omp_get_thread_num();
const int offset = thread_offsets[tid] + thread_counts[tid];
ss_grid.coords[2*offset] = ss_x + static_cast<double>(ss_size_x)/2-0.5;
ss_grid.coords[2*offset + 1] = ss_y + static_cast<double>(ss_size_y)/2-0.5;
ss_grid.mask[j * num_ss_x + i] = offset;
thread_counts[tid]++;
const int offset = ss_grid.mask[j * num_ss_x + i];
if (offset != -1) {
const int ss_x = i * ss_step;
const int ss_y = j * ss_step;
ss_grid.coords[2*offset] =
ss_x + static_cast<double>(ss_size_x)/2 - 0.5;
ss_grid.coords[2*offset + 1] =
ss_y + static_cast<double>(ss_size_y)/2 - 0.5;
}
}
}
Expand Down
Loading