python计算机视觉——第三章图像到图像的映射

ZLiu66

已于 2024-09-01 13:12:21 修改

阅读量187

点赞数 3

文章标签：计算机视觉人工智能

于 2024-08-24 20:37:04 首次发布

本文链接：https://blog.csdn.net/weixin_46002351/article/details/141504564

版权

本章讲解图像之间的变换，以及一些计算变换的实用方法。这些变换可以用于图像扭曲变形和图像配准。

3.1 单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换，这里的平面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如说图像的配准、图像纠正和纹理扭曲，以及创建全景图像等。本质上，单应性变换H，按照下面的方程映射二维中的点：

$\begin{bmatrix}x'\\y'\\w'\end{bmatrix}=\begin{bmatrix}h_1&h_2&h_3\\h_4&h_5&h_6\\h_7&h_8&h_9\end{bmatrix}\begin{bmatrix}x\\y\\w\end{bmatrix}$或$\mathbf{x}^{\prime}=H\mathbf{x}$

对于图像平面内的点，齐次坐标是一个非常有用的表达方式。点的齐次坐标是依赖于其尺度定义的，所以$x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]$，表示的是同一个二维点。单应性矩阵仅依赖尺度定义，具有8个独立的自由度。

创建一个能够实现对点进行归一化和转换齐次坐标的功能：

def normalize(points):
    """
    在其次坐标意义下对点集进行归一化处理
    :param points: 点集
    :return:
    """
    for row in points:
        row /= points[-1]
    return points
 
def make_homog(points):
    """
    将点集转化为齐次坐标表示
    :param points: 点集
    :return:
    """
    return vstack((points, ones((1, points.shape[1]))))

进行点和变换的处理时，按照列优先的原则存储这些点，故n个三维点集将会依次存储为齐次坐标意义下的一个3*n数组。

在这些投影变换中存在一些特别重要的变换如：

仿射变换：

$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}a_1&a_2&t_x\\a_3&a_4&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}$或$x'=\begin{bmatrix}A&t\\0&1\end{bmatrix}x$

w=1，不具有投影变换所具有的强大变形能力。仿射变换包含一个可逆矩阵A和一个平移向量$t[t_{x},t_{y}]$。

相似变换：

$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}s\cos(\theta)&-s\sin(\theta)&t_x\\s\sin(\theta)&s\cos(\theta)&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}$或$x'=\begin{bmatrix}sR&t\\0&1\end{bmatrix}x$

包含尺度变换的二维刚体变换，上式的s向量指定了变换的尺度，R是角度为 $\theta$ 的旋转矩阵，$t=[t_{x},t_{y}]$是一个平移向量，s=1时，该变换甭管保存距离不变，此时变换为刚体变换。

3.1.1 直接线性变换算法

单应性矩阵可以由两幅图像中对于点计算出来，我们已经了解到一个完全射影变换具有8个自由度，而每个对应点可以写出两个方程，对应x和y坐标，所以计算单应性矩阵H要4个对应点对。计算H的方程如下：

$\begin{bmatrix}-x_1&-y_1&-1&0&0&0&x_1x_1'&y_1x_1'&x_1'\\0&0&0&-x_1&-y_1&-1&x_1y'_1&y_1y'_1&y'_1\\-x_2&-y_2&-1&0&0&0&x_2x'_2&y_2x'_2&x'_2\\0&0&0&-x_2&-y_2&-1&x_2y'_2&y_2y'_2&y'_2\\&\vdots&&\vdots&\vdots&&\vdots\end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\\h_4\\h_3\\h_6\\h_7\\h_8\\h_9\end{bmatrix}=\mathbf{0}$

def H_from_points(fp, tp):
    """
    使用使用DLT方法，计算单应性矩阵H，使fp映射到tp，点自动进行归一化
    :param fp:
    :param tp:
    :return:
    """
    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_correrpondences = fp.shape[1]
    A = zeros(2 * nbr_correrpondences, 9)
    for i in range(nbr_correrpondences):
        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]

函数的第一步操作是检查点对两个数组中点的数目是否相同，如果不相同函数将会给出错误提示，对这些点进行归一化操作，使其均值为，方差为，因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题，因此归一化非常重要，之后利用对应点对构造矩阵A，最小二乘解即为矩阵SVD分解后得到矩阵V的最后一行。该行经过变形后得到矩阵H，并对得到的矩阵进行处理和归一化，返回输出。

3.1.2 仿射变换

先前提到过仿射变换是一种重要的投影变换，其具有6个自由度，因此需要3个对应点估计矩阵H，将最后两个元素设置为0，就可以使用DLT算法估计得出。

def Haffine_from_points(fp, tp):
    """
    计算H，仿射变换，使得tp是fp经过仿射变换H得到的
    :param fp: 
    :param tp: 
    :return: 
    """
    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)
    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_cond = dot(C2, tp)
    
    # 因为归一化后的点均值为0，故平移量为0
    A = concatenate((fp_cond[:2], tp_cond[:2]), axis=0)
    U, S, A = linalg.svd(A.T)
    
    # 创建矩阵B,C
    temp = V[:2].T
    B = temp[:2]
    C = temp[2:4]
    temp2 = concatenate((dot(C, linalg.pinv(B)), zeros((2, 1))), axis=1)
    H = vstack((temp2, [0, 0, 1]))
    
    # 反归一化
    H = dot(linalg.inv(C2), dot(H,C1))
    return H / H[2, 2]

3.2　图像扭曲

对图像应用仿射变换就将其称为图像扭曲，扭曲操作可使用Scipy工具包里的ndimage包实现，具体命令为：

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

# 图像扭曲
from scipy import ndimage
from PIL import Image
from numpy import * 
from pylab import *
 
img = array(Image.open('ch03\\3.2\\1.jpg').convert('L'))
gray()
H = array([[1.4, 0.05, -100], [0.05, 1.5, -100], [0, 0, 1]])
img_re = ndimage.affine_transform(img, H[:2, :2], (H[0, 2], H[1, 2]))
subplot(121)
imshow(img), title('original')
subplot(122)
imshow(img_re), title('result')
show()

3.2.1 图像中的图像

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

from PIL import Image
from numpy import *
from pylab import *
import warp

# 仿射扭曲im1 到im2 的例子
im1 = array(Image.open('ch03\\3.2\\ren.jpg').convert('L'))
im2 = array(Image.open('ch03\\3.2\\s.jpg').convert('L'))
# 选定一些目标点
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
im3 = warp.image_in_image(im1,im2,tp)
figure()
gray()
imshow(im3)
axis('equal')
axis('off')
show()

将扭曲的图像和第二幅图像融合，创建 alpha 图像。该图像定义了每个像素从各个图像中获取的像素值成分多少。这里基于以下事实，扭曲的图像是在扭曲区域边界之外以 0 来填充的图像，来创建一个二值的 alpha 图像。严格意义上说，需要在第一幅图像中的潜在 0 像素上加上一个小的数值，或者合理地处理这些 0 像素。

3.2.2 分段仿射扭曲

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准，也可以是非严格配准。为了能够进行图像对比和更精细的图像分析，图像配准是一步非常重要的操作。

配准算法的一般步骤：

特征提取
特征匹配
估计变换模型
图像重采样及变换

特征提取

特征提取是指分别提取两幅图像中共有的图像特征，这种特征是出现在两幅图像中对比列、旋转、平移等变换保持一致性的特征，如线交叉点、物体边缘角点、虚圆闭区域的中心等可提取的特征。特征包括：点、线和面三类。

特征匹配

特征匹配分为两步：

对特征作描述
利用相似度准则进行特征匹配

估计变换模型：

空间变换模型是所有配准技术中需要考虑的一个重要因素，各种配准技术都要建立自己的变换模型，变换空间的选取与图像的变形特性有关。常用的空间变换模型有：刚体变换、仿射变换、投影变换、非线性变换。

图像重采样及变换：

在得到两幅图像的变换参数后，要将输入图像做相应参数的变换，使之与参考图像处于同一坐标系下，则矫正后的输入图像与参考图像可用作后续的图像融合、目标变化检测处理或图像镶嵌；涉及输入图像变换后所得点坐标不一定为整像素数，则应进行插值处理。常用的插值算法有最近领域法，双线性插值法和立方卷积插值法

SIFT特征匹配算法包括两个阶段：

SIFT特征的生成
SIFT特征向量的匹配

SIFT特征向量的生成算法包括四步：

尺度空间极值检测，以初步确定关键点位置和所在尺度。
拟和三维二次函数精确确定位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点。
利用关键点领域像素的梯度方向分布特性为每个关键点指定参数方向，使算子具备旋转不变性。
生成SIFT特征向量。

SIFT特征向量的匹配：对图像1中的某个关键点，找到其与图像2中欧式距离最近的前两个关键点的距离NN和SCN，如果NN/SCN小于某个比例阈值，则接受这一对匹配点。

3.3 创建全景图

在同一位置（即图像的照相机位置相同）拍摄的两幅或者多幅图像是单应性相关的。我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建全景图像。

3.3.1 RANSAC

PANSAC是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC基本的思想是，数据中包含着正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

RANSAC的基本假设是：

数据由“局内点”组成，例如：数据的分布可以用一些模型参数来解释；
“局外点”是不能适应该模型的数据；
除此之外的数据属于噪声。

局外点产生的原因有：噪声的极值；错误的测量方法；对数据的错误假设。

RANSAC也做了以下假设：给定一组（通常很小的）局内点，存在一个可以估计模型参数的过程；而该模型能够解释或者适用于局内点。

RANSAC算法的输入是一组观测数据，一个可以解释或者适应于观测数据的参数化模型，一些可信的参数。 RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点，并用下述方法进行验证：

首先我们先随机假设一小组局内点为初始值。然后用此局内点拟合一个模型，此模型适应于假设的局内点，所有的未知参数都能从假设的局内点计算得出。
用1中得到的模型去测试所有的其它数据，如果某个点适用于估计的模型，认为它也是局内点，将局内点扩充。
如果有足够多的点被归类为假设的局内点，那么估计的模型就足够合理。
然后，用所有假设的局内点去重新估计模型，因为此模型仅仅是在初始的假设的局内点估计的，后续有扩充后，需要更新。
最后，通过估计局内点与模型的错误率来评估模型。

整个这个过程为迭代一次，此过程被重复执行固定的次数，每次产生的模型有两个结局：

要么因为局内点太少，还不如上一次的模型，而被舍弃。
要么因为比现有的模型更好而被选用。

3.3.2 稳健的单应性矩阵估计

我们在任何模型中都可以使用 RANSAC 模块。在使用 RANSAC 模块时，我们只需要在相应 Python 类中实现 fit() 和 get_error() 方法，剩下就是正确地使用 ransac.py 。我们这里使用可能的对应点集来自动找到用于全景图像的单应性矩阵。

from numpy import *
from PIL import Image
from PIL import Image
from pylab import *
from PCV.localdescriptors import sift

import sift

imname = [r'ch03\3.3\Univ' + str(i + 1) + '.jpg' for i in range(5)]
featname = [r'ch03\3.3\Univ' + str(i + 1) + '.sift' for i in range(5)]
im = [array(Image.open(imname[i]).convert('L')) for i in range(5)]
# 提取特征并匹配使用sift算法
l = {}
d = {}

# 循环次数 = 图像数量
for i in range(5):
    sift.process_image(imname[i],featname[i])
    l[i], d[i] = sift.read_features_from_file(featname[i])
 
matches = {}

# 循环次数 = 图像数量 - 1
for i in range(4):
    matches[i] = sift.match(d[i+1],d[i])

 # 循环次数 = 图像数量 - 1
for i in range(4):
    figure()
    gray()
    sift.plot_matches(im[i+1], im[i], l[i+1], l[i], matches[i], show_below=True)

show()

SIFT 是具有很强稳健性的描述子，能够比其他描述子，例如图像块相关的Harris 角点，产生更少的错误的匹配。但是该方法仍然远非完美。

3.3.3 拼接图像

估计出图像间的单应性矩阵（使用RANSAC算法），现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面（否则需要大量变形）。一种方法是创建一个很大的图像，比如图像中全部填充0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充为0，以便为扭曲的图像腾出空间。

全景图像拼接步骤：

特征点的匹配：两张图像要能拼接在一起成为一张图像，就需要这两张图像中存在有重合的部分。通过这些重合的部分使用sift特征点匹配的算法，来寻找到重合部分的特征点。需要注意的是，虽然sift算法比Harris角点的效果更好，但是也会出现错误点，并非完美的匹配方法。合成全景图的第一步是提取并且匹配所有素材图片的局部特征点。$\n$什么是特征点？普遍来讲，一张图片所包含的特征点通常就是周围含有较大信息量的点，而仅通过这些富有特征的局部，基本就可以推测出整张图片。比如说物体的棱角、夜景闪耀的星星，或是图片里的图案和花纹。
图像的匹配：在找到特征点对之后，因为上文提到sift并非所有匹配点都是正确的，这里我们用到了RANSAC这个方法。这个方法的作用是找到一个合理的模型来描述正确的数据并且尽量忽视噪点的影响。接下来就是，找到所有匹配(也就是重叠)的图片部分，连接所有图片之后就可以形成一个基本的全景图了。因为每张图片有可能和其他每张图片有重叠部分，所以匹配全部图片需要差不多匹配图片个数的平方次。不过实际上每两张图片之间只需要那么几个相对精准匹配的点就可以估算出这两张图像里的几何关系。
全景图矫直：矫正拍摄图片时相机的相对3D旋转，主要原因是拍摄图片时相机很可能并不在同一水平线上，并且存在不同程度的倾斜，略过这一步可能导致全景图变成波浪形状。
图像均衡补偿：全局平衡所有图片的光照和色调。
图像频段融合：步骤4之后仍然会存在图像之间衔接边缘、晕影效果（图像的外围部分的亮度或饱和度比中心区域低）、视差效果（因为相机透镜移动导致）。