本文共 10849 字，大约阅读时间需要 36 分钟。

摘要

       作者提出了一种针对全景图的SIFT算法。传统的SIFT算法是基于平面图的，而作者对其进行改进，基本流程和传统SIFT是一样的，但作者把其中平面处理的部分转化为在球面坐标上进行处理。在次基础上提出两种描述子，LSD用于两张全景图的匹配，而LPD用于全景图和平面图的匹配。

主要贡献

• 应用球面尺度空间来替换传统SIFT的平面尺度空间，在下采样过程中考虑到了频率混叠效应（aliasing effect）。
• 提出了LSD（local spherical descrpitor）用于全景图间的匹配。
• 提出了LPD（local plannar descriptor）用于全景图和平面图间的匹配。
• 提出一种在匹配的LPD间的双向映射，用于直接从平面图的一个特征点映射到全景图的对应特征点，反之亦然。（注意并不是全部的映射，而是一对一的）

动机（Motivation： Why SIFT in Spherical Coordinates”）

        为什么使用SIFT，因此它是对尺度和旋转都具有不变性的算法，而把它应用在全景图像中却不如它应用在平面上那样效果良好。

        多数改良的工作集中在，1.把全景图映射到圆柱坐标系，2.把SIFT应用到展开的全景图，3.应用到畸变矫正的全景图。当全景图在成像边缘存在大量信息时，上述方法都会遇到困难；也就是说上述方法仅在局部是有效。而使用球面来描述全景图是自然的，而且不会引入额外的畸变。

原理

        原理介绍的流程基本上按照SIFT的流程来，因此读者阅读之前需要对SIFT有大致的了解。

（一）尺度空间构造

        为了保证尺度不变性，SIFT引入了尺度空间的概念，说白了就是：既然不同的尺度的图像会影响特征点的提取，那么就遍历尽量多的尺度就好了，在不同的尺度都提取一下特征点。如下图：

SIFT构造了不同图片大小的金字塔层（下称octave），每个octave含有多个相同size但经过不同的高斯滤波器（$\sigma$ 不同）处理得到的图片。因此，尺度空间是考虑了两个因素，一个是图像的size，一个是filter的参数。

        而为了把上述过程迁移到球面，需要解决以下问题：

• 如何存储球面和操作球面
• 如何求球面卷积
• 如何求球面高斯滤波器

1. 球面坐标

        首先，定义在球面的数据难以使用和进一步操作，如果能转化成矩阵的形式才方便计算机运算。作者是这样定义球面和平面的映射的（实际上我认为这样定义也是符合成像原理的）：

(a)图中有一个单位球和一个坐标系，球与$X_0$轴相交的N点（北极点）有一个切平面$C^2$，作者定义这个映射为球上任意一点与南极点S的连线与切平面$C^2$相交的点。不妨设球上任意一点为：$\omega \equiv (\theta ,\phi )$，那么如图(b)可见映射到的切面的点为：

$$\Phi (\omega )=2tan\frac{\theta }{2}(cos(\phi ),sin(\phi ))$$

这是一个双射，因此球面和平面（要注意这里的平面不是一般的平面，接下来会说）可以双向转换。一般我们从摄像头拿到的是一张图像也是平面图像，它的行和列规定了像素的位置。这个平面图像显然和球面映射得到的平面图像在定义上是不一样的，因为映射图像的像素是基于球面的，像素的位置是由$(\theta ,\phi )$决定的。但实际上他们是同一个平面，只是用不同的方式表达，如果我们认为成像平面就是切平面，那么我们将得到，原来用行和列表达的图像，现在行和列的含义变为从球面的纬度和经度表达的图像，进一步可以通过映射到达该像素所在的球面坐标。这样就实现了原始图像到球面图像的转换。作者在代码中就是这样做的。

2. 球面卷积

        对于两个定义在球面的函数$f,h\in L^2(S^2)$，球面卷积如下：

$$(f\ast h)(\omega )={\int }_{r\in SO(3)}f(r\eta )h(r^{-1}\omega )dr$$

但是上式很难计算，因此由Driscoll和Healy证明上式可以使用球面傅里叶变换进行高效地计算：

$$\widehat{(f\ast h)}(l,m)=2\pi \sqrt{\frac{4\pi }{2l+1}}\hat{f}(l,m)\hat{h}(l,0)$$

其中，$\widehat{(\cdot )}$表示球面傅里叶变换。

        这个时候，如果f是我们的图像，h为高斯滤波器，那么就可以通过上式构造尺度空间。

3. 球面高斯滤波器

        Bulow给出了球面的高斯滤波器的傅里叶变换表达式如下：

$$\widehat{g^{S^2}}(l,m,\sigma )=\sqrt{\frac{2l+1}{4\pi }}{e}^{\frac{-l(l+1){\sigma }^2}{2}}$$

不妨设定义在球面的图像为$I(\theta ,\phi )$， 那么带入上面的卷积式得到：

$$\widehat{(I\ast g)}(l,m)=2\pi \hat{I}(l,m)e^{\frac{-l(l+1){\sigma }^2}{2}}$$

为了方便描述，定义Octave的每一层为：

$$L^{S^2}(\theta ,\phi ,\sigma )=I(\theta ,\phi )\ast g^{S^2}(\theta ,\phi ,\sigma )$$ 于是SIFT中的DoG就可以计算为： $$\psi (\theta ,\phi ,\sigma )=L^{S^2}(\theta ,\phi ,k\sigma )-L^{S^2}(\theta ,\phi ,\sigma )$$

---

        回答了上述三个问题后，才正式开始描述如何具体地构造尺度空间。首先，对于已经映射在切平面的输入图像（行为纬度，列为经度），作者主观地认为本身具有半像素的方差${\sigma }_N=0.5\ast \pi /H$（H是图像高，这里是为了初始化）。也就是输入图像可以表示为$L^{S^2}(\theta ,\phi ,{\sigma }_N)$。然后开始构造和输入图像size相同的第一个octave。假设每个octave含有S层，每层相邻的图像尺度因子$k=2^{1/S}$，底层的初始尺度为${\sigma }_0$。于是，作者希望每一层的尺度是这样的：

{ σ 0 / k , σ 0 , k σ 0 , . . . , k S − 2 σ 0 } \{\sigma_0/k, \sigma_0, k\sigma_0, ..., k^{S-2}\sigma_0\} {

σ0​/k,σ0​,kσ0​,...,kS−2σ0​}

        现在问题变为，如何从输入的图像，本身设定尺度为${\sigma }_N$得到以下各种尺度呢？我们回顾上面的公式可以发现，空间域的卷积等于频域的乘积在这里也是适用的：

L S 2 ( θ , φ , σ i ) = L S 2 ( θ , φ , σ i − 1 ) ∗ g S 2 ( θ , φ , σ g ) = F − 1 { I ^ ( l , m ) e − l ( l + 1 ) σ i − 1 2 2 ⋅ e − l ( l + 1 ) σ g 2 2 } = F − 1 { I ^ ( l , m ) e − l ( l + 1 ) ( σ i − 1 2 + σ g 2 ) 2 } = L S 2 ( θ , φ , σ i − 1 2 + σ g 2 ) \begin{aligned} L^{S^2}(\theta, \varphi, \sigma_i) &= L^{S^2}(\theta, \varphi, \sigma_{i-1}) * g^{S^2}(\theta, \varphi, \sigma_g) \\ &=F^{-1}\{ \hat{I}(l,m)e^{\frac{-l(l+1)\sigma_{i-1}^2}{2}} \cdot e^{\frac{-l(l+1)\sigma_g^2}{2}}\} \\ &= F^{-1}\{\hat{I}(l,m)e^{\frac{-l(l+1)(\sigma_{i-1}^2+ \sigma_g^2)}{2}} \} \\ &= L^{S^2}(\theta, \varphi, \sqrt{\sigma_{i-1}^2+ \sigma_g^2}) \end{aligned} LS2(θ,φ,σi​)​=LS2(θ,φ,σi−1​)∗gS2(θ,φ,σg​)=F−1{

I^(l,m)e2−l(l+1)σi−12​​⋅e2−l(l+1)σg2​​}=F−1{

I^(l,m)e2−l(l+1)(σi−12​+σg2​)​}=LS2(θ,φ,σi−12​+σg2​ ​)​ 可见，只要合理选择${\sigma }_g$，就可以迭代地求解下一个尺度图像。现在我们有输入图像$L^{S^2}(\theta ,\phi ,{\sigma }_N)$，我们选择${\sigma }_g=\sqrt{({\sigma }_0/k{)}^2-{\sigma }_N^2}$，就可以通过卷积得到第一个octave第一层的图像$L^{S^2}(\theta ,\phi ,{\sigma }_0/k)$，继续使用上面的方法就可以把这个octave的所有层都求出来。

        进一步，我们需要继续求解下一个octave，按照传统的SIFT，会对这个octave的第一层图像尺寸缩小一半（downsample by 2），然后重复上面的操作。然而作者考虑到采样可能造成的频率混叠问题（aliasing effect），尤其在图像操作都是在频域中进行，引入的混叠在傅里叶反变换回去的时候可能使图像失真。因此在每次下采样时，作者引入一个判据来判断是否需要执行抗混叠操作（anti-aliasing）：

$$e^{\frac{-H_n(H_n+1)({\sigma }_0/k{)}^2}{8}}\le e^{-1}$$ 其中，$H_n$为下采样后图像新的高。作者这里说为了保证在下采样后最大频率处指数项保持较小值。个人理解是保证下采样后带宽的一半要小于采样频率（也就是采样定理），从而减轻频率混叠的影响。

上图是满足抗混叠判据的图像没进行抗混叠操作(a)和进行了抗混叠操作(b)的区别，可见(a)图有明显的频率混叠现象，有很多原本没有的横纹。

        到这里就容易理解了，当满足判据，执行抗混叠操作。否则直接下采样，重复之前过程构造octave：

抗混叠操作简单来说就是，先不进行下采样，先增大一倍$\sigma$进行滤波，缩窄图像的带宽，然后再进行下采样。这样虽然能消除混叠，但却使图像模糊，也就是一种trade off——用bluring换aliasing的做法。

（二）极值点提取

        在得到尺度空间后进一步可以构造出DoG。这个时候就可以进行极值点提取了。这个过程和传统SIFT基本保持一致，只有少量参数有所不同。这个过程可以简单总结为“一个提取两个验证”。

        极值的提取简单的说就是，比较当前像素（在DoG上操作）和同层的8邻域和相邻上下两层的9邻域一共26个像素，如果当前像素比这26个像素都大或都小，那就认为是极值点先选上。如图：

        紧接着，需要对候选的极值点进行严格的筛选，只有通过“两个验证”的点才是最终的特征点。

        首先进行对比度验证（contrast conditions）。由于之前找到的极值点并非是真正的极值点（图像是离散的），因此需要这个极值点周围的邻域进行一个二次拟合，寻找一个真正的极值点。这个过程是迭代进行的。先围绕这个点进行一个泰勒展开，只保留到二次项：

然后对这个式子求导并让其等于0，可求出从该点向极值点的偏移量：

$$\widehat{{\delta }_{{\omega }_i}}=-(\frac{{\partial }^2{\psi }^{S^2}}{\partial {\Theta }^2}{)}^{-1}\frac{\partial {\psi }^{S^2}}{\partial \Theta }$$

如果这个偏移量的任意一维大于0.5，那么就对当前极值点对应的那一维加1，认为新的点才是真正的极值点。然后再由这个点开始，重复上述过程。这里迭代的次数最大为5，$\sigma$这个维度是不考虑的。当得到最终的极值点$\widetilde{{\omega }_i}$，检验以下条件：

$$|{\psi }^{S^2}(\widetilde{{\omega }_i})|>\frac{0.02}{k^S{2}^o}$$ 其中，o为所在的octave的id。如果符合这个条件，那么这个极值点就可以保留，否则就丢弃。

        然后进行主曲率验证（ratio of principle curvature）。这一步是用于挑选特征点，强调这个点是因为，粗提取的极值点不仅仅包含了特征点，还有特征边缘。SIFT中提到，Hessian矩阵的特征值正比于平面上该点的主曲率，如下图：

特征值都比较小则是平坦区域，而两个特征值一大一小且差距较大，则是边缘。两个特征值较大且彼此相近，则是特征点。这个时候，需要判断以下条件： $$\frac{trace(H^{S^2}{)}^2}{det(H^{S^2})}<\frac{(r+1{)}^2}{r}$$ 其中H为Hessian矩阵，r取10。只有两个特征值大小相似，左式才有较小值。满足上式则保留为特征点。

（三）LSD描述子

        LSD描述子用于全景图之间的匹配，其过程和SIFT描述子是类似的，修改部分只是为了适应球面坐标系。首先，需要计算这个特征点的角度，它是通过该店的梯度直方图来计算的。以该点为中心，考虑一个$3\sigma \times 3\sigma$的窗口，计算这个窗口内所有点的梯度。这里需要定义两个东西：

1. 球面上的窗口：这个需要用到球面距离的计算公式，因为在不同纬度，$3\sigma$所占列数是不同的。公式如下：

2. 球面上的梯度角：作者定义0度为$\theta$增大的方向，90度为$\phi$增大的方向。

这时每个点的梯度角可以计算为：

$$\alpha (\theta ,\phi ,\sigma )=arctan(\frac{L_{\phi }^{S^2}(\theta ,\phi ,\sigma )}{L_{\theta }^{S^2}(\theta ,\phi ,\sigma )})$$

梯度角划分为36个bin，也就是梯度直方图按每10度划分一个格，而每个格的值是由对应该角度的点的梯度值加权累加得到的。这里由两个注意的地方：

1. 点的梯度值计算：
   $$m(\theta ,\phi ,\sigma )=\sqrt{L_{\phi }^{S^2}(\theta ,\phi ,\sigma )+L_{\theta }^{S^2}(\theta ,\phi ,\sigma )}$$
2. 点的权重：由以该点为中心的标准差为$1.5\sigma$高斯函数分配。

得到直方图后，还会以直方图最大值的点为中心拟合一个二次曲线（为了更准确），求出真正的极值角度。对于其他bin大于0.8倍最大值的也考虑为一个角度，为其创建一个新的特征点，赋予这个角度。

        当特征点求得其对应的角度时，为了实现旋转不变性，现在开始正式计算具有旋转不变性的LSD描述子。首先，以该点为中心划一个$6\sigma \times 6\sigma$的窗口（其实圆形窗口会更准确）。把这个窗口旋转该点的梯度角，在球面上的旋转可以使用Rodrigues公式。旋转后把窗口内分成4x4的窗格，每个窗格内做点的梯度直方图（方法同上，但对角度只划分8个bin）。也就是说每个窗格可以得到一个8维向量，把4x4个窗格的拼起来就得到一个128维的LSD描述子。最后对这个描述子进行归一化。

上图只分了2x2的窗格，但直方图的bin数也是8个。圆形会更加准确！！

（四）LPD描述子

        LPD描述子用于全景图和平面图之间的匹配，由于特征提取一般都是在平面上进行的，LPD描述子使得即使是全景图的场景，已建立好的平面特征数据库依然可用。其实LPD描述子不能说是一种描述子，从代码上更像是一种近似方法，把局部的球面小块近似成平面的方法。我们考虑球面上的点$p_i=({\theta }_i,{\phi }_i)$。上面提到过，传统SIFT在全景图局部是有效的，但需要满足两个条件：1. 小块尺寸较小；2. 小块远离边缘（也就是如果成像面是球面北极点的切面，小块尽量靠近北极点N）。如果以上两个条件满足，容易知道这个球面小块其实和平面小块差别不大，传统的SIFT也能应用到小块上，预先提取的平面SIFT特征也能和这个小块进行匹配。那么，如何使得所有特征点，以其为中心的$6\sigma \times 6\sigma$小块都能有这种最优的近似呢？

        作者使用了一个简单而有效的方法，只要把所有特征点都视为北极点即可，那么它对应的南极点就不是以前那个了，而是坐标为$S^{\prime }(\pi -{\theta }_i,\pi +\phi )$。从这个南极点映射到北极点对应的切面。这个时候，我们需要转换特征点周围小块内的点的坐标，因为他们不再是以前的坐标了。

        如上图，$p^{\prime }$为特征点$p_i$周围小块的一个点。回顾之前提到的求球面两点距离Vincenty公式，$\theta$可以轻易求得：

$$\theta =d(p_i,p^{\prime })$$ 接着需要求$\phi$，其实是需要求$△\phi$，虽然为了理解方便把$p_i$看成是北极点，但实际上在计算的时候仍按照它们的真实坐标。因此现在我们都是求$p^{\prime }$相对于$p_i$的经纬度增量。在Vincenty公式中，容易发现： $$\begin{array}{rl}& A^2+B^2+C^2=1\\ & \Rightarrow tan(d)=\frac{sin(d)}{cos(d)}=\frac{\sqrt{A^2+B^2}}{C}\\ & \Rightarrow cos(d)=C=cos({\theta }_2)cos({\theta }_1)+sin({\theta }_2)sin({\theta }_1)cos(△\phi )\\ & \Rightarrow △\phi =acos(\frac{cos(d)-cos({\theta }_2)cos({\theta }_1)}{sin({\theta }_2)sin({\theta }_1)})\end{array}$$

进而我们也可以求出$\phi$，这个时候转换坐标已经求出。

（五）平面到球面映射

        LPD描述子可以得到球面和平面匹配的小块图像，这两块图像可以进而使用SIFT描述子进行描述和匹配。现在如果知道平面上一个平面的物体可以直接通过一个线性变换H变换到球面，那就可以轻易地从平面的物体分割转为球面的物体分割。

        如上图，作者希望寻找一个3x3的矩阵，把平面的点$p_i^{pl}$转换到三维空间$q_{ij}$（这两个点是匹配的特征点）。其实，这个H矩阵就是单应矩阵。单应矩阵的作用是可以把一个平面上的点转化到另一个平面。因此，在平面图像中，这个物体必须是共面的。通过这个矩阵把平面上的点转换到一个空间上共面的点，然后这个面再投影到单位球面上，得到球面坐标。

        下面求解这个H矩阵。单应矩阵本身有8个自由度，也就是需要8个方程才能求解。如上图，$p_j^{S^2}$必须和$q_{ij}$共线，因此需要满足：$p_j^{S^2}\times H{\widetilde{p_i}}^{pl}=0$（${\widetilde{p_i}}^{pl}$是$p_i^{pl}$的齐次坐标形式）。这个方程会给出两个线性无关的方程组，因此需要4组不共线的匹配特征点才能求解这个H矩阵。由于匹配的特征点数目可能远远超过4组，因此可以使用RANSAC迭代求解最佳的单应矩阵H。

---

参考文献

[1] Cruz-Mota J , Bogdanova I , Paquier B , et al. Scale Invariant Feature Transform on the Sphere: Theory and Applications[J]. International Journal of Computer Vision, 2012, 98(2):217-241.

[2] Lowe D G . Distinctive Image Features from Scale-Invariant Keypoints[J]. International Journal of Computer Vision, 2004, 60(2):91-110.

[3] https://towardsdatascience.com/sift-scale-invariant-feature-transform-c7233dc60f37

[4] https://zhuanlan.zhihu.com/p/74597564

转载地址：http://moklf.baihongyu.com/