计算机视觉学习blog（三）——图像到图像的映射

单应性变换

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

或者 x’=Hx

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

创建homography.py文件加入下面的函数，实现对点进行归一化和转换齐次坐标的功能：

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

单应性矩阵

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

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

或者Ah=0，其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵红，我们可以使用SVD（Singular Value Decomposition奇异值分解）算法找到H的最小二乘解。下面是算法代码，我们也将其加入到 homography.py 文件里：
H_from_points函数

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]

仿射变换

由于仿射变换具有6个自由度，因此我们需要三个对应点来估计矩阵H。通过将最后两个元素设置为0，即h7=h8=0，仿射变换可以用上面的DLT算法估计得出。

下面的函数使用对应点来计算放射变换矩阵，我们继续将其添加到 homography.py 文件中：

Haffine_from_points函数

# 仿射变换
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(8)), zeros((2, 1))), axis=1)
        H = vstack((tmp2, [0, 0, 1]))

        # 反归一化
        H = dot(linalg.inv(C2), dot(H, C1))

        return H / H[2, 2]

图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）

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

# 仿射扭曲
from scipy import ndimage
from numpy import *
from PIL import Image
from pylab import *
im1 = array(Image.open('3.jpg').convert('L'))
H = array([[1.4, 0.05, -100], [0.05, 1.5, -100], [0, 0, 1]])
im2 = ndimage.affine_transform(im1, H[:2, :2], (H[0, 2], H[1, 2]))

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

从运行结果可以看出，原始图像（左）和扭曲后的图像（右）的差别，扭曲后的图像中丢失的像素用零来填充。

图3-1：用仿射变换扭曲图像。原始图像（左图）以及使用ndimage.affine_transform()函数扭曲后的图像（右图）

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

将函数 image_in_image() 添加到 wary.py 文件中。该函数的输入参数为两幅图像和一个坐标。该坐标为将第一幅图像放置到第二幅图像中的角点坐标：

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

# -*- coding: utf-8 -*-

from pcv.geometry import warp, homography
from PIL import  Image
from pylab import *
from scipy import ndimage

# im1仿射到im2的例子
#打开两张图片
im1 = array(Image.open('3.jpg').convert('L'))
im2 = array(Image.open('9.jpg').convert('L'))
# 选定一些目标点，数组tp表示目标区域，三行分别代表y、x和齐次，其中区域按从左上角逆时针旋转,
tp = array([[120,260,260,120],[16,16,305,305],[1,1,1,1]])#将一些目标点存储到数组中
im3 = warp.image_in_image(im1,im2,tp) #使用仿射变换将m1放置到m2
figure()
gray()
subplot(141)
axis('off')
imshow(im1)
subplot(132)
axis('off')
imshow(im2)
subplot(133)
axis('off')
imshow(im3)
show()

运行结果：

图3-2：使用仿射变换将一副图像放置到另一幅图像中

修改参数tp = array([[233, 301, 303, 241],[23, 18, 333, 332],[1,1,1,1]])

坐标值可以通过查看绘制的图像（在Pylab图像中，鼠标的坐标显示在图像底部附近）手工确定，也可以使用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()

从左上角逆时针选取四个点

修改后结果：

图3-3：将图像放置在学校门口的石碑上

以上图片来源集美大学实拍建筑（陈延奎图书馆，集美大学东大门）
https://www.cnblogs.com/Easy-/p/14638459.html