SIFT四部曲之——高斯滤波

本文为原创作品，未经本人同意，禁止转载

欢迎关注我的博客：http://blog.csdn.net/hit2015spring和http://www.cnblogs.com/xujianqing/

或许网络上有各位牛人已经对sift算法进行各种的详解和说明，我（小菜鸟）在翻阅各种资料和对opencv中的代码进行反推之后，终于理解该算法。并记录之，供大家一起交流学习！这个博文主要记录了我的学习历程，或许对你有帮助，或许可以启发你，或许你只是一笑而过！没关系，至少自己总结过。

这篇文章主要是对sift算法的每一个步骤，每一个参数进行说明，并在最后用matlab实现该算法，从理论到代码实现或许需要考虑方方面面，但是它并不是那么的难！

1、算法简介

开始之前要例行一些东西，对sift的简单介绍。如果不想看直接跳到第二部分，没问题！

Sift（Scale-invariant feature transform）尺度不变特征转换，CV界中赫赫有名的算法，由David Lowe（图1的老头）提出，该算法受专利保护，专利权属于英属哥伦比亚大学。

图1 David Lowe

Sift算法可以将一幅图像映射（变换）为一个局部特征向量集；特征向量具有平移、缩放、旋转不变性，同时对光照变化、仿射及投影变换也有一定的不变性。

SIFT算法的特点有：（只是理论，所以看看就好）

1.SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；

2.独特性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；

3.多量性，即使少数的几个物体也可以产生大量的特征向量；

4.高速性，经优化的匹配算法甚至可以达到实时的要求；

5.可扩展性，可以很方便的与其他形式的特征向量进行联合。

SIFT算法可以解决的问题：

目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准目标识别跟踪的性能。而算法在一定程度上可解决：

1.目标的旋转、缩放、平移

2.投影变换

3.光照影响

4.目标遮挡

5.杂物场景

SIFT算法的实质是在不同的尺度空间上查找关键点，并计算出关键点的方向。所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。

SIFT算法可以分解为如下四步：

1.高斯差分（DoG）滤波：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。

2.尺度空间的极值检测和关键点位置确定：对DoG金字塔中的每一层，进行尺度空间的极值检测（极大值和极小值），把每一个极值点作为候选点，在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。

3.关键点方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。

4.构建关键点特征描述符：在每个关键点周围的内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

2.高斯差分（DoG）滤波：

在这个部分开始之前，首先要普及一下一些相关的知识点。

2.1图像尺度空间

在图像信息处理中引入一个名为尺度的参数，通过对图像进行一些变换，获得在多个尺度空间下的图像空间序列，对这些序列可以进行一些特征的提取等操作，可以实现边缘，角点检测和不同分辨率上的特征提取。

对这个尺度空间的理解：它可以模拟人在距离目标由近到远的过程中，目标在视网膜当中形成的图像的过程。尺度越大图像越模糊，相当于我们观察远处物体，这时候关注该物体的轮廓。如下图，我们看到远方只有两个人的外形，并不能看到衣服上的花纹。

尺度空间满足视觉不变性。

• 满足灰度不变性和对比度不变性：当我们用眼睛观察物体时，当物体所处背景的光照条件变化时，视网膜感知图像的亮度水平和对比度是不同的，因此要求尺度空间算子对图像的分析不受图像的灰度水平和对比度变化的影响。
• 满足平移不变性、尺度不变性、欧几里德不变性以及仿射不变性：相对于某一固定坐标系，当观察者和物体之间的相对位置变化时，视网膜所感知的图像的位置、大小、角度和形状是不同的，因此要求尺度空间算子对图像的分析和图像的位置、大小、角度以及仿射变换无关。

一个图像的尺度空间表示在该尺度下，该坐标处的值。

2.2二维高斯函数和高斯模糊

下面这些只是一些基础性的知识，公式并不需要推导，接受就好。

高斯滤波器，使用正态分布计算的一种卷积模板（基本概念，这不懂的话，需要自己入门），利用高斯滤波器和图像进行卷积运算，可对图像进行模糊处理。公式如下（这是一个二维的高斯滤波器）：

二维高斯曲面如下：

图3 高斯二维曲面

其中为正态分布的标准差，在高斯模糊中，它越大，图像越模糊。这里要定义一个模糊半径,

表示模板元素到模板中心的距离。

当然上面只是一个连续的曲面，在对图像进行高斯模糊的过程中需要的是高斯模板，这个模板和图像卷积便可得到高斯模糊图像。正态分布中，在大于3*的范围之外存在的概率占仅0.3%，所以在3*的范围之外，那些像素所起作用基本可以忽略不计了，所以高斯模板只需要计算的大小即可。根据的值计算可以计算出高斯模板。

最后提几条高斯模糊的特性（后面的理解中会被用到）

1. 高斯模糊具有圆对称性，模板是中心对称的
2. 高斯模糊具有线性可分的性质，这样在计算卷积的时候就可以利用一行和一列的两个矩阵和图像进行卷积，可以大大减小计算量。下面的代码就是利用这个性质。
3. 对同一张图片进行连续多次高斯模糊与只用一次大的高斯模糊，可以达到一样的效果。需要满足的是方和根的关系。如两次的模糊值分别为3和4，达到的效果可以只用5就可以。
4. 高斯卷积核被证明是尺度变换唯一的变换核，也是唯一的线性核。

下面是高斯模板的生成代码：

function[g,x] = gaussian_filter(sigma,sample)
sample = 3.5;
if ~exist('sample')
    sample = 7.0/2.0;
end
%设置滤波阈值
n = 2*round(sample * sigma) + 1;

x = 1:n;
x = x-ceil(n/2);
g = exp(-(x.^2)/(2*sigma^2))/(sigma*sqrt(2*pi));
end

 

这里产生的是一行的高斯模板，即运用了对模板进行线性分解的这个性质。当时，该模板为

6.60729028444640e-06	0.000426170838253310	0.00170911194387669	0.000426170838253310	6.60729028444640e-06
0.000426170838253310	0.0274880587288661	0.110237879438303	0.0274880587288661	0.000426170838253310
0.00170911194387669	0.110237879438303	0.442097064144154	0.110237879438303	0.00170911194387669
0.000426170838253310	0.0274880587288661	0.110237879438303	0.0274880587288661	0.000426170838253310
6.60729028444640e-06	0.000426170838253310	0.00170911194387669	0.000426170838253310	0

当时，该高斯模糊的效果如图4

图4高斯模糊效果

2.3高斯金字塔

先说高斯金字塔要得到啥吧：

高斯金字塔主要是为了得到不同尺度的图片，这些图片的尺度必须是连续的，所以要对图片进行高斯滤波。高斯金字塔是一个原始图像，产生几组（octave）每一组中又包含着几层（interval）。如图5：

图5高斯金字塔示意

当然在构建高斯金字塔之前还需要确定的是我们需要构建该金字塔的阶数（o）和每一组的层数（s）。

高斯金字塔的构建主要就是分成两步走

1. 对图像进行不同尺度的高斯模糊（操作上面已经介绍过了）。

   关键：模糊尺度的确定

2. 对高斯金字塔进行降采样

   关键：降采样的母本图片的确定

在高斯金字塔的构建中，图像每一组的大小与相应阶数的对应关系为：（原始图像以512*512为例）

阶数	1	2	3	4	5	6	7	8	9
大小	512	256	128	64	32	16	8	4	2

2.4高斯差分金字塔（这一部分是结论性的知识）

2002年Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其它的特征提取函数（如：梯度，Hessian或Harris角特征）比较，能够产生最稳定的图像特征。

而Lindeberg早在1994年就发现高斯差分函数（Difference of Gaussian 算子）与尺度归一化的高斯拉普拉斯函数非常近似。其中和的关系可以从如下公式推导得到：

利用差分近似代替微分，则有：

其中k-1是个常数，并不影响极值点位置的求取。而高斯拉普拉斯和高斯差分的比较如图6：

图6高斯拉普拉斯和高斯差分

如图所示，红色曲线表示的是高斯差分算子，而蓝色曲线表示的是高斯拉普拉斯算子。使用更高效的高斯差分算子代替拉普拉斯算子进行极值检测，如下：

Ok讲到这里，聪明的读者就应该知道了，学习前面那么多的知识，只是为了对sift特征点的出场做铺垫。在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如图7所示，进行极值检测。这样就能得到sift特征点的候选人，对，只是候选。

图7高斯差分金字塔生成

2.5高斯金字塔生成的细节

上述所有的知识已经把sift关键点候选人的选举办法说清楚了，现在该讲一讲选举过程中所应该注意的4个问题，这样就把这一部分结束。

这两个问题归根到底还是高斯金字塔构建的过程中的4个问题。

1. 金字塔的阶数（o）的确定
2. 每一组层数（s）的确定
3. 每一层的尺度（）的确定
4. 下一组的图片降采样母本的确定

下面一一解答：

1、金字塔的阶数（O）一般为4，也可以根据图像大小来选择，但是要满足下列关系：

分别表示图像的行数和列数。

2、每一阶的层数（S）一般选择5或者6，一般选择6的时候效果最好。在这边就要根据前面的说明，特征点的选举是要在相邻的两层差分金字塔上面进行检测的，所以要得到n个尺度的特征点，就要在层的差分金字塔上检测，（自己画个图就ok了），然而要产生n+2层的差分金字塔，就要有n+3层的高斯金字塔，这样相邻的相减，才能产生n+2层差分。注意：这里的检测都是同阶里面不同层的操作。所以S = n + 3。记住这个n，有用！

3、尺度因子的选择，或许这是本节中最让人头疼的一件事情了。然而在看完众多源码和例程之后我把自己的见解整理如下：

先要理解几个概念：

高斯金字塔的模糊尺度：这个尺度是我们产生模板的尺度

摄像头模糊的尺度：这个尺度是图像被相机镜头模糊后的尺度，一般为固定值，这里定义0.5。

图像的尺度：这个尺度是摄像头模糊尺度和高斯金字塔尺度的合作用，满足方和根的关系。Lowe定义图片的尺度为1.6。

这里插入解答第4个问题，即下一阶的第一层图像是根据上一阶的倒数第三层图像进行降采样得到的。见图8

图8降采样示意图

好了回来继续解答第三个问题，

在一阶的图像内，每一层之间的高斯模糊的尺度因子的比值为 ,

于是同一阶的第s层高斯模糊尺度就变成了,这里面是这一阶的第一层图像的高斯模糊尺度。

根据图7所示，第n+1阶的图像，它的高斯模糊尺度是，这个图像的高斯模糊带到了下一阶的第一层图像中去，于是不同阶相同层的高斯模糊尺度是2倍关系。如图9所示

图9不同层之间的尺度示意

综上可以概括出阶内及阶之间的尺度关系了。

i：表示第i阶的图像

s：表示阶内第s层的图像

在理清楚图像互相之间的关系之后，我们需要的是第一阶，第一层的高斯模糊尺度，这样就能根据上述关系，搞到所有的告示模糊尺度了。还记得两个数吗？0.5和1.6，这时候就派上用场了。

在lowe的论文中，他定义图片的尺度是1.6，被摄像头模糊的尺度是0.5，于是可以算得高斯模糊是

但是这个1.52对我们来说并没有啥意思，lowe为了获得更多的特征点（姑且这么解释吧）他先对原始图像进行扩大一倍。然后呢再对它进行高斯平滑作为第一阶第一层，这时候它的高斯模糊尺度为

于是得到最开始的高斯模糊尺度。这下可以得到所有的高斯模糊尺度了。

到此为止，我们讲完所有理论的工作了，然而在实际实现中，同一阶第s+1层的图片是在第s层图片的基础上进行高斯模糊得到的，还记得那个方和根的关系吧。故，第s+1层图片由第s层图片用一个尺度的模板进行高斯卷积得到。

最后产生高斯金字塔，相邻层相减，便得到高斯差分金字塔。附上这一部分的生成代码：

im = imread('gray_testImage1.jpg');
    im = im2double(im);
 % 设定输入量的默认值
if ~exist('octaves')   %%最大的阶数
    octaves = 4;
 end
 if ~exist('intervals')%%每一阶最大的层数
    intervals = 2;
 end
 if ~exist('object_mask')%%计算模板大小
    object_mask = ones(size(im));
 end
 if size(object_mask) ~= size(im)
    object_mask = ones(size(im));
 end
 if ~exist('contrast_threshold')
    contrast_threshold = 0.02;%%设置去除低对比度特征点阈值大小
 end
 if ~exist('curvature_threshold')
    curvature_threshold = 10.0;%%设置去除边缘特征点阈值大小
 end
 if ~exist('interactive')%%设置迭代次数
    interactive = 1;
 end
interactive = 2;
 
 % 检验输入灰度图像的像素灰度值是否已归一化到[0,1]
 if( (min(im(:)) < 0) | (max(im(:)) > 1) )
    fprintf( 2, 'Warning: image not normalized to [0,1].
' );

 end

 
 % 将输入图像经过高斯平滑处理，采用双线性差值将其扩大一倍.
if interactive >= 1
    fprintf( 2, 'Doubling image size for first octave...
' );
end
 


 tic;
 antialias_sigma = 0.5;           %%高斯平滑参数
 if antialias_sigma == 0         %%不进行平滑操作
    signal = im;
 else
    g = gaussian_filter( antialias_sigma );   %%%产生高斯模板
    if exist('corrsep') == 3
            signal = corrsep( g, g, im );
    else
       signal = conv2( g, g, im, 'same' );  %%高斯卷积
    end
 end
 signal = im;
 [X Y] = meshgrid( 1:0.5:size(signal,2), 1:0.5:size(signal,1) );
 signal = interp2( signal, X, Y, 'linear' );   %%双线性插值扩展图片
 subsample = [0.5]; %  降采样率；

 %下一步是生成高斯和差分高斯(DOG)金字塔，这两个金字塔的数据分别存储在名为gauss_pyr{orient,interval}
 % 和DOG_pyr{orient,interval}的元胞数字中。高斯金字塔含有s+3层，差分高斯金字塔含有s+2层。
 
if interactive >= 1
    fprintf( 2, 'Prebluring image...
' );
end
 

 preblur_sigma = sqrt(sqrt(2)^2 - (2*antialias_sigma)^2);
 if preblur_sigma == 0
    gauss_pyr{1,1} = signal;  %%matlab cell数据类型有待考证
 else
    g = gaussian_filter( preblur_sigma );
    if exist('corrsep') == 3
       gauss_pyr{1,1} = corrsep( g, g, signal );
    else
       gauss_pyr{1,1} = conv2( g, g, signal, 'same' );
    end
 end
 clear signal
 pre_time = toc;
 if interactive >= 1
    fprintf( 2, 'Preprocessing time %.2f seconds.
', pre_time );
 end

 
 initial_sigma = sqrt((2 * antialias_sigma)^2 + preblur_sigma^2);%计算第一层第一阶模糊金字塔的sigma值
 
 %%计算不同阶层的siama
 absolute_sigma = zeros(octaves,intervals + 3);
 absolute_sigma(1,1) = initial_sigma * subsample(1);
 
 %%对生成的金字塔的滤波核大小和标准差进行跟踪
 filter_size = zeros(octaves, intervals+3);
 filter_sigma = zeros(octaves, intervals+3);
 
 tic
 %%计算差分高斯金字塔
 for octave = 1:octaves
     sigma = initial_sigma;
     g = gaussian_filter(sigma);
     filter_size(octave,1) = length(g);
     filter_sigma(octave,1) = sigma;
     DOG_pyr{octave} = zeros(size(gauss_pyr{octave,1},1),size(gauss_pyr{octave,1},2),intervals+2);
     
     %%从第二层计算差分金字塔
     for interval = 2:(intervals + 3)
         sigma_f = sqrt(2^(2/intervals) - 1) * sigma;
         g = gaussian_filter(sigma_f);
         sigma = (2^(1/intervals)) *sigma   %%得到下一个sigma
         absolute_sigma(octave,interval) = sigma * subsample(octave);
         %存储滤波器的核大小及标准差
         filter_size(octave,interval) = length(g);
         filter_sigma(octave,interval) = sigma;
         if exist('corrsep')==3
             gauss_pyr{octave,interval} = corrsep(g,g,gauss_pyr{octave,interval - 1});
         else
             gauss_pyr{octave,interval} = conv2(g,g,gauss_pyr{octave,interval - 1},'same');
         end
         DOG_pyr{octave}(:,:,interval - 1) = gauss_pyr{octave,interval} - gauss_pyr{octave,interval -1};
     end
     if octave < octaves
         sz = size(gauss_pyr{octave,intervals + 1});
         [X,Y]= meshgrid(1:2:sz(2),1:2:sz(1));
         gauss_pyr{octave + 1,1} = interp2(gauss_pyr{octave,intervals + 1},X,Y,'*nearest');
         abslute_sigma(octave+1,1) = absolute_sigma(octave,intervals+1);
         subsample = [subsample subsample(end)*2];
     end
 end
 pyr_time = toc;
 
 
 % 在交互模式下显示高斯金字塔
if interactive >= 2
   sz = zeros(1,2);
   sz(2) = (intervals+3)*size(gauss_pyr{1,1},2);
   for octave = 1:octaves
      sz(1) = sz(1) + size(gauss_pyr{octave,1},1);
   end
   pic = zeros(sz);
   y = 1;
   for octave = 1:octaves
      x = 1;
      sz = size(gauss_pyr{octave,1});
      for interval = 1:(intervals + 3)
                        pic(y:(y+sz(1)-1),x:(x+sz(2)-1)) = gauss_pyr{octave,interval};                         
         x = x + sz(2);
      end
      y = y + sz(1);
   end
   fig = figure;
   clf;
   imshow(pic);
%    resizeImageFig( fig, size(pic), 0.25 );
%    fprintf( 2, 'The gaussian pyramid (0.25 scale).
Press any key to continue.
' );
%    pause;
%    close(fig)
end


% 在交互模式下显示差分高斯金字塔
if interactive >= 2
   sz = zeros(1,2);
   sz(2) = (intervals+2)*size(DOG_pyr{1}(:,:,1),2);
   for octave = 1:octaves
      sz(1) = sz(1) + size(DOG_pyr{octave}(:,:,1),1);
   end
   pic = zeros(sz);
   y = 1;
   for octave = 1:octaves
      x = 1;
      sz = size(DOG_pyr{octave}(:,:,1));
      for interval = 1:(intervals + 2)
                        pic(y:(y+sz(1)-1),x:(x+sz(2)-1)) = DOG_pyr{octave}(:,:,interval);                         
         x = x + sz(2);
      end
      y = y + sz(1);
   end
   fig = figure;
   imshow(pic);
%    clf;
%    imagesc(pic);%%%这个函数可以调整图像显示
%    resizeImageFig( fig, size(pic), 0.25 );
%    fprintf( 2, 'The DOG pyramid (0.25 scale).
Press any key to continue.
' );
%    pause;
%    close(fig)
end

高斯金字塔：

Ok，至此，我们完成第一部分的高斯滤波和金字塔构建工作。

最后，我要贴上几个引导我学习的博客文章链接，对你们表示感谢！尊重作者的劳动成果，知识共享，但不代表着盗取！

未经本人同意禁止转载！

http://blog.csdn.net/zddblog/article/details/7450033