SURF C++代码详细阅读（三）
阅读（二）获取了特征点，接着针对每个特征点，来构造一个SURF描述符。
3 SURF描述符
回顾主函数，1得到特征点后
使用图像和图像的所有特征点向量集合ipts2 构造SURF对象des
接着3 通过函数getDescriptors(upright)得到其描述符。
inline void surfDetDes(IplImage *img, /* image to find Ipoints in */
std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
int single_mem_cpy,
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, single_mem_cpy, octaves, intervals, init_sample, thres);
// Extract interest points and store in vector ipts
clock_t start = clock();
//1得到特征点
fh.getIpoints();
clock_t end = clock();
std::cout<< "Extract keypoints took: " << float(end - start) / CLOCKS_PER_SEC << " seconds" << std::endl;
// Create Surf Descriptor Object
//2构造SURF对象des
Surf des(int_img, ipts);
// Extract the descriptors for the ipts
start = clock();
//3 通过函数getDescriptors(upright)得到其描述符
des.getDescriptors(upright);
end = clock();
std::cout<< "Extract descriptor took: " << float(end - start) / CLOCKS_PER_SEC << " seconds" << std::endl;
// Deallocate the integral image
cvReleaseImage(&int_img);
}
要构造SURF对象，先看SURF类
class Surf {
public:
//! Standard Constructor (img is an integral image)
Surf(IplImage *img, std::vector<Ipoint> &ipts);
//! Describe all features in the supplied vector
void getDescriptors(bool bUpright = false);
private:
//---------------- Private Functions -----------------//
//! Assign the current Ipoint an orientation
void getOrientation(int index);
//! Get the descriptor. See Agrawal ECCV 08
void getDescriptor(bool bUpright = false, int index = 0);
//! Calculate the value of the 2d gaussian at x,y
inline float gaussian(int x, int y, float sig);
inline float gaussian(float x, float y, float sig);
//! Calculate Haar wavelet responses in x and y directions
inline float haarX(int row, int column, int size);
inline float haarY(int row, int column, int size);
//! Get the angle from the +ve x-axis of the vector given by [X Y]
float getAngle(float X, float Y);
//---------------- Private Variables -----------------//
//! Integral image where Ipoints have been detected
IplImage *img;
//! Ipoints vector
IpVec &ipts;
//! Index of current Ipoint in the vector
// int index;
};
3.1 构造函数
其中图像参数使用了赋值运算符初始化，特征点使用了括号赋值初始化（C++特性）
Surf::Surf(IplImage *img, IpVec &ipts)
: ipts(ipts)
{
this->img = img;
}
括号赋值是在给变量分配内存空间的同时对它进行初始化，而赋值运算符赋值则是先分配内存空间，然后再初始化，这里将特征点ipts（引用数据成员）在构造函数调入之后函数体执行之前就进行初始化。
3.2 获取描述符函数（主要）
getDescriptors(bool bUpright = false); 其中使用了SURF类中的高斯滤波、哈尔小波、获取角度函数。
在函数体中，大段描述的是不旋转的描述符，而SURF描述符是考虑旋转的，在得到描述符的邻域矩形计算过程时，要先计算主方向，这里直接执行else。
对所给图像中每一个特征点，得到其主方向getOrientation(i);，再得到其描述符getDescriptor(false, i);这里得到描述符的函数重载为两个参数，false的话是带旋转情况，i是每个特征点索引，便于对每个特征点找到主方向并进行旋转。
else
{
// Main SURF-64 loop assigns orientations and gets descriptors
// #pragma acc parallel loop
for (int i = 0; i < ipts_size; ++i)
{
// Set the Ipoint to be described
// index = i;
// Assign Orientations and extract rotation invariant descriptors
getOrientation(i);
getDescriptor(false, i);
}
}
3.2.1 获取主方向
//! Assign the supplied Ipoint an orientation
void Surf::getOrientation(int index)
{
Ipoint *ipt = &ipts[index];
float gauss = 0.f, scale = ipt->scale;
const int s = fRound(scale), r = fRound(ipt->y), c = fRound(ipt->x);
std::vector<float> resX(109), resY(109), Ang(109);//容量 半径为6的园空间内离散的像素点数109
const int id[] = {6,5,4,3,2,1,0,1,2,3,4,5,6};
int idx = 0;
// calculate haar responses for points within radius of 6*scale
for(int i = -6; i <= 6; ++i)
{
for(int j = -6; j <= 6; ++j)
{
if(i*i + j*j < 36)
{
// gauss = static_cast<float>(gauss33[id[i+6]][id[j+6]]);
gauss = gauss33[id[i+6]][id[j+6]];
resX[idx] = gauss * haarX(r+j*s, c+i*s, 4*s);
resY[idx] = gauss * haarY(r+j*s, c+i*s, 4*s);
Ang[idx] = getAngle(resX[idx], resY[idx]);
++idx;
}
}
}
// calculate the dominant direction
float sumX=0.f, sumY=0.f;
float max=0.f, orientation = 0.f;
float ang1=0.f, ang2=0.f;
//以特征点为中心，张角为π/3的扇形滑动，每次弧度增加0.15，
// loop slides pi/3 window around feature point
for(ang1 = 0; ang1 < 2*pi; ang1+=0.15f) {//每次弧度增加0.15
ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);//覆盖角度pi/3
sumX = sumY = 0.f;
for(unsigned int k = 0; k < Ang.size(); ++k)
{
// 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];
}
}
//计算窗口内的Harr小波响应值dx、dy的累加，选取最大值作为主方向，通过反三角函数getAngle(sumX, sumY)得到主方向角度
// 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)
{
// store largest orientation
max = sumX*sumX + sumY*sumY;
orientation = getAngle(sumX, sumY);
}
}
// assign orientation of the dominant response vector
ipt->orientation = orientation;
}
对每一个特征点，四舍五入fround,得到对应横纵坐标、尺度 r、c、s
接着在像素位置上，进行高斯模糊和哈尔小波运算，哈尔小波计算运用了积分图像。
这里举X方向哈尔小波HarrX例子来看
inline float Surf::haarX(int row, int column, int s)//积分图参数（r+j*s, c+i*s, 4*s）
{
return BoxIntegral(img, row-s/2, column, s, s/2)
-1 * BoxIntegral(img, row-s/2, column-s/2, s, s/2);
}
得到尺度4s的X方向哈尔滤波器
在经过X、Y方向的哈尔小波后，得到像素的差值dx、dy，接着以特征点为中心，张角为π/3的扇形滑动，每次弧度增加0.15，计算窗口内的Harr小波响应值dx、dy的累加，选取累加最大值方向作为主方向，通过反三角函数getAngle(sumX, sumY)得到角度，返回给ipt>orientation。
3.2.1 计算描述符（旋转）
//! Get the modified descriptor. See Agrawal ECCV 08
//! Modified descriptor contributed by Pablo Fernandez
void Surf::getDescriptor(bool bUpright, int index)
{
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)
{
co = 1;
si = 0;
}
else
{
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;
cy += 1.f;
j = j - 4;
ix = i + 5;
jx = j + 5;
xs = fRound(x + ( -jx*scale*si + ix*scale*co));
ys = fRound(y + ( jx*scale*co + ix*scale*si));
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;
}
为了充分利用积分图计算，先使用水平和垂直的哈尔小波滤波器求得dx和dy，然后根据主方向旋转dx和dy与主方向操持一致。
以特征点为中心，在20s邻域矩形划分为4×4个子块，每个子块大小为20s20s，每个子块有5s×5s个像素
每个子块利用尺寸2s哈尔小波滤波器（积分图计算）和高斯滤波器，每个子块统计25个像素的水平方向和垂直方向的哈尔小波特征，得到∑dx、∑|dx|、∑dy、∑|dy|，44个小块具有444=64维数据，这就得到了SURF的64维描述符。
如下图2所示。图中，以特征点为中心，以20s为边长的矩形窗口为特征描述子计算使用的窗口，特征点到矩形边框的线段表示特征点的主方向。
像素点围绕另一个像素点旋转公式：
x1=(x-xb)cosθ - (y-yb)sinθ + xb
y1=(y-yb)cosθ + (x-xb)sinθ + yb
旋转原理可以参考https://guo-pu.blog.csdn.net/article/details/116196097
代码中每个子块中的像素（x-jxs,y-ixs）（图右）围绕中心点（x,y）旋转对应主方向的角度，后四舍五入到像素位，得到旋转后的坐标(xs,ys)（图左）。
xs = fRound(x + ( -jx*scale*si + ix*scale*co));
ys = fRound(y + ( jx*scale*co + ix*scale*si));
最后将特征值desc归一化，得到了SURF特征。
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
