（五）用MSER算法提取特征区域
 分水岭算法，通过逐步水淹创建分水岭，把图像分割成多个区域；MSER（最大稳定外部区域）同样通过提高水位的方法，关注水淹过程中的某段时间内，保持相对稳定的盆地，可发现这些区域对应图像中某些问题的特殊部分。
 5.1 实现
 计算图像MSER的基础类是cv::MSER，它是一个抽象接口，继承自cv::Feature2D类。OpenCV中的所有特征检测类都是从这个类继承的，cv::MSER类的实例可以通过create方法创建。在初始化时指定被检测区域的最小和最大尺寸，以便限制被检测特征的数量，调用方式如下：
 //基本的MSER检测器cv::Ptrcv::MSERptrMSER=cv::MSER::create(5,200,2000);
 参数：局部检测时使用的增量值、允许的最小面积、允许的最大面积
 现在可以通过调用detectRegions方法来获得MSER，指定输入图像和一个相关的输出数据结构，代码如下：
 //点集的容器
 std::vector<std::vectorcv::Point> points;
 //矩形的容器 std::vectorcv::Rect rects;
 //检测MSER特性
 ptrMSER->detectRegions(image,points,rects);
 检测结果放在两个容器中。第一个是区域容器，每个区域用组成它的像素点表示；第二个是矩形容器，每个矩形包围一个区域。为了呈现结果，创建一个空白图像，在图像上用不同的颜色显示检测到的区域（颜色是随机选择的）。实现代码如下：
 
       //创建白色图像
       cv::Mat  output(image,size(),CV_8UC3);
       output= cv::Scalar(255,255,255);
       //OpenCV随机数生成器
       cv::RNG  rng; 
       //针对每个检测到的特征区域，在彩色区域显示MSER
       //反向排序，先显示较大的MSER
       for(std::vector<std::vector<cv::Point>>::reverse_iterator    it=points.rbegin(); it!=points.rend();++it)
       {    
            //生成随机颜色    
            cv::Vec3b  c(rng.uniform(0,254),
                                 rng.uniform(0,254),
                                 rng.uniform(0,254));    
            //针对MSER集合中的每个点    
            for(std::vector<cv::Point>::iterator 
            itPts=it->begin();           itPts=it->end();  ++itPts)   
             {      
                  //不重写MSER的像素      
                  If (output.at<cv::Vec3b>(*itPts)[0]==255)     
                   {         
                       output.at<cv::Vec3b>(*itPts)=c;
                    }
             }
      }
      注意，MSER会形成层叠区域。为了显示全部区域，如果一个较大区域内包含了较小的区域，就不能覆盖它。从结果图中显示可以检测出MSER点集，没有全部区域，但是显示出了有意义区域（如建筑物的窗户等）。 
    5.2 原理
    MSER和分水岭的算法原理相同，即高度为0~255，逐渐淹没图像。图像处理技术中，把高于某个阈值的像素集合成为高度集。
    MSER：随着水位上升，颜色较黑和边界陡峭的区域会形成盆地，且在一段时间内有相对稳定的形状（用水位表示颜色，水位的高低代表像素值的强度），这些稳定的盆地即MSER。
    检测稳定性的方法：计算区域的当前面积以及该区域原先的面积（比当前区域低一个特定值的时候），并比较这两个面积；如果相对变化达到局部最小值，则认识该区域是MSER。 
    注意：
    a 增量值将作为cv::MSER类构造函数的第一个参数，用以测量相对稳定性，默认值为5；     
     b 区域面积必须在预定义的范围内，构造函数中后面两个参数就是允许的最小和最大区域尺寸。      
     c 必须确保MSER是稳定的（第四个参数），即形状的相对变化必须足够小。      
     d 一个稳定区域可以属于另一个更大的区域（称为父区域）；为了确保有效性，一个父MSER和它的子区域必须有足够大的差别，即差异限度，由cv::MSER类构造函数的第五个参数指定。      
     e 最后两个参数（即四、五）都使用了默认值，MSER允许的最大相对变化的默认值为0.25，父MSER与子区域的最小差别的默认值为0.2。要检测MSER，必须对参数规范化，否则难以应对不同环境。
     MSER检测器：
     a 首先输出一个包含像素集的容器，每个像素集构成一个区域；因为我们需要找整个区域的位置，而不是里面的单个像素，所以通常用包含了被检测区域的几何形状表示一个MSER。
     b 检测过程中，输出的第二项是一系列矩形，画出所有矩形就能表示检测的结果。为避免区域之间互相包含，需要提取出所有包含垂直矩形的区域，实现方法就是讲每个矩形的面积与检测到的对应区域进行比较，如果两者一致（此处判断标准是两者比例超过0.6），即为一个区域；代码如下：
 
//提取并显示矩形的MSER
 std::vectorcv::Rect::iterator itr=rects.begin(); std::vector<std::vectorcv::Point>::iterator itp=points.begin();
 for( ; itr!=rects.end(); ++itr, ++itp)
 {
 //检查两者比例
 if (static_cast(itp->size()/itr->area()>0.6)
 {
 cv::rectangle(image, *itr, cv::Saclar(255),2);
 }
 }
 提取到的MSER可用图示。在其他应用程序中，也可采用别的判断标准和显示方法，例如判断依据是检测到的区域不能太细长（将封闭的矩形旋转，计算其宽高比），然后用未旋转的封闭椭圆表示他们，代码如下：
 
     //提取并显示椭圆形的MSER 
     for( std::vector<std::vector<cv::Point>>::iterator  
               it=points.begin();  it!=points.end(); ++it)      
      {          
             //遍历MSER集合中的每个点           
             for (std::vector<cv::Point>::iterator  
                 itPts=it->begin(); itPts!=it->end(); ++itPts)             
                 {                
                        //提取封闭的矩形                 
                        cv::RotateRect  rr=cv::minAreaRect(*it);                
                        //检查椭圆的长宽比               
                         if    (rr.size.height/rr.size.height>0.6||rr.size.height/rr.size.height>1.6)
                         {
                           cv::ellipse(image, rr, cv::Scalar(255), 2);
                         }
                  }
      }
      结果图显示有很多椭圆。注意，对于有父子关系的MSER，表示它们的椭圆通常比较相似。在某些情况下，可以施加一个约束条件，要求椭圆之间的差距不低于某个特定值，以免显示重复的椭圆。

 
 
PNAS | 土壤氮循环微生物功能特征的全球生物地理学
 

 
1 摘要

 
微生物驱动了地球氮循环的大部分，但我们仍然缺乏对执行土壤氮过程的微生物丰度和组成的全球概览。为了解决这一研究差距，我们利用公开可用的土壤宏基因组数据，对微生物氮功能特征的生物地理学进行了描述，定义为八种氮循环途径。生境类型、土壤碳、氮在很大程度上解释了样本的总氮途径相对丰度。相比之下，我们不能确定的驱动氮循环功能基因的物种分类组成。此外，含有某种氮代谢途径的优势属在不同生境类型中的分布普遍相似。在土壤生境发现了一个意料之外的高相对丰度的由细菌携带的途径，异化异硝酸盐还原，其是土壤中一个很少研究的氮循环过程。最后，系统发育分析表明，一些微生物类群似乎是氮循环的specialists（专一）或generalists（广泛）。例如，Deltproteobacteria具有所有的8 种氮循环途径，而蓝细菌主要编码3种氮循环途径。总的来说，这种基于性状的方法为研究全球土壤中微生物多样性和氮循环之间的关系提供了基础。
 
 
 
关键词 硝化作用|固氮|氨同化|宏基因组|异化亚硝酸盐还原
 
 
 
2 研究背景
 
本世纪的一个重大挑战是预测环境变化将如何改变全球生物地球化学循环。生物地理学在这一努力中扮演着重要的角色。环境变化正在改变生物多样性的分布，而生物多样性的分布又是生物地球化学过程的关键驱动因素。从历史上看，生物地理学通过分类学的视角来看待生物多样性，主要解决物种分布问题。然而，关注特性——尤其是那些涉及生态系统过程的特性——可能提供生物多样性模式和生物地球化学之间更清晰的联系。
 
这些概念尤其与微生物有关。微生物催化主要元素的生物地球化学循环,因为他们的绝对丰度占有大量的元素物质。此外,像植物和动物一样,微生物分类组成在空间的变化可以影响生态系统过程。因此，考虑到微生物的特性，应该努力将生物地理模式和生态系统过程联系起来。
 
在这里，我们首次描述了土壤中微生物氮循环特征的全球生物地理模式。微生物通过与大气交换(固氮和反硝化作用)和硝态氮淋失来调节生物有效氮。它们还影响植物吸收可利用氮的形式。与此同时，人类活动通过增加生物圈中活性氮，已经并将继续改变氮循环。在局部尺度上，氮的添加会不断改变土壤和其他生态系统中的微生物组成。因此，微生物性状的分布与当前和未来的氮循环有关。
 

 
图1. 收集宏基因组数据的全球地理位置

 
3 主要结果
 
利用宏基因组学分析服务器(MG-RAST)检索表层土壤的宏基因组数据。在对样本进行序列和数据质量的整理之后，最终获得了365个样本，其中代表了来自10种不同栖息地类型（Habitat type）的118个独特地理位置，覆盖了自然和人为主导的系统(图1）。不同样本的测序深度差异很大，但没有明显偏向于任何特定的栖息地类型(图S1)。为了标准化测序深度，我们分析每N个通路在样本中的相对丰度。无论将通路相对丰度归一化为每个注释序列或每个标记基因检测到的数量(基于30个保守的单克隆基因)，观察到的趋势都是相似的(图S1)。细菌在宏基因组文库中占主导地位，占全部序列的95%，其次是真菌占3%，古菌仅占2%。已知真菌序列在宏基因组文库中的比例低于它们对土壤微生物生物量的贡献(10)。因此，我们集中分析细菌和古菌，只分析真菌的一般趋势。例如，细菌、古菌和真菌总序列的比例在不同生境类型中存在差异(g独立性检验;P<0.001)(图S2)。按栖息地划分，古生菌占所有序列的0.9-11%，在沙漠中发现的比例最高。真菌与细菌序列的比率在温带森林土壤中特别高，正如之前观察到的(28)。
 
3.1 土壤氮代谢途径的相对丰度。
 
在土壤样本中，有0.5%的注释序列与8个N通路中的其中一个相关(图2A)，其中细菌和古菌的每个标记基因平均分别与3.3和4.7个N通路相关。单个通路的相对丰度以几个数量级变化(单因素方差分析P<0.001;F=74.21,df=7)(图2B)。细菌和古菌在N代谢的相对丰度上表现出相似的趋势，除了在古菌中没有亚硝酸盐向铵的异化还原途径。真菌序列只与同化途径有关，包括氨同化、硝酸盐同化、亚硝酸盐同化和亚硝酸盐同化氨。在不同域（domain)中，最常见的途径是氨同化(图2B)。例如，在细菌中，每100万条注释的细菌序列平均检测到280条氨同化途径。相比之下，硝化作用和N固定是最不常见的途径，每百万序列分别检测到6.1次和4.6次。值得注意的是，相对未研究的亚硝酸还原氨途径比这两种途径略常见，平均每百万序列检测9.3次。在所有的土壤样品中，细菌和古菌的N代谢相对丰度绝大多数都是正相关的(图3 A和B)。在细菌和古菌中，反硝化作用与硝酸盐对亚硝酸盐异化还原高度正相关(R2=0.860.97,P≤0.001)。这种关系是可以预料的，因为硝酸还原成亚硝酸盐是完全反硝化过程的第一步;然而，这里我们将这两个步骤分开，因为硝酸还原为亚硝酸也是DNRA的第一步。同样，我们将DNRA分为两种途径:异硝酸还原成亚硝酸盐和异硝酸还原成铵(图2A)。在细菌中，亚硝酸盐对铵途径的同化与其他途径呈负相关。同样，细菌和古菌中氨同化途径的相对丰度与其他N代谢均呈负相关。N固定途径一般与其他途径的相关性较弱或无相关性。
 

 
图2a 氮循环代谢的八种途径

 
3.2 N代谢路径的驱动者。
 
所有氮循环代谢的相对丰度(在所有途径上的总和)在土壤样品中差异很大，初步分析揭示了广泛的生物地理模式。平均而言，在热带森林和以人为主导(牧场、草坪和农田)的土壤中检测到了最高相对丰度，而在寒冷的沙漠中观察到的相对丰度最低(图S3)。N通路总相对丰度也随着纬度的增加呈下降趋势(R2=0.22,P<0.05);图S4)。为了解开这些模式背后的驱动因素，我们进行了多变量回归分析，包括已知的影响微生物丰度和组成的栖息地类型和环境参数。大多数样品没有当地的测量数据;相反，我们从辅助来源估计了这些变量。对于细菌而言，回归模型解释了总N代谢相对丰度中很大比例的可变性(R2=0.58,P<0.001）。栖息地类型对该模型的贡献最大，无论是直接(与总氮途径正相关)还是与土壤碳和N的相互作用。古菌回归模型解释总氮途径相对丰度的变异性低于细菌(R2=0.43,P<0.001)。碳和氮的交互作用对模型的贡献最大，而生境仅通过与温度的交互作用而重要。
 

 
图2b 氮循环代谢途径的相对丰度

 
接下来我们分析了单个N通路相对丰度的驱动因素。由于各通道之间的高协方差(图3 A和B)，我们对每个通道的全相对丰度校正残差拟合回归模型。这些模型解释这一额外变异的能力差异很大(表1)。例如，N固定途径模型分别解释了细菌和古菌样本之间80%和63%的变异(P< 0.001)。相比之下，相同的参数不能解释细菌中亚硝酸盐向铵化途径的差异相对丰度。在重要模式中，生境类型是单个路径相对丰度的重要预测因子(表1)。生境还与降水、温度、土壤氮等其他因素相互作用，影响某些路径的相对丰度。例如，在沙漠中反硝化相对丰度随温度升高而升高，而在热带森林中随温度升高而降低。同样，温带森林氨同化相对丰度随土壤氮增加而增加，热带森林氨同化相对丰度随土壤氮减少。土壤碳似乎是总氮途径相对丰度的主要驱动因素，但不能解释细菌中单个途径相对丰度的差异。在这些模型中包括N沉降的估计数只改善了反硝化模型;反硝化相对丰度随氮沉降量的增加而增加。古菌中单个通路相对丰度的模型通常比细菌的模型解释的变异更少，可能是由于每个样本的序列较少。然而，对于重要的模型来说，单个N通路通常可以用与细菌相同的参数来最好地解释。以生境类型和温度为主要指标。这两个区域内N固定相对丰度的重要预测因子。同样，生境、降水生境和温度生境也导致了古菌和细菌中硝酸盐对亚硝酸盐同化相对丰度的变化。
 

 
图3 氮循环代谢途径之间的相关性

 
3.3 N代谢途径的微生物物种组成和系统发育分布
 
我们发现多种多样的微生物，包括402细菌和53古菌属（Genera)，具有N代谢基因。我们首先研究了同一属（Genera)内通路的相关性(图4和图S5)。我们检测到的10个以上序列的所有属（Genera)均具有氨同化途径。携带反硝化基因的属（Genera)通常也携带前半段反硝化途径基因，将硝酸盐异化为亚硝酸盐。一些属同时包含所有三个途径。事实上，反硝化作用和DNRA最近被证明在同一细菌中存在并发挥作用。我们还检测了只携带异化亚硝酸盐到铵途径(除了氨同化)的属。
 
更确切地说，土壤属和它们所属的门，在它们的氮代谢路径专一化程度上是不同的。蓝藻属似乎是专一性代谢，主要携带同化亚硝酸盐到铵和氮固定途径。相比之下，Delta变形细菌似乎为广泛性代谢，拥有多达6种途径(除了氨同化)。土壤中N固定途径在古菌和细菌中广泛分布(图4和图S5)。变形菌门中发现了丰度最高的固氮微生物，而绿菌门、厚壁菌门和蓝藻门中明显例外(图5A)。大部分序列与N固定属密切相关，可能是缓生根瘤菌和伯克霍德氏菌。其他丰富的属较少被预料到。例如，Azoarcus是一种能降解土壤污染物能力的微生物，而果胶杆菌属(Pectobacterium)主要被认为是一种植物病原菌。事实上,尽管果胶杆菌属编码一组N固定基因，尚不清楚它们是否发挥功能。其次，编码异化亚硝酸盐还原为铵的途径在土壤细菌中也广泛分布(图4)。然而，优势土壤类群局限于两个门，Deltaproteobacteria和Verrucomicrobia(图5B)。已知土壤中有大量的Verrucomicrobia，但它们的生态作用尚不清楚(37,38)。该途径中最丰富的属，Anaeromyxobacter，在农业土壤中很常见，其最近被证明可以进行一种之前未被发现的将N2O还原为N2的非反硝化过程。图S6分析了在土壤样品中编码其他6个N代谢的属的相对丰度。
 

 
 图4 土壤宏基因组中氮循环物种的进化关系

 
3.4 N代谢途径对微生物物种组成的驱动
 
同样的环境因子可以很好地解释N代谢的总体相对丰度，但却不能很好地解释编码这些途径的微生物的分类学组成的变化。对于8种途径，这一模型仅解释了单一N代谢变化的7-19%。然而，就途径相对丰度而言，生境类型是组成的最佳预测因子，解释了高达14%的同化亚硝酸盐到铵途径的组成变化。温度也解释了11%的硝化途径组分变化。所有其他预测因子，包括pH值，最多能解释3%的任何途径的变异。
 
对两种途径的进一步研究证实了栖息地之间微弱的组成差异。携带N固定途径的15个最丰富的属在所有生境中都同样丰富，除了在寒冷的沙漠中(图5A)。编码亚硝酸盐向铵化途径的最丰富的属在不同生境中表现出更大的变异性，但其中只有一个绿菌属(Chlorobium)似乎在湿地上具有特异性(图5B)。
 

 
图5 编码特定代谢途径的物种分布

 
4 参考文献
 
Nelson, Michaeline B.; Martiny, Adam C.; Martiny, Jennifer B. H. (2016). Global biogeography of microbial nitrogen-cycling traits in soil. Proceedings of the National Academy of Sciences, 113(29), 8033–8040. doi:10.1073/pnas.1601070113 
 
猜你喜欢
 
10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
 
系列教程：微生物组入门 Biostar 微生物组  宏基因组
 
专业技能：学术图表 高分文章 生信宝典 不可或缺的人
 
一文读懂：宏基因组 寄生虫益处 进化树
 
必备技能：提问 搜索  Endnote
 
文献阅读 热心肠 SemanticScholar Geenmedical
 
扩增子分析：图表解读 分析流程 统计绘图
 
16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
 
在线工具：16S预测培养基 生信绘图
 
科研经验：云笔记  云协作 公众号
 
编程模板: Shell  R Perl
 
生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  
 
写在后面
 
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

 

 
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

 
点击阅读原文，跳转最新文章目录阅读
 

 
SIFT特征提取与检索
 
1. SIFT算法
1.1 基本概念
1.2 SIFT算法基本原理
1.2.1 特征点
1.2.2 尺度空间
1.2.3 高斯函数
1.2.4 高斯模糊
1.2.5 高斯金字塔
1.2.6 DOG函数
1.2.7 高斯差分金字塔
1.2.8 DOG局部极值检测
1.2.9 关键点方向分配
1.2.10 关键点描述
1.2.11 关键点匹配
   
1.3 SIFT算法实现步骤
  
2. SIFT算法实验过程
2.1 构造小型数据集
2.2 SIFT特征提取
2.2.1 代码实现
2.2.2 实验结果
2.2.3 实验结果分析
   
2.3 SIFT特征匹配
2.2.1 代码实现
2.2.2 实验结果
2.2.3 实验结果分析
   
2.4 SIFT特征检索并排序
2.2.1 代码实现
2.2.2 实验结果
2.2.3 实验结果分析
  
  
3. 匹配地理标记图像
3.1 实验代码
3.2 实验结果
3.3 实验结果分析
3.4 实验遇到的问题及解决
  
4. RANSAC算法
4.1 算法简介
4.2 单应性矩阵
4.3 RANSAC求解单应矩阵
4.2 实验代码
4.3 实验结果及分析
4.3.1 景深单一场景
4.3.2 景深复杂场景
   
4.4 实验总结
4.5 实验中遇到的问题及解决办法
  
5. 实验总结与分析
6. 实验中遇到的问题

 
1. SIFT算法
 
1.1 基本概念
 
尺度不变特征转换(Scale-invariant feature transform或SIFT) 是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。
 1999年David G.Lowe教授总结了基于特征不变技术的检测方法，在图像尺度 空间基础上，提出了对图像缩放、旋转保持不变性的图像局部特征描述算子 －SIFT（尺度不变特征变换），该算法在2004年被加以完善。其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。
 局部影像特征的描述与侦测可以帮助辨识物体，SIFT 特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。
 
SIFT算法可以解决的问题：
 
目标的旋转、缩放、平移（RST）
图像仿射/投影变换（视点viewpoint）
弱光照影响（illumination）
部分目标遮挡（occlusion）
杂物场景（clutter）
噪声
 
1.2 SIFT算法基本原理
 
1.2.1 特征点
 
SIFT中要查找的关键点：
 一些十分突出的点不会因光照、尺度、旋转等因素的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点。既然两幅图像中有相同的景物，那么使用某种方法分别提取各自的稳定点，这些点之间会有相互对应的匹配点。
 
1.2.2 尺度空间
 
尺度空间理论最早于1962年提出，其主要思想是通过 对原始图像进行尺度变换，获得图像多尺度下的空间表示。从而实现边缘、角点检测和不同分辨率上的特征提取，以满足特征点的尺度不变性。
 尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。尺度越大图像越模糊。
 
 
 
一个图像的尺度空间$L(x,y,\sigma )$,定义为一个变化尺度的高斯函数$G(x,y,\sigma )$与原图像$I(x,y)$的卷积：
 $$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)$$ 其中$\ast$表示卷积运算。
 
 
1.2.3 高斯函数
 
$$G(x,y,\sigma )=\frac{1}{2\Pi {\sigma }^2}{e}^{\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}$$
 
其中，$m,n$表示高斯模板的维度(由确定)。$(x,y)$代表图像的像素位置。$\sigma$是尺度空间因子，值越小表示图像被平滑的越少，相应的尺度也就越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。
 
1.2.4 高斯模糊
 
高斯模糊是在Adobe Photoshop等图像处理软件中广泛使用的处理 效果，通常用它来减小图像噪声以及降低细节层次。这种模糊技术生成 的图像的视觉效果是好像经过一个半透明的屏幕观察图像。如下图效果：
 
 
1.2.5 高斯金字塔
 
图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：
 $$n={\log }_{2}min(M,N)-t,t\in [0,{\log }_{2}min(M,N)]$$
 其中$M，N$为原图像的大小,$t$为塔顶图像的最小维数的对数值。
 如，对于大小为512512的图像，金字塔上各层图像的大小如表3.1所示，当塔顶图像为44时，n=7，当塔顶图像为2*2时，n=8。
 
图像大小
512
216
128
64
16
8
4
2
1
金字塔层数
1
2
3
4
5
6
7
8
9
尺度空间在实现时使用高斯金字塔表示,高斯金子塔的构建过程可分为两步：
 ① 对图像做高斯平滑
 ② 对图像做降采样
 为了让尺度体现其连续性，在简单下采样的基础上加上了高斯滤波。，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，每组含有多张(也叫层Interval)图像。另外，降采样时，高斯金字塔上一组图像的初始图像(底层图像)是由前一组图像的倒数第三张图像隔点采样得到的。
 
 高斯图像金字塔共o组、s层， 则有：
 $$\sigma (s)={\sigma }_0{2}^{s/S}$$
 其中，σ为尺度空间坐标； s为sub-level层坐标； ${\sigma }_0$为初始尺度； S为每组层数（一般为3~5）。
 最后可将组内和组间尺度 归为：
 $$2^{i-1}(\sigma ,k\sigma ,k^2\sigma ,...,k^{n-1}\sigma )$$
 $$k=2^{1/S}$$
 其中，i为金字塔组数 n为每一组的层数。
 
1.2.6 DOG函数
 
用更高效的高斯差分算子代替拉普拉斯算子进行极值检测，公式如下：
 $$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)$$
 $$D(x,y,\sigma )=[G(x,y,k\sigma )-G(x,y,\sigma ]\ast I(x,y)=L(x,y,k\sigma )-L(x,y,\sigma )$$
 由上式可知：DoG在计算上只需相邻高斯平滑后图像相减，因此简化了计算。
 
1.2.7 高斯差分金字塔
 
对应DOG算子，需构建DOG金字塔。可以通过高斯差分图像看出图像上的像素值变化情况。（如果没有变化，也就没有特征。特征必须是变化尽可能多的点。） DOG图像描绘的是目标的轮廓。
 在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如下图所示，进行极值检测。
 
 
1.2.8 DOG局部极值检测
 
DoG的局部极值点
 特征点是由DOG空间的局部极值点组成的。为了寻找DoG函数的极值点， 每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。
 中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个 点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。
去除边缘响应
 由于DoG函数在图像边缘有较强的边缘响应，因此需要排除边缘响应。 DoG函数的峰值点在边缘方向有较大的主曲率，而在垂直边缘的方向有 较小的主曲率。主曲率可以通过计算在该点位置尺度的2×2的Hessian矩 阵得到，导数由采样点相邻差来估计：
 $$H=\left[\begin{array}{cc}{D}_{xx}& D_{xy}\\ D_{xy}& D_{yy}\end{array}\right]$$
 $D_{xx}$表示DOG金字塔中某一尺度的图像x方向求导两次。
 $D$的主曲率和$H$的特征值成正比。令$\alpha$ ，$\beta$为特征值，则
 $$\frac{Tr(H{)}^2}{Det(H)}=\frac{(\alpha +\beta {)}^2}{\alpha \beta }$$
 $$det(H)=\alpha \beta$$
 $$trace(H)=\alpha +\beta$$
 该值在两特征值相等时达最小。Lowe论文中建议阈值T为1.2，即$\frac{Tr(H{)}^2}{Det(H)}<T$时保留关键点，反之剔除。
 
1.2.9 关键点方向分配
 
通过尺度不变性求极值点，可以使其具有缩放不变的性质。而利用关键点邻域像素的梯度方向分布特性，可以为每个关键点指定方向参数方向，从而使描述子对图像旋转具有不变性。通过求每个极值点的梯度来为极值点赋予方向。
 
像素点的梯度表示
 $$gradI(x,y)=\left(\begin{array}{c}\frac{\partial I}{\partial x},\frac{\partial I}{\partial y}\end{array}\right)$$
 梯度幅值:
 $$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)=ta{n}^{-1}\left[\begin{array}{c}\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\end{array}\right]$$
方向直方图的生成
 确定关键点的方向采用梯度直方图统计法，统计以关键点为原 点，一定区域内的图像像素点对关键点方向生成所作的贡献。如下图：
 
关键点的主方向与辅方向
 关键点主方向：极值点周围区域梯度直方图的主峰值也是特征点方向 。
 关键点辅方向：在梯度方向直方图中，当存在另一个相当于主峰值80%能量的峰值时，则将这个方向认为是该关键点的辅方向。
 这可以增强匹配的鲁棒性，Lowe的论文指出大概有15%关键点具有 多方向，但这些点对匹配的稳定性至为关键。
 
1.2.10 关键点描述
 
下图是一个SIFT描述子事例。其中描述子由2×2×8维向量表征，也即是 2×2个8方向的方向直方图组成。左图的种子点由8×8单元组成。每一个小格 都代表了特征点邻域所在的尺度空间的一个像素，箭头方向代表了像素梯度方 向，箭头长度代表该像素的幅值。然后在4×4的窗口内计算8个方向的梯度方 向直方图。绘制每个梯度方向的累加可形成一个种子点，如右图所示：一个特 征点由4个种子点的信息所组成。
 
 Lowe实验结果 表明：描述子采用4×4×8＝128维向量表征， 综合效果最优 （不变性与独 特性）。
 
 
1.2.11 关键点匹配
 
分别对模板图（参考图，reference image）和实时图（观测图，observation image）建立关键点描述子集合。目标的识别是通过两点集内关键点描述子的比对来完成。具有128维的关键点描述子的相似 性度量采用欧式距离。
 
 穷举匹配将原图像中的每一个特征点与目标图像中的每一个特征点进行匹配，如下图所示：
 
 模板图中关键点描述子：$R_i=(r_{i1},r_{i2},...,r_{i128})$
 实时图中关键点描述子：$S_i=(s_{i1},s_{i2},...,s_{i128})$
 任意两描述子相似性度量：$d(R_i,S_i)=\sqrt{{\sum }_{j=1}^{128}(r_{ij}-s_{ij}{)}^2}$
 要得到配对的关键点描述子，$d(R_i,S_i)$需满足：
 $$\frac{实时图中距离R_i最近的点S_j}{实时图中距离R_i的次最近的点S_p}＜Threshold$$
 关键点的匹配可以采用穷举法来完成，但是这样耗费的时间太多，一 般都采用Kd树的数据结构来完成搜索。搜索的内容是以目标图像的关 键点为基准，搜索与目标图像的特征点最邻近的原图像特征点和次邻 近的原图像特征点。Kd树是一个平衡二叉树。
 
 
1.3 SIFT算法实现步骤
 
 
 
SIFT算法的实质可以归为在不同尺度空间上查找特征点（关键点）的问题。
 
 
SIFT算法实现特征匹配的流程：
 
提取关键点.
对关键点附加详细的信息(局部特征),即描述符.
通过特征点(附带上特征向量的关 键点)的两两比较找出相互匹配的若干对特征点,建立景物间的对应关系.
 
 
2. SIFT算法实验过程
 
2.1 构造小型数据集
 
针对SIFT特征提取和检测算法，构造了一个小型数据集进行接下来的实验过程，数据集如下图：
 
 
2.2 SIFT特征提取
 
针对数据集的每一张图像实现SIFT特征提取并展示，对比Harris检测的效果。
 
2.2.1 代码实现
 
# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
from PCV.localdescriptors import sift
from PCV.localdescriptors import harris
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
imname='F:\\mr\\JMUclass\\ComputerVison\\CVcode\\test\\test2\\dataset\\1.jpg'
im = array(Image.open(imname).convert('L'))
sift.process_image(imname, 'empire.sift')
l1, d1 = sift.read_features_from_file('empire.sift')
figure()
gray()
subplot(221)
plt.axis('off')
imshow(Image.open(imname))
title(u'原图',fontproperties=font)
subplot(222)
sift.plot_features(im, l1, circle=False)
title(u'SIFT特征',fontproperties=font)
subplot(223)
sift.plot_features(im, l1, circle=True)
title(u'用圆圈表示SIFT特征尺度',fontproperties=font)
# 检测harris角点
harrisim = harris.compute_harris_response(im)
subplot(224)
filtered_coords = harris.get_harris_points(harrisim, 6, 0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
title(u'Harris角点',fontproperties=font)
show()
 
2.2.2 实验结果
 

    
      实验结果1
    

    
      实验结果2
    

    
      实验结果3
    

    
      实验结果4
    

    
      实验结果5
    

    
      实验结果6
    

    
      实验结果7
    

    
      实验结果8
    

    
      实验结果9
    

    
      实验结果10
    

    
      实验结果11
    

    
      实验结果12
    

    
      实验结果13
    

    
      实验结果14
    

    
      实验结果15
    
 
2.2.3 实验结果分析
 
由以上实验结果可得如下分析：
 
SIFT特征提取和检测算法提取到的特征点明显多于Harris角点检测，由于在一定程度上增加了图像的特征点，从后续的特征匹配的工作角度上考虑可知：SIFT算法明显优于Harris算法。在一定程度上这体现了SIFT特征提取与检测算法的优越性。
在实验过程中，由实验结果中的3、5、7、12、13中可以看出Harris角点检测并未检测到任何特征点，而SIFT算法却极大程度地检测出特征点和特征尺度，根据实验结果11、12等可以看出SIFT提取到的特征点几乎覆盖了图像的主体内容，除了白色地墙壁部分其他部分几乎布满了特征点，这体现了SIFT算法的稳定性和极大的鲁棒性。
根据实验结果7、8、9可以分析出SIFT具有多量性，即使少数的几个物体也可以产生大量SIFT特征向量，图像7、8、9结构简单，涵盖的物体较少，但SIFT算法提取出大量的特征点。
从时间效率上来分析，SIFT特征提取和检测消耗了大量的时间，对于一些特征点丰富且像素点较多的图像的检测时间多达5分钟左右，但其提取特征的准确性不可忽略，大可忽略其时间效率影响。可进行进一步分析，对SIFT算法进行简单优化。
 
总结SIFT算法和Harris算法的区别如下：
 
Harris特征点提取的精度低和效率不高,运算慢、抗噪能力差以及在实际应用中存在不必要角点簇，Harris特征点检测算法采用高斯滤波进行平滑,图像存在角点信息丢失与偏移的现象。
SIFT方法对于旋转、尺度缩放、亮度变化保持不变性，而且对视角变化、仿射变换、噪声也保持一定程度的稳定性，特征点的个数和有效点的比例没有要求。
 
2.3 SIFT特征匹配
 
2.2.1 代码实现
 
from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift

if len(sys.argv) >= 3:
    im1f, im2f = sys.argv[1], sys.argv[2]
else:
    im1f = 'F:\\mr\\JMUclass\\ComputerVison\\CVcode\\test\\test2\\1.jpg'
    im2f = 'F:\\mr\\JMUclass\\ComputerVison\\CVcode\\test\\test2\\2.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)
sift.process_image(im2f, 'out_sift_2.txt')
l2, d2 = sift.read_features_from_file('out_sift_2.txt')
subplot(122)
sift.plot_features(im2, l2, circle=False)
# matches = sift.match(d1, d2)
matches = sift.match_twosided(d1, d2)
print ('{} matches'.format(len(matches.nonzero()[0])))
figure()
gray()
sift.plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()
 
2.2.2 实验结果
 
1. 实验结果1
 控制台输出：
 
processed tmp.pgm to out_sift_1.txt
processed tmp.pgm to out_sift_2.txt
4078 matches
[Finished in 1109.8s]
 
可视化匹配效果：
 
 2. 实验结果2
 控制台输出：
 
processed tmp.pgm to out_sift_1.txt
processed tmp.pgm to out_sift_2.txt
679 matches
 
可视化匹配效果：
 
 3. 实验结果3
 控制台输出：
 
processed tmp.pgm to out_sift_1.txt
processed tmp.pgm to out_sift_2.txt
55 matches
 
可视化匹配效果：
 
 
2.2.3 实验结果分析
 
由以上三组实验可得以下分析：
 
在SIFT特征匹配中，针对如上第一个=组实验结果，其运行时间消耗极大，大约执行了约30min，分析可得第一组图像的像素较大，其一张图像共有3024* 4032个像素块，且两张图像的主题一致，其匹配到的特征点共有4078个matches，所以其运行时间消耗大是有原因的。相对来说，第二组和第三组实验中，其像素大小为3024* 3024，相对较小，且其匹配到的特征点的数量是第一组实验的16%和1.2%，所以第二组和第三组的运行时间较快。
针对第一组实验可得：其图像中的建筑物主体几乎完全匹配，由此可以看出SIFT特征匹配的优越表现。
将第一组实验和第二、三组实验进行对比，发现SIFT提取到的特征数量对于特征匹配有一定的影响。第一组实验提取到的特征点多，在后续进行特征匹配的可能性更多，进而获得的matches值较大，而相对于第二、三组实验，其提取到的特征数量较为少，进而得到的matches值也较小。
观察第二组实验结果，其SIFT匹配结果仅为字画的内容匹配较多，而针对于字画的边框并未出现匹配结果，思考可能由于第二组实验中第二张图像对于边框的特征点标注较少，由于其发生旋转后并未寻找到旋转不变性的关键点，从而在匹配过程中得到的关于边框的匹配结果为0。
当特征点不是很多时，经优化的SIFT匹配算法甚至可以达到实时的要求。而且可以很方便的与其他形式的特征向量进行联合。特征点的个数和有效点的比例没有要求。当特征点不是很多时，经优化的SIFT匹配算法甚至可以达到实时的要求。而且可以很方便的与其他形式的特征向量进行联合，SIFT算法也有一些不足。方法通过对特征点构造128维的向量，然后对向量进行匹配，这样图像就得满足足够多的纹理，否则构造出的128维向量区别性就不是太大，容易造成误匹配，极限情况如指纹图像的匹配，星图识别等这类图像特征点周围根本没有什么纹理这时SIFT算法就完全失效了。
 
2.4 SIFT特征检索并排序
 
输入一张新的图像，针对数据集中的每一张图像分别与输入图像进行匹配，每一次匹配计算出其matches的值并记录下来。在数据集中的所有图像匹配完成之后，根据matches的大小排序，选取matches值最大的三张图像作为输出。
 
2.2.1 代码实现
 
#添加必要的库
import os
import numpy as np
import matplotlib.image as mp
from skimage import img_as_ubyte
from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
path = "F:\\mr\\JMUclass\\ComputerVison\\CVcode\\test\\test2\\dataset\\"
filelist = os.listdir(path)  # 打开对应的文件夹
total_num = len(filelist)-1  #得到文件夹中图像的个数
matches_array=np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])
result=np.array([0,0,0])
if len(sys.argv) >= 3:
    im1f, im2f = sys.argv[1], sys.argv[2]
else:
    im1f = 'F:\\mr\\JMUclass\\ComputerVison\\CVcode\\test\\test2\\new.jpg'
im1 = array(Image.open(im1f))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
subplot(2,2,1)
plt.axis('off')
title(u'原图',fontproperties=font)
imshow(Image.open(im1f))
for i in range(total_num):
    im2f = path + str(i + 1) + '.jpg' #拼接图像的读取地址
    im2 = array(Image.open(im2f)) 
    sift.process_image(im2f, 'out_sift_'+str(i+2)+'.txt')
    l2, d2 = sift.read_features_from_file('out_sift_'+str(i+2)+'.txt')
    matches = sift.match_twosided(d1, d2)
    matches_array[i] = len(matches.nonzero()[0])
print (matches_array)
for i in range(3):
	a = np.argmax(matches_array)
	im2f = path + str(a+1) + '.jpg' #拼接图像的读取地址
	subplot(2,2,i+2)
	plt.axis('off')
	mstr='matches：'+str(matches_array[a])
	title(mstr,fontproperties=font)
	imshow(Image.open(im2f))
	matches_array[a]=0
show()
 
2.2.2 实验结果
 
控制台输出：
 
processed tmp.pgm to out_sift_1.txt
processed tmp.pgm to out_sift_2.txt
processed tmp.pgm to out_sift_3.txt
processed tmp.pgm to out_sift_4.txt
processed tmp.pgm to out_sift_5.txt
processed tmp.pgm to out_sift_6.txt
processed tmp.pgm to out_sift_7.txt
processed tmp.pgm to out_sift_8.txt
processed tmp.pgm to out_sift_9.txt
processed tmp.pgm to out_sift_10.txt
processed tmp.pgm to out_sift_11.txt
processed tmp.pgm to out_sift_12.txt
processed tmp.pgm to out_sift_13.txt
processed tmp.pgm to out_sift_14.txt
processed tmp.pgm to out_sift_15.txt
processed tmp.pgm to out_sift_16.txt
[2953 1163  197   41    2    1    5    1    3    2    3    3    0    5   1]
[Finished in 1188.2s]
 
可视化匹配排序结果：
 
 
2.2.3 实验结果分析
 
在输入一张范围涵盖较大场景的图像与数据集中15张图像进行匹配，进而排序，可得观察到其实验结果的排序与我们目视判读的结果一致，可以得出SIFT特征匹配的准确率是相当高的。其实验效果好与高效提取关键点、准确地对关键点附加详细的信息(局部特征)、通过特征点(附带上特征向量的关 键点)的两两比较找出相互匹配的若干对特征点,建立景物间的对应关系有密不可分地关系。
有一个较大的问题在于该实验出结果的速度很慢，考虑如何对SIFT算法进行改进，有了如下思考：sift特征匹配算法128维的特征描述向量使得处理匹配特征点计算庞大，导致时效性不高。考虑将128维的特征描述向量进行降维处理，并且像素的描述范围增大，或许可以有效减少了匹配的运算次数，缩短图像配准时间。
 
3. 匹配地理标记图像
 
3.1 实验代码
 
# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot
import os
os.environ['PATH'] = os.environ['PATH'] + (';C:\\Program Files (x86)\\Graphviz2.38\\bin\\')
""" This is the example graph illustration of matching images from Figure 2-10.
To download the images, see ch2_download_panoramio.py."""
#download_path = "panoimages"  # set this to the path where you downloaded the panoramio images
#path = "/FULLPATH/panoimages/"  # path to save thumbnails (pydot needs the full system path)
download_path = "F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset"  # set this to the path where you downloaded the panoramio images
path = "F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset/"  # path to save thumbnails (pydot needs the full system path)
# list of downloaded filenames
imlist = imtools.get_imlist(download_path)
# print(imlist)
nbr_images = len(imlist)
# extract features
featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    sift.process_image(imname, featlist[i])
matchscores = zeros((nbr_images, nbr_images))
for i in range(nbr_images):
    for j in range(i, nbr_images):  # only compute upper triangle
        print ('comparing ', imlist[i], imlist[j])
        l1, d1 = sift.read_features_from_file(featlist[i])
        l2, d2 = sift.read_features_from_file(featlist[j])
        matches = sift.match_twosided(d1, d2)
        nbr_matches = sum(matches > 0)
        print ('number of matches = ', nbr_matches)
        matchscores[i, j] = nbr_matches
print ("The match scores is: \n", matchscores)
# copy values
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):  # no need to copy diagonal
        matchscores[j, i] = matchscores[i, j]
#可视化
threshold = 2  # min number of matches needed to create link
g = pydot.Dot(graph_type='graph')  # don't want the default directed graph
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):
        if matchscores[i, j] > threshold:
            # first image in pair
            im = Image.open(imlist[i])
            im.thumbnail((100, 100))
            filename = path + str(i) + '.png'
            im.save(filename)  # need temporary files of the right size
            g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))
            # second image in pair
            im = Image.open(imlist[j])
            im.thumbnail((100, 100))
            filename = path + str(j) + '.png'
            im.save(filename)  # need temporary files of the right size
            g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle', image=filename))

            g.add_edge(pydot.Edge(str(i), str(j)))
g.write_png('whitehouse.png')
 
3.2 实验结果
 
控制台输出：
 
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\1.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\10.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\11.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\12.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\13.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\14.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\15.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\2.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\3.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\4.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\5.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\6.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\7.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\8.sift
processed tmp.pgm to F:/mr/JMUclass/ComputerVison/CVcode/test/test2/dataset3\9.sift
The match scores is: 
 [[214.   1.   0.   0.   1.   0.   0.  39.  11.   5.   1.   1.   1.   0.
    0.]
 [  0. 245.  70.   1.   1.   0.   1.   1.   1.   0.   0.   1.   0.   2.
    1.]
 [  0.   0. 353.   0.   0.   2.   1.   0.   0.   0.   0.   1.   0.   4.
    2.]
 [  0.   0.   0. 344.   1.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    1.]
 [  0.   0.   0.   0. 164.  18.  11.   1.   0.   0.   1.   1.   0.   1.
    0.]
 [  0.   0.   0.   0.   0.  81.   8.   0.   0.   0.   2.   1.   2.   0.
    0.]
 [  0.   0.   0.   0.   0.   0.  46.   1.   0.   0.   0.   2.   0.   3.
    0.]
 [  0.   0.   0.   0.   0.   0.   0. 290.   2.   0.   0.   1.   0.   1.
    0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.  99.  31.   0.   0.   0.   1.
    0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0. 202.   0.   0.   0.   0.
    0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0. 178.  33.  42.   0.
    0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0. 125.   7.   2.
    2.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0. 132.   1.
    1.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  68.
    7.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
  158.]]
[Finished in 12.8s]
 
可视化结果：
 
 
3.3 实验结果分析
 
在进行匹配的时候代码中使用嵌套循环实现，即共进行了15*15次匹配。观察输出的匹配分数矩阵，由于读取图片的顺序是错乱的，分别是以1.jpg,10.jpg,11.jpg,12.jpg,13.jpg,14.jpg,15.jpg,2.jpg,3.jpg,4.jpg,5.jpg,6.jpg,7.jpg,8.jpg,9.jpg的顺序，即在观察分数矩阵的时候行和列的顺序应与上述顺序保持一致。
在控制台的输出结果可以看出一个明显的特点：在矩阵的对角线上，匹配分数都很大。原因是：对角线上的匹配分数其实是在每一轮的循环中，自身与自身进行匹配，匹配分数自然很高，所以出现了这种情况。
可视化的是实验结果可以看出，匹配连接的顺序大致正确，其中有两幅画不是同一幅但产生了一条连接线，主要是由于两幅画的主体内容大部分一致，但一定程度上还是存在较小的错误性。
在创建图像缩略图的时候，根据书本中的介绍将数据集中的图像创建为100* 100的图像出现的结果如实验结果中的截图，由于其出现了一张图像有错误，所以考虑将缩略图的大小略增大，增大为200* 200，结果如下图，观察发现部分连线正确，部分连线产生了很大的错误。所以推测缩略图的大小对匹配结果有较大的影响。
 
 
3.4 实验遇到的问题及解决
 
1.出现文件夹读取为空的问题：
 
 将数据集里的图片的后缀改为小写就可以，比如JPG改为jpg
 
2.出现运行很慢的问题：
 
pil_im.thumbnail((100,100))
 
使用该语句完成缩略图的创建即可。
 
3.报错"dot" not found in path.：
 
 解决方法：
 先安装pydot:
 
pip install pydot
 
然后再下载Graphviz：graphviz-2.38.msi
 
 接着安装graphviz，记住默认的安装路径C:\Program Files (x86)\Graphviz2.38
 将Graphviz的bin文件夹的位置配置环境变量path后，添加如下代码：
 
import os
os.environ['PATH'] = os.environ['PATH'] + (';C:\\Program Files (x86)\\Graphviz2.38\\bin\\')
 
测试之后还是报错，后来发现，错误提示是找不到“dot”程序，而不是"dot.exe"程序
 修改self.prog = 'dot’为self.prog = ‘dot.exe’，之后测试，成功运行。
 
4. RANSAC算法
 
RANSAC算法是RANdom SAmple Consensus的缩写，意为随机抽样一致。表面上的意思就是从匹配样本中随机取样，寻找一致的样本点。RANSAC算法是根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，得到有效样本数据的算法。它是在1981年由Fischler和Bolles最先提出。
 
在利用已有算法进行特征点匹配时，常存在的一个问题就是误匹配的问题，这些误匹配的点对匹配的效果产生很大的影响，所以我们需要利用一定的方法剔除误匹配的特征点，在实际应用中，我们常会用到RANSAC算法来消除两两匹配图像的误匹配点，这个算法现在在图像配准以及拼接上得到了广泛的应用。
 
RANSAC算法的核心思想就是在匹配的特征点中随机取4个特征点，通过计算和不断迭代，寻找到最优的参数模型，在这个最优模型中，能匹配上的特征点最多。
 
4.1 算法简介
 
RANSAC算法是一种简单且有效的去除噪声影响，估计模型的一种方法。与普通的去噪算法不同，RANSAC算法是使用尽可能少的点来估计模型参数，然后尽可能的扩大得到的模型参数的影响范围。
 RANSAC算法的具体描述是：给定$N$个数据点组成的集合$P$，假设集合中大多数的点都是可以通过一个模型来产生的，且最少$n$个点（$n<N$）可以拟合出模型的参数，则可以通过以下的迭代方式拟合该参数。
 
对下面的操作执行$k$次：
 
从$P$中随机选择$n$个数据点；
用这$n$个数据点拟合出一个模型$M$；
对$P$中剩余的数据点，计算每个点与模型$M$的距离，距离超过阈值的则认定为局外点，不超过阈值的认定为局内点，并记录该模型$M$所对应的局内点的值$m$；
迭代$k$次以后，选择$m$最大的模型$M$作为拟合的结果。
 
因为在实际应用中$N$的值通常会很大，那么从其中任选$n$个数据点的组合就会很大，如果对所有组合都进行上面的操作运算量就会很大，因此对于$k$的选择就很重要。通常情况下，只要保证模型估计需要的$n$个点都是点的概率足够高即可。因此设$w$为$N$个数据中局内点的比例，$z$为进行$k$次选取后，至少有一次选取的$n$个点都是局内点的概率。则有
 $$z=1-(1-w^n{)}^k$$
 其中$1-w^n$表示一次选取不都是局内点的概率，$(1-w^n{)}^k$表示$k$次选取中没有一次都是局内点的概率。 则有
 $$k=\frac{log(1-z)}{log(1-w^n)}$$
 
这里$z$一般要求满足大于95%即可。
 
4.2 单应性矩阵
 
单应性矩阵描述的是针对同一事物，在不同的视角下拍摄的两幅图像之间的关系。假设这两幅图像之间是透视变换，则单应性矩阵也就是透视变换矩阵$H$定义如下：
 $$H=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& 1\end{array}\right]$$
 则有:
 $$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$
 因此要恢复出$H$中的8个参数，至少需要4对匹配点，过程如下：
 

  

那么就可以每次从所有的匹配点中选出4对，计算单应性矩阵$H$，然后选出内点个数最多的作为最终的结果。计算距离方法如下：
 

  

4.3 RANSAC求解单应矩阵
 
· RANSAC loop：
 
随机选择四对匹配特征
根据DLT计算单应矩阵$H$(唯一解)
对所有匹配点，计算映射误差$\epsilon =\Vert {p_i}^{\prime },H{p}_i\Vert$
根据误差阈值，确定inliers
针对最大inliers集合，重新计算单应矩阵$H$
 
4.2 实验代码
 
该实验设计为先利用SIFT特征匹配对魅族图像进行特征匹配，进而利用RANSAC算法对误匹配点进行剔除，观察实验结果。
 
本次实验设计分为景深复杂和景深单一的两种场景来分析RANSAC算法的实验结果。
 
import cv2
import numpy as np
import random

def compute_fundamental(x1, x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")

    # build matrix for equations
    A = np.zeros((n, 9))
    for i in range(n):
        A[i] = [x1[0, i] * x2[0, i], x1[0, i] * x2[1, i], x1[0, i] * x2[2, i],
                x1[1, i] * x2[0, i], x1[1, i] * x2[1, i], x1[1, i] * x2[2, i],
                x1[2, i] * x2[0, i], x1[2, i] * x2[1, i], x1[2, i] * x2[2, i]]

    # compute linear least square solution
    U, S, V = np.linalg.svd(A)
    F = V[-1].reshape(3, 3)

    # constrain F
    # make rank 2 by zeroing out last singular value
    U, S, V = np.linalg.svd(F)
    S[2] = 0
    F = np.dot(U, np.dot(np.diag(S), V))

    return F / F[2, 2]


def compute_fundamental_normalized(x1, x2):
    """    Computes the fundamental matrix from corresponding points
        (x1,x2 3*n arrays) using the normalized 8 point algorithm. """

    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")

    # normalize image coordinates
    x1 = x1 / x1[2]
    mean_1 = np.mean(x1[:2], axis=1)
    S1 = np.sqrt(2) / np.std(x1[:2])
    T1 = np.array([[S1, 0, -S1 * mean_1[0]], [0, S1, -S1 * mean_1[1]], [0, 0, 1]])
    x1 = np.dot(T1, x1)

    x2 = x2 / x2[2]
    mean_2 = np.mean(x2[:2], axis=1)
    S2 = np.sqrt(2) / np.std(x2[:2])
    T2 = np.array([[S2, 0, -S2 * mean_2[0]], [0, S2, -S2 * mean_2[1]], [0, 0, 1]])
    x2 = np.dot(T2, x2)

    # compute F with the normalized coordinates
    F = compute_fundamental(x1, x2)
    # print (F)
    # reverse normalization
    F = np.dot(T1.T, np.dot(F, T2))

    return F / F[2, 2]

def randSeed(good, num =8 ):
    '''
    :param good: 初始的匹配点对
    :param num: 选择随机选取的点对数量
    :return: 8个点对list
    '''
    eight_point = random.sample(good, num)
    return eight_point

def PointCoordinates(eight_points, keypoints1, keypoints2):
    '''
    :param eight_points: 随机八点
    :param keypoints1: 点坐标
    :param keypoints2: 点坐标
    :return:8个点
    '''
    x1 = []
    x2 = []
    tuple_dim = (1.,)
    for i in eight_points:
        tuple_x1 = keypoints1[i[0].queryIdx].pt + tuple_dim
        tuple_x2 = keypoints2[i[0].trainIdx].pt + tuple_dim
        x1.append(tuple_x1)
        x2.append(tuple_x2)
    return np.array(x1, dtype=float), np.array(x2, dtype=float)


def ransac(good, keypoints1, keypoints2, confidence,iter_num):
    Max_num = 0
    good_F = np.zeros([3,3])
    inlier_points = []
    for i in range(iter_num):
        eight_points = randSeed(good)
        x1,x2 = PointCoordinates(eight_points, keypoints1, keypoints2)
        F = compute_fundamental_normalized(x1.T, x2.T)
        num, ransac_good = inlier(F, good, keypoints1, keypoints2, confidence)
        if num > Max_num:
            Max_num = num
            good_F = F
            inlier_points = ransac_good
    print(Max_num, good_F)
    return Max_num, good_F, inlier_points


def computeReprojError(x1, x2, F):
    """
    计算投影误差
    """
    ww = 1.0/(F[2,0]*x1[0]+F[2,1]*x1[1]+F[2,2])
    dx = (F[0,0]*x1[0]+F[0,1]*x1[1]+F[0,2])*ww - x2[0]
    dy = (F[1,0]*x1[0]+F[1,1]*x1[1]+F[1,2])*ww - x2[1]
    return dx*dx + dy*dy

def inlier(F,good, keypoints1,keypoints2,confidence):
    num = 0
    ransac_good = []
    x1, x2 = PointCoordinates(good, keypoints1, keypoints2)
    for i in range(len(x2)):
        line = F.dot(x1[i].T)
        #在对极几何中极线表达式为[A B C],Ax+By+C=0,  方向向量可以表示为[-B,A]
        line_v = np.array([-line[1], line[0]])
        err = h = np.linalg.norm(np.cross(x2[i,:2], line_v)/np.linalg.norm(line_v))
        # err = computeReprojError(x1[i], x2[i], F)
        if abs(err) < confidence:
            ransac_good.append(good[i])
            num += 1
    return num, ransac_good


if __name__ =='__main__':
    im1 = 'image/33.png'
    im2 = 'image/44.png'

    print(cv2.__version__)
    psd_img_1 = cv2.imread(im1, cv2.IMREAD_COLOR)
    psd_img_2 = cv2.imread(im2, cv2.IMREAD_COLOR)
    # 3) SIFT特征计算
    sift = cv2.xfeatures2d.SIFT_create()
    # find the keypoints and descriptors with SIFT
    kp1, des1 = sift.detectAndCompute(psd_img_1, None)
    kp2, des2 = sift.detectAndCompute(psd_img_2, None)

    # FLANN 参数设计
    match = cv2.BFMatcher()
    matches = match.knnMatch(des1, des2, k=2)

    # Apply ratio test
    # 比值测试，首先获取与 A距离最近的点 B （最近）和 C （次近），
    # 只有当 B/C 小于阀值时（0.75）才被认为是匹配，
    # 因为假设匹配是一一对应的，真正的匹配的理想距离为0
    good = []
    for m, n in matches:
        if m.distance < 0.75 * n.distance:
            good.append([m])
    print(good[0][0])

    print("number of feature points:",len(kp1), len(kp2))
    print(type(kp1[good[0][0].queryIdx].pt))
    print("good match num:{} good match points:".format(len(good)))
    for i in good:
        print(i[0].queryIdx, i[0].trainIdx)


    Max_num, good_F, inlier_points = ransac(good, kp1, kp2, confidence=30, iter_num=500)
    # cv2.drawMatchesKnn expects list of lists as matches.
    # img3 = np.ndarray([2, 2])
    # img3 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, good[:10], img3, flags=2)

    # cv2.drawMatchesKnn expects list of lists as matches.

    img3 = cv2.drawMatchesKnn(psd_img_1,kp1,psd_img_2,kp2,good,None,flags=2)
    img4 = cv2.drawMatchesKnn(psd_img_1,kp1,psd_img_2,kp2,inlier_points,None,flags=2)
    cv2.namedWindow('image1', cv2.WINDOW_NORMAL)
    cv2.namedWindow('image2', cv2.WINDOW_NORMAL)
    cv2.imshow("image1",img3)
    cv2.imshow("image2",img4)
    cv2.waitKey(0)#等待按键按下
    cv2.destroyAllWindows()#清除所有窗口
 
4.3 实验结果及分析
 
4.3.1 景深单一场景
 
控制台输出：
 
3.3.0
<DMatch 00000223D14F2450>
number of feature points: 3629 4108
<class 'tuple'>
good match num:975 good match points:
[51 1314,62 604,73 1314,82 604,99 1314,106 604,141 2859,155 2859,165 3296,169 622,234 1314......]
260 [[ 2.70487221e-07  2.30359993e-06 -3.19941970e-03]
 [-2.24579797e-06 -9.47445524e-09  2.94200972e-03]
 [ 2.83485400e-03 -3.80271280e-03  1.00000000e+00]]
 
可视化结果：
 RANSAC前:
 

  

RANSAC后:
 

  

实验小结：
 
该组实验使用景深较为单一的场景图片进行实验，发现其删除匹配线的效果很明显，在进行SIFT特征匹配时，由于画框的颜色材质与天花板以及场景中床头装饰的颜色材质较为一致，所以发生了较多错乱的匹配线。
在经过RANSAC算法过滤后删除了大多数的错配点，但是仔细观察RANSAC后的实验结果，仍发现一条在画框和天花板装饰之间的粉色匹配线；观察画中特征点的匹配效果得到，画中扇形右下角的特征点与左下角的特征点有一条匹配线，两点并非同一点其周围以及其本身相似程度较大产生了匹配线，而RANSAC算法并没有排除这条匹配先；最后观察画中上半部分原有很多匹配线，而RANSAC后画中上半部分已全无匹配线，可知在RANSAC算法进行中可能删除了大部分匹配正确的匹线。由以上可以得出结论：RANSAC算法并非百分之百准确，仍有优化的空间。
思考：在这组实验中结果并没有那么理想，出现了较多以上叙述的RANSAC算法的筛选不正确性，思考如果增加迭代次数，结果是否会优化。由于以往实验的经验，增加迭代次数并非会带来好的结果，对于这一问题将在之后进行多次实验进行观察。
 
4.3.2 景深复杂场景
 
控制台输出：
 
3.3.0
<DMatch 0000013D09470BD0>
number of feature points: 4544 3785
<class 'tuple'>
good match num:396 good match points:[18 190,19 140,25 312,45 129,60 114,62 139,107 129,116 129,127 114,131 139,133 190,142 175,224 241,240 248,264 277,291 1739,400 1462,408 1073,438 1345,503 569,514 547,543 410,550 569,......]  
228 [[ 4.24486619e-06 -8.04672331e-06 -3.21653120e-03]
 [ 7.58870199e-06  1.42993003e-06 -3.03983563e-03]
 [-3.83888038e-05  5.73794118e-03  1.00000000e+00]]
 
可视化结果：
 RANSAC前:
 

  

RANSAC后:
 

  

实验分析：
 
该组实验使用景深较为丰富的场景图片进行实验，由可视化结果可以看到大多数误匹配的特征点之间的连线被删除，进行RANSAC后的匹配结果更加整齐，并无原匹配结果的错乱复杂。
在这组实验中更好地展现了RANSAC算法的优势，观察建筑物之间的错匹配已成功删除、天空中存在的错匹配也成功删除、街边路灯与建筑物之间的错匹配也成功删除…，由以上可得RANSAC算法求解模型的准确性较高，将匹配点带入求解地模型其中可以较大程度地检测出其不正确性进而删除匹配线。可以得到结论：RANSAC算法经过迭代得到的模型较为准确，模型的准确性决定着删除误匹配的效果。
在该组实验中，程序运行速度很慢长达5min，后来经过思考知道是由于图像本身像素较大，在进行SIFT匹配时速度较慢，按照比例缩小图像的时候运行速度得到改善，同时发现了图像比例不同得到的结果也有所不同。
 
4.4 实验总结
 
RANSAC的优点是它能鲁棒的估计模型参数。例如，它能从包含大量局外点的数据集中估计出高精度的参数。
 RANSAC的缺点是它计算参数的迭代次数没有上限；如果设置迭代次数的上限，得到的结果可能不是最优的结果，甚至可能得到错误的结果。RANSAC只有一定的概率得到可信的模型，概率与迭代次数成正比。RANSAC的另一个缺点是它要求设置跟问题相关的阀值。
 RANSAC只能从特定的数据集中估计出一个模型，如果存在两个（或多个）模型，RANSAC不能找到别的模型。
 
4.5 实验中遇到的问题及解决办法
 
遇到报错：ValueError: Sample larger than population or is negative
 解决方法：
 开始考虑是否为random.py中的函数发生错误，经过观察可得代码中的带的匹配点的数量为0，将0值传入随机函数中发生错误。进而思考为何匹配点数量为0，经过查看图像详细信息发现两张图像的尺寸不同，使用resize()函数将图像尺寸改为一致再进行实验即可成功运行。
遇到报错：
 cv2. error: OpenCV(4.1.2) C: \projects lopencv- python \opencv_ contrib \modules \xfeatures2d\src\sift.cpp:1207: error: ( - 213:The function/ feature is notimplemented) This algorithm is patented and is excluded in this configuration; Set OPENCV ENABLE_ NONFREE CMake option and rebuild the library infunction ’ cv: :xfeatures2d: :SIFT : :create
 解决方法：
 在cmd中运行如下命令即可。
 
pip install opencv-contrib-python==3.3.0.10
 
 
5. 实验总结与分析
 
总体来说，SIFT算法具有以下特性：
 
Sift特征是图像的局部特征，对平移、旋转、尺度缩放、亮度变化、遮挡和噪声等具有良好的不变性，对视觉变化、仿射变换也保持一定程度的稳定性。
SIFT算法独特性好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配。在实验过程中有明确的体现，尤其在SIFT算法与Harris算法的对比中可以看出SIFT算法检测到的信息量丰富。
SIFT算法具有多量性，即使少数的几个物体也可以产生大量Sift特征向量。
SIFT速度相对较快，经查阅可得目前有较多论文对SIFT算法进行了较大程度的优化，经优化的Sift匹配算法甚至可以达到实时的要求。
SIFT算法的可扩展性强，可以很方便的与其他形式的特征向量进行联合。
 
SIFT在图像的不变特征提取方面拥有较大的优势，但仍然存在以下不足：
 
 
实时性不高，主要是由于其运行时间消耗太大，而对于实时性的问题，需要在几毫秒之内或1s之内出结果，而在本次实验中稍微过于复杂的图像进行SIFT匹配时耗时可达30min，这几乎对实时性的问题无法进行很好的解决。
 
 
SIFT算法在进行特征提取的过程中，有时特征点较少。由本次实验中的部分实验结果可以体现。
 
 
SIFT算法对边缘光滑的目标无法准确提取特征点。由于其对图像进行了高斯模糊，所以对于平滑的目标提取其轮廓较为困难。就像上面提到的SIFT特征匹配对两幅字画进行匹配时，只对于图像的主体内容进行匹配，而对字画的边框并无匹配点。
 
 
6. 实验中遇到的问题
 
问题描述如下：
 
OSError: empire.sift not found
 
原因分析：
 报错提示的代码行为l1, d1 = sift.read_features_from_file(‘empire.sift’)，说明没有读取到empire.sift，说明在sift.process_image阶段没有生成sift文件
 
解决方案：
 为了计算图像的SIFT特征，我们需要用到开源工具包VLFeat。下载链接 ：“http://www.vlfeat.org/download/
 注意：这里应该下载0.9.20版本的才可用，我下载的是vlfeat-0.9.20-bin.tar.gz。下载完后解压）。在下载过程中，由于其为外网，下载速度很慢且会出现多次网络错误中断，多次尝试即可解决，如果不行可以联系我把压缩包发送给你。
 
 接下来需要进行的操作步骤：
 
 
把vlfeat文件夹bin下win64中的sift.exe和vl.dll这两个文件复制到项目的文件夹中。
 
 
修改Anaconda文件夹下的PCV（我的PCV位置E:\Anaconda3\Lib\site-packages\PCV）文件夹里面的localdescriptors文件夹中的sift.py文件，使用记事本打开，修改其中的cmmd内的路径为：cmmd = str(r"D:\PythonWork\SIFT\sift.exe “+imagename+” --output="+resultname+" "+params) （路径是你项目文件夹中的sift.exe的路径），一定要记得在括号内加r，否则还是会出错。
 
 
之后就可以运行了。如果在运行过程中提示关于print的错误，记得根据错误提醒的文件夹，去修改相应的print语法，3.5的python的print用法是需要加括号。
 

 
高成远
 
GIS即地理信息系统，是近些年迅速发展起来的一门空间信息分析技术，它是一门综合性学科，结合地理学与地图学以及遥感和计算机科学，已经广泛的应用在不同的领域，是用于输入、存储、查询、分析和显示地理数据的计算机系统。GIS 技术不仅可以对地球表层空间中的地理信息数据进行采集、存储和管理，还能有效地对所得数据进行运算、分析、描述以及显示，从而使其具有更为广泛的实用价值。
 
在开展环境保护工作与建设环境地理信息系统的实践过程中，面对我国环境污染和生态破坏的严峻形势和艰巨繁杂的环境保护任务目标，我国各级环保部门已逐渐认识到及时、准确、高效的获取与处理环境信息的重要性，认识到在环境保护工作中广泛、灵活地应用信息技术的重要性。
 
目前，全国已有20多个省市自治区的环保部门及一百多个市地级环保部门开始使用地理信息系统软件和相应的硬件设施，其中大部分省市已经建立较为完备的环境基础数据库，并在GIS平台基础上开发了城市环境地理信息系统、重点流域水资源管理、环境污染应急预警预报系统等模块，根据目前的反馈来看，GIS平台的应用取得了显著的成效。
 
一、GIS在环境管理中的应用
 
GIS在环境管理中的介入，使得各种环境问题和环境过程描述更加符合实际，并促进了环境监测技术的发展。
 
在GIS和環境管理的实际结合中，可以将3S(GPS、RS、GIS)集成化运用。一是在汽车上装载GPS定位仪等设备，将车辆在移动中或是污染源现场采集的实时数据及时输送到GIS系统；二是针对固定采样点的信息采集，利用GPS将其坐标值导入GIS，然后对管理值进行实时监测并及时更新。
 
三是遥感(RS)技术在环境管理中也有着广泛的运用，环境研究人员通过运用遥感，得到准确的污染源分布、扩散条件、扩散范围等情况。环保部门通过将GIS与遥感技术相结合，对我国西部10多个省市自治区的生态环境状况进行调查，并从中得出了西部地区生态环境的空间分布和空间统计数据。
 
二、GIS在环境监测中的应用
 
环境监测离不开环境信息的采集和处理，而环境信息85%以上与空间位置有关，因此地理信息系统就成为环境监测的有效工具。2.1 GIS在大气环境动态监测中的应用大气环境有着以下特点：它的空间尺度大，人类赖以生存的大气圈有上百公里的厚度；空气在自然环境中有着最好的流动性，地面是其不可逾越的 固体边界 。因此大气环境动态监测最适合用GIS技术进行监测和分析。利用GIS空间分析和数据显示功能，可获得污染物在大气中的浓度分布图，进而可了解污染物的空间分布和超标情况。2.2 GIS在水资源环境监测中的应用 当今，水源污染日趋严重并多以复合型污染为基本特征，因此有必要加强对水资源环境的监测和管理。水资源环境的特点是空间信息量大，而对空间信息的管理与分析正是GIS的优点。GIS用于水资源环境监测，主要是对水质监测数据和空间数据进行科学有效的组织和管理，能够让管理人员方便地对各种空间信息进行查询、修改和编辑等；通过GIS强大的空间分析和图标分析功能，实现对空间和检测数据的分析和专题图的制作，进而为污染治理方案的制订提供有效的信息支持。
 
三、GIS在环境影响评价中的应用
 
环境影响评价是对环境质量优劣的定性定量评价，通过对所有的改、扩、建项目的环境数据进行分析、预测、评价，并给出保护环境的相应措施和对策。
 
GIS技术可以将数据和地图相连起来，建立拓扑关系，进行空间分析、查询以及制作各类专题图，形象、直观地显示环境质量和污染状况；3.1 GIS在项目环境影响评价中的应用利用GIS集成与场地和项目有关的各种数据及用于环境评价的各种模型，进行综合分析、模拟和预测，为环境质量现状进行分析和决策。GIS还具有很强的数据管理、更新和跟踪能力，以此来协助检查和监督环境影响评价单位和工程建设单位履行各自职责，并对环境影响报告书进行事后验证。3.2 GIS 在区域环境影响评价中的应用GIS在区域环境影响评价的作用主要体现在：能有效地管理一个大的地理区域复杂的污染源信息、环境质量信息及其他有关方面的信息，并能统计、分析区域环境影响的变化情况及主要污染源和主要污染物的地理属性和特征等；将区域的污染源数据库和环境特征数据库(如地形、气象等)与各种环境预测模型相关联，采用模型 预测法对区域的环境质量进行预测。
 
四、GIS在环境应急预警中的应用
 
建立环境应急预警，能够对事故风险源的地理位置及其属性、事故敏感区域位置及其属性进行管理，提供污染事故的大气、河流污染扩散的模拟过程和应急方案。GIS在环境应急系统中的作用主要是提供足够空间数据和相关属性数据的快速存取和管理功能，切实提高环境工作人员的管理决策能力。同时，GIS具有提供智能分析及策略的功能。例如，GIS可以通过对污染源位置、污染严重程度进行综合分析，将准确并具有时效性的数据导入到相应的污染应急模型中，便可以生成对应的应急方案，为决策者和应急处置人员提供客观且详实的预案。
 
随着GIS技术的不断发展，相信它在环境领域的运用会进一步得到扩展和完善， GIS系统在环境保护中发挥着越来越重要的作用，它与各种环境管理、分析系统进行有效结合，为环境工作者提供重要的技术指导和建议，有效地提高了环境工作效率，降低工作难度，使环境保护工作更加现代化、科学化。GIS的应用必将为我国的环境保护领域做出巨大的贡献。