序言

??本章講解了影像之間的變換，以及一些計算變換的應用，這些應用可以用于影像扭曲變換和影像配準，

運行環境：

1. python3.8
2. PyCharm2020.3

一、單應性變換

??單應性變換 是將一個平面內的點映射到另一個平面內的二維投影變換，在這里，平面是指影像或者三維中的平面表示，單應性變換具有很強的實用性，比如影像配準，影像糾正和紋理扭曲，以及創建全景影像，我們將頻繁的使用單應性變換，本質上，單應性變換H，按照下面的方程映射二維中的點（齊次坐標意義下）： [ x ′ y ′ w ′ ] = [ h 1 h 2 h 3 h 4 h 5 h 6 h 7 h 8 h 9 ] [ x y z ] \begin{bmatrix} x^{'} \\ y^{'} \\ w^{'}\end{bmatrix}=\begin{bmatrix}h_1 & h_2 & h_3 \\ h_4 & h_5 & h_6 \\ h_7 & h_8 & h_9\end{bmatrix}\begin{bmatrix}x \\ y \\ z\end{bmatrix} ???x′y′w′????=???h1?h4?h7??h2?h5?h8??h3?h6?h9????????xyz????或者 x ′ = H x x^{'}=Hx x′=Hx對于影像平面內（甚至是三維中的點，后面我們會介紹到）的點，齊次坐標是個非常有用的表示方式，點的齊次坐標是依賴于其尺度定義的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一個二維點，因此，單應性矩陣H也僅依賴尺度定義，所以，單應性矩陣具有8個獨立的自由度，我們通常使用w=1來歸一化點，這樣，點具有唯一的影像坐標x和y，這個額外的坐標是的我們可以簡單地使用一個矩陣來表示變換，

??我們創建 homography.py 檔案，用下面的函式可以實作對點進行歸一化和轉換齊次坐標的功能：

from numpy import *

def normallize(points):
    """在齊次坐標意義下，對點集進行歸一化，是最后一行為1"""
    for row in points:
        row /= points[-1]
    return points

def make_homog(points):
    """將點集（dim×n的陣列）轉換為齊次坐標表示"""

    return vstack((points,ones((1, points.shape[1]))))

??進行點和變換的處理時，我們會按照列優先的原則存盤這些點，因此，n個二維點集將會存盤為齊次坐標意義下的一個3×n陣列，這種格式使得矩陣乘法和點的變換操作更加容易，對于其他的例子，比如對于聚類和分類的特征，我們將使用典型的行陣列來存盤資料，

??在這些投影變換中，有一些特別重要的變換，比如，仿射變換： [ x ′ y ′ 1 ] = [ a 1 a 2 t x a 3 a 4 t y 0 0 1 ] [ x y 1 ] \begin{bmatrix}x^{'} \\ y^{'} \\ 1\end{bmatrix}=\begin{bmatrix}a_1 & a_2 & t_x \\ a_3 & a_4 & t_y \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}x \\ y \\ 1\end{bmatrix} ???x′y′1????=???a1?a3?0?a2?a4?0?tx?ty?1???????xy1????或 x ′ = [ A t 0 1 ] x x^{'}=\begin{bmatrix}A & t \\ 0 & 1\end{bmatrix}x x′=[A0?t1?]x保持了w=1，不具有投影變換所具有的強大變形能力，反射變換包括一個可逆矩陣A和一個平移向量t=[t_x,t_y]，仿射變換可以用于很多應用，比如 影像扭曲，

相似變換： [ x ′ y ′ 1 ] = [ s cos ? ( θ ) ? s sin ? ( θ ) t x s sin ? ( θ ) s cos ? ( θ ) t y 0 0 1 ] [ x y 1 ] \begin{bmatrix}x^{'} \\ y^{'} \\ 1\end{bmatrix}=\begin{bmatrix}s\cos(\theta) & -s\sin(\theta) & t_x \\ s\sin(\theta) & s\cos(\theta) & t_y \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}x \\ y \\ 1\end{bmatrix} ???x′y′1????=???scos(θ)ssin(θ)0??ssin(θ)scos(θ)0?tx?ty?1???????xy1????或 x ′ = [ s R t 0 1 ] x x^{'}=\begin{bmatrix}sR & t \\ 0 & 1\end{bmatrix}x x′=[sR0?t1?]x是一個包含尺度變化的二維剛體變換，上式中的向量s指定了變換的尺度，R是角度為θ的旋轉矩陣，t=[t_x,t_y]在這里也是一個平移向量，如果s=1，那么該變換能夠保持距離不變，此時，變換稱為剛體變換，相似變換可以用于很多應用，比如影像配準，

??下面讓我們來一起探討如何設計用于估計單應性矩陣的演算法，然后看一下使用仿射變換進行影像扭曲，使用相似變換進行影像匹配，以及使用完全投影變換進行創建全景影像的一些例子：

1.1 直接線性變換演算法

??單應性矩陣可以有兩幅影像（或者平面）中對應點對計算出來，前面已經提到過，一個完全射影變換具有8個自由度，根據對應點約束，每個對應點對可以寫出兩個方程，分別對應于x和y坐標，因此，計算單應性矩陣H需要4個對應點對，

??DLT（Direct Linear Transformation，直接線性變換）是給定4個點或者更多對應點對矩陣，來計算單應性矩陣H的演算法，將單應性矩陣H作用在對應點上，重新寫出該方程，我們可以得到下面的方程： [ ? x 1 ? y 1 ? 1 0 0 0 x 1 x 1 ′ y 1 x 1 ′ x 1 ′ 0 0 0 ? x 1 ? y 1 ? 1 x 1 y 1 ′ y 1 y 1 ′ y 1 ′ ? x 2 ? y 2 ? 1 0 0 0 x 2 x 2 ′ y 2 x 2 ′ x 2 ′ 0 0 0 ? x 2 ? y 2 ? 1 x 2 y 2 ′ y 2 y 2 ′ y 2 ′ . . . . . . . . . . . . ] [ h 1 h 2 h 3 h 4 h 5 h 6 h 7 h 8 h 9 ] = 0 \begin{bmatrix} -x_1& -y_1&-1& 0&0&0&x_1x_1^{'}&y_1x_1^{'}&x_1^{'} \\ 0&0&0&-x_1&-y_1&-1&x_1y_1{'}&y_1y_1{'}&y_1^{'} \\ -x_2&-y_2&-1&0&0&0&x_2x_2^{'}&y_2x_2^{'}&x_2^{'} \\ 0&0&0&-x_2&-y_2&-1&x_2y_2^{'}&y_2y_2^{'}&y_2^{'} \\ &...&&...&&...&&... \end{bmatrix}\begin{bmatrix}h_1 \\ h_2\\ h_3\\ h_4 \\ h_5\\h_6\\ h_7\\ h_8\\ h_9\end{bmatrix}=0 ????????x1?0?x2?0??y1?0?y2?0...??10?10?0?x1?0?x2?...?0?y1?0?y2??0?10?1...?x1?x1′?x1?y1?′x2?x2′?x2?y2′??y1?x1′?y1?y1?′y2?x2′?y2?y2′?...?x1′?y1′?x2′?y2′????????????????????????h1?h2?h3?h4?h5?h6?h7?h8?h9?????????????????=0或者Ah=0，其中A是一個具有對應點對二倍數量行數的矩陣，將這些對應點對方程的系數堆疊到一個矩陣紅，我們可以使用SVD演算法找到H的最小二乘解，下面是演算法代碼，我們也將其加入到 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]

??上面的函式的第一步操作是檢查點對兩個陣列中點的數目是否相同，如果不相同，函式將會拋出例外資訊，這對于寫出穩健的代碼來說非常有用，

??代碼先對這些點進行歸一化操作，使其均值為0，方差為1，因為演算法的穩定性取決于坐標的表示情況和部分數值計算的問題，所以歸一化操作非常重要，接下來我們使用對應點對來構造矩陣A，最小二乘解即為矩陣SVD分解后所得矩陣V的最后一行，該行經過變換后得到矩陣H，然后對這個矩陣進行處理和歸一化，回傳輸出，

1.2 仿射變換

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

??下面的函式使用對應點來計算放射變換矩陣，我們繼續將其添加到 homography.py 檔案中，

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]

??同樣的，類似于DLT演算法，這些點需要經過預處理和去處理化操作，

二、影像扭曲

??對影像塊應用仿射變換，我們將其稱為影像扭曲（或者 仿射扭曲），該操作不僅經常在計算機圖形學中，而且經常出現在計算機視覺演算法總，扭曲的操作可以使用SciPy工具包中的ndimage包來簡單完成，命令：

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

??使用上面所示的一個線性變換A和一個平移向量b來對影像塊應用放射變換，選項引數 size 可以用來指定輸出影像的大小，

??我們可以運行下列代碼來研究該函式是如何作業的：

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

# 解決中文亂碼
rcParams['font.sans-serif'] = 'SimHei'
rcParams['axes.unicode_minus'] = False

im = array(Image.open('image/image1.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)
title("原始影像")
axis('off')
subplot(122)
imshow(im2)
title("扭曲后的影像")
axis('off')
show()

運行結果：

??從運行結果可以看出，原始影像（左）和扭曲后的影像（右）的差別，扭曲后的影像中丟失的像素用零來填充，

2.1 影像中的影像

??仿射扭曲的一個簡單例子是，將影像或者影像的一部分放置在另一幅影像中，是的他們能夠和指定的區域或者標記物對齊，

??將函式 image_in_image() 添加到 wary.py 檔案中，該函式的輸入引數為兩幅影像和一個坐標，該坐標為將第一幅影像放置到第二幅影像中的角點坐標：

import homography
from scipy import ndimage
from numpy import *

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 = 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

??將扭曲的影像和第二幅影像融合，就創建 alpha 影像，該影像定義了每個像素從各個影像中獲取的像素值成分多少，這里基于以下事實：扭曲的影像是在扭曲區域邊界之外以 0 來填充的影像，來創建一個二值的 alpha 影像，嚴格意義上，需要在第一幅圖象中的潛在 0 像素上加上一個小的數值，或者合理的處理這些 0 像素，

注意：我們使用的影像坐標是齊次坐標意義下的，

使用仿射變換放置影像

??我們試著使用該函式將一幅影像插入另一幅影像，

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

# 將im1仿射扭曲到im2的指定位置
im1 = array(Image.open('image/1.jpg').convert('L'))
im2 = array(Image.open('image/image3.jpg').convert('L'))

figure()
gray()
subplot(121)
imshow(im1)
axis('equal')
axis('off')
subplot(122)
imshow(im2)
axis('equal')
axis('off')

# 選定一些目標點
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

# 呼叫的warp.py的image_in_image函式，從而實作仿射變換
im3 = warp.image_in_image(im1, im2, tp)

figure()
imshow(im3)
axis('equal')
axis('off')
show()

運行結果：

??從運行結果上看，圖片將原本右下角的一個公告去給替換了，因為 tp 中的齊次坐標是對應公告牌的坐標，所以將其放置到該位置，那么，

使用ginput函式獲取 tp 中的齊次坐標值

??其中，tp 中選取的點的坐標值可以通過手工確定，也可以通過 PyLab 類別庫中的 ginput() 函式獲得，可以使用下面的代碼實作點齊次坐標的獲取：

from PIL import Image
from pylab import *
import numpy as np

im = array(Image.open('image/image3.jpg'))
imshow(im)

print('Please click 4 point')
points = ginput(4)

# 轉整數
int_points = np.int_(points)

# 行列互換操作
new_points = list(map(list, zip(*int_points)))

print('you clicked x:', new_points[1])
print('you clicked y:', new_points[0])
show()

運行結果：

選取點的順序	運行結果

??然后就可以用 x 的坐標替換掉 tp 中的第一個引數串列，用 y 替換掉 tp 中的第而個引數串列：

# tp = array([[264, 538, 540, 264], [40, 36, 605, 605], [1, 1, 1, 1]])
# 替換成
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

注意：運行上述代碼選取點時，要按照上圖的給的順序點擊，不然可能會出現圖片位置錯誤、圖片旋轉等一些問題！因為 image_in_image() 函式中的 tp 是按照從左上角逆時針計算的，

??函式 Haffine_from_points() 會回傳給定對應點對的最優仿射變換，在上面的例子中，對應點對為影像和公告牌的角點，如果透視效應比較弱，那么這種方法會回傳很好的結果，

??對于三個點，仿射變換可以將一幅影像進行扭曲，使這三對對應點對可以完美地匹配上，這是因為，仿射變換具有 6 個自由度，三個對應點對可以給出 6 個約束條件（對于這三個對應點對，x 和 y 坐標必須都要匹配），所以，如果你真的打算使用仿射變換將影像放置到公告牌上，可以將影像分成兩個三角形，然后對它們分別進行扭曲影像操作，代碼如下：

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

im1 = array(Image.open('image/1.jpg').convert('L'))
im2 = array(Image.open('image/image3.jpg').convert('L'))

# 選取 im1 角上的一些點
m, n = im1.shape[:2]
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])

# 選定一些目標點
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

# 第一個三角形
tp2 = tp[:, :3]
fp2 = fp[:, :3]

# 計算 H
H = 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 = 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 = 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 = alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im4 = (1 - alpha) * im3 + alpha * im1_t

figure()
gray()
subplot(121)
imshow(im3)
axis('equal')
axis('off')
subplot(122)
imshow(im4)
axis('equal')
axis('off')
show()

運行結果：

??這里我們簡單地為每個三角形創建了 alpha 影像，然后將所有的影像合并起來，該三角形的 alpha 影像可以簡單地通過檢查像素的坐標是否能夠寫成三角形頂點坐標的凸組合來計算得出，如果坐標可以表示成這種形式，那么該像素就位于三角形的內部，上面的例子使用了下面的函式 alpha_for_triangle()，將其添加到 warp.py 檔案中，

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

??這樣顯卡可以極其快速地操作上面的代碼，Python 語言的處理速度比你的顯卡（或者 C/C++ 實作）慢很多，但是對于我們來說已經夠用了，

2.2 分段仿射扭曲

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

為了三角化這些點，我們經常使用狄洛克三角剖分方法，在 Matplotlib（但是不在 PyLab 庫中）中有狄洛克三角剖分，我們可以用下面的方式使用它：

from numpy import *
from matplotlib.pyplot import *
from scipy.spatial import Delaunay

x, y = array(random.standard_normal((2, 100)))
tri = Delaunay(np.c_[x, y]).simplices

figure()
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] # 將第一個點加入到最后
    plot(x[t_ext], y[t_ext], 'r')

plot(x, y, '*')
axis('off')

figure()
plot(x, y, '*')
axis('off')
show()

運行結果：

??運行結果顯示了一些實體點和三角剖分的結果，狄洛克三角剖分選擇一些三角形，使三角剖分中所有三角形的最小角度最大（三角剖分中的邊實際上是泰森圖的對偶圖），函式 delaunay() 有 4 個輸出，其中我們僅需要三角形串列資訊（第三個輸出），在 warp.py 檔案中創建用于三角剖分的函式：

from scipy.spatial import Delaunay

def triangulate_points(x, y):
    """二維點的 Delaunay 三角剖分"""

    # centers, edges, tri, neighbors = Delaunay(x, y)
    tri = Delaunay(np.c_[x, y]).simplices
    return tri

??函式輸出的是一個陣列，該陣列的每一行包含對應陣列 x 和 y 中每個三角形三個點的切片，

??現在將該演算法應用于一個例子，在該例子中，在 5×6 的網格上使用 30 個控制點，將一幅影像扭曲到另一幅影像中的非平坦區域，圖 3-5b 所示的是將一幅影像扭曲到“turning torso”的表面，目標點是使用 ginput() 函式手工選取出來的，將結果保存在 turningtorso_points.txt 檔案中；

??然后我們使用下面簡短的腳本將這些操作統一起來：

 # -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.geometry import warp

# 打開影像，并將其扭曲
fromim = array(Image.open('image/sunset_tree.jpg'))
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('image/turningtorso1.jpg'))
tp = loadtxt('turningtorso1_points.txt') # destination points

figure()
subplot(1, 4, 1)
axis('off')
imshow(im)

# 將點轉換成齊次坐標
fp = array(vstack((y, x, ones((1, len(x))))), 'int')
tp = array(vstack((tp[:, 1], tp[:, 0], ones((1, len(tp))))), 'int')

# 扭曲三角形
im = warp.pw_affine(fromim, im, fp, tp, tri)

# 繪制影像
subplot(1, 4, 2)
axis('off')
imshow(fromim)
warp.plot_mesh(fp[1], fp[0], tri)
subplot(1, 4, 3)
axis('off')
imshow(im)
subplot(1, 4, 4)
axis('off')
imshow(im)
warp.plot_mesh(tp[1], tp[0], tri)
show()

??再通過下面的輔助函式（將其添加到 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')

運行結果：

??使用狄洛克三角剖分標記點進行分段仿射扭曲：第一張為帶有標記物的目標影像；第二張為帶有三角剖分的影像；第三張為扭曲后的影像；第四張為帶有三角剖分的扭曲影像，

如果運行教材上的代碼報錯：ModuleNotFoundError: No module named ‘matplotlib.delaunay’
將 import matplotlib.delaunay as md 改成 from scipy.spatial import Delaunay
并將 triangulate_points() 函式中的 centers, edges, tri, neighbors = Delaunay(x, y) 改成 tri = Delaunay(np.c_[x,y]).simplices

使用狄洛克三角剖分標記點進行分段仿射扭曲使用的 turningtorso1_points.txt 檔案可以從該網站上下載: Programming Computer Vision with Python