From 592f8d8c1bc0ad249de5467aa0d1323acb2714c1 Mon Sep 17 00:00:00 2001 From: Tom Becker Date: Thu, 28 Dec 2017 06:23:11 -0800 Subject: [PATCH] Merge pull request #10232 from TomBecker-BD:hough-many-circles Hough many circles (#10232) * Add Hui's optimization. Merge with latest changes in OpenCV. * Use conditional compilation instead of a runtime flag. * Whitespace. * Create the sequence for the nonzero edge pixels only if using that approach. * Improve performance for finding very large numbers of circles * Return the circles with the larger accumulator values first, as per API documentation. Use a separate step to check distance between circles. Allows circles to be sorted by strength first. Avoids locking in EstimateRadius which was slowing it down. Return centers only if maxRadius == 0 as per API documentation. * Sort the circles so results are deterministic. Otherwise the order of circles with the same strength depends on parallel processing completion order. * Add test for HoughCircles. * Add beads test. * Wrap the non-zero points structure in a common interface so the code can use either a vector or a matrix. * Remove the special case for skipping the radius search if maxRadius==0. * Add performance tests. * Use NULL instead of nullptr. OpenCV should compile with C++98 compiler. * Put test suite name first. Use different test suite names for each test to avoid an error from the test runner. * Address build bot errors and warnings. * Skip radius search if maxRadius < 0. * Dynamically switch to NZPointList when it will be faster than NZPointSet. * Fix compile error: missing 'typename' prior to dependent type name. * Fix compile error: missing 'typename' prior to dependent type name. This time fix it the non C++ 11 way. * Fix compile error: no type named 'const_reference' in 'class cv::NZPointList' * Disable ManySmallCircles tests. Failing on Mac. * Change beads image to JPEG for smaller file size. Try enabling the ManySmallCircles tests again. * Remove ManySmallCircles tests. They are failing on the Mac build. * Fix expectations to check all circles. * Changing case on a case-insensitive file system Step 1: remove the old file names * Changing case on a case-insensitive file system Step 2: add them back with the new names * Fix cmpAccum function to be strictly weak ordered. * Add tests for many small circles. * imgproc(perf): fix HoughCircles tests * imgproc(houghCircles): refactor code - simplify NZPointList - drop broken (de-synchronization of 'current'/'mi' fields) NZPointSet iterator - NZPointSet iterator is replaced to direct area scan - use SIMD intrinsics - avoid std exceptions (build for embedded systems) --- modules/imgproc/include/opencv2/imgproc.hpp | 7 +- modules/imgproc/perf/perf_houghcircles.cpp | 57 ++ modules/imgproc/src/hough.cpp | 639 +++++++++++--------- modules/imgproc/test/test_houghcircles.cpp | 259 ++++++++ 4 files changed, 681 insertions(+), 281 deletions(-) create mode 100644 modules/imgproc/perf/perf_houghcircles.cpp create mode 100644 modules/imgproc/test/test_houghcircles.cpp diff --git a/modules/imgproc/include/opencv2/imgproc.hpp b/modules/imgproc/include/opencv2/imgproc.hpp index 109b547fd5..0d1ba52f58 100644 --- a/modules/imgproc/include/opencv2/imgproc.hpp +++ b/modules/imgproc/include/opencv2/imgproc.hpp @@ -2094,8 +2094,8 @@ Example: : @note Usually the function detects the centers of circles well. However, it may fail to find correct radii. You can assist to the function by specifying the radius range ( minRadius and maxRadius ) if -you know it. Or, you may set maxRadius to 0 to return centers only without radius search, and find the correct -radius using an additional procedure. +you know it. Or, you may set maxRadius to a negative number to return centers only without radius +search, and find the correct radius using an additional procedure. @param image 8-bit, single-channel, grayscale input image. @param circles Output vector of found circles. Each vector is encoded as a 3-element @@ -2114,7 +2114,8 @@ accumulator threshold for the circle centers at the detection stage. The smaller false circles may be detected. Circles, corresponding to the larger accumulator values, will be returned first. @param minRadius Minimum circle radius. -@param maxRadius Maximum circle radius. +@param maxRadius Maximum circle radius. If <= 0, uses the maximum image dimension. If < 0, returns +centers without finding the radius. @sa fitEllipse, minEnclosingCircle */ diff --git a/modules/imgproc/perf/perf_houghcircles.cpp b/modules/imgproc/perf/perf_houghcircles.cpp new file mode 100644 index 0000000000..3f57343c6e --- /dev/null +++ b/modules/imgproc/perf/perf_houghcircles.cpp @@ -0,0 +1,57 @@ +#include "perf_precomp.hpp" +#include "opencv2/imgproc.hpp" +#include "opencv2/imgproc/types_c.h" + +using namespace std; +using namespace cv; +using namespace perf; + +PERF_TEST(PerfHoughCircles, Basic) +{ + string filename = getDataPath("cv/imgproc/stuff.jpg"); + const double dp = 1.0; + double minDist = 20; + double edgeThreshold = 20; + double accumThreshold = 30; + int minRadius = 20; + int maxRadius = 200; + + Mat img = imread(filename, IMREAD_GRAYSCALE); + ASSERT_FALSE(img.empty()) << "Unable to load source image " << filename; + + GaussianBlur(img, img, Size(9, 9), 2, 2); + + vector circles; + declare.in(img); + + TEST_CYCLE() + { + HoughCircles(img, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + } + + SANITY_CHECK_NOTHING(); +} + +PERF_TEST(PerfHoughCircles2, ManySmallCircles) +{ + string filename = getDataPath("cv/imgproc/beads.jpg"); + const double dp = 1.0; + double minDist = 10; + double edgeThreshold = 90; + double accumThreshold = 11; + int minRadius = 7; + int maxRadius = 18; + + Mat img = imread(filename, IMREAD_GRAYSCALE); + ASSERT_FALSE(img.empty()) << "Unable to load source image " << filename; + + vector circles; + declare.in(img); + + TEST_CYCLE() + { + HoughCircles(img, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + } + + SANITY_CHECK_NOTHING(); +} diff --git a/modules/imgproc/src/hough.cpp b/modules/imgproc/src/hough.cpp index f8d168d474..e4c8b3f198 100644 --- a/modules/imgproc/src/hough.cpp +++ b/modules/imgproc/src/hough.cpp @@ -44,6 +44,8 @@ #include "precomp.hpp" #include "opencl_kernels_imgproc.hpp" #include "opencv2/core/hal/intrin.hpp" +#include +#include namespace cv { @@ -885,32 +887,108 @@ void HoughLinesP(InputArray _image, OutputArray _lines, * Circle Detection * \****************************************************************************************/ -struct markedCircle +struct EstimatedCircle { - markedCircle(Vec3f _c, int _idx, int _idxC) : - c(_c), idx(_idx), idxC(_idxC) {} + EstimatedCircle(Vec3f _c, int _accum) : + c(_c), accum(_accum) {} Vec3f c; - int idx, idxC; + int accum; }; -inline bool cmpCircleIndex(const markedCircle &left, const markedCircle &right) +static bool cmpAccum(const EstimatedCircle& left, const EstimatedCircle& right) { - return left.idx > right.idx; + // Compare everything so the order is completely deterministic + // Larger accum first + if (left.accum > right.accum) + return true; + else if (left.accum < right.accum) + return false; + // Larger radius first + else if (left.c[2] > right.c[2]) + return true; + else if (left.c[2] < right.c[2]) + return false; + // Smaller X + else if (left.c[0] < right.c[0]) + return true; + else if (left.c[0] > right.c[0]) + return false; + // Smaller Y + else if (left.c[1] < right.c[1]) + return true; + else if (left.c[1] > right.c[1]) + return false; + // Identical - neither object is less than the other + else + return false; +} + +inline Vec3f GetCircle(const EstimatedCircle& est) +{ + return est.c; } +class NZPointList : public std::vector +{ +private: + NZPointList(const NZPointList& other); // non-copyable + +public: + NZPointList(int reserveSize = 256) + { + reserve(reserveSize); + } +}; + +class NZPointSet +{ +private: + NZPointSet(const NZPointSet& other); // non-copyable + +public: + Mat_ positions; + + NZPointSet(int rows, int cols) : + positions(rows, cols, (uchar)0) + { + } + + void insert(const Point& pt) + { + positions(pt) = 1; + } + + void insert(const NZPointSet& from) + { + cv::bitwise_or(from.positions, positions, positions); + } + + void toList(NZPointList& list) const + { + for (int y = 0; y < positions.rows; y++) + { + const uchar *ptr = positions.ptr(y, 0); + for (int x = 0; x < positions.cols; x++) + { + if (ptr[x]) + { + list.push_back(Point(x, y)); + } + } + } + } +}; + class HoughCirclesAccumInvoker : public ParallelLoopBody { public: HoughCirclesAccumInvoker(const Mat &_edges, const Mat &_dx, const Mat &_dy, int _minRadius, int _maxRadius, float _idp, - std::vector& _accumVec, std::vector& _nz, Mutex& _mtx) : + std::vector& _accumVec, NZPointSet& _nz, Mutex& _mtx) : edges(_edges), dx(_dx), dy(_dy), minRadius(_minRadius), maxRadius(_maxRadius), idp(_idp), accumVec(_accumVec), nz(_nz), mutex(_mtx) { acols = cvCeil(edges.cols * idp), arows = cvCeil(edges.rows * idp); astep = acols + 2; -#if CV_SIMD128 - haveSIMD = hasSIMD128(); -#endif } ~HoughCirclesAccumInvoker() { } @@ -919,8 +997,7 @@ public: { Mat accumLocal = Mat(arows + 2, acols + 2, CV_32SC1, Scalar::all(0)); int *adataLocal = accumLocal.ptr(); - std::vector nzLocal; - nzLocal.reserve(256); + NZPointSet nzLocal(nz.positions.rows, nz.positions.cols); int startRow = boundaries.start; int endRow = boundaries.end; int numCols = edges.cols; @@ -942,7 +1019,6 @@ public: for(; x < numCols; ++x ) { #if CV_SIMD128 - if(haveSIMD) { v_uint8x16 v_zero = v_setzero_u8(); @@ -996,7 +1072,7 @@ _next_step: continue; Point pt = Point(x % edges.cols, y + x / edges.cols); - nzLocal.push_back(pt); + nzLocal.insert(pt); sx = cvRound((vx * idp) * 1024 / mag); sy = cvRound((vy * idp) * 1024 / mag); @@ -1025,9 +1101,11 @@ _next_step: } } - AutoLock lock(mutex); - accumVec.push_back(accumLocal); - nz.insert(nz.end(), nzLocal.begin(), nzLocal.end()); + { // TODO Try using TLSContainers + AutoLock lock(mutex); + accumVec.push_back(accumLocal); + nz.insert(nzLocal); + } } private: @@ -1035,12 +1113,9 @@ private: int minRadius, maxRadius; float idp; std::vector& accumVec; - std::vector& nz; + NZPointSet& nz; int acols, arows, astep; -#if CV_SIMD128 - bool haveSIMD; -#endif Mutex& mutex; }; @@ -1105,294 +1180,301 @@ private: Mutex& _lock; }; +static bool CheckDistance(const std::vector &circles, size_t endIdx, const Vec3f& circle, float minDist2) +{ + bool goodPoint = true; + for (uint j = 0; j < endIdx; ++j) + { + Vec3f pt = circles[j]; + float distX = circle[0] - pt[0], distY = circle[1] - pt[1]; + if (distX * distX + distY * distY < minDist2) + { + goodPoint = false; + break; + } + } + return goodPoint; +} + +static void GetCircleCenters(const std::vector ¢ers, std::vector &circles, int acols, float minDist, float dr) +{ + size_t centerCnt = centers.size(); + float minDist2 = minDist * minDist; + for (size_t i = 0; i < centerCnt; ++i) + { + int center = centers[i]; + int y = center / acols; + int x = center - y * acols; + Vec3f circle = Vec3f((x + 0.5f) * dr, (y + 0.5f) * dr, 0); + + bool goodPoint = CheckDistance(circles, circles.size(), circle, minDist2); + if (goodPoint) + circles.push_back(circle); + } +} + +static void RemoveOverlaps(std::vector& circles, float minDist) +{ + float minDist2 = minDist * minDist; + size_t endIdx = 1; + for (size_t i = 1; i < circles.size(); ++i) + { + Vec3f circle = circles[i]; + if (CheckDistance(circles, endIdx, circle, minDist2)) + { + circles[endIdx] = circle; + ++endIdx; + } + } + circles.resize(endIdx); +} + +template class HoughCircleEstimateRadiusInvoker : public ParallelLoopBody { public: - HoughCircleEstimateRadiusInvoker(const std::vector &_nz, const std::vector &_centers, std::vector &_circles, - int _acols, int _circlesMax, int _accThreshold, int _minRadius, int _maxRadius, - float _minDist, float _dp, Mutex& _mutex) : - nz(_nz), centers(_centers), circles(_circles), acols(_acols), circlesMax(_circlesMax), accThreshold(_accThreshold), - minRadius(_minRadius), maxRadius(_maxRadius), minDist(_minDist), dr(_dp), _lock(_mutex) + HoughCircleEstimateRadiusInvoker(const NZPoints &_nz, int _nzSz, const std::vector &_centers, std::vector &_circlesEst, + int _acols, int _accThreshold, int _minRadius, int _maxRadius, + float _dp, Mutex& _mutex) : + nz(_nz), nzSz(_nzSz), centers(_centers), circlesEst(_circlesEst), acols(_acols), accThreshold(_accThreshold), + minRadius(_minRadius), maxRadius(_maxRadius), dr(_dp), _lock(_mutex) { minRadius2 = (float)minRadius * minRadius; maxRadius2 = (float)maxRadius * maxRadius; - minDist = std::max(dr, minDist); - minDist *= minDist; - nzSz = (int)nz.size(); centerSz = (int)centers.size(); - - iMax = -1; - isMaxCircles = false; - isLastCenter = false; - mc.reserve(64); - loopIdx = std::vector(centerSz + 1, false); -#if CV_SIMD128 - haveSIMD = hasSIMD128(); - if(haveSIMD) - { - v_minRadius2 = v_setall_f32(minRadius2); - v_maxRadius2 = v_setall_f32(maxRadius2); - } -#endif + CV_Assert(nzSz > 0); } ~HoughCircleEstimateRadiusInvoker() {_lock.unlock();} +protected: + inline int filterCircles(const Point2f& curCenter, float* ddata) const; + void operator()(const Range &boundaries) const { - if (isMaxCircles) - return; - + std::vector circlesLocal; const int nBinsPerDr = 10; int nBins = cvRound((maxRadius - minRadius)/dr*nBinsPerDr); - std::vector bins(nBins, 0); - Mat distBuf(1, nzSz, CV_32FC1), distSqrBuf(1, nzSz, CV_32FC1); - float *ddata = distBuf.ptr(); - float *dSqrData = distSqrBuf.ptr(); + AutoBuffer bins(nBins); + AutoBuffer distBuf(nzSz), distSqrtBuf(nzSz); + float *ddata = distBuf; + float *dSqrtData = distSqrtBuf; bool singleThread = (boundaries == Range(0, centerSz)); int i = boundaries.start; - if(boundaries.end == centerSz) - isLastCenter = true; - // For each found possible center // Estimate radius and check support for(; i < boundaries.end; ++i) { - if (isMaxCircles) - return; - int ofs = centers[i]; int y = ofs / acols; int x = ofs - y * acols; //Calculate circle's center in pixels Point2f curCenter = Point2f((x + 0.5f) * dr, (y + 0.5f) * dr); - float rBest = 0; - int j = 0, nzCount = 0, maxCount = 0; - - // Check distance with previously detected valid circles - int curCircleSz = (int)circles.size(); - bool valid = checkDistance(curCenter, 0, curCircleSz); + int nzCount = filterCircles(curCenter, ddata); - if (isMaxCircles) - return; - - if(valid) + int maxCount = 0; + float rBest = 0; + if(nzCount) { -#if CV_SIMD128 - if(haveSIMD) - { - v_float32x4 v_curCenterX = v_setall_f32(curCenter.x); - v_float32x4 v_curCenterY = v_setall_f32(curCenter.y); + Mat_ distMat(1, nzCount, ddata); + Mat_ distSqrtMat(1, nzCount, dSqrtData); + sqrt(distMat, distSqrtMat); - float CV_DECL_ALIGNED(16) rbuf[4]; - int CV_DECL_ALIGNED(16) mbuf[4]; - for(; j <= nzSz - 4; j += 4) - { - v_float32x4 v_nzX, v_nzY; - v_load_deinterleave((const float*)&nz[j], v_nzX, v_nzY); - - v_float32x4 v_x = v_cvt_f32(v_reinterpret_as_s32(v_nzX)); - v_float32x4 v_y = v_cvt_f32(v_reinterpret_as_s32(v_nzY)); - - v_float32x4 v_dx = v_x - v_curCenterX; - v_float32x4 v_dy = v_y - v_curCenterY; - - v_float32x4 v_r2 = (v_dx * v_dx) + (v_dy * v_dy); - v_float32x4 vmask = (v_minRadius2 <= v_r2) & (v_r2 <= v_maxRadius2); - - v_store_aligned(rbuf, v_r2); - v_store_aligned(mbuf, v_reinterpret_as_s32(vmask)); - for(int p = 0; p < 4; p++) - { - if(mbuf[p] < 0) - { - ddata[nzCount] = rbuf[p]; nzCount++; - } - } - } - } -#endif - - // Estimate best radius - for(; j < nzSz; ++j) + memset(bins, 0, sizeof(bins[0])*bins.size()); + for(int k = 0; k < nzCount; k++) { - Point pt = nz[j]; - float _dx = curCenter.x - pt.x, _dy = curCenter.y - pt.y; - float _r2 = _dx * _dx + _dy * _dy; - - if(minRadius2 <= _r2 && _r2 <= maxRadius2) - { - ddata[nzCount] = _r2; - ++nzCount; - } + int bin = std::max(0, std::min(nBins-1, cvRound((dSqrtData[k] - minRadius)/dr*nBinsPerDr))); + bins[bin]++; } - if (isMaxCircles) - return; - - if(nzCount) + for(int j = nBins - 1; j > 0; j--) { - Mat bufRange = distSqrBuf.colRange(Range(0, nzCount)); - sqrt(distBuf.colRange(Range(0, nzCount)), bufRange); - - std::fill(bins.begin(), bins.end(), 0); - for(int k = 0; k < nzCount; k++) - { - int bin = std::max(0, std::min(nBins-1, cvRound((dSqrData[k] - minRadius)/dr*nBinsPerDr))); - bins[bin]++; - } - - if (isMaxCircles) - return; - - for(j = nBins - 1; j > 0; j--) + if(bins[j]) { - if(bins[j]) + int upbin = j; + int curCount = 0; + for(; j > upbin - nBinsPerDr && j >= 0; j--) { - int upbin = j; - int curCount = 0; - for(; j > upbin - nBinsPerDr && j >= 0; j--) - { - curCount += bins[j]; - } - float rCur = (upbin + j)/2.f /nBinsPerDr * dr + minRadius; - if((curCount * rBest >= maxCount * rCur) || - (rBest < FLT_EPSILON && curCount >= maxCount)) - { - rBest = rCur; - maxCount = curCount; - } + curCount += bins[j]; + } + float rCur = (upbin + j)/2.f /nBinsPerDr * dr + minRadius; + if((curCount * rBest >= maxCount * rCur) || + (rBest < FLT_EPSILON && curCount >= maxCount)) + { + rBest = rCur; + maxCount = curCount; } } } } - if(singleThread) + // Check if the circle has enough support + if(maxCount > accThreshold) { - // Check if the circle has enough support - if(maxCount > accThreshold) - { - circles.push_back(Vec3f(curCenter.x, curCenter.y, rBest)); + circlesLocal.push_back(EstimatedCircle(Vec3f(curCenter.x, curCenter.y, rBest), maxCount)); + } + } - if( circles.size() >= (unsigned int)circlesMax ) - return; - } + if(!circlesLocal.empty()) + { + std::sort(circlesLocal.begin(), circlesLocal.end(), cmpAccum); + if(singleThread) + { + std::swap(circlesEst, circlesLocal); } else { - _lock.lock(); - if(isMaxCircles) - { - _lock.unlock(); - return; - } + AutoLock alock(_lock); + if (circlesEst.empty()) + std::swap(circlesEst, circlesLocal); + else + circlesEst.insert(circlesEst.end(), circlesLocal.begin(), circlesLocal.end()); + } + } + } + +private: + const NZPoints &nz; + int nzSz; + const std::vector ¢ers; + std::vector &circlesEst; + int acols, accThreshold, minRadius, maxRadius; + float dr; + int centerSz; + float minRadius2, maxRadius2; + Mutex& _lock; +}; - loopIdx[i] = true; +template<> +inline int HoughCircleEstimateRadiusInvoker::filterCircles(const Point2f& curCenter, float* ddata) const +{ + int nzCount = 0; + const Point* nz_ = &nz[0]; + int j = 0; +#if CV_SIMD128 + { + const v_float32x4 v_minRadius2 = v_setall_f32(minRadius2); + const v_float32x4 v_maxRadius2 = v_setall_f32(maxRadius2); - if( maxCount > accThreshold ) - { - while(loopIdx[iMax + 1]) - ++iMax; + v_float32x4 v_curCenterX = v_setall_f32(curCenter.x); + v_float32x4 v_curCenterY = v_setall_f32(curCenter.y); - // Temporary store circle, index and already checked index for block wise testing - mc.push_back(markedCircle(Vec3f(curCenter.x, curCenter.y, rBest), - i, curCircleSz)); + float CV_DECL_ALIGNED(16) rbuf[4]; + for(; j <= nzSz - 4; j += 4) + { + v_float32x4 v_nzX, v_nzY; + v_load_deinterleave((const float*)&nz_[j], v_nzX, v_nzY); // FIXIT use proper datatype - if(i <= iMax) - { - std::sort(mc.begin(), mc.end(), cmpCircleIndex); - for(int k = (int)mc.size() - 1; k >= 0; --k) - { - if(mc[k].idx <= iMax) - { - if(checkDistance(Point2f(mc[k].c[0], mc[k].c[1]), - mc[k].idxC, (int)circles.size())) - { - circles.push_back(mc[k].c); - if(circles.size() >= (unsigned int)circlesMax) - { - isMaxCircles = true; - _lock.unlock(); - return; - } - } - mc.pop_back(); - } - else - break; - } - } - } + v_float32x4 v_x = v_cvt_f32(v_reinterpret_as_s32(v_nzX)); + v_float32x4 v_y = v_cvt_f32(v_reinterpret_as_s32(v_nzY)); - if(isLastCenter && !mc.empty()) - { - while(loopIdx[iMax + 1]) - ++iMax; + v_float32x4 v_dx = v_x - v_curCenterX; + v_float32x4 v_dy = v_y - v_curCenterY; - if(iMax == centerSz - 1) - { - std::sort(mc.begin(), mc.end(), cmpCircleIndex); - for(int k = (int)mc.size() - 1; k >= 0; --k) - { - if(checkDistance(Point2f(mc[k].c[0], mc[k].c[1]), mc[k].idxC, (int)circles.size())) - { - circles.push_back(mc[k].c); - if(circles.size() >= (unsigned int)circlesMax) - { - isMaxCircles = true; - _lock.unlock(); - return; - } - } - } - } - } - _lock.unlock(); + v_float32x4 v_r2 = (v_dx * v_dx) + (v_dy * v_dy); + v_float32x4 vmask = (v_minRadius2 <= v_r2) & (v_r2 <= v_maxRadius2); + unsigned int mask = v_signmask(vmask); + if (mask) + { + v_store_aligned(rbuf, v_r2); + if (mask & 1) ddata[nzCount++] = rbuf[0]; + if (mask & 2) ddata[nzCount++] = rbuf[1]; + if (mask & 4) ddata[nzCount++] = rbuf[2]; + if (mask & 8) ddata[nzCount++] = rbuf[3]; } } } +#endif -private: - bool checkDistance(Point2f curCenter, int startIdx, int endIdx) const + // Estimate best radius + for(; j < nzSz; ++j) { - // Check distance with previously detected circles - for(int j = startIdx; j < endIdx; ++j ) - { - float dx = circles[j][0] - curCenter.x; - float dy = circles[j][1] - curCenter.y; + const Point pt = nz_[j]; + float _dx = curCenter.x - pt.x, _dy = curCenter.y - pt.y; + float _r2 = _dx * _dx + _dy * _dy; - if( dx * dx + dy * dy < minDist ) - return false; + if(minRadius2 <= _r2 && _r2 <= maxRadius2) + { + ddata[nzCount++] = _r2; } - return true; } + return nzCount; +} - const std::vector &nz; - const std::vector ¢ers; - std::vector &circles; - int acols, circlesMax, accThreshold, minRadius, maxRadius; - float minDist, dr; +template<> +inline int HoughCircleEstimateRadiusInvoker::filterCircles(const Point2f& curCenter, float* ddata) const +{ + int nzCount = 0; + const Mat_& positions = nz.positions; + + const int rOuter = maxRadius + 1; + const Range xOuter = Range(std::max(int(curCenter.x - rOuter), 0), std::min(int(curCenter.x + rOuter), positions.cols)); + const Range yOuter = Range(std::max(int(curCenter.y - rOuter), 0), std::min(int(curCenter.y + rOuter), positions.rows)); #if CV_SIMD128 - bool haveSIMD; - v_float32x4 v_minRadius2, v_maxRadius2; + const int numSIMDPoints = 4; + + const v_float32x4 v_minRadius2 = v_setall_f32(minRadius2); + const v_float32x4 v_maxRadius2 = v_setall_f32(maxRadius2); + const v_float32x4 v_curCenterX_0123 = v_setall_f32(curCenter.x) - v_float32x4(0.0f, 1.0f, 2.0f, 3.0f); #endif - int nzSz, centerSz; - float minRadius2, maxRadius2; - mutable std::vector loopIdx; - mutable std::vector mc; - mutable volatile int iMax; - mutable volatile bool isMaxCircles, isLastCenter; - Mutex& _lock; -}; + for (int y = yOuter.start; y < yOuter.end; y++) + { + const uchar* ptr = positions.ptr(y, 0); + float dy = curCenter.y - y; + float dy2 = dy * dy; + + int x = xOuter.start; +#if CV_SIMD128 + { + const v_float32x4 v_dy2 = v_setall_f32(dy2); + const v_uint32x4 v_zero_u32 = v_setall_u32(0); + float CV_DECL_ALIGNED(16) rbuf[4]; + for (; x <= xOuter.end - 4; x += numSIMDPoints) + { + v_uint32x4 v_mask = v_load_expand_q(ptr + x); + v_mask = v_mask != v_zero_u32; + + v_float32x4 v_x = v_cvt_f32(v_setall_s32(x)); + v_float32x4 v_dx = v_x - v_curCenterX_0123; + + v_float32x4 v_r2 = (v_dx * v_dx) + v_dy2; + v_float32x4 vmask = (v_minRadius2 <= v_r2) & (v_r2 <= v_maxRadius2) & v_reinterpret_as_f32(v_mask); + unsigned int mask = v_signmask(vmask); + if (mask) + { + v_store_aligned(rbuf, v_r2); + if (mask & 1) ddata[nzCount++] = rbuf[0]; + if (mask & 2) ddata[nzCount++] = rbuf[1]; + if (mask & 4) ddata[nzCount++] = rbuf[2]; + if (mask & 8) ddata[nzCount++] = rbuf[3]; + } + } + } +#endif + for (; x < xOuter.end; x++) + { + if (ptr[x]) + { + float _dx = curCenter.x - x; + float _r2 = _dx * _dx + dy2; + if(minRadius2 <= _r2 && _r2 <= maxRadius2) + { + ddata[nzCount++] = _r2; + } + } + } + } + return nzCount; +} static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float dp, float minDist, int minRadius, int maxRadius, int cannyThreshold, - int accThreshold, int maxCircles, int kernelSize ) + int accThreshold, int maxCircles, int kernelSize, bool centersOnly) { CV_Assert(kernelSize == -1 || kernelSize == 3 || kernelSize == 5 || kernelSize == 7); dp = max(dp, 1.f); @@ -1407,19 +1489,20 @@ static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float Mutex mtx; int numThreads = std::max(1, getNumThreads()); std::vector accumVec; - std::vector nz; + NZPointSet nz(_image.rows(), _image.cols()); parallel_for_(Range(0, edges.rows), HoughCirclesAccumInvoker(edges, dx, dy, minRadius, maxRadius, idp, accumVec, nz, mtx), numThreads); - - if(nz.empty()) + int nzSz = cv::countNonZero(nz.positions); + if(nzSz <= 0) return; - Mat accum = accumVec[0].clone(); + Mat accum = accumVec[0]; for(size_t i = 1; i < accumVec.size(); i++) { accum += accumVec[i]; } + accumVec.clear(); std::vector centers; @@ -1437,49 +1520,47 @@ static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float std::vector circles; circles.reserve(256); - - if(maxCircles == 0) + if (centersOnly) + { + // Just get the circle centers + GetCircleCenters(centers, circles, accum.cols, minDist, dp); + } + else { - minDist *= minDist; - for(int i = 0; i < centerCnt; ++i) + std::vector circlesEst; + if (nzSz < maxRadius * maxRadius) { - int _centers = centers[i]; - int y = _centers / accum.cols; - int x = _centers - y * accum.cols; - - bool goodPoint = true; - for(uint j = 0; j < circles.size(); ++j) - { - Vec3f pt = circles[j]; - float distX = x - pt[0], distY = y - pt[1]; - if (distX * distX + distY * distY < minDist) - { - goodPoint = false; break; - } - } - - if(goodPoint) - circles.push_back(Vec3f((x + 0.5f) * dp, (y + 0.5f) * dp, 0)); + // Faster to use a list + NZPointList nzList(nzSz); + nz.toList(nzList); + // One loop iteration per thread if multithreaded. + parallel_for_(Range(0, centerCnt), + HoughCircleEstimateRadiusInvoker(nzList, nzSz, centers, circlesEst, accum.cols, + accThreshold, minRadius, maxRadius, dp, mtx), + numThreads); } - - if(circles.size() > 0) + else { - _circles.create(1, (int)circles.size(), CV_32FC3); - Mat(1, (int)circles.size(), CV_32FC3, &circles[0]).copyTo(_circles.getMat()); - return; + // Faster to use a matrix + // One loop iteration per thread if multithreaded. + parallel_for_(Range(0, centerCnt), + HoughCircleEstimateRadiusInvoker(nz, nzSz, centers, circlesEst, accum.cols, + accThreshold, minRadius, maxRadius, dp, mtx), + numThreads); } - } - // One loop iteration per thread if multithreaded. - parallel_for_(Range(0, centerCnt), - HoughCircleEstimateRadiusInvoker(nz, centers, circles, accum.cols, maxCircles, - accThreshold, minRadius, maxRadius, minDist, dp, mtx), - (numThreads > 1) ? centerCnt : 1); + // Sort by accumulator value + std::sort(circlesEst.begin(), circlesEst.end(), cmpAccum); + std::transform(circlesEst.begin(), circlesEst.end(), std::back_inserter(circles), GetCircle); + RemoveOverlaps(circles, minDist); + } if(circles.size() > 0) { - _circles.create(1, (int)circles.size(), CV_32FC3); - Mat(1, (int)circles.size(), CV_32FC3, &circles[0]).copyTo(_circles.getMat()); + int numCircles = std::min(maxCircles, int(circles.size())); + _circles.create(1, numCircles, CV_32FC3); + Mat(1, numCircles, CV_32FC3, &circles[0]).copyTo(_circles.getMat()); + return; } } @@ -1504,6 +1585,8 @@ static void HoughCircles( InputArray _image, OutputArray _circles, if(maxCircles < 0) maxCircles = INT_MAX; + bool centersOnly = (maxRadius < 0); + if( maxRadius <= 0 ) maxRadius = std::max( _image.rows(), _image.cols() ); else if( maxRadius <= minRadius ) @@ -1514,7 +1597,7 @@ static void HoughCircles( InputArray _image, OutputArray _circles, case CV_HOUGH_GRADIENT: HoughCirclesGradient(_image, _circles, (float)dp, (float)minDist, minRadius, maxRadius, cannyThresh, - accThresh, maxCircles, kernelSize); + accThresh, maxCircles, kernelSize, centersOnly); break; default: CV_Error( Error::StsBadArg, "Unrecognized method id. Actually only CV_HOUGH_GRADIENT is supported." ); diff --git a/modules/imgproc/test/test_houghcircles.cpp b/modules/imgproc/test/test_houghcircles.cpp new file mode 100644 index 0000000000..8fe01bd456 --- /dev/null +++ b/modules/imgproc/test/test_houghcircles.cpp @@ -0,0 +1,259 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2009, Willow Garage Inc., all rights reserved. +// Copyright (C) 2014, Itseez, Inc, all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include "test_precomp.hpp" + +#ifndef DEBUG_IMAGES +#define DEBUG_IMAGES 0 +#endif + +using namespace cv; +using namespace std; + +static string getTestCaseName(const string& picture_name, double minDist, double edgeThreshold, double accumThreshold, int minRadius, int maxRadius) +{ + string results_name = format("circles_%s_%.0f_%.0f_%.0f_%d_%d", + picture_name.c_str(), minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + string temp(results_name); + size_t pos = temp.find_first_of("\\/."); + while (pos != string::npos) { + temp.replace(pos, 1, "_"); + pos = temp.find_first_of("\\/."); + } + return temp; +} + +#if DEBUG_IMAGES +static void highlightCircles(const string& imagePath, const vector& circles, const string& outputImagePath) +{ + Mat imgDebug = imread(imagePath, IMREAD_COLOR); + const Scalar yellow(0, 255, 255); + + for (vector::const_iterator iter = circles.begin(); iter != circles.end(); ++iter) + { + const Vec3f& circle = *iter; + float x = circle[0]; + float y = circle[1]; + float r = max(circle[2], 2.0f); + cv::circle(imgDebug, Point(int(x), int(y)), int(r), yellow); + } + imwrite(outputImagePath, imgDebug); +} +#endif + +typedef std::tr1::tuple Image_MinDist_EdgeThreshold_AccumThreshold_MinRadius_MaxRadius_t; +class HoughCirclesTestFixture : public testing::TestWithParam +{ + string picture_name; + double minDist; + double edgeThreshold; + double accumThreshold; + int minRadius; + int maxRadius; + +public: + HoughCirclesTestFixture() + { + picture_name = std::tr1::get<0>(GetParam()); + minDist = std::tr1::get<1>(GetParam()); + edgeThreshold = std::tr1::get<2>(GetParam()); + accumThreshold = std::tr1::get<3>(GetParam()); + minRadius = std::tr1::get<4>(GetParam()); + maxRadius = std::tr1::get<5>(GetParam()); + } + + HoughCirclesTestFixture(const string& picture, double minD, double edge, double accum, int minR, int maxR) : + picture_name(picture), minDist(minD), edgeThreshold(edge), accumThreshold(accum), minRadius(minR), maxRadius(maxR) + { + } + + void run_test() + { + string test_case_name = getTestCaseName(picture_name, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + string filename = cvtest::TS::ptr()->get_data_path() + picture_name; + Mat src = imread(filename, IMREAD_GRAYSCALE); + EXPECT_FALSE(src.empty()) << "Invalid test image: " << filename; + + GaussianBlur(src, src, Size(9, 9), 2, 2); + + vector circles; + const double dp = 1.0; + HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + + string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; +#if DEBUG_IMAGES + highlightCircles(filename, circles, imgProc + test_case_name + ".png"); +#endif + + string xml = imgProc + "HoughCircles.xml"; + FileStorage fs(xml, FileStorage::READ); + FileNode node = fs[test_case_name]; + if (node.empty()) + { + fs.release(); + fs.open(xml, FileStorage::APPEND); + EXPECT_TRUE(fs.isOpened()) << "Cannot open sanity data file: " << xml; + fs << test_case_name << circles; + fs.release(); + fs.open(xml, FileStorage::READ); + EXPECT_TRUE(fs.isOpened()) << "Cannot open sanity data file: " << xml; + } + + vector exp_circles; + read(fs[test_case_name], exp_circles, vector()); + fs.release(); + + EXPECT_EQ(exp_circles.size(), circles.size()); + } +}; + +TEST_P(HoughCirclesTestFixture, regression) +{ + run_test(); +} + +INSTANTIATE_TEST_CASE_P(ImgProc, HoughCirclesTestFixture, testing::Combine( + // picture_name: + testing::Values("imgproc/stuff.jpg"), + // minDist: + testing::Values(20), + // edgeThreshold: + testing::Values(20), + // accumThreshold: + testing::Values(30), + // minRadius: + testing::Values(20), + // maxRadius: + testing::Values(200) + )); + +TEST(HoughCirclesTest, DefaultMaxRadius) +{ + string picture_name = "imgproc/stuff.jpg"; + const double dp = 1.0; + double minDist = 20; + double edgeThreshold = 20; + double accumThreshold = 30; + int minRadius = 20; + int maxRadius = 0; + + string filename = cvtest::TS::ptr()->get_data_path() + picture_name; + Mat src = imread(filename, IMREAD_GRAYSCALE); + EXPECT_FALSE(src.empty()) << "Invalid test image: " << filename; + + GaussianBlur(src, src, Size(9, 9), 2, 2); + + vector circles; + HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + +#if DEBUG_IMAGES + string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; + highlightCircles(filename, circles, imgProc + "HoughCirclesTest_DefaultMaxRadius.png"); +#endif + + int maxDimension = std::max(src.rows, src.cols); + + EXPECT_GT(circles.size(), size_t(0)) << "Should find at least some circles"; + for (size_t i = 0; i < circles.size(); ++i) + { + EXPECT_GE(circles[i][2], minRadius) << "Radius should be >= minRadius"; + EXPECT_LE(circles[i][2], maxDimension) << "Radius should be <= max image dimension"; + } +} + +TEST(HoughCirclesTest, CentersOnly) +{ + string picture_name = "imgproc/stuff.jpg"; + const double dp = 1.0; + double minDist = 20; + double edgeThreshold = 20; + double accumThreshold = 30; + int minRadius = 20; + int maxRadius = -1; + + string filename = cvtest::TS::ptr()->get_data_path() + picture_name; + Mat src = imread(filename, IMREAD_GRAYSCALE); + EXPECT_FALSE(src.empty()) << "Invalid test image: " << filename; + + GaussianBlur(src, src, Size(9, 9), 2, 2); + + vector circles; + HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + +#if DEBUG_IMAGES + string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; + highlightCircles(filename, circles, imgProc + "HoughCirclesTest_CentersOnly.png"); +#endif + + EXPECT_GT(circles.size(), size_t(0)) << "Should find at least some circles"; + for (size_t i = 0; i < circles.size(); ++i) + { + EXPECT_EQ(circles[i][2], 0.0f) << "Did not ask for radius"; + } +} + +TEST(HoughCirclesTest, ManySmallCircles) +{ + string picture_name = "imgproc/beads.jpg"; + const double dp = 1.0; + double minDist = 10; + double edgeThreshold = 90; + double accumThreshold = 11; + int minRadius = 7; + int maxRadius = 18; + + string filename = cvtest::TS::ptr()->get_data_path() + picture_name; + Mat src = imread(filename, IMREAD_GRAYSCALE); + EXPECT_FALSE(src.empty()) << "Invalid test image: " << filename; + + vector circles; + HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + +#if DEBUG_IMAGES + string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; + string test_case_name = getTestCaseName(picture_name, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + highlightCircles(filename, circles, imgProc + test_case_name + ".png"); +#endif + + EXPECT_GT(circles.size(), size_t(3000)) << "Should find a lot of circles"; +}