RobHess的SIFT代码解析之RANSAC
平臺：win10 x64 +VS 2015專業版 +opencv-2.4.11 + gtk_-bundle_2.24.10_win32
主要參考：1.代碼：RobHess的SIFT源碼：SIFT+KD樹+BBF算法+RANSAC算法
2.書：王永明 王貴錦 《圖像局部不變性特征與描述》
?
?
RobHess的SIFT源碼中的幾個文件說明？
RobHess的SIFT源碼分析：xform.h和xform.c文件
這兩個文件中實現了RANSAC算法(RANdom SAmple Consensus 隨機抽樣一致)。
RANSAC算法可用來篩選兩個圖像間的SIFT特征匹配并計算變換矩陣。
?
請在看這個之前先看完以下內容前五個：
?
SIFT四步驟和特征匹配及篩選：
步驟一：建立尺度空間，即建立高斯差分(DoG)金字塔dog_pyr
步驟二：在尺度空間中檢測極值點，并進行精確定位和篩選創建默認大小的內存存儲器
步驟三：特征點方向賦值，完成此步驟后，每個特征點有三個信息：位置、尺度、方向
步驟四：計算特征描述子
SIFT后特征匹配：KD樹+BBF算法
SIFT后特征匹配后錯誤點篩選：RANSAC算法
?
?
SIFT后特征匹配后錯誤點篩選：RANSAC算法
問題及解答：
?
（1）問題描述：利用RANSAC算法篩選SIFT特征匹配的主要流程是什么？
答：
(1) 從樣本集中隨機抽選一個RANSAC樣本，即4個匹配點對
(2) 根據這4個匹配點對計算變換矩陣M
(3) 根據樣本集，變換矩陣M，和誤差度量函數計算滿足當前變換矩陣的一致集consensus，并返回一致集中元素個數
(4) 根據當前一致集中元素個數判斷是否最優(最大)一致集，若是則更新當前最優一致集
(5) 更新當前錯誤概率p，若p仍大于允許的最小錯誤概率則重復(1)至(4)繼續迭代，直到當前錯誤概率p小于最小錯誤概率
RANSAC算法的迭代終止條件是當前錯誤概率p小于允許的最小錯誤概率p_badxform，一般傳入的值為0.01，即要求正確率達到99%
當前錯誤概率p的計算公式為：p=( 1 - in_fracm)^k，
其中in_frac是當前最優一致集中元素(內點)個數占樣本總數的百分比，表征了當前最優一致集的優劣；
m是計算變換矩陣需要的最小特征點對個數，是固定值，一般是4；k是迭代次數。
所以，內點所占比例in_frac越大，錯誤概率越小；迭代次數k越大，錯誤概率越小，用來保證收斂。
?
（2）問題描述：RobHess的源碼如何實現RANSAC，大概思路是這樣的？
答：（2.1）代碼及說明：——match.c的main函數里
/*下面的觀察看看RANSAC功能如何運作
???? 注意上面的這一行：
???? feat1[i].fwd_match = nbrs[0];
???? 對于RANSAC功能起作用很重要。
? */
? //計算變換參數，利用RANSAC算法篩選匹配點,計算變換矩陣H
? {
??? CvMat* H; //RANSAC算法求出的變換矩陣
??? IplImage* xformed; //xformed臨時拼接圖，即只將圖2變換后的圖
???????? //無論img1和img2的上下順序，H永遠是將feat1中的特征點變換為其匹配點，即將img1中的點變換為img2中的對應點
??? H = ransac_xform( feat1, n1, FEATURE_FWD_MATCH, lsq_homog, 4, 0.01,
?????????????????? ????? homog_xfer_err, 3.0, NULL, NULL );
//feat1：特征點數組；n1：特征點個數；FEATURE_FWD_MATCH：特征點類型，三者之一；lsq_homog：函數指針（2.2.5）；4：在lsq_homog中計算變換矩陣需要最小特征點對個數；0.01：允許的錯誤概率，當前計算出的錯誤概率小于0.01時停止；homog_xfer_err：函數指針（2.2.6.1）；3.0：容錯度（此范圍內認為是內點）
???????? //若能成功計算出變換矩陣，即兩幅圖中有共同區域
???????? if( H )
????? {
??? //為拼接結果圖xformed分配空間,高度為圖2高度
???????? xformed = cvCreateImage( cvGetSize( img2 ), IPL_DEPTH_8U, 3 );
???????? //用變換矩陣H對右圖img1做投影變換(變換后會有坐標右移)，結果放到xformed中
???????? cvWarpPerspective( img1, xformed, H,
??????????????????????????? ?? CV_INTER_LINEAR + CV_WARP_FILL_OUTLIERS,
??????????????????????????? ?? cvScalarAll( 0 ) );
???????? cvNamedWindow( "Xformed", 1 ); //創建窗口
???????? cvShowImage( "Xformed", xformed ); //顯示臨時圖,即只將圖2變換后的圖
???????? cvWaitKey( 0 );
???????? cvReleaseImage( &xformed );
???????? cvReleaseMat( &H ); //釋放變換矩陣H
????? }
}
?
（2.2）ransac_xform代碼及說明：
答：
/*利用RANSAC算法進行特征點篩選，計算出最佳匹配的變換矩陣
參數：
features：特征點數組，只有當mtype類型的匹配點存在時才被用來進行單應性計算
n：特征點個數
mtype：決定使用每個特征點的哪個匹配域進行變換矩陣的計算，應該是FEATURE_MDL_MATCH，
FEATURE_BCK_MATCH，FEATURE_MDL_MATCH中的一個。若是FEATURE_MDL_MATCH，
對應的匹配點對坐標是每個特征點的img_pt域和其匹配點的mdl_pt域，
否則，對應的匹配點對是每個特征點的img_pt域和其匹配點的img_pt域。
xform_fn：函數指針，指向根據輸入的點對進行變換矩陣計算的函數，一般傳入lsq_homog()函數
m：在函數xform_fn中計算變換矩陣需要的最小特征點對個數
p_badxform：允許的錯誤概率，即允許RANSAC算法計算出的變換矩陣錯誤的概率，當前計算出的模型的錯誤概率小于p_badxform時迭代停止
err_fn：函數指針，對于給定的變換矩陣，計算推定的匹配點對之間的變換誤差，一般傳入homog_xfer_err()函數
err_tol：容錯度，對于給定的變換矩陣，在此范圍內的點對被認為是內點
inliers：輸出參數，指針數組，指向計算出的最終的內點集合，若為空，表示沒計算出符合要求的一致集
此數組的內存將在本函數中被分配，使用完后必須在調用出釋放：free(*inliers)
n_in：輸出參數，最終計算出的內點的數目
返回值：RANSAC算法計算出的變換矩陣，若為空，表示出錯或無法計算出可接受矩陣
CvMat* ransac_xform( struct feature* features, int n, int mtype,
?????????????????????????????????? ransac_xform_fn xform_fn, int m, double p_badxform,
?????????????????????????????????? ransac_err_fn err_fn, double err_tol,
struct feature*** inliers, int* n_in )
{
//matched：所有具有mtype類型匹配點的特征點的數組，也就是樣本集
//sample：單個樣本，即4個特征點的數組
//consensus：當前一致集；
//consensus_max：當前最大一致集(即當前的最好結果的一致集)
struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
struct ransac_data* rdata;//每個特征點的feature_data域的ransac數據指針
CvPoint2D64f* pts, * mpts;//每個樣本對應的兩個坐標數組：特征點坐標數組pts和匹配點坐標數組mpts
CvMat* M = NULL;//當前變換矩陣
//p：當前計算出的模型的錯誤概率，當p小于p_badxform時迭代停止
//in_frac：內點數目占樣本總數目的百分比
double p, in_frac = RANSAC_INLIER_FRAC_EST;
//nm：輸入的特征點數組中具有mtype類型匹配點的特征點個數
//in：當前一致集中元素個數
//in_min：一致集中元素個數允許的最小值，保證RANSAC最終計算出的轉換矩陣錯誤的概率小于p_badxform所需的最小內點數目
//in_max：當前最優一致集(最大一致集)中元素的個數
//k：迭代次數，與計算當前模型的錯誤概率有關
int i, nm, in, in_min, in_max = 0, k = 0;
//找到特征點數組features中所有具有mtype類型匹配點的特征點，放到matched數組(樣本集)中，返回值nm是matched數組的元素個數
nm = get_matched_features( features, n, mtype, &matched );
//若找到的具有匹配點的特征點個數小于計算變換矩陣需要的最小特征點對個數，出錯
if( nm < m )
{?? //出錯處理，特征點對個數不足
??? fprintf( stderr, "Warning: not enough matches to compute xform, %s" \
??????? " line %d\n", __FILE__, __LINE__ );
????????????? goto end;
}
srand( time(NULL) );//初始化隨機數生成器
//計算保證RANSAC最終計算出的轉換矩陣錯誤的概率小于p_badxform所需的最小內點數目
in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
//當前計算出的模型的錯誤概率,內點所占比例in_frac越大，錯誤概率越小；迭代次數k越大，錯誤概率越小
p = pow( 1.0 - pow( in_frac, m ), k );
i = 0;
//當前錯誤概率大于輸入的允許錯誤概率p_badxform，繼續迭代
while( p > p_badxform )
?{
//從樣本集matched中隨機抽選一個RANSAC樣本(即一個4個特征點的數組)，放到樣本變量sample中
?? sample = draw_ransac_sample( matched, nm, m );
//從樣本中獲取特征點和其對應匹配點的二維坐標，分別放到輸出參數pts和mpts中
?? extract_corresp_pts( sample, m, mtype, &pts, &mpts );
//調用參數中傳入的函數xform_fn，計算將m個點的數組pts變換為mpts的矩陣，返回變換矩陣給M
M = xform_fn( pts, mpts, m );//一般傳入lsq_homog()函數
if( ! M )//出錯判斷
?? goto iteration_end;
//給定特征點集，變換矩陣，誤差函數，計算出當前一致集consensus，返回一致集中元素個數給in
????????????? in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
?
//若當前一致集大于歷史最優一致集，即當前一致集為最優，則更新最優一致集consensus_max
?????? if( in > in_max )
??????? {
if( consensus_max )//若之前有最優值，釋放其空間
????? ? ?? free( consensus_max );
consensus_max = consensus;//更新最優一致集
in_max = in;//更新最優一致集中元素個數
in_frac = (double)in_max / nm;//最優一致集中元素個數占樣本總個數的百分比
?????? }
else//若當前一致集小于歷史最優一致集，釋放當前一致集
?????????? free( consensus );
?????? cvReleaseMat( &M );
?
iteration_end:
????? release_mem( pts, mpts, sample );
????? p = pow( 1.0 - pow( in_frac, m ), ++k );//更新當前錯誤概率
?????? }
//根據最優一致集計算最終的變換矩陣
//若最優一致集中元素個數大于最低標準，即符合要求
?????? if( in_max >= in_min )
?????? {
//從最優一致集中獲取特征點和其對應匹配點的二維坐標，分別放到輸出參數pts和mpts中
?????? extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
//調用參數中傳入的函數xform_fn，計算將in_max個點的數組pts變換為mpts的矩陣，返回變換矩陣給M
????????????? M = xform_fn( pts, mpts, in_max );
/***********下面會再進行一次迭代**********/
//根據變換矩陣M從樣本集matched中計算出一致集consensus，返回一致集中元素個數給in
?????? in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
????????????? cvReleaseMat( &M );
????????????? release_mem( pts, mpts, consensus_max );
//從一致集中獲取特征點和其對應匹配點的二維坐標，分別放到輸出參數pts和mpts中
????????????? extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
//調用參數中傳入的函數xform_fn，計算將in個點的數組pts變換為mpts的矩陣，返回變換矩陣給M
????????????? M = xform_fn( pts, mpts, in );
????????????? if( inliers )
????????????? {
*inliers = consensus;//將最優一致集賦值給輸出參數：inliers，即內點集合
???????????? consensus = NULL;
????????????? }
????????????? if( n_in )
*n_in = in;//將最優一致集中元素個數賦值給輸出參數：n_in，即內點個數
????????????? release_mem( pts, mpts, consensus );
?????? }
?????? else if( consensus_max )
{?? //沒有計算出符合要求的一致集
????? if( inliers )
???????? *inliers = NULL;
??????? if( n_in )
??????????? *n_in = 0;
????????????? free( consensus_max );
}
//RANSAC算法結束：恢復特征點中被更改的數據域feature_data，并返回變換矩陣M
end:
?????? for( i = 0; i < nm; i++ )
?????? {
//利用宏feat_ransac_data來提取matched[i]中的feature_data成員并轉換為ransac_data格式的指針
??????? rdata = feat_ransac_data( matched[i] );
//恢復feature_data成員的以前的值
?????? matched[i]->feature_data = rdata->orig_feat_data;
free( rdata );//釋放內存
?????? }
?????? free( matched );
return M;//返回求出的變換矩陣M
}
?
（2.2.1）get_matched_features代碼及說明：
答：
?/*找到所有具有mtype類型匹配點的特征點，將他們的指針存儲在matched數組中，
? 并初始化matched數組中每個特征點的feature_data域為ransac_data類型的數據指針
參數：
features：特征點數組
n：特征點個數
mtype：匹配類型
matched：輸出參數，含有mtype類型匹配點的特征點的指針數組
返回值：matched數組中特征點的個數
static int get_matched_features( struct feature* features, int n, int mtype, struct feature*** matched )
{
struct feature** _matched;//輸出數組，具有mtype類型匹配點的特征點指針數組
struct ransac_data* rdata;//ransac_data類型數據指針
?? ? ? int i, m = 0;
?
?? ? _matched = calloc( n, sizeof( struct feature* ) );
?
?? //遍歷輸入的特征點數組
?? ?for( i = 0; i < n; i++ )
{?? //找第i個特征點的mtype類型匹配點，若能正確找到表明此特征點有mtype類型的匹配點，則將其放入_matched數組
if( get_match( features + i, mtype ) )
?? ??? ?{
rdata = malloc( sizeof( struct ransac_data ) );//為ransac_data結構分配空間
memset( rdata, 0, sizeof( struct ransac_data ) );//清零
rdata->orig_feat_data = features[i].feature_data;//保存第i個特征點的feature_data域之前的值
_matched[m] = features + i;//放到_matched數組
_matched[m]->feature_data = rdata;//其feature_data域賦值為ransac_data類型數據的指針
m++;//_matched數組中元素個數
?? ??? ?}
}
*matched = _matched;//輸出參數賦值
return m;//返回值，元素個數
}
問題1：feature結構體里的feature_data在RANSAC算法里是什么？和ransac_data有什么關系？
?答：在SIFT極值點檢測中，是detection_data結構的指針；在k-d樹搜索中，是bbf_data結構的指針；在RANSAC算法中，是ransac_data結構的指針。而在RANSAC算法里ransac_data包括兩部分：
struct ransac_data
{
void* orig_feat_data; //保存此特征點的feature_data域的以前的值
int sampled; //標識位，值為1標識此特征點是否被選為樣本
};
?
?
問題2：RANSAC算法：從樣本集matched中隨機抽選一個RANSAC樣本(即一個4個特征點的數組)，放到樣本變量sample中，四個點怎么隨機選？
答：采用srand和rand()配合使用,而rand()在函數draw_ransac_sample中，使用x = rand() % n;產生隨機下標，而把該點的feature_data成員的sampled （標志位）置為 1;表明此點已選，在for循環找的m（m=4）個點不能重合。
函數說明：srand函數是隨機數發生器的初始化函數。原型：void srand(unsigned int seed);srand和rand()配合使用產生偽隨機數序列。
參考：百度百科——https://baike.baidu.com/item/srand/796881?fr=aladdin
?
問題3：當前計算出的模型的錯誤概率p怎么計算？
答：由公式計算：p = pow( 1.0 - pow( in_frac, m ), k );初始為p=(1-0.25^4)^0 = 1
此公式在書P166（7-4），in_frac為公式中的w，m為N，為迭代次數
?
（2.2.1.1）get_matched_features代碼及說明：
答：
/*根據給定的匹配類型，返回輸入特征點feat的匹配點
參數：
feat：輸入特征點
mtype：匹配類型，是FEATURE_FWD_MATCH，FEATURE_BCK_MATCH，FEATURE_MDL_MATCH之一
返回值：feat的匹配點的指針，若為空表示mtype參數有誤
static __inline struct feature* get_match( struct feature* feat, int mtype )
{
//FEATURE_MDL_MATCH：表明feature結構中的mdl_match域是對應的匹配點
?????? if( mtype == FEATURE_MDL_MATCH )
????????????? return feat->mdl_match;
//FEATURE_BCK_MATCH：表明feature結構中的bck_match域是對應的匹配點
?????? if( mtype == FEATURE_BCK_MATCH )
????????????? return feat->bck_match;
//FEATURE_FWD_MATCH：表明feature結構中的fwd_match域是對應的匹配點
?????? if( mtype == FEATURE_FWD_MATCH )
????????????? return feat->fwd_match;
?????? return NULL;
}
?
（2.2.2）calc_min_inliers代碼及說明：
答：
/*計算保證RANSAC最終計算出的轉換矩陣錯誤的概率小于p_badxform所需的最小內點數目
參數：
n：推定的匹配點對的個數
m：計算模型所需的最小點對個數，初始宏定義值4
p_badsupp：概率，錯誤模型被一個匹配點對支持的概率，初始宏定義值0.1
p_badxform：概率，最終計算出的轉換矩陣是錯誤的的概率，初始宏定義值0.01
返回值：保證RANSAC最終計算出的轉換矩陣錯誤的概率小于p_badxform所需的最小內點數目
static int calc_min_inliers( int n, int m, double p_badsupp, double p_badxform )
{
//根據論文：Chum, O. and Matas, J.? Matching with PROSAC -- Progressive Sample Consensus中的一個公式計算，但是不知道為什么這么計算！
?????? double pi, sum;
?????? int i, j;
?
?????? for( j = m+1; j <= n; j++ )
?????? {
????????????? sum = 0;
????????????? for( i = j; i <= n; i++ )
????????????? {
???????????????????? pi = ( i - m ) * log( p_badsupp ) + ( n - i + m ) * log( 1.0 - p_badsupp ) +
??????????????????????????? log_factorial( n - m ) - log_factorial( i - m ) -
??????????????????????????? log_factorial( n - i );
???????????????????? /*
???????????????????? ?* Last three terms above are equivalent to log( n-m choose i-m )
???????????????????? ?*/
???????????????????? sum += exp( pi );
????????????? }
????????????? if( sum < p_badxform )
???????????????????? break;
?????? }
?????? return j;
}
問題：論文：Chum, O. and Matas, J.? Matching with PROSAC -- Progressive Sample Consensus中公式？代碼如何計算？
?答：公式如下：
代碼中兩邊取對數了，具體看我的演算紙：
?
（2.2.2.1）log_factorial代碼及說明：
答：
//計算n的階乘的自然對數 log（n!）
static __inline double log_factorial( int n )
{
? double f = 0;
? int i;
? for( i = 1; i <= n; i++ )
f += log( i );
?
? return f;
}
具體看我的演算紙：
?
（2.2.3）draw_ransac_sample代碼及說明：
答：
?
/*從給定的特征點集中隨機抽選一個RANSAC樣本(即一個4個特征點的數組)
參數：
features：作為樣本集的特征點數組
n：features中元素個數
m：單個樣本的尺寸，這里是4(至少需要4個點來計算變換矩陣)
返回值：一個指針數組，其元素指向被選為樣本的特征點，被選為樣本的特征點的feature_data域的sampled被設為1
static struct feature** draw_ransac_sample( struct feature** features, int n, int m )
{
struct feature** sample, * feat;//sample：被選為樣本的點的數組
?????? struct ransac_data* rdata;
?????? int i, x;
?
//將所有特征點的feature_data域的sampled值都初始化為0，即未被選為樣本點
?????? for( i = 0; i < n; i++ )
?????? {
//利用宏feat_ransac_data來提取參數中的feature_data成員并轉換為ransac_data格式的指針
rdata = feat_ransac_data( features[i] );
rdata->sampled = 0;//sampled值設為0
?????? }
?
sample = calloc( m, sizeof( struct feature* ) );//為樣本數組分配空間
?
//隨機抽取m個特征點作為一個樣本，將其指針保存在sample數組中
?????? for( i = 0; i < m; i++ )
?????? {
????????????? do
????????????? {
x = rand() % n;//隨機下標
???????????????????? feat = features[x];
rdata = feat_ransac_data( feat );//獲得feature_data成員并轉換為ransac_data格式的指針
????????????? }
while( rdata->sampled );//若抽取的特征點的sampled值為1，繼續選取；否則停止，將其作為樣本中的一個點
//sampled為標志位，初始為0，但（抽后邊m-1次）之前已經抽過，為避免重復抽到，這樣保證抽到了相異的四個點
sample[i] = feat;//放入sample數組
rdata->sampled = 1;//該點的feature_data成員的sampled域值設為1
?????? }
?
return sample;//返回隨機選取的樣本
}
?
問題1：如何保證抽到的四個點隨機且不同？
答：隨機通過 srand和rand()函數配合使用，不同，通過標志位sampled
?
（2.2.4）extract_corresp_pts代碼及說明：
答：
/*從特征點數組中獲取特征點和其對應匹配點的二維坐標，分別放到輸出參數pts和mpts中
參數：
features：特征點數組，將從其中抽取坐標點和其匹配點，此數組中所有特征點都應具有mtype類型的匹配點
n：feantures中特征點個數，n=4
mtype：匹配類型，若是FEATURE_MDL_MATCH，對應的匹配點對坐標是每個特征點的img_pt域和其匹配點的mdl_pt域， 否則，對應的匹配點對是每個特征點的img_pt域和其匹配點的img_pt域。三者之一
pts：輸出參數，從特征點數組features中獲取的二維坐標數組
mpts：輸出參數，從特征點數組features的對應匹配點中獲取的二維坐標數組
static void extract_corresp_pts( struct feature** features, int n, int mtype,
??????????????????????????????????????????????????????? ?CvPoint2D64f** pts, CvPoint2D64f** mpts )
{
struct feature* match;//每個特征點對應的匹配點
?????? CvPoint2D64f* _pts, * _mpts;
?????? int i;
?
_pts = calloc( n, sizeof( CvPoint2D64f ) );//特征點的坐標數組，分配空間
_mpts = calloc( n, sizeof( CvPoint2D64f ) );//對應匹配點的坐標數組，發配空間
?
//匹配類型是FEATURE_MDL_MATCH，匹配點的坐標是mdl_pt域
?????? if( mtype == FEATURE_MDL_MATCH )
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;//特征點的坐標
_mpts[i] = match->mdl_pt;//對應匹配點的坐標
????????????? }
//匹配類型不是FEATURE_MDL_MATCH，匹配點的坐標是img_pt域
?????? else
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;//特征點的坐標
_mpts[i] = match->img_pt;//對應匹配點的坐標
????????????? }
?
????????????? *pts = _pts;
????????????? *mpts = _mpts;
}
?
（2.2.5）xform_fn，一般傳入lsq_homog()函數代碼及說明：
答：
/* 根據4對坐標點計算最小二乘平面單應性變換矩陣
參數：
pts：坐標點數組
mpts：對應點數組，pts[i]與mpts[i]一一對應
n：pts和mpts數組中點的個數，pts和mpts中點的個數必須相同，一般是4
返回值：一個3*3的變換矩陣，將pts中的每一個點轉換為mpts中的對應點，返回值為空表示失敗
CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
{
?????? CvMat* H, * A, * B, X;
double x[9];//數組x中的元素就是變換矩陣H中的值
?????? int i;
?
//輸入點對個數不夠4
?????? if( n < 4 )
?????? {
????????????? fprintf( stderr, "Warning: too few points in lsq_homog(), %s line %d\n",
???????????????????? __FILE__, __LINE__ );
????????????? return NULL;
?????? }
?
//將變換矩陣H展開到一個8維列向量X中，使得AX=B，這樣只需一次解線性方程組即可求出X，然后再根據X恢復H
?????? /* set up matrices so we can unstack homography into X; AX = B */
A = cvCreateMat( 2*n, 8, CV_64FC1 );//創建2n*8的矩陣，一般是8*8
B = cvCreateMat( 2*n, 1, CV_64FC1 );//創建2n*1的矩陣，一般是8*1
X = cvMat( 8, 1, CV_64FC1, x );//創建8*1的矩陣，指定數據為x
H = cvCreateMat(3, 3, CV_64FC1);//創建3*3的矩陣
cvZero( A );//將A清零
?
//由于是展開計算，需要根據原來的矩陣計算法則重新分配矩陣A和B的值的排列
?????? for( i = 0; i < n; i++ )
?????? {
cvmSet( A, i, 0, pts[i].x );//設置矩陣A的i行0列的值為pts[i].x
????????????? cvmSet( A, i+n, 3, pts[i].x );
????????????? cvmSet( A, i, 1, pts[i].y );
????????????? cvmSet( A, i+n, 4, pts[i].y );
????????????? cvmSet( A, i, 2, 1.0 );
????????????? cvmSet( A, i+n, 5, 1.0 );
????????????? cvmSet( A, i, 6, -pts[i].x * mpts[i].x );
????????????? cvmSet( A, i, 7, -pts[i].y * mpts[i].x );
????????????? cvmSet( A, i+n, 6, -pts[i].x * mpts[i].y );
????????????? cvmSet( A, i+n, 7, -pts[i].y * mpts[i].y );
????????????? cvmSet( B, i, 0, mpts[i].x );
????????????? cvmSet( B, i+n, 0, mpts[i].y );
?????? }
?
//調用OpenCV函數，解線性方程組
cvSolve( A, B, &X, CV_SVD );//求X，使得AX=B
x[8] = 1.0;//變換矩陣的[3][3]位置的值為固定值1
?????? X = cvMat( 3, 3, CV_64FC1, x );
cvConvert( &X, H );//將數組轉換為矩陣
?
?????? cvReleaseMat( &A );
?????? cvReleaseMat( &B );
?????? return H;
}
?
問題1：cvMat函數？X = cvMat( 8, 1, CV_64FC1, x );//創建8*1的矩陣，指定數據為x？
答：
初始化矩陣:
double a[] = { 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12 };
CvMat Ma=cvMat(3, 4, CV_64FC1, a);
參看：1.opencv中矩陣計算的一些函數：https://www.cnblogs.com/qiaozhoulin/p/4556199.html
2.CvMat用法詳解：https://blog.csdn.net/zx3517288/article/details/51760541
?
問題2：四個點計算單應性變換矩陣？AX = B中A和B如何構建？
答：在齊次坐標中，假設一點p(xi,yi,1)經過H矩陣的變換變為p‘(xi',yi',1)，即 p' = H*p，通常，對于透視變換，H矩陣有8個自由度，這樣至少需要4對特征點對求解。4個特征點對可以建立8個方程。那么對于有n對特征點的情況(超定方程)，解p' = H*p方程組可以轉化為對齊次方程組Ax = 0 的求解。而對 Ax = 0 的求解轉化為 min ||Ax||2 的非線性優化問題（超定方程，通過最小二乘擬合得到近似解）。
對于某一點(xi,yi)，其變換可表述為 p' = H*p，代入展開可得：
(1)
那么可得：
(2)
進一步變換為：
(3)
這樣便可構造系數矩陣：
(4)
通過系數矩陣我們可以構造出齊次線性方程組(Ax = 0)：
(5)
即：
(6)
對于(6)這樣的超定方程求解，可以通過最小二乘的方式求解。通過對系數矩陣A求取特征值和特征向量得到。通過以下方式獲得最小二乘解：
[V,D] = eig(A'*A) (7)
其中D是特征值對角矩陣(特征值沿主對角線降序)，V是對應D特征值的特征向量(列向量)組成的特征矩陣,A'表示A的轉置。其最小二乘解為V(1)，即系數矩陣A最小特征值對應的特征向量就是超定方程組Ax = 0的最小二乘解。
至此，H矩陣已經求取，后續可以通過隨機采樣一致性(RANSC)進行精選，或者通過LM進行優化。
參看：原理參看：四點求解單應性矩陣：https://blog.csdn.net/hudaliquan/article/details/52121832
結果參看：特征點匹配——使用基礎矩陣、單應性矩陣的RANSAC算法去除誤匹配點對：https://blog.csdn.net/lhanchao/article/details/52849446
?
問題3：計算了單應性變換矩陣H后如何使用？
答：每次從所有的匹配點中選出4對，計算單應性矩陣H，然后選出內點個數最多的作為最終的結果。計算距離方法如下：
參看：特征點匹配——使用基礎矩陣、單應性矩陣的RANSAC算法去除誤匹配點對：https://blog.csdn.net/lhanchao/article/details/52849446
?
（2.2.6）find_consensus代碼及說明：
答：
/*對于給定的模型和錯誤度量函數，從特征點集和中找出一致集
參數：
features：特征點集合，其中的特征點都具有mtype類型的匹配點
n：特征點的個數
mtype：匹配類型，若是FEATURE_MDL_MATCH，對應的匹配點對坐標是每個特征點的img_pt域和其匹配點的mdl_pt域，否則，對應的匹配點對是每個特征點的img_pt域和其匹配點的img_pt域。
M：給定的模型，即一個變換矩陣
err_fn：錯誤度量函數，對于給定的變換矩陣，計算推定的匹配點對之間的變換誤差
err_tol：容錯度，用來衡量err_fn的返回值，小于err_tol的被加入一致集
consensus：輸出參數，一致集，即一致集中特征點構成的數組
返回值：一致集中特征點的個數
static int find_consensus( struct feature** features, int n, int mtype,
????????????????????????????????????????? ?? CvMat* M, ransac_err_fn err_fn, double err_tol,
????????????????????????????????????????? ?? struct feature*** consensus )
{
struct feature** _consensus;//一致集
struct feature* match;//每個特征點對應的匹配點
CvPoint2D64f pt, mpt;//pt:當前特征點的坐標，mpt：當前特征點的匹配點的坐標
double err;//變換誤差
?????? int i, in = 0;
?
_consensus = calloc( n, sizeof( struct feature* ) );//給一致集分配空間
?
//匹配類型是FEATURE_MDL_MATCH，匹配點的坐標是mdl_pt域
?????? if( mtype == FEATURE_MDL_MATCH )
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
pt = features[i]->img_pt;//特征點的坐標
mpt = match->mdl_pt;//對應匹配點的坐標
err = err_fn( pt, mpt, M );//計算"pt經M變換后的點"和mpt之間的距離的平方，即變換誤差
if( err <= err_tol )//若變換誤差小于容錯度，將其加入一致集
??????????????????????????? _consensus[in++] = features[i];
????????????? }
//匹配類型不是FEATURE_MDL_MATCH，匹配點的坐標是img_pt域
?????? else
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
pt = features[i]->img_pt;//特征點的坐標
mpt = match->img_pt;//對應匹配點的坐標
err = err_fn( pt, mpt, M );//計算"pt經M變換后的點"和mpt之間的距離的平方，即變換誤差
if( err <= err_tol )//若變換誤差小于容錯度，將其加入一致集
??????????????????????????? _consensus[in++] = features[i];
????????????? }
?????? *consensus = _consensus;
return in;//返回一致集中元素個數
}
?
（2.2.6.1）homog_xfer_err代碼及說明：
答：
/*對于給定的單應性矩陣H，計算輸入點pt精H變換后的點與其匹配點mpt之間的誤差
例如：給定點x，其對應點x'，單應性矩陣H，則計算x'與Hx之間的距離的平方，d(x', Hx)^2
參數：
pt：一個點
mpt：pt的對應匹配點
H：單應性矩陣
返回值：轉換誤差
*/
double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H )
{
CvPoint2D64f xpt = persp_xform_pt( pt, H );//pt經變換矩陣H變換后的點xpt,即H乘以x對應的向量
return sqrt( dist_sq_2D( xpt, mpt ) );//兩點間距離的平方
}
?
（2.2.6.1.1）persp_xform_pt代碼及說明：
答：
/*計算點pt經透視變換后的點，即給定一點pt和透視變換矩陣T，計算變換后的點
給定點(x,y)，變換矩陣M，計算[x',y',w']^T = M * [x,y,1]^T(^T表示轉置)，
則變換后的點是(u,v) = (x'/w', y'/w')
注意：仿射變換是透視變換的特例
參數：
pt：一個二維點
T：透視變換矩陣
返回值：pt經透視變換后的點
*/
CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T )
{
//XY：點pt對應的3*1列向量，UV：pt變換后的點對應的3*1列向量
?????? CvMat XY, UV;
double xy[3] = { pt.x, pt.y, 1.0 }, uv[3] = { 0 };//對應的數據
CvPoint2D64f rslt;//結果
?
//初始化矩陣頭
?????? cvInitMatHeader( &XY, 3, 1, CV_64FC1, xy, CV_AUTOSTEP );
?????? cvInitMatHeader( &UV, 3, 1, CV_64FC1, uv, CV_AUTOSTEP );
cvMatMul( T, &XY, &UV );//計算矩陣乘法，T*XY,結果放在UV中
rslt = cvPoint2D64f( uv[0] / uv[2], uv[1] / uv[2] );//得到轉換后的點
?
?????? return rslt;
}
?
問題1：為什么要進行透視變換？公式是怎么樣的？
答：透視變換能保持“直線性”，即原圖像里面的直線，經透視變換后仍為直線。
透視變換矩陣變換公式為：
其中透視變換矩陣：
要移動的點，即源目標點為：
另外定點，即移動到的目標點為：
這是一個從二維空間變換到三維空間的轉換，因為圖像在二維平面，故除以Z，? （X';Y';Z'）表示圖像上的點：
令，展開上面公式，得到一個點的情況：
?
4個點可以得到8個方程，即可解出A。
參看：圖像透視變換原理及實現——https://blog.csdn.net/cuixing001/article/details/80261189
?
問題2：cvMatMul函數怎么使用？
答：#define cvMatMul(src1,src2,dst); 表示dst=src1*src2
參看：cvGEMM、cvMatMul和cvMatMulAdd的定義——https://blog.csdn.net/liulianfanjianshi/article/details/11737921
?
（2.2.6.1.2）dist_sq_2D代碼及說明：
答：
//計算兩點之間距離的平方
double dist_sq_2D( CvPoint2D64f p1, CvPoint2D64f p2 )
{
? double x_diff = p1.x - p2.x;
? double y_diff = p1.y - p2.y;
?
? return x_diff * x_diff + y_diff * y_diff;? //(x1-x2)^2+(y1-y2)^2
}
?
問題：如何計算兩點間的距離——歐式距離？
參看：百度百科（歐幾里得度量）——https://baike.baidu.com/item/%E6%AC%A7%E5%87%A0%E9%87%8C%E5%BE%97%E5%BA%A6%E9%87%8F/1274107?fr=aladdin
?
?
（2.2.6）release_mem代碼及說明：
答：
/*釋放輸入參數的存儲空間
static __inline void release_mem( CvPoint2D64f* pts1, CvPoint2D64f* pts2,
??????????????????????????????????????????????????????? ? struct feature** features )
{
?????? free( pts1 );
?????? free( pts2 );
?????? if( features )
????????????? free( features );
}
?
（2.2.7）extract_corresp_pts代碼及說明：
答：
/*從特征點數組中獲取特征點和其對應匹配點的二維坐標，分別放到輸出參數pts和mpts中
參數：
features：特征點數組，將從其中抽取坐標點和其匹配點，此數組中所有特征點都應具有mtype類型的匹配點
n：feantures中特征點個數，n=4
mtype：匹配類型，若是FEATURE_MDL_MATCH，對應的匹配點對坐標是每個特征點的img_pt域和其匹配點的mdl_pt域，否則，對應的匹配點對是每個特征點的img_pt域和其匹配點的img_pt域。三者之一
pts：輸出參數，從特征點數組features中獲取的二維坐標數組
mpts：輸出參數，從特征點數組features的對應匹配點中獲取的二維坐標數組
static void extract_corresp_pts( struct feature** features, int n, int mtype,
??????????????????????????????????????????????????????? ?CvPoint2D64f** pts, CvPoint2D64f** mpts )
{
struct feature* match;//每個特征點對應的匹配點
?????? CvPoint2D64f* _pts, * _mpts;
?????? int i;
?
_pts = calloc( n, sizeof( CvPoint2D64f ) );//特征點的坐標數組
_mpts = calloc( n, sizeof( CvPoint2D64f ) );//對應匹配點的坐標數組
?
//匹配類型是FEATURE_MDL_MATCH，匹配點的坐標是mdl_pt域
?????? if( mtype == FEATURE_MDL_MATCH )
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;//特征點的坐標
_mpts[i] = match->mdl_pt;//對應匹配點的坐標
????????????? }
//匹配類型不是FEATURE_MDL_MATCH，匹配點的坐標是img_pt域
?????? else
????????????? for( i = 0; i < n; i++ )
????????????? {
//根據給定的匹配類型，返回輸入特征點的匹配點
???????????????????? match = get_match( features[i], mtype );
???????????????????? if( ! match )
??????????????????????????? fatal_error( "feature does not have match of type %d, %s line %d",
???????????????????????????????????????????????? mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;//特征點的坐標
_mpts[i] = match->img_pt;//對應匹配點的坐標
????????????? }
?
????????????? *pts = _pts;
????????????? *mpts = _mpts;
}
?
（2.3）cvWarpPerspective函數說明：
答：
說明：密集透視變換
參看：1）第六章 - 圖像變換 - 圖像拉伸、收縮、扭曲、旋轉[2] - 透視變換(cvWarpPerspective)——https://blog.csdn.net/hitwengqi/article/details/6890220
2）【OpenCV3】透視變換——cv::getPerspectiveTransform()與cv::warpPerspective()詳解——https://blog.csdn.net/guduruyu/article/details/72518340
?
?
?附：
1.CvPoint函數說明？
參看：opencv的基本數據類型CvPoint，CvSize，CvRect和CvScalar——https://blog.csdn.net/gdut2015go/article/details/46301821
?
?
轉載于:https://www.cnblogs.com/Alliswell-WP/p/SIFT_code5.html
與50位技術專家面對面20年技術見證，附贈技術全景圖總結
以上是生活随笔為你收集整理的RobHess的SIFT代码解析之RANSAC的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： python守护线程t.setDaemo
- 下一篇： git的忽略文件语法规范