計算機視覺-影像到影像的映射(單應性變換、仿射變換、仿射扭曲)

1.原理

單應性變換（Homography）

齊次坐標系

齊次坐標系( x , y , w x,y,w x,y,w)與常見的三維空間坐標系( x , y , z x,y,z x,y,z)不同，只有兩個自由度,其中 w w w（ w w w>0）對應坐標 x x x和 y y y的縮放尺度:

當 w w w=1與 w w w=0時：

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

：齊次坐標系在計算機影像將3維物體投影到2維平面中起到很大的作用，

單應性變換

單應性變換是將一個平面內的點映射到另一個平面內的二維投影變換，（平面是指影像或者三維中的平面表面），對應的變換矩陣稱為單應性矩陣，

用矩陣的形式可以理解為：

其中( x l x_l xl?, y l y_l yl?)是Left view圖片上的點， ( x r x_r xr?, y r y_r yr?)是Right view圖片上對應的點,H為單應性矩陣 ，使w=1來歸一化點，

每一組匹配點( x i x_i xi?, y i y_i yi?)與( x i ′ x_i' xi′?, y i ′ y_i' yi′?)可以得出：

上式可以表示為：

變換上式可以得到：

寫成矩陣 A H = 0 AH=0 AH=0形式：

根據對應點約束，每個對應點對可以寫出兩個方程，分別對應于 x x x和 y y y坐標，

由下式可以看出單應性矩陣H與aH其實完全一樣（其中a≠0），

即點( x i x_i xi?, y i y_i yi?)無論經過 H H H還是 a H aH aH映射，變化后都是 ( x i ′ x_i' xi′?, y i ′ y_i' yi′?)，
如果使令a=1/ h 33 h_{33} h33?，那么有:

可以看出單應性矩陣 H H H雖然有9個未知數，但只有8個自由度，

由齊次坐標系也可以看出，點的齊次坐標是依賴于其尺度定義的，所以， x = [ x , y , w ] = [ α x , α y , α w ] = [ x / w , y / w , 1 ] x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1] x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1]都表示同一個二維點，因此，單應性矩陣 H H H也僅依賴尺度定義，所以，單應性矩陣具有 8 個獨立的自由度，

根據對應點約束，每個對應點對可以寫出兩個方程，分別對應于 x x x和 y y y坐標，因此，計算單應性矩陣 H H H需要４個對應點對，

SVD(奇異值分解）

DLT（Direct Linear Transformation，直接線性變換）是給定４個或者更多對應點對矩陣，來計算單應性矩陣 H H H的演算法，將單應性矩陣 H H H作用在對應點對上，重新寫出該方程，我們可以得到下面的方程：

我們顯然很難用直接線性的方法來求得H的解，可以用SVD演算法找到h的最小二乘解，具體參考奇異值分解

演算法代碼

def H_from_points(fp,tp):
    """使用線性DLT方法，計算單應性矩陣H，使fp映射到tp，點自動進行歸一化"""
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
    # 對點進行歸一化（對數值計算很重要）
    # ---映射起始點---
    m = mean(fp[:2],axis=1)
    maxstd = max(std(fp[:2],axis=1)) + 1e-9
    C1 = diag([1/maxstd,1/maxstd,1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = m[1]/maxstd
    fp = dot(C1,fp)
    #--- 映射對應點----
    m = mean(tp[:2],axis=1)
    maxstd = max(std(tp[:2],axis=1)) +1e-9
    
    C2 = diag([1/maxstd,1/maxstd,1])
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    
    tp = dot(C2,tp)
    
    #創建用于線性方法的矩陣，對于每個對應對，在矩陣中會出現兩行數值
    nbr_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[0][i]]
    U,S,V = linalg.svd(A)
    H = V[8].reshape((3,3))
    
    #反歸一化
    H = dot(linalg.inv(C2),dot(H,C1))
    
    #歸一化，然后回傳
    return H / H[2,2]

1)檢查fp tp兩張影像矩陣中行列的數目是否相同，如果不相同，函式將會拋出例外資訊，
2)對這些點進行歸一化操作，使其均值為 0，方差為 1，因為演算法的穩定性取決于坐標的表示情況和部分數值計算的問題，所以歸一化操作非常重要，
3)使用對應點來構造矩陣A，最小二乘解即為矩陣SVD分解后得矩陣V的最后一行，改行經過變形后得到矩陣H，
4)對這個矩陣進行反歸一化處理，
5)對矩陣H進行反歸一化，然后回傳，

仿射變換（affine）

通俗的來說，仿射變換就是線性變換+平移，變換前是直線的，變換后依然是直線且直線比例保持不變，

不通俗的說，仿射變換是一種二維坐標到二維坐標之間的線性變換（相同平面），它保持了二維圖形的“平直性”（直線經過變換之后依然是直線）和“平行性”（二維圖形之間的相對位置關系保持不變，平行線依然是平行線，且直線上點的位置順序不變），但是角度會改變，任意的仿射變換都能表示為乘以一個矩陣(線性變換)，再加上一個向量 (平移) 的形式，

由于仿射變換具有 6 個自由度，因此我們需要三個對應點對來估計矩陣 H H H，通過將最后兩個元素設定為 0，即 h7=h8=0，仿射變換可以用上面的 DLT 演算法估計得出，

演算法代碼

def Haffine_from_points(fp, tp):
    """計算H仿射變換，使得tp是fp經過仿射變換H得到的"""

    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    # 對點進行歸一化（對數值計算很重要）
    # --- 映射起始點 ---
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1)) + 1e-9
    C1 = diag([1/maxstd, 1/maxstd, 1]) 
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp_cond = dot(C1,fp)
    
    # --- 映射對應點 ---
    m = mean(tp[:2], axis=1)
    C2 = C1.copy() # 兩個點集，必須都進行相同的縮放
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp_cond = dot(C2,tp)
    
    # 因為歸一化后點的均值為0，所以平移量為0
    A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = linalg.svd(A.T)
    
    # 如Hartley和Zisserman著的Multiplr View Geometry In Computer,Scond Edition所示，
    # 創建矩陣B和C
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]
    
    tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1) 
    H = vstack((tmp2,[0,0,1]))
    
    # 反歸一化
    H = dot(linalg.inv(C2),dot(H,C1))
    
    return H / H[2,2]

2.影像扭曲（仿射扭曲）

對影像塊應用仿射變換，我們將其稱為影像扭曲（或者仿射扭曲）
扭曲操作可以使用SciPy 工具包中的 ndimage 包來簡單完成，命令為：

transformed_im = ndimage.affine_transform(im,A,b,size)

使用如上所示的一個線性變換 A 和一個平移向量 b 來對影像塊應用仿射變換，選項引數 size 可以用來指定輸出影像的大小，默認輸出影像設定為和原始影像同樣大小，

測驗代碼

# -*- coding=utf-8 -*-
# name: nan chen
# date: 2021/4/08 11:11

from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image

im = array(Image.open(r'D:\project3image\shangda001.jpg').convert('L'))
H = array([[1.4, 0.05, -100], [0.05, 1.5, -100], [0, 0, 1]])
im2 = ndimage.affine_transform(im, H[:2, :2], (H[0, 2], H[1, 2]))

gray()
subplot(121)
imshow(im)
axis('off')
subplot(122)
imshow(im2)
axis('off')
show()

運行結果

2.1影像中的影像

目標：通過仿射扭曲變換將影像放置到另一幅影像中，使得它們能夠和指定的區域或者標記物對齊

代碼

手動定位坐標點來調整影像中的位置過于繁瑣

定位獲取點坐標

import cv2
img = cv2.imread(r'D:\project3image\ggp.jpg')


def on_EVENT_LBUTTONDOWN(event, x, y, flags, param):
    if event == cv2.EVENT_LBUTTONDOWN:
        xy = "%d,%d" % (x, y)
        cv2.circle(img, (x, y), 1, (255, 0, 0), thickness=-1)
        cv2.putText(img, xy, (x, y), cv2.FONT_HERSHEY_PLAIN,
                    1.0, (0, 255, 0), thickness=1)
        cv2.imshow("image", img)


cv2.namedWindow("image")
cv2.setMouseCallback("image", on_EVENT_LBUTTONDOWN)
cv2.imshow("image", img)

while (True):
    try:
        cv2.waitKey(100)
    except Exception:
        cv2.destroyAllWindows()
        break

cv2.waitKey(0)
cv2.destroyAllWindows()

結果：

完全影像的仿射扭曲

# -*- coding=utf-8 -*-
# name: nan chen
# date: 2021/4/6 10:36
from PCV.geometry import warp, homography
from PIL import Image
from pylab import *

# 兩張圖片
im1 = array(Image.open(r'D:\project3image\jmu001.jpg').convert('L'))
im2 = array(Image.open(r'D:\project3image\ggp.jpg').convert('L'))

# 設定映射的目標點
tp = array([[45, 308, 319, 53], [254, 247, 591, 599], [1, 1, 1, 1]])

# 使用仿射變換將im1放置在im2上，使im1影像的角和tp盡可能的靠近
im3 = warp.image_in_image(im1, im2, tp)

# 將影像灰度顯示
figure()
gray()
subplot(131)
axis('off')
imshow(im1)
subplot(132)
axis('off')
imshow(im2)
subplot(133)
axis('off')
imshow(im3)
show()

使用兩個三角形的仿射彎曲

# -*- coding=utf-8 -*-
# name: nan chen
# date: 2021/4/6 10:36
from PCV.geometry import warp, homography
from PIL import Image
from pylab import *
from scipy import ndimage

# 兩張圖片
im1 = array(Image.open(r'D:\project3image\jmu001.jpg').convert('L'))
im2 = array(Image.open(r'D:\project3image\ggp.jpg').convert('L'))
gray()
# 設定映射的目標點
tp = array([[45, 308, 319, 53], [254, 247, 591, 599], [1, 1, 1, 1]])

# 選定 im1 角上的一些點
m, n = im1.shape[:2]
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])
# 第一個三角形
tp2 = tp[:, :3]
fp2 = fp[:, :3]
# 計算 H
H = homography.Haffine_from_points(tp2, fp2)
im1_t = ndimage.affine_transform(im1, H[:2, :2], (H[0, 2], H[1, 2]), im2.shape[:2])
# 三角形的 alpha
alpha = warp.alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im3 = (1 - alpha) * im2 + alpha * im1_t
# 第二個三角形
tp2 = tp[:, [0, 2, 3]]
fp2 = fp[:, [0, 2, 3]]
# 計算 H
H = homography.Haffine_from_points(tp2, fp2)
im1_t = ndimage.affine_transform(im1, H[:2, :2],(H[0, 2], H[1, 2]), im2.shape[:2])
# 三角形的 alpha 影像
alpha = warp.alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im4 = (1 - alpha) * im3 + alpha * im1_t
imshow(im4)
axis('equal')
axis('off')
show()

運行結果

完全影像的仿射扭曲

細節

使用兩個三角形的仿射彎曲

細節

結果分析

第一張圖中的p1為集美大學尚大樓的灰度影像，p2為帶有廣告牌的灰度影像，p3是通過完全影像的仿射扭曲將集美大學尚大樓映射至與p2的廣告牌位置對齊，p4為使用兩個三角形的仿射彎曲，通過修改tp = array([[45, 308, 319, 53], [254, 247, 591, 599], [1, 1, 1, 1]])這一行的坐標點可以改變p1映射在p2中的位置，其中前四個數字代表四個角點的縱坐標，中間四個數字代表四個角點的橫坐標，四個角點的順序為從左上角開始按照逆時針方向排序，最后四個數字代表四個角點的α通道，四個1就表示四個角點的透明度均為不透明，將上述用代碼獲取的坐標修改之后便可以得到上面的效果，通過下面的細節圖可以明顯看出完全影像的仿射扭曲變換后映射到廣告牌上的尚大樓影像邊緣不太光滑，使用包含兩個三角形的仿射變換的效果更好，

2.2分段仿射扭曲

從上面的例子可以看出，三角形影像塊的仿射扭曲可以完成角點的精確匹配，三角形影像塊越多，則匹配的越精確，分段仿射扭曲是通過給定任意影像的標記點，通過將這些點進行三角剖分，然后使用仿射扭曲來扭曲每個三角形，然后將影像和另一幅影像的對應標記點扭曲對應，

三角剖分函式

# 三角剖分的函式
def triangulate_points(x, y):
    """二維點的 Delaunay 三角剖分"""
    tri = Delaunay(np.c_[x, y]).simplices
    return tri

分段仿射影像扭曲的通用扭曲函式

def pw_affine(fromim,toim,fp,tp,tri):
    """ Warp triangular patches from an image.
        fromim = image to warp 
        toim = destination image
        fp = from points in hom. coordinates
        tp = to points in hom.  coordinates
        tri = triangulation. """
                
    im = toim.copy()
    
    # check if image is grayscale or color
    is_color = len(fromim.shape) == 3
    
    # create image to warp to (needed if iterate colors)
    im_t = zeros(im.shape, 'uint8') 
    
    for t in tri:
        # compute affine transformation
        H = homography.Haffine_from_points(tp[:,t],fp[:,t])
        
        if is_color:
            for col in range(fromim.shape[2]):
                im_t[:,:,col] = ndimage.affine_transform(
                    fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        else:
            im_t = ndimage.affine_transform(
                    fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        
        # alpha for triangle
        alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
        
        # add triangle to image
        im[alpha>0] = im_t[alpha>0]
        
    return im

繪制三角形函式

def plot_mesh(x,y,tri):
    """ Plot triangles. """ 
    
    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]] # add first point to end
        plot(x[t_ext],y[t_ext],'r')

測驗代碼

上述的三個函式都可以在PCV庫中warp.py代碼中找到，可以匯入包直接呼叫，下面的代碼是對上面的函式進行呼叫來測驗分段仿射扭曲匹配的結果，
代碼：

from PCV.geometry import warp
from PIL import Image
from pylab import *

# 打開影像，并將其扭曲
fromim = array(Image.open(r'D:\project3image\jmu001.jpg').convert('L'))
x, y = meshgrid(range(5), range(6))
x = (fromim.shape[1] / 4) * x.flatten()
y = (fromim.shape[0] / 5) * y.flatten()
# 三角剖分
tri = warp.triangulate_points(x, y)
# 打開影像
im = array(Image.open(r'D:\project3image\ggp.jpg').convert('L'))
gray()
imshow(im)
# 手工選取目標點
tp = plt.ginput(30)
for i in range(0, len(tp)):
    tp[i] = list(tp[i])
    tp[i][0] = int(tp[i][0])
    tp[i][1] = int(tp[i][1])
tp = array(tp)

# 將點轉換成齊次坐標
fp = vstack((y, x, ones((1, len(x)))))
tp = vstack((tp[:, 1], tp[:, 0], ones((1, len(tp)))))
# 扭曲三角形
im = warp.pw_affine(fromim, im, fp, tp, tri)
# 繪制影像
figure()
imshow(im)
# 繪制三角形
warp.plot_mesh(tp[1], tp[0], tri)
axis('off')
show()

結果分析

手工選取目標點：（必須按照從左到右 從上往下的順序整齊的選取目標的，方可得到想要的結果，選取的目標點越準確，匹配的結果越好，如果是隨意選取的話將得到錯亂的結果）

由上圖可以看出多個三角形的選取匹配后的結果要比兩個三角形匹配的結果要好，理論上來說只要選取的點足夠精確，可以達到很好的效果，

3.實驗中遇到的問題

（1）運行代碼時遇到ModuleNotFoundError: No module named ‘matplotlib.delaunay’
解決方法：
將import matplotlib.delaunay as md
改為from scipy.spatial import Delaunay
后將wrap.py中的triangulate_points函式中的陳述句替換為

tri = Delaunay(np.c_[x,y]).simplices

（2）直接運行教材中分段仿射扭曲代碼，由于沒有txt檔案，若直接用陣列取值，會出現下圖的錯誤，

解決方法：使用ginput() 函式手工選取，

參考：

知乎-單應性變換