用Harris Corner Detector进行图像拼接

我们这里的图像拼接主要有以下几个步骤:
一、用Harris corner detector检测两幅或多幅图像中的关键点。
二、在这些关键点为中心的领域内构造描述子，并通过描述子找到不同图像之间的 关键点配对
三、用RANSAC筛选鲁棒性高的的关键点配对，并通过关键点配对求得关键点之间的转换矩阵
四、选定一张图像为参考图像，将其他图像通过转换矩阵转换到参考图像的坐标系，并对将每幅图像调整为进行填充处理，以将每幅图像调整为相同大小
五、对填充后的图像进行融合，并对重叠区域进行平滑处理

一、Harris Corner Detector
Harris Corner Detector的出发点为：如果以一个角点为中心构造一个区域，那么这个区域朝任何一个方向移动都会有较大的亮度值变化。如下图所示

如上面左图所示，如果区域在一个亮度没有变化的区域，亮度值当然没有什么改变；如中间图所示，如果区域中只有一条直线，总能找到找到一个方向使得区域内亮度值无显著变化，这个方向沿着这条这条直线方向；如右图所示，如果是角点，在任何方向上亮度都有显著变化。我们假设一个领域朝$(u,v)$方向移动，则移动前后的亮度值变化为：
$$E(u,v)=\sum _{x,y}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$$
这里的$x,y$是某领域内的像素，求和符号为对该领域内所有的像素点进行求和。这里的权重$w$考虑到不同像素点的贡献不同，通常可以取高斯分布或者1。现在如果领域中心为角点，则无论$u,v$取何值，$E$都具有较大的值。这里我们使用多变量函数的泰勒展开:
$$f(x+u,y+v)=f(x,y)+u{f}_x(x,y)+v{f}_y(x,y)+\frac{1}{2!}[u^2{f}_{xx}(x,y)+uv{f}_{xy}(x,y)+v^2{f}_{yy}(x,y)]$$
$$\frac{1}{3!}[u^3{f}_{xxx}(x,y)+u^{2}v{f}_{xxy}(x,y)+u{v}^2{f}_{xyy}(x,y)+v^3{f}_{yyy}(x,y)]+....$$
我们利用上面展开式将$I(x+u,y+v)$展开到一阶,并代入$E(x+u,y+v)$：
$$E(u,v)=\sum _{x,y}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$$
$$=\sum _{x,y}w(x,y)[u{I}_x+v{I}_y+I(x,y)-I(x,y){]}^2=\sum _{x,y}w(x,y)[u{I}_x+v{I}_y{]}^2$$
$$=\sum _{x,y}w(x,y)[u^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2]$$
$$=\sum _{x,y}w(x,y)\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$$
$$=\left[\begin{array}{cc}u& v\end{array}\right](\sum _{x,y}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right])\left[\begin{array}{c}u\\ v\end{array}\right]$$
令$$M=\sum _{x,y}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$$
由于M为实对称矩阵，所以可以写成如下形式：
$$M=P\Sigma P^{-1}=P\Sigma P^T$$
其中$\Sigma =\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]$为特征值矩阵，$P$为特征值对应特征向量矩阵。所以
$$E(u,v)=\left[\begin{array}{cc}u& v\end{array}\right]P\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]P^T{\left[\begin{array}{cc}u& v\end{array}\right]}^T$$
$$=\left[\begin{array}{cc}{u}^{,}& v^{,}\end{array}\right]\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]{\left[\begin{array}{cc}{u}^{，}& v^{,}\end{array}\right]}^T$$
$$=\frac{(u^{，}{)}^2}{\frac{1}{{\lambda }_1}}+\frac{(v^{，}{)}^2}{\frac{1}{{\lambda }_2}}$$
所以当${\lambda }_1$和${\lambda }_2$都很大时，$E(u,v)$也会很大。由于$detM={\lambda }_1{\lambda }_2$,$tranceM={\lambda }_1+{\lambda }_2$,所以可以构造变量$R=detM-k(tranceM{)}^2$。$R$越大，为角点的可能性越大，$k$为可调参数。下面为Harris 角点检测的代码：

import numpy as np
from skimage import filters

def gaussian_kernel(size,sigma):
    gaussian_kernel=np.zeros((size,size))
    for i in range(size):
        for j in range(size):
            x = i - (size-1)/2
            y = j - (size-1)/2
            gaussian_kernel[i,j]=(1/(2*np.pi*sigma**2))*np.exp(-(x**2 + y**2) / (2*sigma**2))
    return gaussian_kernel

def conv(image,kernel):
    m,n = image.shape
    kernel_m,kernel_n = kernel.shape
    image_pad = np.pad(image,((kernel_m//2,kernel_m//2),(kernel_n//2 , kernel_n//2)),'constant')
    result = np.zeros((m,n))
    for i in range(m):
        for j in range(n):
            value = np.sum(image_pad[i:i+kernel_m,j:j+kernel_n]*kernel)
            result[i,j]=value
    return result

def harris_corners(image,window_size=3,k=0.04,window_type=0):
    if window_type==0:
        window=np.ones((window_size,window_size))
    if window_type==1:
        window = gaussian_kernel(window_size,1)
    m,n = image.shape
    dx = filters.sobel_v(image)
    dy = filters.sobel_h(image)
    dx_dx = dx * dx
    dy_dy = dy * dy
    dx_dy = dx * dy
    w_dx_dx = conv(dx_dx,window)
    w_dy_dy = conv(dy_dy,window)
    w_dx_dy = conv(dx_dy,window)
    reponse = np.zeros((m,n))
    for i in range(m):
        for j in range(n):
            M=np.array([[w_dx_dx[i,j],w_dx_dy[i,j]],[w_dx_dy[i,j],w_dy_dy[i,j]]])
            R = np.linalg.det(M)-k*(np.trace(M))**2
            reponse[i,j] = R
    return reponse

第一个函数为构造高斯滤波器的函数。第二个函数为卷积函数。第三个函数为角点检测函数，image为需要检测的图像；window_size为检测角点时以检测点为中心构造的领域尺寸，这里默认为3；参数$k$默认为0.04；由于权重$w(x,y)$可能是高斯函数也可能为1，所以我们这里定义参数window_type，window_type=0时权重取高斯函数，window_type=1时为权重取1。
下面查看检测效果:

from skimage import io
import matplotlib.pyplot as plt
from skimage.feature import corner_peaks
#import utils
import numpy as np
                               
plt.rcParams['figure.figsize'] = (15.0, 12.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

#Harris Corners Detector
image = io.imread('sudoku.png',as_gray=True)

reponse = harris_corners(image,window_size=5,k=0.04,window_type=0)

plt.figure(figsize=(10,10))
plt.subplot(211)
plt.imshow(image)

plt.subplot(212)
plt.imshow(reponse)
plt.show()

corners = corner_peaks(reponse,threshold_rel=0.01)

plt.imshow(image)
plt.scatter(corners[:,1],corners[:,0],marker='x')
plt.axis('off')
plt.title('Detected Corner')
plt.show()

得到:


第一张图片为原图；第二张为用以上Harris角点检测函数得到的响应图像，越亮的地方对应的$R$值越大，越可能是角点；第三张图像是用skimage.feature.corner_peaks函数定义阈值threshold_rel得到的角点检测结果，响应图像中阈值threshold_rel大于0.01，则认为是角点。
这里我们在应用在需要拼接的两张图片上：

#Describing and Matching Keypoints
image1 = io.imread('uttower1.jpg',as_gray=True)
image2 = io.imread('uttower2.jpg',as_gray=True)

keypoint1 = corner_peaks(harris_corners(image1),threshold_rel=0.05,exclude_border=8)
keypoint2 = corner_peaks(harris_corners(image2),threshold_rel=0.05,exclude_border=8)

plt.subplot(211)
plt.imshow(image1)
plt.scatter(keypoint1[:,1],keypoint1[:,0],marker='x')
plt.axis('off')
plt.title('Detected keyponts in image1')

plt.subplot(212)
plt.imshow(image2)
plt.scatter(keypoint2[:,1],keypoint2[:,0],marker='x')
plt.axis('off')
plt.title('Detected keyponts in image2')
plt.show()

得到输出：

二、为关键点构造描述子并匹配
下面我们为两张图片检测到的关键点为中心构建一个区域计算描述子，然后遍历图片1每一个关键点，计算其与图片2每一个关键点的距离，如果最小距离与次小距离的比值小于某个阈值，则认为最小距离对应的图1关键点和图2关键点相对应。这里构建描述子需要描述子函数，我们可以采用区域的HOG描述子，也可以简单地将区域内像素转化为标准正态分布，然后将其展开，作为描述子。
下面为构建描述子函数以及描述子相匹配的函数：

def keypoint_description(image,keypoint,desc_func,patch_size=16):
    keypoint_desc = []
    for i,point in enumerate(keypoint):
        x,y = point
        patch = image[x-patch_size//2:x+int(np.ceil(patch_size/2)),y-patch_size//2:y+int(np.ceil(patch_size/2))]
        description = desc_func(patch)
        keypoint_desc.append(description)
    return np.array(keypoint_desc)

def description_matches(desc1,desc2,threshold=0.5):
    distance_array = cdist(desc1,desc2)
    matches = []
    i=0
    for each_distance_list in distance_array:
        arg_list = np.argsort(each_distance_list)
        index1 = arg_list[0]
        index2 = arg_list[1]
        if each_distance_list[index1] / each_distance_list[index2] <= threshold:
            matches.append([i,index1])
        i+=1
    return np.array(matches)

在第一个函数keypoint_description中参数image为我们想要处理的图片，keypoint为在步骤1中通过Harris角点检测得到的的关键点的位置坐标，desc_func为描述子计算函数，patch_size为以关键点为中心的区域的尺寸。这里的尺寸可以是偶数也可以是奇数，如：对于下面将像素点化为标准正态分布构建的简单描述子，尺寸为偶数；而对于HOG则尺寸一般为偶数。
有了函数keypoint_description，就可以分别计算图1关键点和图2关键点的描述子，计算图1描述子和图2描述子的距离，距离较小的两个关键点我们认为是相对应的。
下面为上面提到的两种描述子，HOG具体见Histogram of Oriented Gradient方向梯度直方图(HOG)：

def simple_descriptor(patch):
    ave = np.mean(patch)
    std = np.std(patch)
    if std==0:
        std=1
    result_patch = (patch - ave) / std
    return result_patch.flatten()

def hog_description(patch,cell_size=(8,8)):
    if patch.shape[0] % cell_size[0]!=0 or patch.shape[1] % cell_size[1]!=0:
        return 'The size of patch and cell don\'t match'
    n_bins=9
    degree_per_bins=20
    Gx = filters.sobel_v(patch)
    Gy = filters.sobel_h(patch)
    G = np.sqrt(Gx**2 + Gy**2)
    theta = (np.arctan2(Gy,Gx) * 180 / np.pi) % 180
    G_as_cells = view_as_blocks(G,block_shape=cell_size)
    theta_as_cells = view_as_blocks(theta,block_shape=cell_size)
    H = G_as_cells.shape[0]
    W = G_as_cells.shape[1]
    bins_accumulator = np.zeros((H,W,n_bins))
    for i in range(H):
        for j in range(W):
            theta_cell = theta_as_cells[i,j,:,:]
            G_cell = G_as_cells[i,j,:,:]
            for p in range(theta_cell.shape[0]):
                for q in range(theta_cell.shape[1]):
                    theta_value = theta_cell[p,q]
                    G_value = G_cell[p,q]
                    num_bins = int(theta_value // degree_per_bins)
                    k= int(theta_value % degree_per_bins)
                    bins_accumulator[i,j,num_bins % n_bins] += (degree_per_bins - k) / degree_per_bins\
                    * G_value
                    bins_accumulator[i,j,(num_bins+1) % n_bins] += k / degree_per_bins * G_value
    Hog_list = []
    for x in range(H-1):
        for y in range(W-1):
            block_description = bins_accumulator[x:x+2,y:y+2]
            block_description = block_description / np.sqrt(np.sum(block_description**2))
            Hog_list.append(block_description)
    return np.array(Hog_list).flatten()

下面我们分别用两种描述子函数去匹配图片1和图片2中的关键点：
1.simple_descriptor

from panorama import simple_descriptor, keypoint_description ,description_matches
from utils import plot_matches

desc1 = keypoint_description(image1,keypoint1,desc_func=simple_descriptor,patch_size=5)
desc2 = keypoint_description(image2,keypoint2,desc_func=simple_descriptor,patch_size=5)
matches = description_matches(desc1,desc2,threshold=0.7)

fig,ax = plt.subplots(1,1,figsize=(15,12))

ax.axis('off')
plot_matches(ax,image1,image2,keypoint1,keypoint2,matches)
plt.show()

输出：

2.hog_description

from panorama import hog_description
desc1 = keypoint_description(image1,keypoint1,hog_description,patch_size=16)
desc2 = keypoint_description(image2,keypoint2,hog_description,patch_size=16)
hog_matches = description_matches(desc1,desc2,threshold=0.7)
fig,ax = plt.subplots(1,1,figsize=(15,12))
ax.axis('off')
plot_matches(ax,image1,image2,keypoint1,keypoint2,hog_matches)
plt.show()

输出：

这里的plot_matches为将两张图相对应的关键点用直线相连接:

def plot_matches(ax,image1,image2,keypoint1,keypoint2,matches):
    H1,W1 = image1.shape
    H2,W2 = image2.shape
    if H1>H2:
        new_image2 = np.zeros((H1,W2))
        new_image2[:H2,:] = image2
        image2 = new_image2
    if H1<H2:
        new_image1 = np.zeros((H2,W1))
        new_image2[:H1,:]=image1
        image1 = new_image1
    image = np.concatenate((image1,image2),axis=1)
    ax.scatter(keypoint1[:,1],keypoint1[:,0],facecolors='none',edgecolors='k')
    ax.scatter(keypoint2[:,1]+image1.shape[1],keypoint2[:,0],facecolors='none',edgecolors='k')
    ax.imshow(image,interpolation='nearest',cmap='gray')
    for one_match in matches:
        index1 = one_match[0]
        index2 = one_match[1]
        color = np.random.rand(3)
        ax.plot((keypoint1[index1,1],keypoint2[index2,1] + image1.shape[1]),
                (keypoint1[index1,0],keypoint2[index2,0]),'-',color=color)

三、用RANSAC筛选鲁棒性高的的关键点配对，并通过关键点配对求得关键点之间的转换矩阵
由于我们从上面两个步骤得到的关键点匹配并不准确，我们从面得到的匹配结果也看出了这一点。为了去除那些不准确的匹配，我们利用RANSAC算法并得到两张图片的转换矩阵。
RANSAC(RANdom SAmple Consensus)算法的主要思想为在拟合一个图形(如圆，直线等)，我们并不使用所有的点，而是随机选择特定数量的数据去拟合(如圆需要3个点确定，直线2个)，得到图形；计算得到的图形与所有点的距离，距离小于某一个阈值的点我们称为内点，否则称为外点。然后利用所有内点再次拟合得到图形。这一过程需要迭代的进行多次，选择内点最多的情况再次拟合。其主要有4步(以拟合直线为例):

而如何利用RANSAC去筛选匹配度较高的匹配？我们已经得到一系列关键点匹配(N个)，我们这里随机选择n个匹配(n<N)，用选择的n个匹配计算图2到图1的转换矩阵；再利用得到的转换矩阵去转换图2的关键点，如果转换得到的图1关键点和实际的关键点差距不大，则认为是内点。这里所有内点对应的匹配实际上为鲁棒性较高的匹配。再利用所有内点计算转换矩阵，这样我们得到了图2到图1的转换矩阵。这个过程同样需要多次迭代进行。
RANSAC算法代码如下:

def fit_affine_matrix(p1,p2):
    assert (p1.shape[0]==p2.shape[0]),'The number of p1 and p2 are different'
    p1=np.hstack((p1,np.ones((p1.shape[0],1))))
    p2=np.hstack((p2,np.ones((p2.shape[0],1))))
    H = np.linalg.pinv(p2) @ p1
    H[:,2]=np.array([0,0,1])
    return H

def ransac(keypoint1,keypoint2,matches,n_iters=200,threshold=20):
    N=matches.shape[0]
    match_keypoints1 = np.hstack((keypoint1[matches[:,0]],np.ones((N,1))))
    match_keypoints2 = np.hstack((keypoint2[matches[:,1]],np.ones((N,1))))
    n_samples=int(N*0.2)
    n_max = 0
    for i in range(n_iters):
        random_index = np.random.choice(N,n_samples,replace=False)
        p1_choice = match_keypoints1[random_index]
        p2_choice = match_keypoints2[random_index]
        H_choice = np.linalg.pinv(p2_choice) @ p1_choice
        H_choice[:,2] = np.array([0,0,1])
        p1_test = match_keypoints2 @ H_choice
        diff = np.sum((match_keypoints1[:,:2]-p1_test[:,:2])**2,axis=1)
        index=np.where(diff<=threshold)[0]
        n_index = index.shape[0]
        if n_index>n_max:
            H=H_choice
            robust_matches=matches[index]
            n_max=n_index
    return H,robust_matches

其中第一个函数为将求关键点矩阵p2到p1的转换矩阵。第二个函数为RANSAC算法，keypoint1和keypoint2分别为图片1和图片2的关键点；matches的维度为N×2，每一行为keypoint1和keypoint2相匹配的关键点的索引；n_iters为随机选择过程的迭代次数；threshold为阈值，如果matches中对应的关键点2通过转换矩阵得到的关键点1和实际关键点1的距离小于此阈值，则认为是内点。输出H为转换矩阵，robust_matches为通过RANSAC筛选得到的鲁棒性较高的匹配。下面我们用RANSAC算法得到筛选后的匹配结果：

from panorama import ransac


H, robust_matches = ransac(keypoint1, keypoint2, hog_matches, threshold=1)

# Plot matches
fig, ax = plt.subplots(1, 1, figsize=(15, 12))
plot_matches(ax, image1, image2, keypoint1, keypoint2, robust_matches)
plt.axis('off')
plt.show()

得到：

四、选定一张图像为参考图像，将其他图像通过转换矩阵转换到参考图像的坐标系，并对将每幅图像调整为进行填充处理，以将每幅图像调整为相同大小

下面我们需要将需要拼接的图片调整为平移并调整相同大小，这样图片就可通过直接叠加图片以得到拼接图片。这里求平移量和调整后图片大小的函数get_output_space和将图片仿射变换函数warp_image：

def get_output_space(image_ref,images,transforms):
    H_ref , W_ref = image_ref.shape
    corner_ref = np.array([[0,0,1],[H_ref,0,1],[0,W_ref,1],[H_ref,W_ref,1]])
    all_corners=[corner_ref]
    if len(images) != len(transforms):
        print('The size of images and transforms does\'t match')
    for i in range(len(images)):
        H,W = images[i].shape
        corner = np.array([[0,0,1],[H,0,1],[0,W,1],[H,W,1]]) @ transforms[i]
        all_corners.append(corner)
    all_corners = np.vstack(all_corners)
    max_corner = np.max(all_corners,axis=0)
    min_corner = np.min(all_corners,axis=0)
    out_space = np.ceil((max_corner - min_corner)[:2]).astype(int)
    offset = min_corner[:2]
    return out_space,offset
        

def warp_image(image, H, output_shape, offset):
     H_invT = np.linalg.inv(H.T)
     matrix = H_invT[:2,:2]
     o = offset+H_invT[:2,2]
     image_warped = affine_transform(image,matrix,o,output_shape,cval=-1)
     return image_warped

这里的get_output_space函数中参数image_ref为参考图片，所有拼接图片都转换到image_ref所对应的坐标系中；images为所有将要转换的图像列表；H为images转换到image_ref的转换矩阵的列表。out_space为所有图片调整到的尺寸大小，offset为所有图片的平移量。这里由于扩大了图像，我们将扩大的图像的空白区域的像素调整为-1。这里并没有调整为0的原因为原图片实际上有可能具有为0的像素。有了out_space和offset，我们就可以利用H去将图片image转换到对应的坐标系中，并调整图片大小并平移。这样图片就可以直接叠加，得到拼接图片。

output_shape, offset = get_output_space(image1, [image2], [H])

image1_warped = warp_image(image1,np.eye(3),output_shape,offset)
image1_mask = (image1_warped != -1)
image1_warped[~image1_mask]=0

image2_warped = warp_image(image2,H,output_shape,offset)
image2_mask = (image2_warped != -1)
image2_warped[~image2_mask]=0

plt.figure(figsize=(15,12))
plt.subplot(121)
plt.imshow(image1_warped)
plt.subplot(122)
plt.imshow(image2_warped)
plt.show()

merged = image1_warped + image2_warped

overlap = np.maximum(image1_mask*1+image2_mask,1)
merged = merged / overlap

plt.figure(figsize=(15,12))
plt.imshow(merged)
plt.show()

输出为：


第一张图片分别是拼接图片装换后调整大小平移后的结果。第二张图片为第一张图片融合后的结果。

五、对填充后的图像进行融合，并对重叠区域进行平滑处理
这里平滑处理主要是先是确定图像的边缘确定重叠区域，再改变重叠区域的像素大小，具体改变方法为：让重叠区域靠近图片1的地方的图片1对应的像素权重较大，图片1对应的像素权重较小，从左到右均匀的增加重叠区域的图2像素点的权重，减小图1像素点的权重，即线性融合，代码如下：

def linear_blend(image1_warped,image2_warped):
    merged = image1_warped + image2_warped
    H , W = image1_warped.shape
    image1_mask = (image1_warped!=0)
    image2_mask = (image2_warped!=0)
    left_margin = np.argmax(image2_mask,axis=1)
    right_margin = W - np.argmax(np.fliplr(image1_mask),axis=1)
    for i in range(H):
        k = right_margin[i] - left_margin[i]
        for j in range(k):
            alpha = j / (k - 1)
            merged[i,left_margin[i]+j] = (1-alpha) * image1_warped[i,left_margin[i]+j]+\
            alpha * image2_warped[i,left_margin[i]+j]
    return merged

应用到经过warp_image函数处理的得到的image1_warped和image2_warped中：

from panorama import linear_blend
merged = linear_blend(image1_warped,image2_warped)
plt.figure(figsize=(15,12))
plt.imshow(merged)
plt.show()

得到：

六、拼接多幅图像
上面我们只拼接了两幅图像，下面我们拼接多幅图像。实际上凭借多幅图像和拼接两幅图像原理相同，只不过在线性融合时。需要先融合图片1和图片2，再用图片1和图片2融合得到的图片去和图片3融合，以此类推。代码如下：

def stitch_multiple_images(images,desc_func=simple_descriptor,patch_size=5):
    keypoints_list = []
    for image in images:
        keypoints= corner_peaks(harris_corners(image))
        keypoints_list.append(keypoints)
    descriptions = []
    for i,keypoints in enumerate(keypoints_list):
        desc = keypoint_description(images[i],keypoints,desc_func,patch_size=patch_size)
        descriptions.append(desc)
    matches_list=[]
    for i in range(len(images)-1):
        matches = description_matches(descriptions[i],descriptions[i+1],threshold=0.7)
        matches_list.append(matches)
    H_list=[]
    for i in range(len(images)-1):
        H,robust_matches = ransac(keypoints_list[i],keypoints_list[i+1],matches_list[i],n_iters=200,threshold=1)
        H_list.append(H)
        matches_list.append(robust_matches)
    n_images = len(images)
    n_ref = n_images//2
    image_ref = images[n_ref]
    a = images.copy()
    images.pop(n_ref)
    images_rest = images.copy()
    images = a.copy()
    H2ref_list=[]
    H_prior = np.eye(3)
    for i in range(n_ref):
        H_next = np.linalg.inv(H_list[n_ref-1-i]) @ H_prior
        H2ref_list.insert(0,H_next)
        H_prior = H_next.copy()
    H_prior = np.eye(3)
    for i in range(n_ref,n_images-1):
         H_next = H_list[i] @ H_prior
         H2ref_list.append(H_next)
         H_prior=H_next.copy()
    output_space,offset = get_output_space(image_ref,images_rest,H2ref_list)
    H2ref_list.insert(n_ref,np.eye(3))
    warps = []
    for i in range(len(images)):
        warp = warp_image(images[i],H2ref_list[i],output_space,offset)
        warp_mask = (warp != -1)
        warp[~warp_mask]=0
        warps.append(warp)
    prior_image = warps[0]
    for i in range(1,len(images)):
        blend_image = linear_blend(prior_image,warps[i])
        prior_image = blend_image
    return blend_image

这里images是所有按顺序排列的待拼接图片的列表。在实现过程中，我们选中间的图片作为参考图片，所有其他图片都转换到其对应的坐标系中。下面我们应用到实际应用中去：

from panorama import stitch_multiple_images
image1 = io.imread('yosemite1.jpg',as_gray=True)
image2 = io.imread('yosemite2.jpg',as_gray=True)
image3 = io.imread('yosemite3.jpg',as_gray=True)
image4 = io.imread('yosemite4.jpg',as_gray=True)
plt.subplot(411)
plt.imshow(image1)
plt.subplot(412)
plt.imshow(image2)
plt.subplot(413)
plt.imshow(image3)
plt.subplot(414)
plt.imshow(image4)
plt.show()
images = [image1,image2,image3,image4]
blend_image= stitch_multiple_images(images,desc_func=simple_descriptor,patch_size=5)
plt.imshow(blend_image)
plt.show()

得到：