python计算机视觉-图像处理基础章节第三章之根据仿射或单应性变换实现图像的扭曲，映射，融合

---

前言

了解一些基本的概念

图像映射的类型:

• 平移(translation)
• 旋转(rotation)
• 尺度变换(scale)
• 仿射(affine)
• 透视映射(Perspective)

• 刚体变换：平移+旋转，只改变物体位置，不改变物体形
• 仿射变换：改变物体位置和形状，但是保持“平直性”
• 投影变换：彻底改变物体位置和形状

---

一、单应性（Homography）

单应性变换是讲一个平面内的点映射到另一个平面内的二维投影变换

单应性变换的可以用到哪些地方？

• 图像配准
• 图像纠正
• 纹理扭曲
• 创建全景图像

单应性是什么？
用 无镜头畸变 的相机从不同位置拍摄 同一平面物体 的图像之间存在单应性，可以用 投影变换 表示
单应性矩阵H，按照下面的方程映射二维中的点（齐次坐标意义下）:

点的齐次坐标是什么？
点的齐次坐标是依赖于其尺度定义的，
所以， x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1] 都表示同一个二维点。

单应矩阵的自由度是多少呢？
如果给定一个单应H={h_ij}，给它的元素乘上同一个数a，得到的的单应aH和H作用相同，因为新单应无非把齐次点x1变成了齐次点ax1，都是一回事。
因此我们可以把a换成1/h22，那么H就变成了只有8个自由元素的矩阵。

需要多少个点对求解这个H呢？
需要4个点对（对应8个方程，去解H中的8个未知数)求解这个H

单应性矩阵 H 仅依赖尺度定义，所以，单应性矩阵具有 8 个独立自由度。
通常会使第三个值为一来归一化点，即 ：a 点坐标（x , y , z ) 转换成 a’ （x , y , 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]))))

二、直接线性变换算法:

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

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

或者Ah=0，其中 A 是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用 SVD（Singular Value Decomposition， 奇异值分解）算法找到 H 的最小二乘解。

SVD（Singular Value Decomposition， 奇异值分解）算法找到 H 的最小二乘解

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]

函数对这些点进行归一化操作,使其均值为 0，方差为 1。然后使用对应点对来构造矩阵 A。
矩阵 SVD 分解后所得矩阵 V 的最后一行为最小二乘解该行经过变形后得到矩阵 H。然后对这个矩阵进行处理和归一化，返回输出。

为什么要进行归一化操作?
因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题。

三、仿射变换

是线性变换和平移变换的叠加，可用于图像扭曲变形和图像配准。

仿射变换的性质：

1. 仿射变换只有6个自由度（对应变换中的6个系数），因此，仿射变换后互相平行的直线仍是互相平行，三角形映射后也仍是三角形，但却不能保证将四边形以上的多边形映射为等边数的多边形。
2. 仿射变换的乘积和逆变换仍是仿射变换。
3. 仿射变换能够实现平移、旋转、缩放等几何变换。

仿射变换是线性变换后进行平移变换（其实也是齐次空间的线性变换），使用齐次坐标使得仿射变换可以以统一的矩阵形式进行表示。

下面的函数使用对应点对来计算仿射变换矩阵，将其添加到 homograph.py 文件中：

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 包来简单完成。
命令： transformed_im = ndimage.affine_transform(im,A,b,size)

A-线性变换，b-平移向量，size-指定输出图像的大小(默认输出图像设置为和原始图像同样大小)

代码演示:

# -*- codeing =utf-8 -*-
# @Time : 2021/4/6 10:32
# @Author : ArLin
# @File : demo1图像扭曲.py
# @Software: PyCharm
# coding: utf-8
from scipy import ndimage
from PIL import Image
from pylab import *


im = array(Image.open(r'06.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]))


figure()
gray()
subplot(121)
imshow(im)
subplot(122)
imshow(im2)
show()

结果展示:

左图为原图右图为使用 ndimage.affine_transform() 函数扭曲后的图像
通过对比可以得出H针对只有一个平面的图像视觉上来看更加扭曲，
而对于多个平面的图像扭曲的程度比较微弱

五、图像中的图像

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

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

代码演示:

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

我们简单地为每个三角形创建了 alpha 图像，然后将所有的图像合并起来。
该三角形的 alpha 图像可以简单地通过检查像素的坐标是否能够写成三角形顶点坐标的凸 组合来计算得出 1 。如果坐标可以表示成这种形式，那么该像素就位于三角形的内部。

alpha_for_triangle() 函数的代码如下:

def alpha_for_triangle(points,m,n):
#对于由角点定义的三角形创建大小为（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
    #返回alpha的值
    return alpha

5.1 完全仿射变换

使用完全仿射变换将一幅图像放置到另一幅图像中
代码演示:

# -*- codeing =utf-8 -*-
# @Time : 2021/4/6 10:39
# @Author : ArLin
# @File : demo2.py
# @Software: PyCharm
# -*- coding: utf-8 -*-
# -*- coding: utf-8 -*-


from PIL import Image
from pylab import *
from numpy import *
from PCV.geometry import warp


im1 = array(Image.open('04.jpg').convert('L'))
im2 = array(Image.open('05.jpg').convert('L'))


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


# 选定一些目标点
tp = array([[191, 504, 503, 190], [139, 132, 888, 887], [1, 1, 1, 1]])


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

结果展示:



p1映射到p2的目标位置
tp = array([[191, 504, 503, 190], [139, 132, 888, 887], [1, 1, 1, 1]])
从左上角开始按逆时针方向排序
[191, 504, 503, 190]代表四个角点的纵坐标
[139, 132, 888, 887],代表四个角点的横坐标
最后四个数字代表四个角点的α通道，四个1就表示四个角点的透明度为不透明

这些坐标值是通过查看绘制的图像（在 PyLab 图像中，鼠标的坐标显示在图像底部附近）手工确定的。当然，也可以用 PyLab 类 库中的 ginput() 函数获得。

5.2 包含两个三角形的仿射变换

使用包含两个三角形的仿射变换可以很好地将图像完全放置到公告牌 上 对于三个点，仿射变换可以将一幅图像进行扭曲，使这三对对应点对可以完美地匹配上。

这是因为，仿射变换具有 6 个自由度，三个对应点对可以给出 6 个约束条件 （对于这三个对应点对，x 和 y 坐标必须都要匹配）。所以，如果你真的打算使用仿 射变换将图像放置到公告牌上，可以将图像分成两个三角形，然后对它们分别进行扭曲图像操作

代码演示:

# -*- codeing =utf-8 -*-
# @Time : 2021/4/10 14:50
# @Author : ArLin
# @File : demo5.py
# @Software: PyCharm


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


# 两张图片
im1 = array(Image.open(r'09.jpg').convert('L'))
im2 = array(Image.open(r'08.jpg').convert('L'))
gray()
# 设置映射的目标点
tp = array([[191, 504, 503, 190], [139, 132, 888, 887], [1, 1, 1, 1]])


# 选定 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
imshow(im4)
axis('equal')
axis('off')
show()

结果展示:

对比上下图(上图为完全仿射变换细节图，下图为两个三角形的仿射变换细节图)包含两个三角形的仿射变换可以很好地将图像完全放置到公告牌。

六、分段仿射扭曲

从上面的例子可以看出，三角形图像块的仿射扭曲可以完成角点的精确匹配。三角形图像块越多，则匹配的越精确。分段仿射扭曲是通过给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，然后将图像和另一幅图像的对应标记点扭曲对应。

三角剖分函数

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

分段仿射图像扭曲的通用扭曲函数

def pw_affine(fromim,toim,fp,tp,tri):
    """ Warp triangular patches from an image.
        fromim = image to warp
        toim = destination image
        fp = from points in hom. coordinates
        tp = to points in hom.  coordinates
        tri = triangulation. """
                
    im = toim.copy()
    
    # check if image is grayscale or color
    is_color = len(fromim.shape) == 3
    
    # create image to warp to (needed if iterate colors)
    im_t = zeros(im.shape, 'uint8')
    
    for t in tri:
        # compute affine transformation
        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 for triangle
        alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
        
        # add triangle to image
        im[alpha>0] = im_t[alpha>0]
        
    return im

绘制三角形函数

def plot_mesh(x,y,tri):    
""" Plot triangles. """        
for t in tri:        
         t_ext = [t[0], t[1], t[2], t[0]]         
         plot(x[t_ext],y[t_ext],'r')     

多个三角形的选取匹配后的结果比两个三角形匹配的结果要精确，而且相对于取点，如果取点均匀有序的话，效果更佳，理论上来说也是这样。(但是由于本人的手笨，也只能展示到这个程度了)