SIFT演算法特徵描述子構建---關鍵點定位原理及程式碼
阿新 • • 發佈:2019-01-01
0.引言
sift針對區域性特徵進行特徵提取,在尺度空間尋找極值點,提取位置,尺度,旋轉不變數,生成特徵描述子。
總共分四個步驟:
step2 關鍵點/極值點提取
2.1 關鍵點位置初步探查
生成DOG金字塔後,要找到DOG空間中的區域性極值點。第一步首先根據DOG空間內,各點響應大小檢測和篩選極值點。
選擇策略如下:
中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點
比較,以確保在尺度空間和二維影象空間都檢測到極值點。
由於要在相鄰尺度進行比較,高斯差分金子塔S+2層,頭尾兩層無法計算,因此,定義的INTERVALS層數就是指檢測極值點的層數。
bool isExtremum(int x,int y,const vector<Mat>& dog_pyr,int index)//x為列,y為行
{
pixel_t* data = (pixel_t*)dog_pyr[index].data; //當前層資料的指標
pixel_t val = *(data + dog_pyr[index].cols*y + x); //當前點畫素
//檢查是否為極大值
if (val>0)
{
for (int i = -1; i < 2; i++)
{
int col_cur= dog_pyr[index+i].cols; //當前層列數
for (int j = -1; j < 2; j++)
{
for (int k = -1; k < 2; k++)
{
if (val<*((pixel_t*)dog_pyr[index+i].data+col_cur*(y+j)+x+k))
{
return false;
}
}
}
}
}//極大值
//檢查是否為極小值
else
{
for (int i = -1; i < 2; i++)
{
int col_cur = dog_pyr[index + i].cols; //當前層列數
for (int j = -1; j < 2; j++)
{
for (int k = -1; k < 2; k++)
{
if (val > *((pixel_t*)dog_pyr[index + i].data + col_cur*(y + j) + x + k))
{
return false;
}
}
}
}
}//極小值
return true;
}
當然這樣產生離散空間的極值點並不全都是穩定的特徵點,因為某些極值點響應較弱(對比度低)
,同時DOG運算元會產生較強的邊緣響應
。為了解決這兩個問題,需要對極值點進行插值和篩選。
2.2 關鍵點插值定位修正和邊緣響應去除
2.2.1位置修正及去除低響應
插值需要將離散的影象表示式,表達為連續的函式。以下通過擬合三維二次函式(二階的泰勒展開式)來精確確定關鍵點的位置和尺度。
其中是位置和尺度三個變數的向量。
二次函式對稱軸就是極值點偏移量:
此處,為了求此公式,需要提前定義好計算逆矩陣
,二階偏導
,一階偏導
的函式。
/*************************************************************************************************************************
*模組說明:
* 求偏差量預備工作 【1】返回第index層的金字塔在y行x列的值
*功能說明:
* PryAt()---At:在該層y行x列的值
**************************************************************************************************************************/
pixel_t PyrAt(const vector<Mat>& pyr, int index, int x, int y)
{
pixel_t *data = (pixel_t*)pyr[index].data;
pixel_t val = *(data + y*pyr[index].cols + x);
return val;
}
#define At(index, x, y) (PyrAt(dog_pyr, (index), (x), (y)))
/*************************************************************************************************************************
*模組說明:
* 求偏差量預備工作,【2】對x,y,sigma的一階偏導存入dx[3]
*功能說明:
* 使用節點的一階中心差分近似導數, 取最小h=1時,最精確
* f'(xi)=(f(xi+h)-f(xi-h))/2h
* f'(yi)=(f(yi+h)-f(yi-h))/2h
* f'(si)=(f(si+h)-f(si-h))/2h
**************************************************************************************************************************/
void DerivativeOf3D(int x,int y,const vector<Mat>& dog_pyr,int index,double *dx)
{
dx[0] = (double)(At(index, x + 1, y) - At(index, x - 1, y)) / 2.0;
dx[1] = (double)(At(index, x, y + 1) - At(index, x, y - 1)) / 2.0;
dx[2] = (double)(At(index + 1, x, y) - At(index - 1, x, y)) / 2.0;
}
/*************************************************************************************************************************
*模組說明:
* 求偏差量預備工作,【3】對x,y,sigma的二階偏導存入Hat(3*3),對稱陣
* | Dxx Dxy Dxs |
* | Dyx Dxy Dys |
* | Dsx Dsy Dss |
*功能說明:
* 使用節點的一階中心差分近似導數, 取最小h=1時,最精確
* Dxx=f''(xi)=(f(xi+h)+f(xi-h)-2*f(xi))/h^2
* Dxy=f''(yi)=(f(yi+h)+f(yi-h)-2*f(yi))/h^2
* Dss=f''(si)=(f(si+h)+f(si-h)-2*f(si)))/h^2
*
* Dyx=Dxy=f''(xi)(yi)=(f(xi+h,yi+h)+f(xi-h,yi-h)-f(xi+h,yi-h)-f(xi-h,yi+h))/4*h
**************************************************************************************************************************/
#define Hat(i,j) (*(H+3*(i)+(j)))
void Hessian3D(int x, int y, const vector<Mat>& dog_pyr, int index, double*H)
{
Hat(0, 0) = At(index, x + 1, y) + At(index, x - 1, y) - 2 * At(index, x, y);//Dxx
Hat(1, 1) = At(index, x, y + 1) + At(index, x, y - 1) - 2 * At(index, x, y);//Dyy
Hat(1, 1) = At(index + 1, x, y) + At(index - 1, x, y) - 2 * At(index, x, y);//Dss
Hat(1, 0) = Hat(0, 1) = (At(index, x + 1, y + 1) + At(index, x - 1, y - 1) - At(index, x + 1, y - 1) - At(index, x - 1, y + 1)) / 4.0;//Dxy
Hat(2, 0) = Hat(0, 2) = (At(index + 1, x + 1, y) + At(index - 1, x - 1, y) - At(index - 1, x + 1, y) - At(index + 1, x - 1, y)) / 4.0;//Dxs
Hat(2, 1) = Hat(1, 2) = (At(index + 1, x, y + 1) + At(index - 1, x, y - 1) - At(index + 1, x, y - 1) - At(index - 1, x, y + 1)) / 4.0;//Dsy
}
/*************************************************************************************************************************
*模組說明:
* 三階矩陣求逆
**************************************************************************************************************************/
#define HIat(i, j) (*(H_inve+3*(i)+j))
bool Inverse3D(const double *H,double *H_inve)
{
//行列式
double A =Hat(0, 0)*Hat(1, 1)*Hat(2, 2)
+ Hat(0, 1)*Hat(1, 2)*Hat(2, 0)
+ Hat(0, 2)*Hat(1, 0)*Hat(2, 1)
- Hat(0, 0)*Hat(1, 2)*Hat(2, 1)
- Hat(0, 1)*Hat(1, 0)*Hat(2, 2)
- Hat(0, 2)*Hat(1, 1)*Hat(2, 0);
if (fabs(A)<1e-10) return false;
HIat(0, 0) = Hat(1, 1) * Hat(2, 2) - Hat(2, 1)*Hat(1, 2);
HIat(0, 1) = -(Hat(0, 1) * Hat(2, 2) - Hat(2, 1) * Hat(0, 2));
HIat(0, 2) = Hat(0, 1) * Hat(1, 2) - Hat(0, 2)*Hat(1, 1);
//伴隨矩陣第二行
HIat(1, 0) = Hat(1, 2) * Hat(2, 0) - Hat(2, 2)*Hat(1, 0);
HIat(1, 1) = -(Hat(0, 2) * Hat(2, 0) - Hat(0, 0) * Hat(2, 2));
HIat(1, 2) = Hat(0, 2) * Hat(1, 0) - Hat(0, 0)*Hat(1, 2);
//伴隨矩陣第三行
HIat(2, 0) = Hat(1, 0) * Hat(2, 1) - Hat(1, 1)*Hat(2, 0);
HIat(2, 1) = -(Hat(0, 0) * Hat(2, 1) - Hat(0, 1) * Hat(2, 0));
HIat(2, 2) = Hat(0, 0) * Hat(1, 1) - Hat(0, 1)*Hat(1, 0);
for (int i = 0; i < 9; i++)
{
*(H_inve + i) /= A;
}
return true;
}
這樣就可以計算偏差量了。
/*************************************************************************************************************************
*模組說明:計算x,y,sigma 的偏差量x^、y^、sigma^
*
*模組內容 1.x^=-(Dxx*dx+Dxy*dy+Dxs*ds) DOG函式擬合函式的逆矩陣的二階和擬合函式一階(為什麼這個公式)
**************************************************************************************************************************/
void GetoffsetX(int x, int y, const vector<Mat>& dog_pyr, int index, double* offset_x)
{
//求一階偏導
double dx[3];
DerivativeOf3D(x, y, dog_pyr, index, dx);
//求二階偏導
double H[9], H_inve[9] = {0};
Hessian3D(x, y, dog_pyr, index, H);
if (Inverse3D(H, H_inve))
{
for (int i = 0; i < 3; i++)
{
offset_x[i] = 0.0;
for (int j = 0; j < 3; j++)
{
offset_x[i] += H_inve[i * 3 + j] * dx[j];
}
offset_x[i]=-offset_x[i];
}
}
else
cout << "無海森逆矩陣" << endl;
得到偏移量後,首先判斷x,y,sigma任意維度的偏移量是否超過0.5,其次論文中擬合的迭代次數為5次。另外,對比度小於0.03的 極值點也要去除。
極值點處的值計算公式為:
//2.D(x^)=D+0.5*D'^T(x)*x^=D()+0.5*(Dx*offset[0] +Dy* offset[1]+Ds* offset[2])
//計算D(x^)
double GetFabsDx(int x,int y,const vector<Mat>& dog_pyr,int index,const double *offset_x)
{
double dx[3];
DerivativeOf3D(x, y, dog_pyr, index, dx);
double term = 0.0;
for (int i = 0; i < 3; i++)
{
term += dx[i] * offset_x[i];
}
pixel_t val = At(index, x, y);
return fabs(val+0.5*term);
}
整個篩選並儲存關鍵點如下:
/*************************************************************************************************************************
*模組說明:
* 模組四的第二步:修正極值點,刪除不穩定的點
*功能說明:
* 1--根據高斯差分函式產生的極值點並不全都是穩定的特徵點,因為某些極值點的響應較弱,而且DOG運算元會產生較強的邊緣響應
* 2--以上方法檢測到的極值點是離散空間的極值點,下面通過擬合三維二次函式來精確定位關鍵點的位置和尺度,同時去除對比度
* 低和不穩定的邊緣響應點(因為DOG運算元會產生較強的邊緣響應),以增強匹配的穩定性、提高抗噪聲的能力。
* 3--修正極值點,刪除不穩定點,|D(x)| < 0.03 Lowe 2004
**************************************************************************************************************************/
Keypoint* InterpolationExtremum(int x ,int y, const vector<Mat>& dog_pyr, int index, int octave, int interval, double dxthreshold)
{
double offset_x[3] = {0};
const Mat &mat = dog_pyr[index];
int idx = index;
//int intvl=interval;
int i = 0;
for (; i < MAX_INTERPOLATION_STEPS; i++)
{
//求取偏移量
GetoffsetX(x, y, dog_pyr, idx, offset_x);
//如果offset_x 的任一維度大於0.5,插值中心已經偏移到它的鄰近點上
if (fabs(offset_x[0]) < 0.5 && fabs(offset_x[1]) < 0.5 && fabs(offset_x[2]) < 0.5) break;
////將找到的極值點對應成畫素(整數),用周圍點代替
x += cvRound(offset_x[0]);
y += cvRound(offset_x[1]);
interval += cvRound(offset_x[2]);
idx = octave* (INTERVALS+2)+ interval;
//超出影象層或者影象邊界的範圍
if (interval<1||interval>INTERVALS||x>mat.cols-2||x<2|| y>mat.rows - 2 || y<2) return NULL;
}
//超出所設定的迭代次數
if (i>=MAX_INTERPOLATION_STEPS) return NULL;
//|D(x^)| < 0.03取經驗值
if (GetFabsDx(x, y, dog_pyr, idx, offset_x) < dxthreshold / INTERVALS) return NULL;
//則屬於關鍵點
Keypoint* keypoint = new Keypoint;
keypoint->x = x;
keypoint->y = y;
//keypoint->val=At()
keypoint->offset_x = offset_x[0];
keypoint->offset_y = offset_x[1];
keypoint->offset_interval = offset_x[2];
keypoint->octave = octave;
keypoint->interval = interval;
keypoint->dx = (x + offset_x[0])*pow(2.0, octave); //縮放成原影象大小
keypoint->dy = (y + offset_x[1])*pow(2.0, octave);
return keypoint;
}
2.2.2邊緣響應去除
為了剔除邊緣響應點,需要讓下比值小於一定的閾值
)成立時將關鍵點保留,反之剔除。
其中,在Lowe的文章中,取r=10。在Lowe的文章中,取r=10。
bool passEdgeResponse(int x,int y,vector<Mat>& dog_pyr,int index,double r)
{
//求二階偏導
double H[9];
Hessian3D(x, y, dog_pyr, index, H);
double Tr_h = H[0] + H[4];//Dxx+Dyy
double Det_h = H[0] * H[4] - H[1] * H[1];//Dxx * Dyy - Dxy * Dxy
if (Det_h<=0) return false;
if (Tr_h*Tr_h / Det_h < (r + 1)*(r + 1) / r) return true;
return false;
}
2.2.3 關鍵點生成儲存
通過呼叫上述過程函式,將關鍵點位置尺度資訊儲存。
關鍵點資料結構:
struct Keypoint
{
int octave; //【1】關鍵點所在組
int interval; //【2】關鍵點所在層
double offset_interval; //【3】調整後層的增量(是3/2嗎,怎麼是double?)
int x, y; //【4】x,y座標,根據octave和interval可取的層內影象(是指所在組層?)
double scale; //【5】空間尺度座標scale = sigma0*pow(2.0, o+s/S)
double dx, dy; //【6】特徵點座標,該座標被縮放成原影象大小
double offset_x, offset_y;
//============================================================
//特徵描述子
//============================================================
double octave_scale; //【1】offset_i
double ori; //【2】方向
double descriptor[FEATURE_ELEMENT_LENTH];//【3】特徵描述子
int descr_length;
double val; //【4】極值
};
儲存函式:
void DetectionLocalExtrema(vector<Mat>& dog_pyr,vector<Keypoint>& extrema,int intervals)
{
int octaves = (int)dog_pyr.size() / (intervals+2); //DOG組數
double thresh=0.5*DXTHRESHHOLD/intervals;
for (int o = 0; o < octaves; o++)
{
for (int i = 1; i < intervals+1; i++)
{
int index = o*(intervals + 2) + i; //當前層的索引序號
pixel_t *data =(pixel_t*) dog_pyr[index].data; //當前層資料首地址
//對每層進行提取特徵點
for (int y = IMG_BORDER; y < dog_pyr[index].rows-IMG_BORDER; y++)
{
for (int x = IMG_BORDER; x < dog_pyr[index].cols-IMG_BORDER; x++)
{
pixel_t val = *(data + y*dog_pyr[index].cols + x);
//當大於閾值,並且時是26箇中極值點時
if (fabs(val)>thresh&&isExtremum(x,y,dog_pyr,index))
{
Keypoint *extrum = InterpolationExtremum(x, y, dog_pyr, index, o, i);
if (extrum)
{
if (passEdgeResponse(extrum->x,extrum->y,dog_pyr,index))
{
extrum->val = At(index,extrum->x,extrum->y);
extrema.push_back(*extrum);
}
delete extrum;
}
}
}
}
}
}
}
2.2.4 尺度
void CalculateScale(vector<Keypoint>& features, double sigma, int intervals)
{
double intvl = 0.0;
for (int i = 0; i < features.size(); i++)
{
intvl = features[i].interval + features[i].offset_interval;
features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
}
}
//對擴大的影象特徵縮放,對應到原圖大小(非-1層)
void HalfFeatures(vector<Keypoint>& features)
{
for (int i = 0; i < features.size(); i++)
{
features[i].dx /= 2;
features[i].dy /= 2;
features[i].scale /= 2;
}
}