图像到图像的映射
- 一、单应性变换矩阵
- 1. 直接线性变换算法
- 2. 仿射变换矩阵
- 二、图像扭曲
- 三、图像中的图像
- 四、使用三角形仿射弯曲效果
一、单应性变换矩阵
单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表示。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，我们将频繁的使用单应性变换。本质上，单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：
或者
对于图像平面内（甚至是三维中的点，后面我们会介绍到）的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。我们通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标是的我们可以简单地使用一个矩阵来表示变换。
下面让我们来一起探讨如何设计用于估计单应性矩阵的算法，然后看一下使用仿射变换进行图像扭曲，使用相似变换进行图像匹配，以及使用完全投影变换进行创建全景图像的一些例子：
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]2. 仿射变换矩阵
仿射变换是一种二维坐标到二维坐标之间的线性变换，保持二维图形的“平直性”（译注：straightness，即变换后直线还是直线不会打弯，圆弧还是圆弧）和“平行性”（译注：parallelness，其实是指保二维图形间的相对位置关系不变，平行线还是平行线，相交直线的交角不变。）。
简单说来，仿射变换，就是允许图形任意倾斜，而且允许图形在两个方向上任意伸缩的变换（如下图)。
仿射变换，可以保持原来的线共点、点共线的关系不变，保持原来相互平行的线仍然平行，保持原来的中点仍然是中点，保持原来在一直线上几段线段之间的比例关系不变。
但是，仿射变换不能保持原来的线段长度不变，也不能保持原来的夹角角度不变。
下面是计算仿射矩阵H的方法。其中tp是变换后的坐标，fp是变换前的坐标，通过计算H，使得tp是fp通过仿射变换矩阵H得到的，然后返回H。
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]二、图像扭曲
对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常在计算机图形学中，而且经常出现在计算机视觉算法总。扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。
图像扭曲代码：
from array import array
from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
im = array(Image.open('../data/f.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()运行结果：
三、图像中的图像
将一张图像放到另一张图像中，这是仿射扭曲一个比较简单的应用。通过图像扭曲的特点，它能够使一张图像的4个顶点和另一幅图像的指定位置对齐。
方法image_in_image代码：
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
H = homography.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调用image_in_image方法：
from PCV.geometry.warp import image_in_image
from array import array
from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
im1 = array(Image.open('../data/sift/xx02.jpg').convert('L'))
im2 = array(Image.open('../data/sift/xx10.jpg').convert('L'))
gray()
subplot(131)
imshow(im1)
axis('equal')
axis('off')
subplot(132)
imshow(im2)
axis('equal')
axis('off')
# 选定一些目标点
tp = array([[264, 800, 800, 264], [40, 40, 400, 400], [1, 1, 1, 1]])
im3 = image_in_image(im1, im2, tp)
subplot(133)
imshow(im3)
axis('equal')
axis('off')
show()运行结果：
四、使用三角形仿射弯曲效果
使用三角形仿射，将第一张图从左上角到右下角分割成两个三角形，然后先将左下角的三角形先贴到公告牌上，再将右上角的三角形贴到公告牌上，这样就能够使得公告牌被贴满。
三角形仿射代码：
from numpy import *
from matplotlib.pyplot import *
from numpy import array
from scipy import ndimage
from PIL import Image
import alpha
import affinity
from JiSuanJiShiJue.shiyan3 import warp
im1 = array(Image.open('../data/sift/photo9.jpg').convert('L'))
im2 = array(Image.open('../data/img.png').convert('L'))
# 选定im1角上的一些点
print(im1.shape)
m, n = im1.shape[:2]
# 中山纪念馆的4个顶点坐标
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])
# 映射倒公告牌中的坐标
tp = array([[50, 288, 275, 33], [60, 50, 480, 475], [1, 1, 1, 1]])
# 第一个三角形
tp2 = tp[:, :3] # 前三个点
fp2 = fp[:, :3]
# 计算H
H = affinity.Haffine_from_points(tp2, fp2)
im1_t = ndimage.affine_transform(im1, H[:2, :2],
(H[0, 2], H[1, 2]), im2.shape[:2])
# 三角形的alpha
alpha = alpha.alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im3 = (1 - alpha) * im2 + alpha * im1_t
# 第二个三角形
tp2 = tp[:, [0, 2, 3]] # 第1、3、4个坐标
fp2 = fp[:, [0, 2, 3]]
# 计算H
H = affinity.Haffine_from_points(tp2, fp2)
im1_t = ndimage.affine_transform(im1, H[:2, :2],
(H[0, 2], H[0, 2]), im2.shape[:2])
# 三角形的alpha图像
alpha = warp.alpha_for_triangle(tp2, im2.shape[0], im2.shape[1])
im4 = (1 - alpha) * im3 + alpha * im1_t
figure()
gray()
imshow(im4)
axis("equal")
# axis("off")
show()这里简单的为每个三角形创建了alpha图像，然后将所有图像（三角形）合并起来，三角形alpha图像能够简单的通过检查像素的坐标是否能够写成三角形顶点坐标的凸组合来计算得出。通过alpha_for_triangle函数来完成上述工作，代码如下：
def alpha_for_triangle(points, m, n):
"""对于带有由points定义角点的三角形，创建大小为（m,n）的alpha图
（在归一化的齐次坐标意义下）"""
alpha = np.zeros((m, n))
for i in range(min(points[0]), max(points[0])):
for j in range(min(points[1]), max(points[1])):
x = np.linalg.solve(points, [i, j, 1])
if min(x) > 0:
alpha[i, j] = 1
return alpha运行结果：
从图像上可以看到，相同的坐标，使用三角仿射后，使得图像能够贴满公告牌。我们可以看到图像从左上角到右下角被分割成了两块三角形，由于坐标设置的不好，使得两个三角形拼接后发生了位置偏差的问题。
