diff --git a/modules/imgproc/include/opencv2/imgproc.hpp b/modules/imgproc/include/opencv2/imgproc.hpp index 221ee7833e..3d4bb60f57 100644 --- a/modules/imgproc/include/opencv2/imgproc.hpp +++ b/modules/imgproc/include/opencv2/imgproc.hpp @@ -473,7 +473,8 @@ enum HoughModes { /** multi-scale variant of the classical Hough transform. The lines are encoded the same way as HOUGH_STANDARD. */ HOUGH_MULTI_SCALE = 2, - HOUGH_GRADIENT = 3 //!< basically *21HT*, described in @cite Yuen90 + HOUGH_GRADIENT = 3, //!< basically *21HT*, described in @cite Yuen90 + HOUGH_GRADIENT_ALT = 4, //!< variation of HOUGH_GRADIENT to get better accuracy }; //! Variants of Line Segment %Detector @@ -2096,28 +2097,37 @@ 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 a negative number to return centers only without radius -search, and find the correct radius using an additional procedure. +you know it. Or, in the case of #HOUGH_GRADIENT method you may set maxRadius to a negative number +to return centers only without radius search, and find the correct radius using an additional procedure. + +It also helps to smooth image a bit unless it's already soft. For example, +GaussianBlur() with 7x7 kernel and 1.5x1.5 sigma or similar blurring may help. @param image 8-bit, single-channel, grayscale input image. @param circles Output vector of found circles. Each vector is encoded as 3 or 4 element floating-point vector \f$(x, y, radius)\f$ or \f$(x, y, radius, votes)\f$ . -@param method Detection method, see #HoughModes. Currently, the only implemented method is #HOUGH_GRADIENT +@param method Detection method, see #HoughModes. The available methods are #HOUGH_GRADIENT and #HOUGH_GRADIENT_ALT. @param dp Inverse ratio of the accumulator resolution to the image resolution. For example, if dp=1 , the accumulator has the same resolution as the input image. If dp=2 , the accumulator has -half as big width and height. +half as big width and height. For #HOUGH_GRADIENT_ALT the recommended value is dp=1.5, +unless some small very circles need to be detected. @param minDist Minimum distance between the centers of the detected circles. If the parameter is too small, multiple neighbor circles may be falsely detected in addition to a true one. If it is too large, some circles may be missed. -@param param1 First method-specific parameter. In case of #HOUGH_GRADIENT , it is the higher -threshold of the two passed to the Canny edge detector (the lower one is twice smaller). -@param param2 Second method-specific parameter. In case of #HOUGH_GRADIENT , it is the +@param param1 First method-specific parameter. In case of #HOUGH_GRADIENT and #HOUGH_GRADIENT_ALT, +it is the higher threshold of the two passed to the Canny edge detector (the lower one is twice smaller). +Note that #HOUGH_GRADIENT_ALT uses #Scharr algorithm to compute image derivatives, so the threshold value +shough normally be higher, such as 300 or normally exposed and contrasty images. +@param param2 Second method-specific parameter. In case of #HOUGH_GRADIENT, it is the accumulator threshold for the circle centers at the detection stage. The smaller it is, the more false circles may be detected. Circles, corresponding to the larger accumulator values, will be -returned first. +returned first. In the case of #HOUGH_GRADIENT_ALT algorithm, this is the circle "perfectness" measure. +The closer it to 1, the better shaped circles algorithm selects. In most cases 0.9 should be fine. +If you want get better detection of small circles, you may decrease it to 0.85, 0.8 or even less. +But then also try to limit the search range [minRadius, maxRadius] to avoid many false circles. @param minRadius Minimum circle radius. -@param maxRadius Maximum circle radius. If <= 0, uses the maximum image dimension. If < 0, returns -centers without finding the radius. +@param maxRadius Maximum circle radius. If <= 0, uses the maximum image dimension. If < 0, #HOUGH_GRADIENT returns +centers without finding the radius. #HOUGH_GRADIENT_ALT always computes circle radiuses. @sa fitEllipse, minEnclosingCircle */ diff --git a/modules/imgproc/src/hough.cpp b/modules/imgproc/src/hough.cpp index bbb6dfcbdb..79f34521e8 100644 --- a/modules/imgproc/src/hough.cpp +++ b/modules/imgproc/src/hough.cpp @@ -1005,6 +1005,7 @@ void HoughLinesPointSet( InputArray _point, OutputArray _lines, int lines_max, i struct EstimatedCircle { + EstimatedCircle() { accum = 0; } EstimatedCircle(Vec3f _c, int _accum) : c(_c), accum(_accum) {} Vec3f c; @@ -1710,6 +1711,530 @@ static void HoughCirclesGradient(InputArray _image, OutputArray _circles, } } +static int circle_popcnt(uint64 val) +{ + #ifdef CV_POPCNT_U64 + return CV_POPCNT_U64(val); + #else + val -= (val >> 1) & 0x5555555555555555ULL; + val = (val & 0x3333333333333333ULL) + ((val >> 2) & 0x3333333333333333ULL); + val = (val + (val >> 4)) & 0x0f0f0f0f0f0f0f0fULL; + return (int)((val * 0x0101010101010101ULL) >> 56); + #endif +} + +// The structure describes the circle "candidate" that is composed by one or more circle-like arcs. +// * rw - the sum of radiuses multiplied by the corresponding arc lengths. +// * weight - the arc length (the number of pixels it contains) +// * mask - bit mask of 64 elements that shows the coverage of the whole 0..360 degrees angular range of the circle. +// The mask of all 1's means that the whole circle is completely covered. 0's show the uncovered segments. +struct CircleData +{ + CircleData() { rw = 0; weight = 0; mask = 0; } + double rw; + int weight; + uint64 mask; +}; + +enum +{ + HOUGH_CIRCLES_ALT_BLOCK_SIZE = 10, + HOUGH_CIRCLES_ALT_MAX_CLUSTERS = 10 +}; + +static void HoughCirclesAlt( const Mat& img, std::vector& circles, double dp, double rdMinDist, + double minRadius, double maxRadius, double cannyThreshold, double minCos2 ) +{ + const int MIN_COUNT = 10; + const int RAY_FP_BITS = 10; + const int RAY_FP_SCALE = 1 << RAY_FP_BITS; + const int ACCUM_FP_BITS = 6; + const int RAY_SHIFT2 = ACCUM_FP_BITS/2; + const int ACCUM_ALPHA_ONE = 1 << RAY_SHIFT2; + const int ACCUM_ALPHA_MASK = ACCUM_ALPHA_ONE - 1; + const int RAY_SHIFT1 = RAY_FP_BITS - RAY_SHIFT2; + const int RAY_DELTA1 = 1 << (RAY_SHIFT1 - 1); + + const double ARC_DELTA = 80; + const double ARC_EPS = 0.03; + const double CIRCLE_AREA_OFFSET = 4000; + const double ARC2CLUSTER_EPS = 0.06; + const double CLUSTER_MERGE_EPS = 0.075; + const double FINAL_MERGE_DIST_EPS = 0.01; + const double FINAL_MERGE_AREA_EPS = CLUSTER_MERGE_EPS; + + if( maxRadius <= 0 ) + maxRadius = std::min(img.cols, img.rows)*0.5; + if( minRadius > maxRadius ) + std::swap(minRadius, maxRadius); + maxRadius = std::min(maxRadius, std::min(img.cols, img.rows)*0.5); + maxRadius = std::max(maxRadius, 1.); + minRadius = std::max(minRadius, 1.); + minRadius = std::min(minRadius, maxRadius); + cannyThreshold = std::max(cannyThreshold, 1.); + dp = std::max(dp, 1.); + + Mat Dx, Dy, edges; + Scharr(img, Dx, CV_16S, 1, 0); + Scharr(img, Dy, CV_16S, 0, 1); + Canny(Dx, Dy, edges, cannyThreshold/2, cannyThreshold, true); + Mat mask(img.rows + 2, img.cols + 2, CV_8U, Scalar::all(0)); + double idp = 1./dp; + int minR = cvFloor(minRadius*idp); + int maxR = cvCeil(maxRadius*idp); + int acols = cvRound(img.cols*idp); + int arows = cvRound(img.rows*idp); + Mat accum(arows + 1, acols + 1, CV_32S, Scalar::all(0)); + int* adata = accum.ptr(); + int astep = (int)accum.step1(); + minR = std::max(minR, 1); + maxR = std::max(maxR, 1); + + const uchar* edgeData = edges.ptr(); + int estep = (int)edges.step1(); + const short* dxData = Dx.ptr(); + const short* dyData = Dy.ptr(); + int dxystep = (int)Dx.step1(); + uchar* mdata = mask.ptr(); + int mstep = (int)mask.step1(); + + circles.clear(); + std::vector nz; + + std::vector stack; + const int n33[][2] = {{-1, -1}, {-1, 0}, {-1, 1}, {0, 1}, {1, 1}, {1, 0}, {1, -1}, {0, -1}}; + + for( int x = 0; x < mask.cols; x++ ) mdata[x] = mdata[(mask.rows-1)*mstep + x] = (uchar)1; + for( int y = 0; y < mask.rows; y++ ) mdata[y*mstep] = mdata[y*mstep + mask.cols-1] = (uchar)1; + mdata += mstep + 1; + + for( int y = 0; y < edges.rows; y++ ) + { + for( int x = 0; x < edges.cols; x++ ) + { + if(!edgeData[y*estep + x] || mdata[y*mstep + x]) + continue; + + mdata[y*mstep + x] = (uchar)1; + stack.push_back(Point(x, y)); + bool backtrace_mode = false; + + do + { + Point p = stack.back(); + stack.pop_back(); + int vx = dxData[p.y*dxystep + p.x]; + int vy = dyData[p.y*dxystep + p.x]; + + float mag = std::sqrt((float)vx*vx+(float)vy*vy); + nz.push_back(Vec4f((float)p.x, (float)p.y, (float)vx, (float)vy)); + CV_Assert(mdata[p.y*mstep + p.x] == 1); + + int sx = cvRound(vx * RAY_FP_SCALE / mag); + int sy = cvRound(vy * RAY_FP_SCALE / mag); + + int x0 = cvRound((p.x * idp) * RAY_FP_SCALE); + int y0 = cvRound((p.y * idp) * RAY_FP_SCALE); + + // Step from min_radius to max_radius in both directions of the gradient + for(int k1 = 0; k1 < 2; k1++ ) + { + int x1 = x0 + minR * sx; + int y1 = y0 + minR * sy; + + for(int r = minR; r <= maxR; x1 += sx, y1 += sy, r++ ) + { + int x2a = (x1 + RAY_DELTA1) >> RAY_SHIFT1, y2a = (y1 + RAY_DELTA1) >> RAY_SHIFT1; + int x2 = x2a >> RAY_SHIFT2, y2 = y2a >> RAY_SHIFT2; + if( (unsigned)x2 >= (unsigned)acols || + (unsigned)y2 >= (unsigned)arows ) + break; + + // instead of giving everything to the computed pixel of the accumulator, + // do a weighted update of 4 neighbor (2x2) pixels using bilinear interpolation. + // we do it to reduce the aliasing effect, even though it's slower + int* ptr = adata + y2*astep + x2; + int a = (x2a & ACCUM_ALPHA_MASK), b = (y2a & ACCUM_ALPHA_MASK); + ptr[0] += (ACCUM_ALPHA_ONE - a)*(ACCUM_ALPHA_ONE - b); + ptr[1] += a*(ACCUM_ALPHA_ONE - b); + ptr[astep] += (ACCUM_ALPHA_ONE - a)*b; + ptr[astep+1] += a*b; + } + + sx = -sx; sy = -sy; + } + + int neighbors = 0; + for( int k = 0; k < 8; k++ ) + { + int dy = n33[k][0], dx = n33[k][1]; + int y_ = p.y + dy, x_ = p.x + dx; + if( mdata[y_*mstep + x_] || !edgeData[y_*estep + x_]) + continue; + mdata[y_*mstep + x_] = (uchar)1; + stack.push_back(Point(x_, y_)); + neighbors++; + } + + if( neighbors == 0 ) + { + if( backtrace_mode ) + nz.pop_back(); + backtrace_mode = true; + } + else + backtrace_mode = false; + } while(!stack.empty()); + // insert a special "stop marker" in the end of each + // connected component to make sure we + // finalize and analyze the arc segment + nz.push_back(Vec4f(0.f, 0.f, 0.f, 0.f)); + } + } + + if( nz.empty() ) + return; + + // use dilation with massive ((rdMinDisp/dp)*2+1) x ((rdMinDisp/dp)*2+1) kernel. + // this trick helps us quickly find the local maxima of accumulator value + // that are at least within the specified distance from each other. + Mat accum_f, accum_max; + accum.convertTo(accum_f, CV_32F); + int niters = std::max(cvCeil(rdMinDist*idp), 1); + dilate(accum_f, accum_max, Mat(), Point(-1, -1), niters, BORDER_CONSTANT, Scalar::all(0)); + std::vector centers; + + // find the possible circle centers + for( int y = 0; y < arows; y++ ) + { + const float* adataf = accum_f.ptr(y); + const float* amaxdata = accum_max.ptr(y); + int left = -1; + for( int x = 0; x < acols; x++ ) + { + if(adataf[x] == amaxdata[x] && adataf[x] > adataf[x+astep]) + { + if(left < 0) left=x; + } + else if(left >= 0) + { + float cx = (float)((left + x - 1)*dp*0.5f); + float cy = (float)(y*dp); + centers.push_back(Point2f(cx, cy)); + left = -1; + } + } + } + + if(centers.empty()) + return; + + float minR2 = (float)(minRadius*minRadius); + float maxR2 = (float)(maxRadius*maxRadius); + int nstripes = (int)((centers.size() + HOUGH_CIRCLES_ALT_BLOCK_SIZE-1)/HOUGH_CIRCLES_ALT_BLOCK_SIZE); + const int nnz = (int)nz.size(); + Mutex cmutex; + + // Check each possible pair (edge_pixel[i], circle_center[j]). + // For each circle form the clusters to identify possible radius values. + // Several clusters (up to 10) are maintained to help to filter out false alarms and + // to support the concentric circle cases. + + // inside parallel for we process the next "HOUGH_CIRCLES_ALT_BLOCK_SIZE" circles + parallel_for_(Range(0, nstripes), [&](const Range& r) + { + CircleData cdata[HOUGH_CIRCLES_ALT_BLOCK_SIZE*HOUGH_CIRCLES_ALT_MAX_CLUSTERS]; + CircleData arc[HOUGH_CIRCLES_ALT_BLOCK_SIZE]; + int prev_idx[HOUGH_CIRCLES_ALT_BLOCK_SIZE]; + + std::vector local_circles; + for(int j0 = r.start*HOUGH_CIRCLES_ALT_BLOCK_SIZE; j0 < r.end*HOUGH_CIRCLES_ALT_BLOCK_SIZE; j0 += HOUGH_CIRCLES_ALT_BLOCK_SIZE) + { + const Vec4f* nzdata = &nz[0]; + const Point2f* cc = ¢ers[j0]; + int nc = std::min((int)(centers.size() - j0), (int)HOUGH_CIRCLES_ALT_BLOCK_SIZE); + if(nc <= 0) break; + + // reset the statistics about the clusters + for( int j = 0; j < nc; j++ ) + { + for( int k = 0; k < HOUGH_CIRCLES_ALT_BLOCK_SIZE; k++ ) + cdata[j*HOUGH_CIRCLES_ALT_MAX_CLUSTERS + k] = CircleData(); + arc[j] = CircleData(); + arc[j].weight = 1; // avoid division by zero + prev_idx[j] = -2; // we compare the current index "i" with prev_idx[j]+1 + // to check whether we are still at the current Canny + // connected component. so we initially set it to -2 + // to make sure that the initial check gives "false". + } + + for( int i = 0; i < nnz; i++ ) + { + Vec4f v = nzdata[i]; + float x = v[0], y = v[1], vx = v[2], vy = v[3], mag2 = vx*vx + vy*vy; + bool stop_marker = x == 0.f && y == 0.f && vx == 0.f && vy == 0.f; + + for( int j = 0; j < nc; j++ ) + { + float cx = cc[j].x, cy = cc[j].y; + float dx = x - cx, dy = y - cy; + float rij2 = dx*dx + dy*dy; + // check that i-th pixel is within the specified distance range from the center + if( (rij2 > maxR2 || rij2 < minR2) && i < nnz-1 ) continue; + float dv = dx*vx + dy*vy; + // check that the line segment connecting the edge pixel and the center and + // the gradient at the edge pixel are almost collinear + if( (double)dv*dv < (double)minCos2*mag2*rij2 && i < nnz-1 ) continue; + float rij = std::sqrt(rij2); + + CircleData& arc_j = arc[j]; + double r_arc = arc_j.rw/arc_j.weight; + int di0 = 0; + int prev = prev_idx[j]; + prev_idx[j] = i; + + // update the arc statistics if it still looks like an arc + if( std::abs(rij - r_arc) < (r_arc + ARC_DELTA)*ARC_EPS && prev+1 == i && !stop_marker ) + { + arc_j.rw += rij; + arc_j.weight++; + di0 = 1; + r_arc = arc_j.rw/arc_j.weight; + if( i < nnz -1 ) + continue; + } + + // otherwise (or in the very end) store the arc in the cluster collection, + // if the arc is long enough. + if( arc_j.weight >= MIN_COUNT && arc_j.weight >= r_arc*0.15 ) + { + // before doing it, compute the angular range coverage (the mask). + uint64 mval = 0; + for( int di = 0; di < arc_j.weight; di++ ) + { + int i1 = prev + di0 - di; + Vec4f u = nz[i1]; + float x1 = u[0], y1 = u[1]; + float dx1 = x1 - cx, dy1 = y1 - cy; + float af = fastAtan2(dy1, dx1)*(64.f/360.f); + int a = (cvFloor(af) & 63); + int b = (a + 1) & 63; + af -= a; + // this is another protection from aliasing effects + if( af <= 0.25f ) + mval |= (uint64)1 << a; + else if( af > 0.75f ) + mval |= (uint64)1 << b; + else + mval |= ((uint64)1 << a) | ((uint64)1 << b); + } + + double min_eps = DBL_MAX; + int min_mval = (int)(sizeof(mval)*8+1); + int k = 0, best_k = -1, subst_k = -1; + CircleData* cdata_j = &cdata[j*HOUGH_CIRCLES_ALT_MAX_CLUSTERS]; + + for( ; k < HOUGH_CIRCLES_ALT_MAX_CLUSTERS; k++ ) + { + CircleData& cjk = cdata_j[k]; + if( cjk.weight == 0 ) + break; // it means that there is no more valid clusters + double rk = cjk.rw/cjk.weight; + // Compute and use the weighted "cluster with arc" area instead of + // just cluster area or just arc area or their sum. This is because the cluster can + // be small and the arc can be big, or vice versa. Weighted area is more robust. + double r2avg = (rk*rk*cjk.weight + r_arc*r_arc*arc_j.weight)/(cjk.weight + arc_j.weight); + // It seems to be more robust to compare circle areas (without "pi" scale) + // instead of radiuses. When we compare radiuses, when depending on the ALPHA, + // different big circles are merged too easily, or different small circles stay different. + if( std::abs(rk*rk - r_arc*r_arc) < (r2avg + CIRCLE_AREA_OFFSET)*ARC2CLUSTER_EPS ) + { + double eps = std::abs(rk - r_arc)/rk; + if( eps < min_eps ) + { + min_eps = eps; + best_k = k; + } + } + else + { + // Select the cluster with the worst angular coverage. + // We use the angular coverage instead of the arc weight + // in order to protect real small circles + // from "fake" bigger circles with bigger "support". + int pcnt = circle_popcnt(cjk.mask); + if( pcnt < min_mval ) + { + min_mval = pcnt; + subst_k = k; + } + } + } + + if( best_k >= 0 ) // if found the match, merge the arc into the cluster + { + CircleData& cjk = cdata_j[best_k]; + cjk.rw += arc_j.rw; + cjk.weight += arc_j.weight; + cjk.mask |= mval; + } + else + { + if( k < HOUGH_CIRCLES_ALT_MAX_CLUSTERS ) + subst_k = k; // if we have empty space, just add the new cluster, do not throw anything + CircleData& cjk0 = cdata_j[subst_k]; + + // here was the code that attempts to merge the thrown-away cluster with others, + // but apparently it does not have any noticeable effect, + // so we removed it for the sake of simplicity ... + + // add the new cluster + cjk0.rw = arc_j.rw; + cjk0.weight = arc_j.weight; + cjk0.mask = mval; + } + } + // reset the arc statistics. + arc_j.rw = stop_marker ? 0. : rij; + arc_j.weight = 1; + // do not clean arc_j.mval, because we do not alter it. + } + } + + // now merge the final clusters for each particular circle center (cx, cy) + for( int j = 0; j < nc; j++ ) + { + CircleData* cdata_j = &cdata[j*HOUGH_CIRCLES_ALT_MAX_CLUSTERS]; + float cx = cc[j].x, cy = cc[j].y; + + for( int k = 0; k < HOUGH_CIRCLES_ALT_MAX_CLUSTERS; k++ ) + { + CircleData& cjk = cdata_j[k]; + if( cjk.weight == 0 ) + continue; + + // Let in only more or less significant clusters. + // Small clusters more likely correspond to a noise + // (otherwise they would grew more substantial during the + // cluster construction phase). + // Processing those noisy clusters takes time and + // potentially decreases accuracy of computed radiuses + // of good clusters. + double rjk = cjk.rw/cjk.weight; + if( cjk.weight < rjk || circle_popcnt(cjk.mask) < 15 ) + cjk.weight = 0; + } + + // extensive O(nclusters^2) cluster merge algorithm, but since the number + // of clusters is limited with a modest constant HOUGH_CIRCLES_ALT_MAX_CLUSTERS, + // it's still O(1) algorithm :) + for( int k = 0; k < HOUGH_CIRCLES_ALT_MAX_CLUSTERS; k++ ) + { + CircleData& cjk = cdata_j[k]; + if( cjk.weight == 0 ) + continue; + double rk = cjk.rw/cjk.weight; + + int l = k+1; + for( ; l < HOUGH_CIRCLES_ALT_MAX_CLUSTERS; l++ ) + { + CircleData& cjl = cdata_j[l]; + if( l == k || cjl.weight == 0 ) + continue; + double rl = cjl.rw/cjl.weight; + // Here we use a simple sum of areas (without "pi" scale) instead of weighted + // sum just for simplicity and potentially for better accuracy. + if( std::abs(rk*rk - rl*rl) < (rk*rk + rl*rl + CIRCLE_AREA_OFFSET)*CLUSTER_MERGE_EPS) + { + cjk.rw += cjl.rw; + cjk.weight += cjl.weight; + cjk.mask |= cjl.mask; + rk = cjk.rw/cjk.weight; + cjl.weight = 0; + l = -1; // try to merge other clusters again with the updated k-th cluster + } + } + } + + for( int k = 0; k < HOUGH_CIRCLES_ALT_MAX_CLUSTERS; k++ ) + { + CircleData& cjk = cdata_j[k]; + if( cjk.weight == 0 ) + continue; + double rk = cjk.rw/cjk.weight; + uint64 mask_jk = cjk.mask, mask_jk0 = (mask_jk + 1) ^ mask_jk; + int count = 0, count0 = -1, runlen = 0, max_runlen = 0; + int prev_bit = 0; + for( int b = 0; b < 64; b++, mask_jk >>= 1, mask_jk0 >>= 1 ) + { + int bit_k = (mask_jk & 1) != 0; + count += bit_k; + count0 += (mask_jk0 & 1) != 0; + if(bit_k == prev_bit) { runlen++; continue; } + if(prev_bit == 1) + max_runlen = std::max(max_runlen, runlen); + runlen = 1; + prev_bit = bit_k; + } + if( prev_bit == 1) + max_runlen = std::max(max_runlen, runlen + (count < 64 ? count0 : 0)); + + // Those constants are the results of fine-tuning. + // Basically, by lowering thresholds more real circles, as well as fake circles, are accepted. + // By raising the thresholds you get less real circles and less false alarms. + // A better and more safe way to obtain better detection results is to regulate + // [minRadius, maxRadius] range and to play with minCos2 parameter. + // May be some classifier can be trained that takes the weight, + // circle radius and the bit mask as inputs and produces the verdict. + bool accepted = (cjk.weight >= rk*3 && count >= 35 && max_runlen >= 20) || count >= 55; + //if(debug) + //printf("[%c]. cx=%.1f, cy=%.1f, r=%.1f, weight=%d, count=%d, max_runlen=%d, mask=%016llx\n", + // (accepted ? '+' : '-'), cx, cy, rk, cjk.weight, count, max_runlen, cjk.mask); + + if( accepted ) + local_circles.push_back(EstimatedCircle(Vec3f(cx, cy, (float)rk), cjk.weight)); + } + } + } + if(!local_circles.empty()) + { + cmutex.lock(); + std::copy(local_circles.begin(), local_circles.end(), std::back_inserter(circles)); + cmutex.unlock(); + } + }); + + // The final circle merge procedure. + // This is O(ncircles^2) algorithm + // and it can take a long time in some specific scenarious. + // But most of the time it's very fast. + size_t i0 = 0, nc = circles.size(); + for( size_t i = 0; i < nc; i++ ) + { + if( circles[i].accum == 0 ) continue; + EstimatedCircle& ci = circles[i0] = circles[i]; + for( size_t j = i+1; j < nc; j++ ) + { + EstimatedCircle cj = circles[j]; + if( cj.accum == 0 ) continue; + float dx = ci.c[0] - cj.c[0], dy = ci.c[1] - cj.c[1]; + float r2 = dx*dx + dy*dy; + float rs = ci.c[2] + cj.c[2]; + if( r2 > rs*rs*FINAL_MERGE_DIST_EPS) + continue; + if( std::abs(ci.c[2]*ci.c[2] - cj.c[2]*cj.c[2]) < + (ci.c[2]*ci.c[2] + cj.c[2]*cj.c[2] + CIRCLE_AREA_OFFSET)*FINAL_MERGE_AREA_EPS ) + { + int wi = ci.accum, wj = cj.accum; + if( wi < wj ) std::swap(ci, cj); + circles[j].accum = 0; + } + } + i0++; + } + circles.resize(i0); +} + static void HoughCircles( InputArray _image, OutputArray _circles, int method, double dp, double minDist, double param1, double param2, @@ -1728,26 +2253,29 @@ static void HoughCircles( InputArray _image, OutputArray _circles, CV_Assert(!_image.empty() && _image.type() == CV_8UC1 && (_image.isMat() || _image.isUMat())); CV_Assert(_circles.isMat() || _circles.isVector()); - if( dp <= 0 || minDist <= 0 || param1 <= 0 || param2 <= 0) - CV_Error( Error::StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" ); - - int cannyThresh = cvRound(param1), accThresh = cvRound(param2), kernelSize = cvRound(param3); - - minRadius = std::max(0, minRadius); - - if(maxCircles < 0) - maxCircles = INT_MAX; - - bool centersOnly = (maxRadius < 0); - - if( maxRadius <= 0 ) - maxRadius = std::max( _image.rows(), _image.cols() ); - else if( maxRadius <= minRadius ) - maxRadius = minRadius + 2; + if( dp <= 0 || minDist <= 0 || param1 <= 0) + CV_Error( Error::StsOutOfRange, "dp, min_dist and canny_threshold must be all positive numbers" ); switch( method ) { - case CV_HOUGH_GRADIENT: + case HOUGH_GRADIENT: + { + int cannyThresh = cvRound(param1), accThresh = cvRound(param2), kernelSize = cvRound(param3); + minRadius = std::max(0, minRadius); + + if( param2 <= 0 ) + CV_Error( Error::StsOutOfRange, "acc_threshold must be a positive number" ); + + if(maxCircles < 0) + maxCircles = INT_MAX; + + bool centersOnly = (maxRadius < 0); + + if( maxRadius <= 0 ) + maxRadius = std::max( _image.rows(), _image.cols() ); + else if( maxRadius <= minRadius ) + maxRadius = minRadius + 2; + if (type == CV_32FC3) HoughCirclesGradient(_image, _circles, (float)dp, (float)minDist, minRadius, maxRadius, cannyThresh, @@ -1758,6 +2286,33 @@ static void HoughCircles( InputArray _image, OutputArray _circles, accThresh, maxCircles, kernelSize, centersOnly); else CV_Error(Error::StsError, "Internal error"); + } + break; + case HOUGH_GRADIENT_ALT: + { + std::vector circles; + Mat image = _image.getMat(); + HoughCirclesAlt(image, circles, dp, minDist, minRadius, maxRadius, param1, param2); + std::sort(circles.begin(), circles.end(), cmpAccum); + size_t i, ncircles = circles.size(); + + if( type == CV_32FC4 ) + { + std::vector cw(ncircles); + for( i = 0; i < ncircles; i++ ) + cw[i] = GetCircle4f(circles[i]); + Mat(cw).copyTo(_circles); + } + else if( type == CV_32FC3 ) + { + std::vector cwow(ncircles); + for( i = 0; i < ncircles; i++ ) + cwow[i] = GetCircle(circles[i]); + Mat(cwow).copyTo(_circles); + } + else + CV_Error(Error::StsError, "Internal error"); + } 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 index dded1b3ea8..de92800976 100644 --- a/modules/imgproc/test/test_houghcircles.cpp +++ b/modules/imgproc/test/test_houghcircles.cpp @@ -178,26 +178,34 @@ INSTANTIATE_TEST_CASE_P(ImgProc, HoughCirclesTestFixture, testing::Combine( testing::Values(200) )); -TEST(HoughCirclesTest, DefaultMaxRadius) + +class HoughCirclesTest : public testing::TestWithParam +{ +protected: + HoughModes method; +public: + HoughCirclesTest() { method = GetParam(); } +}; + +TEST_P(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); + double dp = 1.0; + double minDist = 20.0; + double edgeThreshold = 20.0; + double param2 = method == HOUGH_GRADIENT_ALT ? 0.9 : 30.; + int minRadius = method == HOUGH_GRADIENT_ALT ? 10 : 20; + int maxRadius = 0; + vector circles; vector circles4f; - HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); - HoughCircles(src, circles4f, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + HoughCircles(src, circles, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); + HoughCircles(src, circles4f, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); #if DEBUG_IMAGES string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; @@ -206,7 +214,15 @@ TEST(HoughCirclesTest, DefaultMaxRadius) int maxDimension = std::max(src.rows, src.cols); - EXPECT_GT(circles.size(), size_t(0)) << "Should find at least some circles"; + if(method == HOUGH_GRADIENT_ALT) + { + EXPECT_EQ(circles.size(), size_t(3)) << "Should find 3 circles"; + } + else + { + 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"; @@ -214,60 +230,80 @@ TEST(HoughCirclesTest, DefaultMaxRadius) } } -TEST(HoughCirclesTest, CentersOnly) +TEST_P(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); + double dp = 1.0; + double minDist = 20.0; + double edgeThreshold = 20.0; + double param2 = method == HOUGH_GRADIENT_ALT ? 0.9 : 30.; + int minRadius = method == HOUGH_GRADIENT_ALT ? 10 : 20; + int maxRadius = -1; vector circles; vector circles4f; - HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); - HoughCircles(src, circles4f, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + + HoughCircles(src, circles, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); + HoughCircles(src, circles4f, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); #if DEBUG_IMAGES string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; - highlightCircles(filename, circles, imgProc + "HoughCirclesTest_CentersOnly.png"); + highlightCircles(filename, circles, imgProc + "HoughCirclesTest_DefaultMaxRadius.png"); #endif - EXPECT_GT(circles.size(), size_t(0)) << "Should find at least some circles"; + if(method == HOUGH_GRADIENT_ALT) + { + EXPECT_EQ(circles.size(), size_t(3)) << "Should find 3 circles"; + } + else + { + 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"; + if( method == HOUGH_GRADIENT ) + { + EXPECT_EQ(circles[i][2], 0.0f) << "Did not ask for radius"; + } EXPECT_EQ(circles[i][0], circles4f[i][0]); EXPECT_EQ(circles[i][1], circles4f[i][1]); EXPECT_EQ(circles[i][2], circles4f[i][2]); } } -TEST(HoughCirclesTest, ManySmallCircles) +TEST_P(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; + const double dp = method == HOUGH_GRADIENT_ALT ? 1.5 : 1.0; + double minDist = 10; + double edgeThreshold = 90; + double accumThreshold = 11; + double minCos2 = 0.85; + double param2 = method == HOUGH_GRADIENT_ALT ? minCos2 : accumThreshold; + int minRadius = 7; + int maxRadius = 18; + int ncircles_min = method == HOUGH_GRADIENT_ALT ? 2000 : 3000; + + Mat src_smooth; + if( method == HOUGH_GRADIENT_ALT ) + GaussianBlur(src, src_smooth, Size(7, 7), 1.5, 1.5); + else + src.copyTo(src_smooth); vector circles; vector circles4f; - HoughCircles(src, circles, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); - HoughCircles(src, circles4f, CV_HOUGH_GRADIENT, dp, minDist, edgeThreshold, accumThreshold, minRadius, maxRadius); + HoughCircles(src_smooth, circles, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); + HoughCircles(src_smooth, circles4f, method, dp, minDist, edgeThreshold, param2, minRadius, maxRadius); #if DEBUG_IMAGES string imgProc = string(cvtest::TS::ptr()->get_data_path()) + "imgproc/"; @@ -275,8 +311,11 @@ TEST(HoughCirclesTest, ManySmallCircles) highlightCircles(filename, circles, imgProc + test_case_name + ".png"); #endif - EXPECT_GT(circles.size(), size_t(3000)) << "Should find a lot of circles"; + EXPECT_GT(circles.size(), size_t(ncircles_min)) << "Should find a lot of circles"; EXPECT_EQ(circles.size(), circles4f.size()); } +INSTANTIATE_TEST_CASE_P(HoughGradient, HoughCirclesTest, testing::Values(HOUGH_GRADIENT)); +INSTANTIATE_TEST_CASE_P(HoughGradientAlt, HoughCirclesTest, testing::Values(HOUGH_GRADIENT_ALT)); + }} // namespace