小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看

搜索
分享

blue_lg的surf c++源碼 解析

 學(xué)海無涯GL 2012-09-11

blue_lg的surf c++源碼 解析

說是surf，其實原算法是基于Notes on the OpenSURF Library這篇文檔。

下載地址：http://opensurf1./files/OpenSURFcpp.zip

首先定義 #define PROCEDURE 1

緊接著進入main()主函數(shù)：

int main(void)
{
if (PROCEDURE == 1) return mainImage();//靜態(tài)圖像的特征點監(jiān)測
if (PROCEDURE == 2) return mainVideo();//通過攝像頭實時監(jiān)測，提取特征點
if (PROCEDURE == 3) return mainMatch();//圖像匹配算法，在視頻序列里尋找一個已知物體圖像
if (PROCEDURE == 4) return mainMotionPoints();//視頻中行為監(jiān)測
if (PROCEDURE == 5) return mainStaticMatch();//靜態(tài)圖像間的特征點匹配
if (PROCEDURE == 6) return mainKmeans();//靜態(tài)圖像的k-means聚類
}

我們以監(jiān)測靜態(tài)圖像特征點（即1）為例了解surf算法如何提取特征點。

int mainImage(void)
{
// Declare Ipoints and other stuff
IpVec ipts;
IplImage *img=cvLoadImage("imgs/sf.jpg");
// Detect and describe interest points in the image
clock_t start = clock();
surfDetDes(img, ipts, false, 5, 4, 2, 0.0004f);//URF描述子特征提取實現(xiàn)函數(shù)
clock_t end = clock();
std::cout<< "OpenSURF found: " << ipts.size() << " interest points" << std::endl;
std::cout<< "OpenSURF took: " << float(end - start) / CLOCKS_PER_SEC << " seconds" << std::endl;
// Draw the detected points
drawIpoints(img, ipts);
// Display the result
showImage(img);
return 0;
}

EXP:

IpVec的定義在ipoint.h里，typedef std::vector<Ipoint> IpVec; 也就是說，IpVec是一個vector向量，每個成員是一個Ipoint。而Ipoint是一個類，也在ipoint.h里，作者將它按照結(jié)構(gòu)體使用，都是public。

clock()的使用是為了測試程序運行的時間，一個記錄起始時間，一個記錄結(jié)束時間，最終將兩者只差輸出，即surfDetDes()特征提取所需要的時間。

drawIpoints()函數(shù)是在img基礎(chǔ)上增加特征點的描述，說白了，就是添加特征點的函數(shù)。

那么，現(xiàn)在進入到surfDetDes()函數(shù)內(nèi)部來探個究竟吧。

surfDetDes定義在surflib.h里面：

//! Library function builds vector of described interest points
inline void surfDetDes(IplImage *img, /* image to find Ipoints in */
std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
bool upright = false, /* run in rotation invariant mode? */
int octaves = OCTAVES, /* number of octaves to calculate */
int intervals = INTERVALS, /* number of intervals per octave */
int init_sample = INIT_SAMPLE, /* initial sampling step */
float thres = THRES /* blob response threshold */)
{
// Create integral-image representation of the image
IplImage *int_img = Integral(img);
// Create Fast Hessian Object
FastHessian fh(int_img, ipts, octaves, intervals, init_sample, thres);
// Extract interest points and store in vector ipts
fh.getIpoints();
// Create Surf Descriptor Object
Surf des(int_img, ipts);
// Extract the descriptors for the ipts
des.getDescriptors(upright);
// Deallocate the integral image
cvReleaseImage(&int_img);
}

　　總的來說，surfDetDes()里面規(guī)劃了具體的步驟：

1.獲取積分圖像init_img;

2.計算Hessian矩陣特征fh;

3.利用fh提取特征點;

4.創(chuàng)建surf特征描述符，并和特征點貫聯(lián);

5.釋放積分圖像。

①積分圖像的實現(xiàn)

首先在Integral.cpp里面找到Integral()，如下：

IplImage *Integral(IplImage *source)
{
// convert the image to single channel 32f
IplImage *img = getGray(source);
IplImage *int_img = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);
// set up variables for data access
int height = img->height;
int width = img->width;
int step = img->widthStep/sizeof(float);
float *data = (float *) img->imageData;
float *i_data = (float *) int_img->imageData;
// first row only
float rs = 0.0f;
for(int j=0; j<width; j++)
{
rs += data[j];
i_data[j] = rs;
}
// remaining cells are sum above and to the left
for(int i=1; i<height; ++i)
{
rs = 0.0f;
for(int j=0; j<width; ++j)
{
rs += data[i*step+j];
i_data[i*step+j] = rs + i_data[(i-1)*step+j];
}
}
// release the gray image
cvReleaseImage(&img);
// return the integral image
return int_img;
}

1.　　首先將原輸入轉(zhuǎn)化為灰度圖像，并創(chuàng)建一個大小等于灰度圖像gray-image的圖像數(shù)組--積分圖像int_img。

2.　　獲取圖像的信息，比如大?。ǜ遠(yuǎn)eight和寬width）以及gray-image和積分圖像int_img的數(shù)據(jù)首地址data && i_data。(注意此時數(shù)據(jù)類型為float)

3.　　首先計算第一行像素的積分值，相當(dāng)于一維數(shù)據(jù)的累加。其他數(shù)據(jù)通過迭代計算獲取，i_data[i*step+j] = rs + i_data[(i-1)*step+j]，若當(dāng)前點為（i0,j0）,其中rs就為第（i+1）行（即i=i0）行j<=j0之前所有像素值和。 如下所示：

[其中黑色為當(dāng)前點i_data[i*step+j]，綠色為i_data[(i-1)*step+j]，而rs=橫放著的和黑色點同行的那塊矩形框?qū)?yīng)的區(qū)域像素值之和]

4.　　釋放灰度圖像，并返回積分圖像。

相關(guān)函數(shù)integral.h中的BoxIntegral().

inline float BoxIntegral(IplImage *img, int row, int col, int rows, int cols)

其中，幾個參數(shù)意思分別為源圖像，row,col為A點的坐標(biāo)值，rows和cols分別為高和寬。

利用上面的積分圖像計算 A B 這樣的box區(qū)域里面所有像素點的灰度值之和。S=int_img(D)+int_img(A)-int_img(B)-int_img(C).

C D

②Hessian矩陣特征的計算

FastHessian，計算hessian矩陣的類，它的定義在fasthessian.h里，實現(xiàn)在fasthessian.cpp里。

class FastHessian {
public:
//! Constructor without image
FastHessian(std::vector<Ipoint> &ipts,
const int octaves = OCTAVES,
const int intervals = INTERVALS,
const int init_sample = INIT_SAMPLE,
const float thres = THRES);
//! Constructor with image
FastHessian(IplImage *img,
std::vector<Ipoint> &ipts,
const int octaves = OCTAVES,
const int intervals = INTERVALS,
const int init_sample = INIT_SAMPLE,
const float thres = THRES);
//! Destructor
~FastHessian();
//! Save the parameters
void saveParameters(const int octaves,
const int intervals,
const int init_sample,
const float thres);
//! Set or re-set the integral image source
void setIntImage(IplImage *img);
//! Find the image features and write into vector of features
void getIpoints();
private:
//---------------- Private Functions -----------------//
//! Build map of DoH responses
void buildResponseMap();
//! Calculate DoH responses for supplied layer
void buildResponseLayer(ResponseLayer *r);
//! 3x3x3 Extrema test
int isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);
//! Interpolation functions - adapted from Lowe's SIFT implementation
void interpolateExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);
void interpolateStep(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b,
double* xi, double* xr, double* xc );
CvMat* deriv3D(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);
CvMat* hessian3D(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);
//---------------- Private Variables -----------------//
//! Pointer to the integral Image, and its attributes
IplImage *img;
int i_width, i_height;
//! Reference to vector of features passed from outside
std::vector<Ipoint> &ipts;
//! Response stack of determinant of hessian values
std::vector<ResponseLayer *> responseMap;
//! Number of Octaves
int octaves;
//! Number of Intervals per octave
int intervals;
//! Initial sampling step for Ipoint detection
int init_sample;
//! Threshold value for blob resonses
float thresh;
};

　　在public里面定義了兩種構(gòu)造函數(shù)分別對應(yīng)有無源圖像這一項參數(shù)，緊接著還定義了析構(gòu)函數(shù)~FastHessian等函數(shù)。下面在fasthessian.cpp對這些函數(shù)的實現(xiàn)一一解釋：

兩個構(gòu)造函數(shù)都是調(diào)用了saveParameters(octaves, intervals, init_sample, thresh)設(shè)置構(gòu)造金字塔的參數(shù)，而帶圖像的構(gòu)造函數(shù)另外多加了一句setIntImage(img)用來設(shè)置當(dāng)前圖像。

//! Save the parameters
void FastHessian::saveParameters(const int octaves, const int intervals,
const int init_sample, const float thresh)
{
// Initialise variables with bounds-checked values
this->octaves =
(octaves > 0 && octaves <= 4 ? octaves : OCTAVES);
this->intervals =
(intervals > 0 && intervals <= 4 ? intervals : INTERVALS);
this->init_sample =
(init_sample > 0 && init_sample <= 6 ? init_sample : INIT_SAMPLE);
this->thresh = (thresh >= 0 ? thresh : THRES);
}
//! Set or re-set the integral image source
void FastHessian::setIntImage(IplImage *img)
{
// Change the source image
this->img = img;
i_height = img->height;
i_width = img->width;
}

　　由于在h頭文件中已設(shè)置

static const int OCTAVES = 5;//組數(shù)
static const int INTERVALS = 4;//每組層數(shù)
static const float THRES = 0.0004f;//閾值
static const int INIT_SAMPLE = 2;//初始采樣因子

　 所以　saveParameters的作用就是調(diào)整參數(shù)，以防過大或過小。

FastHessian::getIpoints()提取興趣點：

//! Find the image features and write into vector of features
void FastHessian::getIpoints()
{
// filter index map
static const int filter_map [OCTAVES][INTERVALS] = {{0,1,2,3}, {1,3,4,5}, {3,5,6,7}, {5,7,8,9}, {7,9,10,11}};
// Clear the vector of exisiting ipts
ipts.clear();
// Build the response map
buildResponseMap();
// Get the response layers
...<BR> ...
}

　　首先初始化filter_map，清空標(biāo)記特征點的ipts結(jié)構(gòu)體。

創(chuàng)建高斯平滑層函數(shù)參數(shù)ResponseMap()，大小與論文所給完全一致，

　　　　　　　　// Oct1: 9, 15, 21, 27
　　　　　　　　 // Oct2: 15, 27, 39, 51
　　　　　　　　 // Oct3: 27, 51, 75, 99
　　　　　　　　 // Oct4: 51, 99, 147,195
　　　　　　　　 // Oct5: 99, 195,291,387

這些都是每組模板的大小，每組間隔遞增，6,12,24,48,96 。ResponseMap這個結(jié)構(gòu)體向量包含4個參數(shù)ResponseLayer(int width, int height, int step, int filter)定義在responsemap.h里面，其中width和height等于實際圖像大小除以step(step初始值為2)，而filter則是濾波器半徑。

然后使用buildResponseLayer(responseMap[i])對圖像處理后將數(shù)據(jù)存放在responses和laplacian兩個數(shù)組里面。

void FastHessian::buildResponseLayer(ResponseLayer *rl)
{
float *responses = rl->responses; // response storage
unsigned char *laplacian = rl->laplacian; // laplacian sign storage
int step = rl->step; // step size for this filter 濾波器尺度因子
int b = (rl->filter - 1) / 2; // border for this filter 濾波器邊界
int l = rl->filter / 3; // lobe for this filter (filter size / 3)
int w = rl->filter; // filter size 濾波器大小
float inverse_area = 1.f/(w*w); // normalisation factor 標(biāo)準(zhǔn)化因子
float Dxx, Dyy, Dxy;
for(int r, c, ar = 0, index = 0; ar < rl->height; ++ar)
{
for(int ac = 0; ac < rl->width; ++ac, index++)
{
// get the image coordinates
r = ar * step;
c = ac * step;
// Compute response components
Dxx = BoxIntegral(img, r - l + 1, c - b, 2*l - 1, w)
- BoxIntegral(img, r - l + 1, c - l / 2, 2*l - 1, l)*3;
Dyy = BoxIntegral(img, r - b, c - l + 1, w, 2*l - 1)
- BoxIntegral(img, r - l / 2, c - l + 1, l, 2*l - 1)*3;
Dxy = + BoxIntegral(img, r - l, c + 1, l, l)
+ BoxIntegral(img, r + 1, c - l, l, l)
- BoxIntegral(img, r - l, c - l, l, l)
- BoxIntegral(img, r + 1, c + 1, l, l);
// Normalise the filter responses with respect to their size
Dxx *= inverse_area;
Dyy *= inverse_area;
Dxy *= inverse_area;
// Get the determinant of hessian response & laplacian sign
responses[index] = (Dxx * Dyy - 0.81f * Dxy * Dxy);
laplacian[index] = (Dxx + Dyy >= 0 ? 1 : 0);
#ifdef RL_DEBUG
// create list of the image coords for each response
rl->coords.push_back(std::make_pair<int,int>(r,c));
#endif
}
}
}

　　其中計算Dxy和Dyy的示意圖如下，其他的Dxx(Dyy的轉(zhuǎn)置)讀者自己參考?！敬藭rw=9,l=9/3=3,對應(yīng)于右圖的總計算區(qū)域高度和寬度2*l-1】

圓點為當(dāng)前點，將紅色方形區(qū)域1內(nèi)的積分值減去綠色方形2區(qū)域內(nèi)的積分值=Dxy*w2 綠色方形區(qū)域2內(nèi)的積分值減去2*紅色色方形區(qū)域1內(nèi)的積分值=Dyy*w2 (==用一整塊區(qū)域-3*紅色區(qū)域）

最后將計算后的結(jié)果存放到ResponseLayer里面的response和laplacian一維數(shù)組里，數(shù)組的大小即為ResponseLayer->width * ResponseLayer->width。

這樣就計算出了每一層的所有像素點的det(Happrox)=Dxx*Dyy-(0.9*Dxy)2,下面開始判斷當(dāng)前點是否是極值點。

③根據(jù)上一步計算每層的每個點的det(Happrox)，判斷極值點。

在fasthessian.cpp里面找到getIpoints()，下面開始抽取每組（共octaves組）相鄰的3層（共有4層，每組有2個這樣層的滿足）：

ResponseLayer *b, *m, *t;
for (int o = 0; o < octaves; ++o) for (int i = 0; i <= 1; ++i)
{
b = responseMap.at(filter_map[o][i]);
m = responseMap.at(filter_map[o][i+1]);
t = responseMap.at(filter_map[o][i+2]);
// loop over middle response layer at density of the most
// sparse layer (always top), to find maxima across scale and space<BR> //總是在最上層去尋找極大值點
for (int r = 0; r < t->height; ++r)
{
for (int c = 0; c < t->width; ++c)
{
if (isExtremum(r, c, t, m, b))
{
interpolateExtremum(r, c, t, m, b);
}
}
}
}

　　在判斷的時候用到了isExtremum()和interpolateExtremum()子函數(shù)，在當(dāng)前cpp內(nèi)找到它們，并分析。

//! Non Maximal Suppression function
int FastHessian::isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)
{
// bounds check
int layerBorder = (t->filter + 1) / (2 * t->step);//邊界值，若當(dāng)前點與邊界的距離小于此值，則無法確定鄰域，顧認(rèn)為當(dāng)前點不是極值點
if (r <= layerBorder || r >= t->height - layerBorder || c <= layerBorder || c >= t->width - layerBorder)
return 0;
// check the candidate point in the middle layer is above thresh
float candidate = m->getResponse(r, c, t);
if (candidate < thresh) //小于某一閾值也認(rèn)定為不是極值點
return 0;
for (int rr = -1; rr <=1; ++rr)
{
for (int cc = -1; cc <=1; ++cc)
{
// if any response in 3x3x3 is greater candidate not maximum<BR> //只要3x3x3鄰域存在一點大于當(dāng)前點的response值，只可認(rèn)定該點非極值點
if (
t->getResponse(r+rr, c+cc) >= candidate ||
((rr != 0 || cc != 0) && m->getResponse(r+rr, c+cc, t) >= candidate) ||
b->getResponse(r+rr, c+cc, t) >= candidate
) //優(yōu)先級&&>||
return 0;
}
}
return 1;
}

　　大家務(wù)必注意，由于各層大小不一致，顧在比較大小的時候，要統(tǒng)一到統(tǒng)一尺度下，在getResponse()里面有所體現(xiàn)，即先計算兩者尺度之比，比如尺度之比為2，最上層當(dāng)前點為（20，25），那么需要比較的緊鄰下層點為（40,50）而不是下層的（20,25）。再看若當(dāng)前點為極值點，那么調(diào)用interpolateExtremum():

//! Interpolate scale-space extrema to subpixel accuracy to form an image feature.
void FastHessian::interpolateExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)
{
// get the step distance between filters
// check the middle filter is mid way between top and bottom
int filterStep = (m->filter - b->filter);
assert(filterStep > 0 && t->filter - m->filter == m->filter - b->filter);
// Get the offsets to the actual location of the extremum
double xi = 0, xr = 0, xc = 0;
interpolateStep(r, c, t, m, b, &xi, &xr, &xc );
// If point is sufficiently close to the actual extremum
if( fabs( xi ) < 0.5f && fabs( xr ) < 0.5f && fabs( xc ) < 0.5f )
{
Ipoint ipt;
ipt.x = static_cast<float>((c + xc) * t->step);
ipt.y = static_cast<float>((r + xr) * t->step);
ipt.scale = static_cast<float>((0.1333f) * (m->filter + xi * filterStep));
ipt.laplacian = static_cast<int>(m->getLaplacian(r,c,t));
ipts.push_back(ipt);
}
}

　　本函數(shù)實現(xiàn)功能，利用插值精確計算極值點在原圖像中的位置，并保存在ipt中。

打開interpolateStep，其中deriv3D是求當(dāng)前點的3的方向上的一階導(dǎo)，hessian3D是求當(dāng)前點的3維hessian二階導(dǎo)，最后計算出3個方向的偏差值xi,xr,xc .

void FastHessian::interpolateStep(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b,
double* xi, double* xr, double* xc )
{
CvMat* dD, * H, * H_inv, X;
double x[3] = { 0 };
dD = deriv3D( r, c, t, m, b );
H = hessian3D( r, c, t, m, b );
H_inv = cvCreateMat( 3, 3, CV_64FC1 );
cvInvert( H, H_inv, CV_SVD );//用svd法求逆矩陣H_inv
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );//新建X
cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );//X=-dD*H_inv
cvReleaseMat( &dD );
cvReleaseMat( &H );
cvReleaseMat( &H_inv );
*xi = x[2];
*xr = x[1];
*xc = x[0];
}

　　下面開始在特征點周圍提取特征描述符

④ 提取特征描述符

轉(zhuǎn)到surf.cpp里尋找getDescriptors()，由于upright初始化設(shè)置為false（為特征點分配主方向，并旋轉(zhuǎn)找到鄰域提取特征描述符），若為true,則直接提取鄰域特征描述符。

//! Describe all features in the supplied vector
void Surf::getDescriptors(bool upright)
{
// Check there are Ipoints to be described
if (!ipts.size()) return;
// Get the size of the vector for fixed loop bounds
int ipts_size = (int)ipts.size();
if (upright)
{
// U-SURF loop just gets descriptors
for (int i = 0; i < ipts_size; ++i)
{
// Set the Ipoint to be described
index = i;
// Extract upright (i.e. not rotation invariant) descriptors
getDescriptor(true);
}
}
else
{
// Main SURF-64 loop assigns orientations and gets descriptors
for (int i = 0; i < ipts_size; ++i)
{
// Set the Ipoint to be described
index = i;
// Assign Orientations and extract rotation invariant descriptors
getOrientation();
getDescriptor(false);
}
}
}

　　其實兩者區(qū)別在于提取特征描述符之前判斷upright的時候是否需要多調(diào)用一句getOrientation()就是了。前者不考慮旋轉(zhuǎn)，兩幅圖只是尺度有異，而后者還考慮了旋轉(zhuǎn)不變性，更加全面。

幾個地方：

1.如論文提出的采取r=6的圓域，共計包含109個像素點，大家可以算算，在此不多解釋了。

2.像素點的方向由harrx和harry計算得到,相當(dāng)于梯度公式里面的dx和dy，然后利用getAngle得到θ（分布在0~2∏，而不是簡單的tan-1(harrx/harry)取值在0~∏）。

//! Calculate Haar wavelet responses in x direction
inline float Surf::haarX(int row, int column, int s)
{
return BoxIntegral(img, row-s/2, column, s, s/2)
-1 * BoxIntegral(img, row-s/2, column-s/2, s, s/2);
}
//-------------------------------------------------------
//! Calculate Haar wavelet responses in y direction
inline float Surf::haarY(int row, int column, int s)
{
return BoxIntegral(img, row, column-s/2, s/2, s)
-1 * BoxIntegral(img, row-s/2, column-s/2, s/2, s);
}
float Surf::getAngle(float X, float Y)//計算每個點的方向角
{
if(X > 0 && Y >= 0)
return atan(Y/X);
if(X < 0 && Y >= 0)
return pi - atan(-Y/X);
if(X < 0 && Y < 0)
return pi + atan(Y/X);
if(X > 0 && Y < 0)
return 2*pi - atan(-Y/X);
return 0;
}

　　圖示如右： harr （黑色代表-1，白色為1，即白色區(qū)域積分值減去黑色區(qū)域積分值）

在getOrienation()里面resx和resy分別是上面計算出來的harrx和harry分別乘上高斯模板系數(shù)gauss25。

void Surf::getOrientation()
{
......
//將像素點按方向角分配到6個寬為60度的區(qū)間去，比如說可以將80度分配到ang1=90度，因為90-30<80<90+30
// loop slides pi/3 window around feature point
for(ang1 = 0; ang1 < 2*pi; ang1+=0.15f) {
ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);//保證ang1+60 不會大于360度
sumX = sumY = 0.f;
for(unsigned int k = 0; k < Ang.size(); ++k)//對于半徑為6的圓鄰域，這里面的Ang.size()=109
{
// get angle from the x-axis of the sample point
const float & ang = Ang[k];
// determine whether the point is within the window
if (ang1 < ang2 && ang1 < ang && ang < ang2)
{
sumX+=resX[k];
sumY+=resY[k];
}
else if (ang2 < ang1 &&
((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*pi) ))
{
sumX+=resX[k];
sumY+=resY[k];
}
}
// if the vector produced from this window is longer than all
// previous vectors then this forms the new dominant direction
if (sumX*sumX + sumY*sumY > max)//尋找一個ang1使得角度在此區(qū)間的長度最大 ，然后計算出累加的resx和resy，得出的方向角
{
// store largest orientation
max = sumX*sumX + sumY*sumY;
orientation = getAngle(sumX, sumY);
}
}
// assign orientation of the dominant response vector
ipt->orientation = orientation;
}

　　好了，終于要開始提取特征描述符了哈~.~

void Surf::getDescriptor(bool bUpright)
{
int y, x, sample_x, sample_y, count=0;
int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
float scale, *desc, dx, dy, mdx, mdy, co, si;
float gauss_s1 = 0.f, gauss_s2 = 0.f;
float rx = 0.f, ry = 0.f, rrx = 0.f, rry = 0.f, len = 0.f;
float cx = -0.5f, cy = 0.f; //Subregion centers for the 4x4 gaussian weighting
Ipoint *ipt = &ipts[index];
scale = ipt->scale;//尺度
x = fRound(ipt->x);
y = fRound(ipt->y);//空間位置
desc = ipt->descriptor;
if (bUpright)
{//不需旋轉(zhuǎn)的情況
co = 1;
si = 0;
}
else
{//需要旋轉(zhuǎn)調(diào)整選取鄰域的情況
co = cos(ipt->orientation);
si = sin(ipt->orientation);
}
i = -8;
//Calculate descriptor for this interest point
while(i < 12)
{
j = -8;
i = i-4;
cx += 1.f;
cy = -0.5f;
while(j < 12)
{
dx=dy=mdx=mdy=0.f;//特征向量的構(gòu)成
cy += 1.f;
j = j - 4;
ix = i + 5;//ix=i0+1,這里面i0和j0分別取得的值為-8，-3，2，7
jx = j + 5;//iy=j0+1
xs = fRound(x + ( -jx*scale*si + ix*scale*co));
ys = fRound(y + ( jx*scale*co + ix*scale*si));
//fRound為4舍五入算法，最近鄰插值尋找旋轉(zhuǎn)對應(yīng)點
for (int k = i; k < i + 9; ++k)
{
for (int l = j; l < j + 9; ++l)
{
//Get coords of sample point on the rotated axis
sample_x = fRound(x + (-l*scale*si + k*scale*co));
sample_y = fRound(y + ( l*scale*co + k*scale*si));
//Get the gaussian weighted x and y responses
gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);
rx = haarX(sample_y, sample_x, 2*fRound(scale));
ry = haarY(sample_y, sample_x, 2*fRound(scale));
//Get the gaussian weighted x and y responses on rotated axis
rrx = gauss_s1*(-rx*si + ry*co);
rry = gauss_s1*(rx*co + ry*si);
dx += rrx;
dy += rry;
mdx += fabs(rrx);
mdy += fabs(rry);
}
}
//Add the values to the descriptor vector
gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);
desc[count++] = dx*gauss_s2;
desc[count++] = dy*gauss_s2;
desc[count++] = mdx*gauss_s2;
desc[count++] = mdy*gauss_s2;
len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy) * gauss_s2*gauss_s2;
j += 9;
}
i += 9;
}
//Convert to Unit Vector 特征向量歸一化
len = sqrt(len);
for(int i = 0; i < 64; ++i)
desc[i] /= len;
}

　　其中i和j分別取的值為-8，-3，2，7，很明顯i,j確定的鄰域為7-（-8）+1=16，16x16的鄰域，旋轉(zhuǎn)對應(yīng)的在原圖像的點位置為

xs = fRound(x + ( -jx*scale*si + ix*scale*co)); ys = fRound(y + ( jx*scale*co + ix*scale*si));

co = cos(ipt->orientation); si = sin(ipt->orientation); 【ipt->orientation為特征點的方向角】

共有 4x4個子塊（9x9=81個像素點），每個子塊分別計算了其中16個dx,dy,|dx|,|dy|之和（當(dāng)然還要考慮高斯濾波權(quán)重系數(shù)），則最終的特征描述符為4x4x4=64維向量。

main.cpp內(nèi)mainImage函數(shù)內(nèi)部drawIpoints(img, ipts)就不用再做解釋了吧。

搞了一天，終于搞完了~~

來個截圖~~

謝謝大家~

轉(zhuǎn)載請注明blue_lghttp://www.cnblogs.com/blue-lg/archive/2012/07/20/2600436.html

    本站是提供個人知識管理的網(wǎng)絡(luò)存儲空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點。請注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購買等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請點擊一鍵舉報。
    轉(zhuǎn)藏 分享 獻花（0）

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多