《Scale Invariant Feature Transform on the Sphere: Theory and Applications》论文阅读和源码理解（一）
摘要
作者提出了一种针对全景图的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
X_0
X0轴相交的N点（北极点）有一个切平面
C
2
C^2
C2，作者定义这个映射为球上任意一点与南极点S的连线与切平面
C
2
C^2
C2相交的点。不妨设球上任意一点为：
ω
≡
(
θ
,
φ
)
\omega \equiv (\theta, \varphi)
ω≡(θ,φ)，那么如图(b)可见映射到的切面的点为：
Φ
(
ω
)
=
2
t
a
n
θ
2
(
c
o
s
(
φ
)
,
s
i
n
(
φ
)
)
\Phi(\omega) = 2tan\frac{\theta}{2}(cos(\varphi), sin(\varphi))
Φ(ω)=2tan2θ(cos(φ),sin(φ))
这是一个双射，因此球面和平面（要注意这里的平面不是一般的平面，接下来会说）可以双向转换。一般我们从摄像头拿到的是一张图像也是平面图像，它的行和列规定了像素的位置。这个平面图像显然和球面映射得到的平面图像在定义上是不一样的，因为映射图像的像素是基于球面的，像素的位置是由 ( θ , φ ) (\theta, \varphi) (θ,φ)决定的。但实际上他们是同一个平面，只是用不同的方式表达，如果我们认为成像平面就是切平面，那么我们将得到，原来用行和列表达的图像，现在行和列的含义变为从球面的纬度和经度表达的图像，进一步可以通过映射到达该像素所在的球面坐标。这样就实现了原始图像到球面图像的转换。作者在代码中就是这样做的。
2. 球面卷积
对于两个定义在球面的函数
f
,
h
∈
L
2
(
S
2
)
f, h \in L^2(S^2)
f,h∈L2(S2)，球面卷积如下：
(
f
∗
h
)
(
ω
)
=
∫
r
∈
S
O
(
3
)
f
(
r
η
)
h
(
r
−
1
ω
)
d
r
(f*h)(\omega) = \int_{r\in SO(3)}f(r\eta)h(r^{-1}\omega)dr
(f∗h)(ω)=∫r∈SO(3)f(rη)h(r−1ω)dr
但是上式很难计算，因此由Driscoll和Healy证明上式可以使用球面傅里叶变换进行高效地计算：
(
f
∗
h
)
^
(
l
,
m
)
=
2
π
4
π
2
l
+
1
f
^
(
l
,
m
)
h
^
(
l
,
0
)
\hat{(f*h)}(l,m) = 2\pi \sqrt{\frac{4\pi}{2l+1}}\hat{f}(l,m) \hat{h}(l,0)
(f∗h)^(l,m)=2π2l+14π
f^(l,m)h^(l,0)
其中， ( ⋅ ) ^ \hat{(\cdot)} (⋅)^表示球面傅里叶变换。
这个时候，如果f是我们的图像，h为高斯滤波器，那么就可以通过上式构造尺度空间。
3. 球面高斯滤波器
Bulow给出了球面的高斯滤波器的傅里叶变换表达式如下：
g S 2 ^ ( l , m , σ ) = 2 l + 1 4 π e − l ( l + 1 ) σ 2 2 \hat{g^{S^2}}(l,m,\sigma) = \sqrt{\frac{2l+1}{4\pi}}e^{\frac{-l(l+1)\sigma^2}{2}} gS2^(l,m,σ)=4π2l+1 e2−l(l+1)σ2
不妨设定义在球面的图像为
I
(
θ
,
φ
)
I(\theta, \varphi)
I(θ,φ)， 那么带入上面的卷积式得到：
(
I
∗
g
)
^
(
l
,
m
)
=
2
π
I
^
(
l
,
m
)
e
−
l
(
l
+
1
)
σ
2
2
\hat{(I*g)}(l,m) = 2\pi \hat{I}(l,m)e^{\frac{-l(l+1)\sigma^2}{2}}
(I∗g)^(l,m)=2πI^(l,m)e2−l(l+1)σ2
为了方便描述，定义Octave的每一层为：
L
S
2
(
θ
,
φ
,
σ
)
=
I
(
θ
,
φ
)
∗
g
S
2
(
θ
,
φ
,
σ
)
L^{S^2}(\theta, \varphi, \sigma) = I(\theta, \varphi) * g^{S^2}(\theta, \varphi, \sigma)
LS2(θ,φ,σ)=I(θ,φ)∗gS2(θ,φ,σ)
于是SIFT中的DoG就可以计算为：
ψ
(
θ
,
φ
,
σ
)
=
L
S
2
(
θ
,
φ
,
k
σ
)
−
L
S
2
(
θ
,
φ
,
σ
)
\psi(\theta, \varphi, \sigma) = L^{S^2}(\theta, \varphi, k\sigma) - L^{S^2}(\theta, \varphi, \sigma)
ψ(θ,φ,σ)=LS2(θ,φ,kσ)−LS2(θ,φ,σ)
回答了上述三个问题后，才正式开始描述如何具体地构造尺度空间。首先，对于已经映射在切平面的输入图像（行为纬度，列为经度），作者主观地认为本身具有半像素的方差
σ
N
=
0.5
∗
π
/
H
\sigma_N = 0.5*\pi/H
σN=0.5∗π/H（H是图像高，这里是为了初始化）。也就是输入图像可以表示为
L
S
2
(
θ
,
φ
,
σ
N
)
L^{S^2}(\theta, \varphi, \sigma_N)
LS2(θ,φ,σN)。然后开始构造和输入图像size相同的第一个octave。假设每个octave含有S层，每层相邻的图像尺度因子
k
=
2
1
/
S
k = 2^{1/S}
k=21/S，底层的初始尺度为
σ
0
\sigma_0
σ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}
现在问题变为，如何从输入的图像，本身设定尺度为
σ
N
\sigma_N
σ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
)
可见，只要合理选择
σ
g
\sigma_g
σg，就可以迭代地求解下一个尺度图像。现在我们有输入图像
L
S
2
(
θ
,
φ
,
σ
N
)
L^{S^2}(\theta, \varphi, \sigma_N)
LS2(θ,φ,σN)，我们选择
σ
g
=
(
σ
0
/
k
)
2
−
σ
N
2
\sigma_g = \sqrt{(\sigma_0/k)^2 - \sigma_N^2}
σg=(σ0/k)2−σN2
，就可以通过卷积得到第一个octave第一层的图像
L
S
2
(
θ
,
φ
,
σ
0
/
k
)
L^{S^2}(\theta, \varphi, \sigma_0/k)
LS2(θ,φ,σ0/k)，继续使用上面的方法就可以把这个octave的所有层都求出来。
进一步，我们需要继续求解下一个octave，按照传统的SIFT，会对这个octave的第一层图像尺寸缩小一半（downsample by 2），然后重复上面的操作。然而作者考虑到采样可能造成的频率混叠问题（aliasing effect），尤其在图像操作都是在频域中进行，引入的混叠在傅里叶反变换回去的时候可能使图像失真。因此在每次下采样时，作者引入一个判据来判断是否需要执行抗混叠操作（anti-aliasing）：
e
−
H
n
(
H
n
+
1
)
(
σ
0
/
k
)
2
8
≤
e
−
1
e^{\frac{-H_n(H_n+1)(\sigma_0/k)^2}{8}} \le e^{-1}
e8−Hn(Hn+1)(σ0/k)2≤e−1
其中，
H
n
H_n
Hn为下采样后图像新的高。作者这里说为了保证在下采样后最大频率处指数项保持较小值。个人理解是保证下采样后带宽的一半要小于采样频率（也就是采样定理），从而减轻频率混叠的影响。
上图是满足抗混叠判据的图像没进行抗混叠操作(a)和进行了抗混叠操作(b)的区别，可见(a)图有明显的频率混叠现象，有很多原本没有的横纹。
到这里就容易理解了，当满足判据，执行抗混叠操作。否则直接下采样，重复之前过程构造octave：
抗混叠操作简单来说就是，先不进行下采样，先增大一倍
σ
\sigma
σ进行滤波，缩窄图像的带宽，然后再进行下采样。这样虽然能消除混叠，但却使图像模糊，也就是一种trade off——用bluring换aliasing的做法。
（二）极值点提取
在得到尺度空间后进一步可以构造出DoG。这个时候就可以进行极值点提取了。这个过程和传统SIFT基本保持一致，只有少量参数有所不同。这个过程可以简单总结为“一个提取两个验证”。
极值的提取简单的说就是，比较当前像素（在DoG上操作）和同层的8邻域和相邻上下两层的9邻域一共26个像素，如果当前像素比这26个像素都大或都小，那就认为是极值点先选上。如图：
紧接着，需要对候选的极值点进行严格的筛选，只有通过“两个验证”的点才是最终的特征点。
首先进行对比度验证（contrast conditions）。由于之前找到的极值点并非是真正的极值点（图像是离散的），因此需要这个极值点周围的邻域进行一个二次拟合，寻找一个真正的极值点。这个过程是迭代进行的。先围绕这个点进行一个泰勒展开，只保留到二次项：
然后对这个式子求导并让其等于0，可求出从该点向极值点的偏移量：
δ
ω
i
^
=
−
(
∂
2
ψ
S
2
∂
Θ
2
)
−
1
∂
ψ
S
2
∂
Θ
\hat{\delta_{\omega_i}} = -(\frac{\partial^2 \psi^{S^2}}{\partial \Theta^2})^{-1} \frac{\partial \psi^{S^2}}{\partial \Theta}
δωi^=−(∂Θ2∂2ψS2)−1∂Θ∂ψS2
如果这个偏移量的任意一维大于0.5，那么就对当前极值点对应的那一维加1，认为新的点才是真正的极值点。然后再由这个点开始，重复上述过程。这里迭代的次数最大为5，
σ
\sigma
σ这个维度是不考虑的。当得到最终的极值点
ω
i
~
\tilde{\omega_i}
ωi~，检验以下条件：
∣
ψ
S
2
(
ω
i
~
)
∣
>
0.02
k
S
2
o
|\psi^{S^2}(\tilde{\omega_i})| > \frac{0.02}{k^S 2^o}
∣ψS2(ωi~)∣>kS2o0.02
其中，o为所在的octave的id。如果符合这个条件，那么这个极值点就可以保留，否则就丢弃。
然后进行主曲率验证（ratio of principle curvature）。这一步是用于挑选特征点，强调这个点是因为，粗提取的极值点不仅仅包含了特征点，还有特征边缘。SIFT中提到，Hessian矩阵的特征值正比于平面上该点的主曲率，如下图：
特征值都比较小则是平坦区域，而两个特征值一大一小且差距较大，则是边缘。两个特征值较大且彼此相近，则是特征点。这个时候，需要判断以下条件：
t
r
a
c
e
(
H
S
2
)
2
d
e
t
(
H
S
2
)
<
(
r
+
1
)
2
r
\frac{trace(H^{S^2})^2}{det(H^{S^2})} < \frac{(r+1)^2}{r}
det(HS2)trace(HS2)2<r(r+1)2
其中H为Hessian矩阵，r取10。只有两个特征值大小相似，左式才有较小值。满足上式则保留为特征点。
（三）LSD描述子
LSD描述子用于全景图之间的匹配，其过程和SIFT描述子是类似的，修改部分只是为了适应球面坐标系。首先，需要计算这个特征点的角度，它是通过该店的梯度直方图来计算的。以该点为中心，考虑一个 3 σ × 3 σ 3\sigma \times 3\sigma 3σ×3σ的窗口，计算这个窗口内所有点的梯度。这里需要定义两个东西：
- 球面上的窗口：这个需要用到球面距离的计算公式，因为在不同纬度， 3 σ 3\sigma 3σ所占列数是不同的。公式如下：
- 球面上的梯度角：作者定义0度为 θ \theta θ增大的方向，90度为 φ \varphi φ增大的方向。
这时每个点的梯度角可以计算为：
α
(
θ
,
φ
,
σ
)
=
a
r
c
t
a
n
(
L
φ
S
2
(
θ
,
φ
,
σ
)
L
θ
S
2
(
θ
,
φ
,
σ
)
)
\alpha(\theta, \varphi, \sigma) = arctan(\frac{L_{\varphi}^{S^2}(\theta, \varphi, \sigma)}{L_{\theta}^{S^2}(\theta, \varphi, \sigma)})
α(θ,φ,σ)=arctan(LθS2(θ,φ,σ)LφS2(θ,φ,σ))
梯度角划分为36个bin，也就是梯度直方图按每10度划分一个格，而每个格的值是由对应该角度的点的梯度值加权累加得到的。这里由两个注意的地方：
- 点的梯度值计算：
m ( θ , φ , σ ) = L φ S 2 ( θ , φ , σ ) + L θ S 2 ( θ , φ , σ ) m(\theta, \varphi, \sigma) = \sqrt{L_{\varphi}^{S^2}(\theta, \varphi, \sigma) + L_{\theta}^{S^2}(\theta, \varphi, \sigma)} m(θ,φ,σ)=LφS2(θ,φ,σ)+LθS2(θ,φ,σ) - 点的权重：由以该点为中心的标准差为 1.5 σ 1.5\sigma 1.5σ高斯函数分配。
得到直方图后，还会以直方图最大值的点为中心拟合一个二次曲线（为了更准确），求出真正的极值角度。对于其他bin大于0.8倍最大值的也考虑为一个角度，为其创建一个新的特征点，赋予这个角度。
当特征点求得其对应的角度时，为了实现旋转不变性，现在开始正式计算具有旋转不变性的LSD描述子。首先，以该点为中心划一个 6 σ × 6 σ 6\sigma \times 6\sigma 6σ×6σ的窗口（其实圆形窗口会更准确）。把这个窗口旋转该点的梯度角，在球面上的旋转可以使用Rodrigues公式。旋转后把窗口内分成4x4的窗格，每个窗格内做点的梯度直方图（方法同上，但对角度只划分8个bin）。也就是说每个窗格可以得到一个8维向量，把4x4个窗格的拼起来就得到一个128维的LSD描述子。最后对这个描述子进行归一化。
上图只分了2x2的窗格，但直方图的bin数也是8个。圆形会更加准确！！
（四）LPD描述子
LPD描述子用于全景图和平面图之间的匹配，由于特征提取一般都是在平面上进行的，LPD描述子使得即使是全景图的场景，已建立好的平面特征数据库依然可用。其实LPD描述子不能说是一种描述子，从代码上更像是一种近似方法，把局部的球面小块近似成平面的方法。我们考虑球面上的点 p i = ( θ i , φ i ) p_i = (\theta_i, \varphi_i) pi=(θi,φi)。上面提到过，传统SIFT在全景图局部是有效的，但需要满足两个条件：1. 小块尺寸较小；2. 小块远离边缘（也就是如果成像面是球面北极点的切面，小块尽量靠近北极点N）。如果以上两个条件满足，容易知道这个球面小块其实和平面小块差别不大，传统的SIFT也能应用到小块上，预先提取的平面SIFT特征也能和这个小块进行匹配。那么，如何使得所有特征点，以其为中心的 6 σ × 6 σ 6\sigma \times 6\sigma 6σ×6σ小块都能有这种最优的近似呢？
作者使用了一个简单而有效的方法，只要把所有特征点都视为北极点即可，那么它对应的南极点就不是以前那个了，而是坐标为 S ′ ( π − θ i , π + φ ) S'(\pi - \theta_i, \pi + \varphi) S′(π−θi,π+φ)。从这个南极点映射到北极点对应的切面。这个时候，我们需要转换特征点周围小块内的点的坐标，因为他们不再是以前的坐标了。
如上图，
p
′
p'
p′为特征点
p
i
p_i
pi周围小块的一个点。回顾之前提到的求球面两点距离Vincenty公式，
θ
\theta
θ可以轻易求得：
θ
=
d
(
p
i
,
p
′
)
\theta = d(p_i, p')
θ=d(pi,p′)
接着需要求
φ
\varphi
φ，其实是需要求
△
φ
\triangle \varphi
△φ，虽然为了理解方便把
p
i
p_i
pi看成是北极点，但实际上在计算的时候仍按照它们的真实坐标。因此现在我们都是求
p
′
p'
p′相对于
p
i
p_i
pi的经纬度增量。在Vincenty公式中，容易发现：
A
2
+
B
2
+
C
2
=
1
⇒
t
a
n
(
d
)
=
s
i
n
(
d
)
c
o
s
(
d
)
=
A
2
+
B
2
C
⇒
c
o
s
(
d
)
=
C
=
c
o
s
(
θ
2
)
c
o
s
(
θ
1
)
+
s
i
n
(
θ
2
)
s
i
n
(
θ
1
)
c
o
s
(
△
φ
)
⇒
△
φ
=
a
c
o
s
(
c
o
s
(
d
)
−
c
o
s
(
θ
2
)
c
o
s
(
θ
1
)
s
i
n
(
θ
2
)
s
i
n
(
θ
1
)
)
\begin{aligned} &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(\triangle \varphi) \\ &\Rightarrow \triangle \varphi = acos(\frac{cos(d) - cos(\theta_2)cos(\theta_1)}{sin(\theta_2)sin(\theta_1)} ) \end{aligned}
A2+B2+C2=1⇒tan(d)=cos(d)sin(d)=CA2+B2
⇒cos(d)=C=cos(θ2)cos(θ1)+sin(θ2)sin(θ1)cos(△φ)⇒△φ=acos(sin(θ2)sin(θ1)cos(d)−cos(θ2)cos(θ1))
进而我们也可以求出 φ \varphi φ，这个时候转换坐标已经求出。
（五）平面到球面映射
LPD描述子可以得到球面和平面匹配的小块图像，这两块图像可以进而使用SIFT描述子进行描述和匹配。现在如果知道平面上一个平面的物体可以直接通过一个线性变换H变换到球面，那就可以轻易地从平面的物体分割转为球面的物体分割。
如上图，作者希望寻找一个3x3的矩阵，把平面的点 p i p l p_i^{pl} pipl转换到三维空间 q i j q_{ij} qij（这两个点是匹配的特征点）。其实，这个H矩阵就是单应矩阵。单应矩阵的作用是可以把一个平面上的点转化到另一个平面。因此，在平面图像中，这个物体必须是共面的。通过这个矩阵把平面上的点转换到一个空间上共面的点，然后这个面再投影到单位球面上，得到球面坐标。
下面求解这个H矩阵。单应矩阵本身有8个*度，也就是需要8个方程才能求解。如上图， p j S 2 p_j^{S^2} pjS2必须和 q i j q_{ij} qij共线，因此需要满足： p j S 2 × H p i ~ p l = 0 p_j^{S^2} \times H\tilde{p_i}^{pl} = 0 pjS2×Hpi~pl=0（ p i ~ p l \tilde{p_i}^{pl} pi~pl是 p i p l p_i^{pl} pipl的齐次坐标形式）。这个方程会给出两个线性无关的方程组，因此需要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