2021-04-11

1、单应性变换

单应性变换，可简单理解为用来描述物体在世界坐标系和像素坐标系之间的位置映射关系，对应的变换矩阵H称为单应性矩阵。对于由于不同的图像变换方式，H矩阵的取值会有所不同。假设两张图像中的对应点对对齐次坐标为（x’,y’,1）和（x,y,1），单应矩阵H定义：
于是有：

即 x=Hx’，将矩阵展开后有3个等式，将第3个等式代入前两个等式中可得：
可见一个点对应了两个等式的存在，其中每个等式未知数有8个，即8个自由度。

对于8个未知数的H的计算：
将上述等式乘以分母展开，得到如下：

整理得到，得到：

假如已知两幅图片中对应的N个点对（特征点匹配对），可以得到如下线性方程组：

写成矩阵形式：

由于单应矩阵H包含了||H||=1约束，因此根据上图的线性方程，8自由度的H我们至少需要4对对应的点才能计算出单应矩阵，然而，在真实的应用场景中，计算的点对中都会包含噪声，比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象，如果只使用4个点对来计算单应矩阵，那会出现很大的误差。因此，为了使得计算更精确，一般都会使用可以用直接线性变换算法算法计算H，将单应性矩阵H作用于对应点对上，重新写出一个齐次方程Ah=0，其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，可以使用奇异值分解找到H的最小二乘解。
H_from_points(fp,tp)函数算法代码：

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]

2、仿射变换

仿射变换就是图像的线性变换加上平移，用一幅图表示，就是 ：

由 image1 到 image2 的转换经过了三个操作，旋转 (线性变换)，缩放操作(线性变换)，平移 (向量加)。任意的仿射变换都能表示为乘以一个矩阵(线性变换)，再加上一个向量 (平移) 的形式。单应性变换有8个自由度，仿射有6个自由度：

其中令 h31=h32=0，有6个自由度，因此我们需要三个对应点来估计矩阵H。与单应性变换求H矩阵一样，仿射变换参数的求取也可采用上面的DLT算法估计H，在真实的应用场景中，计算的点对中都会包含噪声，比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象。

算法代码：Haffine_form_points(fp,tp)函数：

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]

3、仿射变换的应用

1、图像扭曲

1.代码：

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

im = array(Image.open('D:\\PhotoTest\\11.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')
show()

2.实现结果

2、图像映射融合（将一副图像融合进另一幅图像）

1.代码：

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

# 图像11.jpg到图像12.jpg的扭曲

im1 = array(Image.open(r'D:\PhotoTest\11.jpg').convert('L'))
im2 = array(Image.open(r'D:\PhotoTest\12.jpg').convert('L'))
# tp变量是图像变换目标坐标，它按照从左上角逆时针的顺序，定义第一幅图像的四个角的点在第二幅中的位置。
tp = array([[180,474,474,180],[340,340,735,735],[1,1,1,1]])

im3 = warp.image_in_image(im1,im2,tp)
figure()
gray()
subplot(141)
axis('off')
imshow(im1)
subplot(142)
axis('off')
imshow(im2)
subplot(143)
axis('off')
imshow(im3)

# 设置点tp在图片12.jpg的位置
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]
# compute 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 = 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]]
# compute 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 = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im3 + alpha*im1_t
subplot(144)
imshow(im4)
axis('off')
show()

2.实现结果

4、遇到的问题

导入书中代码时，运行图像融合代码时候发现的错误，出现：
几经查找，在一个博主上的博客找到了解决方法附上博文链接https://blog.csdn.net/weixin_42648848/article/details/88667243（解决方法）：
1.找到导入错误的包

2.将 import matplotlib.delaunay as md 改成 from scipy.spatial import Delaunay

3.将找到函数里把triangulate_points(x,y)里面的代码替换成 tri = Delaunay(np.c_[x,y]).simplices