图像到图像的映射
单应性变换
基本概念
单应性变换描述了平面上的两幅图像之间的映射关系。举个例子,缩放图像和旋转等都是单应性变换。
从数学上来讲,在齐次坐标意义下,单应性变换的变换如下:
[
x
′
y
′
1
]
=
[
h
11
h
12
h
13
h
21
h
22
h
23
h
31
h
32
h
33
]
[
x
y
1
]
\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&h_{33}\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}
x′y′1
=
h11h21h31h12h22h32h13h23h33
xy1
或者说
x
′
=
H
x
x' = Hx
x′=Hx
这里的齐次坐标是依赖于尺度定义的(也就是上式中的1)。正常情况下定义点的坐标只需要定义(x,y)即可,但是此处我们增加了尺度定义可以更方便的定义一些变换。
比方说平移变换和尺度变换原本需要两步来做:
平移变换:
T
=
[
1
0
t
x
0
1
t
y
0
0
1
]
T=\begin{bmatrix}1&0&t_x\\0&1&t_y\\0&0&1\end{bmatrix}
T=
100010txty1
尺度变换:
S
=
[
s
0
0
0
s
0
0
0
1
]
S=\begin{bmatrix}s&0&0\\0&s&0\\0&0&1\end{bmatrix}
S=
s000s0001
现在可以直接合并:
[
x
′
y
′
1
]
=
[
s
0
t
x
0
s
t
y
0
0
1
]
[
x
y
1
]
\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}s&0&t_x\\0&s&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}
x′y′1
=
s000s0txty1
xy1
此外补充一点:单应性变换矩阵 H H H中的 h 33 h_{33} h33通常定义为1,因此这个矩阵只有8个独立的自由度,因此只需要4个点对就可以解出原方程(也就是说,假如我们已知两张图上的四个点对,我们就可以知道这两张图是怎么变换的了)
直接线性变换算法
DLT(Direct Linear Transformation,直接线性变换)是给定一些点对,来计算单应性矩阵 H 的算法。
我们前面提到过,只需要4个点对即可求出变换矩阵,具体解法此处不再赘述,直接给出代码:
def H_from_points(fp, tp):
""" Find homography H, such that fp is mapped to tp
using the linear DLT method. Points are conditioned
automatically. """
if fp.shape != tp.shape:
raise RuntimeError('number of points do not match')
# condition points (important for numerical reasons)
# --from points--
m = np.mean(fp[:2], axis=1)
maxstd = max(np.std(fp[:2], axis=1)) + 1e-9
C1 = np.diag([1/maxstd, 1/maxstd, 1])
C1[0][2] = -m[0]/maxstd
C1[1][2] = -m[1]/maxstd
fp = np.dot(C1, fp)
# --to points--
m = np.mean(tp[:2], axis=1)
maxstd = max(np.std(tp[:2], axis=1)) + 1e-9
C2 = np.diag([1/maxstd, 1/maxstd, 1])
C2[0][2] = -m[0]/maxstd
C2[1][2] = -m[1]/maxstd
tp = np.dot(C2, tp)
# create matrix for linear method, 2 rows for each correspondence pair
nbr_correspondences = fp.shape[1]
A = np.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 = np.linalg.svd(A)
H = V[8].reshape((3, 3))
# decondition
H = np.dot(np.linalg.inv(C2), np.dot(H, C1))
return H/H[2, 2]
仿射变换
仿射变换(Affine Transformation)是指在二维平面上进行的一类几何变换,它保持了点的共线性、平行性和比例关系。在仿射变换中,原始图像上的所有平行线在变换后仍然保持平行,并且直线上的点在变换后仍然共线。
其主要分类有:
- 平移
- 旋转
- 缩放
- 错切
简单来说,仿射变换就是允许图像任意倾斜,而且允许图形在两个方向上任意伸缩变换,但是,不能保证原来的线段长度不变,也不能保证原来的夹角角度不变。
下面直接上代码:
def Haffine_from_points(fp, tp):
""" Find H, affine transformation, such that
tp is affine transf of fp. """
if fp.shape != tp.shape:
raise RuntimeError('number of points do not match')
# condition points
# --from points--
m = np.mean(fp[:2], axis=1)
maxstd = max(np.std(fp[:2], axis=1)) + 1e-9
C1 = np.diag([1/maxstd, 1/maxstd, 1])
C1[0][2] = -m[0]/maxstd
C1[1][2] = -m[1]/maxstd
fp_cond = np.dot(C1, fp)
# --to points--
m = np.mean(tp[:2], axis=1)
C2 = C1.copy() # must use same scaling for both point sets
C2[0][2] = -m[0]/maxstd
C2[1][2] = -m[1]/maxstd
tp_cond = np.dot(C2, tp)
# conditioned points have mean zero, so translation is zero
A = np.concatenate((fp_cond[:2], tp_cond[:2]), axis=0)
U, S, V = np.linalg.svd(A.T)
# create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
tmp = V[:2].T
B = tmp[:2]
C = tmp[2:4]
tmp2 = np.concatenate((np.dot(C, np.linalg.pinv(B)), np.zeros((2, 1))), axis=1)
H = np.vstack((tmp2, [0, 0, 1]))
# decondition
H = np.dot(np.linalg.inv(C2), np.dot(H, C1))
return H / H[2, 2]
Haffine_from_points 与H_from_points 非常相似,这两个函数的主要区别在于它们计算的变换类型不同。H_from_points 函数计算的是单应性矩阵(homography),而 Haffine_from_points 函数计算的是仿射变换矩阵(affine transformation)。单应性矩阵是一种更为一般的变换,可以将一个平面上的点映射到另一个平面上的点,而仿射变换则是单应性矩阵的一种特殊情况,它只能进行平移、旋转、缩放和错切这些基本变换。
具体到代码层面的话,H_from_points 函数通过使用线性 DLT 方法计算单应性矩阵,而 Haffine_from_points 函数则通过使用 SVD 分解计算仿射变换矩阵。此外,这两个函数在对输入点进行归一化时使用的变换矩阵也不同,因为单应性矩阵和仿射变换矩阵的计算方式不同。
图像扭曲
在本节中,我们将主要介绍仿射变换是如何处理图像的。上面的代码仅用于展示如何实现,本节中我们将采用封装好的方法来操作。
import numpy as np
import cv2
import matplotlib.pyplot as plt
# 仿射变换
img = cv2.imread('img.jpg')
rows, cols, ch = img.shape
pts1 = np.float32([[50, 50], [200, 50], [50, 200]])
pts2 = np.float32([[10, 100], [200, 50], [100, 250]])
M = cv2.getAffineTransform(pts1, pts2)
dst = cv2.warpAffine(img, M, (cols, rows))
# 转换通道
b, g, r = cv2.split(img)
img = cv2.merge((r, g, b))
b, g, r = cv2.split(dst)
dst = cv2.merge((r, g, b))
# 规定大小
plt.subplot(121), plt.imshow(img), plt.title('Input')
plt.subplot(122), plt.imshow(dst), plt.title('Output')
plt.show()
运行结果如下图所示:
可以看到输出这里图像被扭曲了,黑色部分是由于被0填充。从代码中我们也可以很清楚的看到,这个仿射变换是由三个点对求得的。
图像中的图像
仿射扭曲的一个简单例子是,将图像或者图像的一部分放置在另一幅图像中,使得它们能够和指定的区域或者标记物对齐。
这么讲可能不够直观,看书上给的例子:
代码:
import numpy as np
import cv2
import matplotlib.pyplot as plt
# 将im1仿射扭曲到im2的指定位置
def warp_image_in_image(im1, im2, tp, tp1):
# 计算仿射变换矩阵
H = cv2.getAffineTransform(tp, tp1)
# 仿射变换
im1_warp = cv2.warpAffine(im1, H, (im2.shape[1], im2.shape[0]))
# 将im1扭曲到im2的指定位置
im2_warp = im2.copy()
im2_warp[im1_warp>0] = im1_warp[im1_warp>0]
return im2_warp
# 读取图像
im1 = cv2.imread('im1.jpg')
im2 = cv2.imread('im2.jpg')
# 设置仿射变换的目标点
tp = np.array([[0, 0], [0, im1.shape[0]], [im1.shape[1], 0]], np.float32)
tp1 = np.array([[0, 0], [0, im2.shape[0]/2], [im2.shape[1]/2, 0]], np.float32)
# 将im1扭曲到im2的指定位置
im_warp = warp_image_in_image(im1, im2, tp, tp1)
# 显示结果
plt.subplot(131)
plt.imshow(im1[:,:,::-1])
plt.subplot(132)
plt.imshow(im2[:,:,::-1])
plt.subplot(133)
plt.imshow(im_warp[:,:,::-1])
plt.show()
运行结果:
到目前为止,这个东西似乎还并不是那么有用,不过如果这两张图是有对应点对的话,那么是不是可以去拼接两张图呢?我觉得是完全可以的。
比方说如果有个图片中不动的电视机,那么我们完全可以通过这种手段来修改电视机上的内容,想象空间还是很大的。
图像配准
图像配准是对图像进行变换,使变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准,也可以是非严格配准。为了能够进行图像对比和更精细的图像分析,图像配准是一步非常重要的操作。
上面这一段话是书中的原文,但是比较难以理解,这里来举一个例子吧:现在我们有很多很多张人脸的图片,但是他们在图片上的坐标不同,比方说第一张在图片右侧,第二张在图片左侧,而且图片的尺寸也不相同。但是我们想要制作一个标准的数据集,这个时候就需要对图片进行裁剪旋转等一系列操作,使得图片大小一致,人脸在图片中的大小一致,方位大致相同(这几个条件未必是必须的,根据任务灵活处理)
更进一步,如果我们有两张关于同一个物体的图片,他们的角度相同,大小不同。我们通过一些操作,比如放大之类的操作让其看起来大小一致,这也可以算是一种图像配准。示例如下,把人脸在图中的位置给统一了:
图像配准一般需要经过以下步骤:
- 特征提取
- 特征匹配
- 估计变换模型
- 图像重采样以及变换
这几步我们通过前三章的学习已经基本学完了,像是特征提取与特征匹配在第二章中都已经接触过,如何变换在本章也已经介绍过,如果再做实验也无非是二者组合,因此本段不再进行实验,直接进入到下一节。
创建全景图
在同一位置(即图像的照相机位置相同)拍摄的两幅或者多幅图像是单应性相关的 (也即是平面到平面的变换)。我们经常使用该约束将很多图像缝补起来,拼成一个大的图像来创建全景图像。
RANSAC
RANSAC 是“RANdom SAmple Consensus”(随机一致性采样)的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型,例如点集之间的单应性矩阵,RANSAC 基本的思想是,数据中包含正确的点和噪声点,合理的模型应该能够在描述正确数据点的同时摒弃噪声点。
其基本步骤是:
- 从数据集中随机选择一个最小的样本集,并使用这些样本点来估计模型的参数。
- 计算所有数据点与估计的模型之间的误差,并找到在给定阈值内的内点(inliers),即与模型拟合得比较好的数据点。
- 如果内点数量大于阈值(一般情况下,是事先设定的一个阈值),则使用内点重新估计模型参数。
- 重复步骤1-3,迭代固定次数或直到满足停止条件为止。
- 最终,使用所有内点(或所有迭代中获得的内点最多的模型)来估计最终的模型参数。
创建全景图
其基本步骤如下:
- 加载图像
- 特征检测与匹配
- 计算透视变换矩阵(RANSAC算法在这里起作用,选择最正确的拟合)
- 图像拼接
- 显示全景图
下面给出了创建全景图的完整示例代码:
import numpy as np
import cv2
import matplotlib.pyplot as plt
# 加载图像
image1 = cv2.imread('m.jpg')
image2 = cv2.imread('l.jpg')
# 使用SIFT特征检测器找到关键点和特征描述符
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(image1, None)
kp2, des2 = sift.detectAndCompute(image2, None)
# 使用FLANN匹配器进行特征匹配
FLANN_INDEX_KDTREE = 1
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des1, des2, k=2)
# 应用比值测试来剔除错误匹配
good_matches = []
for m, n in matches:
if m.distance < 0.7 * n.distance:
good_matches.append(m)
# 获取关键点对应的坐标
src_pts = np.float32([kp1[m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)
dst_pts = np.float32([kp2[m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)
# 计算透视变换矩阵
M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)
# 使用透视变换矩阵进行图像拼接
result = cv2.warpPerspective(image1, M, (image1.shape[1] + image2.shape[1], image1.shape[0]))
result[0:image2.shape[0], 0:image2.shape[1]] = image2
# 显示全景图
plt.imshow(cv2.cvtColor(result, cv2.COLOR_BGR2RGB))
# plt.axis('off')
plt.show()
运行结果如下:
原始图像此处不再给出,可以看到上图中有明显的色差,这应该是两张图拼接在了一起的原因。同时右侧的黑色部分是因为我们选择将两张图片扭曲到一张黑色基地的图片上的原因,除此之外,横坐标4000附近的直线是倾斜的,这应该是由于单映射变换时将原本的图片扭曲了的原因。
小结
这一章主要学习了单应性变换,即怎么将图片按照我们想要的方式去扭曲,这里其实是有非常多能做的事情,比方说我可以拿去换脸这样子,选取几个点然后把另一个人的脸去覆盖掉。但是这是非常传统的方法了,深度学习似乎并不是采用的这种方法,深度学习是直接重绘;那么,这种方法能否与深度学习结合呢?比方说将这个函数提供给神经网络去调用这样子,感觉不应该没人做,估计是我看的论文太少了吧。