Python计算机视觉 第3章-图像到图像的映射

Python计算机视觉 第3章-图像到图像的映射

3.1　单应性变换

单应性变换（Homography）是计算机视觉中非常重要的一种几何变换，它用于将一个平面内的点映射到另一个平面内。具体来说，单应性变换可以描述一个图像在摄像机视角变化、平面移动或旋转时，如何从一个视角变换到另一个视角。

这种变换在多个应用场景中非常有用，比如：

1. 图像配准：将不同视角或不同时间拍摄的图像对齐，找到它们之间的对应关系。
2. 图像校正：修正由于摄像机角度或透视导致的图像扭曲，使图像看起来更平整。
3. 纹理扭曲：将一个平面的纹理准确地映射到另一个平面上。
4. 全景图像创建：将多个图像拼接成一个大的全景图像。

单应性变换的频繁使用，尤其是在涉及多个视角或需要精确对齐图像的情况下，能够显著提升算法的鲁棒性和精度。在项目中，理解和正确应用单应性变换是处理图像和三维几何信息的关键技能。

单应性变换（Homography）将二维平面上的点映射到另一个平面上的点，在齐次坐标（homogeneous coordinates）下，这种映射可以通过以下方程来表示：

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}\right)\mathbf{H}\cdot \left(\begin{array}{c}x\\ y\\ w\end{array}\right)$$

其中，单应性矩阵 $\mathbf{H}$ 为：

$$\mathbf{H}=\left(\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right)$$

经过单应性变换后的目标点的常规二维坐标 $(x^{\prime },y^{\prime })$ 为：

$$x^{\prime }=\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}}$$

$$y^{\prime }=\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}}$$

通过这些公式，你可以描述平面间的各种变换，比如旋转、缩放、平移、透视变换等。

3.1.1　直接线性变换算法

单应性矩阵可以由两幅图像（或者平面）中对应点对计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要４个对应点对。

DLT（Direct Linear Transformation，直接线性变换）是给定４个或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点对上，重新写出该方程，我们可以得到下面的方程：

$$\left[\begin{array}{ccccccccc}-x_1& -y_1& -1& 0& 0& 0& x_1{x}_1^{\prime }& y_1{x}_1^{\prime }& x_1^{\prime }\\ 0& 0& 0& -x_1& -y_1& -1& x_1{y}_1^{\prime }& y_1{y}_1^{\prime }& y_1^{\prime }\\ -x_2& -y_2& -1& 0& 0& 0& x_2{x}_2^{\prime }& y_2{x}_2^{\prime }& x_2^{\prime }\\ 0& 0& 0& -x_2& -y_2& -1& x_2{y}_2^{\prime }& y_2{y}_2^{\prime }& y_2^{\prime }\\ & ⋮& & ⋮& ⋮& ⋮& ⋮& ⋮& \end{array}\right]\left[\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\\ h_9\end{array}\right]=0$$

或者Ah=0，其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用SVD（Singular Value Decomposition，奇异值分解）算法找到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[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]

上面函数的第一步操作是检查点对的两个数组中点的数目是否相同。如果不相同，函数将会抛出异常信息。这对于写出稳健的代码来说非常有用。但是，为了使得代码例子更简单、更容易理解，我们在本书中仅在很少的例子中使用异常处理技巧。

3.1.2　仿射变换

由于仿射变换具有6个自由度，因此我们需要三个对应点对来估计矩阵H。通过将最后两个元素设置为0，即$h_7=h_8=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 著的 Multiple View Geometry in Computer, Second 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算法，这些点需要经过预处理和去处理化操作。

3.2　图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常应用在计算机图形学中，而且经常出现在计算机视觉算法中。扭曲操作可以使用SciPy工具包中的ndimage包来简单完成。

以下为实验代码：

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from skimage import data, color

# 读取示例图像
image = color.rgb2gray(data.astronaut())

# 定义仿射变换矩阵
# 例如，这里是一个旋转矩阵和一个平移矩阵的组合
affine_matrix = np.array([
    [1.2, 0.2, -30],  # x轴的缩放和旋转，以及平移
    [0.1, 1.2, 20],   # y轴的缩放和旋转，以及平移
    [0, 0, 1]         # 齐次坐标的归一化因子
])

# 对图像应用仿射变换
transformed_image = ndimage.affine_transform(
    image,
    affine_matrix[:2, :2],  # 2x2 仿射矩阵
    offset=affine_matrix[:2, 2],  # 平移偏移
    mode='reflect'  # 边界处理模式
)

# 显示原始和变换后的图像
plt.figure(figsize=(10, 5))

plt.subplot(1, 2, 1)
plt.title('Original Image')
plt.imshow(image, cmap='gray')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.title('Transformed Image')
plt.imshow(transformed_image, cmap='gray')
plt.axis('off')

plt.show()

实验图1 图像扭曲处理结果

3.2.1　图像中的图像

仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from skimage import io, color


def image_in_image(background, overlay, position):
    """
    将图像 overlay 放置到图像 background 中的指定位置。

    :param background: 背景图像
    :param overlay: 要放置的图像
    :param position: 放置 overlay 的坐标 (x, y) 元组
    :return: 带有 overlay 的背景图像
    """
    # 确保 overlay 图像的尺寸
    h, w = overlay.shape[:2]

    # 生成仿射变换矩阵，将 overlay 图像的四个角点对齐到背景图像中的指定区域
    src_points = np.array([
        [0, 0],  # overlay 的左上角
        [w, 0],  # overlay 的右上角
        [w, h],  # overlay 的右下角
        [0, h]  # overlay 的左下角
    ])

    dst_points = np.array([
        [position[0], position[1]],  # 背景图像中放置位置的左上角
        [position[0] + w, position[1]],  # 背景图像中放置位置的右上角
        [position[0] + w, position[1] + h],  # 背景图像中放置位置的右下角
        [position[0], position[1] + h]  # 背景图像中放置位置的左下角
    ])

    # 构建矩阵 A 和向量 b 以求解仿射变换矩阵
    A = []
    b = []
    for i in range(4):
        A.append([src_points[i][0], src_points[i][1], 1, 0, 0, 0, -dst_points[i][0] * src_points[i][0],
                  -dst_points[i][0] * src_points[i][1]])
        A.append([0, 0, 0, src_points[i][0], src_points[i][1], 1, -dst_points[i][1] * src_points[i][0],
                  -dst_points[i][1] * src_points[i][1]])
        b.append(dst_points[i][0])
        b.append(dst_points[i][1])

    A = np.array(A)
    b = np.array(b)

    # 通过最小二乘法求解仿射变换矩阵
    h = np.linalg.lstsq(A, b, rcond=None)[0]
    H = np.append(h, [1]).reshape(3, 3)

    # 将 overlay 图像进行仿射变换
    transformed_overlay = ndimage.affine_transform(
        overlay,
        H[:2, :2],
        offset=H[:2, 2],
        output_shape=background.shape,
        mode='constant',
        cval=0
    )

    # 合并图像
    mask = (transformed_overlay > 0).astype(float)
    result = background.copy()
    result = result * (1 - mask) + transformed_overlay * mask

    return result


# 示例使用
if __name__ == "__main__":
    # 读取内置示例图像
    background = color.rgb2gray(io.imread('img.png'))  # 背景图像
    overlay = color.rgb2gray(io.imread('python.png'))  # 要放置的图像
    position = (100, 100)  # 放置位置（x, y）

    # 应用函数
    result_image = image_in_image(background, overlay, position)

    # 显示结果
    plt.figure(figsize=(10, 5))
    plt.subplot(1, 3, 1)
    plt.title('Background Image')
    plt.imshow(background, cmap='gray')
    plt.axis('off')

    plt.subplot(1, 3, 2)
    plt.title('Overlay Image')
    plt.imshow(overlay, cmap='gray')
    plt.axis('off')

    plt.subplot(1, 3, 3)
    plt.title('Result Image')
    plt.imshow(result_image, cmap='gray')
    plt.axis('off')

    plt.show()

实验图2 处理结果

3.2.2　分段仿射扭曲

对应点对集合之间最常用的扭曲方式：分段仿射扭曲。给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这些都是最基本的操作。
为了三角化这些点，我们经常使用狄洛克三角剖分方法。在Matplotlib（但是不在PyLab 库中）中有狄洛克三角剖分，我们可以用下面的方式使用它：

以下是实验代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage
from scipy.spatial import Delaunay
 
x,y = np.array(np.random.standard_normal((2,100)))
tri = Delaunay(np.c_[x, y]).simplices
plt.figure() 
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
    plt.plot(x[t_ext],y[t_ext],'r')
plt.plot(x,y,'*')
plt.axis('off')
plt.show()

实验图3 处理结果

3.2.3　图像配准

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准，也可以是非严格配准。为了能够进行图像对比和更精细的图像分析，图像配准是一步非常重要的操作。

3.3　创建全景图

在同一位置（即图像的照相机位置相同）拍摄的两幅或者多幅图像是单应性相关的（如图3-9所示）。我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建全景图像。

图3-9 瑞典隆德主要大学建筑的５幅图像。这些图像都是从同一个视点拍摄的

3.3.1　RANSAC

RANSAC是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

RANSAC的标准例子：用一条直线拟合带有噪声数据的点集。简单的最小二乘在该例子中可能会失效，但是RANSAC能够挑选出正确的点，然后获取能够正确拟合的直线。

图3-10 使用RANSAC算法用一条直线来拟合包含噪声数据点集

3.3.2　拼接图像

估计出图像间的单应性矩阵（使用RANSAC算法），现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面（否则，需要进行大量变形）。一种方法是创建一个很大的图像，比如图像中全部填充0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充０，以便为扭曲的图像腾出空间。以下为示例代码：

def panorama(H, fromim, toim, padding=2400, delta=2400):
    """ 使用单应性矩阵 H（使用 RANSAC 健壮性估计得出），协调两幅图像，创建水平全景图像。结果为一幅和 toim 具有相同高度的图像。padding 指定填充像素的数目，delta 指定额外的平移量。 """
    
    # 检查图像是灰度图像，还是彩色图像
    is_color = len(fromim.shape) == 3

    # 用于 geometric_transform() 的单应性变换
    def transf(p):
        p2 = np.dot(H, [p[0], p[1], 1])
        return (p2[0] / p2[2], p2[1] / p2[2])

    if H[1, 2] < 0:  # fromim 在右边
        print('warp - right')

        # 变换 fromim
        if is_color:
            # 在目标图像的右边填充 0
            toim_t = np.hstack((toim, np.zeros((toim.shape[0], padding, 3))))
            fromim_t = np.zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(
                    fromim[:, :, col], transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的右边填充 0
            toim_t = np.hstack((toim, np.zeros((toim.shape[0], padding))))
            fromim_t = ndimage.geometric_transform(fromim, transf, (toim.shape[0], toim.shape[1] + padding))
    
    else:  # fromim 在左边
        print('warp - left')

        # 为了补偿填充效果，在左边加入平移量
        H_delta = np.array([[1, 0, 0], [0, 1, -delta], [0, 0, 1]])
        H = np.dot(H, H_delta)

        # 变换 fromim
        if is_color:
            # 在目标图像的左边填充 0
            toim_t = np.hstack((np.zeros((toim.shape[0], padding, 3)), toim))
            fromim_t = np.zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(
                    fromim[:, :, col], transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的左边填充 0
            toim_t = np.hstack((np.zeros((toim.shape[0], padding)), toim))
            fromim_t = ndimage.geometric_transform(fromim, transf, (toim.shape[0], toim.shape[1] + padding))

    # 协调后返回（将 fromim 放置在 toim 上）
    if is_color:
        # 所有非黑色像素
        alpha = ((fromim_t[:, :, 0] > 0) | (fromim_t[:, :, 1] > 0) | (fromim_t[:, :, 2] > 0))
        for col in range(3):
            toim_t[:, :, col] = fromim_t[:, :, col] * alpha + toim_t[:, :, col] * ~alpha
    else:
        alpha = (fromim_t > 0)
        toim_t = fromim_t * alpha + toim_t * ~alpha

    return toim_t

对于通用的geometric_transform() 函数，我们需要指定能够描述像素到像素间映射的函数。在这个例子中，transf()函数就是该指定的函数。该函数通过将像素和H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，我们可以决定应该将该图像填补到左边还是右边。当该图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边”情况中，需要在单应性矩阵中加入平移。简单起见，我们同样使用0像素的技巧来寻找alpha图。现在在图像中使用该操作，函数如下所示：

# 扭曲图像
delta = 2000  # 用于填充和平移

# 读取图像
im1 = np.array(Image.open(imname[1]))
im2 = np.array(Image.open(imname[2]))

# 图像拼接
im_12 = warp.panorama(H_12, im1, im2, delta, delta)

im1 = np.array(Image.open(imname[0]))
im_02 = warp.panorama(np.dot(H_12, H_01), im1, im_12, delta, delta)

im1 = np.array(Image.open(imname[3]))
im_32 = warp.panorama(H_32, im1, im_02, delta, delta)

im1 = np.array(Image.open(imname[j + 1]))  # 确保 imname[j + 1] 是一个有效的索引
im_42 = warp.panorama(np.dot(H_32, H_43), im1, im_32, delta, 2 * delta)