OpenCV霍夫變換系列（前篇）-經典霍夫線變換
阿新 • • 發佈：2019-01-02
前言：最近新來的的我的大學室友（現在也是我的學弟）在研究霍夫線變換，我之前只是知道這玩意可以拿來做直線檢測，並沒有深入研究，那既然提到了，還是按照我們的老規矩，原理，示例以及OpenCV這一套流程走下來。
菜鳥一枚，好多寫的不好，有點囉嗦，見諒
主要參考部落格：
還有一些博主未給出連結，也衷心表示感謝！
霍夫變換（Hough）
1.基本原理：
一條直線可由兩個點A=(X1,Y1)和B=(X2,Y2)確定(笛卡爾座標)
另一方面，
對應的變換可以通過圖形直觀表示：
變換後的空間成為霍夫空間。即：笛卡爾座標系中一條直線，對應霍夫空間的一個點 。
反過來同樣成立（霍夫空間的一條直線，對應笛卡爾座標系的一個點）：
再來看看A、B兩個點，對應霍夫空間的情形：
一步步來，再看一下三個點共線的情況：
可以看出如果笛卡爾座標系的點共線，這些點在霍夫空間對應的直線交於一點：這也是必然，共線只有一種取值可能。
如果不止一條直線呢？再看看多個點的情況（有兩條直線）：
其實（3，2）與（4，1）也可以組成直線，只不過它有兩個點確定，而圖中A、B兩點是由三條直線匯成，這也是霍夫變換的後處理的基本方式：選擇由儘可能多直線匯成的點。
看看，霍夫空間：選擇由三條交匯直線確定的點（中間圖），對應的笛卡爾座標系的直線（右圖）。
到這裡問題似乎解決了，已經完成了霍夫變換的求解，但是如果像下圖這種情況呢？
k=∞是不方便表示的，而且q怎麼取值呢，這樣不是辦法。因此考慮將笛卡爾座標系換為：極座標表示。
備註：很多文章一來就是極座標，感覺沒有這樣看著順暢，舒服。
在極座標系下，其實是一樣的：極座標的點--->霍夫空間的直線，這個地方要注意，這個必須注意，只有垂直的時候才可能是
一一對應關係。
只不過霍夫空間不再是[k,q]的引數，而是[ρ,θ]的引數，給出對比圖：
是不是就一目瞭然了？
有一個離散化的過程，本質上就是這樣（這張圖太深動形象了，怪不得後面原始碼分析就有博主總是在說格子數）
交點怎麼求解呢？細化成座標形式，取整後將交點對應的座標進行累加，最後找到數值最大的點就是求解的[ρ,θ]，也就求解出了直線 。(OpenCV不是這樣做的，但是也用到了累加的思想)
下面給出OpenCV中HoughLines的霍夫變換的直線檢測步驟:
1、對邊緣影象進行霍夫空間變換；
2、在4鄰域內找到霍夫空間變換的極大值；
3、對這些極大值安裝由大到小順序進行排序，極大值越大，越有可能是直線；
這裡說明一下極大值的含義，指的是霍夫空間中相交於某一點的曲線的數目的極大值，不要理解為數學上某一條曲線的極大值
4、輸出直線。
2. Samples
OpenCV中的霍夫直線檢測的函式為HoughLines
還有一個改進版本的HoughLinesP函式（統計概論霍夫直線檢測）
函式原型分別如下：
//函式HoughLines的原型為：
void HoughLines(InputArray image,OutputArray lines, double rho, double theta, int threshold, double srn=0,double stn=0 )
/*
image為輸入影象，要求是單通道的二值影象
lines為輸出直線向量，兩個元素的向量（ρ，θ）代表一條直線，ρ是從原點（影象的左上角）的距離，θ是直線的角度（單位是弧度），0表示垂直線，π/2表示水平線
rho為距離解析度
theta為角度解析度
threshold為閾值，在步驟2中，只有大於該值的點才有可能被當作極大值，即至少有多少條正弦曲線交於一點才被認為是直線
srn和stn在多尺度霍夫變換的時候才會使用，在這裡我們只研究經典霍夫變換的原始碼
*///函式HoughLinesP的原型為：
void HoughLinesP(InputArray image, OutputArray lines, double rho, double theta, int threshold, double minLineLength=0, double maxLineGap=0 )
/*
第一個引數，InputArray型別的image，輸入影象，即源影象，需為8位的單通道二進位制影象，可以將任意的源圖載入進來後由函式修改成此格式後，再填在這裡。
第二個引數，InputArray型別的lines，經過呼叫HoughLinesP函式後後儲存了檢測到的線條的輸出向量，每一條線由具有四個元素的向量(x_1,y_1, x_2, y_2） 表示，其中，(x_1, y_1)和(x_2, y_2) 是是每個檢測到的線段的結束點。
第三個引數，double型別的rho，以畫素為單位的距離精度。另一種形容方式是直線搜尋時的進步尺寸的單位半徑。
第四個引數，double型別的theta，以弧度為單位的角度精度。另一種形容方式是直線搜尋時的進步尺寸的單位角度。
第五個引數，int型別的threshold，累加平面的閾值引數，即識別某部分為圖中的一條直線時它在累加平面中必須達到的值。大於閾值threshold的線段才可以被檢測通過並返回到結果中。
第六個引數，double型別的minLineLength，有預設值0，表示最低線段的長度，比這個設定引數短的線段就不能被顯現出來。
第七個引數，double型別的maxLineGap，有預設值0，允許將同一行點與點之間連線起來的最大的距離。
*/由於函式HoughLines只輸出直線所對應的角度和距離，所以在進行直線檢測的時候還要把其轉換為直角座標系下的資料，另外輸入影象還必須是邊緣影象，下面就是具體的例項：
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>
#include <iostream>
#include <Windows.h>
using namespace cv;
using namespace std;
Mat src, dst, cdst;
const char* filename = "1.jpg";
const char *winname="hough";
const char *trackbarname="1.houghlines\n2.houghlinep\n3.houghcircle";
int savevalue=100,houghtype;
const int maxthreshold=150;
void help()
{
cout << "\nThis program demonstrates line finding with the Hough transform.\n";
}
void choisehoughlines()
{
vector<Vec2f> lines;
Canny(src, dst,50, 200, 3);
cvtColor(dst, cdst, CV_GRAY2BGR);
HoughLines(dst, lines, 1, CV_PI/180, savevalue+10, 0, 0 );
for( size_t i = 0; i < lines.size(); i++ )
{
float rho = lines[i][0], theta = lines[i][1];
Point pt1, pt2;
double a = cos(theta), b = sin(theta);
double x0 = a*rho, y0 = b*rho;
pt1.x = cvRound(x0 + 1000*(-b));
pt1.y = cvRound(y0 + 1000*(a));
pt2.x = cvRound(x0 - 1000*(-b));
pt2.y = cvRound(y0 - 1000*(a));
line( cdst, pt1, pt2, Scalar(0,0,255), 1, CV_AA);
float angle = theta / CV_PI *180;
//cout<<"line "<< i << ": "<<"rho: "<<rho<<"theta: "<<theta
//<<"angle: "<< angle<<endl;
}
imshow(winname, cdst);
}
void choisehoughlinep()
{
vector<Vec4i> lines;
Canny(src, dst,50, 200, 3);
cvtColor(dst, cdst, CV_GRAY2BGR);
HoughLinesP(dst, lines, 1, CV_PI/180, savevalue+10, savevalue+10, 10 );
for( size_t i = 0; i < lines.size(); i++ )
{
Vec4i l = lines[i];
line( cdst, Point(l[0], l[1]), Point(l[2], l[3]), Scalar(0,0,255), 1, CV_AA);
}
imshow(winname, cdst);
}
void choisehoughcircle()
{
vector<Vec3f> circles;
cvtColor(src, cdst, CV_GRAY2BGR);
/// Apply the Hough Transform to find the circles
HoughCircles( src, circles, CV_HOUGH_GRADIENT, 1, dst.rows/10, 200, savevalue+10, 0, 0 );
/// Draw the circles detected
printf("%d",circles.size());
for( size_t i = 0; i < circles.size(); i++ )
{
Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
int radius = cvRound(circles[i][2]);
// circle center
circle( cdst, center, 3, Scalar(0,255,0), -1, 8, 0 );
// circle outline
circle( cdst, center, radius, Scalar(0,0,255), 3, 8, 0 );
}
imshow(winname, cdst);
}
void choice(int,void *)
{
switch (houghtype)
{
case 0:choisehoughlines();break;
case 1:choisehoughlinep();break;
//case 2:choisehoughcircle();break;
}
}
int main(int argc, char** argv)
{
SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN |FOREGROUND_INTENSITY);
//system("color 3E");
help();
src = imread(filename, 0);
if(src.empty())
{
help();
cout << "can not open " << filename << endl;
return -1;
}
namedWindow(winname);
createTrackbar(trackbarname,winname,&houghtype,2,choice);
createTrackbar("thresholdvalue",winname,&savevalue,maxthreshold,choice);
imshow("source", src);
choice(0,0);
while(char(waitKey())!= 'q');
return 0;
}
這是一個綜合示例，我把霍夫圓變換註釋掉了，那個放在下一篇單獨講解，其中輸出經典霍夫線變換的角度的程式碼我也註釋掉了，
需要的可以自己加上。
對於經典的線變換計算得到的兩點的座標為（ρcosθ-1000sinθ，ρsinθ+1000cosθ），（ρcosθ+1000sinθ，ρsinθ-1000cosθ），我個人覺得
計算過程如下：
結論：
houghlines的計算效率比較低O(n*n*m)，耗時較長，而且沒有檢測出直線的端點。
統計概論霍夫直線檢測houghlinesP是一個改進，不僅執行效率較高，而且能檢測到直線的兩個端點。
思想：
先隨機檢測出一部分直線，然後將直線上點的排查掉，再進行其他直線的檢測
1. 首先僅統計影象中非零點的個數，對於已經確認是某條直線上的點就不再變換了。
2. 對所以有非零點逐個變換到霍夫空間
a. 並累加到霍夫統計表(影象)中，並統計最大值
b. 最大值與閾值比較，小於閾值，則繼續下一個點的變換
c. 若大於閾值，則有一個新的直線段要產生了
d. 計算直線上線段的端點、長度，如果符合條件，則儲存此線段，並mark這個線段上的點不參與其他線段檢測的變換。
3.原始碼分析
HoughLines函式是在sources/modules/imgproc/src/hough.cpp檔案中定義的：
void cv::HoughLines( InputArray _image, OutputArray _lines,
double rho, double theta, int threshold,
double srn, double stn )
{
//申請一段記憶體，用於儲存霍夫變換後所檢測到的直線
Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
//提取輸入影象矩陣
Mat image = _image.getMat();
CvMat c_image = image;
//呼叫由C語音寫出的霍夫變換的函式
CvSeq* seq = cvHoughLines2( &c_image, storage, srn == 0 && stn == 0 ?
CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,
rho, theta, threshold, srn, stn );
//把由cvHoughLines2函式得到的霍夫變換直線序列轉換為陣列
seqToMat(seq, _lines);
}CV_IMPL CvSeq*
cvHoughLines2( CvArr* src_image, void* lineStorage, int method,
double rho, double theta, int threshold,
double param1, double param2 )
{
CvSeq* result = 0;
CvMat stub, *img = (CvMat*)src_image;
CvMat* mat = 0;
CvSeq* lines = 0;
CvSeq lines_header;
CvSeqBlock lines_block;
int lineType, elemSize;
int linesMax = INT_MAX; //輸出最多直線的數量，設為無窮多
int iparam1, iparam2;
img = cvGetMat( img, &stub );
//確保輸入影象是8位單通道
if( !CV_IS_MASK_ARR(img))
CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
//記憶體空間申請成功，以備輸出之用
if( !lineStorage )
CV_Error( CV_StsNullPtr, "NULL destination" );
//確保rho，theta，threshold這三個引數大於0
if( rho <= 0 || theta <= 0 || threshold <= 0 )
CV_Error( CV_StsOutOfRange, "rho, theta and threshold must be positive" );
if( method != CV_HOUGH_PROBABILISTIC ) //經典霍夫變換和多尺度霍夫變換
{
//輸出直線的型別為32位浮點雙通道，即ρ和θ兩個浮點型變數
lineType = CV_32FC2;
elemSize = sizeof(float)*2;
}
else //概率霍夫變換
{
//輸出直線的型別為32位有符號整型4通道，即兩個畫素點的4個座標
lineType = CV_32SC4;
elemSize = sizeof(int)*4;
}
//判斷lineStorage的型別，經分析lineStorage只可能是STORAGE型別
if( CV_IS_STORAGE( lineStorage ))
{
//在lineStorage記憶體中建立一個序列，用於儲存霍夫變換的直線，lines為該序列的指標
lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage );
}
else if( CV_IS_MAT( lineStorage ))
{
mat = (CvMat*)lineStorage;
if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )
CV_Error( CV_StsBadArg,
"The destination matrix should be continuous and have a single row or a single column" );
if( CV_MAT_TYPE( mat->type ) != lineType )
CV_Error( CV_StsBadArg,
"The destination matrix data type is inappropriate, see the manual" );
lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,
mat->rows + mat->cols - 1, &lines_header, &lines_block );
linesMax = lines->total;
cvClearSeq( lines );
}
else
CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );
iparam1 = cvRound(param1);
iparam2 = cvRound(param2);
switch( method )
{
case CV_HOUGH_STANDARD: //經典霍夫變換
icvHoughLinesStandard( img, (float)rho,
(float)theta, threshold, lines, linesMax );
break;
case CV_HOUGH_MULTI_SCALE: //多尺度霍夫變換
icvHoughLinesSDiv( img, (float)rho, (float)theta,
threshold, iparam1, iparam2, lines, linesMax );
break;
case CV_HOUGH_PROBABILISTIC: //概率霍夫變換
icvHoughLinesProbabilistic( img, (float)rho, (float)theta,
threshold, iparam1, iparam2, lines, linesMax );
break;
default:
CV_Error( CV_StsBadArg, "Unrecognized method id" );
}
//在前面判斷lineStorage型別時，已經確定它是STORAGE型別，因此沒有進入else if內，也就是沒有對mat賦值，所以在這裡進入的是else
if( mat )
{
if( mat->cols > mat->rows )
mat->cols = lines->total;
else
mat->rows = lines->total;
}
else
result = lines;
//返回lines序列的指標，該序列儲存有霍夫變換的直線
return result;
}
並且其中有一個地方還是不明白,已經標註出來，希望有人解釋。
static void icvHoughLinesStandard( const CvMat* img, float rho, float theta,
int threshold, CvSeq *lines, int lineMax)
{
cv::AutoBuffer<int> _accum, _sort_buf;
cv::AutoBuffer<float> _tabSin, _tabCos;
const uchar* image;
int step, width, height;
int numangle, numrho;
int total=0;
float ang;
int r, n;
int i, j;
float irho = 1 / rho;
double scale;
//再次確保輸入影象的正確性
CV_Assert( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1);
image = img->data.ptr; //得到影象的指標
step = img->step; //得到影象的步長
width = img->cols; //得到影象的寬
height = img->rows; //得到影象的高
//由角度和距離的解析度得到角度和距離的數量，即霍夫變換后角度和距離的個數
numangle = cvRound(CV_PI / theta); // 霍夫空間,角度方向的大小
numrho = cvRound(((width + height)*2 + 1) / rho); //r的空間範圍
/*
allocator類是一個模板類，定義在標頭檔案memory中，用於記憶體的分配、釋放、管理，它幫助我們將記憶體分配和物件構造分離開來。具體地說，allocator類將記憶體的分配和物件的構造解耦，分別用allocate和construct兩個函式完成，同樣將記憶體的釋放和物件的析構銷燬解耦，分別用deallocate和destroy函式完成。
*/
//為累加器陣列分配記憶體空間
//該累加器陣列其實就是霍夫空間，它是用一維陣列表示二維空間
_accum.allocate((numangle+2) * (numrho + 2));
//為排序陣列分配記憶體空間
_sort_buf.allocate(numangle * numrho);
//為正弦和餘弦列表分配記憶體空間
_tabSin.allocate(numangle);
_tabCos.allocate(numangle);
//分別定義上述記憶體空間的地址指標
int *accum = _accum, *sort_buf = _sort_buf;
int *tabSin = _tabSin, *tabCos = _tabCos;
//累加器陣列清零
memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );
//為避免重複運算，事先計算好sinθi/ρ和cosθi/ρ
for( ang = 0, n = 0; n < numangle; ang += theta, n++ ) //計算正弦曲線的準備工作,查表
{
tabSin[n] = (float)(sin(ang) * irho);
tabCos[n] = (float)(cos(ang) * irho);
}
//stage 1. fill accumulator
//執行步驟1，逐點進行霍夫空間變換，並把結果放入累加器陣列內
for( i = 0 ; i < height; i++)
for( j = 0; j < width; j++)
{
//只對影象的非零值處理，即只對影象的邊緣畫素進行霍夫變換
if( image[i * step + j] != 0 ) //將每個非零點,轉換為霍夫空間的離散正弦曲線,並統計。
for( n = 0; n < numangle; n++ )
{
//根據公式: ρ = xcosθ + ysinθ
//cvRound()函式：四捨五入
r = cvRound( j * tabCos[n] + i * tabSin[n])
//numrho是ρ的最大值，或者說最大取值範圍
r += (numrho - 1) / 2; //這一步是真的看不太明白？？？
//(另一位博主的解釋)哈哈,這裡讓我想了好久,為什麼這樣安排呢? 可能是因為θ、ρ正負問題 ，但我感覺解釋不通
//r表示的是距離，n表示的是角點，在累加器內找到它們所對應的位置（即霍夫空間內的位置），其值加1
accum[(n+1) * (numrho+2)+ r + 1]++;
/*
最初我也是下面這樣理解的，覺得比較好理解，直觀，但是是一維陣列
哪來的行與列額!
n+1是為了第一行空出來
// numrho+2 是總共的列數,這裡實際意思應該是半徑的所有可能取值，加2是為了防止越界，但是說成列數我覺得也沒錯，並且我覺得這樣比較好理解
//r+1 是為了把第一列空出來
//因為程式後面要比較前一行與前一列的資料, 這樣省的越界
因此空出首行和首列*/
}
}
// stage 2. find local maximums
// 執行步驟2，找到區域性極大值，即非極大值抑制
// 霍夫空間,區域性最大點,採用四領域判斷,比較。(也可以使8鄰域或者更大的方式),如果不判斷區域性最大值,同時選用次大值與最大值,就可能會是兩個相鄰的直線,但實際是一條直線。選用最大值,也是去除離散的近似計算帶來的誤差,或合併近似曲線。
for( r = 0 ; r < numrho; r++ )
for( n = 0; n < numangle; n++ )
{
//得到當前值在累加器陣列的位置
int base = (n+1)*(numrho+2) + r + 1;
//得到計數值，並以它為基準，看看它是不是區域性極大值
if( accum[base] > threshold &&
accum[base] > accum[base - 1] && accum[base] >= accum[base+1] &&
accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2)
//把極大值位置存入排序陣列內——sort_buf
sort_buf[total++] = base;
}
//stage 3. sort the detected lines by accumulator value
//執行步驟3，對儲存在sort_buf陣列內的累加器的資料按由大到小的順序進行排序
icvHoughSortDescent32s( sort_buf, total, accum );
/*OpenCV中自帶了一個排序函式，名稱為：
void icvHoughSortDescent32s(int *sequence , int sum , int*data)，引數解釋：
第三個引數：陣列的首地址
第二個引數：要排序的資料總數目
第一個引數：此陣列存放data中要參與排序的資料的序號
而且這個排序演算法改變的只是sequence[]陣列中的元素，源資料data[]未發生絲毫改變。
*/
// stage 4. store the first min(total, linesMax ) lines to the output buffer,輸出直線
lineMax = MIN(lineMax, total); //linesMax是輸入引數，表示最多輸出幾條直線
//事先定義一個尺度
scale = 1./(numrho+2);
for(i=0; i<linesMax; i++) // 依據霍夫空間解析度,計算直線的實際r,theta引數
{
//CvLinePolar 直線的資料結構
//CvLinePolar結構在該檔案的前面被定義
CvLinePolar line;
//idx為極大值在累加器陣列的位置
int idx = sort_buf[i]; //找到索引(下標)
//分離出該極大值在霍夫空間中的位置
int n = cvFloor(idx*scale) - 1; //向下取整
int r = idx - (n+1)*(numrho+2) - 1;
//最終得到極大值所對應的角度和距離
line.rho = (r - (numrho - 1)*0.5f)*rho;
line.angle = n * theta;
//儲存到序列內
cvSeqPush( lines, &line ); //用序列存放多條直線
}
}