struct SURFInvoker
{enum{ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20};// Parametersconst Mat* img;const Mat* sum;vector<KeyPoint>* keypoints;Mat* descriptors;bool extended;bool upright;// Pre-calculated valuesint nOriSamples;vector<Point> apt; // 特征点周围用于描述方向的邻域的点vector<float> aptw; // 描述 方向时的 高斯 权vector<float> DW;SURFInvoker(const Mat& _img, const Mat& _sum,vector<KeyPoint>& _keypoints, Mat& _descriptors,bool _extended, bool _upright){keypoints = &_keypoints;descriptors = &_descriptors;img = &_img;sum = &_sum;extended = _extended;upright = _upright;// 用于描述特征点的 方向的 邻域大小： 12*sigma+1 (sigma =1.2) 因为高斯加权的核的参数为2sigma// nOriSampleBound为 矩形框内点的个数const int nOriSampleBound = (2 * ORI_RADIUS + 1)*(2 * ORI_RADIUS + 1); // 这里把s近似为1 ORI_DADIUS = 6s// 分配大小
        apt.resize(nOriSampleBound);aptw.resize(nOriSampleBound);DW.resize(PATCH_SZ*PATCH_SZ); // PATHC_SZ为特征描述子的 区域大小 20s(s 这里初始为1了)/* 计算特征点方向用的 高斯分布 权值与坐标 */Mat G_ori = getGaussianKernel(2 * ORI_RADIUS + 1, SURF_ORI_SIGMA, CV_32F); // SURF_ORI_SIGMA = 1.2 *2 =2.5nOriSamples = 0;for (int i = -ORI_RADIUS; i <= ORI_RADIUS; i++){for (int j = -ORI_RADIUS; j <= ORI_RADIUS; j++){if (i*i + j*j <= ORI_RADIUS*ORI_RADIUS) // 限制在圆形区域内
                {apt[nOriSamples] = cvPoint(i, j);// 下面这里有个坐标转换，因为i,j都是从-ORI_RADIUS开始的。aptw[nOriSamples++] = G_ori.at<float>(i + ORI_RADIUS, 0) * G_ori.at<float>(j + ORI_RADIUS, 0);}}}CV_Assert(nOriSamples <= nOriSampleBound); // nOriSamples为圆形区域内的点，nOriSampleBound是正方形区域的点/* 用于特征点描述子的高斯 权值 */Mat G_desc = getGaussianKernel(PATCH_SZ, SURF_DESC_SIGMA, CV_32F); // 用于生成特征描述子的 高斯加权 sigma = 3.3s (s初取1)for (int i = 0; i < PATCH_SZ; i++){for (int j = 0; j < PATCH_SZ; j++)DW[i*PATCH_SZ + j] = G_desc.at<float>(i, 0) * G_desc.at<float>(j, 0);}/* x与y方向上的 Harr小波,参数为4s */const int NX = 2, NY = 2;const int dx_s[NX][5] = { { 0, 0, 2, 4, -1 }, { 2, 0, 4, 4, 1 } };const int dy_s[NY][5] = { { 0, 0, 4, 2, 1 }, { 0, 2, 4, 4, -1 } };float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound]; // 用于计算特生点主方向uchar PATCH[PATCH_SZ + 1][PATCH_SZ + 1];float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ]; // 20s * 20s区域的 梯度值CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);Mat _patch(PATCH_SZ + 1, PATCH_SZ + 1, CV_8U, PATCH);int dsize = extended ? 128 : 64;int k, k1 = 0, k2 = (int)(*keypoints).size();// k2为Harr小波的 模板尺寸float maxSize = 0;for (k = k1; k < k2; k++){maxSize = std::max(maxSize, (*keypoints)[k].size);}// maxSize*1.2/9 表示最大的尺度 sint imaxSize = std::max(cvCeil((PATCH_SZ + 1)*maxSize*1.2f / 9.0f), 1);Ptr<CvMat> winbuf = cvCreateMat(1, imaxSize*imaxSize, CV_8U);for (k = k1; k < k2; k++){int i, j, kk, nangle;float* vec;SurfHF dx_t[NX], dy_t[NY];KeyPoint& kp = (*keypoints)[k];float size = kp.size;Point2f center = kp.pt;/* s是当前层的尺度参数 1.2是第一层的参数，9是第一层的模板大小*/float s = size*1.2f / 9.0f;/* grad_wav_size是 harr梯度模板的大小 边长为 4s */int grad_wav_size = 2 * cvRound(2 * s);if (sum->rows < grad_wav_size || sum->cols < grad_wav_size){/* when grad_wav_size is too big,* the sampling of gradient will be meaningless* mark keypoint for deletion. */kp.size = -1;continue;}float descriptor_dir = 360.f - 90.f;if (upright == 0){// 这一步 是计算梯度值，先将harr模板放大，再根据积分图计算，与前面求D_x，D_y一致类似resizeHaarPattern(dx_s, dx_t, NX, 4, grad_wav_size, sum->cols);resizeHaarPattern(dy_s, dy_t, NY, 4, grad_wav_size, sum->cols);for (kk = 0, nangle = 0; kk < nOriSamples; kk++){int x = cvRound(center.x + apt[kk].x*s - (float)(grad_wav_size - 1) / 2);int y = cvRound(center.y + apt[kk].y*s - (float)(grad_wav_size - 1) / 2);if (y < 0 || y >= sum->rows - grad_wav_size ||x < 0 || x >= sum->cols - grad_wav_size)continue;const int* ptr = &sum->at<int>(y, x);float vx = calcHaarPattern(ptr, dx_t, 2);float vy = calcHaarPattern(ptr, dy_t, 2);X[nangle] = vx*aptw[kk];Y[nangle] = vy*aptw[kk];nangle++;}if (nangle == 0){// No gradient could be sampled because the keypoint is too// near too one or more of the sides of the image. As we// therefore cannot find a dominant direction, we skip this// keypoint and mark it for later deletion from the sequence.kp.size = -1;continue;}matX.cols = matY.cols = _angle.cols = nangle;// 计算邻域内每个点的 梯度角度cvCartToPolar(&matX, &matY, 0, &_angle, 1);float bestx = 0, besty = 0, descriptor_mod = 0;for (i = 0; i < 360; i += SURF_ORI_SEARCH_INC) // SURF_ORI_SEARCH_INC 为扇形区域扫描的步长
                {float sumx = 0, sumy = 0, temp_mod;for (j = 0; j < nangle; j++){// d是 分析到的那个点与 现在主方向的偏度int d = std::abs(cvRound(angle[j]) - i);if (d < ORI_WIN / 2 || d > 360 - ORI_WIN / 2){sumx += X[j];sumy += Y[j];}}temp_mod = sumx*sumx + sumy*sumy;// descriptor_mod 是最大峰值if (temp_mod > descriptor_mod){descriptor_mod = temp_mod;bestx = sumx;besty = sumy;}}descriptor_dir = fastAtan2(-besty, bestx);}kp.angle = descriptor_dir;if (!descriptors || !descriptors->data)continue;/* 用特征点周围20*s为边长的邻域 计算特征描述子 */int win_size = (int)((PATCH_SZ + 1)*s);CV_Assert(winbuf->cols >= win_size*win_size);Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);if (!upright){descriptor_dir *= (float)(CV_PI / 180); // 特征点的主方向 弧度值float sin_dir = -std::sin(descriptor_dir); //  - sin dirfloat cos_dir = std::cos(descriptor_dir);float win_offset = -(float)(win_size - 1) / 2;float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;uchar* WIN = win.data;int ncols1 = img->cols - 1, nrows1 = img->rows - 1;size_t imgstep = img->step;for (i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir){double pixel_x = start_x;double pixel_y = start_y;for (j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir){int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);if ((unsigned)ix < (unsigned)ncols1 &&(unsigned)iy < (unsigned)nrows1){float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);const uchar* imgptr = &img->at<uchar>(iy, ix);WIN[i*win_size + j] = (uchar)cvRound(imgptr[0] * (1.f - a)*(1.f - b) +imgptr[1] * a*(1.f - b) +imgptr[imgstep] * (1.f - a)*b +imgptr[imgstep + 1] * a*b);}else{int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);WIN[i*win_size + j] = img->at<uchar>(y, x);}}}}else{float win_offset = -(float)(win_size - 1) / 2;int start_x = cvRound(center.x + win_offset);int start_y = cvRound(center.y - win_offset);uchar* WIN = win.data;for (i = 0; i < win_size; i++, start_x++){int pixel_x = start_x;int pixel_y = start_y;for (j = 0; j < win_size; j++, pixel_y--){int x = MAX(pixel_x, 0);int y = MAX(pixel_y, 0);x = MIN(x, img->cols - 1);y = MIN(y, img->rows - 1);WIN[i*win_size + j] = img->at<uchar>(y, x);}}}// Scale the window to size PATCH_SZ so each pixel's size is s. This// makes calculating the gradients with wavelets of size 2s easyresize(win, _patch, _patch.size(), 0, 0, INTER_AREA);// Calculate gradients in x and y with wavelets of size 2sfor (i = 0; i < PATCH_SZ; i++)for (j = 0; j < PATCH_SZ; j++){float dw = DW[i*PATCH_SZ + j]; // 高斯加权系数float vx = (PATCH[i][j + 1] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i + 1][j])*dw;float vy = (PATCH[i + 1][j] - PATCH[i][j] + PATCH[i + 1][j + 1] - PATCH[i][j + 1])*dw;DX[i][j] = vx;DY[i][j] = vy;}// Construct the descriptorvec = descriptors->ptr<float>(k);for (kk = 0; kk < dsize; kk++)vec[kk] = 0;double square_mag = 0;if (extended){// 128维描述子，考虑dx与dy的正负号for (i = 0; i < 4; i++)for (j = 0; j < 4; j++){// 每个方块内是一个5s * 5s的区域，每个方法由8个特征描述for (int y = i * 5; y < i * 5 + 5; y++){for (int x = j * 5; x < j * 5 + 5; x++){float tx = DX[y][x], ty = DY[y][x];if (ty >= 0){vec[0] += tx;vec[1] += (float)fabs(tx);}else {vec[2] += tx;vec[3] += (float)fabs(tx);}if (tx >= 0){vec[4] += ty;vec[5] += (float)fabs(ty);}else {vec[6] += ty;vec[7] += (float)fabs(ty);}}}for (kk = 0; kk < 8; kk++)square_mag += vec[kk] * vec[kk];vec += 8;}}else{// 64位描述子for (i = 0; i < 4; i++)for (j = 0; j < 4; j++){for (int y = i * 5; y < i * 5 + 5; y++){for (int x = j * 5; x < j * 5 + 5; x++){float tx = DX[y][x], ty = DY[y][x];vec[0] += tx; vec[1] += ty;vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);}}for (kk = 0; kk < 4; kk++)square_mag += vec[kk] * vec[kk];vec += 4;}}// 归一化 描述子 以满足 光照不变性vec = descriptors->ptr<float>(k);float scale = (float)(1. / (sqrt(square_mag) + DBL_EPSILON));for (kk = 0; kk < dsize; kk++)vec[kk] *= scale;}}
};