0


【OpenCV 例程 300篇】241. 尺度不变特征变换(SIFT)

『youcans 的 OpenCV 例程300篇 - 总目录』

【youcans 的 OpenCV 例程 300篇】241. 尺度不变特征变换(SIFT)

6.4 尺度不变特征变换(SIFT)

6.4.1 简介

尺度不变特征转换算法(Scale-invariant feature transform,SIFT)是图像处理中经典的局部特征描述算法,广泛应用于物体识别、动作识别、影像配准、影像追踪和 3D 建模。

当不同图像尺度系统、方向相似时,可以用角检测和 MSER 作为整体图像特征。但存在尺度变化、方向旋转、光照变化和视点变化时,则需要使用 SIFT 方法提取图像中的不变特征。

SIFT 算法检测与描述影像中的局部性特征,在空间尺度中寻找极值点,并提取出其位置、尺度、旋转不变量。

SIFT 算法的实质是在不同的尺度空间上查找关键点(特征点),计算特征点的大小、方向、尺度信息,利 用这些信息组成关键点对特征点进行描述的问题。SIFT 算法查找的关键点都是高度显著且容易获取的“稳定”特征点,如角点、边缘点、暗区的亮点以及亮区的暗点等,这些特征与大小、旋转无关,对于光线、噪声、视角改变的鲁棒性也很高。

SIFT 算法的特点是:

(1)SIFT 特征不仅具有尺度不变性,而且具有旋转不变性、亮度不变性,对视角变化、仿射变换、噪声也保持一定程度的稳定性;
(2)独特性好,信息量丰富,适用于在海量的特征数据库中进行快速、准确的匹配;
(3)多量性,即使少数几个物体也可以产生大量的 SIFT 特征向量;
(4)高速性,优化的 SIFT 匹配算法可以满足实时要求;
(5)可扩展性,方便与其它形式的特征向量相结合。

SIFT 算法比较复杂,Lowe 将其分解为四个基本步骤:

(1)尺度空间极值检测:通过高斯差分金字塔识别潜在的对于尺度和旋转不变的兴趣点。
(2)关键点的精确定位:通过模型拟合确定关键点位置和尺度。
(3)确定关键点的方向:基于图像局部的梯度方向,确定每个关键点的一个或多个方向。
(4)关键点描述:在关键点的邻域内,在选定的尺度上测量图像局部的梯度。

6.4.2 数学方法

1、构造尺度空间

高斯差分金字塔(DoG pyramid) 是 SIFT 检测的基础。

物体在不同的观测尺度下有不同的表现形态。人在距离目标由近到远时,各尺度图像的模糊程度逐渐变大,尺度越大图像越模糊。图像高斯金字塔是一种图像的尺度空间,模拟人眼看到物体的远近程度以及模糊程度。

图像的尺度空间表达就是图像在所有尺度下的描述。尺度空间表达是由不同高斯核平滑卷积得到,在所有尺度上具有相同的分辨率,而高斯金字塔的多分辨率表达在每层分辨率减少。

高斯核是唯一可以产生多尺度空间的核。一个图像的尺度空间

    L
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    )
   
  
  
   L(x,y,\sigma)
  
 
L(x,y,σ) 定义为原始图像 

 
  
   
    I
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    )
   
  
  
   I(x,y)
  
 
I(x,y) 与一个可变尺度二维高斯函数 

 
  
   
    G
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    )
   
  
  
   G(x,y,\sigma)
  
 
G(x,y,σ) 的卷积运算。

尺度空间:

     L
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     =
    
    
     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     ⋆
    
    
     I
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     )
    
   
   
     L(x, y, \sigma) = G(x,y,\sigma) \star I(x,y) 
   
  
 L(x,y,σ)=G(x,y,σ)⋆I(x,y)

二维空间高斯函数:

     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     =
    
    
     
      1
     
     
      
       2
      
      
       π
      
      
       
        σ
       
       
        2
       
      
     
    
    
     e
    
    
     x
    
    
     
      p
     
     
      
       −
      
      
       (
      
      
       
        x
       
       
        2
       
      
      
       +
      
      
       
        y
       
       
        2
       
      
      
       )
      
      
       /
      
      
       2
      
      
       
        σ
       
       
        2
       
      
     
    
   
   
     G(x,y,\sigma) = \frac{1}{2\pi \sigma ^2} exp^{-(x^2+y^2) /2 \sigma ^2} 
   
  
 G(x,y,σ)=2πσ21​exp−(x2+y2)/2σ2

输入图像

    I
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    )
   
  
  
   I(x,y)
  
 
I(x,y) 依次与标准差为 

 
  
   
    σ
   
   
    ,
   
   
    k
   
   
    σ
   
   
    ,
   
   
    
     k
    
    
     2
    
   
   
    σ
   
   
    .
   
   
    .
   
   
    .
   
  
  
   \sigma,k \sigma,k^2 \sigma...
  
 
σ,kσ,k2σ... 的高斯核卷积,生成一组由常量因子 k 分隔的高斯平滑图像。

SIFT 将尺度空间细分为倍频程(八音度),每个倍频程对应于

    σ
   
  
  
   \sigma
  
 
σ 的加倍。倍频程是指音乐中的频率之比为 2 的两个频段之间的间隔,每个频段的上限频率是下限频率的 2倍,如 31.25、62.5、125、250、500、1K …。

SIFT 将每个倍频程细分为 n 个区间。对于行高为 M、列宽为 N 的图像,可以求得图像高斯金字塔的组数 O、每组的层数 S:

     O
    
    
     =
    
    
     [
    
    
     l
    
    
     o
    
    
     
      g
     
     
      2
     
    
    
     m
    
    
     i
    
    
     n
    
    
     (
    
    
     M
    
    
     ,
    
    
     N
    
    
     )
    
    
     ]
    
    
     −
    
    
     3
    
    
    
     S
    
    
     =
    
    
     n
    
    
     +
    
    
     3
    
   
   
     O = [log_2 min(M,N)] - 3\\ S = n + 3 
   
  
 O=[log2​min(M,N)]−3S=n+3

则:

     σ
    
    
     (
    
    
     o
    
    
     ,
    
    
     r
    
    
     )
    
    
     =
    
    
     
      σ
     
     
      0
     
    
    
     
      2
     
     
      
       o
      
      
       +
      
      
       r
      
      
       /
      
      
       n
      
     
    
    
     , 
    
    
     o
    
    
     ∈
    
    
     [
    
    
     0
    
    
     ,
    
    
     .
    
    
     .
    
    
     .
    
    
     O
    
    
     −
    
    
     1
    
    
     ]
    
    
     ,
    
    
     r
    
    
     ∈
    
    
     [
    
    
     0
    
    
     ,
    
    
     .
    
    
     .
    
    
     .
    
    
     n
    
    
     +
    
    
     2
    
    
     ]
    
   
   
     \sigma (o,r) = \sigma _0 2^{o+r/n},\ o \in [0,... O-1], r \in [0,...n+2] 
   
  
 σ(o,r)=σ0​2o+r/n, o∈[0,...O−1],r∈[0,...n+2]

式中 o 为组 Octave 的索引号,r 为组内各层的索引号,

    σ
   
   
    (
   
   
    o
   
   
    ,
   
   
    r
   
   
    )
   
  
  
   \sigma(o,r)
  
 
σ(o,r) 为各组各层中图像的高斯模糊系数,初值被设为 

 
  
   
    
     σ
    
    
     0
    
   
   
    =
   
   
    1.6
   
  
  
   \sigma _0 = 1.6
  
 
σ0​=1.6 (考虑相机采样带来的模糊,修正为 

 
  
   
    
     σ
    
    
     0
    
   
   
    =
   
   
    
     
      1.
     
     
      
       6
      
      
       2
      
     
     
      −
     
     
      0.
     
     
      
       5
      
      
       2
      
     
    
   
   
    =
   
   
    1.52
   
  
  
   \sigma _0 = \sqrt{1.6^2 - 0.5^2} = 1.52
  
 
σ0​=1.62−0.52​=1.52 )。

2、检测局部极值

关键点是由高斯差分空间的局部极值点组成的。

关键点的初步探查是通过对同一组内各相邻层之间的高斯差分图像进行局部最大值搜索,在空间位置和尺度空间定位局部特征点。为了寻找局部极值点,每一个检测点要与其同尺度的 8 个相邻点及上下相邻尺度对应的 9*2 个点进行比较,以确保检测到尺度空间和图像空间的极值点。

SIFT 检测一个倍频程中两幅相邻尺度空间图像的高斯差的机制,然后与对应于这个倍频程的输入图像进行卷积。在实际计算中,对两个不同高斯尺度空间的图像相减,就得到高斯差分图像

    D
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    )
   
  
  
   D(x,y,\sigma)
  
 
D(x,y,σ):

 
  
   
    
     D
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     =
    
    
     [
    
    
     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     k
    
    
     σ
    
    
     )
    
    
     −
    
    
     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     ]
    
    
     ⋆
    
    
     I
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     )
    
    
    
     =
    
    
     L
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     k
    
    
     σ
    
    
     )
    
    
     −
    
    
     L
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
   
   
     D(x,y,\sigma) = [G(x,y,k\sigma) - G(x,y,\sigma)] \star I(x,y)\\ = L(x,y,k\sigma) - L(x,y,\sigma) 
   
  
 D(x,y,σ)=[G(x,y,kσ)−G(x,y,σ)]⋆I(x,y)=L(x,y,kσ)−L(x,y,σ)

Lindberg 证明尺度空间中的尺度不变性要求使用

     σ
    
    
     2
    
   
  
  
   \sigma ^2
  
 
σ2 对高斯差分算子进行归一化:

 
  
   
    
     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     k
    
    
     σ
    
    
     )
    
    
     −
    
    
     G
    
    
     (
    
    
     x
    
    
     ,
    
    
     y
    
    
     ,
    
    
     σ
    
    
     )
    
    
     ≈
    
    
     (
    
    
     k
    
    
     −
    
    
     1
    
    
     )
    
    
     
      σ
     
     
      2
     
    
    
     
      ▽
     
     
      2
     
    
    
     G
    
   
   
     G(x,y,k\sigma) - G(x,y,\sigma) \approx (k-1) \sigma ^2 \triangledown ^2G 
   
  
 G(x,y,kσ)−G(x,y,σ)≈(k−1)σ2▽2G

在这里插入图片描述

构造尺度空间和检测局部机制的具体实现步骤如下:

(1)对原图进行不同尺度的高斯模糊,得到一组图片,称为一个倍频程(Octave)。每个倍频程内的图片尺寸相同,高斯模糊系数

    σ
   
  
  
   \sigma
  
 
σ 不同;

(2)对第一个倍频程内的图片进行差分操作,得到一组差分数据;
(3)对第一个倍频程中最上面的图片进行高斯降采样,得到尺寸缩小的图片,再对缩小图片进行不同尺度的高斯模糊,得到一组尺寸缩小的图片。
(4)对第二个倍频程内的图片进行差分操作,得到一组差分数据;
(5)重复以上过程,得到各组差分数据。

3、极值点(Key points)的精确定位

高斯差分金字塔在尺度空间和像素点都是离散的,检测到的极值点位置不太准确。高斯差分算子会产生较强的边缘响应,某些极值点对比度较低或稳定性较低。

通过三维二次拟合函数可以求出亚像素位置精度的极值点,精确确定关键点的位置和尺度,并去除低对比度的关键点和不稳定的边缘响应点,以增强匹配稳定性、提高抗噪声能力。

对高斯差分图像

    D
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    )
   
  
  
   D(x,y,\sigma)
  
 
D(x,y,σ) 进行泰勒级数展开:

 
  
   
    
     D
    
    
     (
    
    
     x
    
    
     )
    
    
     =
    
    
     D
    
    
     +
    
    
     (
    
    
     
      
       ∂
      
      
       D
      
     
     
      
       ∂
      
      
       x
      
     
    
    
     
      )
     
     
      T
     
    
    
     x
    
    
     +
    
    
     
      1
     
     
      2
     
    
    
     
      x
     
     
      T
     
    
    
     
      
       
        ∂
       
       
        2
       
      
      
       D
      
     
     
      
       ∂
      
      
       
        x
       
       
        2
       
      
     
    
    
     x
    
   
   
     D(x) = D + (\frac{\partial D}{\partial x})^T x + \frac{1}{2} x^T \frac{\partial ^2 D}{\partial x^2} x 
   
  
 D(x)=D+(∂x∂D​)Tx+21​xT∂x2∂2D​x

对上式求导后,求出极值点的精确位置为:

      x
     
     
      ^
     
    
    
     =
    
    
     −
    
    
     
      H
     
     
      
       −
      
      
       1
      
     
    
    
     (
    
    
     ▽
    
    
     D
    
    
     )
    
    
    
     ▽
    
    
     D
    
    
     =
    
    
     
      
       ∂
      
      
       D
      
     
     
      
       ∂
      
      
       x
      
     
    
    
     =
    
    
     
      [
     
     
      
       
        
         
          
           D
          
          
           x
          
         
        
       
      
      
       
        
         
          
           D
          
          
           y
          
         
        
       
      
      
       
        
         
          
           D
          
          
           σ
          
         
        
       
      
     
     
      ]
     
    
    
    
     
      H
     
     
      
       x
      
      
       y
      
      
       σ
      
     
    
    
     
      [
     
     
      
       
        
         
          
           D
          
          
           
            x
           
           
            x
           
          
         
        
       
       
        
         
          
           D
          
          
           
            x
           
           
            y
           
          
         
        
       
       
        
         
          
           D
          
          
           
            x
           
           
            σ
           
          
         
        
       
      
      
       
        
         
          
           D
          
          
           
            x
           
           
            y
           
          
         
        
       
       
        
         
          
           D
          
          
           
            y
           
           
            y
           
          
         
        
       
       
        
         
          
           D
          
          
           
            y
           
           
            σ
           
          
         
        
       
      
      
       
        
         
          
           D
          
          
           
            σ
           
           
            x
           
          
         
        
       
       
        
         
          
           D
          
          
           
            σ
           
           
            y
           
          
         
        
       
       
        
         
          
           D
          
          
           
            σ
           
           
            σ
           
          
         
        
       
      
     
     
      ]
     
    
   
   
     \hat{x} = -H^{-1}(\triangledown D) \\ \triangledown D = \frac{\partial D}{\partial x} = \begin{bmatrix} D_{x}\\ D_{y}\\ D_{\sigma} \end{bmatrix}\\ H_{xy\sigma}\begin{bmatrix} D_{xx} & D_{xy} & D_{x \sigma}\\ D_{xy} & D_{yy} & D_{y \sigma}\\ D_{\sigma x} & D_{\sigma y} & D_{\sigma \sigma}\\ \end{bmatrix} 
   
  
 x^=−H−1(▽D)▽D=∂x∂D​=⎣⎡​Dx​Dy​Dσ​​⎦⎤​Hxyσ​⎣⎡​Dxx​Dxy​Dσx​​Dxy​Dyy​Dσy​​Dxσ​Dyσ​Dσσ​​⎦⎤​

利用已知的离散空间点插值得到的连续空间极值点的方法,称为子像素插值(Sub-pixel interpolation)。

剔除低对比度的极值点:

SIFT 使用极值位置的函数值

    D
   
   
    (
   
   
    
     x
    
    
     ^
    
   
   
    )
   
  
  
   D(\hat{x})
  
 
D(x^) 来剔除低对比度的不稳定极值。Lowe 提出应剔除 

 
  
   
    ∣
   
   
    D
   
   
    (
   
   
    
     x
    
    
     ^
    
   
   
    )
   
   
    ∣
   
   
    <
   
   
    0.03
   
  
  
   |D(\hat{x})| < 0.03
  
 
∣D(x^)∣<0.03 的极值点。

 
  
   
    
     D
    
    
     (
    
    
     
      x
     
     
      ^
     
    
    
     )
    
    
     =
    
    
     D
    
    
     +
    
    
     
      1
     
     
      2
     
    
    
     (
    
    
     ▽
    
    
     D
    
    
     
      )
     
     
      T
     
    
    
     
      x
     
     
      ^
     
    
   
   
     D(\hat{x}) = D + \frac{1}{2} (\triangledown D)^T \hat{x} 
   
  
 D(x^)=D+21​(▽D)Tx^

消除不稳定的边缘效应:

考虑边缘点与角点的差异,边在边缘方向具有很大的主曲率,而在垂直边缘的方向具有较小的主曲率。

在图像空间计算特征点的二维 Hessian 矩阵:

          H
         
         
          
           x
          
          
           y
          
         
        
        
         =
        
       
      
     
     
      
       
        
        
         
          [
         
         
          
           
            
             
              
               D
              
              
               
                x
               
               
                x
               
              
             
            
           
           
            
             
              
               D
              
              
               
                x
               
               
                y
               
              
             
            
           
          
          
           
            
             
              
               D
              
              
               
                x
               
               
                y
               
              
             
            
           
           
            
             
              
               D
              
              
               
                y
               
               
                y
               
              
             
            
           
          
         
         
          ]
         
        
       
      
     
    
    
     
      
       
        
         T
        
        
         r
        
        
         (
        
        
         
          H
         
         
          
           x
          
          
           y
          
         
        
        
         )
        
       
      
     
     
      
       
        
        
         =
        
        
         
          D
         
         
          
           x
          
          
           x
          
         
        
        
         +
        
        
         
          D
         
         
          
           x
          
          
           y
          
         
        
        
         =
        
        
         α
        
        
         +
        
        
         β
        
       
      
     
    
    
     
      
       
        
         D
        
        
         e
        
        
         t
        
        
         (
        
        
         
          H
         
         
          
           x
          
          
           y
          
         
        
        
         )
        
       
      
     
     
      
       
        
        
         =
        
        
         
          D
         
         
          
           x
          
          
           x
          
         
        
        
         
          D
         
         
          
           y
          
          
           y
          
         
        
        
         −
        
        
         (
        
        
         
          D
         
         
          
           x
          
          
           y
          
         
        
        
         
          )
         
         
          2
         
        
        
         =
        
        
         α
        
        
         β
        
       
      
     
    
   
   
     \begin{aligned} H_{xy} =& \begin{bmatrix} D_{xx} & D_{xy}\\ D_{xy} & D_{yy} \end{bmatrix}\\ Tr(H_{xy}) &= D_{xx} + D_{xy} = \alpha + \beta\\ Det(H_{xy}) &= D_{xx} D_{yy} - (D_{xy})^2 = \alpha \beta \end{aligned} 
   
  
 Hxy​=Tr(Hxy​)Det(Hxy​)​[Dxx​Dxy​​Dxy​Dyy​​]=Dxx​+Dxy​=α+β=Dxx​Dyy​−(Dxy​)2=αβ​

以 r 表示较大特征值

    α
   
  
  
   \alpha
  
 
α 与较小特征值 

 
  
   
    β
   
  
  
   \beta
  
 
β 之比,则比值 

 
  
   
    (
   
   
    r
   
   
    +
   
   
    1
   
   
    
     )
    
    
     2
    
   
   
    /
   
   
    r
   
  
  
   (r+1)^2 /r
  
 
(r+1)2/r 在特征值相等时为最小,随着 r 的增大而增大。

为了检查主曲率是否低于某个阈值 r,只要检查是否满足:

       [
      
      
       T
      
      
       r
      
      
       (
      
      
       
        H
       
       
        
         x
        
        
         y
        
       
      
      
       )
      
      
       
        ]
       
       
        2
       
      
     
     
      
       D
      
      
       e
      
      
       t
      
      
       (
      
      
       
        H
       
       
        
         x
        
        
         y
        
       
      
      
       )
      
     
    
    
     <
    
    
     
      
       (
      
      
       r
      
      
       +
      
      
       1
      
      
       
        )
       
       
        2
       
      
     
     
      r
     
    
   
   
     \frac{[Tr(H_{xy})]^2}{Det(H_{xy})} < \frac{(r+1)^2}{r} 
   
  
 Det(Hxy​)[Tr(Hxy​)]2​<r(r+1)2​

Lowe 提出取 r=10,剔除不满足上式的关键点,可以消除不稳定的边缘效应。

4、关键点的方向

根据图像的局部特征为每个关键点分配一个基准方向,就可以使描述子具有旋转不变性,以便用于图像匹配。

SIFT 使用关键点的邻域像素的梯度分布特性来确定其方向参数,再用梯度直方图求关键点局部结构的稳定方向。

首先用关键点的尺度选择最接近的高斯平滑图像

    L
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    )
   
  
  
   L(x,y,\sigma)
  
 
L(x,y,σ),从而以尺度不变的方式计算所有关键点的方向。对每个关键点,计算以关键点为中心、以 

 
  
   
    (
   
   
    3
   
   
    ∗
   
   
    1.5
   
   
    σ
   
   
    ,
   
   
    3
   
   
    ∗
   
   
    1.5
   
   
    σ
   
   
    )
   
  
  
   (3*1.5\sigma, 3*1.5\sigma)
  
 
(3∗1.5σ,3∗1.5σ) 为半径的邻域图像的幅值 

 
  
   
    M
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    )
   
  
  
   M(x,y)
  
 
M(x,y) 和方向 

 
  
   
    θ
   
   
    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    )
   
  
  
   \theta (x,y)
  
 
θ(x,y):

 
  
   
    
     
      
       
        
         M
        
        
         (
        
        
         x
        
        
         ,
        
        
         y
        
        
         )
        
       
      
     
     
      
       
        
        
         =
        
        
         
          
           [
          
          
           L
          
          
           (
          
          
           x
          
          
           +
          
          
           1
          
          
           ,
          
          
           y
          
          
           )
          
          
           −
          
          
           L
          
          
           (
          
          
           x
          
          
           −
          
          
           1
          
          
           ,
          
          
           y
          
          
           )
          
          
           
            ]
           
           
            2
           
          
          
           +
          
          
           [
          
          
           L
          
          
           (
          
          
           x
          
          
           ,
          
          
           y
          
          
           +
          
          
           1
          
          
           )
          
          
           −
          
          
           L
          
          
           (
          
          
           x
          
          
           ,
          
          
           y
          
          
           −
          
          
           1
          
          
           )
          
          
           
            ]
           
           
            2
           
          
         
        
       
      
     
    
    
     
      
       
        
         θ
        
        
         (
        
        
         x
        
        
         ,
        
        
         y
        
        
         )
        
       
      
     
     
      
       
        
        
         =
        
        
         a
        
        
         r
        
        
         c
        
        
         t
        
        
         a
        
        
         n
        
        
         
          
           L
          
          
           (
          
          
           x
          
          
           ,
          
          
           y
          
          
           +
          
          
           1
          
          
           )
          
          
           −
          
          
           L
          
          
           (
          
          
           x
          
          
           ,
          
          
           y
          
          
           −
          
          
           1
          
          
           )
          
         
         
          
           L
          
          
           (
          
          
           x
          
          
           +
          
          
           1
          
          
           ,
          
          
           y
          
          
           )
          
          
           −
          
          
           L
          
          
           (
          
          
           x
          
          
           −
          
          
           1
          
          
           ,
          
          
           y
          
          
           )
          
         
        
       
      
     
    
   
   
     \begin{aligned} M(x,y) &= \sqrt {[L(x+1,y)-L(x-1,y)]^2 + [L(x,y+1) - L(x,y-1)]^2}\\ \theta (x,y) &= arctan \frac{L(x,y+1) - L(x,y-1)}{L(x+1,y) - L(x-1,y)} \end{aligned} 
   
  
 M(x,y)θ(x,y)​=[L(x+1,y)−L(x−1,y)]2+[L(x,y+1)−L(x,y−1)]2​=arctanL(x+1,y)−L(x−1,y)L(x,y+1)−L(x,y−1)​​

使用直方图统计关键点的邻域像素对应的梯度方向和幅值,将直方图中的最大峰值作为特征点的主方向。梯度直方图的横轴是梯度方向的角度,纵轴是梯度方向对应梯度幅值的累加。梯度方向的范围是 0~360 度,分为 10 个柱(或8个柱)。为了减少突变的影响,可以对直方图进行平滑。

得到特征点的主方向后,对于每个特征点可以得到包括位置、尺度和方向的三个信息

    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    ,
   
   
    θ
   
   
    )
   
  
  
   (x, y, \sigma, \theta)
  
 
(x,y,σ,θ),由此可以确定一个 SIFT 特征区域。SIFT 特征区域由三个值表示,中心表示特征点位置,半径表示特征点的尺度,箭头表示主方向。

5、关键点的描述

通过以上步骤,对于每个特征点可以得到包括位置、尺度和方向的信息

    (
   
   
    x
   
   
    ,
   
   
    y
   
   
    ,
   
   
    σ
   
   
    ,
   
   
    θ
   
   
    )
   
  
  
   (x, y, \sigma, \theta)
  
 
(x,y,σ,θ),由此可以确定一个 SIFT 特征区域。接下来需要用一组向量来描述关键点,构造特征点描述子。特征点描述符不仅包含特征点,还包含特征点周围对其有贡献的像素点。

Lowe 使用在关键点尺度空间内 44 的邻域中计算 8 个方向的梯度,共 44*8=128 维向量作为描述子,具体步骤如下:

(1)确定计算描述子所需的图像区域

梯度方向直方图由关键点所在尺度的高斯图像计算产生。将关键点附近的邻域窗口划分为 d*d 个子区域,每个子区域作为一个种子点, 每个种子点有 8 个方向。考虑到旋转因素,所需窗口的半径 r 为:

     r
    
    
     =
    
    
     [
    
    
     3
    
    
     
      2
     
    
    
      
    
    
     σ
    
    
     (
    
    
     d
    
    
     +
    
    
     1
    
    
     )
    
    
     +
    
    
     1
    
    
     ]
    
    
     /
    
    
     2
    
    
      ,
    
    
     d
    
    
     =
    
    
     4
    
   
   
     r = [3 \sqrt{2} \ \sigma (d+1)+1]/2 \ ,d=4 
   
  
 r=[32​ σ(d+1)+1]/2 ,d=4

(2)将坐标旋转到关键点方向

将坐标轴旋转为关键点的主方向,以实现描述子的旋转不变性。旋转后邻域窗口内采样点的新坐标为:

      (
     
     
      
       
        x
       
       
        ′
       
      
      
       
        y
       
       
        ′
       
      
     
     
      )
     
    
    
     =
    
    
     
      (
     
     
      
       
        
         
          
           c
          
          
           o
          
          
           s
          
          
           θ
          
         
        
       
       
        
         
          
           −
          
          
           s
          
          
           i
          
          
           n
          
          
           θ
          
         
        
       
      
      
       
        
         
          
           s
          
          
           i
          
          
           n
          
          
           θ
          
         
        
       
       
        
         
          
           c
          
          
           o
          
          
           s
          
          
           θ
          
         
        
       
      
     
     
      )
     
    
    
     
      (
     
     
      
       x
      
      
       y
      
     
     
      )
     
    
    
     ,
    
    
    
     x
    
    
     ,
    
    
     y
    
    
     ∈
    
    
     [
    
    
     −
    
    
     r
    
    
     ,
    
    
     r
    
    
     ]
    
   
   
     \binom{x'}{y'}= \begin{pmatrix} cos \theta & -sin \theta\\ sin \theta & cos \theta \end{pmatrix} \binom{x}{y} ,\quad x,y \in [-r,r] 
   
  
 (y′x′​)=(cosθsinθ​−sinθcosθ​)(yx​),x,y∈[−r,r]

在这里插入图片描述

(3)将邻域内的采样点分配到对应的子区域,将子区域内的梯度值分配到 8 个方向,计算其权值。

(4)插值计算每个种子点的 8 个方向的梯度。

(5)统计 44 个子区域、8 个方向的梯度信息,形成 44*8=128 维的向量作为描述该关键点的特征描述向量。

(6)将特征向量进行归一化处理,并设置阈值门限截断较大的梯度值,以消除光照变化的影响。

(7)按特征点的尺度对特征描述向量进行排序。

6、特征点的匹配:

SIFT 算法的基本应用是基于特征提取和描述的结果进行特征匹配。特征点的匹配是通过对两幅图像的特征点集合内的关键点描述子的相似性比对来实现的。

分别对模板图/基准图(Reference image)和观测图/目标图(Observation image)建立关键点描述子集合,采用欧式距离作为 128维的关键点描述向量的相似性度量。

记模板图中关键点描述子为:

     R
    
    
     i
    
   
   
    =
   
   
    (
   
   
    
     r
    
    
     
      i
     
     
      1
     
    
   
   
    ,
   
   
    .
   
   
    .
   
   
    .
   
   
    ,
   
   
    
     r
    
    
     
      i
     
     
      128
     
    
   
   
    )
   
  
  
   R_i = (r_{i1},...,r_{i128})
  
 
Ri​=(ri1​,...,ri128​),观测图图中关键点描述子为:

 
  
   
    
     S
    
    
     i
    
   
   
    =
   
   
    (
   
   
    
     s
    
    
     
      i
     
     
      1
     
    
   
   
    ,
   
   
    .
   
   
    .
   
   
    .
   
   
    ,
   
   
    
     s
    
    
     
      i
     
     
      128
     
    
   
   
    )
   
  
  
   S_i = (s_{i1},...,s_{i128})
  
 
Si​=(si1​,...,si128​),任意两个描述子的相似性度量值为:

 
  
   
    
     d
    
    
     (
    
    
     
      R
     
     
      i
     
    
    
     ,
    
    
     
      S
     
     
      i
     
    
    
     )
    
    
     =
    
    
     
      
       
        ∑
       
       
        
         j
        
        
         =
        
        
         1
        
       
       
        128
       
      
      
       (
      
      
       
        r
       
       
        
         i
        
        
         j
        
       
      
      
       −
      
      
       
        s
       
       
        
         i
        
        
         j
        
       
      
      
       
        )
       
       
        2
       
      
     
    
   
   
     d(R_i,S_i) = \sqrt{\sum_{j=1}^{128}(r_{ij}-s_{ij})^2} 
   
  
 d(Ri​,Si​)=j=1∑128​(rij​−sij​)2​

    (
   
   
    
     R
    
    
     i
    
   
   
    ,
   
   
    
     S
    
    
     i
    
   
   
    )
   
  
  
   (R_i,S_i)
  
 
(Ri​,Si​) 满足阈值条件时,则认为 

 
  
   
    (
   
   
    
     R
    
    
     i
    
   
   
    ,
   
   
    
     S
    
    
     i
    
   
   
    )
   
  
  
   (R_i,S_i)
  
 
(Ri​,Si​) 是配对的关键点描述子。

关键点的匹配可以采用穷举法实现,但一般采用二叉树算法来搜索以提高算法效率。以目标图像的关键点为基准,在模板图像中搜索与最邻近的和次邻近的特征点。

Lowe 提出,比较最近邻距离与次近邻距离之比 ratio,小于某个阈值的认为是正确匹配 ,推荐 0.4~0.6。

运用最近邻点比次近邻点进行限制,可大大提高特征点匹配的正确率

6.4.3 实现算法

6.4.4 OpenCV 中的 SIFT 类

OpenCV 提供了丰富的特征检测算法,而且继承了 cv::Feature2D 类,采用了统一的定义和封装。如:AffineFeature、AgastFeatureDetector、AKAZE、BRISK、FastFeatureDetector、GFTTDetector、KAZE、MSER、ORB、SimpleBlobDetector、SIFT、SURF 等。
OpenCV 中提供 cv::SIFT 类实现 SIFT方法。cv::SIFT 类继承了 cv::Feature2D 类,通过 create 静态方法创建。

SIFT类封装了 SIFT 提取关键点和计算描述符的参数,类的声明在 include/opencv2/ features2d .hpp 文件中。

详细内容请参见: 【OpenCV 说明文档 (https://docs.opencv.org/5.x/d7/d60/classcv_1_1SIFT.html)】

1、构造 SIFT 对象:

在 Python 语言中,OpenCV 提供了 SIFT 类的接口函数 cv.SIFT.create() 实例化 SIFT 类。

cv.SIFT.create([, nfeatures[, nOctaveLayers[, contrastThreshold[, edgeThreshold[, sigma]]]]]) → retval
cv.SIFT.create(nfeatures, nOctaveLayers, contrastThreshold, edgeThreshold, sigma, descriptorType) → retval
cv.SIFT_create([, nfeatures[, nOctaveLayers[, contrastThreshold[, edgeThreshold[, sigma]]]]]) → retval
cv.SIFT_create(nfeatures, nOctaveLayers, contrastThreshold, edgeThreshold, sigma, descriptorType) → retval

参数说明:

  • nfeatures:保留的最佳特征的数量,按算法测度值大小排序
  • nOctaveLayers:每个倍频程的层数,根据图像分辨率自动计算得到
  • contrastThreshold:对比度阈值,用于滤除低对比度的弱特征,默认值 0.04
  • edgeThreshold:用于过滤边缘特征的阈值,默认值 10
  • sigma:高斯模糊系数 σ \sigma σ 的初值,默认值 1.6
  • descriptorType:描述符的类型,仅支持 CV_32F 和 CV_8U

注意事项:

  • 函数 cv.SIFT.create() 实例化 SIFT 类,构造一个 SIFT 对象。
  • 描述符的类型 descriptorType 仅支持 CV_32F 和 CV_8U
  • 函数 cv.SIFT.getDefaultName 返回字符串的标识符。当对象保存到文件或字符串时,该字符串用作顶级 xml/yml 节点的标记。

2、检测特征点:

sift.detect(image,[, mask]) → keypoints 

参数说明:

  • sift:实例化的 SIFT 对象
  • image:输入图像,单通道
  • mask:掩模图像,8 位单通道,指定查找关键点的区域,可选项
  • keypoints :检测到的关键点,是一个特殊的数据结构,包括以下属性: - Point2f pt:坐标- float size:关键点的邻域直径- float angle:关键点的方向,取值为 [0,360)- float response:关键点的响应强度- int octave:关键点所在的图像金字塔的组- int class_id:关键点的 id 编号

3、绘制特征点:

cv.drawKeypoint(image, keypoints, outImage[, color[, flags]]) → outImage 

参数说明:

  • image:输入图像
  • keypoints:输入图像,单通道
  • outimage:输出图像,
  • color:绘制关键点的颜色
  • flags:绘制关键点特征的选择 - cv.DRAW_MATCHES_FLAGS_DEFAULT,默认值,对每个关键点仅绘制中心点,不包括围绕关键点的圆和大小、方向- cv.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS,对于每个关键点绘制表示关键点大小和方向的圆- cv.DRAW_MATCHES_FLAGS_DRAW_OVER_OUTIMG,在输入图像上绘制匹配关系- cv.DRAW_MATCHES_FLAGS_NOT_DRAW_SINGLE_POINTS,不绘制单个关键点

4、其它成员函数:

cv.Feature2D.compute(images, keypoints[, descriptors]) → keypoints, descriptors
cv.Feature2D.detect(image[, mask]) → keypoints
cv.Feature2D.detectAndCompute(image, mask[, descriptors[, useProvidedKeypoints]]) → keypoints, descriptors
cv.Feature2D.defaultNorm() → retval
cv.Feature2D.descriptorType() → retval
cv.Feature2D.empty() → retval
cv.Feature2D.getDefaultName() → retval
cv.Feature2D.read(fileName) → None
cv.Feature2D.write(fileName) → None

本节只列出相关函数供参考,具体使用将在后续内容中介绍。

关于专利问题的说明:

由于 SIFT 和 SURF 都是专有和专利算法,在学术研究中可以免费使用的,但在商业应用中需要获得授权许可。

因此,OpenCV3、OpenCV4 在默认安装中删除了 SIFT 和 SURF 等专利算法,并将其转移到 “OpenCV_contrib” 包中。为了使用 “OpenCV_contrib” 包的算法,在 OpenCV/C++ 中,需要启用 OpenCV contrib 支持,从源代码编译和安装 OpenCV。在 OpenCV/Python 中,需要通过 pip 安装 OpenCV contrib python 包。

但在 OpenCV3/Python 的一些版本中,由于 OpenCV 用

OPENCV_ENABLE_NONFREE

对编译进行了限制,使用 SIFT 和 SURF 算法时仍然会出现错误:

"module ‘cv2.cv2’ has no attribute ‘xfeatures2d’ ”

或者:

cv2.error: OpenCV(3.4.8) C:\projects\opencv-python\opencv_contrib\modules\xfeatures2d\src\sift.cpp: 1207: error: (-213:The function/feature is not implemented) This algorithm is patented and is excluded in this configuration; Set OPENCV_ENABLE_NONFREE CMake option and rebuild the library in function ‘cv::xfeatures2d::SIFT::create’

卸载之前的 OpenCV 和包,将 OpenCV 版本退到 3.4.2 即可解决。

2020年,由于 SIFT 专利到期, OpenCV 4.4 已经将 SIFT 移到主库,可以正常使用。

例程 14.22:特征检测之 SIFT 角点检测

# 14.22 特征检测之尺度不变特征变换 (SIFT)
    img = cv2.imread("../images/imgRoof1.png", flags=1)# (400, 600, 3)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)# 比较算法:Harris 检测角点
    dst = cv2.cornerHarris(gray,2,3, k=0.04)
    corners = np.column_stack(np.where(dst>0.1*dst.max()))# 筛选并返回角点坐标 (y,x)
    corners = corners.astype(np.int)# (128, 2), 检测到的角点的点集 (y,x)
    imgHarris = np.copy(img)for point in corners:# 注意坐标次序
        cv2.circle(imgHarris,(point[1], point[0]),4,(0,0,255),1)# 在点 (x,y) 处画圆# SIFT 关键点检测# sift = cv2.xfeatures2d.SIFT_create()
    sift = cv2.SIFT.create()# sift实例化对象
    kp = sift.detect(gray,None)# 关键点检测,kp 为关键点信息(包括方向)# kp, des = sift.compute(gray, None)  # 关键点检测,des 为特征描述向量 (923, 128)
    img1 = cv2.drawKeypoints(img, kp,None)# 只绘制关键点位置
    img2 = cv2.drawKeypoints(img, kp,None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)# 关键点大小和方向

    plt.figure(figsize=(9,6))
    plt.subplot(221), plt.axis('off'), plt.title("Origin")
    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    plt.subplot(222), plt.axis('off'), plt.title("Harris corners")
    plt.imshow(cv2.cvtColor(imgHarris, cv2.COLOR_BGR2RGB))
    plt.subplot(223), plt.axis('off'), plt.title("SIFT keypoints")
    plt.imshow(cv2.cvtColor(img1, cv2.COLOR_BGR2RGB))
    plt.subplot(224), plt.axis('off'), plt.title("SIFT keypoints")
    plt.imshow(cv2.cvtColor(img2, cv2.COLOR_BGR2RGB))
    plt.tight_layout()
    plt.show()

在这里插入图片描述

【本节完】

版权声明:
youcans@xupt 原创作品,转载必须标注原文链接:(https://blog.csdn.net/youcans/article/details/125828053)
Copyright 2022 youcans, XUPT
Crated:2022-7-25

  1. OpenCV 中的 Harris 角点检测
  2. Harris 角点检测之精确定位(cornerSubPix)
  3. OpenCV 中的 Shi-Tomas 角点检测
  4. 尺度不变特征变换(SIFT)

本文转载自: https://blog.csdn.net/youcans/article/details/126018654
版权归原作者 YouCans 所有, 如有侵权,请联系我们删除。

“【OpenCV 例程 300篇】241. 尺度不变特征变换(SIFT)”的评论:

还没有评论