From 2dd3d1371a9975248b1c236f4a94295a7ac5777d Mon Sep 17 00:00:00 2001 From: pratham-mcw Date: Mon, 16 Mar 2026 16:21:03 +0530 Subject: [PATCH] Merge pull request #28636 from pratham-mcw:imgproc_distance_transform_opt imgproc: add SIMD support for distance transform function #28636 - This PR adds OpenCV SIMD intrinsics-based optimizations to the distance transform functions for improved performance. - The optimized implementation uses vectorized operations to accelerate the forward and backward passes of distance computation. - In x64 architecture, distance transform function benefit from IPP-based optimized implementations. However, on ARM64 platforms, the execution falls back to scalar implementation, which results in lower performance. - After introducing these changes, the distance transform functions showed noticeable performance improvements on Windows-ARM64. image - [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 --- modules/imgproc/src/distransform.cpp | 486 +++++++++++++++++++++------ 1 file changed, 384 insertions(+), 102 deletions(-) diff --git a/modules/imgproc/src/distransform.cpp b/modules/imgproc/src/distransform.cpp index 66be0e8ea5..3885bc5320 100644 --- a/modules/imgproc/src/distransform.cpp +++ b/modules/imgproc/src/distransform.cpp @@ -70,7 +70,7 @@ static void distanceTransform_3x3( const Mat& _src, Mat& _temp, Mat& _dist, const float* metrics ) { const int BORDER = 1; - int i, j; + const unsigned int HV_DIST = CV_FLT_TO_FIX( metrics[0], DIST_SHIFT ); const unsigned int DIAG_DIST = CV_FLT_TO_FIX( metrics[1], DIST_SHIFT ); const unsigned int DIST_MAX = UINT_MAX - DIAG_DIST; @@ -79,65 +79,190 @@ distanceTransform_3x3( const Mat& _src, Mat& _temp, Mat& _dist, const float* met const uchar* src = _src.ptr(); int* temp = _temp.ptr(); float* dist = _dist.ptr(_dist.rows - 1); - int srcstep = (int)(_src.step/sizeof(src[0])); - int step = (int)(_temp.step/sizeof(temp[0])); - int dststep = (int)(_dist.step/sizeof(dist[0])); - Size size = _src.size(); + + const int srcstep = (int)_src.step; + const int step = (int)(_temp.step / sizeof(int)); + const int dststep = (int)(_dist.step / sizeof(float)); + + const int width = _src.cols; + const int height = _src.rows; initTopBottom( _temp, BORDER, DIST_MAX ); // forward pass - unsigned int* tmp = (unsigned int*)(temp + BORDER*step) + BORDER; - const uchar* s = src; - for( i = 0; i < size.height; i++ ) - { - for( j = 0; j < BORDER; j++ ) - tmp[-j-1] = tmp[size.width + j] = DIST_MAX; + unsigned int* row = (unsigned int*)(temp + BORDER * step) + BORDER; + const uchar* srow = src; - for( j = 0; j < size.width; j++ ) + for (int i = 0; i < height; i++) + { + unsigned int* cur = row; + const uchar* s = srow; + const unsigned int* top = (const unsigned int*)(cur - step); + + // set horizontal border sentinels + cur[-1] = DIST_MAX; + cur[width] = DIST_MAX; + + unsigned int left = cur[-1]; + +#if defined(CV_SIMD) + const int BLOCK = 8; + int j = 0; + const v_uint32 v_diag = vx_setall(DIAG_DIST); + const v_uint32 v_hv = vx_setall(HV_DIST); + + for (; j <= width - BLOCK; j += BLOCK) { - if( !s[j] ) - tmp[j] = 0; - else + v_uint32 tl0 = vx_load(top + j - 1); + v_uint32 t0 = vx_load(top + j); + v_uint32 tr0 = vx_load(top + j + 1); + v_uint32 tl1 = vx_load(top + j + 3); + v_uint32 t1 = vx_load(top + j + 4); + v_uint32 tr1 = vx_load(top + j + 5); + tl0 = v_add(tl0, v_diag); + t0 = v_add(t0, v_hv); + tr0 = v_add(tr0, v_diag); + tl1 = v_add(tl1, v_diag); + t1 = v_add(t1, v_hv); + tr1 = v_add(tr1, v_diag); + v_uint32 m0 = v_min(tl0, t0); + m0 = v_min(m0, tr0); + v_uint32 m1 = v_min(tl1, t1); + m1 = v_min(m1, tr1); + v_store(cur + j, m0); + v_store(cur + j + 4, m1); + + for (int k = 0; k < BLOCK; k++) { - unsigned int t0 = tmp[j-step-1] + DIAG_DIST; - unsigned int t = tmp[j-step] + HV_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step+1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-1] + HV_DIST; - if( t0 > t ) t0 = t; - tmp[j] = (t0 > DIST_MAX) ? DIST_MAX : t0; + int idx = j + k; + if (!s[idx]) + { + cur[idx] = 0; + left = 0; + } + else + { + unsigned int m = cur[idx]; + unsigned int lm = left + HV_DIST; + if (lm < m) m = lm; + if (m > DIST_MAX) m = DIST_MAX; + cur[idx] = m; + left = m; + } } } - tmp += step; - s += srcstep; + + for (; j < width; j++) + { + if (!s[j]) + { + cur[j] = 0; + left = 0; + } + else + { + unsigned int t0 = top[j - 1] + DIAG_DIST; + unsigned int t1 = top[j] + HV_DIST; + unsigned int t2 = top[j + 1] + DIAG_DIST; + unsigned int t3 = left + HV_DIST; + + unsigned int m = t0 < t1 ? t0 : t1; + m = m < t2 ? m : t2; + m = m < t3 ? m : t3; + + cur[j] = m > DIST_MAX ? DIST_MAX : m; + left = cur[j]; + } + } +#else + for (int j = 0; j < width; j++) + { + if (!s[j]) + { + cur[j] = 0; + left = 0; + } + else + { + unsigned int t0 = top[j - 1] + DIAG_DIST; + unsigned int t1 = top[j] + HV_DIST; + unsigned int t2 = top[j + 1] + DIAG_DIST; + unsigned int t3 = left + HV_DIST; + + unsigned int m = t0 < t1 ? t0 : t1; + m = m < t2 ? m : t2; + m = m < t3 ? m : t3; + + cur[j] = m > DIST_MAX ? DIST_MAX : m; + left = cur[j]; + } + } +#endif + + row += step; + srow += srcstep; } // backward pass - float* d = (float*)dist; - for( i = size.height - 1; i >= 0; i-- ) + float* drow = dist; + unsigned int* cur = row - step; +#if defined(CV_SIMD) + const v_float32 scale_v = vx_setall(scale); +#endif + for (int i = height - 1; i >= 0; i--) { - tmp -= step; + unsigned int* bottom = cur + step; - for( j = size.width - 1; j >= 0; j-- ) + unsigned int right = cur[width]; + for (int j = width - 1; j >= 0; j--) { - unsigned int t0 = tmp[j]; - if( t0 > HV_DIST ) + unsigned int t0 = cur[j]; + + if (t0 > HV_DIST) { - unsigned int t = tmp[j+step+1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step] + HV_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step-1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+1] + HV_DIST; - if( t0 > t ) t0 = t; - tmp[j] = t0; + unsigned int b1 = bottom[j + 1] + DIAG_DIST; + unsigned int b2 = bottom[j] + HV_DIST; + unsigned int b3 = bottom[j - 1] + DIAG_DIST; + unsigned int r = right + HV_DIST; + + unsigned int t = b1 < b2 ? b1 : b2; + t = t < b3 ? t : b3; + t = t < r ? t : r; + + if (t < t0) + { + t0 = t; + cur[j] = t0; + } } - d[j] = (float)(t0 * scale); + + right = cur[j]; } - d -= dststep; + +#if defined(CV_SIMD) + int j = 0; + const uint32_t* tptr = (const uint32_t*)cur; + float* dptr = drow; + + for (; j <= width - 8; j += 8) + { + v_uint32 v0 = vx_load(tptr + j); + v_uint32 v1 = vx_load(tptr + j + 4); + v_float32 f0 = v_mul(v_cvt_f32(v0), scale_v); + v_float32 f1 = v_mul(v_cvt_f32(v1), scale_v); + v_store(dptr + j, f0); + v_store(dptr + j + 4, f1); + } + + for (; j < width; j++) + dptr[j] = (float)(tptr[j] * scale); +#else + for (int j = 0; j < width; j++) + drow[j] = (float)(cur[j] * scale); +#endif + + cur -= step; + drow -= dststep; } } @@ -146,91 +271,248 @@ static void distanceTransform_5x5( const Mat& _src, Mat& _temp, Mat& _dist, const float* metrics ) { const int BORDER = 2; - int i, j; + const unsigned int HV_DIST = CV_FLT_TO_FIX( metrics[0], DIST_SHIFT ); const unsigned int DIAG_DIST = CV_FLT_TO_FIX( metrics[1], DIST_SHIFT ); const unsigned int LONG_DIST = CV_FLT_TO_FIX( metrics[2], DIST_SHIFT ); const unsigned int DIST_MAX = UINT_MAX - LONG_DIST; const float scale = 1.f/(1 << DIST_SHIFT); - const uchar* src = _src.ptr(); + const uchar* src = _src.ptr(); int* temp = _temp.ptr(); float* dist = _dist.ptr(_dist.rows - 1); - int srcstep = (int)(_src.step/sizeof(src[0])); - int step = (int)(_temp.step/sizeof(temp[0])); - int dststep = (int)(_dist.step/sizeof(dist[0])); - Size size = _src.size(); - initTopBottom( _temp, BORDER, DIST_MAX ); + const int srcstep = (int)_src.step; + const int step = (int)(_temp.step / sizeof(int)); + const int dststep = (int)(_dist.step / sizeof(float)); + + const int width = _src.cols; + const int height = _src.rows; + + initTopBottom(_temp, BORDER, DIST_MAX); // forward pass - unsigned int* tmp = (unsigned int*)(temp + BORDER*step) + BORDER; - const uchar* s = src; - for( i = 0; i < size.height; i++ ) - { - for( j = 0; j < BORDER; j++ ) - tmp[-j-1] = tmp[size.width + j] = DIST_MAX; + unsigned int* row = (unsigned int*)(temp + BORDER * step) + BORDER; + const uchar* srow = src; - for( j = 0; j < size.width; j++ ) + for (int i = 0; i < height; i++) + { + unsigned int* cur = row; + const uchar* s = srow; + const unsigned int* top1 = cur - step; + const unsigned int* top2 = cur - step * 2; + cur[-1] = cur[-2] = DIST_MAX; + cur[width] = cur[width + 1] = DIST_MAX; +#if defined(CV_SIMD) + unsigned int left1 = DIST_MAX; + unsigned int left2 = DIST_MAX; + const int BLOCK = 8; + int j = 0; + const v_uint32 v_hv = vx_setall(HV_DIST); + const v_uint32 v_diag = vx_setall(DIAG_DIST); + const v_uint32 v_long = vx_setall(LONG_DIST); + const v_uint32 v_max = vx_setall(DIST_MAX); + + for (; j <= width - BLOCK; j += BLOCK) { - if( !s[j] ) - tmp[j] = 0; - else + v_uint32 t2m1_0 = vx_load(top2 + j - 1); + v_uint32 t2p1_0 = vx_load(top2 + j + 1); + v_uint32 t1m2_0 = vx_load(top1 + j - 2); + v_uint32 t1m1_0 = vx_load(top1 + j - 1); + v_uint32 t1_0 = vx_load(top1 + j); + v_uint32 t1p1_0 = vx_load(top1 + j + 1); + v_uint32 t1p2_0 = vx_load(top1 + j + 2); + + t2m1_0 = v_add(t2m1_0, v_long); + t2p1_0 = v_add(t2p1_0, v_long); + t1m2_0 = v_add(t1m2_0, v_long); + t1m1_0 = v_add(t1m1_0, v_diag); + t1_0 = v_add(t1_0, v_hv); + t1p1_0 = v_add(t1p1_0, v_diag); + t1p2_0 = v_add(t1p2_0, v_long); + + v_uint32 m0 = v_min(t2m1_0, t2p1_0); + v_uint32 m0_1 = v_min(t1m2_0, t1m1_0); + v_uint32 m0_2 = v_min(t1p1_0, t1p2_0); + m0 = v_min(m0, m0_1); + m0_2 = v_min(m0_2, t1_0); + m0 = v_min(m0, m0_2); + m0 = v_min(m0, v_max); + + v_store(cur + j, m0); + + v_uint32 t2m1_1 = vx_load(top2 + j + 3); + v_uint32 t2p1_1 = vx_load(top2 + j + 5); + v_uint32 t1m2_1 = vx_load(top1 + j + 2); + v_uint32 t1m1_1 = vx_load(top1 + j + 3); + v_uint32 t1_1 = vx_load(top1 + j + 4); + v_uint32 t1p1_1 = vx_load(top1 + j + 5); + v_uint32 t1p2_1 = vx_load(top1 + j + 6); + + t2m1_1 = v_add(t2m1_1, v_long); + t2p1_1 = v_add(t2p1_1, v_long); + t1m2_1 = v_add(t1m2_1, v_long); + t1m1_1 = v_add(t1m1_1, v_diag); + t1_1 = v_add(t1_1, v_hv); + t1p1_1 = v_add(t1p1_1, v_diag); + t1p2_1 = v_add(t1p2_1, v_long); + + v_uint32 m1 = v_min(t2m1_1, t2p1_1); + v_uint32 m1_1 = v_min(t1m2_1, t1m1_1); + v_uint32 m1_2 = v_min(t1p1_1, t1p2_1); + m1 = v_min(m1, m1_1); + m1_2 = v_min(m1_2, t1_1); + m1 = v_min(m1, m1_2); + m1 = v_min(m1, v_max); + + v_store(cur + j + 4, m1); + + for (int k = 0; k < BLOCK; k++) { - unsigned int t0 = tmp[j-step*2-1] + LONG_DIST; - unsigned int t = tmp[j-step*2+1] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step-2] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step-1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step] + HV_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step+1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-step+2] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j-1] + HV_DIST; - if( t0 > t ) t0 = t; - tmp[j] = (t0 > DIST_MAX) ? DIST_MAX : t0; + int idx = j + k; + if (!s[idx]) + { + cur[idx] = 0; + left1 = 0; + left2 = 0; + } + else + { + unsigned int m = cur[idx]; + unsigned int lm = left1 + HV_DIST; + if (lm < m) m = lm; + if (m > DIST_MAX) m = DIST_MAX; + cur[idx] = m; + left2 = left1; + left1 = m; + } } } - tmp += step; - s += srcstep; + + for (; j < width; j++) + { + if (!s[j]) + { + cur[j] = 0; + left1 = left2 = 0; + } + else + { + unsigned int t0 = top2[j - 1] + LONG_DIST; + unsigned int t = top2[j + 1] + LONG_DIST; + if (t < t0) t0 = t; + t = top1[j - 2] + LONG_DIST; if (t < t0) t0 = t; + t = top1[j - 1] + DIAG_DIST; if (t < t0) t0 = t; + t = top1[j] + HV_DIST; if (t < t0) t0 = t; + t = top1[j + 1] + DIAG_DIST; if (t < t0) t0 = t; + t = top1[j + 2] + LONG_DIST; if (t < t0) t0 = t; + t = left1 + HV_DIST; if (t < t0) t0 = t; + + cur[j] = t0 > DIST_MAX ? DIST_MAX : t0; + left2 = left1; + left1 = cur[j]; + } + } +#else + unsigned int left = DIST_MAX; + for (int j = 0; j < width; j++) + { + if (!s[j]) + { + cur[j] = 0; + left = 0; + } + else + { + unsigned int t0 = top2[j - 1] + LONG_DIST; + unsigned int t = top2[j + 1] + LONG_DIST; + if (t < t0) t0 = t; + t = top1[j - 2] + LONG_DIST; + if (t < t0) t0 = t; + t = top1[j - 1] + DIAG_DIST; + if (t < t0) t0 = t; + t = top1[j] + HV_DIST; + if (t < t0) t0 = t; + t = top1[j + 1] + DIAG_DIST; + if (t < t0) t0 = t; + t = top1[j + 2] + LONG_DIST; + if (t < t0) t0 = t; + t = left + HV_DIST; + if (t < t0) t0 = t; + cur[j] = t0 > DIST_MAX ? DIST_MAX : t0; + left = cur[j]; + } + } +#endif + row += step; + srow += srcstep; } // backward pass - float* d = (float*)dist; - for( i = size.height - 1; i >= 0; i-- ) - { - tmp -= step; + float* drow = dist; + unsigned int* cur = row - step; - for( j = size.width - 1; j >= 0; j-- ) +#if defined(CV_SIMD) + const v_float32 scale_v = vx_setall(scale); +#endif + + for (int i = height - 1; i >= 0; i--) + { + unsigned int* bot1 = cur + step; + unsigned int* bot2 = cur + step * 2; + unsigned int right = DIST_MAX; + + for (int j = width - 1; j >= 0; j--) { - unsigned int t0 = tmp[j]; + unsigned int t0 = cur[j]; if( t0 > HV_DIST ) { - unsigned int t = tmp[j+step*2+1] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step*2-1] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step+2] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step+1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step] + HV_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step-1] + DIAG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+step-2] + LONG_DIST; - if( t0 > t ) t0 = t; - t = tmp[j+1] + HV_DIST; - if( t0 > t ) t0 = t; - tmp[j] = t0; + unsigned int t = bot2[j + 1] + LONG_DIST; + if (t < t0) t0 = t; + t = bot2[j - 1] + LONG_DIST; + if (t < t0) t0 = t; + t = bot1[j + 2] + LONG_DIST; + if (t < t0) t0 = t; + t = bot1[j + 1] + DIAG_DIST; + if (t < t0) t0 = t; + t = bot1[j] + HV_DIST; + if (t < t0) t0 = t; + t = bot1[j - 1] + DIAG_DIST; + if (t < t0) t0 = t; + t = bot1[j - 2] + LONG_DIST; + if (t < t0) t0 = t; + t = right + HV_DIST; + if (t < t0) t0 = t; + + cur[j] = t0; } - d[j] = (float)(t0 * scale); + right = cur[j]; } - d -= dststep; + +#if defined(CV_SIMD) + int j = 0; + const uint32_t* tptr = (const uint32_t*)cur; + float* dptr = drow; + + for (; j <= width - 8; j += 8) + { + v_uint32 v0 = vx_load(tptr + j); + v_uint32 v1 = vx_load(tptr + j + 4); + v_float32 f0 = v_mul(v_cvt_f32(v0), scale_v); + v_float32 f1 = v_mul(v_cvt_f32(v1), scale_v); + v_store(dptr + j, f0); + v_store(dptr + j + 4, f1); + } + + for (; j < width; j++) + dptr[j] = (float)(tptr[j] * scale); +#else + for (int j = 0; j < width; j++) + drow[j] = (float)(cur[j] * scale); +#endif + + cur -= step; + drow -= dststep; } }