GaussianFilter

https://blog.csdn.net/zxpddfg/article/details/45912561

二维高斯滤波，是图像处理中非常常见的操作。操作的核心是使用一个从高斯分布中采样得到核和原图像做卷积操作。 设高斯核窗口尺寸为\((2\omega+1) \times (2\omega + 1)\),高斯分布的标准差为\(\sigma\).则高斯核可以用矩阵表示：

\[ \mathbf G = \begin{bmatrix} G(-w, -w) &\dots &G(-w, 0) &\dots &G(-w, w)\\ \vdots &{} &\vdots &{} &\vdots\\ G(0, -w) &\dots &G(0, 0) &\dots &G(0, w)\\ \vdots &{} &\vdots &{} &\vdots\\ G(w, -w) &\dots &G(w, 0) &\dots &G(w, w) \end{bmatrix} \]

其中：

\[ G(u,v) = \dfrac{1}{S} \exp\left(-\dfrac{u^2}{2\sigma^2}-\dfrac{v^2}{2\sigma^2}\right) \]

其中\(u\)表示行，\(v\)表示列，\(u, v \in \{-w, -w + 1, \dots, w - 1, w\}\)，\(S\)是归一化常数。

\[ S = \sum_{u = -w}^{w}\sum_{v = -w}^{w} \exp\left(-\dfrac{u^2}{2\sigma^2}-\dfrac{v^2}{2\sigma^2}\right) \]

由于高斯分布概率密度函数非零区间主要分布在\((-3\sigma, 3\sigma)\)，所以一般取\(w \approx 3\sigma\)。

对于图片的第\(i\)行第\(j\)列，使用高斯核计算如下：

\[ Y(i,j) = \sum_{u = -w}^{w} \sum_{v = -w}^{w}X(i+u, j+v)G(u, v) \]

可以看到需要\((2w + 1)\times(2w + 1)\)次乘法以及\((2w + 1)\times(2w + 1) - 1\)次加法，时间复杂度为\(O(w^2)\)。

注意，高斯核的表达式是可分离的。令：

\[ g(x) = \exp\left(-\dfrac{x^2}{2\sigma^2}\right) \]

则：

\[ G(u, v) = \dfrac{1}{S}g(u)g(v) \]

那么高斯核矩阵可以改写为归一化常熟乘以一个行向量乘以一个列向量形式：

\[ \begin{align} \mathbf G &= \dfrac{1}{S} \begin{bmatrix} g(-w)g(-w) &\dots &g(-w)g(0) &\dots &g(-w)g(w) \\ \vdots &\ &\vdots &\ &\vdots \\ g(0)g(-w) &\dots &g(0)g(0) &\dots &g(0)g(w)\\ \vdots &\ &\vdots &\ &\vdots \\ g(w)g(-w) &\dots &g(w)g(0) &\dots &g(w)g(w) \end{bmatrix}\\ &=\dfrac{1}{S} \begin{bmatrix} g(-w)\\ \vdots \\ g(0) \\ \vdots \\ g(w) \end{bmatrix} \times \begin{bmatrix} g(-w) \dots g(0) \dots g(w) \end{bmatrix} \end{align} \]

其中：

\[ \begin{align} S &= \sum_{u = -w}^{w}\sum_{v = -w}^{w}g(u)g(v) =\sum_{u = -w}^{w}\left[\sum_{v = -w}^{w}g(v)\right]g(u) \\ &=\left[\sum_{u = -w}^{w}g(u)\right] \times \left[\sum_{v = -w}^{w}g(v)\right] = S' \times S' \end{align} \]

其中\(S'\)为一维高斯核的归一化系数的倒数，得：

\[ \mathbf G = \mathbf G_1 \times \mathbf G_2 \]

\[ \begin{align} \mathbf G_1 &= \dfrac{1}{S'}\left[g(-w) \dots g(0) \dots g(w)\right]^{\mathrm T}\\ \mathbf G_2 &= \dfrac{1}{S'}\left[g(-w) \dots g(0) \dots g(w)\right] \end{align} \]

可见\(\mathbf G\)分离成两个向量乘积形式。因此对于图片的第\(i\)行第\(j\)列，使用高斯核计算如下：

\[ \begin{align} Y(i,j) &= \sum_{u = -w}^{w} \sum_{v = -w}^{w}X(i+u, j+v)G(u, v) \\ &=\sum_{u = -w}^{w} \sum_{v = -w}^{w}X(i+u, j+v)\dfrac{1}{S} g(u)g(v)\\ &=\sum_{u = -w}^{w} \sum_{v = -w}^{w}X(i+u, j+v)\dfrac{1}{S'}g(u)\dfrac{1}{S'}g(v)\\ &= \sum_{u = -w}^{w}\left[ \sum_{v = -w}^{w}X(i+u, j+v)\dfrac{1}{S'}g(v)\right]\dfrac{1}{S'}g(u)\\ &= \sum_{u = -w}^{w} Z(i + u)\dfrac{1}{S'} g(u) \end{align} \]

上面的式子表明，为了获得最终的高斯滤波的结果，可以先用横向一维高斯核\(\mathbf{G_2}\)与输入图片\(\mathbf{X}\)进行计算，得到中间结果\(\mathbf{Z}\)。再用纵向一维高斯核\(\mathbf{G_1}\)与中间结果\(\mathbf{Z}\)进行计算，得到输出\(\mathbf{Y}\)。需要\((w\omega +2)\)次乘法和\(4\omega\)次加法，时间复杂度为\(O(\omega)\)​.

使用vulkan compute shader实现如下：

#version 450
layout(local_size_x_id = 0, local_size_y_id = 1, local_size_z_id = 2) in;
layout(binding = 0, r32f) uniform readonly image2D input_image;
layout(binding = 1, r32f) uniform image2D output_image;
layout(push_constant) uniform KernelType {
    layout(offset = 0) int size;
    layout(offset = 4) float value[15];
} kernel;

void main(){
    int kerelsize = kernel.size - 1;
    float value = 0., weight = 0.;
    //ivec2 index(gl_GlobalInvocationID.x, gl_GlobalInvocationID.y);
    for (int pos = -kerelsize; pos <= kerelsize; ++pos){
        float w = kernel.value[abs(pos)];
        vec4 pixel = imageLoad(input_image, ivec2(gl_GlobalInvocationID.x + pos, gl_GlobalInvocationID.y));
        value += pixel.r * w; weight += w;
    }
    imageStore(output_image, ivec2(gl_GlobalInvocationID.xy), vec4(value/weight, 0., 0., 0.));
}

尺度空间-高斯金字塔

https://www.cnblogs.com/JiePro/p/sift_1.html

https://www.cnblogs.com/starfire86/p/5735061.html

https://liuchang.men/2020/02/24/SIFT%E7%AE%97%E6%B3%95%E6%B7%B1%E5%85%A5%E7%90%86%E8%A7%A3/

在说高斯金字塔之前，我们先来说一下人的眼睛，我们人眼对世界的感知有两种特性：一是近大远小：同一物体，近处看时感觉比较大，远处看时感觉比较小；二是"模糊"：更准确说应该是"粗细"，我们看近处，可以看到物体的细节(人会觉得比较清楚)，比如一片树叶，近看可以看到该树叶的纹理，远处看只能看到该片的大概轮廓(人会觉得比较模糊). 从频率的角度出发，图像的细节(比如纹理，轮廓等)代表图像的高频成分，图像较平滑区域表示图像的低频成分.

图像高斯金字塔实际上是一种图像的尺度空间(分线性和非线性空间，此处仅讨论线性空间)，尺度的概念用来模拟观察者距离物体的远近程度，在模拟物体远近的同时，还得考虑物体的粗细程序.

综上，图像的尺度空间是模拟人眼看到物体的远近程度以及模糊程度.

其中使用下采样模拟图像的远近程度，采用高斯核对图像进行平滑处理模拟图像的模糊程度。

\[ \begin{align} & O = [log_2 min(M,N)] - 3 \\ & S = n + 3 \\ & k = 2^{1/n} \\ & \sigma(o, r) = \sigma_0 2^{o + \frac{r}{n}} = \sigma_0 2^o {k}^r\\ & o \in [0, 1, \dots, O - 1], r \in [0, 1, \dots, n + 2] \end{align} \]

其中\(M\)为图像行高，\(N\)为图像列宽，\(O\)为高斯金字塔组数，\(n\)为待提取图像特征的图像数（高斯差分求极值的层），\(S\)为图像高斯金字塔每层的层数，\(o\)为组索引号，\(r\)为层索引号，\(\sigma(o,r)\)为对应图像高斯模糊参数。\(\sigma_0\)为高斯模糊初始值，David G.Lowe 教授刚开始设置为1.6，考虑相机实际已对图像进行\(\sigma = 0.5\)的模糊处理,故：

\[ \sigma_0 = \sqrt{1.6^2 - 0.5^2} = 1.52 \]

通过上式即可计算出对应高斯金字塔模糊系数。

假设\(O = 4,S = 6\)，可得：

	第0层	第1层	第2层	第3层	第4层	第5层
第0组	\(\sigma(0,0) \\= \sigma_0 2^{0 + 0/3}\)​​	\(\sigma(0,0) \\= \sigma_0 2^{0 + 1/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{0 + 2/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{0 + 3/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{0 + 4/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{0 + 5/3}\)​
第1组	\(\sigma(0,0) \\= \sigma_0 2^{1 + 0/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{1 + 1/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{1 + 2/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{1 + 3/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{1 + 4/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{1 + 5/3}\)​
第2组	\(\sigma(0,0) \\= \sigma_0 2^{2 + 0/3}\)​​	\(\sigma(0,0) \\= \sigma_0 2^{2 + 1/3}\)​​	\(\sigma(0,0) \\= \sigma_0 2^{2 + 2/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{2 + 3/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{2 + 4/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{2 + 5/3}\)​
第3组	\(\sigma(0,0) \\= \sigma_0 2^{3 + 0/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{3 + 1/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{3 + 2/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{3 + 3/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{4 + 4/3}\)​	\(\sigma(0,0) \\= \sigma_0 2^{4 + 5/3}\)​

在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性，生成高斯金字塔每组有\(S + 3\)层图像。DOG金字塔每组有\(S+2\)层图像。见上图，为\(S=3\)的情况，由上一问题可知，高斯尺度空间中倒数第三幅尺度与下一octave第一幅的尺度相同，由图中红色矩形中的尺度对应为DoG中极值检测的图像，将其各层红色矩形框中的尺度依次排列，即可发现其为以\(k=2^{1/S} = 2^{1/3}\)为等比的连续尺度。所以极值检测是在一个连续变化的尺度空间中进行的。

特征检测

---

References

[1] https://en.wikipedia.org/wiki/Separable_filter [2] https://people.eecs.berkeley.edu/~malik/cs294/lowe-ijcv04.pdf [3] http://www.ipol.im/pub/art/2011/my-asift/ [4] HartSift: A High-Accuracy and Real-Time SIFT based on GPU