{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第8章 图像拼接"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当我们游览祖国的大好河山时，我们往往一边感慨山川的奇伟峻丽，一边遗憾自己没有设备去拍一张全景图，而只能拍一些不同视角下分离的照片。通过这一章的学习，或许我们就可以自己实现全景图的生成。\n",
    "\n",
    "图像拼接 (Stitching)，是指将多张拍摄同一场景的有重叠的图像（可能是不同时间、不同视角或者不同传感器获得的）拼接成一幅无缝的大尺寸全景图的技术。图像拼接在运动检测和跟踪、增强现实、分辨率增强、视频压缩和图像稳定等计算机视觉任务上有很多的应用。图 8-1 和图 8-2 展示了一个图像拼接的例子。接下来，就让我们一起学习如何实现图像拼接。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-d564a1d560b5bfb7d3a967b260e63a32_1440w.jpg\" width=800>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-1 同一场景在不同视角下分离的图像。</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-b06869d874d36b527144f68673eb91a8_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-2 图像拼接。</div>\n",
    "</center>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.2 图像变换"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大家仔细观察图 8-1 中的一系列图像与图 8-2 中的全景图，他们在几何上有什么相似性吗？没错，图 8-2 中的全景图可以通过图 8-1 中的子图经过缩放旋转投影等一些操作之后再拼接组合得到，对图像进行的这些操作我们统称为“图像变换”（Image Warping）。接下来，我们先简单介绍几种常用的图像变换。\n",
    "\n",
    "1. 平移：将图像$\\mathbf{I}$沿着$[t_x, t_y]$平移，设图像上一点$\\mathbf{p}=[i,j]$经过平移后得到新的点为$\\mathbf{p}^{\\prime}=[i^{\\prime}, j^{\\prime}]$，我们有$\\mathbf{p}^{\\prime} = \\mathbf{p} + \\mathbf{t}$，其中，$\\mathbf{t} = [t_x, t_y]$为平移方向；\n",
    "\n",
    "2. 旋转：我们将图像沿着图像左上角起始点$[0,0]$逆时针旋转$\\theta$角度，则图像上一点$\\mathbf{p}=[i,j]$旋转之后对应的新的点坐标为$\\mathbf{p}^{\\prime} = \\mathbf{R}\\mathbf{p}$，其中，$\\mathbf{R} = \\begin{bmatrix}\n",
    "\\cos\\theta &-\\sin\\theta \\\\\n",
    "\\sin\\theta & \\cos\\theta\n",
    "\\end{bmatrix}$\n",
    "表示旋转矩阵；\n",
    "\n",
    "3. 缩放：以图像左上角起始点$[0,0]$为缩放中心对图像沿$x$轴缩放$s_x$倍，沿$y$轴缩放$s_y$倍，则图像上一点$\\mathbf{p}=[i,j]$缩放之后对应的点坐标为$\\mathbf{p}^{\\prime} = \\mathbf{S}\\mathbf{p}$。这里，$\\mathbf{S} = \\begin{bmatrix}\n",
    "s_x & 0\\\\\n",
    "0 & s_y\n",
    "\\end{bmatrix}$是缩放矩阵；\n",
    "\n",
    "4. 对称：将图像沿着$y$轴镜像对称，我们不难得到图像上一点$\\mathbf{p}=[i,j]$对称之后对应的点坐标为$\\mathbf{p}^{\\prime} = [-i,j]$，那么$\\mathbf{p}^{\\prime} = \\begin{bmatrix}\n",
    "-1 & 0\\\\\n",
    "0 & 1\n",
    "\\end{bmatrix}\\mathbf{p}$；如果将图像沿直线$y=x$对称，则$\\mathbf{p}^{\\prime} = [j,i]$，即$\\mathbf{p}^{\\prime} = \\begin{bmatrix}\n",
    "0 & 1\\\\\n",
    "1 & 0\n",
    "\\end{bmatrix}\\mathbf{p}$。\n",
    "\n",
    "我们不难发现，上述2、3、4中图像变换前后的对应关系是线性的，即可以写成$\\mathbf{p}^{\\prime} = \\mathbf{T}\\mathbf{p}$的矩阵形式，而平移却不可以用这样的线性变换来表示。考虑到矩阵乘法的线性变换形式易于计算，那么我们能不能也将平移用矩阵乘法形式表示呢？\n",
    "\n",
    "当然可以。只不过，我们需要把常用的笛卡尔坐标系拓展至齐次坐标系（Homogenenous Coordinates）。在笛卡尔坐标系下，对于一个向量$\\mathbf{v}$，令$\\lambda$是一个常数，我们知道一般$\\mathbf{v}\\neq\\lambda\\mathbf{v}$；而在齐次坐标系下，总有$\\mathbf{v}=\\lambda\\mathbf{v}$。那么笛卡尔坐标系和齐次坐标系如何转换呢？如图 8-3 所示，我们在$x$轴、$y$轴的基础上增加$w$轴，并将点$\\mathbf{p}=(i,j)$表示为$\\mathbf{p}=(i,j,1)$，我们把这样的坐标系称为齐次坐标系。对于齐次坐标系中任意一点$(x,y,w)$，其对应的笛卡尔坐标系中的坐标为$(x/w, y/w)$。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-90f2a36dceb54cadabf0c197672d63c4_1440w.jpg\" width=300>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-3 齐次坐标系</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "在这基础上，我们便可以把平移前后的点坐标表示为：\n",
    "\n",
    "$$\\begin{bmatrix}i^{\\prime}\\\\j^{\\prime}\\\\1 \\end{bmatrix}= \\begin{bmatrix}\n",
    "1 & 0 & t_x \\\\\n",
    "0 & 1 & t_y \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix} \\begin{bmatrix} i \\\\j \\\\1 \\end{bmatrix}$$\n",
    "\n",
    "同样，对于旋转、缩放，我们也可以将其依次表示为：\n",
    "\n",
    "旋转：$\\begin{bmatrix}i^{\\prime}\\\\j^{\\prime}\\\\1 \\end{bmatrix}= \\begin{bmatrix}\n",
    "\\cos\\theta &-\\sin\\theta & 0 \\\\\n",
    "\\sin\\theta & \\cos\\theta & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix} \\begin{bmatrix} i \\\\j \\\\1 \\end{bmatrix}$，\n",
    "缩放：$\\begin{bmatrix}i^{\\prime}\\\\j^{\\prime}\\\\1 \\end{bmatrix}= \\begin{bmatrix}\n",
    "s_x & 0 & 0 \\\\\n",
    "0 & s_y & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix} \\begin{bmatrix} i \\\\j \\\\1 \\end{bmatrix}$。\n",
    "\n",
    "至此，在齐次坐标系下，上述的图像变换操作都可以用一个$3\\times 3$的矩阵$\\mathbf{T}$来表示。不知道细心的同学观察到没有，平移、旋转、缩放与对称操作对应的矩阵$\\mathbf{T}$的最后一行都是$[0,0,1]$。我们把形如$\\mathbf{T}=\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}$对应的图像变换称为仿射变换。那变换矩阵$\\mathbf{T}$的最后一行不是$[0,0,1]$，而是任意的$[g,h,1]$，图像会发生什么样的变化呢？\n",
    "\n",
    "我们可以通过实验得到类似于图 8-4 的结果，很明显，在$\\mathbf{T}$最后一行的值不再是$[0,0,1]$的情况下，图像发生了一定程度的扭曲与畸变，我们把形如$\\mathbf{T}=\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "g & h & 1\n",
    "\\end{bmatrix}$变换操作称作投影（Homography）。投影的变换矩阵形式是最通用的，所以投影包含了上述所有图像变换。不同视角下的平面目标和固定相机位置下旋转相机角度，这两种情况下拍摄得到的图像都满足投影变换。通常我们对同一场景拍摄多张图像的情况总可以近似为上述两种情况，所以我们可以通过投影变换来实现全景图的拼接。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-8a1f6066c5b8533dc930e93ea06e95f9_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-4 投影变换。</div>\n",
    "</center>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 图像拼接"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.1 计算变换矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像拼接是要将一组图像拼接成一张大的全景图。这个任务中最基本的问题就是如何拼接两张图像。给定两张图像，要将它们拼接在一起，就需要计算它们之间的变换矩阵。要计算变换矩阵$\\mathbf{T}$，就需要找到两张图像之间的匹配关系。还记得我们在上一章中学习过的SIFT特征点检测算法和SIFT特征描述吗？通过SIFT特征点检测算法找到两张图像上的特征点，然后基于SIFT特征描述，就可以确定这些特征点之间的匹配关系，如图 8-5 所示。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pica.zhimg.com/80/v2-12e218c65547fcfe8a00fde248cf800f_1440w.jpg?source=d16d100b\" width=700>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-5 基于特征点的图像匹配。 </div>\n",
    "</center>\n",
    "\n",
    "设$\\mathbf{p}=[x,y]$和$\\mathbf{p}^{\\prime}=[x^{\\prime},y^{\\prime}]$是两张图像中匹配的一对特征点，在齐次坐标系中表示为$\\mathbf{p} = [x,y,1]^\\mathsf{T}, \\mathbf{p}^{\\prime} = [x^{\\prime},y^{\\prime},1]^\\mathsf{T}$，则我们有$\\mathbf{p}^{\\prime}=\\mathbf{T}\\mathbf{p}$，即\n",
    "\n",
    "$$\\begin{bmatrix}x^{\\prime} \\\\ y^{\\prime}\\\\1 \\end{bmatrix}= \\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "g & h & i\n",
    "\\end{bmatrix} \\begin{bmatrix} x \\\\ y \\\\ 1 \\end{bmatrix}$$\n",
    "\n",
    "我们可以得到$x,y$和$x^{\\prime},y^{\\prime}$之间的关系为：\n",
    "\n",
    "$$\n",
    "x^{\\prime} = \\frac{ax+by+c}{gx+hy+i},\n",
    "y^{\\prime} = \\frac{dx+ey+f}{gx+hy+i}\n",
    "$$\n",
    "\n",
    "为了更方便的求解矩阵$\\mathbf{T}$中的未知量，我们将上式写成\n",
    "\n",
    "$$\n",
    "x^{\\prime}(gx+hy+i) = ax+by+c$$$$\n",
    "y^{\\prime}(gx+hy+i) = dx+ey+f\n",
    "$$\n",
    "\n",
    "即\n",
    "\n",
    "$$\\begin{bmatrix}x & y & 1 & 0 & 0 & 0 &-x^{\\prime}x& -x^{\\prime}y & -x^{\\prime} \\\\ 0 & 0 & 0 & x & y & 1 & -y^{\\prime}x & -y^{\\prime}y & -y^{\\prime}\\end{bmatrix}  \\begin{bmatrix} a \\\\ b \\\\ c \\\\ d \\\\ e \\\\ f \\\\ g \\\\ h \\\\ i \\end{bmatrix} = \\begin{bmatrix} 0 \\\\ 0 \\end{bmatrix}$$\n",
    "\n",
    "假设两张图像共有$n$对匹配的特征点，我们可以得到以下方程：\n",
    "\n",
    "$$\\begin{bmatrix}x_1 & y_1 & 1 & 0 & 0 & 0 &-x_1^{\\prime}x_1& -x_1^{\\prime}y_1 & -x_1^{\\prime} \\\\ 0 & 0 & 0 & x_1 & y_1 & 1 & -y_1^{\\prime}x_1 & -y_1^{\\prime}y_1 & -y_1^{\\prime} \\\\ ... \\\\x_n & y_n & 1 & 0 & 0 & 0 &-x_n^{\\prime}x_n& -x_n^{\\prime}y_n & -x_n^{\\prime} \\\\ 0 & 0 & 0 & x_n & y_n & 1 & -y_n^{\\prime}x_n & -y_n^{\\prime}y_n & -y_n^{\\prime}\\end{bmatrix}  \\begin{bmatrix} a \\\\ b \\\\ c \\\\ d \\\\ e \\\\ f \\\\ g \\\\ h \\\\ i \\end{bmatrix} = \\begin{bmatrix} 0 \\\\ 0 \\\\...\\\\ 0\\end{bmatrix}_{2n}$$\n",
    "\n",
    "令矩阵$\\mathbf{A} = \\begin{bmatrix}x_1 & y_1 & 1 & 0 & 0 & 0 &-x_1^{\\prime}x_1& -x_1^{\\prime}y_1 & -x_1^{\\prime} \\\\ 0 & 0 & 0 & x_1 & y_1 & 1 & -y_1^{\\prime}x_1 & -y_1^{\\prime}y_1 & -y_1^{\\prime} \\\\ ... \\\\x_n & y_n & 1 & 0 & 0 & 0 &-x_n^{\\prime}x_n& -x_n^{\\prime}y_n & -x_n^{\\prime} \\\\ 0 & 0 & 0 & x_n & y_n & 1 & -y_n^{\\prime}x_n & -y_n^{\\prime}y_n & -y_n^{\\prime}\\end{bmatrix}_{2n\\times 9}$, $\\mathbf{t} = [a,b,c,d,e,f,g,h,i]^\\mathsf{T}$，上述等式可以简写成\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{t}=\\mathbf{0}$$\n",
    "\n",
    "一般来说，我们能够得到的特征点对的数量远多于方程未知量的个数，且由于有的特征点对的位置存在误差，上式可能没有唯一闭式解。因此，我们往往通过优化下述问题来求解$\\mathbf{t}$:\n",
    "\n",
    "$$\\min \\|\\mathbf{At}-\\mathbf{0}\\|_F^2, {\\rm s.t.} \\|\\mathbf{t}\\|^2=1$$\n",
    "\n",
    "因此，我们只需要求解上述方程得到$\\mathbf{t}$，便可得到变换矩阵$\\mathbf{T}$。注意到变换矩阵$\\mathbf{T}$是定义在齐次坐标系下的，所以尺度对其无影响。故而，我们只需要求解单位向量$\\hat{\\mathbf{t}}$，即为矩阵$\\mathbf{A}^\\mathsf{T}\\mathbf{A}$的最小特征值对应的特征向量（证明过程详见附录）。我们把图像拼接的过程总结为以下三步：\n",
    "\n",
    "1. 计算两幅图像中的特征点；\n",
    "2. 将两张图像的特征点进行匹配；\n",
    "3. 根据匹配的特征点对计算图像变换矩阵。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.2 利用RANSAC算法去除误匹配"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当我们利用SIFT进行特征匹配时，有些时候可能会出现图 8-6 的情况。图 8-6 中的左图对应右图中绿色圆圈内，但我们利用SIFT匹配特征点时，会将左图中部分特征点匹配到右图绿色圆圈之外的特征点（如红色区域）。这样一些特征点匹配是错误的匹配，应该被移除，从而保证变换矩阵计算的鲁棒性。那么我们应该如何移除错误的匹配点对呢？\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-cc484a49132c92e7ceeb0be3cc0e4eb4_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-6 基于SIFT特征描述得到的两张图像中的匹配特征点对。</div>\n",
    "</center>\n",
    "\n",
    "我们先放下匹配这个问题，从一个更简单的例子入手。图 8-7a 显示了一组数据点，我们希望用一条直线来拟合这些数据点。但这些数据点的分布并不是一条直线的分布，所以其中肯定有所谓的异常样本，即异常的数据点。然而，我们并不知道哪些数据点是正常的样本，哪些数据点是异常的样本（outlier），因此我们并不知道应该用哪部分点去拟合直线。但是，从直觉出发，我们会认为图 8-7c 中的直线比图 8-7b 中的直线拟合效果更好，因为有更多的数据点分布在绿色直线的附近。深入思考，正常的数据点是大多数的具有相似分布规律的数据点，而异常数据点的分布则杂乱无章。因此，正常的样本总是有很多一致的同类，称为“inlier”，异常的点和正常的样本不一致，称为“outlier”。如果我们用任意两个数据点拟合一条直线，如果有越多的数据点与这条直线一致，即分布在其附近，则这条直线具有越多的与其一致的inlier，这条直线也就越有可能是我们希望拟合的直线，而与这条直线不符的数据点则是outlier。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-c468a014945893ded1c7f1bb8f313c1e_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-7 直线拟合。（a）用于直线拟合的数据点。（b）与红色直线一致的数据点较少。（c）与绿色直线一致的数据点很多。 </div>\n",
    "</center>\n",
    "\n",
    "随机样本一致算法RANSAC（Random SAmple Consensus）就是利用上述思想从可能含有异常值的观测数据集中拟合数学模型，是一种鲁棒参数估计的迭代方法。RANSAC算法的步骤如下：\n",
    "\n",
    "1. 随机选择$s$个样本（其中，$s$是可以用于拟合模型的最小样本量）；\n",
    "2. 用这些样本拟合模型；\n",
    "3. 统计与模型一致的样本（inlier）个数；\n",
    "4. 上述过程重复N次；\n",
    "5. 选择具有最多一致样本（inlier）个数的模型作为最终拟合的结果。\n",
    "\n",
    "对于直线拟合这问题，一个样本是指一个数据点，所以这里$s=2$，因为两个数据点确定一条直线。回到图像特征点对匹配问题这个问题，这样的一个样本是指一对匹配的特征点对。同为inlier的匹配特征点对应具有一致的几何变换关系，而作为outlier的误匹配特征点对，它们之间则各有各的不同。对于不同的图像变换，$s$大小不同。比如对于平移变换，$s=1$，因为只需要一对匹配点对就可以确定平移变换的参数；而对于投影变换，$s=4$，大家想想为什么？\n",
    "\n",
    "图 8-8 展示了使用RANSAC算法拟合平移变换的一次试验。在左图的试验中，我们首先选择黄色箭头代表的平移拟合模型，那么与该平移模型一致的inlier（蓝色箭头与黄色箭头）个数为5，而与之不一致的outlier（红色箭头）个数为2。在下一次尝试中（右图），我们再随机选择数据（黄色箭头）拟合对应的模型，并统计inlier的个数（只有一个）。最终，我们会选择inlier最多的模型（在图 8-8 中为蓝色箭头的方向）作为最终拟合的结果。而那些匹配错误的特征点对（如红色线）则没有被用到，即这些outlier被“去除”。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-b363dd13843dd12c05eb36f8f979fec1_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-8 利用RANSAC算法剔除匹配异常的特征点对 </div>\n",
    "</center>\n",
    "\n",
    "\n",
    "\n",
    "需要说明的是，RANSAC算法不能保证得到最优解，因为它只在一定的概率下去逼近最优解。不过随着执行更多的迭代，逼近最优解的概率会增加。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.3 图像变换与缝合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像拼接的最后一步是将输入图像变换并缝合到一张图像中。对于两张图像$\\mathbf{A}$和$\\mathbf{B}$，我们在已经检测出对应的特征点对，并利用RANSAC算法计算得到变换矩阵$\\mathbf{T}$之后，我们将图像$\\mathbf{B}$整体变换为$\\mathbf{TB}$。在这之后，我们将转换后的图像$\\mathbf{B}$（即$\\mathbf{TB}$）与图像$\\mathbf{A}$在重叠部分上的像素值求取平均值，以优化图像缝合的边界。如此，我们便可得到最终缝合好的拼接图像。\n",
    "\n",
    "\n",
    "综上所述，我们把图像拼接的全过程总结为以下四步：\n",
    "1. 计算两幅图像中的特征点；\n",
    "2. 将两张图像的特征点进行匹配；\n",
    "3. 根据匹配的特征点对利用RANSAC算法计算图像变换矩阵；\n",
    "4. 将图像进行拼接。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 编程实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一节中，我们将编程实现图像拼接。我们先导入上一章中编写的SIFT相关函数，并在此基础按照8.3中讲述的流程编写实现RANSAC算法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:14:27.878595Z",
     "iopub.status.busy": "2022-09-16T03:14:27.878175Z",
     "iopub.status.idle": "2022-09-16T03:14:27.904993Z",
     "shell.execute_reply": "2022-09-16T03:14:27.903759Z",
     "shell.execute_reply.started": "2022-09-16T03:14:27.878559Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import *\n",
    "\n",
    "# 利用RANSAC算法剔除无效点\n",
    "def compute_affine_xform(corners1, corners2, matches):\n",
    "    \n",
    "    # 设置最大迭代次数\n",
    "    iteration = 50\n",
    "    M_list = []\n",
    "    \n",
    "    # 存放每一次模型的inlier个数\n",
    "    inlier_num_list = []\n",
    "    for _ in range(iteration):\n",
    "        \n",
    "        # 随机选择4个特征点对来拟合模型\n",
    "        sample_index = random.sample(range(len(matches)), 4)\n",
    "        x1_s, y1_s = corners1[matches[sample_index[0]][0]]\n",
    "        x1_t, y1_t = corners2[matches[sample_index[0]][1]]\n",
    "        x2_s, y2_s = corners1[matches[sample_index[1]][0]]\n",
    "        x2_t, y2_t = corners2[matches[sample_index[1]][1]]\n",
    "        x3_s, y3_s = corners1[matches[sample_index[2]][0]]\n",
    "        x3_t, y3_t = corners2[matches[sample_index[2]][1]]\n",
    "        x4_s, y4_s = corners1[matches[sample_index[3]][0]]\n",
    "        x4_t, y4_t = corners2[matches[sample_index[3]][1]]\n",
    "        \n",
    "        # 编写矩阵A\n",
    "        A = np.array([[x1_s, y1_s, 1, 0, 0, 0, -x1_t*x1_s, -x1_t*y1_s, -x1_t],\n",
    "                      [0, 0, 0, x1_s, y1_s, 1, -y1_t*x1_s, -y1_t*y1_s, -y1_t],\n",
    "                      [x2_s, y2_s, 1, 0, 0, 0, -x2_t*x2_s, -x2_t*y2_s, -x2_t],\n",
    "                      [0, 0, 0, x2_s, y2_s, 1, -y2_t*x2_s, -y2_t*y2_s, -y2_t],\n",
    "                      [x3_s, y3_s, 1, 0, 0, 0, -x3_t*x3_s, -x3_t*y3_s, -x3_t],\n",
    "                      [0, 0, 0, x3_s, y3_s, 1, -y3_t*x3_s, -y3_t*y3_s, -y3_t],\n",
    "                      [x4_s, y4_s, 1, 0, 0, 0, -x4_t*x4_s, -x4_t*y4_s, -x4_t],\n",
    "                      [0, 0, 0, x4_s, y4_s, 1, -y4_t*x4_s, -y4_t*y4_s, -y4_t]\n",
    "                      ])\n",
    "        \n",
    "        # 求解A的特征向量\n",
    "        _,_, v = np.linalg.svd(A)\n",
    "        \n",
    "        # 取最小特征值对应的特征向量作为最终的结果\n",
    "        M = np.reshape(v[-1], (3, 3))\n",
    "        \n",
    "        inlier_num = 0\n",
    "        # 统计inlier个数\n",
    "        for (index1, index2) in matches:\n",
    "            # coord是齐次坐标系的坐标\n",
    "            coord1 = [corners1[index1][0], corners1[index1][1], 1]\n",
    "            coord2 = [corners2[index2][0], corners2[index2][1], 1]\n",
    "            # 计算将coord1齐次变换之后的坐标\n",
    "            mapcoor = np.dot(M, coord1)\n",
    "            # 将齐次坐标系中的w置为1\n",
    "            mapcoor = mapcoor / mapcoor[-1]\n",
    "            if np.linalg.norm(coord2 - mapcoor) < 5:\n",
    "                inlier_num += 1\n",
    "                \n",
    "        # 将inlier个数和对应的变换矩阵M记录在列表里\n",
    "        M_list.append(M)\n",
    "        inlier_num_list.append(inlier_num)\n",
    "    \n",
    "    # 获取列表中inlier值最大的元素的下标\n",
    "    best_index = np.argmax(inlier_num_list)\n",
    "    # 获取对应的变换矩阵\n",
    "    xform = M_list[best_index].astype(np.float64)\n",
    "    \n",
    "    # 统计属于模型outlier的特征点对\n",
    "    outlier_labels = []\n",
    "    for (index1, index2) in matches:\n",
    "        coord1 = [corners1[index1][0], corners1[index1][1], 1]\n",
    "        coord2 = [corners2[index2][0], corners2[index2][1], 1]\n",
    "        mapcoor = np.dot(xform, coord1)\n",
    "        mapcoor = mapcoor/mapcoor[-1]\n",
    "        if np.linalg.norm(coord2 - mapcoor) < 12:\n",
    "            outlier_labels.append(1)\n",
    "        else:\n",
    "            outlier_labels.append(0)\n",
    "    \n",
    "    return xform, outlier_labels  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了在绘制时能够体现出异常点对，我们改写第7章中的draw_matches()函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_matches(image1, image2, corners1, corners2, matches, outliers=None):\n",
    "    # 获得两张图像的分辨率\n",
    "    h1, w1 = image1.shape\n",
    "    h2, w2 = image2.shape\n",
    "    hres = 0\n",
    "    if h1 >= h2:\n",
    "        hres = int((h1 - h2) / 2)\n",
    "    \n",
    "        # 将两张图像拼接成高度一致的图像，方便进行特征点对的比对\n",
    "        match_image = np.zeros((h1, w1 + w2, 3), np.uint8)\n",
    "\n",
    "        # 对R、G、B图像分别处理\n",
    "        for i in range(3):\n",
    "            match_image[: h1, : w1, i] = image1\n",
    "            match_image[hres: hres + h2, w1: w1 + w2, i] = image2\n",
    "\n",
    "        for i in range(len(matches)):\n",
    "            m = matches[i]\n",
    "            # 获得匹配的特征点对在图中的坐标\n",
    "            pt1 = (int(corners1[m[0]][0]), int(corners1[m[0]][1]))\n",
    "            pt2 = (int(corners2[m[1]][0] + w1), int(corners2[m[1]][1] + hres))\n",
    "            # 将其圈出\n",
    "            cv2.circle(match_image, pt1, 1, (0,255,0), 2)\n",
    "            cv2.circle(match_image, (pt2[0], pt2[1]), 1, (0,255,0), 2)\n",
    "            # 画线相连\n",
    "            if outliers:\n",
    "                cv2.line(match_image, pt1, pt2, (255, 0, 0))\n",
    "            else:\n",
    "                cv2.line(match_image, pt1, pt2, (0, 0, 255))\n",
    "    else:\n",
    "        hres = int((h2 - h1) / 2)\n",
    "        \n",
    "         # 将两张图像拼接成高度一致的图像，方便进行特征点对的比对\n",
    "        match_image = np.zeros((h2, w1 + w2, 3), np.uint8)\n",
    "\n",
    "        # 对R、G、B图像分别处理\n",
    "        for i in range(3):\n",
    "            match_image[hres: hres + h1, : w1, i] = image1\n",
    "            match_image[: h2, w1: w1 + w2, i] = image2\n",
    "\n",
    "        for i in range(len(matches)):\n",
    "            m = matches[i]\n",
    "            # 获得匹配的特征点对在图中的坐标\n",
    "            pt1 = (int(corners1[m[0]][0]), int(corners1[m[0]][1] + hres))\n",
    "            pt2 = (int(corners2[m[1]][0] + w1), int(corners2[m[1]][1]))\n",
    "            # 将其圈出\n",
    "            cv2.circle(match_image, pt1, 1, (0,255,0), 2)\n",
    "            cv2.circle(match_image, (pt2[0], pt2[1]), 1, (0,255,0), 2)\n",
    "            # 画线相连\n",
    "            if outliers:\n",
    "                cv2.line(match_image, pt1, pt2, (255, 0, 0))\n",
    "            else:\n",
    "                cv2.line(match_image, pt1, pt2, (0, 0, 255))\n",
    "        \n",
    "    return match_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之后，我们编写图像变换与缝合的相关函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 图像变换与缝合\n",
    "def stitch_images(image1, image2, xform):\n",
    "    # 将图像按照变换矩阵的方向进行变换\n",
    "    new_image = cv2.warpPerspective(image1, xform, (image2.shape[1],image2.shape[0]))\n",
    "    image = new_image.copy()\n",
    "    \n",
    "    # 获取变换后的图像位置坐标\n",
    "    h, w, _ = image1.shape\n",
    "    pts = np.float32([[0, 0],\n",
    "                      [0, h - 1],\n",
    "                      [w - 1, h - 1],\n",
    "                      [w - 1, 0]]).reshape(-1, 1, 2)\n",
    "    dst = cv2.perspectiveTransform(pts, xform)\n",
    "    \n",
    "    # 画出缝合线\n",
    "    image = cv2.polylines(image2, [np.int32(dst)], True, 255, 1, cv2.LINE_AA)\n",
    "    \n",
    "    # 将重叠部分的图像求平均\n",
    "    image = np.where(new_image>0, 0.5*new_image+0.5*image, image)\n",
    "    return image.astype(np.uint8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "导入两张图像$\\mathbf{A}$和$\\mathbf{B}$，我们先计算图像的变换矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:14:27.906774Z",
     "iopub.status.busy": "2022-09-16T03:14:27.906345Z",
     "iopub.status.idle": "2022-09-16T03:15:33.710850Z",
     "shell.execute_reply": "2022-09-16T03:15:33.709748Z",
     "shell.execute_reply.started": "2022-09-16T03:14:27.906730Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 162052 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 读取图像\n",
    "img1 = cv2.imread('t1.png', cv2.IMREAD_COLOR)\n",
    "img2 = cv2.imread('t2.png', cv2.IMREAD_COLOR)\n",
    "\n",
    "gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) / 255.0\n",
    "gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY) / 255.0\n",
    "\n",
    "g_images1, corners1, scales1, orientations1, layers1 = detect_blobs(gray1)\n",
    "g_images2, corners2, scales2, orientations2, layers2 = detect_blobs(gray2)\n",
    "\n",
    "# 计算两张图像的特征点\n",
    "descriptors1 = compute_descriptors(g_images1, corners1, scales1, orientations1, layers1)\n",
    "descriptors2 = compute_descriptors(g_images2, corners2, scales2, orientations2, layers2)\n",
    "\n",
    "# 匹配两张图像中的特征点\n",
    "matches = match_descriptors(descriptors1, descriptors2)\n",
    "\n",
    "image1 = cv2.imread('t1.png', cv2.IMREAD_GRAYSCALE)\n",
    "image2 = cv2.imread('t2.png', cv2.IMREAD_GRAYSCALE)\n",
    "\n",
    "# xform为变换矩阵，outlier_labels为模型的outlier\n",
    "xform, outlier_labels = compute_affine_xform(corners1, corners2, matches) \n",
    "# 展示两张图像特征点的匹配\n",
    "match_image = draw_matches(image1, image2, corners1, corners2, matches,outliers=outlier_labels)\n",
    "\n",
    "plt.title('Match')\n",
    "plt.imshow(match_image)\n",
    "plt.axis('off')\n",
    "plt.plot()    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之后，我们将其中一张图像进行变换，并与另一张图像拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:15:33.713310Z",
     "iopub.status.busy": "2022-09-16T03:15:33.712886Z",
     "iopub.status.idle": "2022-09-16T03:15:33.951886Z",
     "shell.execute_reply": "2022-09-16T03:15:33.950889Z",
     "shell.execute_reply.started": "2022-09-16T03:15:33.713272Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 366444 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "image = stitch_images(img1, img2, xform)\n",
    "\n",
    "plt.title('Stich')\n",
    "plt.imshow(image.astype(np.uint8))\n",
    "plt.axis('off')\n",
    "plt.plot()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.5 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本章中，我们学习了基于SIFT特征检测与RANSAC算法的图像拼接方法。图像拼接在实际的应用场景很广，比如无人机航拍，遥感图像等等。这些场景的视角通常很大，需要图像拼接的技术将多次拍摄的图像拼接成一张大图，继而进行后续的图像处理。NASA的火星勇气号探测车也通过图像拼接技术拍摄火星上的全景图。大家如果有兴趣，可以去查阅Matthies所著的《Computer Vision on Mars》一书。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 附录"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用最小二乘法求解\n",
    "$$\n",
    "\\begin{array}{c}\n",
    "\\min _{\\mathbf{p}}\\|A \\mathbf{p}\\|^{2}~~\\text{s.t.}~~\\left \\| \\bar{p} \\right \\| _{2} = 1\n",
    "\\end{array}\n",
    "$$\n",
    "已知:\n",
    "$$\n",
    "\\begin{array}{c}\n",
    "\\|A \\mathbf{p}\\|^{2}=(A \\mathbf{p})^{T}(A \\mathbf{p})=\\mathbf{p}^{T} A^{T} A \\mathbf{p}\\quad ~~\\text{s.t.}~~\\quad\\|\\mathbf{p}\\|^{2}=\\mathbf{p}^{T} \\mathbf{p}=1\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "那么我们可以构建一个损失函数: $\\mathbf{L}_{P}$ 来迭代找到满足方程组的解。\n",
    "$$\n",
    "\\begin{array}{c}\n",
    "\\min _{\\mathbf{p}}\\left\\{L(\\mathbf{p})=\\mathbf{p}^{T} A^{T} A \\mathbf{p}+\\lambda\\left(\\mathbf{p}^{T} \\mathbf{p}-1\\right)\\right\\}\n",
    "\\end{array}\n",
    "$$\n",
    "直接对$P$求导，那么可以得到:\n",
    "$$\n",
    "\\begin{array}{c}\n",
    "A^{T} A \\mathbf{p} + \\lambda\\mathbf{p}=0\n",
    "\\end{array}\n",
    "$$\n",
    "那么显然我们可以得到：矩阵$A^{T}A$的最小特征值$\\lambda$对应的特征向量$\\mathbf{p}$是此方程组的近似解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.6 习题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 习题8.1：为什么我们需要至少四对特征点对才能拟合一个仿射变换？\n",
    "\n",
    "#### 习题8.2：经过本章的学习，相信你已经掌握如何实现两张图像的拼接。那么，如何将多张图像拼接成一张大的全景图像呢？请给出解决思路并编程实现。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
