第三章 图像到图像的映射

本章讲解图像之间的变换，以及一些计算变换的实用方法。这些变换可以用于图像扭曲变形和图像配准。

3.1　单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。本质上，按照下面的方程映射二维中的点（齐次坐标意义下）:$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}\right]=\left[\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}\right]\left[\begin{array}{c}x\\ y\\ w\end{array}\right]\text{或}{\mathbf{x}}^{\prime }=H\mathbf{x}$$

下面的函数可以实现对点进行归一化和转换齐次坐标的功能

import numpy as np


def normalize(points):
    """ 在齐次坐标意义下，对点集进行归一化，使最后一行为 1 """
    for row in points:
        row /= points[-1]
    return points


def make_homog(points):
    """ 将点集（dim×n 的数组）转换为齐次坐标表示 """
    return np.vstack((points, np.ones((1, points.shape[1]))))

复习一下数字图像处理中的放射变换

$$\begin{gathered} 仿射变换 \\ \begin{array}{ccc}变换名称& 仿射矩阵& 坐标公式\\ 恒等变换& \left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]& \begin{array}{c}x=v\\ y=w\end{array}\\ 尺度变换& \left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& 1\end{array}\right]& \begin{array}{c}x=c_{x}v\\ y=c_{y}w\end{array}\\ 旋转变换& \left[\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]& \begin{array}{c}x=v\cos \theta -w\sin \theta \\ y=v\sin \theta +w\cos \theta \end{array}\\ 平移变换& \left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ t_x& t_y& 1\end{array}\right]& \begin{array}{c}x=v+t_x\\ y=w+t_y\end{array}\\ 垂直偏移变换& \left[\begin{array}{ccc}1& 0& 0\\ s_v& 1& 0\\ 0& 0& 1\end{array}\right]& \begin{array}{c}x=v+s_{v}w\\ y=w\end{array}\\ 水平偏移变换& \left[\begin{array}{ccc}1& s_h& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]& \begin{array}{c}x=v\\ y=s_{h}v+w\end{array}\end{array} \end{gathered}$$

import math
import cv2
import numpy as np
import matplotlib.pyplot as plt


def plt_show(*args):
    for ttt in range(len(args)):
        img = args[ttt]
        if (len(img.shape) == 3):
            img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
        elif (len(img.shape) == 2):
            img = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB)
        plt.subplot(231+ttt), plt.imshow(img)

img = cv2.imread('cv.png', 0)

# cv2.imwrite('cv.png',img)
n, m = img.shape
n += 100
m += 100
points = []
for i in range(len(img)):
    for j in range(len(img[0])):
        if img[i][j]:
            points.append([i, j])
angle = 3.14/6
mats = [np.array([1, 0, 0, 0, 1, 0, 0, 0, 1]).reshape(3, 3),  # 恒等变换
        np.array([1.2, 0, 0, 0, 1.2, 0, 0, 0, 1]).reshape(3, 3),  # 尺度变换
        np.array([math.cos(angle), math.sin(angle), 0,  # 旋转变换
                 -math.sin(angle), math.cos(angle), 0,
                 0, 0, 1]).reshape(3, 3),
        np.array([1, 0, 0, 0, 1, 0, 30, 50, 1]).reshape(3, 3),  # 平移变换
        np.array([1, 0, 0, 1.1, 1, 0, 0, 0, 1]).reshape(3, 3),  # 垂直变换
        np.array([1, 1.1, 0, 0, 1, 0, 0, 0, 1]).reshape(3, 3),  # 水平变换
        ]
imgs = []

for mat in mats:
    img_t = np.zeros([n, m])
    for v, w in points:
        x, y, z = np.dot([v, w, 1], mat)
        try:  # 避免新坐标越界
            img_t[int(x)][int(y)] = 254
        except:
            pass
    imgs.append(img_t.astype(np.uint8))
plt_show(imgs[0], imgs[1], imgs[2], imgs[3], imgs[4], imgs[5])

上图分别是：原图、放大、旋转、平移、垂直偏移、水平偏移。这些是基本的仿射变换。

3.1.1　直接线性变换算法

单应性矩阵可以由两幅图像（或者平面）中对应点对计算出来。一个完全射影变换具有 8 个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于 x 和 y 坐标。因此，计算单应性矩阵 H 需要４个对应点对。DLT（Direct Linear Transformation，直接线性变换）是给定４个或者更多对应点对矩阵，来计算单应性矩阵 H 的算法。我们可以得到下面的方程：$$\left[\begin{array}{ccccccccc}-x_1& -y_1& -1& 0& 0& 0& x_1{x}_1^{\prime }& y_1{x}_1^{\prime }& x_1^{\prime }\\ 0& 0& 0& -x_1& -y_1& -1& x_1{y}_1^{\prime }& y_1{y}_1^{\prime }& y_1^{\prime }\\ -x_2& -y_2& -1& 0& 0& 0& x_2{x}_2^{\prime }& y_2{x}_2^{\prime }& x_2^{\prime }\\ 0& 0& 0& -x_2& -y_2& -1& x_2{y}_2^{\prime }& y_2{y}_2^{\prime }& y_2^{\prime }\\ & ⋮& & ⋮& ⋮& ⋮& ⋮& ⋮& \end{array}\right]\left[\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_5\\ h_6\\ h_7\\ h_8\\ h_9\end{array}\right]=0$$ 即 Ah=0。

我们可以使用 SVD（Singular Value Decomposition，奇异值分解）算法找到 H 的最小二乘解。下面是该算法的代码

import numpy as np


def H_from_points(fp, tp):
    """ 使用线性 DLT 方法，计算单应性矩阵 H，使 fp 映射到 tp。点自动进行归一化 """
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
    # 对点进行归一化（对数值计算很重要）
    # --- 映射起始点 ---
    m = np.mean(fp[:2], axis=1)
    maxstd = np.max(np.std(fp[:2], axis=1)) + 1e-9
    C1 = np.diag([1/maxstd, 1/maxstd, 1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp = np.dot(C1, fp)
    # --- 映射对应点 ---
    m = np.mean(tp[:2], axis=1)
    maxstd = np.max(np.std(tp[:2], axis=1)) + 1e-9
    C2 = np.diag([1/maxstd, 1/maxstd, 1])
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp = np.dot(C2, tp)
    # 创建用于线性方法的矩阵，对于每个对应对，在矩阵中会出现两行数值
    nbr_correspondences = fp.shape[1]
    A = np.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 = np.linalg.svd(A)
    H = V[8].reshape((3, 3))
    # 反归一化
    H = np.dot(np.linalg.inv(C2), np.dot(H, C1))
    # 归一化，然后返回
    return H / H[2, 2]

因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题，所以归一化操作非常重要。接下来我们使用对应点对来构造矩阵 A。最小二乘解即为矩阵 SVD 分解后所得矩阵 V 的最后一行。该行经过变形后得到矩阵 H。然后对这个矩阵进行处理和归一化，返回输出

3.1.2　仿射变换

由于仿射变换具有 6 个自由度，因此我们需要三个对应点对来估计矩阵 H。通过将最后两个元素设置为 0，即 h7=h8=0，仿射变换可以用上面的 DLT 算法估计得出。这里我们将使用不同的方法来计算单应性矩阵 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 = np.mean(fp[:2], axis=1)
    maxstd = np.max(np.std(fp[:2], axis=1)) + 1e-9
    C1 = np.diag([1/maxstd, 1/maxstd, 1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp_cond = np.dot(C1, fp)
    # --- 映射对应点 ---
    m = np.mean(tp[:2], axis=1)
    C2 = C1.copy()  # 两个点集，必须都进行相同的缩放
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp_cond = np.dot(C2, tp)
    # 因为归一化后点的均值为 0，所以平移量为 0
    A = np.concatenate((fp_cond[:2], tp_cond[:2]), axis=0)
    U, S, V = np.linalg.svd(A.T)
    # 如 Hartley 和 Zisserman 著的Multiple View Geometry in Computer, Scond Edition 所示，
    # 创建矩阵 B 和 C
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]
    tmp2 = np.concatenate((np.dot(C, np.linalg.pinv(B)), np.zeros((2, 1))), axis=1)
    H = np.vstack((tmp2, [0, 0, 1]))
    # 反归一化
    H = np.dot(np.linalg.inv(C2), np.dot(H, C1))
    return H / H[2, 2]

3.2　图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常应用在计算机图形学中，而且经常出现在计算机视觉算法中。扭曲操作可以使用
SciPy 工具包中的 ndimage 包来简单完成。命令：

transformed_im = ndimage.affine_transform(im,A,b,size)

from scipy import ndimage
from PIL import Image  

img = cv2.cvtColor(cv2.imread('am.png'), cv2.COLOR_BGR2RGB)
H = np.array([[1.4,0.2,-100],[0.1,1.5,-100],[0,0,1]])

im2 = ndimage.affine_transform( cv2.cvtColor(img, cv2.COLOR_BGR2GRAY),  H[:2,:2],(H[0,2],H[1,2]))

plt.gray()
plt.figure(figsize=(8,4))

plt.subplot(1,2,1)
plt.imshow(img)
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(im2)
plt.axis('off')

plt.show()

<Figure size 640x480 with 0 Axes>

3.2.1　图像中的图像

仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。
将函数 image_in_image() 。该函数的输入参数为两幅图像和一个坐标。该坐标为将第一幅图像放置到第二幅图像中的角点坐标：

该图像定义了每个像素从各个图像中获取的像素值成分多少。这里我们基于以下事实，扭曲的图像是在扭曲区域边界之外以 0 来填充的图像，来创建一个二值的 alpha 图像。

def image_in_image(im1, im2, tp):
    """ 使用仿射变换将 im1 放置在 im2 上，使 im1 图像的角和 tp 尽可能的靠近
    tp 是齐次表示的，并且是按照从左上角逆时针计算的 """
    # 扭曲的点
    m, n = im1.shape[:2]
    fp = np.array([[0, m, m, 0], [0, 0, n, n], [0, 0, 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 到 im2 的例子'
im1 = cv2.cvtColor(cv2.imread('am.png'), cv2.COLOR_BGR2GRAY)
im2 = cv2.cvtColor(cv2.imread('letter.png'), cv2.COLOR_BGR2GRAY)

# 选定一些目标点
tp = np.array([[80,600,300,600],[40,40,600,600],[1,1,1,1]])

im3 = image_in_image(im1, im2, tp)
plt.gray()
plt.figure(figsize=(10,5))


plt.subplot(131)
plt.imshow(im1)
plt.axis('off')

plt.subplot(132)
plt.imshow(im2)
plt.axis('off')

plt.subplot(133)
plt.imshow(im3)
plt.axis('off')

plt.show()

<Figure size 640x480 with 0 Axes>

这效果很差，代码有些调不过来，改成用cv2实现。

def update(img,a,b,t):
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if abs(img[i][j][0]-a[0])<=t and abs(img[i][j][1]-a[1])<=t and abs(img[i][j][2]-a[2])<=t :
                img[i][j][0],img[i][j][1],img[i][j][2]=b[0],b[1],b[2]
    return img

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取im1和im2
im1 = cv2.cvtColor(cv2.imread('am.png'), cv2.COLOR_BGR2RGB)
im1=update(im1,[255,255,255],[177,165,123],10)
im2 = cv2.cvtColor(cv2.imread('letter.png'), cv2.COLOR_BGR2RGB)
# 获取im1和im2的尺寸
height1, width1,_ = im1.shape
height2, width2,_ = im2.shape

# 指定im2中的四个点和它们在im1中的对应点
points_im2 = np.array([[219, 160], [1116, 134], [144, 1526], [1214, 1528]], dtype=np.float32)
points_im1 = np.array([[0, 0], [width1, 0], [0, height1], [width1, height1]], dtype=np.float32)

# 计算透视变换矩阵M
M = cv2.getPerspectiveTransform(points_im1, points_im2)

# 使用透视变换将im1拉伸
im1_stretched = cv2.warpPerspective(im1, M, (width2, height2))

# 创建一个新的图像im3，将im1_stretched放入其中
im3 = np.copy(im2)
alpha = (im1_stretched > 0)
im3 = (1 - alpha) * im2 + alpha * im1_stretched

# 显示im1_stretched和im3
plt.gray()
plt.figure(figsize=(20, 8))

plt.subplot(131)
plt.imshow(im1, vmin=0, vmax=255)
plt.axis('off')
plt.title('im1')

plt.subplot(132)
plt.imshow(im2, vmin=0, vmax=255)
plt.axis('off')
plt.title('im2')

plt.subplot(133)
plt.imshow(im3, vmin=0, vmax=255)
plt.axis('off')
plt.title('im3')

plt.show()

<Figure size 640x480 with 0 Axes>

这就完美地将im1放入im2，效果很好

3.2.2　分段仿射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配。让我们看一下对应点对集合之间最常用的扭曲方式：分段仿射扭曲。给定任意图像的标记
点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这些都是最基本的操作。为了三角化这些点，我们经常使用狄洛克三角剖分方法。在 Matplotlib中有狄洛克三角剖分，我们可以用下面的方式使用它：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri

# 生成服从标准正态分布的随机数据
x, y = np.random.normal(0, 1, (2, 100))

# 创建绘图
plt.figure(figsize=(15,7))
plt.subplot(121)
plt.plot(x,y, '*')
plt.axis('off')

plt.subplot(122)
triang = tri.Triangulation(x, y)
plt.triplot(triang, linestyle='-', color='r') 
plt.plot(x, y, '*')
plt.axis('off')

plt.show()

上图是狄洛克三角剖分示例。

3.2.3　图像配准

没有找到 jkface.xml 文件

3.3　创建全景图

在同一位置(即图像的照相机位置相同)拍摄的两幅或者多幅图像是单应性相关的 (如图 3-9 所示)。 我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建全景图像

3.3.1　RANSAC

RANSAC 是“RANdom SAmple Consensus”(随机一致性采样)的缩写。该方法是 用来找到正确模型来拟合带有噪声数据的迭代方法。 给定一个模型，例如点集之间的单应性矩阵，RANSAC 基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

3.3.2　稳健的单应性矩阵估计

我们在任何模型中都可以使用 RANSAC 模块。在使用 RANSAC 模块时， 我们只需 要在相应 Python 类中实现 fit() 和 get_error() 方法，剩下就是正确地使用 ransac.py。 我们这里使用可能的对应点集来自动找到用于全景图像的单应性矩阵。

现显示图片

imname = ['Univ'+str(i+1)+'.jpg' for i in range(5)]
plt.figure(figsize=(15,8))
for i in range(5):
    img=cv2.cvtColor(cv2.imread(imname[i]), cv2.COLOR_BGR2RGB)
    plt.subplot(2,3,i+1)
    plt.axis('off')
    plt.imshow(img)

之后使用 SIFT 特征自动找到匹配对应

import cv2

sift = cv2.SIFT_create()

l = {}
d = {}

for i in range(5):
    image = cv2.imread(imname[i], cv2.IMREAD_GRAYSCALE)
    keypoints = sift.detect(image, None)
    keypoints, descriptors = sift.compute(image, keypoints)
    l[i], d[i] = keypoints, descriptors

matches = {}
for i in range(4):
    bf = cv2.BFMatcher()
    matches[i] = bf.knnMatch(d[i+1], d[i], k=2)

'''homograph.py'''
from numpy import *
from scipy import ndimage
    

class RansacModel(object):
    """ Class for testing homography fit with ransac.py from
        http://www.scipy.org/Cookbook/RANSAC"""
    
    def __init__(self,debug=False):
        self.debug = debug
        
    def fit(self, data):
        """ Fit homography to four selected correspondences. """
        
        # transpose to fit H_from_points()
        data = data.T
        
        # from points
        fp = data[:3,:4]
        # target points
        tp = data[3:,:4]
        
        # fit homography and return
        return H_from_points(fp,tp)
    
    def get_error( self, data, H):
        """ Apply homography to all correspondences, 
            return error for each transformed point. """
        
        data = data.T
        
        # from points
        fp = data[:3]
        # target points
        tp = data[3:]
        
        # transform fp
        fp_transformed = dot(H,fp)
        
        # normalize hom. coordinates
        fp_transformed = normalize(fp_transformed)
                
        # return error per point
        return sqrt( sum((tp-fp_transformed)**2,axis=0) )
        

def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
    """ Robust estimation of homography H from point 
        correspondences using RANSAC (ransac.py from
        http://www.scipy.org/Cookbook/RANSAC).
        
        input: fp,tp (3*n arrays) points in hom. coordinates. """
    
    from PCV.tools import ransac
    
    # group corresponding points
    data = vstack((fp,tp))
    
    # compute H and return
    H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)
    return H,ransac_data['inliers']


def H_from_points(fp,tp):
    """ Find homography H, such that fp is mapped to tp
        using the linear DLT method. Points are conditioned
        automatically. """
    
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    # condition points (important for numerical reasons)
    # --from points--
    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)
    
    # --to points--
    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)
    
    # create matrix for linear method, 2 rows for each correspondence pair
    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))    
    
    # decondition
    H = dot(linalg.inv(C2),dot(H,C1))
    
    # normalize and return
    return H / H[2,2]


def Haffine_from_points(fp,tp):
    """ Find H, affine transformation, such that 
        tp is affine transf of fp. """
    
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    # condition points
    # --from points--
    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)
    
    # --to points--
    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)
    
    # conditioned points have mean zero, so translation is zero
    A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = linalg.svd(A.T)
    
    # create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
    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 normalize(points):
    """ Normalize a collection of points in 
        homogeneous coordinates so that last row = 1. """

    for row in points:
        row /= points[-1]
    return points
    
    
def make_homog(points):
    """ Convert a set of points (dim*n array) to 
        homogeneous coordinates. """
        
    return vstack((points,ones((1,points.shape[1])))) 
 

# 将匹配转换成齐次坐标点的函数
def convert_points(j):
    ndx = [i[0].queryIdx for i in matches[j]]
    fp = np.array([l[j+1][i].pt for i in ndx])
    ndx2 = [i[0].trainIdx for i in matches[j]]
    tp = np.array([l[j][i].pt for i in ndx2])
    return fp, tp

# 估计单应性矩阵
model = RansacModel()

fp,tp = convert_points(0)
H_01 = H_from_ransac(fp,tp,model)[0] # im0 到 im1 的单应性矩阵

fp,tp = convert_points(1)
H_12 = H_from_ransac(fp,tp,model)[0] # im1 到 im2 的单应性矩阵


tp,fp = convert_points(2) # 注意:点是反序的
H_32 = H_from_ransac(fp,tp,model)[0] # im3 到 im2 的单应性矩阵

tp,fp = convert_points(3) # 注意:点是反序的
H_43 = H_from_ransac(fp,tp,model)[0] # im4 到 im3 的单应性矩阵

3.3.3　拼接图像

def panorama(H,fromim,toim,padding=2400,delta=2400):
    """ 使用单应性矩阵 H(使用 RANSAC 健壮性估计得出)，协调两幅图像，创建水平全景图像。结果
        为一幅和 toim 具有相同高度的图像。padding 指定填充像素的数目，delta 指定额外的平移量
    """

    # 检查图像是灰度图像，还是彩色图像

    is_color = len(fromim.shape) == 3
    # 用于 geometric_transform() 的单应性变换

    def transf(p):
        p2 = np.dot(H,[p[0],p[1],1])
        return (p2[0]/p2[2],p2[1]/p2[2])

    if H[1,2]<0: # fromim 在右边
        print( 'warp - right')
        # 变换 fromim
        if is_color:
            # 在目标图像的右边填充 0
            toim_t = np.hstack((toim,np.zeros((toim.shape[0],padding,3))))
            fromim_t = np.zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
            for col in range(3):
                fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
                    transf,(toim.shape[0],toim.shape[1]+padding))

        else:
            # 在目标图像的右边填充 0
            toim_t = hstack((toim,np.zeros((toim.shape[0],padding))))
            fromim_t = ndimage.geometric_transform(fromim,transf,
                        (toim.shape[0],toim.shape[1]+padding))

    else:
        print( 'warp - left')
        # 为了补偿填充效果，在左边加入平移量
        H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
        H = np.dot(H,H_delta)
        # fromim 变换
        if is_color:
            # 在目标图像的左边填充 0
            toim_t = hstack((zeros((toim.shape[0],padding,3)),toim))
            fromim_t = zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
            for col in range(3):
                fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
                    transf,(toim.shape[0],toim.shape[1]+padding))

        else:
            # 在目标图像的左边填充 0
            toim_t = np.hstack((np.zeros((toim.shape[0],padding)),toim))
            fromim_t = ndimage.geometric_transform(fromim,
                transf,(toim.shape[0],toim.shape[1]+padding))

            # 协调后返回(将 fromim 放置在 toim 上)
            if is_color:
                # 所有非黑色像素
                alpha = ((fromim_t[:,:,0] * fromim_t[:,:,1] * fromim_t[:,:,2] ) > 0)
                for col in range(3):
                    toim_t[:,:,col] = fromim_t[:,:,col]*alpha + toim_t[:,:,col]*(1-alpha)

            else:
                alpha = (fromim_t > 0)
                toim_t = fromim_t*alpha + toim_t*(1-alpha)

    return toim_t

# 扭曲图像
delta = 2000 # 用于填充和平移
im1 = np.array(cv2.cvtColor(cv2.imread(imname[1]), cv2.COLOR_BGR2RGB))
im2 = np.array(cv2.cvtColor(cv2.imread(imname[2]), cv2.COLOR_BGR2RGB))
im_12 = panorama(H_12,im1,im2,delta,delta)

im1 = np.array(cv2.cvtColor(cv2.imread(imname[0]), cv2.COLOR_BGR2RGB))
im_02 = panorama(dot(H_12,H_01),im1,im_12,delta,delta)

im1 = np.array(cv2.cvtColor(cv2.imread(imname[3]), cv2.COLOR_BGR2RGB))
im_32 = panorama(H_32,im1,im_02,delta,delta)

im1 = array(cv2.cvtColor(cv2.imread(imname[j+1]), cv2.COLOR_BGR2RGB))
im_42 = warp.panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)

由于图像曝光不同，在单个图像的边界上存在边缘效应，但总体效果还行。