【懒懒的计算机视觉笔记之单应性变换】

最近一直在学习Python计算机视觉编程中图像到图像之间的映射，这些变化可以用于图像扭曲变形和图像配准。所谓单应性变换就是将一个平面内的点映射到另一个平面内的二维投影变换。单应性变换具有很强的实用性，比如图像配准、图像纠正和纹理扭曲，以及创建全景图像等。其实，单应性变换H，是按照下面的公式映射二维中的点（其次坐标意义下）：

                                                                


齐次坐标是在实际应用中非常有用的表达式，点的齐次坐标依赖其尺度定义，所以x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都可以表示同一个二维点，所以单应性矩阵H也仅仅依赖尺度定义，单应性矩阵H有8个独立的自由度。我们通常使用w=1来归一化点，这样点具有唯一的图像坐标x和y。
下面的代码实现了对点进行归一化和转换成齐次坐标：

# coding=UTF-8
from PIL import Image
from numpy import *
from pylab import *

def make_homog(points):
'''
:param points:
:return:
'''
return vstack((points, ones((1, points.shape[1]))))

def normalize(points):
'''
在齐次坐标下进行归一化处理
:param points:
:return:
'''
for row in points:
row /= points[-1]
return points单应性矩阵可以通过对应点计算出来，一个完全的单应性矩阵有8个自由度，所以需要4个对应点对计算单应性矩阵，下面我们使用SVD算法找到H的最小二乘解：
def Hfrompoints(fp, tp):
'''
求解单应性矩阵，使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[0] / 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[0] / maxstd
tp = dot(C1, tp)

# 求解 H
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]