图像到图像的映射

一、映射变换

1.单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。本质上，单应性变换 H，按照下面的方程映射二维中的点（齐次坐标意义下）：

对于图像平面内（甚至是三维中的点，后面我们会介绍到）的点， 齐次坐
标是个非常有用的表示方式。

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]))))

↑ 建立齐次坐标

单应性矩阵主要解决以下两个问题：

• 表述真实世界中一个平面与对应它图像的透视变换
• 从通过透视变换实现图像从一种视图变换到另外一种视图

2.直接线性变换算法

单应性矩阵可以由两幅图像中对应点对计算出来，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H至少需要4个对应点对。

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

DLT实现:

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_from_points()

先对点进行归一化处理，使其均值为0，方差为1。因为坐标的表示情况和部分数值的计算问题决定了算法的稳定性，所以归一化很重要。
之后，就可以通过对应点对来构造矩阵A。最小二乘解即是矩阵SVD分解后的矩阵V的最后一行。该行经过变形后就可以得到矩阵H。然后对矩阵进行处理和归一化，返回输出。

3. 仿射变换

由于仿射变换具有 6 个自由度，因此我们需要三个对应点对来估计矩阵 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]

↑ Haffine_from_points()

二、图像扭曲

1.ndimage的仿射变换

from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
import pylab
import cv2 as cv

im = array(Image.open('thebeatles.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')
pylab.show()
cv.imwrite('img_trans.jpg',im2)

↑ img_trans.py

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

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

上面的函数使用一个线性变换 A 和一个平移向量 b 来对图像块应用仿射变换。选项参数 size 可以用来指定输出图像的大小，默认输出图像设置为和原始图像同样大小。

2.图像扭曲

将输出结果和原始图像对比可以看到，经过仿射变换后的扭曲图像丢失掉的像素用零来填充(也就是图像周围那圈黑色)

三、图像中的图像img_in_img

1.完全图像的仿射扭曲

from scipy import ndimage
import homography


def image_in_image(im1, im2, tp):
    """使用仿射变换将im1放置在im2上，使im1图像的角和tp尽可能的靠近
        tp是齐次表示的，并且是按照从左上角逆时针计算的"""

    # 扭曲的点
    m, n = im1.shape[:2]
    fp = homography.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

↑ warp.py

将扭曲的图像和第二幅图像融合，我们就创建了 alpha 图像。严格意义上说，我们需要在第一幅图像中的潜在 0 像素上加上一个小的数值，或者合理地处理这些 0 像素。

from numpy import *
from matplotlib.pyplot import *
from PIL import Image
import wrap
import cv2 as cv

im1 = array(Image.open('Abbey_road.jpg').convert('L'))
im2 = array(Image.open('AD_board.jpg').convert('L'))

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

# 选定一些目标点 [[y的坐标][x的坐标][1,1,1,1]]
tp = array([[44,338,333,55], [112,112,524,518], [1, 1, 1, 1]])

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

cv.imwrite('imginimg.jpg',im3)

↑ img_in_img.py

这里的tp坐标是事先用OpenCV选好的点坐标，这里是鼠标事件获取点坐标的代码。

结果：

但是，放大图3我们就会发现图1贴在广告绿幕的时候不能完全匹配广告牌范围的四个角：

↑左上角

↑右下角

所以，我们在进行仿射变换时，可以使用一个小技巧。

2.复数三角形仿射

对于三个点，仿射变换可以将一幅图像进行扭曲，使得三对对应点可以完美匹配上。这是因为，仿射变换具有6个自由度，三个对应点对可以给出6个约束条件。所以，当我们使用仿射变换将图像投射到广告牌上的时候，可以把图像分成两个三角形，然后对他们呢分别进行图像扭曲。

# 三角形图像仿射
# 选定 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
figure()
gray()
imshow(im4)
axis('equal')
axis('off')
show()

subplot(133)
gray()
imshow(im3)
axis('equal')
axis('off')
show()

↑这段代码替换img_in_img.py中的wrap.image_in_image()

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

↑这段代码添加到warp.py



↑两个三角形进行仿射映射
使用两个三角形进行仿射变换，还带有一些透视效应。如上图，包含两个三角形的仿射变换将图像映射到广告位，图像的贴合程度有肉眼可见的提高。

五、分段仿射扭曲

我们发现，三角形图像块的仿射扭曲可以完成角点的精确匹配。所以，给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。这就是分段仿射扭曲。

1.三角剖分

from scipy.spatial import Delaunay
import numpy as np


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

↑ triangulate_points()

这里会有一个十分严重的问题：
ModuleNotFoundError: No module named ‘matplotlib.delaunay’
这是因为书上的例子比较陈旧，所以才会有这样的问题。
这里是→完整解决方法

2. 分段仿射变换

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, 'uint8')
    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])
        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

↑ pw_affine()全称 piecewise affine 部分仿射

在这里，将仿射图像分成多个三角形后，对于每个小三角，仿射变换是唯一确定的，所以还要用到homography.py中的Haffine_from_points()来得到仿射变换矩阵。

3.多个三角形的分段仿射变换效果

1.运行代码：

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

# 打开将要扭曲的图像，并将其扭曲
fromim = array(Image.open(r'img\mao.png'))
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'img\yintong.jpg'))
imshow(im)

tp = loadtxt(r'data\coordinate.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()
imshow(im)
warp.plot_mesh(tp[1],tp[0],tri)
axis('off')
show()

↑ img_triangle_trans.py

在img_tirangle_rans.py中可以通过改变

x,y = meshgrid(range(5),range(6))

的range和shape在目标图像中选取更多的点，仿射图像会细分成更多块，能应对更复杂的情况。

PS：coordinate.txt 是提前在目标图像上选取目标点时，写入的像素坐标，
这是选点的过程：

点这里→看实现的代码

2.图像

↑ 将要仿射的图像

↑ 目标图像

↑ 大功告成 猫猫引桐楼(゜∀。)