python 图像到图像的映射

单应性变换

  单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表示。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，我们将频繁的使用单应性变换。本质上，单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：

  对于图像平面内（甚至是三维中的点，后面我们会介绍到）的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。我们通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标是的我们可以简单地使用一个矩阵来表示变换。

下面的函数可以实现对点进行归一化和转换齐次坐标的功能：

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数组。这种格式使得矩阵乘法和点的变换操作更加容易。对于其他的例子，比如对于聚类和分类的特征，我们将使用典型的行数组来存储数据。

  在这些投影变换中，有一些特别重要的变换。比如，仿射变换：

保持了w=1，不具有投影变换所具有的强大变形能力，仿射变换包括一个可逆矩阵A和一个平移向量t=[t_x,t_y]。仿射变换可以用于很多应用，比如图像扭曲。

线性变换

线性变换需满足的要求：

• 变换前是直线，变换后还是直线
• 直线保持比例不变
• 变换前是原点的，变换后依然是原点

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]

仿射变换

仿射变换需满足的条件：

• 变换前是直线，变换后还是直线
• 直线的比例保持不变

  从上面的条件可以发现，仿射变换和线性变换相比少了一个原点不变的条件，所以仿射变换相当于线性变换加平移。
  计算仿射矩阵H的方法，其中tp是变换后的坐标，fp是变换前的坐标，通过计算H，使得tp是fp通过仿射变换矩阵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]

图像扭曲

  仿射变换可应用于图像扭曲，当仿射变换应用于图像扭曲时，变换矩阵定义下图公式：

代码：

from numpy import *
from matplotlib.pyplot import *
from numpy.core._multiarray_umath import array
from scipy import ndimage
from PIL import Image

im = array(Image.open('img/dog.jpg').convert('L'))
H = array([[1.4, 0.3, -100], [0.3, 2.0, -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.3, -100], [0.3, 2.0, -100], [0, 0, 1]]，则s cosθ=1.4，-s sinθ=0.3，tx=-100，s sinθ=0.3，scosθ=2.0，ty=-100。

图像中的图像

  定义方法image_in_image：将一副图像放到另一幅图像中。

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

调用image_in_image方法，将猫放入狗的图片中

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

# 猫和狗
im1 = array(Image.open('img/cat.jpg').convert('L'))
im2 = array(Image.open('img/dog.jpg').convert('L'))

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

# 选定一些目标点
tp = array([[264, 538, 540, 264], [50, 46, 505, 505], [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，接着使用图像扭曲方法将两幅图像组合在一起。

三角形仿射

以下例子尝试将猫图像放到公路边的广告栏中。

  假如我们直接将猫图像放到公路边的广告栏中，可能会出现下面这种情况，底部和右侧填充不完整，细节处理不好。

  使用三角形仿射，图像将会以左上-右下分成两个三角形，再分别将这两个三角形放到指定坐标围成的区域中。

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

代码：

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

im1 = array(Image.open('img/cat.jpg').convert('L'))
im2 = array(Image.open('img/03.jpeg').convert('L'))

# 选定im1角上的一些点
m, n = im1.shape[:2]

fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])
tp = array([[228, 815, 787, 269],[102, 115, 465, 470] , [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()

运行结果：

可以看到使用三角形仿射，图像原本底部和右侧的细节问题已经解决了。

分段仿射扭曲

使用更多的三角形，就能实现更精确的匹配。

• 三角剖分的函数triangulate_points，获得三角形

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

• 通过得到的三角形对图像进行扭曲拼接，该函数为pw_affine

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

• 在拼接后的图片中绘制三角形，该函数为plot_mesh

# 绘制三角形
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')

• 手动抓取目标点代码部分

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

代码：

from array import array
from matplotlib import pyplot as plt
from numpy import *
from matplotlib.pyplot import *
from PIL import Image
from numpy import meshgrid, loadtxt, vstack, ones

import warp

# 打开图像，并将其扭曲
fromim = array(Image.open('img/cat.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/03.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 = np.array(tp)
# 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()

运行结果：

  从运行结果中我们可以看到，使用多个三角形，拼接的效果比两个三角形高，两个三角形的拼接效果可能会由于坐标选取的不够好而出现一点错位，比如上个例子中的左上角区域。

  目标点的抓取需要使用ginput函数手动抓取，抓取30个点，按照从上到下，从左到右的顺序依次抓取，抓取的点的位置尽量排列整齐，每行每列的点数需相同，这里例子中是采用了六行五列，运行结果效果比较好，如果随意抓取，可能就会出现图片的各个三角形之间错乱的情况。