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
This commit is contained in:
lamm45 2023-06-19 08:11:01 -04:00 committed by GitHub
parent 67fd2d02ef
commit ddcbd2cc26
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 59 additions and 14 deletions

View File

@ -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<uchar> _buf(std::max(m*2*sizeof(float) + (m*3+1)*sizeof(int), n*2*sizeof(float)));
cv::AutoBuffer<uchar> _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;

View File

@ -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<uchar>(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<float>(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