cvHoughCircle 解析

查志強發表於2016-04-08

【原文：http://blog.csdn.net/mysee1989/article/details/12405949】

在網上看了一篇關於cvHoughCircle函式的解析，轉載過來看看。學習Opencv，不少人只能呼叫函式而不能作出適配或者修改，這是一種遺憾，貌似很少有這方面的部落格，不知道有人推薦一些關於OpenCV原始碼解析的資料沒有？

一、基本原理

opencv中圓識別的基本原理如下：

1、canny運算元求影象的單畫素二值化邊緣

2、假設我們需要找半徑為R的所有圓，則對於邊緣圖中的每一個邊緣點，該邊緣點的切線的法線方向上(正負兩個方向)，尋找到該邊緣點距離為R的點，將該點的計數加1(初始化所有點的計數都是0)

3、找到計數值大於門限值的點，即圓心所在的點

 二、程式碼分析

程式碼在/modules\imgproc\src\hough.cpp檔案icvHoughCirclesGradient函式中

 

[cpp] view plaincopy
  1. static void  
  2. icvHoughCirclesGradient( CvMat* img, float dp, float min_dist,  
  3.                          int min_radius, int max_radius,  
  4.                          int canny_threshold, int acc_threshold,  
  5.                          CvSeq* circles, int circles_max )  
  6. {  
  7. //引數：  
  8. //img: 輸入影象  
  9. //dp: 識別精度,1.0表示按照原圖精度  
  10. //min_dist: 圓心點位置識別精度  
  11. //min_radius: 所需要找的圓的最小半徑  
  12. //max_radius：所需要找的圓的最大半徑  
  13. //canny_threshold：canny運算元的高閥值  
  14. //acc_threshold：累加器閥值，計數大於改閥值的點即被認為是可能的圓心  
  15. //circles: 儲存找到的符合條件的所有圓  
  16. //circles_max: 最多需要的找到的圓的個數  
  17.   
  18.     const int SHIFT = 10, ONE = 1 << SHIFT, R_THRESH = 30;  
  19.     cv::Ptr<CvMat> dx, dy;  
  20.     cv::Ptr<CvMat> edges, accum, dist_buf;  
  21.     std::vector<int> sort_buf;  
  22.     cv::Ptr<CvMemStorage> storage;  
  23.       
  24.     int x, y, i, j, k, center_count, nz_count;  
  25.     float min_radius2 = (float)min_radius*min_radius;  
  26.     float max_radius2 = (float)max_radius*max_radius;  
  27.     int rows, cols, arows, acols;  
  28.     int astep, *adata;  
  29.     float* ddata;  
  30.     CvSeq *nz, *centers;  
  31.     float idp, dr;  
  32.     CvSeqReader reader;  
  33.       
  34.     //canny運算元求單畫素二值化邊緣，儲存在edges變數中  
  35.     edges = cvCreateMat( img->rows, img->cols, CV_8UC1 );  
  36.     cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 );  
  37.       
  38.     //sobel運算元求水平和垂直方向的邊緣，用於計算邊緣點的法線方向  
  39.     dx = cvCreateMat( img->rows, img->cols, CV_16SC1 );  
  40.     dy = cvCreateMat( img->rows, img->cols, CV_16SC1 );  
  41.     cvSobel( img, dx, 1, 0, 3 );  
  42.     cvSobel( img, dy, 0, 1, 3 );  
  43.       
  44.     //dp表示識別精度  
  45.     if( dp < 1.f )  
  46.         dp = 1.f;  
  47.     idp = 1.f/dp;  
  48.       
  49.     //accum用作累加器，包含影象中每一個點的計數。影象中每一個點都有一個計數，點的計數表示每一個canny邊緣點法線方向上，  
  50.     //到該點距離為R的邊緣點的個數，初始化為0  
  51.     accum = cvCreateMat( cvCeil(img->rows*idp)+2, cvCeil(img->cols*idp)+2, CV_32SC1 );  
  52.     cvZero(accum);  
  53.       
  54.     storage = cvCreateMemStorage();  
  55.     nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage );  
  56.       
  57.     //centers用於儲存可能的圓心點  
  58.     centers = cvCreateSeq( CV_32SC1, sizeof(CvSeq), sizeof(int), storage );  
  59.       
  60.     rows = img->rows;  
  61.     cols = img->cols;  
  62.     arows = accum->rows - 2;  
  63.     acols = accum->cols - 2;  
  64.     adata = accum->data.i;  
  65.     astep = accum->step/sizeof(adata[0]);  
  66.       
  67.     //以下這個迴圈用於獲取所有可能的圓邊緣點，儲存在nz中，同時設定  
  68.     //累加器中的值  
  69.     for( y = 0; y < rows; y++ )  
  70.     {  
  71.         const uchar* edges_row = edges->data.ptr + y*edges->step;  
  72.         const short* dx_row = (const short*)(dx->data.ptr + y*dx->step);  
  73.         const short* dy_row = (const short*)(dy->data.ptr + y*dy->step);  
  74.           
  75.         for( x = 0; x < cols; x++ )  
  76.         {  
  77.             float vx, vy;  
  78.             int sx, sy, x0, y0, x1, y1, r, k;  
  79.             CvPoint pt;  
  80.       
  81.             vx = dx_row[x];  
  82.             vy = dy_row[x];  
  83.       
  84.             if( !edges_row[x] || (vx == 0 && vy == 0) )  
  85.                 continue;  
  86.       
  87.             float mag = sqrt(vx*vx+vy*vy);  
  88.             assert( mag >= 1 );  
  89.               
  90.             //sx表示cos, sy表示sin  
  91.             sx = cvRound((vx*idp)*ONE/mag);  
  92.             sy = cvRound((vy*idp)*ONE/mag);  
  93.               
  94.             x0 = cvRound((x*idp)*ONE);  
  95.             y0 = cvRound((y*idp)*ONE);  
  96.               
  97.             //迴圈兩次表示需要計算兩個方向，法線方向和法線的反方向  
  98.             for( k = 0; k < 2; k++ )  
  99.             {  
  100.                 //半徑方向的水平增量和垂直增量  
  101.                 x1 = x0 + min_radius * sx;  
  102.                 y1 = y0 + min_radius * sy;  
  103.                   
  104.                 //在法線方向和反方向上，距離邊緣點的距離為輸入的最大半徑和最小半徑範圍內找點  
  105.                 //每找到一個點，該點的累加器計數就加1  
  106.                 for( r = min_radius; r <= max_radius; x1 += sx, y1 += sy, r++ )  
  107.                 {  
  108.                     int x2 = x1 >> SHIFT, y2 = y1 >> SHIFT;  
  109.                     if( (unsigned)x2 >= (unsigned)acols || (unsigned)y2 >= (unsigned)arows )  
  110.                         break;  
  111.                     adata[y2*astep + x2]++;  
  112.                 }  
  113.                 //反方向  
  114.                 sx = -sx; sy = -sy;  
  115.             }  
  116.               
  117.             //儲存可能的圓邊緣點  
  118.             pt.x = x; pt.y = y;  
  119.             cvSeqPush( nz, &pt );  
  120.         }  
  121.     }  
  122.   
  123.     nz_count = nz->total;  
  124.     if( !nz_count )  
  125.         return;  
  126.       
  127.     //累加器中，計數大於閥值的點，被認為可能的圓心點。因為計算各點計數過程中，距離有誤差，所以  
  128.     //在與閥值比較時，該點計數先要與4鄰域內的各個點的計數比較，最大者才能和閥值比較。可能的圓心  
  129.     //點儲存在centers中  
  130.     for( y = 1; y < arows - 1; y++ )  
  131.     {  
  132.         for( x = 1; x < acols - 1; x++ )  
  133.         {  
  134.             int base = y*(acols+2) + x;  
  135.             if( adata[base] > acc_threshold &&  
  136.                 adata[base] > adata[base-1] && adata[base] > adata[base+1] &&  
  137.                 adata[base] > adata[base-acols-2] && adata[base] > adata[base+acols+2] )  
  138.             cvSeqPush(centers, &base);  
  139.         }  
  140.     }  
  141.   
  142.     center_count = centers->total;  
  143.     if( !center_count )  
  144.         return;  
  145.   
  146.     sort_buf.resize( MAX(center_count,nz_count) );  
  147.       
  148.     //連結串列結構的certers轉化成連續儲存結構sort_buf  
  149.     cvCvtSeqToArray( centers, &sort_buf[0] );  
  150.       
  151.     //經過icvHoughSortDescent32s函式後，以sort_buf中元素作為adata陣列下標,   
  152.     //adata中的元素降序排列, 即adata[sort_buf[0]]是adata所有元素中最大的,   
  153.     //adata[sort_buf[center_count-1]]是所有元素中最小的  
  154.     icvHoughSortDescent32s( &sort_buf[0], center_count, adata );  
  155.       
  156.     cvClearSeq( centers );  
  157.       
  158.     //經過排序後的元素，重新以連結串列形式儲存到centers中  
  159.     cvSeqPushMulti( centers, &sort_buf[0], center_count );  
  160.       
  161.     dist_buf = cvCreateMat( 1, nz_count, CV_32FC1 );  
  162.     ddata = dist_buf->data.fl;  
  163.       
  164.     dr = dp;  
  165.     min_dist = MAX( min_dist, dp );  
  166.     min_dist *= min_dist;  
  167.       
  168.     //對於每一個可能的圓心點，計算所有邊緣點到該圓心點的距離。由於centers中的  
  169.     //元素已經經過前面排序，所以累加器計數最大的可能圓心點最先進行下面的操作  
  170.     for( i = 0; i < centers->total; i++ )  
  171.     {  
  172.         int ofs = *(int*)cvGetSeqElem( centers, i );  
  173.         y = ofs/(acols+2) - 1;  
  174.         x = ofs - (y+1)*(acols+2) - 1;  
  175.         float cx = (float)(x*dp), cy = (float)(y*dp);  
  176.         float start_dist, dist_sum;  
  177.         float r_best = 0, c[3];  
  178.         int max_count = R_THRESH;  
  179.           
  180.         //如果該可能的圓心點和已經確認的圓心點的距離小於閥值，則表示  
  181.         //這個圓心點和已經確認的圓心點是同一個點  
  182.         for( j = 0; j < circles->total; j++ )  
  183.         {  
  184.             float* c = (float*)cvGetSeqElem( circles, j );  
  185.             if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )  
  186.                 break;  
  187.         }  
  188.   
  189.         if( j < circles->total )  
  190.             continue;  
  191.           
  192.         cvStartReadSeq( nz, &reader );  
  193.           
  194.         //求所有邊緣點到當前圓心點的距離，符合條件的距離值儲存在ddata中  
  195.         for( j = k = 0; j < nz_count; j++ )  
  196.         {  
  197.             CvPoint pt;  
  198.             float _dx, _dy, _r2;  
  199.             CV_READ_SEQ_ELEM( pt, reader );  
  200.             _dx = cx - pt.x; _dy = cy - pt.y;  
  201.             _r2 = _dx*_dx + _dy*_dy;  
  202.             if(min_radius2 <= _r2 && _r2 <= max_radius2 )  
  203.             {  
  204.                 ddata[k] = _r2;  
  205.                 sort_buf[k] = k;  
  206.                 k++;  
  207.             }  
  208.         }  
  209.   
  210.         int nz_count1 = k, start_idx = nz_count1 - 1;  
  211.         if( nz_count1 == 0 )  
  212.             continue;  
  213.         dist_buf->cols = nz_count1;  
  214.         cvPow( dist_buf, dist_buf, 0.5 );  
  215.           
  216.         //經過如下處理後，以sort_buf中元素作為ddata陣列下標,ddata中的元素降序排列,   
  217.         //即ddata[sort_buf[0]]是ddata所有元素中最大的, ddata[sort_buf[nz_count1-1]]  
  218.         //是所有元素中最小的  
  219.         icvHoughSortDescent32s( &sort_buf[0], nz_count1, (int*)ddata );  
  220.           
  221.         //對所有的距離值做處理，求出最可能圓半徑值，max_count為到圓心的距離為最可能半徑值的點的個數  
  222.         dist_sum = start_dist = ddata[sort_buf[nz_count1-1]];  
  223.         for( j = nz_count1 - 2; j >= 0; j-- )  
  224.         {  
  225.             float d = ddata[sort_buf[j]];         
  226.             if( d > max_radius )  
  227.                 break;  
  228.               
  229.             if( d - start_dist > dr )  
  230.             {  
  231.                 float r_cur = ddata[sort_buf[(j + start_idx)/2]];  
  232.                 if( (start_idx - j)*r_best >= max_count*r_cur ||  
  233.                     (r_best < FLT_EPSILON && start_idx - j >= max_count) )   
  234.                 {  
  235.                     r_best = r_cur;  
  236.                     max_count = start_idx - j;  
  237.                 }  
  238.                 start_dist = d;  
  239.                 start_idx = j;  
  240.                 dist_sum = 0;  
  241.             }  
  242.             dist_sum += d;  
  243.         }  
  244.         //max_count大於閥值，表示這幾個邊緣點構成一個圓  
  245.         if( max_count > R_THRESH )  
  246.         {  
  247.             c[0] = cx;  
  248.             c[1] = cy;  
  249.             c[2] = (float)r_best;  
  250.             cvSeqPush( circles, c );  
  251.             if( circles->total > circles_max )  
  252.             return;  
  253.         }  
  254.     }  
  255. }  
頂
0
踩

相關文章