计算机视觉课程5图像配准

1、模式识别计算机视觉图像配准申抒含中国自动化计算机视觉课程结构图12. 识别13. 人脑工程8-11.三维重建6. 运动估计7. 目标跟踪5.图像配准4.图像分割2. 特征检测3. 特征匹配图像什么是图像配准图像配准(Image registration)就是将不同时间、不同传感器或不同条件下获取的两幅或多幅图像进行匹配、叠加的过程。图像配准的种类刚体配准非刚体配准多模态配准图像配准的应用PanoramaBrown and Lowe, 2007图像配准的应用PhotomontageAgarwala, Dontcheva, Agrawala. 2004图像配准的应用Poisson image ed

2、itingPerez, Gangnet, and Blake, 2003图像配准的基本思想图像变换W：22源图像S目标图像T寻找变换W(u,v)，将源图像S一个变换图像，使其与目标图像T一致。图像配准的基本思想图像逆变换W-1： 22源图像S目标图像T为使变换图像完整，可以通过逆变换W-1将每一个像素回源图像S来构造变换图像。注意：并非所变换W都可逆图像配准中的变换模型平移（2D Translation）éu' ùé1ù éuù010txú ê úêúêêv

3、' ú = ê0ty ú êv ú1 úû êë1 úûê1 úêë0ëû图像点的齐次坐标度：2结构保持性：姿态图像配准中的变换模型欧式（2D Euclidean）éu' ùécosq-sinq cosq0ù éuùtxú ê úêúêêv' ú = ês

4、inqty ú êv úê1 úêë1 úû êë1 úû0ëû度：3结构保持性：长度图像配准中的变换模型相似（2D Similarity）éu' ùés×cosq-s×sinqs×cosq 0ù éuùtxú ê úêúêêv' ú = ês×

5、;sinqty ú êv ú1 úû êë1 úûê1 úêë0ëû度：4结构保持性：夹角图像配准中的变换模型仿射（2D Affine）éu' ùéaù éuùaa000102ú êêúêúv=a'aaú êv úêúê101112ê1 

6、50;êë 0úû êë1 úû01ëû度：6结构保持性：平行性图像配准中的变换模型摄影（2D Projective），也称（Homography）éu' ùé pù éuùpp000102ú êêúêúv=p'ppú êv úú ê1 úêúê101112ê1

7、úê pp1ë20û ë ûëû21度：8结构保持性：直线图像配准中的变换模型小孔相机模型éuùé f0f0ùêêëúúúû3ê 00úêêë 0ú1úû3´1图像配准中的变换模型éuùé f0ù0fêêëúúú

8、1;3ê 00úêêë 0ú1úû31é f0f00ùK = ê 00ú相机矩阵êêë 0ú1úû小孔相机模型由(K, R, C)完全决定（未考虑相机畸变参数）图像配准中的变换模型（Homography）R (C - C ) nTH = K (R R+-1-1 212nT X) K2211(X, n)éu' ùé pù éuùpp000102

9、0; êêúêúv=p'ppú êv úú ê1 úêúê101112ê1 úê pp1ëû ë ûëû2021HH(K2, R2, C2)(K1, R1, C1)度：8通过空间平面诱导的，图像配准中的变换模型（Homography）特殊形式R (C - C ) nT(X, n)H = K2 (R2R1+-1-1 212nT X) K1C1= C2R-1K-1

10、H = K R2211与空间平面(X, n)无关(K1, R1, C1)(K2, R2, C2)相机纯旋转的，度：3-5（根据相机焦距是否已知），常用于全景图拼接。图像配准中的变换模型非刚体变换通过网格形变进行非刚体配准度：任意变换模型的求解基于特征 q 变换W(. ;x) q22源图像S目标图像Tå(q«q')2e (x) =q '- W(q; x)变换模型的求解基于像素源图像S目标图像T变换W(. ;x)22e (x) = å T(q) -S(W-1(q; x)2qÎSSSD (Sum of Square Differences)变换

11、图像变换模型的求解基于特征点 VS基于像素基于特征点基于像素优点 变换模型计算只依赖特征点位置，对图像外观变化有一定鲁棒性； 能够计算宽基线时的变换； 使用了所有图像数据； 求解精度通常更高； 不需要特征检测和匹配；缺点 没有使用所有图像数据； 需要进行特征检测和匹配。 容易受到图像外观变化影响； 一般只适用于窄基线时的变换。变换模型的求解线性最小二乘如何求解(tx, ty)？(tx, ty)变换模型的求解线性最小二乘 (u1, v1) (u1, v1) éu'i ùtx ù éui ùé1010- u ù0

12、9; étùéu'é1êv' ú = ê0ú êv úxú =y ûiitêêv'úê01úêi úy ú êúêêë0i- vtëû ëëii ûêë1 úû1 úû êë1 ú

13、51;方程数量通常要远多于未知数数量， (tx, ty)为最小二乘解。变换模型的求解线性最小二乘x=(tx, ty) 为如下线性最小二乘的解：2éu'1 - u1ùúúúé10ù1úêv' - vê0êêêêú ét11ùú - êxú êminú ëty ûêú2tÎêu'n - un 

14、50;ê1êë00ú1úûêëv'n - vnúû2- bminAx2n´22´12n´12xÎ变换模型的求解线性最小二乘2minAx- bx为线性最小二乘的解：2xÎ2定义目标函数： e (x) =Ax- b22+- 2 xT AT bAxb¶e最优解x应满足：Ax = b= 2 AT Ax- 2 AT b = 0¶ x其中：A = (AT A)-1 AT因此：x = Ab称为A的伪逆变换模型的求解非线性最小二乘m&

15、#229;i=12Tmina x- bi线性最小二乘：inxÎx的线性函数må2minf (x) - b非线性最小二乘：iinxÎi=1x的非线性函数变换模型的求解非线性最小二乘mmin å2®,i = 1,., mn1fi :fi (x) - biwithnxÎi=1写成紧凑形式：2®minF(x) - bnmF: withnxÎæöæö÷÷F(x) = ç f(x) ÷= ç bbç÷ç(m&

16、#180;1)jjç÷ç÷èøèø变换模型的求解非线性最小二乘min e (x)2e (x) =F(x) - bnxÎ求解非线性最小二乘的：非线性最小二乘是非线性优化的一类特殊形式；一般非线性优化的梯度下降法（1阶）、阶）；（2、LMevenberg-非线性最小二乘的Marquardt）。变换模型的求解非线性最小二乘梯度下降法（Gradient Descent）将(x)在x处一阶Taylor展开：e (x+) = e (x) + g(x)T +o(2 )g(x) = Ñe (x)则当|0时，(x

17、) 的变化值可以表示为：De = e (x+) - e (x)= gT cosq=g为x处梯度与向量间的夹角。显然， =时(x)下降最快。变换模型的求解非线性最小二乘梯度下降法（Gradient Descent）梯度下降法的基本思想是将自变量向负梯度方向移动：x ¬ x-a g(x)为移动步长，最优的可以通过一维搜索获得：min e (x-a g(x)a ³0寻找最优的 计算代价很高，因此常使用一些简单的近似算法，如折半查找法（从 =1开始不断将 折半）。变换模型的求解非线性最小二乘梯度下降法（Gradient Descent）2： e (x) =F(x) - b对于最小二

18、乘梯度方向为： g(x) = J(x)T r(x)Gradient Descent方向为：(x) = -J(x)T r(x)= ¶ F(x)观测状态的Jacobian矩阵J(x)其中：m´n¶ xr(x) = F(x) - b残差向量变换模型的求解非线性最小二乘（Newton Method）梯度下降法是一阶收敛算法，收敛较慢。为加快收敛速度，可将(x)在x处Taylor展开：1e (x+) » e (x) + g(x) +H(x)TT2g(x) 为n×1的梯度向量，H(x) 为n×n的Hessian矩阵。将(x+)对求导，并令导数为0

19、，可得使(x+)最小化的：g(x)+ H(x)=0(x)= - H(x)-1 g(x)变换模型的求解非线性最小二乘（Newton Method）2： e (x) =F(x) - b对于最小二乘梯度向量为：g(x) = J(x)Tr(x)Hessian矩阵为：H(x) = J(x)TJ(x) + S(x)J(x) = ¶ F(x)¶ xr(x) = F(x) - bS(x)包含(x)的导数变换模型的求解非线性最小二乘（Newton Method）2： e (x) =F(x) - b对于最小二乘Newton方向为： (x) = -(J(x)TJ(x) + S(x)-1 J(x)

20、T r(x)是收敛算法，收敛速度比梯度下降法快；需要求Hesian矩阵，即目标函数的导数。变换模型的求解非线性最小二乘（Gauss-Newton Method）最小二乘的一种有效解法是，避免了中对Hessian矩阵的计算。2： e (x) =F(x) - b对于最小二乘将F(x)在x处一阶Taylor展开：F(x+) » F(x) + J(x)J(x) = ÑF(x)将上式代入(x)可得：e (x+ ) »F(x) + J(x) -b 2变换模型的求解非线性最小二乘（Gauss-Newton Method）F(x) + J(x) -b 2e (x+ ) »

21、;r(x) = F(x) - b由可得：2e (x+ ) »r(x) + J(x)显然，如果要求(x+) 最小化，则r(x) + J(x) = 0则Gauss-Newton方向为：(x) = - J(x)：r(x)= -(J(x)T J(x)-1 J(x)Tr(x)变换模型的求解非线性最小二乘（Gauss-Newton Method）与的比较：(x) = -(J(x)T J(x) + S(x)-1 J(x)Tr(x)Newton方向为：Gauss-Newton方向为：(x) = -(J(x)TJ(x)-1 J(x)T r(x)S(x)可以忽略的条件：2中F(x)适度非线性e (x)

22、=F(x) - b最小二乘残差r(x)比较小大多数计算机视觉都满足这两个条件。变换模型的求解非线性最小二乘阻尼（Damped Gauss-Newton Method）的特点：的收敛速度比梯度下降法快的多（超线性收敛收敛） ；，近似不能保证收敛（梯度下降法和收敛）。：(x,) = -(J(x)T J(x) + I)-1 J(x)T r(x)阻尼lim (x,)则成为®0lim (x,)则成为梯度下降法®¥变换模型的求解非线性最小二乘Levenberg-Marquardt Method (LM法)LM法通过启发式在每一步动态调整：(x,) = -(J(x)T J(x)

23、 + I)-1 J(x)T r(x)如果误差减少，则令0.1如果误差增大，则令10LM法是一种启发式的阻尼泛使用。，在计算机视觉中广变换模型的求解非线性最小二乘一阶和收敛算法求解非线性最小二乘比较：多数情况下我们会使用阻尼（如LM法）来处理计算机视觉中的最小二乘。图像配准实例：全景图拼接M. Brown and D. G. Lowe, Automatic Panoramic Image Stitching using Invariant Features, International Journal of Computer Vision. 74(1): 59-73, 2007.纯旋转拍摄无序图

24、像全景图像图像配准实例：全景图拼接算法基本流程：1.2.3.4.5.6.特征点检测与匹配；图像匹配与；变换模型参数计算；全景图矫正； 增益补偿；多频段光度融合。图像配准实例：全景图拼接步骤1：特征点检测与匹配检测每幅图像中的SIFT特征点（共n幅图像）； 在每幅图像上建立k-d树；每幅图像中每个特征点与其他图像匹配（只和每幅图像k-d 树中4个最近邻特征点计算匹配值）。回忆一下k-d树：个非线性最小二乘，为了简化操作，放松为一个普（ 8通度）。(K1, R1, C)(K2, R2, C)图像配准实例：全景图拼接步骤2：图像匹配与相机纯旋转拍摄，因此图像间变换为一个5度的变换（一般变换为8度）H

25、 = K R R-1K-122115度的变换图像配准实例：全景图拼接步骤2：图像匹配与相机纯旋转拍摄，因此图像间变换为一个5度的变换变换为8（一般度）éu' ùé pù éuùpp000102ú êêúêúv=p'ppú êv úú ê1 úêúê101112ê1 úê pp1ëû ë ûë&

26、#251;2021HH包含8个未知量，可以通过4组对应点求得。(K1, R1, C)(K , R , C)22图像配准实例：全景图拼接步骤2：图像匹配与H，只需选择对每一幅图像，不需要和其他所有图像计算匹配点最多的6幅图像进行计算。匹配点中不可避免的点，因此使用RANSAC估计H。外RANSAC计算H流程：1.2.3.4.随机选4对匹配点8个线性方程组求H；计算H的一致集（所有符合H的匹配点步骤1-2循环500次；的集合）；最大一致集中的所有匹配点性最小二乘法求解H。超定线性方程组，通过线图像配准实例：全景图拼接步骤2：图像匹配与假定两幅图像中正确匹配点所占比例为p=0.5，则经过n次RANS

27、AC后找到正确里r=4）：H的概率为（r为最小数据点个数，这p(H is correct) = 1- (1- ( p)r )n当n=500时，正确的H未被找到的概率为1×10-14图像配准实例：全景图拼接步骤2：图像匹配与变换H，下面对每一对匹配图像，我们获得了它们之间的需要这对图像是否确实匹配？的基本依据：源图像经过变换后与目标图像产生一定面积的重叠区域，记这一区域内特征点总数为nf，记这一区域中符合H的特征点（内点）数量为ni，如果ni /nf足够大，则这两幅图像确实匹配。图像配准实例：全景图拼接步骤2：图像匹配与定义指标m，m=1表示图像匹配，m=0表示图像不匹配； 定义指标f

28、 (i)，f (i) =1表示第i个特征点为内点，f (i) =0表示第i个特征点为外点；分布（Bernoulli distribution）：则总内点数符合nf !-n )(1:n )(1- p )(n| m = 1) = B(n ; n, p ) =np( fpffiiif111n !(n- n )!ifnf !i-n )(1:n )(1- p )(n| m = 0) = B(n ; n , p ) =np( fpffiiif000n !(n- n )!ifip1和p0分别为图像匹配和不匹配时特征点为内点的概率图像配准实例：全景图拼接步骤2：图像匹配与根据公式，图像正确匹配的后验概率为：(

29、1:n f)| m = 1) p(m = 1)p( f(1:n )p(m = 1| f) =f(1:n f ) )1p( f=(1:n f )| m = 0) p(m = 0)p( f1+(1:n )| m = 1) p(m = 1)p( ff(1:nf) ) > p当后验概率 p(m = 1| f时认为两幅图像正确匹配。min图像配准实例：全景图拼接步骤2：图像匹配与等价于当下面不等式成立时，两幅图像正确匹配：B(ni ; nf , p1) p(m = 1)1>1B(ni ; nf , p0 ) p(m = 0)-1pmin给定参数：p1=0.6，p0=0.1，p(m=1)=10

30、-6，pmin=0.999，上述判定条件可以简化为：ni > a + b nf其中：=8.0，=0.3图像配准实例：全景图拼接步骤2：图像匹配与使用匹配条件剔除掉错误匹配图像对后，可以获得图像匹配的连接集，即输入图像可以分成几幅全景图。图像配准实例：全景图拼接步骤3：变换模型参数计算特征点匹配中使用RANSAC得到的H是根据两两图像计算获得的，为保证匹配点的全局一致性，在变换模型参数计调整算法（Bundle Adjustment）。算中使用相机纯旋转拍摄时对应一个5度的变换，令u k和u k为ij图像i和j中第k组匹配特征点，则将u k通过变换到图像i上j的位置为：= K R R-1K-

31、1 ukpkijiijjj变换误差为：= uk - pkrkijiij图像配准实例：全景图拼接步骤3：变换模型参数计算误差函数定义为所有图像中所有匹配点的变换误差平方和：ne = å ååh( rk )iji=1 jÎI (i) kÎF (i, j )其中I(i)表示所有与图像i匹配的图像，F(i,j)表示图像i与图像j 中的特征匹配点集合，h()为Huber误差函数(Huber, 1981)：ìï| x |2 ,if | x |< sh(x) =íss ,if | x |³ s| x | -2&#

32、239;î2在优化循环初期使用=，在后期使用=2。图像配准实例：全景图拼接步骤3：变换模型参数计算误差函数定义为所有图像中所有匹配点的变换误差平方和：ne = å ååh( rk )iji=1 jÎI (i) kÎF (i, j )误差函数的自变量包含每个相机的焦距和旋转矩阵，共4n个参数（n为图像数量），是一个非线性最小二乘Levenberg-Marquardt (LM)算法对其进行求解。，使用图像配准实例：全景图拼接步骤4：全景图矫正在上一步骤中，我们已经计算了每个相机的参数（焦距f和姿态R），下面需要考虑如何把它们拼接起来。平面

33、坐标：最简单的拼接是将某一副图像作为参考图像（如图像i），通过Hij=K R R -1K -1 (j=1,n)计算，并将图iijj像j (j=1,n)都转到图像i所在平面完成拼接。图像配准实例：全景图拼接步骤4：全景图矫正平面坐标拼接在源图像和目标图像间旋转角度较大时，图像会被严重拉伸，影响拼接效果。图像配准实例：全景图拼接步骤4：全景图矫正解决图像过度拉伸的是将图像变换到一个圆柱面上，这个圆柱中轴线通过相机光心，底面为圆。图像配准实例：全景图拼接步骤4：全景图矫正uu ' = s tan-1fvv ' = s+ fu22(u,v)：图像平面坐标(u,v)：柱面展开平面坐标柱面

34、坐标：使用柱面坐标系进行拼接时，每个图像点变换为圆柱面上一点：I(u,v)I(,h)I(u,v)柱面坐标适合纯水平方向旋转的全景图拼接，如果在竖直方向也有旋转，则可能导致竖直方向旋转角大的图像被过度拉伸。图像配准实例：全景图拼接步骤4：全景图矫正解决竖直方向图像过度拉伸的是将图像变换到一个单位球面上，这个球的球心和相机光心重合。图像配准实例：全景图拼接步骤4：全景图矫正uu ' = s tan-1fvv ' = s tan-1+ fu22(u,v)：图像平面坐标(u,v)：球面展开平面坐标球面坐标：使用球面坐标系进行拼接时，每个图像点变换为球面上一点：I(u,v)I(,)I(u

35、,v)在后续步骤中，全景图的拼接均使用球面坐标。图像配准实例：全景图拼接步骤4：全景图矫正图像变换到球面坐标系时，根据球面坐标系主轴方向不同，最终获得的全景图会有很大不同。图像配准实例：全景图拼接步骤4：全景图矫正图像变换到球面坐标系时，根据球面坐标系主轴方向不同，最终获得的全景图会有很大不同。图像配准实例：全景图拼接步骤4：全景图矫正图像变换到球面坐标系时，根据球面坐标系主轴方向不同，最终获得的全景图会有很大不同。因此，将图像变换到球面坐标系时需要确定合适的主轴方向，主轴方向与图像拍摄时的相机姿态直接相关。图像配准实例：全景图拼接步骤4：全景图矫正全景图拍摄时，相机可能发生水平旋转、竖直旋转

36、，但一般极少发生绕光轴旋转，因此相机X轴通常会保持在一个平面，这一平面的法向u可以作为一个理想的主轴方向。图像配准实例：全景图拼接步骤4：全景图矫正全景图拍摄时，相机可能发生水平旋转、竖直旋转，但一般极少发生绕光轴旋转，因此相机X轴通常会保持在一个平面，这一平面的法向u可以作为一个理想的主轴方向。我们期望得到的主轴方向u与所有图像坐标系X轴垂直，因此：u = 0XTi=1,ni写成最小二乘形式：nmin åi=1XTui图像配准实例：全景图拼接步骤4：全景图矫正，最优解u满足：这是一个线性最小二乘æöu = 0nÑånçå=

37、 0XTXTuX÷iiièøi=1i=1nu为3×3矩阵åXi=1XT最小特征值对应的特征向量。ii以u作为主轴方向确定球面坐标系，根据相机参数将所有图像变换到球面上I(u,v)I(,h)I(u,v)，获得全景图。图像配准实例：全景图拼接步骤4：全景图矫正未校正全景图校正后全景图图像配准实例：全景图拼接步骤5：增益补偿在前面的步骤中已经确定了全景图拼接的几何结构参数（相机焦距、相机姿态、全景图坐标），已经可以将图像拼接为全景图像，但这样直接拼接时可能产生光度不一致的现象。直接拼接结果图像配准实例：全景图拼接步骤5：增益补偿为保证拼接后重叠区域的

38、光度一致性，需要对每幅图像进行增益补偿。× g=图像配准实例：全景图拼接步骤5：增益补偿为使全景图全局光度一致，我们在所有图像上定义全局光度误差函数：(g I- g j I ji )2æç2 ö()1- gnne(g ,., g ) = 1 åå N÷iji+iç÷ss1nij2N2g2ç÷i=1j =1èøn为图像数量，gi为第i幅图像的增益，Nij为图像i和j重叠区域像素数量，Iij为图像i在图像i与j重叠区域的亮度均值，N为归一化 光度误差方差（设为10），g

39、为增益方差（设为0.1）。最小化e(.)是一个二次优化优增益g1gn，可通过超定线性方程组求得最图像配准实例：全景图拼接步骤5：增益补偿直接拼接结果增益补偿结果图像配准实例：全景图拼接步骤6：多频段光度融合由于光晕、光心移动、配准误差、镜头畸变等因素的影响，即使经过增益补偿后，全景图中可能某些光度不一致的区域，需要进一步进行。图像重叠区域可能包含多幅图像，一种有效的光度方式是将重叠区域像素用多幅图像进行融合。图像配准实例：全景图拼接步骤6：多频段光度融合令Ii(u,v)表示第i幅图像上的像素点(u,v)，Wi (u,v)=(u) (v)表示像素点(u,v)的融合权重，(.)的值在图像中心为1，在图像边缘为0，中心和边缘之间线性递减。将Ii(u,v)和Wi (u,v)均转换到球面坐标系Ii(,)和Wi (,)，则全景图像中每一个像素点可以通过下式进行光度融合：åni=1Ii (q ,f) Wi (q ,f)Ilinear (q ,f) =åni=1Wi (q ,f)这一融合称为线性光度融合。图像配准实例：全景图拼接步骤6：多频段光度融合如果在图像配准过程中出现微小误差，那么线性光度融合在图像高频部分可能出现模糊，因此一