From ffa38e1b74654a78c17f65a992adeb60953b1d53 Mon Sep 17 00:00:00 2001 From: Ding zhehao <2915288232@qq.com> Date: Wed, 10 Jun 2026 23:23:35 +0800 Subject: [PATCH] Merge pull request #29194 from Enilrats:opt-emd-performance imgproc: optimize EMD (Earth Mover's Distance) solver performance #29194 ### imgproc: optimize EMD solver performance using O(V) spanning tree traversal --- This PR significantly optimizes the performance of the Earth Mover's Distance (`cv::EMD`) solver in `modules/imgproc/src/emd_new.cpp`. Specifically, it refactors the dual-variable calculation in `EMDSolver::findBasicVars()` from a naive $O((N+M)^2)$ linked-list scanning approach to an optimized $O(N+M)$ BFS tree traversal utilizing existing adjacency lists, leading to a massive speedup. --- #### Technical Details & Core Bottleneck Fixed 1. **Algorithmic Complexity Reduction in `findBasicVars()`:** - **Before**: The original implementation solved the dual variables $u_i$ and $v_j$ by traversing the entire unmarked rows (`u0_head`) or columns (`v0_head`) linked-lists and invoking `getIsX(i, j)` inside nested loops to find connected basic variables. This resulted in an $O((N+M)^2)$ complexity per simplex iteration. For a scale of $2000 \times 2000$, this performed up to $16,000,000$ operations per iteration. - **After**: Since `EMDSolver` already maintains the adjacency lists of the basic variables tree (`rows_x` and `cols_x`), we can traverse the spanning tree in linear time. This PR implements a dual-queue BFS tree traversal. The complexity per iteration is drastically reduced to $O(N+M)$, performing at most $4000$ operations per iteration. 2. **Cache Locality & Pointer-Chasing Elimination:** - Replaced pointer-chasing on dynamically-allocated linked lists with contiguous, stack-allocated array queues (`cv::AutoBuffer`), significantly improving CPU L1/L2 cache hit rates and enabling hardware prefetching. --- #### Performance Benchmarks Below is the benchmark comparison evaluated on a standard CPU. #### Test 1: dims = 64 | Scale ($N, M$) | Original EMD (ms) | Optimized EMD (This PR) | Speedup | | :--- | :--- | :--- | :--- | | **100** | 5.473 | 3.419 | **1.60x** | | **500** | 381.871 | 244.355 | **1.56x** | | **1000** | 1893.053 | 1369.016 | **1.38x** | | **2000** | 11387.792 | 8331.221 | **1.37x** | #### Test 2: dims = 3 | Scale ($N, M$) | Original EMD (ms) | Optimized EMD (This PR) | Speedup | | :--- | :--- | :--- | :--- | | **100** | 4.433 | 3.042 | **1.46x** | | **500** | 365.762 | 259.735 | **1.41x** | | **1000** | 1989.400 | 1421.952 | **1.40x** | | **2000** | 12731.836 | 7952.210 | **1.60x** | *(Note: The exact performance figures may vary slightly depending on the compiler and test machine.)* --- #### Verification - All existing tests in `opencv_test_imgproc` (including EMD tests) pass successfully. No regressions were introduced. --- ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [ ] There is a reference to the original bug report and related work - [ ] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [ ] The feature is well documented and sample code can be built with the project CMake --- modules/imgproc/perf/perf_emd.cpp | 37 ++++++++ modules/imgproc/src/emd_new.cpp | 143 +++++++++++------------------- 2 files changed, 88 insertions(+), 92 deletions(-) create mode 100644 modules/imgproc/perf/perf_emd.cpp diff --git a/modules/imgproc/perf/perf_emd.cpp b/modules/imgproc/perf/perf_emd.cpp new file mode 100644 index 0000000000..1952a58867 --- /dev/null +++ b/modules/imgproc/perf/perf_emd.cpp @@ -0,0 +1,37 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html + +#include "perf_precomp.hpp" +#include + +namespace opencv_test { namespace { + +typedef tuple EMD_Size_Dim_t; +typedef perf::TestBaseWithParam EMD_Fixture; + +PERF_TEST_P(EMD_Fixture, L1_Distance, testing::Combine( + testing::Values(100, 500, 1000), + testing::Values(3, 64) +)) +{ + int size = get<0>(GetParam()); + int dims = get<1>(GetParam()); + + Mat sign1(size, dims + 1, CV_32FC1); + Mat sign2(size, dims + 1, CV_32FC1); + + theRNG().fill(sign1, RNG::UNIFORM, 0.1, 1.0); + theRNG().fill(sign2, RNG::UNIFORM, 0.1, 1.0); + + declare.in(sign1, sign2); + + TEST_CYCLE() + { + cv::EMD(sign1, sign2, cv::DIST_L1); + } + + SANITY_CHECK_NOTHING(); +} + +}} \ No newline at end of file diff --git a/modules/imgproc/src/emd_new.cpp b/modules/imgproc/src/emd_new.cpp index 59a6b24342..81528d01b0 100644 --- a/modules/imgproc/src/emd_new.cpp +++ b/modules/imgproc/src/emd_new.cpp @@ -442,119 +442,78 @@ double EMDSolver::calcFlow(Mat* flow_) const int EMDSolver::findBasicVars() const { - int i, j; - int u_cfound, v_cfound; - Node1D u0_head, u1_head, *cur_u, *prev_u; - Node1D v0_head, v1_head, *cur_v, *prev_v; - bool found; - CV_Assert(u != 0 && v != 0); - /* initialize the rows list (u) and the columns list (v) */ - u0_head.next = u; - for (i = 0; i < ssize; i++) - { - u[i].next = u + i + 1; - } - u[ssize - 1].next = 0; - u1_head.next = 0; + // 1. Initialize status flags using contiguous memory to eliminate pointer chasing + AutoBuffer computed_buf(ssize + dsize); + char* row_computed = computed_buf.data(); + char* col_computed = computed_buf.data() + ssize; + memset(row_computed, 0, ssize + dsize); - v0_head.next = ssize > 1 ? v + 1 : 0; - for (i = 1; i < dsize; i++) - { - v[i].next = v + i + 1; - } - v[dsize - 1].next = 0; - v1_head.next = 0; + // 2. Create BFS queues + AutoBuffer queue_buf(ssize + dsize); + int* row_queue = queue_buf.data(); + int* col_queue = queue_buf.data() + ssize; - /* there are ssize+dsize variables but only ssize+dsize-1 independent equations, - so set v[0]=0 */ + int row_head = 0, row_tail = 0; + int col_head = 0, col_tail = 0; + + // Initial condition: enqueue column 0 as the root node and set its value to 0 v[0].val = 0; - v1_head.next = v; - v1_head.next->next = 0; + col_computed[0] = true; + col_queue[col_tail++] = 0; - /* loop until all variables are found */ - u_cfound = v_cfound = 0; - while (u_cfound < ssize || v_cfound < dsize) + int u_cfound = 0; + int v_cfound = 1; + + // 3. Dual-queue interactive BFS traversal over the spanning tree (Time Complexity: O(N + M)) + while (row_head < row_tail || col_head < col_tail) { - found = false; - if (v_cfound < dsize) + // Process currently marked columns to update their connected rows + while (col_head < col_tail) { - /* loop over all marked columns */ - prev_v = &v1_head; - cur_v = v1_head.next; - found = found || (cur_v != 0); - for (; cur_v != 0; cur_v = cur_v->next) - { - float cur_v_val = cur_v->val; + int j = col_queue[col_head++]; + float cur_v_val = v[j].val; - j = (int)(cur_v - v); - /* find the variables in column j */ - prev_u = &u0_head; - for (cur_u = u0_head.next; cur_u != 0;) + // Use adjacency list cols_x to directly access rows connected to column j, avoiding full scans + for (Node2D* xp = cols_x[j]; xp != 0; xp = xp->next[1]) + { + int i = xp->i; + if (!row_computed[i]) { - i = (int)(cur_u - u); - if (getIsX(i, j)) - { - /* compute u[i] */ - cur_u->val = getCost(i, j) - cur_v_val; - /* ...and add it to the marked list */ - prev_u->next = cur_u->next; - cur_u->next = u1_head.next; - u1_head.next = cur_u; - cur_u = prev_u->next; - } - else - { - prev_u = cur_u; - cur_u = cur_u->next; - } + u[i].val = getCost(i, j) - cur_v_val; + row_computed[i] = true; + row_queue[row_tail++] = i; // Enqueue the newly resolved row + u_cfound++; } - prev_v->next = cur_v->next; - v_cfound++; } } - if (u_cfound < ssize) + // Process currently marked rows to update their connected columns + while (row_head < row_tail) { - /* loop over all marked rows */ - prev_u = &u1_head; - cur_u = u1_head.next; - found = found || (cur_u != 0); - for (; cur_u != 0; cur_u = cur_u->next) + int i = row_queue[row_head++]; + float cur_u_val = u[i].val; + + // Use adjacency list rows_x to directly access columns connected to row i + for (Node2D* xp = rows_x[i]; xp != 0; xp = xp->next[0]) { - float cur_u_val = cur_u->val; - i = (int)(cur_u - u); - /* find the variables in rows i */ - prev_v = &v0_head; - for (cur_v = v0_head.next; cur_v != 0;) + int j = xp->j; + if (!col_computed[j]) { - j = (int)(cur_v - v); - if (getIsX(i, j)) - { - /* compute v[j] */ - cur_v->val = getCost(i, j) - cur_u_val; - /* ...and add it to the marked list */ - prev_v->next = cur_v->next; - cur_v->next = v1_head.next; - v1_head.next = cur_v; - cur_v = prev_v->next; - } - else - { - prev_v = cur_v; - cur_v = cur_v->next; - } + v[j].val = getCost(i, j) - cur_u_val; + col_computed[j] = true; + col_queue[col_tail++] = j; // Enqueue the newly resolved column + v_cfound++; } - prev_u->next = cur_u->next; - u_cfound++; } } - - if (!found) - return -1; } + // If the number of traversed nodes is insufficient, the graph is disconnected and the spanning tree is incomplete + if (u_cfound < ssize || v_cfound < dsize) + return -1; + return 0; } @@ -1008,4 +967,4 @@ float cv::wrapperEMD(InputArray _sign1, OutputArray _flow) { return EMD(_sign1, _sign2, distType, _cost, lowerBound.get(), _flow); -} +} \ No newline at end of file