SIFT系列01--SIFT算法OpenCV源码

#include "precomp.hpp"

#include <iostream>

#include <stdarg.h>

namespace cv

{

// default number of sampled intervals per octave

static const int SIFT_INTVLS = 3;

// default sigma for initial gaussian smoothing

static const float SIFT_SIGMA = 1.6f;

// default threshold on keypoint contrast |D(x)|

static const float SIFT_CONTR_THR = 0.04f;

// default threshold on keypoint ratio of principle curvatures

static const float SIFT_CURV_THR = 10.f;

// double image size before pyramid construction?

static const bool SIFT_IMG_DBL = true;

// default width of descriptor histogram array

static const int SIFT_DESCR_WIDTH = 4;

// default number of bins per histogram in descriptor array

static const int SIFT_DESCR_HIST_BINS = 8;

// assumed gaussian blur for input image

static const float SIFT_INIT_SIGMA = 0.5f;

// width of border in which to ignore keypoints

static const int SIFT_IMG_BORDER = 5;

// maximum steps of keypoint interpolation before failure

static const int SIFT_MAX_INTERP_STEPS = 5;

// default number of bins in histogram for orientation assignment

static const int SIFT_ORI_HIST_BINS = 36;

// determines gaussian sigma for orientation assignment

static const float SIFT_ORI_SIG_FCTR = 1.5f;

// determines the radius of the region used in orientation assignment

static const float SIFT_ORI_RADIUS = 3 * SIFT_ORI_SIG_FCTR;

// orientation magnitude relative to max that results in new feature

static const float SIFT_ORI_PEAK_RATIO = 0.8f;

// determines the size of a single descriptor orientation histogram

static const float SIFT_DESCR_SCL_FCTR = 3.f;

// threshold on magnitude of elements of descriptor vector

static const float SIFT_DESCR_MAG_THR = 0.2f;

// factor used to convert floating-point descriptor to unsigned char

static const float SIFT_INT_DESCR_FCTR = 512.f;

static const int SIFT_FIXPT_SCALE = 48;

//获取原始图像，把读取的图像长宽变为原来的2倍

static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma )

{

      Mat gray, gray_fpt;

      if( img.channels() == 3 || img.channels() == 4 )

            cvtColor(img, gray, COLOR_BGR2GRAY);

      else

            img.copyTo(gray);

      gray.convertTo(gray_fpt, CV_16S, SIFT_FIXPT_SCALE, 0);

      float sig_diff;

      if( doubleImageSize )

      {

            sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );

            Mat dbl;

            resize(gray_fpt, dbl, Size(gray.cols*2, gray.rows*2), 0, 0, INTER_LINEAR);

            GaussianBlur(dbl, dbl, Size(), sig_diff, sig_diff);

            return dbl;

      }

      else

      {

            sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );

            GaussianBlur(gray_fpt, gray_fpt, Size(), sig_diff, sig_diff);

            return gray_fpt;

      }

}

//构建nOctaves组高斯金字塔，每组有nOctaveLayers +3层，nOctaveLayers应该是每组Dog要用到的层数

void SIFT::buildGaussianPyramid( const Mat& base, vector& pyr, int nOctaves ) const

{

      vector sig(nOctaveLayers + 3);

      pyr.resize(nOctaves*(nOctaveLayers + 3));   //pyr保存所有组所有层

      // precompute Gaussian sigmas using the following formula:

      //   \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2

//计算每层的尺度因子sig[i]

      sig[0] = sigma;                        //第0层尺度为sigma

      double k = pow( 2., 1. / nOctaveLayers );

      for( int i = 1; i < nOctaveLayers + 3; i++ )

      {

            double sig_prev = pow(k, (double)(i-1))*sigma;

            double sig_total = sig_prev*k;

            sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);

      }

      for( int o = 0; o < nOctaves; o++ )

      {

            for( int i = 0; i < nOctaveLayers + 3; i++ )

            {

                  Mat& dst = pyr[o*(nOctaveLayers + 3) + i];

//第0组第0层为base层，即原始图像

                  if( o == 0   &&   i == 0 )

                        dst = base;

                  // base of new octave is halved image from end of previous octave

//高斯金字塔的新组(new octave)的第0幅为上一组的第nOctaveLayers幅下采样得到，采样步长为2

                  else if( i == 0 )

                  {

                        const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];

                        resize(src, dst, Size(src.cols/2, src.rows/2),

                                   0, 0, INTER_NEAREST);

                  }

// 每一组的第i幅图像是由第i-1幅图像用sig[i]高斯模糊得到，即得到尺度为sig[i]的图像

                  else

                  {

                        const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1];

                        GaussianBlur(src, dst, Size(), sig[i], sig[i]);

                  }

            }

      }

}

//构建Dog金字塔

void SIFT::buildDoGPyramid( const vector& gpyr, vector& dogpyr ) const

{

      int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3);    //nOctaves表示组的个数

      dogpyr.resize( nOctaves*(nOctaveLayers + 2) );      //保存所有组的Dog图像

//每组相邻两幅图像相减，获取Dog图像

      for( int o = 0; o < nOctaves; o++ )

      {

            for( int i = 0; i < nOctaveLayers + 2; i++ )

            {

                  const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i];

                  const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1];

                  Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i];

                  subtract(src2, src1, dst, noArray(), CV_16S);    //两幅图像相减

            }

      }

}

//计算某一个特征点的周围区域梯度方向直方图

// Computes a gradient orientation histogram at a specified pixel

static float calcOrientationHist( const Mat& img, Point pt, int radius,

                                                   float sigma, float* hist, int n )

{

//len表示以2*radius+1为半径的圆(因为点是离散的，其实为正方形)的像素个数

      int i, j, k, len = (radius*2+1)*(radius*2+1);

      float expf_scale = -1.f/(2.f * sigma * sigma);

      AutoBuffer buf(len*4 + n+4);

      float *X = buf, *Y = X + len, *Mag = X, *Ori = Y + len, *W = Ori + len;

      float* temphist = W + len + 2;

      for( i = 0; i < n; i++ )

            temphist[i] = 0.f;

      for( i = -radius, k = 0; i <= radius; i++ )

      {

            int y = pt.y + i;

            if( y <= 0 || y >= img.rows - 1 )

                  continue;

            for( j = -radius; j <= radius; j++ )

            {

                  int x = pt.x + j;

                  if( x <= 0 || x >= img.cols - 1 )

                        continue;

                  float dx = (float)(img.at(y, x+1) - img.at(y, x-1));

                  float dy = (float)(img.at(y-1, x) - img.at(y+1, x));

                  X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale;

                  k++;

            }

      }

      len = k;

      // compute gradient values, orientations and the weights over the pixel neighborhood

      exp(W, W, len);

      fastAtan2(Y, X, Ori, len, true);    //计算向量角度

      magnitude(X, Y, Mag, len);            //计算梯度

      for( k = 0; k < len; k++ )

      {

            int bin = cvRound((n/360.f)*Ori[k]);

            if( bin >= n )

                  bin -= n;

            if( bin < 0 )

                  bin += n;

            temphist[bin] += W[k]*Mag[k];

      }

      // smooth the histogram

//平滑

      temphist[-1] = temphist[n-1];

      temphist[-2] = temphist[n-2];

      temphist[n] = temphist[0];

      temphist[n+1] = temphist[1];

      for( i = 0; i < n; i++ )

      {

            hist[i] = (temphist[i-2] + temphist[i+2])*(1.f/16.f) +

                  (temphist[i-1] + temphist[i+1])*(4.f/16.f) +

                  temphist[i]*(6.f/16.f);

      }

      float maxval = hist[0];

      for( i = 1; i < n; i++ )

            maxval = std::max(maxval, hist[i]);

      return maxval;       //返回直方图中最大值

}

//

// Interpolates a scale-space extremum's location and scale to subpixel

// accuracy to form an image feature.   Rejects features with low contrast.

// Based on Section 4 of Lowe's paper.

//   dog_pyr:dog金字塔;

//   kpt:关键点;

//   octv:组序号

//   layer: dog层序号

//   r: 行号; c:列号

//   nOctaveLayers:dog中要用到的层数，为3

//   contrastThreshold:对比度阈值=0.04

//    edgeThreshold:边界阈值=10

//   sigma: 尺度因子

static bool adjustLocalExtrema( const vector& dog_pyr, KeyPoint& kpt, int octv,

                                                int& layer, int& r, int& c, int nOctaveLayers,

                                                float contrastThreshold, float edgeThreshold, float sigma )

{

      const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE);    //SIFT_FIXPT_SCALE=48

      const float deriv_scale = img_scale*0.5f;       

      const float second_deriv_scale = img_scale;     

      const float cross_deriv_scale = img_scale*0.25f;     

      float xi=0, xr=0, xc=0, contr;

      int i = 0;

//循环5次，对坐标进行5次修正。

//SIFT_MAX_INTERP_STEPS=5，表示 maximum steps of keypoint interpolation before failure

      for( ; i < SIFT_MAX_INTERP_STEPS; i++ )

      {

            int idx = octv*(nOctaveLayers+2) + layer;

            const Mat& img = dog_pyr[idx];

            const Mat& prev = dog_pyr[idx-1];

            const Mat& next = dog_pyr[idx+1];

//计算一阶偏导数，通过临近点差分求得

            Vec3f dD((img.at(r, c+1) - img.at(r, c-1))*deriv_scale,

                          (img.at(r+1, c) - img.at(r-1, c))*deriv_scale,

                          (next.at(r, c) - prev.at(r, c))*deriv_scale);

//计算二阶偏导数，通过临近点差分求得

            float v2 = (float)img.at(r, c)*2;

            float dxx = (img.at(r, c+1) + img.at(r, c-1) - v2)*second_deriv_scale;

            float dyy = (img.at(r+1, c) + img.at(r-1, c) - v2)*second_deriv_scale;

            float dss = (next.at(r, c) + prev.at(r, c) - v2)*second_deriv_scale;

            float dxy = (img.at(r+1, c+1) - img.at(r+1, c-1) -

                                img.at(r-1, c+1) + img.at(r-1, c-1))*cross_deriv_scale;

            float dxs = (next.at(r, c+1) - next.at(r, c-1) -

                                prev.at(r, c+1) + prev.at(r, c-1))*cross_deriv_scale;

            float dys = (next.at(r+1, c) - next.at(r-1, c) -

                                prev.at(r+1, c) + prev.at(r-1, c))*cross_deriv_scale;

//二阶偏导数矩阵

            Matx33f H(dxx, dxy, dxs,

                           dxy, dyy, dys,

                           dxs, dys, dss);

//令尺度方程的泰勒展开导数为0，求解出偏移量X

            Vec3f X = H.solve(dD, DECOMP_LU);

            xi = -X[2];         //层偏移？？不应该是尺度偏移吗，层偏移即尺度偏移

            xr = -X[1];       //行偏移

            xc = -X[0];      //列偏移

//如果求解出的偏移量均小于0.5，则退出循环，说明该关键点的选取是正确的

            if( std::abs( xi ) < 0.5f   &&   std::abs( xr ) < 0.5f   &&   std::abs( xc ) < 0.5f )

                  break;

//求解出的偏移量均大于0.5，则将原坐标加上求出来的偏移量，得到更精确的坐标

            c += cvRound( xc );

            r += cvRound( xr );

            layer += cvRound( xi );      //？？？？

//层数或行列超出边界则退出

            if( layer < 1 || layer > nOctaveLayers ||

                 c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER   ||

                 r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER )

                  return false;

      }

    

//保证插值收敛，SIFT_MAX_INTERP_STEPS=5

      if( i >= SIFT_MAX_INTERP_STEPS )

            return false;

      {

            int idx = octv*(nOctaveLayers+2) + layer;   //求出修正后的层序号

            const Mat& img = dog_pyr[idx];

            const Mat& prev = dog_pyr[idx-1];

            const Mat& next = dog_pyr[idx+1];

            Matx31f dD((img.at(r, c+1) - img.at(r, c-1))*deriv_scale,

                             (img.at(r+1, c) - img.at(r-1, c))*deriv_scale,

                             (next.at(r, c) - prev.at(r, c))*deriv_scale);

//一阶偏导向量与上面求出的偏移量向量求点积

            float t = dD.dot(Matx31f(xc, xr, xi));

//修正后的坐标带入泰勒展开式，得到的结果小于0.04则抛弃该点

            contr = img.at(r, c)*img_scale + t * 0.5f;     

            if( std::abs( contr ) * nOctaveLayers < contrastThreshold )

                  return false;

          

//利用Hessian矩阵的迹和行列式计算该关键点的主曲率的比值

            float v2 = img.at(r, c)*2.f;

            float dxx = (img.at(r, c+1) + img.at(r, c-1) - v2)*second_deriv_scale;

            float dyy = (img.at(r+1, c) + img.at(r-1, c) - v2)*second_deriv_scale;

            float dxy = (img.at(r+1, c+1) - img.at(r+1, c-1) -

                                img.at(r-1, c+1) + img.at(r-1, c-1)) * cross_deriv_scale;

            float tr = dxx + dyy;

            float det = dxx * dyy - dxy * dxy;

//edgeThreshold=10,去除边缘点

            if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )

                  return false;

      }

//精确特征点在原图像上的位置

      kpt.pt.x = (c + xc) * (1 << octv);      //高斯金字塔坐标根据组数扩大相应的倍数

      kpt.pt.y = (r + xr) * (1 << octv);

      kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16);    //特征点被检测出时所在的金字塔组

      kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2;    //关键点邻域直径

      return true;

}

//

// Detects features at extrema in DoG scale space.   Bad features are discarded

// based on contrast and ratio of principal curvatures.

//找到尺度空间Dog图像中的极值点，即候选关键点

void SIFT::findScaleSpaceExtrema( const vector& gauss_pyr, const vector& dog_pyr,

                                                   vector& keypoints ) const

{

      int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);      //组的个数

// The contrast threshold used to filter out weak features in semi-uniform  

// (low-contrast) regions. The larger the threshold, the less features are produced by the detector.  

//低 对比度的阈值， contrastThreshold默认为0.04  

      int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);

      const int n = SIFT_ORI_HIST_BINS;   //直方图bin的个数=36，每个10度

      float hist[n];

      KeyPoint kpt;

      keypoints.clear();

      for( int o = 0; o < nOctaves; o++ )

            for( int i = 1; i <= nOctaveLayers; i++ )   // nOctaveLayers表示每层Dog图像个数-2，即最终用到的层数

            {

                  int idx = o*(nOctaveLayers+2)+i;             

                  const Mat& img = dog_pyr[idx];                //获取该层Dog图像，序号从1开始。第0层和最后一层不用

                  const Mat& prev = dog_pyr[idx-1];      //上一幅Dog图像

                  const Mat& next = dog_pyr[idx+1];    //下一幅Dog图像

                  int step = (int)img.step1();

                  int rows = img.rows, cols = img.cols;

//SIFT_IMG_BORDER=5，边界5个像素的距离

                  for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)     

                  {

//获取3幅相邻的Dog图像的行指针

                        const short* currptr = img.ptr(r);

                        const short* prevptr = prev.ptr(r);

                        const short* nextptr = next.ptr(r);

                        for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)

                        {

//获取该层图像r行c列的像素值

                              int val = currptr[c];

                              // find local extrema with pixel accuracy

//与周围26个点比较，极大或极小值则为局部极值点

                              if( std::abs(val) > threshold &&

                                   ((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&

                                      val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] &&

                                      val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] &&

                                      val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] &&

                                      val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] &&

                                      val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] &&

                                      val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] &&

                                      val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] &&

                                      val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||

                                    (val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&

                                      val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] &&

                                      val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] &&

                                      val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] &&

                                      val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] &&

                                      val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] &&

                                      val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] &&

                                      val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] &&

                                      val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1])))

                              {

//找到极值点之后，在Dog中调整局部极值点

                                    int r1 = r, c1 = c, layer = i;

                                    if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1,

                                                                        nOctaveLayers, (float)contrastThreshold,

                                                                        (float)edgeThreshold, (float)sigma) )

                                          continue;

//计算该特征点周围的方向直方图，并返回直方图最大值

                                    float scl_octv = kpt.size*0.5f/(1 << o);    //获取该特征点的尺度

//方向直方图的计算是在该尺度的高斯金字塔图像中计算的，不是在Dog图像，也不是在原图

                                    float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3) + layer],

                                                                                      Point(c1, r1),

                                                                                      cvRound(SIFT_ORI_RADIUS * scl_octv),   //直方图统计半径：3*1.5*σ，SIFT_ORI_RADIUS=3*1.5

                                                                                      SIFT_ORI_SIG_FCTR * scl_octv,

                                                                                      hist, n); 

                                    float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO);    //辅方向为0.8*主方向最大值，SIFT_ORI_PEAK_RATIO=0.8f

                                    for( int j = 0; j < n; j++ )

                                    {

                                          int l = j > 0 ? j - 1 : n - 1;

                                          int r2 = j < n-1 ? j + 1 : 0;

                                          if( hist[j] > hist[l]   &&   hist[j] > hist[r2]   &&   hist[j] >= mag_thr )

                                          {

                                                float bin = j + 0.5f * (hist[l]-hist[r2]) / (hist[l] - 2*hist[j] + hist[r2]);

                                                bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;

                                                kpt.angle = (float)((360.f/n) * bin);      //得到关键点的方向

                                                keypoints.push_back(kpt);       //这里保存的特征点具有位置，尺度和方向3个信息

                                          }

                                    }

                              }

                        }

                  }

            }

}

//img:特征点所在金字塔图像

//ptf:特征点在金字塔中的坐标

//ori:特征点角度

//scl:特征点邻域半径

//d=4

//n=8

//dst:第i个特征点描述符指针

static void calcSIFTDescriptor( const Mat& img, Point2f ptf, float ori, float scl,

                                               int d, int n, float* dst )

{

      Point pt(cvRound(ptf.x), cvRound(ptf.y));   //坐标点取整

      float cos_t = cosf(ori*(float)(CV_PI/180));   //余弦值

      float sin_t = sinf(ori*(float)(CV_PI/180));    //正弦值

      float bins_per_rad = n / 360.f;

      float exp_scale = -1.f/(d * d * 0.5f);

      float hist_width = SIFT_DESCR_SCL_FCTR * scl;

      int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);

      cos_t /= hist_width;

      sin_t /= hist_width;

      int i, j, k, len = (radius*2+1)*(radius*2+1), histlen = (d+2)*(d+2)*(n+2);

      int rows = img.rows, cols = img.cols;

      AutoBuffer buf(len*6 + histlen);

      float *X = buf, *Y = X + len, *Mag = Y, *Ori = Mag + len, *W = Ori + len;

      float *RBin = W + len, *CBin = RBin + len, *hist = CBin + len;

//初始化直方图

      for( i = 0; i < d+2; i++ )

      {

            for( j = 0; j < d+2; j++ )

                  for( k = 0; k < n+2; k++ )

                        hist[(i*(d+2) + j)*(n+2) + k] = 0.;

      }

      for( i = -radius, k = 0; i <= radius; i++ )

            for( j = -radius; j <= radius; j++ )

            {

                

                  float c_rot = j * cos_t - i * sin_t;

                  float r_rot = j * sin_t + i * cos_t;

                  float rbin = r_rot + d/2 - 0.5f;

                  float cbin = c_rot + d/2 - 0.5f;

                  int r = pt.y + i, c = pt.x + j;

                  if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&

                       r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )

                  {

                        float dx = (float)(img.at(r, c+1) - img.at(r, c-1));

                        float dy = (float)(img.at(r-1, c) - img.at(r+1, c));

                        X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;

                        W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;

                        k++;

                  }

            }

      len = k;

      fastAtan2(Y, X, Ori, len, true);

      magnitude(X, Y, Mag, len);

      exp(W, W, len);

      for( k = 0; k < len; k++ )

      {

            float rbin = RBin[k], cbin = CBin[k];

            float obin = (Ori[k] - ori)*bins_per_rad;

            float mag = Mag[k]*W[k];

            int r0 = cvFloor( rbin );

            int c0 = cvFloor( cbin );

            int o0 = cvFloor( obin );

            rbin -= r0;

            cbin -= c0;

            obin -= o0;

            if( o0 < 0 )

                  o0 += n;

            if( o0 >= n )

                  o0 -= n;

            // histogram update using tri-linear interpolation

            float v_r1 = mag*rbin, v_r0 = mag - v_r1;

            float v_rc11 = v_r1*cbin, v_rc10 = v_r1 - v_rc11;

            float v_rc01 = v_r0*cbin, v_rc00 = v_r0 - v_rc01;

            float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 - v_rco111;

            float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 - v_rco101;

            float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 - v_rco011;

            float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 - v_rco001;

            int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;

            hist[idx] += v_rco000;

            hist[idx+1] += v_rco001;

            hist[idx+(n+2)] += v_rco010;

            hist[idx+(n+3)] += v_rco011;

            hist[idx+(d+2)*(n+2)] += v_rco100;

            hist[idx+(d+2)*(n+2)+1] += v_rco101;

            hist[idx+(d+3)*(n+2)] += v_rco110;

            hist[idx+(d+3)*(n+2)+1] += v_rco111;

      }

      // finalize histogram, since the orientation histograms are circular

      for( i = 0; i < d; i++ )

            for( j = 0; j < d; j++ )

            {

                  int idx = ((i+1)*(d+2) + (j+1))*(n+2);

                  hist[idx] += hist[idx+n];

                  hist[idx+1] += hist[idx+n+1];

                  for( k = 0; k < n; k++ )

                        dst[(i*d + j)*n + k] = hist[idx+k];

            }

      // copy histogram to the descriptor,

      // apply hysteresis thresholding

      // and scale the result, so that it can be easily converted

      // to byte array

      float nrm2 = 0;

      len = d*d*n;

      for( k = 0; k < len; k++ )

            nrm2 += dst[k]*dst[k];

      float thr = std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;

      for( i = 0, nrm2 = 0; i < k; i++ )

      {

            float val = std::min(dst[i], thr);

            dst[i] = val;

            nrm2 += val*val;

      }

      nrm2 = SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2), FLT_EPSILON);

      for( k = 0; k < len; k++ )

      {

            dst[k] = saturate_cast(dst[k]*nrm2);

      }

}

static void calcDescriptors(const vector& gpyr, const vector& keypoints,

                                          Mat& descriptors, int nOctaveLayers )

{

//SIFT_DESCR_WIDTH=4，描述直方图的宽度

//SIFT_DESCR_HIST_BINS=8

      int d = SIFT_DESCR_WIDTH, n = SIFT_DESCR_HIST_BINS;

      for( size_t i = 0; i < keypoints.size(); i++ )

      {

            KeyPoint kpt = keypoints[i];

            int octv=kpt.octave & 255, layer=(kpt.octave >> 8) & 255;    //该特征点所在的组序号和层序号

            float scale = 1.f/(1 << octv);         //缩放倍数

            float size=kpt.size*scale;                //该特征点所在组的图像尺寸

            Point2f ptf(kpt.pt.x*scale, kpt.pt.y*scale);    //该特征点在金字塔组中的坐标

            const Mat& img = gpyr[octv*(nOctaveLayers + 3) + layer];    //该点所在的金字塔图像

            calcSIFTDescriptor(img, ptf, kpt.angle, size*0.5f, d, n, descriptors.ptr((int)i));

      }

}

//

//构造函数，用来初始化成员变量

SIFT::SIFT( int _nfeatures, int _nOctaveLayers,

                 double _contrastThreshold, double _edgeThreshold, double _sigma )

      : nfeatures(_nfeatures), nOctaveLayers(_nOctaveLayers),

      contrastThreshold(_contrastThreshold), edgeThreshold(_edgeThreshold), sigma(_sigma)

{

}

int SIFT::descriptorSize() const

{

      return SIFT_DESCR_WIDTH*SIFT_DESCR_WIDTH*SIFT_DESCR_HIST_BINS;    //4*4*8

}

int SIFT::descriptorType() const

{

      return CV_32F;

}

void SIFT::operator()(InputArray _image, InputArray _mask,

                                 vector& keypoints) const

{

      (*this)(_image, _mask, keypoints, noArray());

}

void SIFT::operator()(InputArray _image, InputArray _mask,

                                 vector& keypoints,

                                 OutputArray _descriptors,

                                 bool useProvidedKeypoints) const

{

      Mat image = _image.getMat(), mask = _mask.getMat();

      if( image.empty() || image.depth() != CV_8U )

            CV_Error( CV_StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" );

      if( !mask.empty() && mask.type() != CV_8UC1 )

            CV_Error( CV_StsBadArg, "mask has incorrect type (!=CV_8UC1)" );

//把读取的图像长宽变为2倍，得到初始图像base。

      Mat base = createInitialImage(image, false, (float)sigma);  

//gpyr:高斯金字塔

//dogpyr:dog金字塔

      vector gpyr, dogpyr;

//高斯金字塔的组数

      int nOctaves = cvRound(log( (double)std::min( base.cols, base.rows ) ) / log(2.) - 2);

      //double t, tf = getTickFrequency();

      //t = (double)getTickCount();

      buildGaussianPyramid(base, gpyr, nOctaves);   //构造高斯金字塔，共nOctaves组

      buildDoGPyramid(gpyr, dogpyr);                //根据高斯金字塔构造dog金字塔

      //t = (double)getTickCount() - t;

      //printf("pyramid construction time: %g\n", t*1000./tf);

//useProvidedKeypoints=true表示使用用户提供的特征点，默认为flase;

      if( !useProvidedKeypoints )

      {

            //t = (double)getTickCount();

// //获取尺度空间极值，包含关键点精确定位，分配方向

            findScaleSpaceExtrema(gpyr, dogpyr, keypoints);  

//去除重复的关键点

            KeyPointsFilter::removeDuplicated( keypoints );

// Remove keypoints from some image by mask for pixels of this image.

// filter keypoints by mask

            if( !mask.empty() )

                  KeyPointsFilter::runByPixelsMask( keypoints, mask );

//保留指定数目的关键点

            if( nfeatures > 0 )

                  KeyPointsFilter::retainBest(keypoints, nfeatures);

            //t = (double)getTickCount() - t;

            //printf("keypoint detection time: %g\n", t*1000./tf);

      }

      else

      {

            // filter keypoints by mask

            //KeyPointsFilter::runByPixelsMask( keypoints, mask );

      }

//计算特征描述符

      if( _descriptors.needed() )

      {

            //t = (double)getTickCount();

            int dsize = descriptorSize();

            _descriptors.create((int)keypoints.size(), dsize, CV_32F);

            Mat descriptors = _descriptors.getMat();

            calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers);

            //t = (double)getTickCount() - t;

            //printf("descriptor extraction time: %g\n", t*1000./tf);

      }

}

void SIFT::detectImpl( const Mat& image, vector& keypoints, const Mat& mask) const

{

      (*this)(image, mask, keypoints, noArray());

}

void SIFT::computeImpl( const Mat& image, vector& keypoints, Mat& descriptors) const

{

      (*this)(image, Mat(), keypoints, descriptors, true);

}

}