From 18bd9b26e1179a4ae18675b7f271d445e0e76d4a Mon Sep 17 00:00:00 2001 From: JoelPhys Date: Tue, 28 Jul 2026 14:27:30 +0100 Subject: [PATCH 1/4] small bugfix. inactive subsets are skipped for MULTIWINDOW method. --- src/pyvale/dic/cpp/dicmultiwindow_only.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/pyvale/dic/cpp/dicmultiwindow_only.cpp b/src/pyvale/dic/cpp/dicmultiwindow_only.cpp index c1192788..05651ede 100644 --- a/src/pyvale/dic/cpp/dicmultiwindow_only.cpp +++ b/src/pyvale/dic/cpp/dicmultiwindow_only.cpp @@ -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]; From c548074727a5e43d71d7776f11261acb7c98c892 Mon Sep 17 00:00:00 2001 From: JoelPhys Date: Fri, 28 Aug 2026 08:53:39 +0100 Subject: [PATCH 2/4] cleared pointers for partial multiwindow. --- src/pyvale/dic/cpp/dicmain.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pyvale/dic/cpp/dicmain.cpp b/src/pyvale/dic/cpp/dicmain.cpp index 08b9c4d1..484c8cb4 100644 --- a/src/pyvale/dic/cpp/dicmain.cpp +++ b/src/pyvale/dic/cpp/dicmain.cpp @@ -191,9 +191,9 @@ void engine(const py::array_t& 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 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; @@ -354,9 +354,9 @@ void engine(const py::array_t& 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 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; From d36356e3a851da4f02b3c58dc68c49542054bd27 Mon Sep 17 00:00:00 2001 From: JoelPhys Date: Fri, 28 Aug 2026 08:54:21 +0100 Subject: [PATCH 3/4] small var name changes in dicmultiwindow_util to be consistent with other parts of code. --- src/pyvale/dic/cpp/dicmultiwindow_util.cpp | 43 ++++++++++++---------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/pyvale/dic/cpp/dicmultiwindow_util.cpp b/src/pyvale/dic/cpp/dicmultiwindow_util.cpp index 20d91f86..b584c129 100644 --- a/src/pyvale/dic/cpp/dicmultiwindow_util.cpp +++ b/src/pyvale/dic/cpp/dicmultiwindow_util.cpp @@ -14,6 +14,8 @@ #include #include #include +#include +#include // Common Header files #include "../../common_cpp/progressbar.hpp" @@ -63,7 +65,7 @@ void multiwindow_init_partial(std::vector &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, @@ -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 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; }; @@ -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> 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(std::floor(cx / prev_step)); + const int idx_y = static_cast(std::floor(cy / prev_step)); // range of neighbour search int min_x = std::max(0,idx_x-5); @@ -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); @@ -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; @@ -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); @@ -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) { @@ -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; From c20976b21ef93a8b8620bc65a7a28a51066d7930 Mon Sep 17 00:00:00 2001 From: JoelPhys Date: Fri, 28 Aug 2026 08:56:15 +0100 Subject: [PATCH 4/4] Make subset indexing deterministic across OpenMP threads. --- src/pyvale/dic/cpp/dicsubset.cpp | 88 ++++++++------------------------ 1 file changed, 21 insertions(+), 67 deletions(-) diff --git a/src/pyvale/dic/cpp/dicsubset.cpp b/src/pyvale/dic/cpp/dicsubset.cpp index 9ba4b823..7758a4b4 100644 --- a/src/pyvale/dic/cpp/dicsubset.cpp +++ b/src/pyvale/dic/cpp/dicsubset.cpp @@ -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 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 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++) { @@ -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(valid); } } - // Compute prefix sum to get offsets - std::vector 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(ss_size_x)/2-0.5; - ss_grid.coords[2*offset + 1] = ss_y + static_cast(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(ss_size_x)/2 - 0.5; + ss_grid.coords[2*offset + 1] = + ss_y + static_cast(ss_size_y)/2 - 0.5; } } }