单应性变换
单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里,平面是指图像或者三维中的平面表面。
如果场景是平面,或者近似平面,或者低视差时,我们能应用单应性矩阵(homography)
如下图所示的平面的两幅图像。红点表示两幅图像中的相同物理点,我们称之为对应点。这里显示了四种不同颜色的四个对应点即红色,绿色,黄色和橙色。
单应性变换具有很强的实用性,比如图像配准,图像纠正和纹理扭曲,以及创建全景图像,我们将频繁的使用单应性变换。本质上,单应性变换H,按照下面的方程映射二维中的点(齐次坐标意义下):
或者 x′=Hx
对于图像平面内(甚至是三维中的点,后面我们会介绍到)的点,齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的,所以,x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此,单应性矩阵H也仅依赖尺度定义,所以,单应性矩阵具有8个独立的自由度。我们通常使用w=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]))))
进行点和变换的处理时。我们会按照列优先的原则存储这些点。因此,n个二维点集将会存储为齐次坐标意义下的一个3×n数组。这种格式使得矩阵乘法和点的变换操作更加容易。对于其他的例子,比如对于聚类和分类的特征,我们将使用典型的行数组来存储数据。
在这些投影变换中,有一些特别重要的变换。比如,仿射变换:
或者
保持了w=1,不具有投影变换所具有的强大变形能力,反射变换包括一个可逆矩阵A和一个平移向量t=[tx,ty]。仿射变换可以用于很多应用,比如图像扭曲。
图注:仿射变换(正方形-平行四边形)
相似变换:
或者
是一个包含尺度变化的二维刚体变换。上式中的向量s指定了变换的尺度,R是角度为θ的旋转矩阵,t=[tx,ty]在这里也是一个平移向量。如果s=1,那么该变换能够保持距离不变。此时,变换称为刚体变换。相似变换可以用于很多应用,比如图像配准。
下面让我们来一起探讨如何设计用于估计单应性矩阵的算法,然后看一下使用仿射变换进行图像扭曲,使用相似变换进行图像匹配,以及使用完全投影变换进行创建全景图像的一些例子:
1.直接线性变换算法
直接线性变换解法是建立像点的“坐标仪坐标”和相应物点的物方空间坐标直接的线性关系的解法。
直接线性变换解法的特点:
- 不归心、不定项
- 不需要内外方位元素的起始值
- 物方空间需布置一组控制点
- 特别适合于处理非量测相机所摄影像
- 本质是一种空间后交-欠交解法。
单应性矩阵可以由两幅图像(或者平面)中对应点对计算出来。前面已经提到过,一个完全射影变换具有8个自由度。根据对应点约束,每个对应点对可以写出两个方程,分别对应于x和y坐标。因此,计算单应性矩阵H需要4个对应点对。
DLT(Direct Linear Transformation,直接线性变换)是给定4个点或者更多对应点对矩阵,来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点上,重新写出该方程,我们可以得到下面的方程:
或者Ah=0,其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵红,我们可以使用SVD算法找到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。然后对这个矩阵进行处理和归一化,返回输出。
2.仿射变换
由于仿射变换具有6个自由度,因此我们需要三个对应点来估计矩阵H。通过将最后两个元素设置为0,即h7=h8=0,仿射变换可以用上面的DLT算法估计得出。
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]
图像扭曲
对图像块应用仿射变换,我们将其称为图像扭曲(或者仿射扭曲)。
仿射变换其实就是将图片上每一个像素点按照一定的规律映射到新的位置上,实际上就是求解新的x,y的过程,这跟矩阵仿射变换非常像,仿射变换的实质其实就是将原图像的三个点映射到目标图片三个新的位置上,这三个位置对于原图片来说位于左上角,左下角以及右上角,通过把原图片上三个点映射到目标图片三个新的位置上。
扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。
代码实现
from numpy import *
from matplotlib.pyplot import *
from numpy import array
from scipy import ndimage
from PIL import Image
im = array(Image.open('C:\\Users\\Administrator\\PycharmProjects\\pythonProject\\VC\\img\\listen.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()
运行结果
1.图像中的图像
仿射扭曲的一个简单例子是,将图像或者图像的一部分放置在另一幅图像中,是的他们能够和指定的区域或者标记物对齐。
代码实现
from numpy import *
from matplotlib.pyplot import *
from numpy import array, dot, concatenate, zeros, diag
from numpy.ma import vstack
from scipy import ndimage, linalg
from PIL import Image
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]
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 = 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\\Administrator\\PycharmProjects\\pythonProject\\VC\\img\\miao.jpg').convert('L'))
im2 = array(Image.open('C:\\Users\\Administrator\\PycharmProjects\\pythonProject\\VC\\img\\2.jpg').convert('L'))
gray()
subplot(131)
imshow(im1)
axis('equal')
axis('off')
subplot(132)
imshow(im2)
axis('equal')
axis('off')
# 选定一些目标点
tp = array([[264, 538, 540, 264], [40, 36, 605, 605], [1, 1, 1, 1]])
im3 = image_in_image(im1, im2, tp)
subplot(133)
imshow(im3)
axis('equal')
axis('off')
show()