From ddcbd2cc268f6293bb636a41c1aa19dfb3c181ad Mon Sep 17 00:00:00 2001 From: lamm45 <96844552+lamm45@users.noreply.github.com> Date: Mon, 19 Jun 2023 08:11:01 -0400 Subject: [PATCH] Merge pull request #22798 from lamm45:distransform-large Fix distransform to work with large images #22798 This attempts to fix the following bug which was caused by storing squares of large integers into 32-bit floating point variables: https://github.com/opencv/opencv/issues/22732 ### 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 - [x] There is a reference to the original bug report and related work - [x] 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/src/distransform.cpp | 31 +++++++------- .../imgproc/test/test_distancetransform.cpp | 42 +++++++++++++++++++ 2 files changed, 59 insertions(+), 14 deletions(-) diff --git a/modules/imgproc/src/distransform.cpp b/modules/imgproc/src/distransform.cpp index a2c5a0d5fd..dd22d4f50b 100644 --- a/modules/imgproc/src/distransform.cpp +++ b/modules/imgproc/src/distransform.cpp @@ -448,7 +448,7 @@ static void getDistanceTransformMask( int maskType, float *metrics ) struct DTColumnInvoker : ParallelLoopBody { - DTColumnInvoker( const Mat* _src, Mat* _dst, const int* _sat_tab, const float* _sqr_tab) + DTColumnInvoker( const Mat* _src, Mat* _dst, const int* _sat_tab, const int* _sqr_tab) { src = _src; dst = _dst; @@ -481,7 +481,7 @@ struct DTColumnInvoker : ParallelLoopBody { dist = dist + 1 - sat_tab[dist - d[j]]; d[j] = dist; - dptr[0] = sqr_tab[dist]; + dptr[0] = (float)sqr_tab[dist]; } } } @@ -489,12 +489,12 @@ struct DTColumnInvoker : ParallelLoopBody const Mat* src; Mat* dst; const int* sat_tab; - const float* sqr_tab; + const int* sqr_tab; }; struct DTRowInvoker : ParallelLoopBody { - DTRowInvoker( Mat* _dst, const float* _sqr_tab, const float* _inv_tab ) + DTRowInvoker( Mat* _dst, const int* _sqr_tab, const float* _inv_tab ) { dst = _dst; sqr_tab = _sqr_tab; @@ -529,7 +529,7 @@ struct DTRowInvoker : ParallelLoopBody for(;;k--) { p = v[k]; - float s = (fq + sqr_tab[q] - d[p] - sqr_tab[p])*inv_tab[q - p]; + float s = (fq - d[p] + (sqr_tab[q]-sqr_tab[p]))*inv_tab[q - p]; if( s > z[k] ) { k++; @@ -552,28 +552,28 @@ struct DTRowInvoker : ParallelLoopBody } Mat* dst; - const float* sqr_tab; + const int* sqr_tab; const float* inv_tab; }; static void trueDistTrans( const Mat& src, Mat& dst ) { - const float inf = 1e15f; + const int inf = INT_MAX; CV_Assert( src.size() == dst.size() ); CV_Assert( src.type() == CV_8UC1 && dst.type() == CV_32FC1 ); int i, m = src.rows, n = src.cols; - cv::AutoBuffer _buf(std::max(m*2*sizeof(float) + (m*3+1)*sizeof(int), n*2*sizeof(float))); + cv::AutoBuffer _buf(std::max(m*2*sizeof(int) + (m*3+1)*sizeof(int), n*2*sizeof(float))); // stage 1: compute 1d distance transform of each column - float* sqr_tab = (float*)_buf.data(); + int* sqr_tab = (int*)_buf.data(); int* sat_tab = cv::alignPtr((int*)(sqr_tab + m*2), sizeof(int)); int shift = m*2; for( i = 0; i < m; i++ ) - sqr_tab[i] = (float)(i*i); + sqr_tab[i] = i*i; for( i = m; i < m*2; i++ ) sqr_tab[i] = inf; for( i = 0; i < shift; i++ ) @@ -584,13 +584,14 @@ trueDistTrans( const Mat& src, Mat& dst ) cv::parallel_for_(cv::Range(0, n), cv::DTColumnInvoker(&src, &dst, sat_tab, sqr_tab), src.total()/(double)(1<<16)); // stage 2: compute modified distance transform for each row - float* inv_tab = sqr_tab + n; + float* inv_tab = (float*)sqr_tab + n; - inv_tab[0] = sqr_tab[0] = 0.f; + inv_tab[0] = 0.f; + sqr_tab[0] = 0; for( i = 1; i < n; i++ ) { inv_tab[i] = (float)(0.5/i); - sqr_tab[i] = (float)(i*i); + sqr_tab[i] = i*i; } cv::parallel_for_(cv::Range(0, m), cv::DTRowInvoker(&dst, sqr_tab, inv_tab)); @@ -752,7 +753,9 @@ void cv::distanceTransform( InputArray _src, OutputArray _dst, OutputArray _labe CV_IPP_CHECK() { #if IPP_DISABLE_PERF_TRUE_DIST_MT - if(cv::getNumThreads()<=1 || (src.total()<(int)(1<<14))) + // IPP uses floats, but 4097 cannot be squared into a float + if((cv::getNumThreads()<=1 || (src.total()<(int)(1<<14))) && + src.rows < 4097 && src.cols < 4097) #endif { IppStatus status; diff --git a/modules/imgproc/test/test_distancetransform.cpp b/modules/imgproc/test/test_distancetransform.cpp index 32fc7aa2c3..742595631a 100644 --- a/modules/imgproc/test/test_distancetransform.cpp +++ b/modules/imgproc/test/test_distancetransform.cpp @@ -302,4 +302,46 @@ BIGDATA_TEST(Imgproc_DistanceTransform, large_image_12218) EXPECT_EQ(nz, (size.height*size.width / 2)); } +TEST(Imgproc_DistanceTransform, wide_image_22732) +{ + Mat src = Mat::zeros(1, 4099, CV_8U); // 4099 or larger used to be bad + Mat dist(src.rows, src.cols, CV_32F); + distanceTransform(src, dist, DIST_L2, DIST_MASK_PRECISE, CV_32F); + int nz = countNonZero(dist); + EXPECT_EQ(nz, 0); +} + +TEST(Imgproc_DistanceTransform, large_square_22732) +{ + Mat src = Mat::zeros(8000, 8005, CV_8U), dist; + distanceTransform(src, dist, DIST_L2, DIST_MASK_PRECISE, CV_32F); + int nz = countNonZero(dist); + EXPECT_EQ(dist.size(), src.size()); + EXPECT_EQ(dist.type(), CV_32F); + EXPECT_EQ(nz, 0); + + Point p0(src.cols-1, src.rows-1); + src.setTo(1); + src.at(p0) = 0; + distanceTransform(src, dist, DIST_L2, DIST_MASK_PRECISE, CV_32F); + EXPECT_EQ(dist.size(), src.size()); + EXPECT_EQ(dist.type(), CV_32F); + bool first = true; + int nerrs = 0; + for (int y = 0; y < dist.rows; y++) + for (int x = 0; x < dist.cols; x++) { + float d = dist.at(y, x); + double dx = (double)(x - p0.x), dy = (double)(y - p0.y); + float d0 = (float)sqrt(dx*dx + dy*dy); + if (std::abs(d0 - d) > 1) { + if (first) { + printf("y=%d, x=%d. dist_ref=%.2f, dist=%.2f\n", y, x, d0, d); + first = false; + } + nerrs++; + } + } + EXPECT_EQ(0, nerrs) << "reference distance map is different from computed one at " << nerrs << " pixels\n"; +} + }} // namespace