计算机视觉之全景拼接

全景拼接,顾名思义就是将同一场景下几张图片拼接在一起,形成一张全貌图。现在的手机相机普遍有这一个功能,今天我们就来学习一下这一有趣的东西,世界如此之大,无奇不有。

一、图像拼接原理

1.提取图像的特征和匹配(SIFT算法)

2.将匹配转化成齐次坐标

3.估计单应性矩阵(RANSAC算法)

4.拼接图像

二、全景拼接实现

1.程序代码

2.homography.make_homog()

3.homography.H_from_ransac()

4.warp.panorama()

三、结果分析

1.全景拼接

2.景深差小的情况

2.景深差大的情况

一、图像拼接原理

首先,我们先了解一下全景拼接的原理。

全景拼接的基础流程如下:

(1)针对同一场景拍摄系列图像

(2)提取图像的特征和匹配

(3)将匹配转化成齐次坐标点

(4)估计单应性矩阵

(5)拼接图像

接下来我们就这几部分分别讲解,当然,第一个步骤拍照,就不赘述了,地球人都会的

1.提取图像的特征和匹配(SIFT算法)

本文运用SIFT算法来提取图像的特征和匹配,是因为SIFT算法对于旋转和尺度均具有不变性,并且对于噪声、视角变化和光照变化具有良好的鲁棒性。关于SIFT算法,博主在之前的博客已经详细讲述了,这里就不重复了,详见下方链接:

https://blog.csdn.net/qq_40369926/article/details/88597406

2.将匹配转化成齐次坐标

将匹配转化成齐次坐标是为了矩阵乘法和点的操作更容易,因为后面计算的单应性矩阵是齐次的。

为什么要引入齐次坐标呢?这是因为齐次坐标可以将图像的各种变化统一成矩阵乘法,详见下方链接:

https://blog.csdn.net/saltriver/article/details/79680364
 

3.估计单应性矩阵(RANSAC算法)

单应性变换是将一个平面内的点映射到另一平面内的二维投影变换,而两张图像变换的对应关系就叫做单应性矩阵,值得注意的是单应性矩阵是齐次的,单应性变换用数学表示如:

其中,h_2_2为1。

因为h_2_2为1,所以单应性矩阵H有8个自由度,求解图像之间的单应性变换其实就是求解单应性矩阵H中的这8个自由度。单应性矩阵可以通过SIFT找到的特征对应点对计算出来。一个完全映射变换需要8个自由度。根据对应点越是,每个对应点可以写出两个方程,分别对应于xy坐标,因此,计算单应性矩阵H需要4个对应点对。

我们RANSAC(RANdom ASmple Consensus,随机一致性采样)来估计单应性矩阵。这其实是DLT(Direct Linear Transformation,直接线性变换)的优化。DLT算法是将给定的4个或更多个对应点方程的系数堆叠到一个矩阵中,使用SVD(奇异值分解)算法求得H的最小二乘解。但是,值得注意的是并非所有的对应点对都是正确的,如下图所示:

这样一口气把所有对应点堆叠在一起求H,很容易导致结果错误,也就是说DLT算法对噪声的鲁棒性并不强。而RANSAC算法经常用来在边矩阵的求解中剔除噪声。RANSAC的基本思想是,数据包含正确的点和噪声点,合理的模型应该能在描述正确数据点的同时摒弃噪声点。RANSAC的求解过程如下:

(1)随机选择4对匹配特征对应点对

(2)根据DLT计算单应性矩阵H(唯一解)

(3)对所有匹配点,计算映射误差(对每个对应点使用该单应性矩阵,然后返回平方距离之和)

(4)根据误差阈值,确定正常值inliers(小于阈值的点看做正确点,否则认为是错误的点)

重复步骤(1)-(4)

(5)针对最大inliers集合,重新计算单应性矩阵H

4.拼接图像

估计出单应性矩阵H后,我们需要将所有图像扭曲到一个公共的图像平面上。通常,这里的公共平面为中心平面(否则,需要进行大量变形)。一种方法是创建一个很大的图像,比如将平面中全部填充为0,使其和中心图像平行,然后将所有的图像扭曲到上面。其基本步骤如下:

(1)实现像素间的映射(计算像素和和单应性矩阵H相乘,然后对齐次坐标进行归一化)

(2)判断图像填补位置(查看H中的平移量,左边为正,右边为负)

(3)在单应性矩阵中加入平移量,进行alpha图填充

二、全景拼接实现

1.程序代码

from pylab import *
from numpy import *
from PIL import Image

from PCV.geometry import homography, warp
from PCV.localdescriptors import sift

"""
这是3.3节中的全景示例。
"""

# 设置数据文件夹的路径

#featname = ['D:/test/Univ'+str(i+1)+'.sift' for i in range(5)] 
#imname = ['D:/test/Univ'+str(i+1)+'.jpg' for i in range(5)]

# 提取特征和匹配
l = {}
d = {}
for i in range(5): 
    sift.process_image(imname[i],featname[i])#处理图像生成sift并保存在文件中
    l[i],d[i] = sift.read_features_from_file(featname[i])#读取特征符并以矩阵形式返回

matches = {}
for i in range(4):
    matches[i] = sift.match(d[i+1],d[i])#匹配特征符

# 可视化匹配
for i in range(4):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i+1]))
    figure()
    sift.plot_matches(im2,im1,l[i+1],l[i],matches[i],show_below=True)


# 将匹配转化成齐次坐标点的函数
def convert_points(j):
    ndx = matches[j].nonzero()[0]
    fp = homography.make_homog(l[j+1][ndx,:2].T) 
    ndx2 = [int(matches[j][i]) for i in ndx]
    tp = homography.make_homog(l[j][ndx2,:2].T) #将点集转化为齐次坐标表示
    
    # 转化x和y  -  TODO这应该移动到其他地方
    fp = vstack([fp[1],fp[0],fp[2]])
    tp = vstack([tp[1],tp[0],tp[2]])
    return fp,tp


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

fp,tp = convert_points(1)
H_12 = homography.H_from_ransac(fp,tp,model)[0] #m1到m2的单应性矩阵

fp,tp = convert_points(0)
H_01 = homography.H_from_ransac(fp,tp,model)[0] #m0到m1的单应性矩阵

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

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


# 扭曲图像
delta = 2000 # 用于填充和平移

im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12,im1,im2,delta,delta)

im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12,H_01),im1,im_12,delta,delta)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32,im1,im_02,delta,delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)


figure()
imshow(array(im_42, "uint8"))
axis('off')
show()

下面介绍程序中调用的关键函数,大家需要特别注意啦!

2.homography.make_homog()

该函数是将点集转化为齐次坐标表示

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

3.homography.H_from_ransac()

该函数用于估计单应性矩阵

class RansacModel(object):
    """ 用于测试单应性矩阵的类,其中单应性矩阵是由网站
        http://www.scipy.org/Cookbook/RANSAC上的ransac.py计算出来的"""
    
    def __init__(self,debug=False):
        self.debug = debug
        
    def fit(self, data):
        """计算选取的4个对应的单应性矩阵 """
        
        #将其转置,来调用H_from_points()计算单应性矩阵
        data = data.T
        
        # 映射的起始点
        fp = data[:3,:4]
        #映射的目标点
        tp = data[3:,:4]
        
        #计算单应性矩阵然后返回
        return H_from_points(fp,tp)
    
    def get_error( self, data, H):
        """对所有的对应计算单应性矩阵,然后对每个变换后的点,返回响应的误差"""
        
        data = data.T
        
        #映射的起始点
        fp = data[:3]
        #映射的目标点
        tp = data[3:]
        
        # 变换fp
        fp_transformed = dot(H,fp)
        
        # 归一化齐次坐标
        fp_transformed = normalize(fp_transformed)
                
        # 返回每个点的误差
        return sqrt( sum((tp-fp_transformed)**2,axis=0) )
        

def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
    """ 使用RANSAC稳健性聚集点对应间的单应性矩阵H(ransac.py为从
        http://www.scipy.org/Cookbook/RANSAC下载的版本).
        
        输入: 齐次坐标表示的点fp,tp (3*n数组) """
    
    from PCV.tools import ransac
    
    # 对应点组
    data = vstack((fp,tp))
    
    # 计算H,并返回
    H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)
    return H,ransac_data['inliers']

4.warp.panorama()

该函数用于图像扭曲。

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 = 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 = hstack((toim,zeros((toim.shape[0],padding,3))))
            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 = hstack((toim,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 = 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 = hstack((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

三、结果分析

1.全景拼接

测试图:

结果图:

从结果图可以看出,我们可以基本实现全景拼接,但是由于图像之间的曝光度不同,我们可以清楚看到右侧存在图像边缘效应。这个是我有待解决的问题。针对这个问题,我们可以采用一些方法将图像强度进行归一化,并对平移进行平滑场景转换。

接下来我们要针对不同景深差进行分析

2.景深差小的情况

测试图:                                                                                         结果图:

       

结果图中左边的图为拼接后的图像,右侧为测试图中左边的图像扭曲后的图像。我们可以清楚的看出,当景深落差小时,图片的扭曲程度较小,对应点可以基本匹配上,不会产生重影或者错位,图像拼接的效果良好。

2.景深差大的情况

测试图:                                                                           结果图:

     

结果图中左边的图为拼接后的图像,右侧为测试图中左边的图像扭曲后的图像。我们可以清楚的看出,当景深落差大时,图片的扭曲程度很大,有时可以甚至完全看不清图像,对应点难以匹配上,会产生严重错位,图像拼接的效果不好。

 

  • 4
    点赞
  • 19
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值