{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第8章 图像拼接"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当我们游览祖国的大好河山时，我们往往一边感慨山川的奇伟峻丽，一边遗憾自己没法用设备拍一张全景图，而只能拍一些离散的照片。实际上，当你学习完这一章的内容之后，或许你就可以自己实现全景图的生成。\n",
    "\n",
    "\n",
    "图像拼接 (Stitching)，是指将数张描绘同一场景的有重叠的图像（可能是不同时间、不同视角或者不同传感器获得的）拼接成一幅无缝的全景图或者高分辨率图像的技术。图像拼接在运动检测和跟踪、增强现实、分辨率增强、视频压缩和图像稳定等机器视觉领域有很大的应用。如图 8-1 所示，图像拼接的任务是将类似于图中的一系列离散的图像拼接成一张完整的图像，如图 8-2 所示。接下来，就让我们一起学习图像拼接的相关知识吧。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"1.png\" 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=\"2.png\" 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": [
    "在正式介绍图像拼接之前，我们先介绍图像变换的相关工作。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2.1 平移"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像$\\mathbf{I}$上一点$\\mathbf{p}=(i,j)$，平移后的点为$\\mathbf{p}^{\\prime} = (i^{\\prime}, j^{\\prime})$，我们有\n",
    "\n",
    "$$\\mathbf{p}^{\\prime} = \\mathbf{p} + \\mathbf{t}$$\n",
    "\n",
    "其中，$\\mathbf{t} = (t_x, t_y)$为平移方向。平移前后图像的形状不发生变化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2.2 旋转"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将图像$\\mathbf{I}$上一点$\\mathbf{p}=(i,j)$逆时针旋转$\\theta$角度，则对应的新的点坐标为：\n",
    "\n",
    "$$\\mathbf{p}^{\\prime} = \\mathbf{R}\\mathbf{p}$$\n",
    "\n",
    "其中，\n",
    "\n",
    "$$\\mathbf{R} = \\begin{bmatrix}\n",
    "\\cos\\theta &-\\sin\\theta \\\\\n",
    "\\sin\\theta & \\cos\\theta\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "是旋转矩阵。旋转前后图像的长度不发生变化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2.3 相似"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将图像$\\mathbf{I}$上一点$\\mathbf{p}=(i,j)$进行相似变换，则对应的新的点坐标为：\n",
    "\n",
    "$$\\mathbf{p}^{\\prime} = s\\mathbf{R}\\mathbf{p} + \\mathbf{t}$$\n",
    "\n",
    "其中，$s$为常数，$\\mathbf{t}$为平移量。图像相似变换前后角度不变"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2.4 仿射变换"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将图像$\\mathbf{I}$上一点$\\mathbf{p}=(i,j)$进行仿射变换，则对应的新的点坐标为：\n",
    "\n",
    "$$\\mathbf{p}^{\\prime} = \\begin{bmatrix}\n",
    "a_{11} & a_{12} & a_{13} \\\\\n",
    "a_{21} & a_{22} & a_{23}\n",
    "\\end{bmatrix}[\\mathbf{p}, 1]^\\mathbf{T}$$\n",
    "\n",
    "其中，$s$为常数，$\\mathbf{t}$为平移量。仿射变换前后图像中的平行关系不变"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2.5 投影"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将图像$\\mathbf{I}$上一点$\\mathbf{p}=(i,j)$投影到三维空间中，则：\n",
    "\n",
    "$$\\mathbf{p}^{\\prime} = \\mathbf{H}[\\mathbf{p}, 1]^\\mathbf{T}$$\n",
    "\n",
    "这里，$\\mathbf{H}$是一个$3\\times3$的投影矩阵。投影前后图像中的直线仍然为直线。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 图像拼接"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像拼接是将同一场景的多个重叠图像拼接成较大的图像的一种方法，在医学成像、计算机视觉、卫星数据、军事目标自动识别等领域具有重要意义。图像拼接的输出是两个输入图像的并集。图像拼接通常用到以下四个步骤：\n",
    "\n",
    "1. 图像特征提取\n",
    "2. 图像配准\n",
    "3. RANSAC算法计算单应矩阵\n",
    "4. 图像变形与融合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.1 图像特征提取"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们如果要匹配的两个输入图像中的元素，首先需要提取图像的特征点。还记得我们在上一章学过的SIFT算法吗？它是尺度不变的特征点检测算法，可用于识别图像中的相似目标。SIFT算法既具有旋转不变性，又具有尺度不变性，因此非常适合于高分辨率图像中的目标检测。我们可以利用SIFT算法进行图像特征的提取。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.2 图像配准"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在两个输入图像$\\mathbf{I}_1,\\mathbf{I}_2$的特征点被检测出来之后，我们需要以某种方式将它们关联起来，可以通过计算归一化互相关或互信息的方法来确定其对应关系。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.2.1 归一化互相关（Normalized Cross Correlation，NCC）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们在模板匹配中学习过互相关与标准化互相关的相关概念，这里我们进一步介绍另一种形式的互相关————归一化互相关。\n",
    "\n",
    "互相关的工作原理是将第一幅图像中每个点周围的像素窗口与第二幅图像中每个点周围的像素窗口关联起来，将双向相关性最大的点作为对应的“对”。我们定义归一化互相关如下：\n",
    "\n",
    "$$NCC(i,j) = \\frac{\\sum_{m}\\sum_n [\\mathbf{I}_1[m,n] - \\overline{\\mathbf{I}}_1] \\cdot [\\mathbf{I}_2[m+i, n+j] - \\overline{\\mathbf{I}}_2]}{\\sqrt{\\sum_{m}\\sum_n[\\mathbf{I}_1[m,n] - \\overline{\\mathbf{I}}_1]^2}\\sqrt{\\sum_{m}\\sum_n[\\mathbf{I}_2[m+i, n+j] - \\overline{\\mathbf{I}}_2]^2}}$$\n",
    "\n",
    "其中，$\\overline{\\mathbf{I}}_1,\\overline{\\mathbf{I}}_2$是窗口的平均值图像，即：\n",
    "\n",
    "$$\\overline{\\mathbf{I}}_1 = \\frac{1}{MN}\\sum_{i=1}^M\\sum_{j=1}^N \\mathbf{I}_1(i,j)$$\n",
    "\n",
    "$$\\overline{\\mathbf{I}}_2 = \\frac{1}{MN}\\sum_{i=1}^M\\sum_{j=1}^N \\mathbf{I}_2(i,j)$$\n",
    "\n",
    "NCC系数的范围为$[-1,1]$。我们令$u=(i,j)$，NCC峰值相对应的位移参数$u$表示为两个图像之间的几何变换。此方法的优点是计算简单，但是速度特别慢。此外，此类算法要求源图像之间必须有显著的重叠。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.2.2 互信息（Mutual Information，MI）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "互信息测量基于两个图像之间共享信息数量的相似性，互信息越大，说明彼此之间越相似。两个图像$\\mathbf{I}_1,\\mathbf{I}_2$之间的互信息用熵表示为：\n",
    "\n",
    "$$MI(\\mathbf{I}_1,\\mathbf{I}_2) = E(\\mathbf{I}_1) + E(\\mathbf{I}_2) - E(\\mathbf{I}_1\\mathbf{I}_2)$$\n",
    "\n",
    "其中，$E(\\mathbf{I}_1)$和$E(\\mathbf{I}_2)$分别是图像$\\mathbf{I}_1$和$\\mathbf{I}_2$的熵，$E(\\mathbf{I}_1\\mathbf{I}_2)$是两个图像的联合熵。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3.3 单应矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "单应矩阵估计是图像拼接的第三步。在单应矩阵估计中，不属于重叠区域的不需要的特征点被删除。我们采用RANSAC算法进行单应。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.3.1 RANSAC 算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随机样本一致算法RANSAC(random sample consensus)从可能含有异常值的观测数据集中拟合数学模型，是一种鲁棒参数估计的迭代方法。该算法是不确定性的，因为它只在一定的概率下产生一个合理的结果，随着执行更多的迭代，这个概率会增加。RANSAC算法用于在存在大量可用数据外行的情况下以鲁棒的方式拟合模型，在计算机视觉中有许多应用。\n",
    "\n",
    "算法的流程可以简要地描述如下：\n",
    "\n",
    "1. 从数据样本$S$中随机选取$K$个点，得到假设子集$S_K$，然后使用$S_K$优化模型参数$M(\\theta)$，得到初始模型；\n",
    "2. 对所有数据样本$S$应用模型$M(\\theta)$，根据每个数据样本的误差$||M(\\theta, x_i)-y_i||$来所有数据样本划分为内群点（inliner points）和离群点（outliner points），其中误差大于超参数$\\epsilon$的称为离群点$S_{outliner}$，反之则为内群点$S_{inliner}$；\n",
    "3. 如果内群点的占比$R = \\frac{S_{inliner}}{S}$大于超参数阈值$\\tau$，则使用所有的内群点$S_{inliner}$来重新优化模型参数，并将得到的模型作为最终模型，终止迭代；\n",
    "4. 否则使用内群点$S_{inliner}$优化模型，并重复执行步骤 2，直到条件3满足，或者达到迭代次数上限停止。\n",
    "\n",
    "图 8-3 是使用RANSAC拟合数据的一个流程。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"3.png\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图8-3 RANSAC示意</div>\n",
    "</center>\n",
    "\n",
    "可以看到 RANSAC 算法的性能由超参数$N,K,\\epsilon,\\tau$来决定，通常来讲，迭代次数必须设置的足够大，才能使得至少有一次数据采样$S_K$不包含任何离群点的概率$p$（一般来说$p=0.99$）满足我们的使用要求，以此来避免离群点的影响。我们可以通过简单地推导来得到迭代上限的近似估计:\n",
    "\n",
    "假设每次随机选择一个点，该点为内群点的概率为$p_{i}$，则其为离群点的概率为$p_o = 1-p_i$ ，则在$N$次迭代，每次迭代选取$K$个点的情况下，有\n",
    "\n",
    "$$1-p = (1-p_i^K)^N$$\n",
    "\n",
    "从而估计出\n",
    "\n",
    "$$N = \\frac{\\log(1-p)}{\\log(1-(1-p_i)^K)}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.3.2 计算单应矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有了两组相关点，接下来就需要建立两组点的转换关系，也就是图像变换关系。单应矩阵是两个空间之间的映射，常用于表示同一场景的两个图像之间的对应关系。我们可以利用单应矩阵匹配大部分相关的特征点，并实现图像投影，使一张图通过投影和另一张图重合。\n",
    "\n",
    "假设两个图像的匹配点为$\\mathbf{X} = [x,y]^\\mathbf{T}, \\mathbf{X}^{\\prime} = [x^{\\prime},y^{\\prime}]^\\mathbf{T}$\n",
    "\n",
    "我们有\n",
    "\n",
    "$$\\mathbf{X}^{\\prime} = \\mathbf{H}[\\mathbf{X},1]$$\n",
    "\n",
    "由于$\\mathbf{X}$和$\\mathbf{X}^{\\prime}$两向量共线，我们有\n",
    "\n",
    "$$[\\mathbf{X}^{\\prime},1]^\\mathsf{T}\\mathbf{H}[\\mathbf{{X}},1] = 0$$\n",
    "\n",
    "令$\\mathbf{H} = \\begin{bmatrix}\n",
    "h_{11} & h_{12} & h_{13} \\\\\n",
    "h_{21} & h_{22} & h_{23} \\\\\n",
    "h_{31} & h_{32} & 1\n",
    "\\end{bmatrix}$为一个含有8个参数的变换矩阵，由此我们知道四点可以确定$\\mathbf{H}$，即：\n",
    "\n",
    "\\begin{equation}\n",
    "\\left[\\begin{array}{c}\n",
    "x_{1}, y_{1}, 1,0,0,0,-x_{1}^{\\prime} x_{1},-x_{1}^{\\prime} y_{1} \\\\\n",
    "0,0,0, x_{1}, y_{1}, 1,-y_{1}^{\\prime} x_{1},-y_{1}^{\\prime} y_{1} \\\\\n",
    "x_{2}, y_{2}, 1,0,0,0,-x_{2}^{\\prime} x_{2},-x_{2}^{\\prime} y_{2} \\\\\n",
    "0,0,0, x_{2}, y_{2}, 1,-y_{2}^{\\prime} x_{2},-y_{2}^{\\prime} y_{2} \\\\\n",
    "x_{3}, y_{3}, 1,0,0,0,-x_{3}^{\\prime} x_{3},-x_{3}^{\\prime} y_{3} \\\\\n",
    "0,0,0, x_{3}, y_{3}, 1,-y_{3}^{\\prime} x_{3},-y_{3}^{\\prime} y_{3} \\\\\n",
    "x_{4}, y_{4}, 1,0,0,0,-x_{4}^{\\prime} x_{4},-x_{4}^{\\prime} y_{4} \\\\\n",
    "0,0,0, x_{4}, y_{4}, 1,-y_{4}^{\\prime} x_{4},-y_{4}^{\\prime} y_{4}\n",
    "\\end{array}\\right]\\left[\\begin{array}{l}\n",
    "h_{12} \\\\\n",
    "h_{21} \\\\\n",
    "h_{22} \\\\\n",
    "h_{23} \\\\\n",
    "h_{31} \\\\\n",
    "h_{32}\n",
    "\\end{array}\\right]=\\left[\\begin{array}{c}\n",
    "x_{1}^{\\prime} \\\\\n",
    "y_{1}^{\\prime} \\\\\n",
    "x_{2}^{\\prime} \\\\\n",
    "y_{2}^{\\prime} \\\\\n",
    "x_{3}^{\\prime} \\\\\n",
    "y_{3}^{\\prime} \\\\\n",
    "x_{4}^{\\prime} \\\\\n",
    "y_{4}^{\\prime}\n",
    "\\end{array}\\right]\n",
    "\\end{equation}\n",
    "\n",
    "那么就可以每次从所有的匹配点中选出4对，计算单应性矩阵H，然后选出内点个数最多的作为最终的结果。其中，特征点的误差定义为：\n",
    "\n",
    "$$Dist(\\mathbf{HX}、\\mathbf{X}^{\\prime}) = \\left\\|\\left[\\begin{array}{c}\n",
    "x^{\\prime} \\\\\n",
    "y^{\\prime} \\\\\n",
    "1\n",
    "\\end{array}\\right]-\\mathbf{H}\\left[\\begin{array}{c}\n",
    "x \\\\\n",
    "y \\\\\n",
    "1\n",
    "\\end{array}\\right]\\right\\|$$\n",
    "\n",
    "我们利用用RANSAC方法估算$\\mathbf{H}$的过程如下：\n",
    "\n",
    "1. 首先检测两张图像的特征点；\n",
    "2. 在特征点之间应用归一化互相关，收集相关性足够高的对；\n",
    "3. 使用匹配的对计算矩阵$\\mathbf{H}$；\n",
    "4. 选择与单应性一致的配对。如果有$Dist(\\mathbf{HX}、\\mathbf{X}^{\\prime}) <\\epsilon$，则点对$(\\mathbf{X},\\mathbf{X}^{\\prime})$被认为与单应性H一致\n",
    "5. 重复34步，直到足够多的点对满足$\\mathbf{H}$；\n",
    "6. 使用所有满足条件的点对，通过公式重新计算$\\mathbf{H}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.4 图像变形与融合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后一步是将所有输入图像变形并融合到一个符合的输出图像中。我们可以简单地将所有输入的图像变形到一个平面上，这个平面名为复合全景平面。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.4.1 图像变形步骤"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 计算每个输入图像的变形图像坐标范围，得到输出图像大小。我们可以很容易地通过映射每个源图像的四个角确定输出图像的大小。最后，需要计算指定参考图像原点相对于输出全景图的偏移量的偏移量；\n",
    "2. 使用上述的变形，将每个输入图像的像素映射到参考图像定义的平面。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3.4.2 图形融合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当图像变形完成并展示在一张图像上之后，我们需要在重叠区域融合像素颜色以避免接缝。最简单的可用形式是使用羽化，它使用加权平均颜色值融合重叠的像素。我们通常使用$\\alpha$因子（通常称为$\\alpha$通道），它在中心像素处的值为1，在与边界像素线性递减后变为0。当输出拼接图像中至少有两幅重叠图像时，我们将使用如下的$\\alpha$值来计算其中一个像素处的颜色.\n",
    "\n",
    "假设两个图像$\\mathbf{I}_1,\\mathbf{I}_2$在输出图像中重叠，每个像素点$(x,y)$在图像中表示为$\\mathbf{I}_i(x,y) = (\\alpha_iR, \\alpha_iG, \\alpha_iB, \\alpha_i)$，其中$(R,G,B)$是像素的颜色值。我们将在缝合后的输出图像中计算$(x,y)$的像素值为：\n",
    "\n",
    "$$\\frac{[(\\alpha_1R, \\alpha_1G, \\alpha_1B, \\alpha_1)+[(\\alpha_2R, \\alpha_2G, \\alpha_2B, \\alpha_2)]}{\\alpha_1+\\alpha_2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 编程实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一节中，我们将编程实现图像拼接。我们先导入上一章中编写的SIFT相关函数，并在此基础上利用RANSAC算法实现图像拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import *\n",
    "def draw_matches(image1, image2, corners1, corners2, matches,\n",
    "    outlier_labels=None):\n",
    "    \"\"\"Draw matched corners between images.\n",
    "    \n",
    "    Args:\n",
    "    - matches (list of 2-tuples)\n",
    "    - image1 (3D uint8 array): A color image having shape (H1, W1, 3).\n",
    "    - image2 (3D uint8 array): A color image having shape (H2, W2, 3).\n",
    "    - corners1 (list of 2-tuples)\n",
    "    - corners2 (list of 2-tuples)\n",
    "    - outlier_labels (list of bool)\n",
    "    \n",
    "    Returns:\n",
    "    - match_image (3D uint8 array): A color image having shape\n",
    "        (max(H1, H2), W1 + W2, 3).\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",
    "    match_image = np.zeros((h1, w1 + w2, 3), np.uint8)\n",
    "\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",
    "        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",
    "        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",
    "        if outlier_labels[i] == 1:\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\n",
    "\n",
    "\n",
    "# 利用RANSAC算法剔除无效点\n",
    "def compute_affine_xform(corners1, corners2, matches):\n",
    "    \"\"\"Compute affine transformation given matched feature locations.\n",
    "    \n",
    "    Args:\n",
    "    - corners1 (list of 2-tuples)\n",
    "    - corners1 (list of 2-tuples)\n",
    "    - matches (list of 2-tuples)\n",
    "    \n",
    "    Returns:\n",
    "    - xform (2D float64 array): A 3x3 matrix representing the affine\n",
    "        transformation that maps coordinates in image1 to the corresponding\n",
    "        coordinates in image2.\n",
    "    - outlier_labels (list of bool): A list of Boolean values indicating whether\n",
    "        the corresponding match in `matches` is an outlier or not. For example,\n",
    "        if `matches[42]` is determined as an outlier match after RANSAC, then\n",
    "        `outlier_labels[42]` should have value `True`.\n",
    "    \"\"\"\n",
    "    iteration = 50\n",
    "    M_list = []\n",
    "    inlier_num_list = []\n",
    "    for _ in range(iteration):\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 = 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",
    "        _,_, v = np.linalg.svd(A)\n",
    "        M = np.reshape(v[-1], (3, 3))\n",
    "        \n",
    "        #M = t.reshape((3,3))\n",
    "        inlier_num = 0\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(M, coord1)\n",
    "            mapcoor = mapcoor/mapcoor[-1]\n",
    "            if np.linalg.norm(coord2 - mapcoor) < 5:\n",
    "                inlier_num += 1\n",
    "        M_list.append(M)\n",
    "        inlier_num_list.append(inlier_num)\n",
    "        \n",
    "    best_index = np.argmax(inlier_num_list)\n",
    "    xform = M_list[best_index].astype(np.float64)\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  \n",
    "\n",
    "\n",
    "def stitch_images(image1, image2, xform):\n",
    "    \"\"\"Stitch two matched images given the transformation between them.\n",
    "    \n",
    "    Args:\n",
    "    - image1 (3D uint8 array): A color image.\n",
    "    - image2 (3D uint8 array): A color image.\n",
    "    - xform (2D float64 array): A 3x3 matrix representing the transformation\n",
    "        between image1 and image2. This transformation should map coordinates\n",
    "        in image1 to the corresponding coordinates in image2.\n",
    "    \n",
    "    Returns:\n",
    "    - image_stitched (3D uint8 array)\n",
    "    \"\"\"\n",
    "    new_image = cv2.warpPerspective(image1, xform, (image2.shape[1],image2.shape[0]))\n",
    "    image = new_image.copy()\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",
    "    image = cv2.polylines(image2, [np.int32(dst)], True, 255, 1, cv2.LINE_AA)\n",
    "    image = np.where(new_image>0, 0.5*new_image+0.5*image, image)\n",
    "    return image.astype(np.uint8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 78868 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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",
    "descriptors1 = compute_descriptors(g_images1, corners1, scales1, orientations1, layers1)\n",
    "descriptors2 = compute_descriptors(g_images2, corners2, scales2, orientations2, layers2)\n",
    "    \n",
    "matches = match_descriptors(descriptors1, descriptors2)\n",
    "image1 = cv2.imread('t1.png', cv2.IMREAD_GRAYSCALE)\n",
    "image2 = cv2.imread('t2.png', cv2.IMREAD_GRAYSCALE)\n",
    "\n",
    "xform, outlier_labels = compute_affine_xform(corners1, corners2, matches)    \n",
    "match_image = draw_matches(image1, image2, corners1, corners2, matches,outlier_labels=outlier_labels)\n",
    "\n",
    "plt.title('Match')\n",
    "plt.imshow(match_image)\n",
    "plt.axis('off')\n",
    "plt.plot()    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 162932 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "image1 = cv2.imread('t1.png')\n",
    "image2 = cv2.imread('t2.png')\n",
    "    \n",
    "image = stitch_images(image1, image2, xform)\n",
    "\n",
    "plt.title('Stich')\n",
    "plt.imshow(image.astype(np.uint8))\n",
    "plt.axis('off')\n",
    "plt.plot()    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.5 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本章中，我们学习了利用如何进行图像拼接。我们首先利用SIFT算法提取出需要拼接的图像中的特征点，在此基础上使用RANSAC算法剔除无效点，之后进行匹配与图像变形融合，得到最终的拼接图像。由于篇幅的限制我们只做了两张图像的拼接，你也可以自己动手尝试多张图像的拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}
