Python計算機視覺——圖像到圖像的映射

文章目錄

• Python計算機視覺——圖像到圖像的映射
• • 寫在前面
  • 1 單應性變換
  • • 1.1 直接線性變換算法
    • 1.2 仿射變換
  • 2 圖像扭曲
  • • 2.1 圖像中的圖像
    • 2.2 分段仿射扭曲
    • 2.3 圖像配準
  • 3 創建全景圖
  • • 3.1 RANSAC
    • 3.2 拼接圖像

寫在前面

• 剛體變換：平移+旋轉，只改變物體位置，不改變物體形
• 仿射變換：改變物體位置和形狀，但是保持“平直性”
• 投影變換：徹底改變物體位置和形狀

1 單應性變換

單應性變換是將一個平面內的點映射到另一個平面內的二維投影變換。在這里，平 面是指圖像或者三維中的平面表面。單應性變換具有很強的實用性，比如圖像配準、 圖像糾正和紋理扭曲，以及創建全景圖像。本質上，單應性變換H，按照下面的方程映射二維中的點（齊次坐標意義下）：
$$?\begin{array}{l}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}?=?\begin{array}{lll}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}??\begin{array}{l}x\\ y\\ w\end{array}?\text{?或?}{\mathbf{x}}^{\prime }=\mathbf{H}\mathbf{x}$$
也就是原圖中的點為x，經過經過矩陣H的變換，就能變換為$x^{\prime }$等。點的齊次坐標是依賴于其尺度定義的，所以， x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1] 都表示同一個二維點。因此，單應性矩陣 H也僅 依賴尺度定義，所以，單應性矩陣具有 8 個獨立的自由度。我們通常使用 w=1 來歸 一化點，這樣，點具有唯一的圖像坐標 x 和 y。這個額外的坐標使得我們可以簡單地使用一個矩陣來表示變換。

1.1 直接線性變換算法

單應性矩陣可以由兩幅圖像中對應點對計算出來。且一個完全射影變換具有 8 個自由度。根據對應點約束，每個對應點對可以寫出兩個 方程，分別對應于 x 和 y 坐標。因此，計算單應性矩陣H需要４個對應點對。所謂的DLT（Direct Linear Transformation，直接線性變換）是給定４個或者更多對應點對 矩陣，來計算單應性矩陣H的算法。參數求解過程如下：
$$?\begin{array}{c}w{x}^{\prime }\\ w{y}^{\prime }\\ w\end{array}?=?\begin{array}{lll}{h}_{00}& h_{01}& h_{02}\\ h_{10}& h_{11}& h_{12}\\ h_{20}& h_{21}& h_{22}\end{array}??\begin{array}{l}x\\ y\\ 1\end{array}?$$
其中$w$=1，且
$$\begin{array}{r}{x}_i^{\prime }=\frac{h_{00}{x}_i+h_{01}{y}_i+h_{02}}{h_{20}{x}_i+h_{21}{y}_i+h_{22}}\\ y_i^{\prime }=\frac{h_{10}{x}_i+h_{11}{y}_i+h_{12}}{h_{20}{x}_i+h_{21}{y}_i+h_{22}}\end{array}$$
得到：
$$\begin{array}{rl}& x_i^{\prime }\left(h_{20}{x}_i+h_{21}{y}_i+h_{22}\right)=h_{00}{x}_i+h_{01}{y}_i+h_{02}\\ & y_i^{\prime }\left(h_{20}{x}_i+h_{21}{y}_i+h_{22}\right)=h_{10}{x}_i+h_{11}{y}_i+h_{12}\end{array}$$

$$\left[\begin{array}{ccccccccc}{x}_i& y_i& 1& 0& 0& 0& ?x_i^{\prime }{x}_i& ?x_i^{\prime }{y}_i& ?x_i^{\prime }\\ 0& 0& 0& x_i& y_i& 1& ?y_i^{\prime }{x}_i& ?y_i^{\prime }{y}_i& ?y_i^{\prime }\end{array}\right]?\begin{array}{l}{h}_{00}\\ h_{01}\\ h_{02}\\ h_{10}\\ h_{11}\\ h_{12}\\ h_{20}\\ h_{21}\\ h_{22}\end{array}?=\left[\begin{array}{l}0\\ 0\end{array}\right]$$

即：
$$\begin{array}{c}?\begin{array}{ccccccccc}{x}_1& y_1& 1& 0& 0& 0& ?x_1^{\prime }{x}_1& ?x_1^{\prime }{y}_1& ?x_1^{\prime }\\ 0& 0& 0& x_1& y_1& 1& ?y_1^{\prime }{x}_1& ?y_1^{\prime }{y}_1& ?y_1^{\prime }\\ \dots \dots & & & & & & & & \\ x_n& y_n& 1& 0& 0& 0& ?x_n^{\prime }{x}_n& ?x_n^{\prime }{y}_n& ?x_n^{\prime }\\ 0& 0& 0& x_n& y_n& 1& ?y_n^{\prime }{x}_n& ?y_n^{\prime }{y}_n& ?y_n^{\prime }\end{array}??\begin{array}{l}{h}_{00}\\ h_{01}\\ h_{02}\\ h_{10}\\ h_{11}\\ h_{12}\\ h_{20}\\ h_{21}\\ h_{22}\end{array}?=?\begin{array}{c}0\\ 0\\ ?\\ 0\\ 0\end{array}?\\ ?\end{array}$$

可記為Ah = 0 ，其中A是一個具有對應點對二倍數量行數的矩陣，且A的維度為2n×9，h的維度為9，0的維度為2n。根據最小二乘解，$\hat{h}=A^{T}A$最小特征值對應的特征向量。代碼如下：

def H_from_points(fp,tp):""" Find homography H, such that fp is mapped to tpusing the linear DLT method. Points are conditionedautomatically. """if fp.shape != tp.shape:raise RuntimeError('number of points do not match')# condition points (important for numerical reasons)# --from points--m = mean(fp[:2], axis=1)maxstd = max(std(fp[:2], axis=1)) + 1e-9C1 = diag([1/maxstd, 1/maxstd, 1]) C1[0][2] = -m[0]/maxstdC1[1][2] = -m[1]/maxstdfp = dot(C1,fp)# --to points--m = mean(tp[:2], axis=1)maxstd = max(std(tp[:2], axis=1)) + 1e-9C2 = diag([1/maxstd, 1/maxstd, 1])C2[0][2] = -m[0]/maxstdC2[1][2] = -m[1]/maxstdtp = dot(C2,tp)# create matrix for linear method, 2 rows for each correspondence pairnbr_correspondences = fp.shape[1]A = zeros((2*nbr_correspondences,9))for i in range(nbr_correspondences):        A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]U,S,V = linalg.svd(A)H = V[8].reshape((3,3))    # deconditionH = dot(linalg.inv(C2),dot(H,C1))# normalize and returnreturn H / H[2,2]

其實可以偽造一些數據來調用此函數，直觀感受下輸入輸出是什么。首先偽造輸入，即兩個隨機矩陣，大小為3×3，然后調用上述函數：

fp = np.random.randint(0, 500, size=(3, 3))
tp = np.random.randint(0, 500, size=(3, 3))
result = H_from_points(fp,tp)

然后result結果如下：

array([[ 2.96402671e+00, -1.27728682e+00, -1.22660693e+01],[ 1.44447830e+00, -3.97930961e-01, -3.30726226e+01],[ 4.88783355e-03, -3.40812578e-03,  1.00000000e+00]])

可以看到數組的最后一個值就是1。

1.2 仿射變換

仿射變換有 6 個自由度，因此需要三個對應點對來估計矩陣H。通過將 最后兩個元素設置為 0，即 h7=h8=0，仿射變換可以用上面的 DLT 算法估計得出。代碼如下：

def Haffine_from_points(fp,tp):""" Find H, affine transformation, such that tp is affine transf of fp. """if fp.shape != tp.shape:raise RuntimeError('number of points do not match')# condition points# --from points--m = mean(fp[:2], axis=1)maxstd = max(std(fp[:2], axis=1)) + 1e-9C1 = diag([1/maxstd, 1/maxstd, 1]) C1[0][2] = -m[0]/maxstdC1[1][2] = -m[1]/maxstdfp_cond = dot(C1,fp)# --to points--m = mean(tp[:2], axis=1)C2 = C1.copy() #must use same scaling for both point setsC2[0][2] = -m[0]/maxstdC2[1][2] = -m[1]/maxstdtp_cond = dot(C2,tp)# conditioned points have mean zero, so translation is zeroA = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)U,S,V = linalg.svd(A.T)# create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.tmp = V[:2].TB = tmp[:2]C = tmp[2:4]tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1) H = vstack((tmp2,[0,0,1]))# deconditionH = dot(linalg.inv(C2),dot(H,C1))return H / H[2,2]

同樣按照上述可看下函數的輸入輸出：

array([[ 860.17452003, -599.25709744, -111.19311744],[ 658.73707606, -458.91440696,  -85.53262921],[   0.        ,    0.        ,    1.        ]])

能看到最后一行的數據為0，0，1，因此也就對應了6個自由度。

2 圖像扭曲

所謂圖像扭曲，就是仿射扭曲，進行的操作即是改變物體位置和形狀，但是保持平直性。可先簡單看下效果：

2.1 圖像中的圖像

圖像扭曲的一個簡單例子是，將圖像或者圖像的一部分放置在另一幅圖像中，使得它們能夠和指定的區域或者標記物對齊。也就是把圖像放在圖像上，圖像中的圖像。函數如下：

def image_in_image(im1,im2,tp):""" Put im1 in im2 with an affine transformationsuch that corners are as close to tp as possible.tp are homogeneous and counter-clockwise from top left. """ # points to warp fromm,n = im1.shape[:2]fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])# compute affine transform and applyH = homography.Haffine_from_points(tp,fp)im1_t = ndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])alpha = (im1_t > 0)return (1-alpha)*im2 + alpha*im1_t

該函數的輸入參數為兩幅圖像和 一個坐標。該坐標為將第一幅圖像放置到第二幅圖像中的角點坐標：

函數 Haffine_from_points() 會返回給定對應點對的最優仿射變換。

2.2 分段仿射扭曲

分段仿射扭曲是對應點對集合之間最常用的扭曲方式。給定任意圖像的標記 點，通過將這些點進行三角剖分，然后使用仿射扭曲來扭曲每個三角形，我們可以 將圖像和另一幅圖像的對應標記點扭曲對應。對于任何圖形和圖像處理庫來說，這 些都是最基本的操作。下面我們來演示一下如何使用 Matplotlib 和 SciPy 來完成該操作。使用的是 Matplotlib中的狄洛克三角剖分：

import numpy as np 
from scipy.spatial import Delaunay
import matplotlib.pyplot as pltx,y = array(random.standard_normal((2,100))) 
tri = Delaunay(np.c_[x,y]).simplicesplt.figure(dpi = 140)
plt.subplot(121)
plt.plot(x,y,'*') 
plt.axis('off')
for t in tri:t_ext = [t[0], t[1], t[2], t[0]] # 將第一個點加入到最后 plot(x[t_ext],y[t_ext],'r')plt.subplot(122)plt.plot(x[t_ext],y[t_ext],'r')
plt.plot(x,y,'*')
plt.axis('off')
plt.show()

上圖顯示了一些實例點和三角剖分的結果。狄洛克三角剖分選擇一些三角形， 使三角剖分中所有三角形的最小角度最大。

2.3 圖像配準

圖像配準是對圖像進行變換，使變換后的圖像能夠在常見的坐標系中對齊。配準可以是嚴格配準，也可以是非嚴格配準。為了能夠進行圖像對比和更精細的圖像分析，圖像配準是一步非常重要的操作。圖像配準具有廣泛的應用，適用于同一個場景中有多張圖像需要進行匹配或疊加。在醫學圖像領域以及衛星圖像分析和光流（optical flow）方面非常普遍。

3 創建全景圖

在同一位置（即圖像的照相機位置相同）拍攝的兩幅或者多幅圖像是單應性相關的，我們可以使用該約束將很多圖像縫補起來，拼成一個大的圖像來創建全景圖。

3.1 RANSAC

RANSAC(RAndom SAmple Consensus,隨機采樣一致)算法是從一組含有“外點”(outliers)的數據中正確估計數學模型參數的迭代算法。“外點”一般指的的數據中的噪聲，比如說匹配中的誤匹配和估計曲線中的離群點。所以，RANSAC也是一種“外點”檢測算法。對于RANSAC算法來說一個基本的假設就是數據是由“內點”和“外點”組成的。“內點”就是組成模型參數的數據，“外點”就是不適合模型的數據。同時RANSAC假設：在給定一組含有少部分“內點”的數據，存在一個程序可以估計出符合“內點”的模型。可用代碼實現：

import numpy as np
import matplotlib.pyplot as plt
import random
import math# 數據量。
SIZE = 50
# 產生數據。np.linspace 返回一個一維數組，SIZE指定數組長度。
# 數組最小值是0，最大值是10。所有元素間隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10fig = plt.figure()
# 畫圖區域分成1行1列。選擇第一塊區域。
ax1 = fig.add_subplot(1,1, 1)
# 標題
ax1.set_title("RANSAC")# 讓散點圖的數據更加隨機并且添加一些噪聲。
random_x = []
random_y = []
# 添加直線隨機噪聲
for i in range(SIZE):random_x.append(X[i] + random.uniform(-0.5, 0.5)) random_y.append(Y[i] + random.uniform(-0.5, 0.5)) 
# 添加隨機噪聲
for i in range(SIZE):random_x.append(random.uniform(0,10))random_y.append(random.uniform(10,40))
RANDOM_X = np.array(random_x) # 散點圖的橫軸。
RANDOM_Y = np.array(random_y) # 散點圖的縱軸。# 畫散點圖。
ax1.scatter(RANDOM_X, RANDOM_Y)
# 橫軸名稱。
ax1.set_xlabel("x")
# 縱軸名稱。
ax1.set_ylabel("y")# 使用RANSAC算法估算模型
# 迭代最大次數，每次得到更好的估計會優化iters的數值
iters = 100000
# 數據和模型之間可接受的差值
sigma = 0.25
# 最好模型的參數估計和內點數目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正確模型的概率
P = 0.99
for i in range(iters):# 隨機在數據中紅選出兩個點去求解模型sample_index = random.sample(range(SIZE * 2),2)x_1 = RANDOM_X[sample_index[0]]x_2 = RANDOM_X[sample_index[1]]y_1 = RANDOM_Y[sample_index[0]]y_2 = RANDOM_Y[sample_index[1]]# y = ax + b 求解出a，ba = (y_2 - y_1) / (x_2 - x_1)b = y_1 - a * x_1# 算出內點數目total_inlier = 0for index in range(SIZE * 2):y_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1# 判斷當前的模型是否比之前估算的模型好if total_inlier > pretotal:iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))pretotal = total_inlierbest_a = abest_b = b# 判斷是否當前模型已經符合超過一半的點if total_inlier > SIZE:break# 用我們得到的最佳估計畫圖
Y = best_a * RANDOM_X + best_b# 直線圖
ax1.plot(RANDOM_X, Y)plt.show()

普通最小二乘是保守派：在現有數據下，如何實現最優。是從一個整體誤差最小的角度去考慮，盡量誰也不得罪。

RANSAC是改革派：首先假設數據具有某種特性（目的），為了達到目的，適當割舍一些現有的數據。

3.2 拼接圖像

實際情況中，創建全景圖的步驟大致分為三步：圖像獲取、圖像配準、圖像融合。兩張待配準的圖像處于不同的坐標系中，配準的目的就是將他們投影到拼接平面（即同一坐標系下）對齊。拼接的前提是，我們的圖片集中圖片之間要有共視區域，這樣我們才能在各圖像之間提取特征并尋找匹配點，然后在盡量去除錯配點對的情況下，根據匹配點計算圖像間的映射關系，最后將圖像根據映射關系拼接以獲得更寬廣的可視面積，多次進行配準拼接融合步驟，就得到了我們所需的全景圖。

總結

以上是生活随笔為你收集整理的Python计算机视觉——图像到图像的映射的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。