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)
pull/10463/head
Tom Becker 7 years ago committed by Alexander Alekhin
parent d9c0231475
commit 592f8d8c1b
  1. 7
      modules/imgproc/include/opencv2/imgproc.hpp
  2. 57
      modules/imgproc/perf/perf_houghcircles.cpp
  3. 639
      modules/imgproc/src/hough.cpp
  4. 259
      modules/imgproc/test/test_houghcircles.cpp

@ -2094,8 +2094,8 @@ Example: :
@note Usually the function detects the centers of circles well. However, it may fail to find correct @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 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 you know it. Or, you may set maxRadius to a negative number to return centers only without radius
radius using an additional procedure. search, and find the correct radius using an additional procedure.
@param image 8-bit, single-channel, grayscale input image. @param image 8-bit, single-channel, grayscale input image.
@param circles Output vector of found circles. Each vector is encoded as a 3-element @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 false circles may be detected. Circles, corresponding to the larger accumulator values, will be
returned first. returned first.
@param minRadius Minimum circle radius. @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 @sa fitEllipse, minEnclosingCircle
*/ */

@ -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<Vec3f> 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<Vec3f> circles;
declare.in(img);
TEST_CYCLE()
{
HoughCircles(img, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius);
}
SANITY_CHECK_NOTHING();
}

@ -44,6 +44,8 @@
#include "precomp.hpp" #include "precomp.hpp"
#include "opencl_kernels_imgproc.hpp" #include "opencl_kernels_imgproc.hpp"
#include "opencv2/core/hal/intrin.hpp" #include "opencv2/core/hal/intrin.hpp"
#include <algorithm>
#include <iterator>
namespace cv namespace cv
{ {
@ -885,32 +887,108 @@ void HoughLinesP(InputArray _image, OutputArray _lines,
* Circle Detection * * Circle Detection *
\****************************************************************************************/ \****************************************************************************************/
struct markedCircle struct EstimatedCircle
{ {
markedCircle(Vec3f _c, int _idx, int _idxC) : EstimatedCircle(Vec3f _c, int _accum) :
c(_c), idx(_idx), idxC(_idxC) {} c(_c), accum(_accum) {}
Vec3f c; 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<Point>
{
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_<uchar> 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<uchar>(y, 0);
for (int x = 0; x < positions.cols; x++)
{
if (ptr[x])
{
list.push_back(Point(x, y));
}
}
}
}
};
class HoughCirclesAccumInvoker : public ParallelLoopBody class HoughCirclesAccumInvoker : public ParallelLoopBody
{ {
public: public:
HoughCirclesAccumInvoker(const Mat &_edges, const Mat &_dx, const Mat &_dy, int _minRadius, int _maxRadius, float _idp, HoughCirclesAccumInvoker(const Mat &_edges, const Mat &_dx, const Mat &_dy, int _minRadius, int _maxRadius, float _idp,
std::vector<Mat>& _accumVec, std::vector<Point>& _nz, Mutex& _mtx) : std::vector<Mat>& _accumVec, NZPointSet& _nz, Mutex& _mtx) :
edges(_edges), dx(_dx), dy(_dy), minRadius(_minRadius), maxRadius(_maxRadius), idp(_idp), edges(_edges), dx(_dx), dy(_dy), minRadius(_minRadius), maxRadius(_maxRadius), idp(_idp),
accumVec(_accumVec), nz(_nz), mutex(_mtx) accumVec(_accumVec), nz(_nz), mutex(_mtx)
{ {
acols = cvCeil(edges.cols * idp), arows = cvCeil(edges.rows * idp); acols = cvCeil(edges.cols * idp), arows = cvCeil(edges.rows * idp);
astep = acols + 2; astep = acols + 2;
#if CV_SIMD128
haveSIMD = hasSIMD128();
#endif
} }
~HoughCirclesAccumInvoker() { } ~HoughCirclesAccumInvoker() { }
@ -919,8 +997,7 @@ public:
{ {
Mat accumLocal = Mat(arows + 2, acols + 2, CV_32SC1, Scalar::all(0)); Mat accumLocal = Mat(arows + 2, acols + 2, CV_32SC1, Scalar::all(0));
int *adataLocal = accumLocal.ptr<int>(); int *adataLocal = accumLocal.ptr<int>();
std::vector<Point> nzLocal; NZPointSet nzLocal(nz.positions.rows, nz.positions.cols);
nzLocal.reserve(256);
int startRow = boundaries.start; int startRow = boundaries.start;
int endRow = boundaries.end; int endRow = boundaries.end;
int numCols = edges.cols; int numCols = edges.cols;
@ -942,7 +1019,6 @@ public:
for(; x < numCols; ++x ) for(; x < numCols; ++x )
{ {
#if CV_SIMD128 #if CV_SIMD128
if(haveSIMD)
{ {
v_uint8x16 v_zero = v_setzero_u8(); v_uint8x16 v_zero = v_setzero_u8();
@ -996,7 +1072,7 @@ _next_step:
continue; continue;
Point pt = Point(x % edges.cols, y + x / edges.cols); Point pt = Point(x % edges.cols, y + x / edges.cols);
nzLocal.push_back(pt); nzLocal.insert(pt);
sx = cvRound((vx * idp) * 1024 / mag); sx = cvRound((vx * idp) * 1024 / mag);
sy = cvRound((vy * idp) * 1024 / mag); sy = cvRound((vy * idp) * 1024 / mag);
@ -1025,9 +1101,11 @@ _next_step:
} }
} }
AutoLock lock(mutex); { // TODO Try using TLSContainers
accumVec.push_back(accumLocal); AutoLock lock(mutex);
nz.insert(nz.end(), nzLocal.begin(), nzLocal.end()); accumVec.push_back(accumLocal);
nz.insert(nzLocal);
}
} }
private: private:
@ -1035,12 +1113,9 @@ private:
int minRadius, maxRadius; int minRadius, maxRadius;
float idp; float idp;
std::vector<Mat>& accumVec; std::vector<Mat>& accumVec;
std::vector<Point>& nz; NZPointSet& nz;
int acols, arows, astep; int acols, arows, astep;
#if CV_SIMD128
bool haveSIMD;
#endif
Mutex& mutex; Mutex& mutex;
}; };
@ -1105,294 +1180,301 @@ private:
Mutex& _lock; Mutex& _lock;
}; };
static bool CheckDistance(const std::vector<Vec3f> &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<int> &centers, std::vector<Vec3f> &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<Vec3f>& 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 NZPoints>
class HoughCircleEstimateRadiusInvoker : public ParallelLoopBody class HoughCircleEstimateRadiusInvoker : public ParallelLoopBody
{ {
public: public:
HoughCircleEstimateRadiusInvoker(const std::vector<Point> &_nz, const std::vector<int> &_centers, std::vector<Vec3f> &_circles, HoughCircleEstimateRadiusInvoker(const NZPoints &_nz, int _nzSz, const std::vector<int> &_centers, std::vector<EstimatedCircle> &_circlesEst,
int _acols, int _circlesMax, int _accThreshold, int _minRadius, int _maxRadius, int _acols, int _accThreshold, int _minRadius, int _maxRadius,
float _minDist, float _dp, Mutex& _mutex) : float _dp, Mutex& _mutex) :
nz(_nz), centers(_centers), circles(_circles), acols(_acols), circlesMax(_circlesMax), accThreshold(_accThreshold), nz(_nz), nzSz(_nzSz), centers(_centers), circlesEst(_circlesEst), acols(_acols), accThreshold(_accThreshold),
minRadius(_minRadius), maxRadius(_maxRadius), minDist(_minDist), dr(_dp), _lock(_mutex) minRadius(_minRadius), maxRadius(_maxRadius), dr(_dp), _lock(_mutex)
{ {
minRadius2 = (float)minRadius * minRadius; minRadius2 = (float)minRadius * minRadius;
maxRadius2 = (float)maxRadius * maxRadius; maxRadius2 = (float)maxRadius * maxRadius;
minDist = std::max(dr, minDist);
minDist *= minDist;
nzSz = (int)nz.size();
centerSz = (int)centers.size(); centerSz = (int)centers.size();
CV_Assert(nzSz > 0);
iMax = -1;
isMaxCircles = false;
isLastCenter = false;
mc.reserve(64);
loopIdx = std::vector<bool>(centerSz + 1, false);
#if CV_SIMD128
haveSIMD = hasSIMD128();
if(haveSIMD)
{
v_minRadius2 = v_setall_f32(minRadius2);
v_maxRadius2 = v_setall_f32(maxRadius2);
}
#endif
} }
~HoughCircleEstimateRadiusInvoker() {_lock.unlock();} ~HoughCircleEstimateRadiusInvoker() {_lock.unlock();}
protected:
inline int filterCircles(const Point2f& curCenter, float* ddata) const;
void operator()(const Range &boundaries) const void operator()(const Range &boundaries) const
{ {
if (isMaxCircles) std::vector<EstimatedCircle> circlesLocal;
return;
const int nBinsPerDr = 10; const int nBinsPerDr = 10;
int nBins = cvRound((maxRadius - minRadius)/dr*nBinsPerDr); int nBins = cvRound((maxRadius - minRadius)/dr*nBinsPerDr);
std::vector<int> bins(nBins, 0); AutoBuffer<int> bins(nBins);
Mat distBuf(1, nzSz, CV_32FC1), distSqrBuf(1, nzSz, CV_32FC1); AutoBuffer<float> distBuf(nzSz), distSqrtBuf(nzSz);
float *ddata = distBuf.ptr<float>(); float *ddata = distBuf;
float *dSqrData = distSqrBuf.ptr<float>(); float *dSqrtData = distSqrtBuf;
bool singleThread = (boundaries == Range(0, centerSz)); bool singleThread = (boundaries == Range(0, centerSz));
int i = boundaries.start; int i = boundaries.start;
if(boundaries.end == centerSz)
isLastCenter = true;
// For each found possible center // For each found possible center
// Estimate radius and check support // Estimate radius and check support
for(; i < boundaries.end; ++i) for(; i < boundaries.end; ++i)
{ {
if (isMaxCircles)
return;
int ofs = centers[i]; int ofs = centers[i];
int y = ofs / acols; int y = ofs / acols;
int x = ofs - y * acols; int x = ofs - y * acols;
//Calculate circle's center in pixels //Calculate circle's center in pixels
Point2f curCenter = Point2f((x + 0.5f) * dr, (y + 0.5f) * dr); Point2f curCenter = Point2f((x + 0.5f) * dr, (y + 0.5f) * dr);
float rBest = 0; int nzCount = filterCircles(curCenter, ddata);
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);
if (isMaxCircles) int maxCount = 0;
return; float rBest = 0;
if(nzCount)
if(valid)
{ {
#if CV_SIMD128 Mat_<float> distMat(1, nzCount, ddata);
if(haveSIMD) Mat_<float> distSqrtMat(1, nzCount, dSqrtData);
{ sqrt(distMat, distSqrtMat);
v_float32x4 v_curCenterX = v_setall_f32(curCenter.x);
v_float32x4 v_curCenterY = v_setall_f32(curCenter.y);
float CV_DECL_ALIGNED(16) rbuf[4]; memset(bins, 0, sizeof(bins[0])*bins.size());
int CV_DECL_ALIGNED(16) mbuf[4]; for(int k = 0; k < nzCount; k++)
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)
{ {
Point pt = nz[j]; int bin = std::max(0, std::min(nBins-1, cvRound((dSqrtData[k] - minRadius)/dr*nBinsPerDr)));
float _dx = curCenter.x - pt.x, _dy = curCenter.y - pt.y; bins[bin]++;
float _r2 = _dx * _dx + _dy * _dy;
if(minRadius2 <= _r2 && _r2 <= maxRadius2)
{
ddata[nzCount] = _r2;
++nzCount;
}
} }
if (isMaxCircles) for(int j = nBins - 1; j > 0; j--)
return;
if(nzCount)
{ {
Mat bufRange = distSqrBuf.colRange(Range(0, nzCount)); if(bins[j])
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]) int upbin = j;
int curCount = 0;
for(; j > upbin - nBinsPerDr && j >= 0; j--)
{ {
int upbin = j; curCount += bins[j];
int curCount = 0; }
for(; j > upbin - nBinsPerDr && j >= 0; j--) float rCur = (upbin + j)/2.f /nBinsPerDr * dr + minRadius;
{ if((curCount * rBest >= maxCount * rCur) ||
curCount += bins[j]; (rBest < FLT_EPSILON && curCount >= maxCount))
} {
float rCur = (upbin + j)/2.f /nBinsPerDr * dr + minRadius; rBest = rCur;
if((curCount * rBest >= maxCount * rCur) || maxCount = curCount;
(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 circlesLocal.push_back(EstimatedCircle(Vec3f(curCenter.x, curCenter.y, rBest), maxCount));
if(maxCount > accThreshold) }
{ }
circles.push_back(Vec3f(curCenter.x, curCenter.y, rBest));
if( circles.size() >= (unsigned int)circlesMax ) if(!circlesLocal.empty())
return; {
} std::sort(circlesLocal.begin(), circlesLocal.end(), cmpAccum);
if(singleThread)
{
std::swap(circlesEst, circlesLocal);
} }
else else
{ {
_lock.lock(); AutoLock alock(_lock);
if(isMaxCircles) if (circlesEst.empty())
{ std::swap(circlesEst, circlesLocal);
_lock.unlock(); else
return; circlesEst.insert(circlesEst.end(), circlesLocal.begin(), circlesLocal.end());
} }
}
}
private:
const NZPoints &nz;
int nzSz;
const std::vector<int> &centers;
std::vector<EstimatedCircle> &circlesEst;
int acols, accThreshold, minRadius, maxRadius;
float dr;
int centerSz;
float minRadius2, maxRadius2;
Mutex& _lock;
};
loopIdx[i] = true; template<>
inline int HoughCircleEstimateRadiusInvoker<NZPointList>::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 ) v_float32x4 v_curCenterX = v_setall_f32(curCenter.x);
{ v_float32x4 v_curCenterY = v_setall_f32(curCenter.y);
while(loopIdx[iMax + 1])
++iMax;
// Temporary store circle, index and already checked index for block wise testing float CV_DECL_ALIGNED(16) rbuf[4];
mc.push_back(markedCircle(Vec3f(curCenter.x, curCenter.y, rBest), for(; j <= nzSz - 4; j += 4)
i, curCircleSz)); {
v_float32x4 v_nzX, v_nzY;
v_load_deinterleave((const float*)&nz_[j], v_nzX, v_nzY); // FIXIT use proper datatype
if(i <= iMax) 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));
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;
}
}
}
if(isLastCenter && !mc.empty()) v_float32x4 v_dx = v_x - v_curCenterX;
{ v_float32x4 v_dy = v_y - v_curCenterY;
while(loopIdx[iMax + 1])
++iMax;
if(iMax == centerSz - 1) v_float32x4 v_r2 = (v_dx * v_dx) + (v_dy * v_dy);
{ v_float32x4 vmask = (v_minRadius2 <= v_r2) & (v_r2 <= v_maxRadius2);
std::sort(mc.begin(), mc.end(), cmpCircleIndex); unsigned int mask = v_signmask(vmask);
for(int k = (int)mc.size() - 1; k >= 0; --k) if (mask)
{ {
if(checkDistance(Point2f(mc[k].c[0], mc[k].c[1]), mc[k].idxC, (int)circles.size())) v_store_aligned(rbuf, v_r2);
{ if (mask & 1) ddata[nzCount++] = rbuf[0];
circles.push_back(mc[k].c); if (mask & 2) ddata[nzCount++] = rbuf[1];
if(circles.size() >= (unsigned int)circlesMax) if (mask & 4) ddata[nzCount++] = rbuf[2];
{ if (mask & 8) ddata[nzCount++] = rbuf[3];
isMaxCircles = true;
_lock.unlock();
return;
}
}
}
}
}
_lock.unlock();
} }
} }
} }
#endif
private: // Estimate best radius
bool checkDistance(Point2f curCenter, int startIdx, int endIdx) const for(; j < nzSz; ++j)
{ {
// Check distance with previously detected circles const Point pt = nz_[j];
for(int j = startIdx; j < endIdx; ++j ) float _dx = curCenter.x - pt.x, _dy = curCenter.y - pt.y;
{ float _r2 = _dx * _dx + _dy * _dy;
float dx = circles[j][0] - curCenter.x;
float dy = circles[j][1] - curCenter.y;
if( dx * dx + dy * dy < minDist ) if(minRadius2 <= _r2 && _r2 <= maxRadius2)
return false; {
ddata[nzCount++] = _r2;
} }
return true;
} }
return nzCount;
}
const std::vector<Point> &nz; template<>
const std::vector<int> &centers; inline int HoughCircleEstimateRadiusInvoker<NZPointSet>::filterCircles(const Point2f& curCenter, float* ddata) const
std::vector<Vec3f> &circles; {
int acols, circlesMax, accThreshold, minRadius, maxRadius; int nzCount = 0;
float minDist, dr; const Mat_<uchar>& 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 #if CV_SIMD128
bool haveSIMD; const int numSIMDPoints = 4;
v_float32x4 v_minRadius2, v_maxRadius2;
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 #endif
int nzSz, centerSz;
float minRadius2, maxRadius2;
mutable std::vector<bool> loopIdx; for (int y = yOuter.start; y < yOuter.end; y++)
mutable std::vector<markedCircle> mc; {
mutable volatile int iMax; const uchar* ptr = positions.ptr(y, 0);
mutable volatile bool isMaxCircles, isLastCenter; float dy = curCenter.y - y;
Mutex& _lock; 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, static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float dp, float minDist,
int minRadius, int maxRadius, int cannyThreshold, 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); CV_Assert(kernelSize == -1 || kernelSize == 3 || kernelSize == 5 || kernelSize == 7);
dp = max(dp, 1.f); dp = max(dp, 1.f);
@ -1407,19 +1489,20 @@ static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float
Mutex mtx; Mutex mtx;
int numThreads = std::max(1, getNumThreads()); int numThreads = std::max(1, getNumThreads());
std::vector<Mat> accumVec; std::vector<Mat> accumVec;
std::vector<Point> nz; NZPointSet nz(_image.rows(), _image.cols());
parallel_for_(Range(0, edges.rows), parallel_for_(Range(0, edges.rows),
HoughCirclesAccumInvoker(edges, dx, dy, minRadius, maxRadius, idp, accumVec, nz, mtx), HoughCirclesAccumInvoker(edges, dx, dy, minRadius, maxRadius, idp, accumVec, nz, mtx),
numThreads); numThreads);
int nzSz = cv::countNonZero(nz.positions);
if(nz.empty()) if(nzSz <= 0)
return; return;
Mat accum = accumVec[0].clone(); Mat accum = accumVec[0];
for(size_t i = 1; i < accumVec.size(); i++) for(size_t i = 1; i < accumVec.size(); i++)
{ {
accum += accumVec[i]; accum += accumVec[i];
} }
accumVec.clear();
std::vector<int> centers; std::vector<int> centers;
@ -1437,49 +1520,47 @@ static void HoughCirclesGradient(InputArray _image, OutputArray _circles, float
std::vector<Vec3f> circles; std::vector<Vec3f> circles;
circles.reserve(256); circles.reserve(256);
if (centersOnly)
if(maxCircles == 0) {
// Just get the circle centers
GetCircleCenters(centers, circles, accum.cols, minDist, dp);
}
else
{ {
minDist *= minDist; std::vector<EstimatedCircle> circlesEst;
for(int i = 0; i < centerCnt; ++i) if (nzSz < maxRadius * maxRadius)
{ {
int _centers = centers[i]; // Faster to use a list
int y = _centers / accum.cols; NZPointList nzList(nzSz);
int x = _centers - y * accum.cols; nz.toList(nzList);
// One loop iteration per thread if multithreaded.
bool goodPoint = true; parallel_for_(Range(0, centerCnt),
for(uint j = 0; j < circles.size(); ++j) HoughCircleEstimateRadiusInvoker<NZPointList>(nzList, nzSz, centers, circlesEst, accum.cols,
{ accThreshold, minRadius, maxRadius, dp, mtx),
Vec3f pt = circles[j]; numThreads);
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));
} }
else
if(circles.size() > 0)
{ {
_circles.create(1, (int)circles.size(), CV_32FC3); // Faster to use a matrix
Mat(1, (int)circles.size(), CV_32FC3, &circles[0]).copyTo(_circles.getMat()); // One loop iteration per thread if multithreaded.
return; parallel_for_(Range(0, centerCnt),
HoughCircleEstimateRadiusInvoker<NZPointSet>(nz, nzSz, centers, circlesEst, accum.cols,
accThreshold, minRadius, maxRadius, dp, mtx),
numThreads);
} }
}
// One loop iteration per thread if multithreaded. // Sort by accumulator value
parallel_for_(Range(0, centerCnt), std::sort(circlesEst.begin(), circlesEst.end(), cmpAccum);
HoughCircleEstimateRadiusInvoker(nz, centers, circles, accum.cols, maxCircles, std::transform(circlesEst.begin(), circlesEst.end(), std::back_inserter(circles), GetCircle);
accThreshold, minRadius, maxRadius, minDist, dp, mtx), RemoveOverlaps(circles, minDist);
(numThreads > 1) ? centerCnt : 1); }
if(circles.size() > 0) if(circles.size() > 0)
{ {
_circles.create(1, (int)circles.size(), CV_32FC3); int numCircles = std::min(maxCircles, int(circles.size()));
Mat(1, (int)circles.size(), CV_32FC3, &circles[0]).copyTo(_circles.getMat()); _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) if(maxCircles < 0)
maxCircles = INT_MAX; maxCircles = INT_MAX;
bool centersOnly = (maxRadius < 0);
if( maxRadius <= 0 ) if( maxRadius <= 0 )
maxRadius = std::max( _image.rows(), _image.cols() ); maxRadius = std::max( _image.rows(), _image.cols() );
else if( maxRadius <= minRadius ) else if( maxRadius <= minRadius )
@ -1514,7 +1597,7 @@ static void HoughCircles( InputArray _image, OutputArray _circles,
case CV_HOUGH_GRADIENT: case CV_HOUGH_GRADIENT:
HoughCirclesGradient(_image, _circles, (float)dp, (float)minDist, HoughCirclesGradient(_image, _circles, (float)dp, (float)minDist,
minRadius, maxRadius, cannyThresh, minRadius, maxRadius, cannyThresh,
accThresh, maxCircles, kernelSize); accThresh, maxCircles, kernelSize, centersOnly);
break; break;
default: default:
CV_Error( Error::StsBadArg, "Unrecognized method id. Actually only CV_HOUGH_GRADIENT is supported." ); CV_Error( Error::StsBadArg, "Unrecognized method id. Actually only CV_HOUGH_GRADIENT is supported." );

@ -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<Vec3f>& circles, const string& outputImagePath)
{
Mat imgDebug = imread(imagePath, IMREAD_COLOR);
const Scalar yellow(0, 255, 255);
for (vector<Vec3f>::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<string, double, double, double, int, int> Image_MinDist_EdgeThreshold_AccumThreshold_MinRadius_MaxRadius_t;
class HoughCirclesTestFixture : public testing::TestWithParam<Image_MinDist_EdgeThreshold_AccumThreshold_MinRadius_MaxRadius_t>
{
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<Vec3f> 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<Vec3f> exp_circles;
read(fs[test_case_name], exp_circles, vector<Vec3f>());
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<Vec3f> 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<Vec3f> 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<Vec3f> 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";
}
Loading…
Cancel
Save