【python】图像映射:单应性变换与图像扭曲
- 单应性变换(Homography)
- 图像扭曲(仿射变换)
- 图中图
- 分段仿射扭曲
单应性变换(Homography)
单应性变换(Homography)即一个平面中的点到另一个平面的映射关系,如下图所示,使用单应性矩阵,将不同角度拍摄的图像图1和图2变换到同一平面。
而实现单应性变化的重点就是对单应性矩阵H的求解。如下图所示,其中x1y1表示图像一的矩阵数据,x2y2表示图像二的矩阵数据,而H则为单应性矩阵,通过求解H即可得出两图像间的映射关系。
在一般情况下,矩阵A不可逆,若需要求得t的值(下图中t即为单应性矩阵)则需要做一定的数学变换
如下公式所示,进行转换后,通过计算最小二乘解求得最小特征值对应的特征向量
再通过DLT(Direct Linear Transform)算法解决该最小二乘问题,如下图,通过DTL算法得A = UDVt, 而h (单应性矩阵)即为V的最后一行
单应性变换的结果诶下图,图中将图1与图2通过单应性变换进行了图片融合
图像扭曲(仿射变换)
图像扭曲即仿射变换。仿射变换变化包括缩放、平移、旋转、反射等,,原来的直线仿射变换后还是直线,原来的平行线经过仿射变换之后还是平行线。
仿射变换中一些性质会保持不变:
1.凸性
2.共线性:若几个点变换前在一条线上,则仿射变换后仍然在一条线上
3.平行性:若两条线变换前平行,则变换后仍然平行
4.共线比例不变性:变换前一条线上两条线段的比例,在变换后比例不变
from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
im = array(Image.open('11.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()
图中图
如下图所示,将一图融入了二图中,虽然图片有些扭曲但这就是该算法所需要呈现的效果。图像中的图像即为仿射扭曲的一个简单例子,它将图像或者图像的一部分放置在另一幅图像中,使得他们能够和指定的区域或者标记物对齐。为了能够实现更加精准的匹配,我们可以使用分段仿射扭曲。
from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
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('31.jpg').convert('L'))
im2 = array(Image.open('32.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()
分段仿射扭曲
为了能够实现图像中角点的精确匹配,我们通常使用分段仿射扭曲。即首先通过三角化图像块再使用仿射扭曲来扭曲每个三角形,其中的核心部分即为图像块的三角化,可以通过狄洛克三角剖分实现,其具体如下:
from PIL import Image
from pylab import *
import numpy as np
from scipy.spatial import Delaunay
x,y = array(np.random.standard_normal((2,100)))
"""centers, edges, tri, neighbors = md.delaunay(x, y"""
tri = Delaunay(np.c_[x,y]).simplices
figure()
for t in tri:
t_ext = [t[0],t[1],t[2],t[0]]
plot(x[t_ext], y[t_ext], 'g')
plot(x,y,'*')
axis('off')
show()