1. 拼接算法现状#首先先需要说明的是不同于我们之前介绍的 CV 和 NLP在这些领域中DL 方法已经成为绝对主流传统机器学习和统计学方法更多只作为思想上的参考。其根本原因还是我们之前提到过的DL 本质上是数据驱动的方法而这些领域恰恰拥有海量数据从而使模型能够依靠规模持续提升性能。而在高光谱拼接这个细分领域中情况有所不同可以说传统方法仍然非常强势但与此同时也已经出现了大量基于 DL 的相关模型与技术可以预见未来 DL 方法的占比一定会持续提升。造成这一现状的原因大致可以归结为以下几点高光谱图像本身采集成本高、设备昂贵、波段标准不统一因此公开数据集相对稀缺。拼接问题往往不存在唯一严格正确的结果导致高质量监督标签难以构建。拼接问题本质上仍然属于几何配准与优化问题而传统数学方法在这一领域已经发展得非常成熟。高光谱数据与传感器、光照、成像系统等物理过程密切相关因此物理建模天然具有较强适用性。因此虽然在后续我们会逐步引入 DL 相关方法但这不代表近几篇的数学方法过时事实上即使是在这个领域的最新论文中依然能够频繁看到传统数学算法的身影。2. Harris 角点探测Harris Corner Detector#Harris 起源于 1988 年的论文 A Combined Corner and Edge Detector它可以说影响了整个局部特征检测时代。虽然现在的论文里看不到这种原始方法了但数学不能跳过加减乘除直接去学积分我们可以在后续大部分相关算法中都看到它的思想。简单来说我们可以把图像内容分成下面几类平坦区域往任何方向移动内容几乎没变化比如背景里的天空。边缘沿着边缘方向移动变化不大但垂直边缘方向移动变化很大比如墙的边。角点从这个位置往任何方向移动变化都很大比如窗户角物理上的不是“角所在的那个像素点”。Harris 的检测目标就是第三种各方向都敏感的位置。因此要强调的是Harris 不是边缘检测算法也不是“特殊像素点”检测算法而是局部特征检测算法。这种局部特征检测逻辑在后续会作为拼接算法的一环来找到不同图像中的相同物体。2.1 图像数据中的近似偏导#在展开具体逻辑前首先我们知道要求某个变量发生微小变化带来的变化量就要去求偏导limℎ→0(ℎ,)−(,)ℎ也就是说如果一个位置在多个方向上都存在较大导数值那么它更可能是角点。但这里有一个问题标准求导需要的是连续函数但我们现在只有离散的图像数据以灰度图为例[1234745235891625173286419374195582146693113877]显然我们没法直接计算每个元素在两个方向的导数值因此一种最基础的方法就是用差分法来近似导数比如以 方向为例前向差分≈(1)−()后向差分≈()−(−1)中心差分最常用≈(1)−(−1)2而 Harris 里常用的是信号处理中最基础的Sobel 算子同样以 方向为例我们 之前 用过它来进行边缘检测。[−101−202−101]你可能会觉得怎么突然从公式跳到矩阵了但我们把这个矩阵对应的运算展开你会发现这就是用加入平滑的差分法来近似导数对应图像里的一个 3×3 小块[111213212223313233]其运算就是(−1)11012113(−2)21022223(−1)31032133展开逻辑中间列全零表示不参与左右对比。-2 和 2 是权重表示中心信息更重要。上下两行的元素用来进行平滑降低噪声。值得一提的是如果是 RGB 图像可以对三个通道分别求各方向的近似偏导再融合也可以直接对三通道值加权融合后再计算近似偏导方法并不唯一更多通道也是同理。最终你会发现这个公式的语义是右边强度左边强度右边强度左边强度≈(右边强度)−(左边强度)结果 就表示这个 3×3 小区域在 方向上的变化强度。在这个例子里它就表示22 位置在 方向上的近似偏导。现在对整幅灰度图像分别应用两个方向的 Sobel 算子我们就得到了每个位置在 方向和 方向上的近似偏导即梯度 和 。这里还有一个细节在图像的边界卷积核会“伸出图像外”。因此需要对这一圈数据进行边界处理如抛弃、零填充、边缘复制或镜像填充镜像填充是实际工程中最常见的。2.2 窗口变化函数#我们已经得到了图像中每个位置的梯度 和 下一步问题是如何利用这些梯度判断一个区域是不是角点这里就引出了 Harris 的经典思想不看单个像素而看一个局部窗口在发生微小位移后整体变化有多大。原因首先是单个点的噪声因素过多没有统计稳定性并且很难区分边缘和角点。并且“角点”本身就是应该是局部结构因为物理意义上的角点几乎不可能只占一个像素。最终Harris 的结论是真正稳定的角点特征应该在一个小区域内都体现出明显变化。因此它定义了一个“窗口变化函数”(,)∑(,)∈(,)[(,)−(,)]2具体展开局部窗口在具体实现里尺寸在 3∼21 的范围内都有。(,)窗口内偏移注意是相对坐标不是绝对坐标。(,)整个窗口的平移。(,)权重通常是二维高斯权重。平方项避免正负抵消并强化变化明显的位置。(,) 表示这个窗口移动后与原窗口的差异大小。这里简单补充一下二维高斯权重(,)122exp⁡(−2222)在这里其逻辑就是让“中心像素更重要离中心越远越不重要”。这样会削弱边缘噪声同时不让微小平移过多影响输出。式子里的 是超参数分析就会发现小 权重集中在中心只看“非常局部”对细节敏感。大 权重分布更“平”、边缘部分的权重增加、更稳定但更粗糙。同时在具体工程实现中局部窗口大小 和 设计通常是≈61因为高斯在 ±3 外几乎为 0过大的窗口没有意义过小的窗口又会丢失信息。现在根据对同一个窗口应用不同偏移量的函数结果我们就可以对窗口内容进行区分如果窗口怎么移动都差不多平坦区域。如果只在某个方向变化大边缘。如果任何方向变化都大角点。这里有一个关键问题(,)对于这种“偏移后函数值”的形式我们当然可以规定偏移量为整数然后对每个窗口暴力枚举不同方向的响应值但这种方法无论是在精度还是计算量都十分局限。所以我们要回归数学原理Harris 的下一步是对这个式子进行一阶泰勒展开(,)≈(,)现在我们上一步计算的梯度就再次出现了再代回原式(,)∑(,)∈(,)[]2展开平方(,)∑(,)(22222)到这里我们就使用梯度把图像局部变化写成关于位移 (,) 的解析函数但是到这还没结束我们还要继续推导。2.3 Harris 矩阵#现在 我们进行移项把整个式子整理为矩阵运算形式(,)[][]其中[∑2∑∑∑2]这个矩阵就叫Harris 矩阵它是整个算法的核心∑2 表示区域内水平方向变化强度。∑2 表示区域内竖直方向变化强度。∑ 表示两个方向变化的相关性。你会发现这个式子的形式和之前提到过的协方差矩阵几乎完全一样但它统计的不是数据分布而是“这个局部区域在各方向上的整体梯度分布”。因此Harris 矩阵可以代表一个局部区域里的变化方向结构。而具有这种“局部结构描述工具”性质的矩阵后续一般统称为结构张量Structure Tensor。现在我们知道了通过Harris 矩阵就可以分析窗口在各个方向的变化强度接下来就进入经典内容了2.4 特征值分解和结果分析#因为矩阵 是一个实对称矩阵[]所以它一定可以进行特征值分解。设其两个特征值为:1,2它们的物理意义就是这个局部区域变化最剧烈的两个主方向的变化强度。展开一下逻辑在这里要求窗口在哪个方向变化最大我们要看 (,)[][]也就说我们其实在求相同长度位移得到的最大响应值max||1这部分的推导我们之前在 PCA 里已经展开过了。总之对于对称矩阵其在单位球上的最大值等于最大特征值。而取得最大值的方向就是对应特征向量。最终我们可以以此判断这个窗口的内容区域类型特征值关系局部变化特性几何含义典型例子平坦区域1≈0, 2≈0任意方向移动窗口变化都很小没有明显梯度与方向结构天空、纯色墙面边缘区域1≫2一个方向变化大另一个方向变化小仅存在单一主变化方向墙边、纸张边缘、道路边界角点区域1≫0, 2≫0任意方向移动都会产生明显变化两个主方向都存在强变化窗角、棋盘格交点、纹理交汇处2.5 工程优化和非极大值抑制#还没结束实际计算中我们不可能真的对每个矩阵都求特征值这会带来过大的计算量。因此Harris 又做一步工程优化根据齐次方程存在非零解时行列式为 0 和韦达定理对于矩阵 :[]其两个特征值满足12trace()12det()−2于是 Harris 构造了一个新的响应函数det()−⋅(trace())2其中det() 希望两个特征值都大。(trace())2 会惩罚“单方向特别大”的情况。 是超参数通常 ∈[0.04,0.06]。所以Harris 的最终判定逻辑就是区域类型Harris 响应 平坦区域≈0边缘区域0角点区域≫0最终在代码里角点判定条件就是阈值阈值阈值但到这里还还还还没有结束因为实际应用时会出现一个问题真正的角点附近往往会出现一大片连续高响应区域。原因就是我们最开始说的物理上的角点几乎不可能只在一个窗口或者一个像素上。因此在完成检测后还会再进行非极大值抑制再输出我们之前在检测算法那里展开过这里的基本逻辑也是相同的对于每个候选点取一个局域判断自己是不是其中最大的响应值如果不是就删除。到这里Harris 的完成逻辑才终于结束了。代码流程图如下Harris 虽然计算高效、旋转稳定并首次系统性地利用梯度统计描述局部结构但它仍存在明显局限固定窗口导致缺乏尺度不变性只能检测“哪里是角点”却无法描述角点本身同时对光照变化与噪声较为敏感。