2021-04-11

一丶求解单应性变换矩阵和仿射变换矩阵

1.原理

1.1 单应性变换矩阵

单应性变换就是一个平面到另一个平面的映射关系

如图，两张图片中相同颜色的点叫做corresponding Points，比如两个红点就是一对corresponding points。单应性矩阵(Homography)就是一个从一张图像到另一张图像映射关系的转换矩阵(3*3)。可以由下面的公式来表示：

以图像中的红点为例，可以将单应性变换写成如下形式：

1.2仿射变换矩阵

仿射变换变化包括缩放（Scale、平移(transform)、旋转(rotate)、反射（reflection,对图形照镜子）、错切(shear mapping，感觉像是一个图形的倒影)，原来的直线仿射变换后还是直线，原来的平行线经过仿射变换之后还是平行线，这就是仿射

仿射变换中集合中的一些性质保持不变：
（1）凸性
（2）共线性：若几个点变换前在一条线上，则仿射变换后仍然在一条线上
（3）平行性：若两条线变换前平行，则变换后仍然平行
（4）共线比例不变性：变换前一条线上两条线段的比例，在变换后比例不变

2 仿射变换数学表达
一个集合 XX 的仿射变换为：
f ( x ) = A x + b , x ∈ X f(x)=Ax+b, \quad x\in X
f(x)=Ax+b,x∈X

仿射变换是二维平面中一种重要的变换，在图像图形领域有广泛的应用，在二维图像变换中，一般表达为：

可以视为线性变换R和平移变换T的叠加

2.代码

def normalize(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]))))

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]

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[:]),axis=0)
    U,S,V = linalg.svd(A.T)

    #如Hartley和Zisserman著的。。。所示
    #创建矩阵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]

二丶图像扭曲

1.原理

仿射变换（Affine Transformation 或Affine Map）是一种二维坐标（x, y）到二维坐标（u, v）的线性变换，其数学表达式形式如下：

对应的齐次坐标矩阵表示形式为：

仿射变换保持了二维图形的“平直性”（直线经仿射变换后依然为直线）和“平行性”（直线之间的相对位置关系保持不变，平行线经仿射变换后依然为平行线，且直线上点的位置顺序不会发生变化）。非共线的三对对应点确定一个唯一的仿射变换。

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

2.代码

import matplotlib.pyplot as plt
from PIL import Image

import warp
from numpy import array

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

im1 = array(Image.open('C:\\Users\\Gypsophila\\Desktop\\1.jpg').convert('L'))
im2 = array(Image.open('C:\\Users\\Gypsophila\\Desktop\\2.jpg').convert('L'))

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

# 选定一些目标点
tp = array([[264, 538, 540, 264], [40, 36, 605, 605], [1, 1, 1, 1]])

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

3.结果

三丶图像中的图像

1.原理

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

2.代码

import matplotlib.pyplot as plt
from PIL import Image

import homography
from warp import *
from numpy import array

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
    
im1 = array(Image.open('1.jpg').convert('L'))
im2 = array(Image.open('2.jpg').convert('L'))

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

# 选定一些目标点
tp = array([[264, 538, 540, 264], [40, 36, 605, 605], [1, 1, 1, 1]])

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

3.结果