[人工智能]SIFT python实现以及公式总结

SIFT python实现以及公式总结

算法简介

以下来自百度：
??SIFT由David Lowe在1999年提出，在2004年加以完善 [1-2] 。SIFT在数字图像的特征描述方面当之无愧可称之为最红最火的一种，许多人对SIFT进行了改进，诞生了SIFT的一系列变种。SIFT已经申请了专利。
??SIFT特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。

实验结果对比图

研究sift开始原因是课程作业，在两个星期没日没夜的研究下，翻遍了各大博客，lowe的论文也看了好几遍，终于搞明白了大概过程，特此记录。文中有部分自己实现的源代码，大部分计算过程中的公式我都总结出来了，缺憾是理论部分写的不是很清楚，另一个缺憾是英文水平不够，变量名函数名命名还不够规范。有时间再仔细学理论部分吧，期末周要破大防了orz

实验结果

target和dataset中图片的匹配结果










sift函数调用结果



综合上述对比，target和dataset/3.jpg匹配效果最好，这与实际相符合。对比sift内部函数调用的结果和自己实现的结果，基本吻合，甚至错误匹配比sift更少。对比sift和手写计算出的描述子个数，sift对target生成的描述子个数为1824，匹配个数为427。而自己实现的描述子个数为1829，匹配个数为347，基本持平，侧面说明实现效果较优秀。且对每一张图片生成描述子的时间控制在两分钟以内，不算太慢。

函数实现

SIFT各个函数调用关系

def sift(image):
    print('---------------start sift-----------------')
    image = image.astype('float32')
    image = buildBaseImage(image)
    Octaves = computeOctaves(image.shape)
    GaussPyd = buildGaussPyd(image, Octaves)
    DoGImages = buildDoG(GaussPyd)
    keyPoints = findExtrema(GaussPyd, DoGImages)
    descriptors = buildDescriptors(keyPoints, GaussPyd)
    print('------------------done---------------------')
    return cvt_to_inputImageSize(keyPoints), descriptors

整体关系这张图很形象，取自Alliswell_WP的博客¹

初始化图像

??在建立尺度空间前，Lowe定义原始图像的尺度为0.5，而第0层的尺度为1.6，因此需要先将原始图像扩大一倍，并进行一次高斯模糊使得图像的尺度变成1.6，并能够获取足够多的特征点数量。

??注：想要得到尺度为 ${\sigma }_i$ 的图像，只需要对尺度 ${\sigma }_{i?1}$ 的图像做一次 $\sigma =\sqrt{{\sigma }_i^2?{\sigma }_{i?1}^2}$ 的高斯模糊即可。

代码如下：

def buildBaseImage(img,sigma=1.6):
    img = cv2.resize(img, (0,0), fx=2, fy=2, interpolation = cv2.INTER_LINEAR)
    sigma_diff = sqrt(sigma**2 - 1)
    return cv2.GaussianBlur(img, (0,0), sigmaX=sigma_diff, sigmaY=sigma_diff)

计算高斯金字塔和DoG差分金字塔

??二维图像的尺度空间定义为：
$L(x,y,\sigma )=G(x,y,\sigma )?I(x,y)$

??G为标准高斯函数，$\sigma$的大小决定了图像的精细程度，越小越精细，越大越模糊，模拟了不同尺度下的图像。

??我们需要建立一个共有 Octaves 层，每一层含有 S+3 个图像的金字塔。特征为：

1. 每一层的大小是上一层的一半，且每一层中的所有图像大小相同。
2. 其中金字塔的层数定义为：
   
   $Octaves={\log }_2(min(M,N))?a$
   其中a 一般取1~3，M, N 代表图片的大小。
3. 金字塔中 (O, i) 的图像相当于原图做过一个高斯核为 $\sigma =1.6?2^{O+\frac{i}{S}}$ 的高斯模糊后生成的图像。

这个金字塔即为高斯金字塔。详细公式见图，此处引用Alliswell-WP¹ 博客中的图片

算法思想如下：

??显然，由上述特征以及高斯模糊的叠加性质可以发现，每一层第二张图片开始只需要选择如下数组中的高斯核 [ $k\sigma$ , $k^2\sigma$ , $k^3\sigma$ ,… ] 对第一张图片做高斯模糊即可，而第一张图片为保证高斯空间的连续性，选择上一层中倒数第3张图片直接缩小一倍得到。这样获得的每一张图片根据高斯模糊的叠加性质即符合高斯金字塔的特征。

高斯核数组的生成：

def computeSigma(sigma=1.6, S=3):
    k = 2 ** (1. / S)
    sigma_diff = np.zeros(S+3)
    sigma_diff[0] = sqrt(sigma**2 - 1)
    pre_sigma = sigma 
    for image_index in range(1,S+3):
        cur_sigma = pre_sigma * k
        sigma_diff[image_index] = sqrt(cur_sigma**2 - pre_sigma**2)
        pre_sigma = cur_sigma
    return sigma_diff

按照上述算法思想，算法如下

def buildGaussPyd(image, Octaves, S=3, sigma=1.6):
    GaussPyd = []
    sigma_diff = computeSigma(sigma)
    # 只需要计算第一层的即可，之后每多一层，自然会多乘一次
    for octaves_index in range(Octaves):
        Gauss_images = [image]
        for sigma_kernel in sigma_diff[1:]:
            image = cv2.GaussianBlur(image, (0, 0), sigmaX=sigma_kernel, sigmaY=sigma_kernel)
            Gauss_images.append(image)
        GaussPyd.append(Gauss_images)
        image = Gauss_images[-3]
        image = cv2.resize(image, (int(image.shape[1] / 2),int(image.shape[0] / 2)),interpolation=cv2.INTER_NEAREST)
    return np.array(GaussPyd,dtype=object)

??高斯差分尺度空间定义为：
$D(x,y,\sigma )=L(x,y,k\sigma )?L(x,y,\sigma )$
高斯差分尺度空间由定义即相邻图片作差即可得到

def buildDoG(GaussPyd):
    print("generating DoG...")
    DoGImages = []
    for GaussImages_In_Octaves in GaussPyd:
        DoG_In_Octaves = []
        for i in range(len(GaussImages_In_Octaves)-1):
            DoG_In_Octaves.append(cv2.subtract(GaussImages_In_Octaves[i+1],GaussImages_In_Octaves[i]))
        
        DoGImages.append(DoG_In_Octaves)
    return np.array(DoGImages, dtype=object)

寻找关键点

DoG空间极值点

??寻找关键点的第一步是在DoG金字塔每一组中第二层至倒数第二层中寻找空间极值点，舍弃第一层和倒数第一层。判断是否为空间极值点标准为：将该点和周围8个相邻像素点以及上下两个相邻尺度中2*9个像素点共26个像素点比较，判断是否为极值点。为提高计算的稳定性，舍弃了图像边缘的检索，边缘间隔设定为5.

def isPixelAnExtremum(first_subimage, second_subimage, third_subimage, threshold):
    center_val = second_subimage[1,1]
    if abs(center_val) > threshold:
        # 忽略极值过小的点
        if center_val > 0:
            # 极大值点
            return np.all(center_val >= first_subimage) and \
                    np.all(center_val >= second_subimage) and \
                   np.all(center_val >= third_subimage)
        elif center_val < 0:
            # 极小值点
           return np.all(center_val <= first_subimage) and \
                    np.all(center_val <= second_subimage) and \
                    np.all(center_val <= third_subimage)
    else:
        return False

关键点精确定位与过滤

对于关键点，设image_index为某一层中该图片的编号，我们定义在三维平面上关键点函数 $f(x,y,image\_index)=D(x,y,image\_index)$

??对于DoG空间极值点，需要精确确定关键点的位置和尺度，这里采用牛顿法确定极值点（你说巧不巧，研究的时候正好凸优化学了牛顿法求极值点）。函数取泰勒展开式展开至二阶，如图：




对上述式子求导，令导数等于0，即可得到极值点的偏移量为：




当偏移量在任一维度上的偏移量大于0.5时（即x或y或 image_index），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除。
对于图像像素点中的求导方式，通过如下方式计算：



注： 在后续有关坐标(x, y)的计算中，本文采取了和网上大多数源码一致的方式，即对于numpy内置的矩阵行和列，行作为y，列作为x。

y 
|  1  1   1
|  1  1   P
|  1  1   1
O —— —— —— ——x

如图矩阵，则 P(1, 2) = Mat(2, 1)， 显然x，y和矩阵的行和列下标正好相反！这一点在代码实现中体现较为重要。


牛顿法实现如下

SIFT_MAX_INTERP_STEPS = 5
is_suitable_Extrema = False
atempt = 0
while atempt < SIFT_MAX_INTERP_STEPS:
    first_image, second_image, third_image = DoGImages_In_Octave[image_index-1:image_index+2]
    pixel_cube = np.stack([first_image[i-1:i+2, j-1:j+2],
                        second_image[i-1:i+2, j-1:j+2],
                        third_image[i-1:i+2, j-1:j+2]]).astype('float32') / 255.
    image_shape = first_image.shape
    gradient = computeGradient(pixel_cube)
    Hessian = computeHessian(pixel_cube)
    X_bar = -np.linalg.lstsq(Hessian, gradient, rcond=None)[0]
    if np.all( abs(X_bar) < 0.5 ):
        # 极值点收敛
        is_suitable_Extrema = True
        break
    i += int(round(X_bar[1]))
    j += int(round(X_bar[0]))
    image_index += int(round(X_bar[2]))
    if i < image_border_width or i >= image_shape[0] - image_border_width \
        or j < image_border_width or j >= image_shape[1] - image_border_width \
            or image_index < 1 or image_index > S:
            break
    atempt += 1

对于最终精确定位的极值点，我们还需要去除极值点小于某一个经验值的极值点，即若 $abs(f)<\frac{T}{S}$ ，则去除该极值点。

边缘响应点检测

??除了DoG响应较低的点，还有一些响应较强的点也不是稳定的特征点。DoG对图像中的边缘有较强的响应值，所以落在图像边缘的点也不是稳定的特征点。Lowe通过计算图像的二阶求导矩阵，即Hessain矩阵，得到其较大的特征值 ${\lambda }_{max}$ 和较小的特征值 ${\lambda }_{min}$， 令$\frac{{\lambda }_{max}}{{\lambda }_{min}}=ratio$，若ratio值越大，说明图像在某一个方向梯度值越大， 另一个方向则越小，这恰好是边缘图像的特征；反之，若比例接近于1，则说明，该点附近较为平缓。若想过滤掉边缘点，则可以通过定义一个阈值 r ，当比值大于 r 时过滤掉。为减少计算，Lowe通过如下推导，得到一个判断条件：





若
$\frac{Tr(H{)}^2}{Det(H)}>\frac{(r+1{)}^2}{r}$ ，则判定为边缘点，舍去。Lowe建议取 $r=10$

源代码实现

??综合上述分析，极值点的筛选需要通过四个步骤的筛选和调整，分别为：

1. 判断该像素点是否为周围3x3x3像素空间的极值点
2. 通过牛顿法精确确定该像素点的极值
3. 该极值点归一化后不能小于某一个经验值，默认threhold取0.03，但在opencv中使用了 $\frac{0.04}{S}$
4. 边缘点检测，满足条件 $\frac{Tr(H{)}^2}{Det(H)}<=\frac{(r+1{)}^2}{r}$

当一个极值点通过了上述四个步骤的筛选和调整，则可以认定为一个极值点，在我的实现中，采取了字典来存储该极值点的位置信息(x, y, octave, image_index)，以及该组中的尺度差值 ${\sigma }_{oct}=1.6?2^{\frac{image\_index}{S}}$（参见之前的推导公式），以及实际的尺度 $\sigma ={\sigma }_{oct}?2^{octave}$， 方便后续处理使用。

判断是否为领域极值点

def isPixelAnExtremum(first_subimage, second_subimage, third_subimage, threshold):
    center_val = second_subimage[1,1]
    if abs(center_val) > threshold:
        # 忽略极值过小的点
        if center_val > 0:
            # 极大值点
            return np.all(center_val >= first_subimage) and \
                    np.all(center_val >= second_subimage) and \
                   np.all(center_val >= third_subimage)
        elif center_val < 0:
            # 极小值点
           return np.all(center_val <= first_subimage) and \
                    np.all(center_val <= second_subimage) and \
                    np.all(center_val <= third_subimage)
    else:
        return False

调整极值点 && 过滤过小极值点 && 过滤边缘极值点

def adjustLocalExtrema\
    (DoGImages_In_Octave,i,j,octave_index,image_index,contrast_threshold,sigma,S=3,r=10,image_border_width=5):
    SIFT_MAX_INTERP_STEPS = 5
    is_suitable_Extrema = False
    atempt = 0
    while atempt < SIFT_MAX_INTERP_STEPS:
        first_image, second_image, third_image = DoGImages_In_Octave[image_index-1:image_index+2]
        pixel_cube = np.stack([first_image[i-1:i+2, j-1:j+2],
                            second_image[i-1:i+2, j-1:j+2],
                            third_image[i-1:i+2, j-1:j+2]]).astype('float32') / 255.
        image_shape = first_image.shape
        gradient = computeGradient(pixel_cube)
        Hessian = computeHessian(pixel_cube)
        X_bar = -np.linalg.lstsq(Hessian, gradient, rcond=None)[0]
        if np.all( abs(X_bar) < 0.5 ):
            # 极值点收敛
            is_suitable_Extrema = True
            break
        i += int(round(X_bar[1]))
        j += int(round(X_bar[0]))
        image_index += int(round(X_bar[2]))
        if i < image_border_width or i >= image_shape[0] - image_border_width \
            or j < image_border_width or j >= image_shape[1] - image_border_width \
                or image_index < 1 or image_index > S:
                break
        atempt += 1

    if not is_suitable_Extrema :return None
    adjustExtremum = pixel_cube[1, 1, 1] + 0.5 * np.dot(gradient, X_bar)
    if abs(adjustExtremum)*S >= contrast_threshold:
        Hessian_xy = Hessian[:2,:2]
        Hessian_xy_trace = np.trace(Hessian_xy)
        Hessian_xy_det = np.linalg.det(Hessian_xy)
        if Hessian_xy_det > 0 and r * (Hessian_xy_trace ** 2) < Hessian_xy_det * ((r + 1) ** 2) :
            keyPoint = dict()
            # pt代表原图像中的坐标位置 (x,y) -> (j,i)
            keyPoint['pt'] = ((j+X_bar[0]), (i+X_bar[1]))
            keyPoint['octave_index'] = octave_index
            keyPoint['image_index'] = image_index
            # keyPoint['image_index_bar'] = int(round(X_bar[2]+0.5)*255)
            keyPoint['sigma_oct']  = sigma * (2 ** ((image_index + X_bar[2]) / np.float32(S)))
            keyPoint['scale'] = keyPoint['sigma_oct'] * (2 ** (octave_index + 1))
            keyPoint['functionValueAtUpdatedExtremum'] = abs(adjustExtremum)
            return keyPoint, image_index
    return None

统计所有符合条件的极值点

def findExtrema(GaussPyd, DoGImages, sigma=1.6, S=3, contrast_threshold=0.04):
    print('calc Extrema points...')
    threshold = contrast_threshold / S * 255 # 需要归一化
    keyPoints = []
    image_border_width = 5 # 图像边缘保留
    for octave_index, DoGImages_In_Octave in enumerate(DoGImages):
        for k in range(1, len(DoGImages_In_Octave)-1): # S + 1
            first_image, second_image, third_image = \
             DoGImages_In_Octave[k-1], DoGImages_In_Octave[k], DoGImages_In_Octave[k+1]
            W, H = first_image.shape
            for i in range(image_border_width, W-image_border_width):
                for j in range(image_border_width, H-image_border_width):
                    if isPixelAnExtremum(first_image[i-1:i+2, j-1:j+2], \
                        second_image[i-1:i+2, j-1:j+2], \
                        third_image[i-1:i+2, j-1:j+2], threshold):
                        result = adjustLocalExtrema(DoGImages_In_Octave, i, j, octave_index, k, contrast_threshold, sigma)
                        if result != None:
                            keyPoint, localized_image_index = result
                            keyPoints_with_Oris = computeKeyPointOris(keyPoint, GaussPyd[octave_index, localized_image_index])
                            keyPoints.extend(keyPoints_with_Oris)
    return keyPoints

计算关键点方向

??经过上述过程，我们已经得到了图像中的极值点。接下来，我们需要通过极值点的一系列信息，计算出极值点的主方向Oris，这个方向代表了这个极值点所处的邻域的像素点主要变化方向，反映到图像上则代表着附近颜色变化的主要方向。计算出这个方向后，在后续的特征提取前，若全部旋转到以这个方向为基坐标，则将得到一组满足方向不变性的坐标。假设图像img1, 仅通过旋转得到的图像img2。若将img2中的坐标旋转至img2的主方向上，则理论上来说，将得到和img1完全一致的坐标！
上述分析则是计算关键点方向极为重要的原因之一。

??首先我们可以定义像素点的梯度模值和角度为该像素点箭头方向的长度和角度。计算公式如下：

$$\begin{gathered} m(x,y)=\sqrt{d{x}^2+d{y}^2} \\ theta(x,y)=arctan(\frac{dy}{dx}) \\ dx=L(x+1,y)?L(x?1,y) \\ dy=L(x,y+1)?L(x,y?1) \end{gathered}$$
高斯加权公式如下:
$$w_{i,j}=e^{\frac{?(i^2+j^2)}{2?1.(5\sigma {)}^2}}$$ i，j表示相对于关键点中心的相对坐标。

因此累加到幅度直方图中的值为：
$$m(x,y)?w_{i,j}$$
??Lowe首先采集了关键点所在图像$3\sigma$ 邻域的像素点的所有箭头(Lowe定义$\sigma =1.5{\sigma }_{oct}$)，对所有像素箭头长度分别在360度以10度为间隔，对模值采取高斯分布累加加成，然后取最高值以及高于最高值％80的角度作为主方向和辅助方向，opencv中还要求该方向必须为局部极值点（该角度的左右值都比它小）。是不是给这一连串东西整蒙了，别慌，下面还有详细的步骤
算法步骤及简单解释：

1. 首先建立一个以10度为间隔的，范围为 [0°, 360°] 的直方图。

2. 对于关键点邻域每一个像素点箭头，通过上述公式计算梯度值幅值和角度，计算高斯加权权值，累加到直方图中。

   角度调用numpy的函数 arctan2 而非 arctan 函数，因为第一个函数返回的弧度值范围为(-Pi, Pi)，可以反应到四个象限。注意，这里需要仔细思考角度转化为直方图下标的方式问题，在下一步中类似操作中要保持一致性！否则就会像我一样找错误找了三天都没找明白orz

3. 对直方图进行平滑处理。该方式主要分为两种，第一种是均值平滑，第二种是高斯平滑，都是卷积运算(opencv 源码采取这种方式)。此处不详细赘诉。
   具体参见博客: http://www.mamicode.com/info-detail-447178.html

4. 选择直方图中最大累加值的角度作为主方向，把超过最大角度累加值*peak_ratio的角度作为该关键的辅助方向。

   Lowe的论文指出取peak_ratio=0.8时，大概有15%关键点具有多方向，但这些点对匹配的稳定性至为关键。而对于具有多个方向的关键点，只需要复制一份完全相同的关键点，存储不同的方向即可。
   

5. 对选取的方向进行插值拟合处理，得到更加精确的主方向。差值拟合方法参见: https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
   ??

注：对于该直方图实际上是首尾相接的，因为坐标是连续的，因此在插值的时候，最后一个柱的下一个柱子便是第一个柱。

??在实际实现时，为保持上下函数的一致性，我分别完成了两个工具函数，用于将弧度值转换成梯度直方图下标和插值拟合处理.

def computePeakIndexWithInterpolation(alpha, beta, gamma, peak_index):
    p = 0.5 * (alpha - beta) / (alpha - 2* beta + gamma)
    return peak_index + p

def rad2index(rad, num_bins):
    PI2 = 2 * np.pi
    theta = rad % PI2
    index =  int(floor(num_bins * theta / PI2))
    return index 

计算关键点主方向完整代码如下：

def computeKeyPointOris(keyPoint, Gauss_image, r_factor=3, scl_factor=1.5, width_bins=10, peak_ratio=0.8):
    octave_index = keyPoint['octave_index']
    gauss_shape = Gauss_image.shape 
    # 但在hist中权重w 采用的sigma = 1.5*scale
    sigma_oct = keyPoint['sigma_oct'] 
    r = int(round(r_factor * scl_factor * sigma_oct))
    exp_factor = - 1 / (2 * ((scl_factor * sigma_oct) ** 2))
    num_bins = int(360 // width_bins)
    hist = np.zeros(num_bins)
    smooth_hist = np.zeros(num_bins)
    keyPoints_with_Oris = [] 
    # 计算m theta hist
    for contrast_y in range(-r, r+1):
        y = int(round(keyPoint['pt'][1])) + contrast_y
        if y<= 0 or y >= gauss_shape[0] - 1:continue
        for contrast_x in range(-r, r+1):
            x = int(round(keyPoint['pt'][0])) + contrast_x
            if x<= 0 or x >= gauss_shape[1] - 1:continue
            dx = Gauss_image[y, x+1] - Gauss_image[y, x-1]
            # 注意dy的方向问题！
            dy = Gauss_image[y-1, x] - Gauss_image[y+1, x]
            gradient_m = sqrt(dx*dx + dy*dy)
            # arctan2 可以直接判断是在哪个象限 取值范围为(-Pi, Pi)
            gradient_theta = np.arctan2(dy,dx) # /rad
            # (0, 2 * pi) -> (0, 36)
            hist_index = rad2index(gradient_theta, num_bins)
            hist[hist_index] += gradient_m * np.exp(exp_factor * (contrast_x**2 + contrast_y**2)) 
     
    # 高斯平滑处理
    gauss_weight = np.array([1, 4, 6, 4, 1])
    for i in range(num_bins):
        sub_hist = np.array([hist[i-2], hist[i-1], hist[i], hist[(i+1)%num_bins], hist[(i+2)%num_bins]])
        smooth_hist[i] = np.dot(sub_hist,gauss_weight.T) / 16

    value_max = np.max(smooth_hist)
    for i,value in enumerate(smooth_hist):
        if value_max * peak_ratio <= value and value > smooth_hist[i-1] and value > smooth_hist[(i+1)%num_bins]:
                
                peak_index = computePeakIndexWithInterpolation(smooth_hist[i-1], value, smooth_hist[(i+1)%num_bins], i)
                peak_index %= num_bins
                Oris = 2 * np.pi * peak_index / num_bins
                new_keyPoint = keyPoint.copy()
                new_keyPoint['Oris'] = Oris
                keyPoints_with_Oris.append(new_keyPoint)
    return keyPoints_with_Oris

生成关键点描述子

??经历过上述一系列复杂的操作后，得到了每一个关键点的位置、尺度和主方向信息，在接下来，需要根据关键点生成描述子。首先以关键点为中心，划分成d x d个子区域。(Lowe论文中d取4)每个子区域统计一个以45°为间隔的直方图，直方图高度同样需要高斯分布函数加权累加得到，共8个方向。因此，对于每一个关键点，共有 4*4*8 = 128 个数据，称之为128维向量描述子，为去除光照影响，还需要进行归一化。

??首先，为满足方向不变性，需要将统计的像素点都旋转到以关键点主方向为x轴的坐标系下。虽然读取的像素点是在原本的图像上，但是用于统计直方图的坐标以及角度应该相对于主方向Oris做坐标变换后再记录。坐标和角度旋转到以主方向为x轴的计算公式为：
$$\begin{gathered} x_{rot}=x?cos(Oris)?y?sin(Oris) \\ {y}_{rot}=x?sin(Oris)+y?cos(Oris) \\ thet{a}_{rot}=theta?Oris \end{gathered}$$
其次，我们需要确定以关键点为中心需要遍历的所有像素点范围。以关键点为中心，划分 d x d 个区域，每个区域内有 $m{\sigma }_{oct}$个像素点(m=3, d=4)，考虑到实际计算时需要用到三线性插值法，邻域边长 $a=m{\sigma }_{oct}(d+1)$，再考虑到旋转的因素，最终邻域边长为 $a=m{\sigma }_{oct}(d+1)\sqrt{2}$（像素点个数宜多不宜少）

综合上述两点讨论，描述子的生成步骤如下：

1. 初始化一个直方图hist（window_width， wondow_width， 8）,hist（i，j，theta）表示在子像素区域（i，j）的直方图中方向为theta的柱子。
2. 计算每一个子像素区域的宽度，关键点邻域的大小。
3. 遍历邻域内的每一个像素点，计算旋转至关键点主方向为x轴下的坐标，以及theta角，并计算该像素点处于d x d 中的哪一个子像素区域以及theta位于8个方向中的直方图下标。
4. 和计算关键点主方向时采取的方式类似，计算每个像素点的幅度值乘以高斯分布函数权重，这里高斯核$\sigma =0.5d$
5. 将像素点计算出来的magnitude通过三线性插值累加到对应的子像素区域对应的方向上，即累加到hist(i_bin, j_bin, oris_bin)。

   三线性插值在这里本质上是让每一个处在不同直方图下标中间的像素点按照一定权重同时对多个直方图产生影响。可参见网址：https://blog.csdn.net/Xyc19930930/article/details/77806504

总结公式如下：
$$\begin{gathered} default:m=3,d=4,num\_bins=8 \\ sub\_width=m?{\sigma }_{oct} \\ radius=sub\_width?(d+1)?\sqrt{2}?0.5 \\ {j}_{rot}=j?cos(Oris)?i?sin(Oris) \\ {i}_{rot}=j?sin(Oris)+i?cos(Oris) \\ thet{a}_{rot}=theta?Oris \\ {i}_{bin}=i_{rot}+d?0.5?0.5 \\ {j}_{bin}=j_{rot}+d?0.5?0.5 \\ {w}_{i_{rot},j_{rot}}=\exp (?\frac{(i_{rot}/sub\_width{)}^2+(j_{rot}/sub\_width{)}^2}{2?(0.5d{)}^2}) \\ magnitude=gradient?w_{i_{rot},j_{rot}} \\ thet{a}_{bin}=rad2index(thet{a}_{rot},num\_bins) \end{gathered}$$
三线性插值公式
$$weight=w?d{r}^k?(1?dr{)}^{1?k}?d{c}^m?(1?dc{)}^{1?m}?do?(1?do{)}^{1?n}$$
其中k,m,n为0（像素点的下一个向上取整坐标），或者1（像素点的向下取整坐标）。

描述子门限和归一化

??128维特征向量子生成后，为减少光照变化等影响，需要对他们进行归一化处理。归一化处理后，还需要作描述子门限化。非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值(大于0.2的则就令它等于0.2，小于0.2的则保持不变)。然后再进行一次归一化处理，提高特征的鉴别性。归一化公式为：

$$L_i=\frac{h_i}{\sqrt{{\sum }_{j=1}^{128}{h}_j^2}}$$

思考

实现心得

??高斯金字塔每一层要生成S+3张图像，且每一层第一张照片直接由上一层倒数第三张resize得到，主要原因可能是这样DoG金字塔则为S+2层，而实际寻找关键点时，实际上只在DoG金字塔每一层中间三张图像中寻找，也就说实际用到的只有S层，因此，升采样可能是为了保证寻找极值点用到的S层图像在每一组中满足尺度连续性。

各类函数操作的根本原因

1. 高斯金字塔采用了尺度空间理论，试图在图像领域中模拟人眼观察物体，尺度大则图像细节模糊，尺度小则图像细节特征丰富。
2. 在寻找关键点中，主要目的是寻找图像中十分突出的不会因为光照，旋转，尺度产生较大变化的点，这些点主要包括角点、边缘点、暗区域的亮点以及亮区域的暗点等，这些关键点的提取是后续创建特征描述子以及提高匹配精确度的关键要素。
3. 高斯金字塔使得描述子满足了视角不变性，在主方向上建立的描述子使得描述子满足了旋转不变性，归一化和描述子门限化则使得描述子消除了量纲的差距，减少光照在不同角度对图片特征向量产生的影响。

参考资料

0. Lowe论文（论文里面理论写的很详细，如果要研究论文真的值得一读）
1. 博主Alliswell_WP的sift算法原理详解（他写的超级好，很多其它博客仅介绍原理，且很多公式不清晰甚至是错误的，并且该博客同样给了很多其他参考博客建议）
2. 【OpenCV】SIFT原理与源码分析
   (相信研究sift的人一定看过这个博客)
3. SIFT特征提取分析（128维向量的创建过程很详细）