图像到图像的映射
原理
仿射变换
仿射变换是一种二维坐标到二维坐标之间的线性变换（相同平面），它保持了二维图形的“平直性”（直线经过变换之后依然是直线）和“平行性”（二维图形之间的相对位置关系保持不变，平行线依然是平行线，且直线上点的位置顺序不变），但是角度会改变。在有限维的情况，每个仿射变换可以由一个矩阵A和一个向量b给出，它可以写作A和一个附加的列b。一个仿射变换对应于一个矩阵和一个向量的乘法，而仿射变换的复合对应于普通的矩阵乘法，只要加入一个额外的行到矩阵的底下，这一行全部是0除了最右边是一个1，而列向量的底下要加上一个1。
空间坐标变换经常使用的仿射变换式子如下。其中(x,y)为变换前的坐标，(x’,y’)为变换后的坐标。
由于仿射变换具有6个自由度，因此我们需要三个对应点来计算。
在仿射变换中有两种基本的方法：前向映射和逆向映射，如下图：
以逆向映射为例，对于 g(x’)中的每个像素 x’， 根据变换模型计算相应的映射坐标 x = h-1(x’)， 并将x的像素值赋给g(x’)。前向映射则相反。从原始图像映射到变换图像，赋值的时候需要进行插值运算。通常情况下有三种插值运算：最邻近插值法、双线性插值法、双三次插值法。
单应矩阵
把一张图映射到另一张图上，与前面提到的仿射变换需要变换矩阵T相似，这里我们需要单应矩阵H。我们假设两张图像中的对应点对齐次坐标为(x’,y’,1)和(x,y,1)，单应矩阵H定义为：
注意其中h22=1
则有：
将其展开，我们可以得到：
将其乘以分母得到：
左式展开减去右式，并写成矩阵形式有：
假设我们得到两幅图片中的N个对应点，则有：
如果已经有幅图像，只需要找到原始图像中的任意四个点坐标(其中至少三个点不在同一条直线上)，并且指定目标图像中的四个点，这样通过这八个点，就能求出变换矩阵H。
alpha通道
阿尔法通道是一个8位的灰度通道，该通道用256级灰度来记录图像中的透明度信息，定义透明、不透明和半透明区域，其中白表示不透明，黑表示透明，灰表示半透明。
将第一张图映射在第二张图的部分设为透明，第二张图的其余部分为不透明，即可完成映射。
实现
from PIL import Image
from pylab import *
from scipy import ndimage
def image_in_image(im1,im2,tp):
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
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
def Haffine_from_points(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_cond = dot(C1,fp)
m = 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 = dot(C2,tp)
A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
U,S,V = linalg.svd(A.T)
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]))
# decondition
H = dot(linalg.inv(C2),dot(H,C1))
return H / H[2,2]
def alpha_for_triangle(points,m,n):
alpha = zeros((m,n))
for i in range(min(points[0]),max(points[0])):
for j in range(min(points[1]),max(points[1])):
x = linalg.solve(points,[i,j,1])
if min(x) > 0: #all coefficients positive
alpha[i,j] = 1
return alpha
im1 = array(Image.open('D:\\beatles.jpg').convert('L'))
im2 = array(Image.open('D:\\timg.jpg').convert('L'))
tp = array([[90,259,256,86],[500,498,134,138],[1,1,1,1]])
im3 = image_in_image(im1,im2,tp)
figure()
gray()
subplot(141)
axis('off')
imshow(im1)
subplot(142)
axis('off')
imshow(im2)
subplot(143)
axis('off')
imshow(im3)
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
tp2 = tp[:,:3]
fp2 = fp[:,:3]
H = 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_for_triangle(tp2,im2.shape[0],im2.shape[1])
im3 = (1-alpha)*im2 + alpha*im1_t
tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]
H = 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_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im3 + alpha*im1_t
subplot(144)
imshow(im4)
axis('off')
show()
运行结果如下