單應性變換

1、單應性變換時將一個平面內的點映射到另一個平面內的二維投影變換，下圖中的兩幅圖片，存在對應的點，如第一幅橋頂的A點對應第二幅圖橋頂的A’點，但是這兩幅圖橋頂的點在影像中的左邊都不一樣，要找到這兩個坐標的對應關系，就用到了單應性變換，單應性變換就是一個平面到另一個平面的變換關系

2、假設A(x, y)、A’(x’, y’)，那么x和x’、y和y’就存在一種變換關系，定義變換關系矩陣為H，則A點和A’點的變換關系公式定義為下圖公式，從下圖的關系可以看出，x、y、x’、y’是可以比較容易得到的，所以關鍵問題是如何去求解H矩陣，由于不同的影像變換方式，H矩陣的取值會有所不同，下面主要講講線性變換和仿射變換

線性變換

1、線性變換需要滿足一下三點要求：

• 變換前是直線，變換后還是直線
• 直線保持比例不變
• 變換前是原點的，變換后依然是原點

例如，一個矩形，經過線性變換后（旋轉90度），原點（紅色點）還是在原點，位置沒有偏移，直線還是直線，且直線的比例沒有發生變化

2、將上面的單應性變換矩陣展開后得到x和x’、y和y’的關系如下

3、再將其轉成矩陣后可以表示為下圖的形式，我們可以使用SVD演算法找到H矩陣的最小二乘解

4、直接線性變換演算法（SVD），將下面的方法放入homography.py檔案中

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[1][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、仿射變換需要滿足一下兩點條件：

• 變換前是直線，變換后還是直線
• 直線的比例保持不變

從上面的條件可以發現，仿射變換和線性變換相比少了一個原點不變的條件，所以仿射變換相當于線性變換加平移，例如：（圖片來源網路）

2、仿射變換的變換矩陣定義為下圖公式，H矩陣最后一排h7=h8=0，h9=1，且權重w=1，c和 f 表示x和y的平移量

3、計算仿射矩陣H的方法，將下面的方法放到homography.py檔案中，其中tp是變換后的坐標，fp是變換前的坐標，通過計算H，使得tp是fp通過仿射變換矩陣H得到的，然后回傳H

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]

影像扭曲

1、仿射變換可以用于很多的應用，影像扭曲就是其中一個應用，當仿射變換應用于影像扭曲時，變換矩陣定義下圖公式，其中向量s指定了變換的尺度，如果s=1，那么該變換能夠保持距離不變，此時變換為剛體變換，θ表示旋轉（扭曲）的角度

2、影像扭曲代碼

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

im = array(Image.open('img/jmu01.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()

運行結果

根據代碼，設定了H矩陣值為[[1.4, 0.05, -100], [0.05, 1.5, -100], [0, 0, 1]]，則s cosθ=1.4，-s sinθ=0.05，tx=-100，s sinθ=0.05，scosθ=1.5，ty=-100，從上圖可以看到扭曲區域邊界之外的地方用0像素來填充，來創建一個二值的alpha影像，

影像中的影像

1、將一張影像放到另一張影像中，這是仿射扭曲一個比較簡單的應用，通過影像扭曲的特點，它能夠使一張影像的4個頂點和另一幅影像的指定位置對齊，

2、下面定義方法image_in_image：將一副影像放到另一幅影像中，將該方法放到warp.py檔案中，

def image_in_image(im1, im2, tp):
    """使用仿射變換將im1放置在im2上，使im1影像的角和tp盡可能的靠近
        tp是齊次表示的，并且是按照從左上角逆時針計算的"""
    # 扭曲的點
    m, n = im1.shape[:2]
    fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])

    # 計算仿射變換，并且將其應用于影像im1中
    # 計算仿射變換矩陣H
    H = 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

3、呼叫image_in_image方法，將陳嘉庚雕像放到中山紀念館前面

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

im1 = array(Image.open('img/jmu01.jpg').convert('L'))
im2 = array(Image.open('img/jmu12.jpg').convert('L'))

gray()
subplot(131)
imshow(im1)
axis('equal')
axis('off')
subplot(132)
imshow(im2)
axis('equal')
axis('off')

# 選定4個點的y，x坐標
tp = array([[264, 800, 800, 264], [40, 40, 400, 400], [1, 1, 1, 1]])

im3 = warp.image_in_image(im1, im2, tp)
subplot(133)
imshow(im3)
axis('equal')
axis('off')
show()

在這個案例中，在中山紀念館影像中設定了4個頂點的位置，存放到變數tp，然后傳入imge_in_image方法計算陳嘉庚影像的仿射變換矩陣H，接著使用影像扭曲方法將兩幅影像組合在一起

使用三角形仿射彎曲效果

1、為什么需要三角形仿射？如果不使用三角形仿射，我們把中山紀念館貼到一張公告牌上，那么會出現貼不滿的情況，如下圖，

2、這里簡單說兩個三角形的情況，如果使用三角形仿射，將中山紀念館影像從左上角到右下角分割成兩個三角形，然后先將左下角的三角形先貼到公告牌上，再將右上角的三角形貼到公告牌上，這樣就能夠使得公告牌被貼滿，

3、三角形仿射代碼

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

im1 = array(Image.open('img/jmu12.jpg').convert('L'))
im2 = array(Image.open('img/announce.jpg').convert('L'))

# 選定im1角上的一些點
print(im1.shape)
m, n = im1.shape[:2]

# 中山紀念館的4個頂點坐標
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])
# 映射倒公告牌中的坐標
tp = array([[50, 188, 175, 33], [60, 50, 375, 370], [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]]   # 第1、3、4個坐標
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[0, 2]), im2.shape[:2])

# 三角形的alpha影像
alpha = warp.alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im4 = (1 - alpha) * im3 + alpha * im1_t

figure()
gray()
imshow(im4)
axis("equal")
# axis("off")
show()

這里簡單的為每個三角形創建了alpha影像，然后將所有影像（三角形）合并起來，三角形alpha影像能夠簡單的通過檢查像素的坐標是否能夠寫成三角形頂點坐標的凸組合來計算得出，通過alpha_for_triangle函式來完成上述作業，代碼如下，

def alpha_for_triangle(points, m, n):
    """對于帶有由points定義角點的三角形，創建大小為（m,n）的alpha圖
    （在歸一化的齊次坐標意義下）"""

    alpha = zeros((m, n))
    for i in range(min(points[0]), max(points[0])):
        for j in range(min(points[1]), max(points[1])):
            x = linalg.solve(points, [i, j, 1])
            if min(x) > 0:
                alpha[i, j] = 1
    return alpha

4、三角形仿射運行結果

從影像上可以看到，相同的坐標，使用三角仿射后，使得影像能夠貼滿公告牌，我們可以看到影像從左上角到右下角被分割成了兩塊三角形，由于坐標設定的不好，使得兩個三角形拼接后發生了位置偏差的問題，

分段仿射扭曲

從上面的結果可以知道，將一張影像分成若干個三角形，可以實作角點的精確匹配，分的三角形越多，匹配越精確，

1、首先需要三角剖分的函式triangulate_points，獲得三角形，將其加入warp.py檔案中

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

2、然后是通過得到的三角形對影像進行扭曲拼接，該函式為pw_affine，將其放入warp.py檔案中

def pw_affine(fromim, toim, fp, tp, tri):
    """ 從一幅影像中扭曲矩形影像塊
    fromim= 將要扭曲的影像
    toim= 目標影像
    fp= 齊次坐標表示下，扭曲前的點
    tp= 齊次坐標表示下，扭曲后的點
    tri= 三角剖分 """

    im = toim.copy()

    # 檢查影像是灰度影像還是彩色圖象
    is_color = len(fromim.shape) == 3

    # 創建扭曲后的影像（如果需要對彩色影像的每個顏色通道進行迭代操作，那么有必要這樣做）
    im_t = zeros(im.shape, 'u8')

    for t in tri:
        # 計算仿射變換
        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])
                # im1_t = ndimage.affine_transform(im1, H[:2, :2],
                #                                  (H[0, 2], H[1, 2]), im2.shape[:2])
        else:
            im_t = ndimage.affine_transform(
                fromim, H[:2, :2], (H[0, 2], H[1, 2]), im.shape[:2])
        # 三角形的 alpha
        alpha = alpha_for_triangle(tp[:, t], im.shape[0], im.shape[1])

        # 將三角形加入到影像中
        im[alpha > 0] = im_t[alpha > 0]
    return im

3、然后在拼接后的圖片中繪制三角形的函式plot_mesh，將其放入warp.py檔案中

# 繪制三角形
def plot_mesh(x, y, tri):
    """ 繪制三角形 """
    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]]  # 將第一個點加入到最后
        plot(x[t_ext], y[t_ext], 'r')

4、最后是主程式呼叫上面的函式

from numpy import *
from matplotlib.pyplot import *
from PIL import Image
import warp

# 打開影像，并將其扭曲
fromim = array(Image.open('img/jmu12.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('img/announce.jpg').convert('L'))

tp = loadtxt('turningtorso1_points.txt')  # destination points

# 將點轉換成齊次坐標
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()
gray()
imshow(im)
warp.plot_mesh(tp[1], tp[0], tri)
axis('off')
show()

運行結果

通過上面的結果可以看到，多個三角形的拼接精確度，比兩個三角形的精確度要高，拼接后發生錯位的地方看起來沒有兩個三角形那么別扭，
完成實驗存在的一個問題是，如果傳入的是一幅彩色影像，pw_affine函式處理彩色影像的時候發生了too many indices for array錯誤，隨后查看了影像扭曲affine_transform函式的呼叫方式，發現影像扭曲的處理的是灰度影像，處理彩色影像和灰度影像的方式不一樣，于是將影像改成了灰度便可以運行成功，但是彩色影像的錯誤目前還無法弄清楚是什么原因
需要說明的是，運行結果的右圖中的所有目標點需要通過ginput函式去手動選取，然后將這些點寫到turningtorso1_points.txt檔案中，至少選擇30個點（因為triangulate_points函式得到的三角形的點最大下標是29），而且手動選取這些點的時候，需要從左向右、從上到下的順序進行采點，順序不能錯誤，否則拼接出來的三角形就會發生很大的錯位，同時選取的目標點不能隨意的亂點，需要按照上面的運行結果那樣有規律、整齊的選取，如果隨意亂取目標點，也會造成圖片錯位混亂，以下就是因為隨意選取目標點而導致將中山紀念館放到廣告牌上的時候變得很混亂

以下是選取目標點的代碼，選取目標點的時候ginput函式應該是會等待選取完30個點才會繼續運行代碼，但是實際上在我取完30個點之前，代碼就自動運行然后結束了，這個目前還不知道是什么原因，只能快速選完點，然后寫入txt檔案，最后通過txt檔案進行坐標調整

im = array(Image.open('img/announce.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)
print(tp)