创建全景图

创建全景图（实验四）

RANSAC

随机抽样一致算法（RANdom SAmple Consensus,RANSAC）,采用迭代的方式从一组包含离群的被观测数据中估算出数学模型的参数。RANSAC算法假设数据中包含正确数据和异常数据（或称为噪声）。正确数据记为内点（inliers），异常数据记为外点（outliers）。同时RANSAC也假设，给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。该算法基本思想是，数据中包含正确的点，和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。核心思想是随机性和假设性。其中随机性是根据正确数据出现概率去随机选取抽样数据，根据大数定律，随机性模拟可以近似得到正确结果。假设性是假设选取出的抽样数据都是正确数据，然后用这些正确数据通过问题满足的模型，去计算其他点，然后对这次结果进行一个评分。
  示例
   一个简单的例子是从一组观测数据中找出合适的2维直线。假设观测数据中包含局内点和局外点，其中局内点近似的被直线所通过，而局外点远离于直线。简单的最小二乘法不能找到适应于局内点的直线，原因是最小二乘法尽量去适应包括局外点在内的所有点。相反，RANSAC能得出一个仅仅用局内点计算出模型，并且概率还足够高。
  

稳键的单应性矩阵

我我们使用 RANSAC 算法来求解单应性矩阵，首先需要将下面模型类添加到 homography.py 文件中：

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)

     # 归一化齐次坐标
     # for i in range(3):
     #     fp_transformed[i] /= fp_transformed[2]
     fp_transformed = normalize(fp_transformed)

     # 返回每个点的误差
     return sqrt(sum((tp - fp_transformed) ** 2, axis=0))

可以看到，这个类包含 fit() 方法。该方法仅仅接受由 ransac.py 选择的４个对应点对（data 中的前 4 个点对），然后拟合一个单应性矩阵。其中 RANSAC 算法能够判定哪些点对是正确的，哪些是错误的。在实际中，我们需要在距离上使用一个阈值来决定哪些单应性矩阵是合理的。为了方便使用，将下面的函数添加到 homography.py 文件中：

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']

拼接图象

估计出图像间的单应性矩阵，现在我们需要将所有的图像扭曲到一个公共的图像平面上。这里一个方法是创建一个很大的图像，在里面全部填充0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充0，以便为扭曲的图像腾出空间。

def panorama(H,fromim,toim,padding=2400,delta=2400):
  """ Create horizontal panorama by blending two images 
      using a homography H (preferably estimated using RANSAC).
      The result is an image with the same height as toim. 'padding' 
      specifies number of fill pixels and 'delta' additional translation. """ 
  
  # check if images are grayscale or color
  is_color = len(fromim.shape) == 3
  
  # homography transformation for 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 is to the right
      print('warp - right')
      # transform fromim
      if is_color:
          # pad the destination image with zeros to the right
          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:
          # pad the destination image with zeros to the right
          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')
      # add translation to compensate for padding to the left
      H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
      H = dot(H,H_delta)
      # transform fromim
      if is_color:
          # pad the destination image with zeros to the left
          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:
          # pad the destination image with zeros to the left
          toim_t = hstack((zeros((toim.shape[0],padding)),toim))
          fromim_t = ndimage.geometric_transform(fromim,
                                  transf,(toim.shape[0],toim.shape[1]+padding))
  
  # blend and return (put fromim above toim)
  if is_color:
      # all non black pixels
      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

示例图

分析：在学校取景拍了几张图，最后一张图为前面几张图拼接后的图片，这里蓝笔勾画出的为拼接缝（虽然有点不明显）。这里的拼接效果较好，看不出问题。下面就一个反例。

分析：很明显这里是一个拼接不成功的图，通过分析因为我拍摄的角度的垂直角度有明显变化，接着或许是拍摄的时间段为夜里，这里对应特征点取得其实有很多错误的地方，导致了拼接有明显错误。然后可能是仰视拍摄的原因，效果不太好。如果拍摄的角度是平视同时是白天取得景的话，效果可能会好很多。