單應性Homograph估計

單應性原理被廣泛應用於圖像配准，全景拼接，機器人定位SLAM，AR增強現實等領域。這篇文章從基礎圖像坐標知識系為起點，講解圖像變換與坐標系的關系，介紹單應性矩陣計算方法，並分析深度學習在單應性方向的進展。

圖像變換與平面坐標系的關系

旋轉

將圖形圍繞原點\((0,0)\)逆時針方向旋轉\(\theta\)角，用解析式表示為：

\[\begin{matrix}x'=x\cdot cos\theta -y\cdot sin\theta\\y'=x\cdot sin\theta+y\cdot cos \theta\end{matrix} \]

寫成矩陣乘法形式：

\[\begin{pmatrix}x'\\y'\end{pmatrix}=\begin{bmatrix}cos\theta&-sin\theta\\sin\theta&cos\theta\end{bmatrix}\begin{pmatrix}x\\y\end{pmatrix}=R\begin{pmatrix}x\\y\end{pmatrix} \]

平移

\[\begin{matrix}x'=x+t_x\\y'=y+t_y\end{matrix} \]

\[\begin{pmatrix}x'\\y'\end{pmatrix}=\begin{pmatrix}x\\y\end{pmatrix}+\begin{bmatrix}t_x\\t_y\end{bmatrix} \]

但是現在遇到困難了，平移無法寫成和上面旋轉一樣的矩陣乘法形式。所以引入齊次坐標 \((x,y)\Leftrightarrow\)，再寫成矩陣形式：

\[\begin{pmatrix}x'\\y'\\1\end{pmatrix}=\begin{bmatrix}1&0&t_x\\0&1&t_y\\0&0&1\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix}=\begin{bmatrix}I_{2\times 2}&T_{2\times 1}\\0^T&1\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix} \]

其中\(I_{2\times 2}=\begin{bmatrix}1&0\\0&1\end{bmatrix}\)表示單位矩陣，而\(T_{2\times 1}=\begin{bmatrix}t_x\\t_y\end{bmatrix}\)表示平移向量。

剛體變換

那么就可以把把旋轉和平移統一寫在一個矩陣乘法公式中，即剛體變換：

\[\begin{pmatrix}x'\\y'\\1\end{pmatrix}=\begin{bmatrix}cos\theta&-sin\theta&t_x\\sin\theta&cos\theta&t_y\\0&0&1\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix}=\begin{bmatrix}R_{2\times 2}&T_{2\times 1}\\0^T&1\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix} \]

而旋轉矩陣\(R_{2\times 2}\)是正交矩陣\((RR^T=R^TR=I)\)

仿射變換

\[\begin{pmatrix}x'\\y'\\1\end{pmatrix}=\begin{bmatrix}A_{2\times 2}&T_{2\times 1}\\0^T&1\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix} \]

其中\(A_{2\times 2}=\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}\)可以是任意\(2\times 2\)矩陣，而\(R\)必須是正交矩陣

可以看到，相比剛體變換（旋轉和平移），仿射變換除了改變目標位置，還改變目標的形狀，但是會保持物體的“平直性（如圖形中平行的兩條線變換后依然平行）”。

不同\(A\)和\(T\)矩陣對應的各種基本仿射變換：

投影變換（單應性變換）

\[\begin{pmatrix}x'\\y'\\1\end{pmatrix}=\begin{bmatrix}A_{2\times 2}&T_{2\times 1}\\V^T&s\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix}=H_{3\times 3}\begin{pmatrix}x\\y\\1\end{pmatrix} \]

簡單說，投影變換徹底改變目標的形狀。

總結

1. 剛體變換：平移+旋轉，只改變物體位置，不改變物體形狀
2. 仿射變換：改變物體位置和形狀，但是保持“平直性”（原來平行的邊依然平行）
3. 投影變換：徹底改變物體位置和形狀

我們來看看完整投影變換矩陣各個參數的物理含義：

\[H_{3\times 3}=\begin{bmatrix}A_{2\times 2}&T_{2\times 1}\\V^T&s\end{bmatrix}=\begin{bmatrix}a_{11}&a_{12}&t_x\\a_{21}&a_{22}&t_y\\v_1&v_2&s\end{bmatrix} \]

其中\(A_{2\times 2}\)代表仿射變換參數，\(T_{1\times 1}\)代表平移變換參數。

而\(V^T=[v_1,v_2]\)表示一種“變換后邊緣交點”關系，如：

至於\(s\)則是一個與\(V^T=[v_1.v_2]\)相關的縮放因子。

\[\begin{bmatrix}1&0&t_x\\0&1&t_y\\v_1&v_2&s\end{bmatrix}\begin{pmatrix}x\\y\\1\end{pmatrix}=\begin{pmatrix}x\\y\\v_1x+v_2y+s\end{pmatrix}\Leftrightarrow\begin{pmatrix}\frac{x}{v_1x+v_2y+s}\\\frac{y}{v_1x+v_2y+s}\end{pmatrix} \]

一般情況下都會通過歸一化使得\(s=1\)，原因見下文。

平面坐標系與齊次坐標系

問題來了，齊次坐標到底是什么？

齊次坐標系\((x,y,w)\in \mathbb{P}^3\)與常見的三維空間坐標系\((x,y,z)\in \mathbb{R}^3\)不同，只有兩個自由度：

\[\begin{pmatrix}\frac{x}{w}\\\frac{y}{w}\end{pmatrix}\Leftrightarrow \begin{pmatrix}x\\y\\w\end{pmatrix} \]

而\(w\)（其中\(w>0\)）對應坐標\(x\)和\(y\)的縮放尺度。當\(w=1\)時：

\[\begin{pmatrix}x\\y\end{pmatrix}\Leftrightarrow \begin{pmatrix}x\\y\\1\end{pmatrix} \]

特別的當\(w=0\)時，對應無窮遠：

\[\begin{pmatrix}\text{Inf}\\\text{Inf}\end{pmatrix}\Leftrightarrow \begin{pmatrix}x\\y\\1\end{pmatrix} \]

從二維平面上看，\((x,y,w)\)隨\(w\)的變化在從原點到\((x,y)\)的藍虛線示意的射線上滑動：

單應性變換

單應性是什么？

此處給出單應性不嚴謹的定義：用 [無鏡頭畸變] 的相機從不同位置拍攝 [同一平面物體] 的圖像之間存在單應性，可以用 [投影變換] 表示 。

單應變換就是投影變換，逆投影變換是投影變換的逆

注意：
單應性的嚴格定義與成立條件非常復雜，超出本文范圍。
有需要的朋友請自行查閱《Multiple View Geometry in Computer Vision》書中相關內容。

簡單說就是：

\[\begin{pmatrix}x_l\\y_l\\1\end{pmatrix}=H_{3\times3}\times\begin{pmatrix}x_r\\y_r\\1\end{pmatrix} \]

其中\((x_l,y_l)\)是Left view圖片上的點，\((x_r,y_r)\)是Right view圖片上對應的點。

那么這個\(H_{3\times3}\)單應性矩陣如何求解呢？

從更一般的情況分析，每一組匹配點\((x_i,y_i)\overset{match}{\longrightarrow}(x_i',y_i')\)有以下等式成立：

\[\begin{pmatrix}x_i'\\y_i'\\1\end{pmatrix}=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&h_{33}\end{bmatrix}\begin{pmatrix}x_i\\y_i\\1\end{pmatrix}=\begin{pmatrix}h_{11}x_i+h_{12}y_i+h_{13}\\h_{21}x_i+h_{22}y_i+h_{23}\\h_{31}x_i+h_{32}y_i+h_{33}\end{pmatrix} \]

由平面坐標與齊次坐標對應關系\((\frac{x}{w},\frac{y}{w})\in \mathbb{R}^2\Leftrightarrow (x,y,w)\in \mathbb{P}^3\)，上式可以表示為：

\(\begin{matrix}x'_i=\frac{h_{11}x_i+h_{12}y_i+h_{13}}{h_{31}x_i+h_{32}y_i+h_{33}}\\y'_i=\frac{h_{21}x_i+h_{22}y_i+h_{23}}{h_{31}x_i+h_{32}y_i+h_{33}}\end{matrix}\)

進一步變換為：

\[\begin{matrix}(h_{31}x_i+h_{32}y_i+h_{33})\cdot x'_i=h_{11}x_i+h_{12}y_i+h_{13}\\(h_{31}x_i+h_{32}y_i+h_{33})\cdot y'_i=h_{21}x_i+h_{22}y_i+h_{23}\end{matrix} \]

寫成矩陣\(AX=0\)形式：

\[\begin{bmatrix}x_i&y_i&1&0&0&0&-x_i' x_i&-x_i' y_i&-x_i'\\0&0&0&x_i&y_i&1&-y_i'x_i&-y_i'y_i&-y_i'\end{bmatrix}\begin{bmatrix}h_{11}\\h_{12}\\h_{13}\\h_{21}\\h_{22}\\h_{23}\\h_{31}\\h_{32}\\h_{33}\end{bmatrix}=0 \]

也就是說一組匹配點\((x_i,y_i)\overset{match}{\longrightarrow}(x_i',y_i')\)可以獲得2組方程。

單應性矩陣8自由度

注意觀察：單應性矩陣\(H\)與\(aH\)其實是完全一樣的（其中$a\ne 0 $），例如：

\[\begin{pmatrix}x_i'\\y_i'\\1\end{pmatrix}=aH\begin{pmatrix}x_i\\y_i\\1\end{pmatrix} \]

\[\begin{matrix}x'_i=\frac{ah_{11}x_i+ah_{12}y_i+ah_{13}}{ah_{31}x_i+ah_{32}y_i+ah_{33}}=\frac{h_{11}x_i+h_{12}y_i+h_{13}}{h_{31}x_i+h_{32}y_i+h_{33}}\\y'_i=\frac{ah_{21}x_i+ah_{22}y_i+ah_{23}}{ah_{31}x_i+ah_{32}y_i+ah_{33}}=\frac{h_{21}x_i+h_{22}y_i+h_{23}}{h_{31}x_i+h_{32}y_i+h_{33}}\end{matrix} \]

即點\((x_i,y_i)\)無論經過\(H\)還是\(aH\)映射，變化后都是\((x_i',y_i')\)。

如果使\(a=\frac{1}{h_{33}}\)，那么有：

\[H'=aH=\begin{bmatrix}\frac{h_{11}}{h_{33}}&\frac{h_{12}}{h_{33}}&\frac{h_{13}}{h_{33}}\\\frac{h_{21}}{h_{33}}&\frac{h_{22}}{h_{33}}&\frac{h_{23}}{h_{33}}\\\frac{h_{31}}{h_{33}}&\frac{h_{32}}{h_{33}}&1\end{bmatrix}=\begin{bmatrix}h_{11}'&h_{12}'&h_{13}'\\h_{21}'&h_{22}'&h_{23}'\\h_{31}'&h_{32}'&1\end{bmatrix} \]

所以單應性矩陣\(H\)雖然有9個未知數，但是只有8個自由度。

在求\(H\)時一般添加約束\(h_{33}=1\)（也有用\(\sqrt{h_{11}^2+h_{12}^2+...+h_{33}^2}=1\)約束），所以還有\(h_{11}\sim h_{32}\)共8個未知數。由於一組匹配點\((x_i,y_i)\overset{match}{\longrightarrow}(x_i',y_i')\)對應2組方程，那么只需要\(n=4\)組不共線的匹配點即可求解\(H\)的唯一解。

下圖為xiaomi9拍攝，有鏡頭畸變：

OpenCV已經提供了相關API，代碼和變換結果如下。

import cv2
import numpy as np

im1 = cv2.imread('left.jpg')
im2 = cv2.imread('right.jpg')

src_points = np.array([[581, 297], [1053, 173], [1041, 895], [558, 827]])
dst_points = np.array([[571, 257], [963, 333], [965, 801], [557, 827]])

H, _ = cv2.findHomography(src_points, dst_points)

h, w = im2.shape[:2]

im2_warp = cv2.warpPerspective(im2, H, (w, h))

可以看到：

1. 紅框所在平面上內容基本對齊，但受到鏡頭畸變影響無法完全對齊；
2. 平面外背景物體不符合單應性原理，偏離很大，完全無法對齊。

傳統方法估計單應性矩陣

一般傳統方法估計單應性變換矩陣，需要經過以下4個步驟：

1. 提取每張圖SIFT/SURF/FAST/ORB等特征點
2. 提取每個特征點對應的描述子
3. 通過匹配特征點描述子，找到兩張圖中匹配的特征點對（這里可能存在錯誤匹配）
4. 使用RANSAC算法剔除錯誤匹配
5. 求解方程組，計算Homograph單應性變換矩陣

示例代碼如下：

#coding:utf-8

# This code only tested in OpenCV 3.4.2!
import cv2 
import numpy as np

# 讀取圖片
im1 = cv2.imread('left.jpg')
im2 = cv2.imread('right.jpg')

# 計算SURF特征點和對應的描述子，kp存儲特征點坐標，des存儲對應描述子
surf = cv2.xfeatures2d.SURF_create()
kp1, des1 = surf.detectAndCompute(im1, None)
kp2, des2 = surf.detectAndCompute(im2, None)

# 匹配特征點描述子
bf = cv2.BFMatcher()
matches = bf.knnMatch(des1, des2, k=2)

# 提取匹配較好的特征點
good = []
for m,n in matches:
    if m.distance < 0.7*n.distance:
        good.append(m)

# 通過特征點坐標計算單應性矩陣H
# （findHomography中使用了RANSAC算法剔初錯誤匹配）
src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1,1,2)
dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1,1,2)
H, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)
matchesMask = mask.ravel().tolist()

# 使用單應性矩陣計算變換結果並繪圖
h, w, d = im1.shape
pts = np.float32([[0,0], [0,h-1], [w-1,h-1], [w-1,0]]).reshape(-1,1,2)
dst = cv2.perspectiveTransform(pts, H)
img2 = cv2.polylines(im2, [np.int32(dst)], True, 255, 3, cv2.LINE_AA)

draw_params = dict(matchColor = (0,255,0), # draw matches in green color
                   singlePointColor = None,
                   matchesMask = matchesMask, # draw only inliers
                   flags = 2)

im3 = cv2.drawMatches(im1, kp1, im2, kp2, good, None, **draw_params)

相關內容網上資料較多，這里不再重復造輪子。需要說明，一般情況計算出的匹配的特征點對\((x_i,y_i)\overset{match}{\longrightarrow}(x_i',y_i')\)數量都有\(n \gg 4\)，此時需要解超定方程組（類似於求解線性回歸時數據點的數量遠多於未知數）。

深度學習在單應性方向的進展

HomographyNet（深度學習end2end估計單應性變換矩陣）

HomographyNet是發表在CVPR 2016的一種用深度學習計算單應性變換的網絡，即輸入兩張圖，直接輸出單應性矩陣\(H\)。

在之前的分析中提到，只要有4組\((x_i,y_i)\overset{match}{\longrightarrow}(x_i',y_i')\)匹配點即可計算\(H_{3\times 3}\)的唯一解。

相似的，只要有4組\((\bigtriangleup x_i,\bigtriangleup y_i)\)也可以計算出\(H_{3\times 3}\)的唯一解：

\[\begin{cases}(\bigtriangleup x_1,\bigtriangleup y_1)\\(\bigtriangleup x_2,\bigtriangleup y_2)\\(\bigtriangleup x_3,\bigtriangleup y_3)\\(\bigtriangleup x_4,\bigtriangleup y_4)\end{cases}\rightarrow H=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&1\end{bmatrix} \]

其中\(\bigtriangleup x_i=x_i-x_i'\)且\(\bigtriangleup y_i=y_i-y_i'\)

分析到這里，如果要計算\(H\)​，網絡輸出可以有以下2種情況：

1. Regression：網絡直接輸出\((\bigtriangleup x_1,\bigtriangleup y_1)\sim (\bigtriangleup x_4,\bigtriangleup y_4)\)共8個數值

   這樣設置網絡非常直觀，使用L2損失訓練，測試時直接輸出8個float values，但是沒有置信度confidence。即在使用網絡時，無法知道當前輸出單應性可靠程度。

2. Classification：網絡輸出\((\bigtriangleup x_1,\bigtriangleup y_1)\sim (\bigtriangleup x_4,\bigtriangleup y_4)\)共8個數值的量化值+confidence

   這時將網絡輸出每個\(\bigtriangleup x_i\)和\(\bigtriangleup y_i\)​量化成21個區間，用分類的方法判斷落在哪一個區間。訓練時使用Softmax損失。相比回歸直接輸出數值，量化必然會產生誤差，但是能夠輸出分類置信度評判當前效果好壞，更便於實際應用。

另外HomographyNet訓練時數據生成方式也非常有特色。

1. 首先在隨機\(p\)​位置獲取正方形圖像塊Patch A
2. 然后對正方形4個點進行隨機擾動，同時獲得4組\((\bigtriangleup x_i,\bigtriangleup y_i)\)
3. 再通過4組\((\bigtriangleup x_i,\bigtriangleup y_i)\)計算\(H^{AB}\)
4. 最后將圖像通過\(H^{AB}=(H^{AB})^{-1}\)變換，在變換后圖像\(p\)位置獲取正方形圖像塊Patch B

那么圖像塊A和圖像塊B作為輸入，4組\((\bigtriangleup x_i,\bigtriangleup y_i)\)​作為監督Label，進行訓練

可以看到，在無法提取足夠特征點的弱紋理區域，HomographyNet相比傳統方法確實有一定的優勢：

Spatial Transformer Networks（直接對CNN中的卷積特征進行變換）

其實早在2015年，就已經有對CNN中的特征進行變換的STN結構。

假設有特征層\(U\)，經過卷積變為\(V\)，可以在他們之間插入STN結構。這樣就可以直接學習到從特征\(U\)上的點\((\bigtriangleup x_i^u,\bigtriangleup y_i^u)\)的仿射變換。

\[\begin{pmatrix}x_i^u\\y_i^u\end{pmatrix}=T_\theta(G_i)=A_\theta\cdot \begin{pmatrix}x_i^u\\y_i^u\\1\end{pmatrix}=\begin{bmatrix}\theta_{11}&\theta_{12}&\theta_{13}\\\theta_{21}&\theta_{22}&\theta_{23}\end{bmatrix}\begin{pmatrix}x_i^u\\y_i^u\\1\end{pmatrix} \]

其中\(A_\theta\)對應STN中的仿射變換參數。STN直接在特征維度進行變換，且可以插入輕松任意兩層卷積中。

DELF: DEep Local Features（深度學習提取特征點與描述子）

之前提到傳統方法使用SIFT和Surf等特征點估計單應性。顯然單應性最終估計准確度嚴重依賴於特征點和描述子性能。Google在ICCV 2017提出使用使用深度學習提取特征點。

考慮到篇幅，這里不再展開DELF，請有興趣的讀者自行了解相關內容。

關於OpenCV圖像坐標系的問題

需要說明的是，在上述分析中使用的是\((x,y)\)坐標系，但是在OpenCV等常用圖像庫中往往使用圖像左上角為原點的\((x_{col},y_{row})\)的像素坐標系，會導致OpenCV中的Homograph矩陣與上述推導有一些差異。