{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第7章 SIFT特征检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1 块状区域检测与尺度空间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "仍然回到上一章所提到的图像匹配问题。如果需要去匹配图7-1左上和右上两个子图，要怎么做才能实现呢？是否可以通过上一章介绍的角点来进行匹配呢？注意到图7-1的左上和右上两个子图除了有旋转变化还有尺度变化，而角点对于尺度并不是不变的，因此需要寻找一种新的特征点（关键点），用于匹配复杂变化下的图像。如果以这种特征点为中心能够定位到图像中一些感兴趣（interesting）的区域，比如图7-1左下和右下所示的怪兽的眼睛，并且如果知道这两个区域是匹配的，那么就可以计算出它们之间的旋转和尺度变化，从而可以匹配两张图像。但是图像中并不是每个区域都是可以用于匹配的特征点区域，这样的区域需要满足这样几种性质：\n",
    "\n",
    "<center>\n",
    "      <img style=\"border-radius: 0.3125em;\" src=\"https://pic1.zhimg.com/80/v2-3c70adb701efda7e4365723fd6d4b018_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-1 基于感兴趣的区域特征去匹配图像</div>\n",
    "</center>\n",
    "\n",
    "1. 具有丰富的语义信息；\n",
    "2. 在图像中具有明确的位置（可定位）；\n",
    "3. 具有明确的特征描述，可用于匹配或区分其他感兴趣区域。\n",
    "\n",
    "而对于区域的特征描述，还有以下要求：\n",
    "\n",
    "1. 对图像的平移旋转和放缩具有不变性；\n",
    "2. 对光照变化具有鲁棒性。\n",
    "\n",
    "\n",
    "接下来，根据上面的性质来判断图7-2中的哪些区域是感兴趣的特征点区域。首先，图 7-2 中的边缘是特征点区域吗？我们发现不同位置的边缘区域的特征基本类似，因此无法对其准确的定位，所以边缘区域不是特征点区域。而图 7-2 中被闭合轮廓包围的块状区域（blob）具有确定的位置和大小，以及特有的语义信息，可以用于匹配，所以块状区域是感兴趣的特征点区域。\n",
    "\n",
    "<center>\n",
    "        <img style=\"border-radius: 0.3125em;\" \n",
    "   src=\"https://pic2.zhimg.com/80/v2-f984e685ebc4e625e126455197784ca9_1440w.jpg\" width=700>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-2 图像中的边缘区域和块状区域</div>\n",
    "</center>\n",
    "\n",
    "那么检测图像中的特征点，就需要检测图像中的块状区域，确定其位置、大小和方向，并最终将其表述为一个与大小和方向无关的特征描述。由于块状区域是具有闭合轮廓（边缘）的区域，自然会想到用边缘检测器来试试对其进行检测。图7-3显示了一些不同大小的一维块状区域，试试用二阶微分算子（拉普拉斯算子）对其卷积，并观察效果。注意，在第5章中我们学习过，做边缘检测前都有一步用高斯平滑去噪的过程，这里仍然沿用这一步骤。结合了高斯平滑和二阶微分的算子叫做高斯拉普拉斯算子（Laplacian of Gaussian, LoG）。\n",
    "\n",
    "<center>\n",
    "   <img style=\"border-radius: 0.3125em;\" \n",
    " src=\"https://pic3.zhimg.com/80/v2-2e564f910382e558dfaebbdf7c71e82a_1440w.webp\" width=800>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-3 块状区域对高斯拉普拉斯算子的响应</div>\n",
    "</center>\n",
    "\n",
    "图 7-3 展示了不同大小的一维块状区域与不同标准差的高斯拉普拉斯算子的卷积结果。不难发现，随着高斯拉普拉斯算子的标准差的增大，当该标准差的大小和块状区域大小相当时，它们卷积结果会出现唯一的极值点。此时的标准差大小称为特征尺度$\\sigma^\\ast$，与块状区域大小相同。既然存在这样一个响应极值，是否可以用不同大小标准差的高斯拉普拉斯算子对输入图像进行卷积，然后对卷积结果求取极值就可以找到输入图像的块状区域了呢？然而事情并没有这么简单。注意到前面所说的唯一极值点是在某一特定标准差$\\sigma$（尺度）下$x$空间上的唯一极值点，如果考虑多个尺度，即在整个$(x,\\sigma)$空间上，如图 7-4 上子图所示，该极值并不是一个全局最小值。这使得上述对块状区域检测的方法无法实施。出现这样现象的原因是因为高斯拉普拉斯算子的响应会随着标准差变大而衰减。所以需要对响应做归一化，即对其乘以$\\sigma^2$。这样归一化后的高斯拉普拉斯算子的响应如图 7-4 下子图所示，此时特征尺度下的极值点是整个$(x,\\sigma)$空间上的全局最小值。那么只要检测到$(x,\\sigma)$空间上归一化高斯拉普拉斯算子响应的全局最小值，就可以检测到块状区域，并同时确定其特征尺度。\n",
    "\n",
    "<center>\n",
    "   <img style=\"border-radius: 0.3125em;\" \n",
    " src=\"https://pic3.zhimg.com/80/v2-3d031c91bd563911178af65d108546be_1440w.jpg\" width=800>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-4 高斯拉普拉斯算子的尺度归一化响应</div>\n",
    "</center>\n",
    "\n",
    "至此，我们学习了基于尺度空间的块状区域检测的基本原理。那么，如何将该原理应用到图像中呢？\n",
    "\n",
    "图像的尺度空间$\\mathbf{L}(x,y,\\sigma)$可通过一系列标准差从小到大的高斯函数$\\mathbf{G}(x,y,\\sigma)$与输入图像$\\mathbf{I}(x,y)$卷积得到，即：\n",
    "\n",
    "$$\\mathbf{L}(x,y,\\sigma) = \\mathbf{G}(x,y,\\sigma) * \\mathbf{I}(x,y)$$\n",
    "\n",
    "其中，$(x,y)$是二维空间坐标，$\\sigma$是高斯函数的标准差，也称为尺度空间因子。根据在滤波中学习的知识，我们知道$\\sigma$值越小，即尺度越小，图像模糊的程度越小；$\\sigma$值越大，即尺度越大，图像模糊的程度越大。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。计算图像尺度空间与归一化拉普拉斯算子卷积的结果，根据卷积结果中是否存在极值点来判断当前位置是否是块状区域，即是否存在特征点。如图 7-5 所示，当前位置（屋檐位置），其卷积响应存在极大值，说明当前位置是一个块状区域，该极大值对应的尺度为这个块状区域的特征尺度。如图7-6所示，当前位置，其卷积响应不存在明显的极值，说明当前位置不是一个块状区域，不存在特征点。\n",
    "\n",
    "<center>\n",
    "   <img style=\"border-radius: 0.3125em;\" \n",
    " src=\"https://pic3.zhimg.com/80/v2-897d5db7f80152a8818bc286dd9e0ce6_1440w.webp\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-5 图像尺度空间的归一化拉普拉斯算子响应中存在极值点，说明存在块状区域，即存在特征点。</div>\n",
    "</center>\n",
    "\n",
    "<center>\n",
    "   <img style=\"border-radius: 0.3125em;\" \n",
    " src=\"https://pic3.zhimg.com/80/v2-9dd7fd80bd3e5de696a801f4d59568ee_1440w.webp\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-6 图像尺度空间的归一化拉普拉斯算子响应中不存在明显极值点，说明不存在块状区域，即不存在特征点。</div>\n",
    "</center>\n",
    "\n",
    "如果找到了两张图像中对应的特征点，根据特征点对应的特征尺度的比例，就可以得到两张图像尺度缩放的比例。如图 7-7 中的两张图像，他们之间对应的放缩比例为$\\sigma_1^\\ast/\\sigma_2^\\ast$\n",
    "\n",
    "<center>\n",
    "   <img style=\"border-radius: 0.3125em;\" \n",
    " src=\"https://pic2.zhimg.com/80/v2-3604abee859ae87813475b9ee2c2f1dd_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-7 特征尺度的比例与图像尺度变化的关系。</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "\n",
    "接下来，我们将学习一种高效的基于尺度空间的块状区域检测算法——SIFT。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2 SIFT算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SIFT [1] 的全称是 Scale Invariant Feature Transform（尺度不变特征变换），由加拿大教授David G.Lowe于1999年在国际计算机视觉会议（IEEE International Conference on Computer Vision，ICCV）上提出，是计算机视觉界近二十年来引用率最高的文章之一。SIFT是一种高效的基于尺度空间的块状区域检测算法，包含块状区域的检测及块状区域的描述，即SIFT特征。它的提出极大推动了计算机视觉领域的发展，并被广泛用于视觉识别和三维重建。SIFT算法的主要思想同7.1中描述的一致，即利用尺度空间来检测图像中的块状区域。接下来，将详细介绍SIFT算法。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.1 局部极值点检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在前文中，我们已经知道如何利用尺度空间检测特征点，即利用尺度归一化的高斯拉普拉斯算子$\\sigma^2 \\nabla^2\\mathbf{G}$与图像$\\mathbf{I}$做卷积，寻找响应中极值点对应的$\\sigma$。为了提高计算效率，我们可以使用高斯差分算子（Difference of Gaussian，DOG）$\\mathbf{D}$近似代替表示归一化高斯拉普拉斯算子。DOG算子定义如下：\n",
    "\n",
    "$$\\mathbf{D}(x,y,\\sigma)=\\mathbf{G}(x,y,k\\sigma) - \\mathbf{G}(x,y,\\sigma) $$\n",
    "\n",
    "因为\n",
    "\n",
    "$$\\sigma \\nabla^2\\mathbf{G} = \\frac{\\partial \\mathbf{G}}{\\partial \\sigma} \\approx \\frac{\\mathbf{G}(x,y,k\\sigma) - \\mathbf{G}(x,y,\\sigma)}{k\\sigma - \\sigma} $$\n",
    "\n",
    "因此，我们可以得到：\n",
    "\n",
    "$$\\mathbf{D}(x,y,\\sigma)=\\mathbf{G}(x,y,k\\sigma) - \\mathbf{G}(x,y,\\sigma) \\approx {(k-1)\\sigma^2}\\nabla^2\\mathbf{G}$$\n",
    "\n",
    "其中k-1是个常数，对极值点不会有影响。根据图 7-8 所示，我们也可以发现归一化高斯拉普拉斯算子和高斯差分算子近似相同。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-c541d33c0b2628b47a7cca0f50887f59_1440w.webp\" width=300>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-8 高斯拉普拉斯算子和高斯差分算子的比较</div>\n",
    "</center>\n",
    "\n",
    "用高斯差分算子$\\mathbf{D}(x,y,\\sigma)$对图像$\\mathbf{I}(x,y)$作卷积，得到：\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf{D}(x,y,\\sigma) * \\mathbf{I}(x,y) & = (\\mathbf{G}(x,y,k\\sigma) - \\mathbf{G}(x,y,\\sigma))* \\mathbf{I}(x,y) \\\\\n",
    "& = \\mathbf{L}(x,y,k\\sigma)-\\mathbf{L}(x,y,\\sigma)\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "由此可见，高斯差分算子与图像卷积的结果可以由不同尺度的$\\mathbf{L}$直接相减得到，因此，可以通过如图 7-9 所示的方法来计算高斯差分算子与图像卷积的结果。首先得到一张图像对应的尺度空间，即通过由小到大不同大小标准差（尺度）的高斯函数与输入图像做卷积得到一系列响应图像，并将尺度大小相邻的两幅响应图像做差，即可得到相应的不同大小尺度的高斯差分算子与输入图像卷积的结果。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic4.zhimg.com/80/v2-18fefd08b31842f757a7ad6ea2d43143_1440w.webp\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-9 基于高斯差分算子快速计算图像尺度空间。</div>\n",
    "</center>\n",
    "\n",
    "如前所述，特征点对应于高斯差分算子响应图上的极值点。为了检测响应图上的极值点，需要将响应图上每个点的响应值与它所有相邻点比较。如图 7-10 中间子图所示，需要将中间橙色的待检测点与它在整个$(x,y,\\sigma)$图像尺度空间上的26邻域点的响应值进行比较，包括与它在同一尺度响应图上的8个相邻点以及两个相邻尺度响应图上对应的$9\\times2$个相邻点。如果该点的响应值在26领域点中是最大的，则认为该点是候选特征点。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic4.zhimg.com/80/v2-ddad446b6005bb9666ef523e26840583_1440w.webp\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-10 图像尺度空间上26邻域中的局部极值点为候选特征点。</div>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.2 特征点定位与筛选"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在章节7.2.1中，通过尺度空间中极值点检测，得到候选特征点。注意，“候选”两个字说明有些检测到的特征点并不是真正的特征点，需要被进一步筛选。那么哪些候选特征点需要被去除呢？由于这些特征点是离散的尺度空间中的极值点，而离散空间中的极值点与连续空间中的极值点是有差别的，如图 7-11 所示：\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-c6c065ab1b1609498ff9d2ff5fd5f369_1440w.webp\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-11 离散空间与连续空间极值点的差别 </div>\n",
    "</center>\n",
    "\n",
    "因此，为了从候选特征点中筛选出真正需要的特征点，可以利用插值的方法对其进行过滤。首先对高斯差分算子进行曲线拟合，该算子的泰勒展开为：\n",
    "\n",
    "$$\\mathbf{D}(\\mathbf{x}) = \\mathbf{D}(\\mathbf{0})+ \\frac{\\partial\\mathbf{D}^\\mathsf{T}}{\\partial \\mathbf{x}}\\mathbf{x} +\\frac{1}{2}\\mathbf{x}^\\mathbf{T} \\frac{\\partial^2\\mathbf{D}}{\\partial \\mathbf{x}^2}\\mathbf{x}$$\n",
    "\n",
    "其中，$\\mathbf{x} = (x,y,\\sigma)^\\mathsf{T}$。令方程求导为0，可以得到真实极值点位置相对于当前位置的偏移量：\n",
    "\n",
    "$$\\hat{\\mathbf{x}} = -\\frac{\\partial^2\\mathbf{D}}{\\partial \\mathbf{x}^2}\\frac{\\partial\\mathbf{D}}{\\partial \\mathbf{x}}$$\n",
    "\n",
    "一般认为，当$\\hat{\\mathbf{x}}$在任一维度上的偏移量大于0.5（即$x,y$或$\\sigma$）时，这意味着插值中心已经偏移到它的邻近点上，所以必须更新当前候选特征点的位置，并在新的位置上反复插值直到收敛。在实际操作中这一过程也有可能超出所设定的迭代次数或者更新的极值点的位置超出图像边界的范围，此时这样的点应该删除。\n",
    "\n",
    "代入极值点$\\mathbf{x}^\\ast = \\mathbf{x}+\\hat{\\mathbf{x}}$，方程的值为：\n",
    "\n",
    "$$\\mathbf{D}(\\mathbf{x}^\\ast) = \\mathbf{D}(\\mathbf{0}) +\\frac{1}{2}\\frac{\\partial\\mathbf{D}^\\mathsf{T}}{\\partial \\mathbf{x}}\\hat{\\mathbf{x}}$$\n",
    "\n",
    "我们把$\\mathbf{D}({\\mathbf{x}^\\ast})$也称为该点的对比度。由于过小的点容易受噪声的干扰而变得不稳定，因此将低对比度（小于某个经验值如0.03）的极值点删除。\n",
    "\n",
    "除此之外，还需要消除边缘响应（Edge Response）。还记得在第5章中学习的边缘检测的相关知识吗？边缘处像素点的梯度方向与边缘垂直，在边缘方向上没有分量。这意味着对于边缘的像素点而言，能找到两个方向，一个方向的图像梯度变化最快，而另一个方向最慢。这就是所谓的边缘响应。应该尽量筛掉边缘的候选特征点，因为他们只有一个方向有极值，很难被唯一地描述。\n",
    "\n",
    "那么该如何消除边缘响应呢？还记得在第6章角点检测中，对于图像中某一位置$[i,j]$，可以用角点响应函数确定其是否为边缘。令$\\lambda_1$和$\\lambda_2$是$[i,j]$位置对应的Hessian矩阵$\\mathbf{H}$的两个特征值。当其中一个特征值远大于另一个时，认为该位置是边缘。假设$\\lambda_1=\\gamma\\lambda_2$，有\n",
    "\n",
    "$$\\frac{{\\rm Tr}(\\mathbf{H})^2}{{\\rm Det}(\\mathbf{H})} = \\frac{(\\lambda_1+\\lambda_2)^2}{\\lambda_1 \\lambda_2} = \\frac{(\\gamma+1)^2 \\lambda_2^2}{\\gamma \\lambda_2^2} = \\frac{(\\gamma+1)^2}{\\gamma}$$\n",
    "\n",
    "根据不等式的相关知识我们知道，上述公式的值在$\\gamma=1$（即两个特征值相等）时最小，而随着$\\gamma$的增大而增大。$\\gamma$越大即表示两个特征值差距越大，越有可能是边缘。因此，为了检测候选特征点是否为边缘点，只需要检验下式是否被满足：\n",
    "\n",
    "$$\\frac{{\\rm Tr}(\\mathbf{H})^2}{{\\rm Det}(\\mathbf{H})}<\\frac{(\\gamma+1)^2}{\\gamma}$$\n",
    "\n",
    "根据经验，$\\gamma=10$。这样，经过筛选后的特征点即为最终确定的特征点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.3 特征点方向计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了能够匹配不同旋转角度下的相似目标，特征点应具有旋转不变性，所以需要确定特征点的主方向。那么特征点的主方向如何确定呢？在第5章的边缘检测中我们介绍过图像梯度可以描述图像局部结构的方向，那么也可以基于这个思想去计算特征点的主方向。根据章节7.2.1，可以得到每个特征点的位置以及特征尺度，根据Lowe的建议，以特征点的位置为中心，4.5倍的特征尺度（即1.5倍大小的3$\\sigma$原则）为半径的圆的外接正方形区域即是该特征点对应的块状区域。完成关键点梯度计算后，使用直方图统计关键点邻域内像素的梯度幅值和方向。具体做法是，将360°分为36柱，每10°为一柱，然后在该正方形区域内，将梯度方向在某一个柱内的像素挑选出，并将他们的幅值相加作为柱的高度。因为在正方形区域内像素的梯度幅值对中心像素的贡献是不同的，还需要对幅值进行加权处理，这里采用高斯加权。若当前像素点与中心像素点的偏移为$[i,j]$，则对应的高斯权重$w$为\n",
    "\n",
    "$$w = e^{\\left(-\\frac{1}{2 \\times(1.5\\sigma)^2)}(i^2+j^2 )\\right)}$$\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-d3696e0187b7a7e4214c7544f9124fd6_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-12 直方图统计选取特征点方向 </div>\n",
    "</center>\n",
    "\n",
    "在统计完成各个方向的加权梯度幅值得到直方图$H$之后，为了防止某个某个梯度方向角度因受到噪声的干扰而改变，还需要对梯度方向直方图进行平滑处理，平滑公式为\n",
    "\n",
    "$$H[i] = \\frac{6}{16}H[i] + \\frac{4}{16}(H[i-1]+H[(i+1)]) + \\frac{1}{16}(H[i-2]+H[i+2])$$\n",
    "\n",
    "由于角度是循环的（即$0^\\circ  =360^\\circ$），若上述公式中索引越界，则可以通过圆周循环的方法找到其所对应的$[0^\\circ,360^\\circ]$之间的位置，如$H[38] = H[2]$。平滑处理后直方图的峰值代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。如在图 7-12 的直方图中的正西方向代表了关键点的主方向（为简化，图中只画了八个方向的直方图）。\n",
    "\n",
    "除了给每个特征点分配一个主方向外，还可能会给其分配一个或更多个辅方向。增加辅方向是为了增强图像匹配的鲁棒性，辅方向的定义是，当存在另一个柱体高度大于主方向柱体高度的80%的柱体时，该柱体所代表的方向就是该特征点的辅方向。\n",
    "\n",
    "至此，检测出的图像SIFT特征点含有位置、尺度和方向三种信息。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.4 特征点描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来需要为每个特征点建立一个描述符，即用一组向量表示这个特征点。这个描述符需要满足1）描述符具有独特性，即相似的特征点具有相似的描述符，不相似的特征点具有不同的描述符，从而可以基于描述符对不同图像中的特征点进行匹配；2）描述符不随各种变化而改变，比如光照变化、视角变化等等，从而保证匹配不受这些变化的影响。 \n",
    "\n",
    "为了生成描述符，对于每个特征点，以该点为中心，以半径$r=\\frac{3\\sigma \\sqrt{2}(d+1)}{2}$作圆，确定该圆的外接正方形区域，并将该区域划分成$d\\times d$个子区域。根据经验，$d=4$。如图 7-13 所示，在每个子区域内对8个方向上的图像梯度幅值（8个方向为东，西，南，北，东南，东北，西南，西北）做直方图统计。梯度幅值可以经过高斯加权，也可以不加权。这样每个子区域可以用长度为8的直方图向量表示。最后将16个子区域的向量连接构成的长度为128维的向量就是该特征点的描述符，也称为SIFT特征向量。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-2c29dcd83b5f3090fe4ecd1f63fcc3d2_1440w.webp\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-13 SIFT特征向量的构建。 </div>\n",
    "</center>\n",
    "\n",
    "然而，上述操作得到的SIFT特征向量并不具有旋转不变性，因此，在进行梯度幅值计算之前需要将每个特征点所在的块状区域的横轴方向旋转到主方向上（如图 7-14 所示）。这种旋转并不会改变特征点的邻居像素，也不会改变邻域内像素点梯度的幅值，它仅仅改变了每个像素点的梯度方向。思考一下，为什么通过将块状区域的横轴方向旋转到主方向便可以实现特征向量的旋转不变性？这是因为，对特征点进行特征描述时，需要利用到横纵轴方向进行梯度方向的确定。如果该方向和图像本身相关，那么当图像旋转时，该坐标轴也会随之而旋转，得到的描述符也就不同。当该方向仅仅依赖于块状区域的主方向时，无论对图像进行多大角度的旋转，计算的特征描述符都不会改变。换言之，依赖的坐标系并不是绝对的，而是相对的。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-ee9c7d8cb6e0595b1f855058d5048aed_1440w.webp\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\"> 图7-14 特征点所在的块状区域的横轴旋转到主方向上。 </div>\n",
    "</center>\n",
    "\n",
    "对块状区域中的一点$[i,j]$，将其逆时针旋转$\\alpha$角度，得到的新的坐标点坐标$[i^{'},j^{'}]$为\n",
    "\n",
    "\n",
    "$$[i^{'},j^{'}]^\\mathsf{T} =\\begin{bmatrix}\n",
    "\\cos\\alpha &-\\sin\\alpha \\\\\n",
    "\\sin\\alpha & \\cos\\alpha\n",
    "\\end{bmatrix} [i,j]^\\mathsf{T}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 编程实现\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一节中，我们将通过编程实现SIFT算法。将SIFT的每一步流程都编写成一个函数。首先，我们来实现尺度空间的生成。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-04T10:13:16.674225Z",
     "iopub.status.busy": "2022-09-04T10:13:16.673759Z",
     "iopub.status.idle": "2022-09-04T10:13:17.088035Z",
     "shell.execute_reply": "2022-09-04T10:13:17.086929Z",
     "shell.execute_reply.started": "2022-09-04T10:13:16.674130Z"
    }
   },
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from scipy.ndimage import maximum_filter\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "\n",
    "\n",
    "# 生成指定数目的高斯核，并对图像做高斯滤波\n",
    "def generateGimage(image, sigma, num_layers = 8, k_stride = 2):\n",
    "    sigma_res = np.sqrt(np.max([(sigma ** 2) - ((1) ** 2), 0.01]))\n",
    "    # (0,0)表示卷积核的大小根据sigma来决定\n",
    "    image = cv2.GaussianBlur(image, (0, 0), sigmaX=sigma_res, sigmaY=sigma_res)\n",
    "    \n",
    "    # 生成高斯核\n",
    "    k = 2 ** (1 / k_stride)\n",
    "    gaussian_kernels = np.zeros(num_layers)\n",
    "    # 第一层高斯就是1.6\n",
    "    gaussian_kernels[0] = sigma\n",
    "    \n",
    "    for i in range(1, num_layers):\n",
    "        # 根据高斯的性质，可以将大的高斯核拆分，减少计算量\n",
    "        gaussian_old = k**(i-1) * sigma\n",
    "        gaussian_new = k * gaussian_old\n",
    "        gaussian_kernels[i] = np.sqrt(gaussian_new**2 - gaussian_old**2)\n",
    "    \n",
    "    # 至此，我们已经得到一系列高斯核\n",
    "    \n",
    "    # 进行高斯模糊\n",
    "    gaussian_images = [image]\n",
    "    for kernel in gaussian_kernels:\n",
    "        tmp_image = cv2.GaussianBlur(image, (0, 0), sigmaX=kernel, sigmaY=kernel)\n",
    "        gaussian_images.append(tmp_image)\n",
    "    \n",
    "    # 返回不同高斯核对图像滤波的结果\n",
    "    return np.array(gaussian_images)\n",
    "        \n",
    "    \n",
    "# 生成DoG空间\n",
    "def generateDoGSpace(gaussian_images):\n",
    "    dog_images = []\n",
    "    for img1, img2 in zip(gaussian_images, gaussian_images[1:]):\n",
    "        dog_images.append(img2 - img1)\n",
    "    return dog_images\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在生成图像的尺度空间之后，需要检测图像的特征点，这包含以下两个步骤：1.判断每个位置其是否为局部极值点；2.若为极值点，则对其进行定位与筛选，并计算最终选出的特征点的特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.判别当前像素点是否为局部最大值点\n",
    "def isLocalExtremum(l1, l2, l3, threshold):\n",
    "    # l1，l2，l3是DoG尺度空间中相邻的三层，大小均已被切片成 3*3，即[l1,l2,l3]是一个cube\n",
    "    # 需要确定l2层的中心位置，即(1,1)位置是否为极值点\n",
    "    # threshold是一个设定的阈值，l2[1,1]必须大于阈值才可进行极值点的判定\n",
    "    if l2[1,1] > threshold:\n",
    "        if l2[1,1] > 0:\n",
    "            return np.all(l2[1,1]>=l1) and np.all(l2[1,1]>=l3) and np.sum(l2[1,1]<l2)==0\n",
    "        \n",
    "        elif l2[1,1] < 0:\n",
    "            return np.all(l2[1,1]<=l1) and np.all(l2[1,1]<=l3) and np.sum(l2[1,1]>l2)==0\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2.若为极值点，则对其进行定位与筛选，并计算最终选出的特征点的特征。\n",
    "# 在这一模块中，我们将编写实现计算特征点各个属性\n",
    "# 参考章节7.2.2以及7.2.3\n",
    "\n",
    "\n",
    "# 计算极值点特征\n",
    "def computeBlobAttribute(x, y, layer, dog_images, sigma, threshold, border, num_layers,\n",
    "                         g_image, corners, scales, orientations, layers):\n",
    "    # 计算(x,y)位置的特征点的特征\n",
    "    # 此时假设该位置的点为邻域内的最大值\n",
    "    # 首先需要对其进行筛选，根据牛顿迭代法判定是否为真正的极值点\n",
    "    # 筛除低对比度的特征点\n",
    "    # 消除边缘响应\n",
    "    # 在此之后便可得到特征点的信息\n",
    "    \n",
    "    gamma = 10\n",
    "    image_shape = dog_images[0].shape\n",
    "    out_flag = False\n",
    "    \n",
    "    # 使用牛顿迭代法得到极值点\n",
    "    # 至多更新5次，如果未收敛，则认为该候选点不是极值点\n",
    "    for iter_num in range(5):\n",
    "        # 得到DoG空间中相邻的三层\n",
    "        img1, img2, img3 = dog_images[layer-1:layer+2]\n",
    "        # 对这三层进行切片，得到大小为3*3的三层\n",
    "        cube = np.array([img1[x-1:x+2, y-1:y+2], img2[x-1:x+2, y-1:y+2],img3[x-1:x+2, y-1:y+2]])\n",
    "        \n",
    "        # 分别得到cube的一二阶导数\n",
    "        grad = compute1derivative(cube)\n",
    "        hessian = compute2derivative(cube)\n",
    "        \n",
    "        # 解方程得到牛顿迭代的更新值\n",
    "        update = -np.linalg.lstsq(hessian, grad, rcond=None)[0]\n",
    "        # 如果移动的距离太小，说明当前点里极值已收敛，直接返回当前点即可\n",
    "        if abs(update[0]) < 0.5 and abs(update[1]) < 0.5 and abs(update[2]) < 0.5:\n",
    "            break\n",
    "        # 更新当前点\n",
    "        y += int(round(update[0]))\n",
    "        x += int(round(update[1]))\n",
    "        layer += int(round(update[2]))\n",
    "        # 确保新的cube在DoG空间里\n",
    "        if x < border or x >= image_shape[0] - border or y < border or y >= image_shape[1] - border \\\n",
    "                                                      or layer < 1 or layer > num_layers - 2:\n",
    "            # 若不在空间中，则 out_flag = True\n",
    "            out_flag = True\n",
    "            break\n",
    "        \n",
    "    # 超出scale space或者未不收敛，直接return\n",
    "    if out_flag or iter_num >= 4:\n",
    "        return \n",
    "    \n",
    "    # 使用公式计算极值点的对比度\n",
    "    Extremum = cube[1, 1, 1] + 0.5 * np.dot(grad, update)\n",
    "    # 筛除低对比度的特征点\n",
    "    if np.abs(Extremum) >= threshold:\n",
    "        # 得到xy的hessian矩阵\n",
    "        xy_hessian = hessian[:2, :2]\n",
    "        xy_hessian_trace = np.trace(xy_hessian)\n",
    "        xy_hessian_det = np.linalg.det(xy_hessian)\n",
    "        # 消除边缘响应\n",
    "        if xy_hessian_det > 0 and (xy_hessian_trace ** 2) / xy_hessian_det < ((gamma + 1) ** 2) / gamma:\n",
    "            # 极值点坐标\n",
    "            pt = ((y + update[0]), (x + update[1]))\n",
    "            # 极值点尺度\n",
    "            size = sigma * (2 ** ((layer + update[2])))\n",
    "            # 计算特征点的方向（主方向和辅方向）\n",
    "            # 利用computeOrien函数进行计算\n",
    "            orien_list = computeOrien(pt, size, layer, g_image)\n",
    "            for tmp_orien in orien_list:\n",
    "                # 尺度\n",
    "                scales.append(size)\n",
    "                # 位置\n",
    "                layers.append(layer)\n",
    "                corners.append(pt)\n",
    "                # 方向\n",
    "                orientations.append(tmp_orien)\n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着对上述代码中出现的一些函数进行进一步的实现。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 编写一些数学函数\n",
    "\n",
    "# 计算一阶导数\n",
    "def compute1derivative(cube):\n",
    "    # cube是由DoG空间中相邻的三层组成，其中每层的大小被切片成 3*3\n",
    "    # 需要计算cube正中心位置的梯度\n",
    "    dx = (cube[1, 1, 2] - cube[1, 1, 0]) / 2\n",
    "    dy = (cube[1, 2, 1] - cube[1, 0, 1]) / 2\n",
    "    ds = (cube[2, 1, 1] - cube[0, 1, 1]) / 2\n",
    "    return np.array([dx, dy, ds])\n",
    "\n",
    "\n",
    "# 计算二阶导数\n",
    "def compute2derivative(cube):\n",
    "    # cube是由DoG空间中相邻的三层组成，其中每层的大小被切片成 3*3\n",
    "    # 需要计算cube正中心位置的梯度\n",
    "    \n",
    "    # 根据二阶导数的定义写出各个分量\n",
    "    center = cube[1, 1, 1]\n",
    "    \n",
    "    dxx = cube[1, 1, 2] + cube[1, 1, 0] - 2 * center\n",
    "    dyy = cube[1, 2, 1] + cube[1, 0, 1] - 2 * center\n",
    "    dss = cube[2, 1, 1] + cube[0, 1, 1] - 2 * center\n",
    "\n",
    "    dxy = (cube[1, 2, 2] - cube[1, 2, 0] - cube[1, 0, 2] + cube[1, 0, 0]) / 4\n",
    "    dxs = (cube[2, 1, 2] - cube[2, 1, 0] - cube[0, 1, 2] + cube[0, 1, 0]) / 4\n",
    "    dys = (cube[2, 2, 1] - cube[2, 0, 1] - cube[0, 2, 1] + cube[0, 0, 1]) / 4\n",
    "    \n",
    "    return np.array([[dxx, dxy, dxs], [dxy, dyy, dys],[dxs, dys, dss]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 计算特征点的方向\n",
    "def computeOrien(pt, size, layer, g_image):\n",
    "    # pt为特征点的位置\n",
    "    # size为对应的尺度大小\n",
    "    # layer为特征点所在层级\n",
    "    # g_image为高斯尺度空间的一层，即 g_images[layer]\n",
    "    \n",
    "    # 1.5倍的3sigma原则，决定圆半径的大小\n",
    "    radius = int(round(3 * size * 1.5))\n",
    "    image_shape = g_image.shape\n",
    "    # 设置直方图柱的数量\n",
    "    num_bins = 36\n",
    "    histogram = np.zeros(num_bins)\n",
    "    smooth_histogram = np.zeros(num_bins)\n",
    "    orien_list = []\n",
    "\n",
    "    # 遍历以特征点为中心的块状区域，块状区域的边长为2*radius\n",
    "    for i in range(-radius, radius + 1):\n",
    "        y = int(round(pt[1])) + i\n",
    "        # 判断坐标是否越界\n",
    "        if y > 0 and y < image_shape[0] - 1:\n",
    "            for j in range(-radius, radius + 1):\n",
    "                x = int(round(pt[0])) + j\n",
    "                # 判断坐标是否越界\n",
    "                if x > 0 and x < image_shape[1] - 1:\n",
    "                    # 计算当前位置的dx和dy\n",
    "                    dx = 0.5 * (g_image[y, x + 1] - g_image[y, x - 1])\n",
    "                    dy = 0.5 * (g_image[y + 1, x] - g_image[y - 1, x])\n",
    "                    # 计算当前位置梯度的值\n",
    "                    value = np.sqrt(dx * dx + dy * dy)\n",
    "                    # 计算当前位置梯度的方向\n",
    "                    orien = np.rad2deg(np.arctan2(dy, dx))\n",
    "                    # 高斯加权\n",
    "                    weight = np.exp(-0.5 / ((size*1.5) ** 2) * (i ** 2 + j ** 2))  \n",
    "                    histogram_index = int(round(orien * num_bins / 360.))\n",
    "                    histogram[histogram_index % num_bins] += weight * value\n",
    "\n",
    "    # 对直方图进行高斯平滑\n",
    "    for n in range(num_bins):\n",
    "        smooth_histogram[n] = (6 * histogram[n] + 4 * (histogram[n-1] + histogram[(n+1) % num_bins]) \n",
    "                               + histogram[n-2] + histogram[(n+2) % num_bins]) / 16.\n",
    "    # 选择主方向\n",
    "    orien_max = np.max(smooth_histogram)\n",
    "    orien_local_max = list(i for i in range(len(smooth_histogram)) \n",
    "                           if smooth_histogram[i] > smooth_histogram[i-1] and\n",
    "                           smooth_histogram[i] > smooth_histogram[(i+1)%num_bins])\n",
    "    \n",
    "    # 选择辅方向\n",
    "    for index in orien_local_max:\n",
    "        if smooth_histogram[index] >= 0.8 * orien_max:\n",
    "            orien_list.append(index * 360. / num_bins)\n",
    "    \n",
    "    return orien_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "至此，可以编写确定特征点的函数模块。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 对一张图像进行特征点检测\n",
    "def detect_blobs(image):\n",
    "    sigma = 1.6\n",
    "    num_layers = 4\n",
    "    border = 5\n",
    "    k_stride = 1\n",
    "    \n",
    "    # 生成高斯尺度空间\n",
    "    g_images = generateGimage(image, sigma, num_layers, k_stride)\n",
    "    # 生成DoG空间\n",
    "    dog_images = generateDoGSpace(g_images)\n",
    "        \n",
    "    # 开始寻找blob\n",
    "    threshold = 0.02\n",
    "    corners = []\n",
    "    scales = []\n",
    "    orientations = []\n",
    "    layers = []\n",
    "    \n",
    "    for layer, (image1, image2, image3) in enumerate(zip(dog_images, dog_images[1:], dog_images[2:])):\n",
    "        # 忽略太靠近边缘的点\n",
    "        for x in range(border, image1.shape[0]-border):\n",
    "            for y in range(border, image2.shape[1]-border):\n",
    "                # 检测当前位置是否为局部极值\n",
    "                if isLocalExtremum(image1[x-1:x+2, y-1:y+2], image2[x-1:x+2, y-1:y+2], \n",
    "                                   image3[x-1:x+2, y-1:y+2], threshold):\n",
    "                    # 如果是候选点，则进行进一步定位筛选，并返回其信息\n",
    "                    computeBlobAttribute(x, y, layer+1, dog_images, sigma, threshold, border, \n",
    "                                         num_layers, g_images[layer], corners, scales, orientations, layers)\n",
    "                    \n",
    "    return g_images, corners, scales, orientations, layers  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，对检测到的特征点进行描述符的生成。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIFT特征点描述符的生成\n",
    "def compute_descriptors(g_images, corners, scales, orientations, layers):\n",
    " \n",
    "    if len(corners) != len(scales) or len(corners) != len(orientations):\n",
    "        raise ValueError(\n",
    "          '`corners`, `scales` and `orientations` must all have the same length.')\n",
    "    \n",
    "    descriptors_list = []\n",
    "    \n",
    "    for pt, size, orien, layer in zip(corners, scales, orientations, layers):\n",
    "        \n",
    "        # 读取特征点的各项信息\n",
    "        g_image = g_images[layer]\n",
    "        x, y = np.round(np.array(pt)).astype(np.int32)\n",
    "        orien = 360 - orien\n",
    "        \n",
    "        # 计算块状区域大小\n",
    "        win_s = 3 * size\n",
    "        win_l = int(round(min(2**0.5 * win_s * (4+1) / 2, np.sqrt(g_image.shape[0]**2+g_image.shape[1]**2))))\n",
    "        \n",
    "        # 用list依次存储块状区域内所有点的信息\n",
    "        i_index = []\n",
    "        j_index = []\n",
    "        value_list = []\n",
    "        orien_index = []\n",
    "        \n",
    "        # 三维数组存储16个窗口的8个方向，为防止计算时边界溢出，在行、列的首尾各扩展一次\n",
    "        result_cube = np.zeros((4 + 2, 4 + 2, 8)) \n",
    "        \n",
    "        # 统计一个块状区域内16个子区域的像素梯度直方图\n",
    "        for i in range(-win_l, win_l+1):\n",
    "            for j in range(-win_l, win_l+1):\n",
    "                \n",
    "                # 获得旋转之后的坐标\n",
    "                i_rotate = j * np.sin(np.deg2rad(orien)) + i * np.cos(np.deg2rad(orien))\n",
    "                j_rotate = j * np.cos(np.deg2rad(orien)) - i * np.sin(np.deg2rad(orien))\n",
    "                \n",
    "                # 计算4*4子区域对应的下标\n",
    "                tmp_i = (i_rotate / win_s) + 2 - 0.5\n",
    "                tmp_j = (j_rotate / win_s) + 2 - 0.5\n",
    "                \n",
    "                # 邻域的点在旋转后，仍然处于4×4的区域内\n",
    "                if tmp_i > -1 and tmp_j > -1 and tmp_i < 4 and tmp_j < 4:\n",
    "                    \n",
    "                    # 该特征点在原图像中的位置\n",
    "                    i_inimg = int(round(y + i))\n",
    "                    j_inimg = int(round(x + j))\n",
    "                    \n",
    "                    if i_inimg > 0 and j_inimg > 0 and i_inimg < g_image.shape[0]-1 \\\n",
    "                                                   and j_inimg < g_image.shape[1]-1:\n",
    "                        # 计算梯度大小\n",
    "                        dx = g_image[i_inimg, j_inimg + 1] - g_image[i_inimg, j_inimg - 1]\n",
    "                        dy = g_image[i_inimg - 1, j_inimg] - g_image[i_inimg + 1, j_inimg]\n",
    "                        grad_value = np.sqrt(dx**2 + dy**2)\n",
    "                        # 计算梯度方向\n",
    "                        grad_orien = np.rad2deg(np.arctan2(dy, dx)) % 360\n",
    "                        \n",
    "                        i_index.append(tmp_i)\n",
    "                        j_index.append(tmp_j)\n",
    "                        # 进行高斯加权\n",
    "                        g_weight = np.exp(-1 / 8 * ((i_rotate / win_s) ** 2 + (j_rotate / win_s) ** 2))\n",
    "                        value_list.append(g_weight * grad_value)\n",
    "                        # 将梯度方向投影到8个方向，这里的grad_orien是原图像中的梯度，需要叠加上旋转之后的角度\n",
    "                        orien_index.append((grad_orien - orien) * 8 / 360)\n",
    "        \n",
    "        # 将每个方向的幅值插入进矩阵result_cube中\n",
    "        for i, j, value, orien1 in zip(i_index, j_index, value_list, orien_index):\n",
    "            tirlinearInterpolation(i, j, value, orien1, result_cube)\n",
    "            \n",
    "        descriptor = result_cube[1:-1, 1:-1, :].flatten()  \n",
    "        # 计算描述符的大小\n",
    "        l2norm = np.linalg.norm(descriptor)\n",
    "        # 设定阈值\n",
    "        threshold = l2norm * 0.2\n",
    "        # 将描述符进行截断，将大于阈值的值设定为阈值\n",
    "        descriptor[descriptor > threshold] = threshold\n",
    "        # 归一化，确保描述符具有尺度不变性\n",
    "        descriptor /= l2norm\n",
    "        # 添加描述符\n",
    "        descriptors_list.append(descriptor)\n",
    "        \n",
    "    return descriptors_list\n",
    "\n",
    "\n",
    "# 将每个方向的幅值插入进矩阵result_cube中\n",
    "def tirlinearInterpolation(i, j, value, orien, result_cube):\n",
    "    \n",
    "    # 由于位置坐标(i,j,orien)是一个浮点数，没有办法直接在result_cube中赋值\n",
    "    # 考虑到该浮点位置会对临近的8个整数坐标都有一定的贡献（可以将该浮点坐标想象成一个正方体中的一点，它对正方体的8个顶点都有贡献）\n",
    "    # 我们可以根据三线性插值法，根据该点到各顶点的距离为相邻的8个顶点赋值\n",
    "    # 感兴趣的同学可以自行学习三线性插值法\n",
    "    \n",
    "    # 对坐标先进性量化，将其转换为整数\n",
    "    i_quant = int(np.floor(i))\n",
    "    j_quant = int(np.floor(j))\n",
    "    orien_quant = int(np.floor(orien)) % 8\n",
    "    \n",
    "    # 计算量化前后的偏移量\n",
    "    i_residual = i - i_quant\n",
    "    j_residual =  j - j_quant\n",
    "    orien_residual = (orien - orien_quant) % 8\n",
    "    \n",
    "    # 根据三线性插值法写出当前位置对每个顶点的权重\n",
    "    c0 = (1 - i_residual) * value\n",
    "    c1 = i_residual * value\n",
    "    c11 = c1 * j_residual\n",
    "    c10 = c1 * (1 - j_residual)\n",
    "    c01 = c0 * j_residual\n",
    "    c00 = c0 * (1 - j_residual)\n",
    "    \n",
    "    c111 = c11 * orien_residual\n",
    "    c110 = c11 * (1 - orien_residual)\n",
    "    c101 = c10 * orien_residual\n",
    "    c100 = c10 * (1 - orien_residual)\n",
    "    c011 = c01 * orien_residual\n",
    "    c010 = c01 * (1 - orien_residual)\n",
    "    c001 = c00 * orien_residual\n",
    "    c000 = c00 * (1 - orien_residual)\n",
    "    \n",
    "    # 进行赋值操作\n",
    "    result_cube[i_quant + 1, j_quant + 1, orien_quant] += c000\n",
    "    result_cube[i_quant + 1, j_quant + 1, (orien_quant + 1) % 8] += c001\n",
    "    result_cube[i_quant + 1, j_quant + 2, orien_quant] += c010\n",
    "    result_cube[i_quant + 1, j_quant + 2, (orien_quant + 1) % 8] += c011\n",
    "    result_cube[i_quant + 2, j_quant + 1, orien_quant] += c100\n",
    "    result_cube[i_quant + 2, j_quant + 1, (orien_quant + 1) % 8] += c101\n",
    "    result_cube[i_quant + 2, j_quant + 2, orien_quant] += c110\n",
    "    result_cube[i_quant + 2, j_quant + 2, (orien_quant + 1) % 8] += c111\n",
    "    \n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "至此，我们已经完成SIFT特征的编写。先导入两张图像，寻找他们的特征点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 展示图像中的特征点\n",
    "def draw(img, corners):\n",
    "    # corners为特征点的坐标\n",
    "    img1 = img.copy()\n",
    "    for m in corners:\n",
    "        pt = (int(m[0]), int(m[1]))\n",
    "        cv2.circle(img1, pt, 1, (0,255,0), 2)\n",
    "    return img1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-04T10:13:17.091358Z",
     "iopub.status.busy": "2022-09-04T10:13:17.090633Z",
     "iopub.status.idle": "2022-09-04T10:13:26.846228Z",
     "shell.execute_reply": "2022-09-04T10:13:26.845160Z",
     "shell.execute_reply.started": "2022-09-04T10:13:17.091306Z"
    }
   },
   "outputs": [],
   "source": [
    "# 读取图像\n",
    "img1 = cv2.imread('test1.png', cv2.IMREAD_COLOR)\n",
    "img2 = cv2.imread('test2.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",
    "# 获得两张图像特征点的各项信息\n",
    "g_images1, corners1, scales1, orientations1, layers1 = detect_blobs(gray1)\n",
    "g_images2, corners2, scales2, orientations2, layers2 = detect_blobs(gray2)\n",
    "\n",
    "# 获得两张图像的SIFT特征点\n",
    "descriptors1 = compute_descriptors(g_images1, corners1, scales1, orientations1, layers1)\n",
    "descriptors2 = compute_descriptors(g_images2, corners2, scales2, orientations2, layers2)\n",
    "\n",
    "img1_detect = draw(img1, corners1)\n",
    "img2_detect = draw(img2, corners2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-04T10:13:26.847674Z",
     "iopub.status.busy": "2022-09-04T10:13:26.847344Z",
     "iopub.status.idle": "2022-09-04T10:13:27.273221Z",
     "shell.execute_reply": "2022-09-04T10:13:27.272223Z",
     "shell.execute_reply.started": "2022-09-04T10:13:26.847643Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 273192 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 225196 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 展示两张图像的特征点\n",
    "plt.title('Output1')\n",
    "plt.axis('off')\n",
    "plt.imshow(img1_detect[:, :, ::-1])\n",
    "plt.show()\n",
    "\n",
    "plt.title('Output2')\n",
    "plt.axis('off')\n",
    "plt.imshow(img2_detect[:, :, ::-1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用SIFT算法很容易就能发现图像中的特征点。接下来，需要将寻找出的特征点进行匹配。将两幅图像中的特征点集合分别定义为$\\mathcal{S}$，$\\mathcal{R}$。接着，对于$\\mathcal{S}$中的某一个特征点$\\mathbf{s}_i$，利用$\\ell_2$范数去度量其对应的特征向量$\\mathbf{Vs}_{i}$和集合$\\mathcal{R}$中每个特征点$\\mathbf{r}_j$对应的特征向量$\\mathbf{Vr}_{j}$之间的距离，并确定距离$\\mathbf{s}_i$最近的特征点$\\mathbf{r}_k$和第二近的特征点$\\mathbf{r}_d$。如果\n",
    "\n",
    "$$\\frac{\\ell_2(\\mathbf{s}_i, \\mathbf{r}_k)}{\\ell_2(\\mathbf{s}_i, \\mathbf{r}_d)} < \\theta,$$\n",
    "\n",
    "其中$\\theta$是设定的阈值，便可以认为$\\mathbf{s}_i$与$\\mathbf{r}_k$相匹配。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 特征点匹配\n",
    "def match_descriptors(descriptors1, descriptors2):\n",
    "\n",
    "    max_index = np.zeros((len(descriptors1))) - 1\n",
    "    # 初始化第一和第二近的特征点对数组\n",
    "    maxmatch = np.zeros((len(descriptors1))) + 1e10\n",
    "    secmatch = np.zeros((len(descriptors1))) + 1e10\n",
    "    # 设定阈值\n",
    "    threshold = 0.8\n",
    "\n",
    "    for vec1_index in range(len(descriptors1)):\n",
    "        for vec2_index in range(len(descriptors2)):\n",
    "            # 计算特征点对的距离\n",
    "            distance = np.linalg.norm(descriptors1[vec1_index] - descriptors2[vec2_index])\n",
    "            if distance < maxmatch[vec1_index]:\n",
    "                # 更新当前最匹配的特征点\n",
    "                maxmatch[vec1_index] = distance\n",
    "                # 记录当前特征点对的距离\n",
    "                max_index[vec1_index] = vec2_index\n",
    "                \n",
    "            elif distance < secmatch[vec1_index]:\n",
    "                # 更新第二近特征点对的距离\n",
    "                # 这里只需要更新距离即可，因为不关心哪一个点是第二近的\n",
    "                secmatch[vec1_index] = distance\n",
    "\n",
    "    matches = []\n",
    "    # 返回匹配的特征点对的坐标信息\n",
    "    for i in range(len(descriptors1)):\n",
    "        if maxmatch[i] / secmatch[i] < threshold:\n",
    "            matches.append((i, int(max_index[i])))\n",
    "\n",
    "    return matches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-04T10:13:27.276188Z",
     "iopub.status.busy": "2022-09-04T10:13:27.275695Z",
     "iopub.status.idle": "2022-09-04T10:13:27.294246Z",
     "shell.execute_reply": "2022-09-04T10:13:27.293190Z",
     "shell.execute_reply.started": "2022-09-04T10:13:27.276141Z"
    }
   },
   "outputs": [],
   "source": [
    "def draw_matches(image1, image2, corners1, corners2, matches):\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",
    "            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",
    "            cv2.line(match_image, pt1, pt2, (0, 0, 255))\n",
    "        \n",
    "    return match_image\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-04T10:13:27.296527Z",
     "iopub.status.busy": "2022-09-04T10:13:27.295792Z",
     "iopub.status.idle": "2022-09-04T10:13:27.445171Z",
     "shell.execute_reply": "2022-09-04T10:13:27.444111Z",
     "shell.execute_reply.started": "2022-09-04T10:13:27.296438Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 471.5, 171.5, -0.5)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 111492 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 进行图像匹配\n",
    "matches = match_descriptors(descriptors1, descriptors2)\n",
    "image1 = cv2.imread('test1.png', cv2.IMREAD_GRAYSCALE)\n",
    "image2 = cv2.imread('test2.png', cv2.IMREAD_GRAYSCALE)\n",
    "\n",
    "# 绘制匹配连线\n",
    "match_image = draw_matches(image1, image2, corners1, corners2, matches)\n",
    "\n",
    "plt.title('Match')\n",
    "plt.imshow(match_image)\n",
    "plt.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以很明显的看出，两幅图像中相同的特征点完美的匹配在了一起。再这之后，便可以进行图像拼接等一系列后续操作了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4 小结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一章中，我们学习了SIFT算法的基本原理，并动手实现了SIFT算法。SIFT特征具有多种优点，它稳定性高，对旋转、尺度缩放、亮度变化等能够保持不变；独特性好，信息量丰富，所以SIFT特征可以作为许多计算机视觉任务的基础。在后面的章节中，我们将介绍SIFT的一个实际应用———图像拼接。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.5 习题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 习题 7.1: 请解释SIFT特征点提取的基本原理和步骤，并说明如何选择适当的尺度空间和关键点。\n",
    "\n",
    "#### 习题 7.2: SIFT特征检测算法如何实现尺度不变性？解释尺度不变性的原理和实现方法。\n",
    "\n",
    "#### 习题 7.3: SIFT特征检测算法的主要优点是什么？与其他特征检测算法（如Harris角点检测）相比有什么不同之处？\n",
    "\n",
    "#### 习题 7.4: SIFT特征检测算法对于图像的光照变化是否敏感？如果是，请解释原因并提出一种解决光照变化的方法；如果不是，请解释原因。\n",
    "\n",
    "#### 习题 7.5：SIFT特征点描述时，为什么要把特征点所在的块状区域的横轴旋转到主方向上？这样做的原理是什么？\n",
    "\n",
    "#### 习题 7.6：在SIET中，我们使用“最近点对”和“次近点对“之间的比值来挑选匹配两幅图像中的特征点，以选择最相似的特征点对。然而，当某种纹理在一张图像中反复出现 (如图 7-2 中的块状区域)时，便会给特征点对的匹配带来极大的干扰。你能设计出更加鲁棒的方法来匹配特征点对吗?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.6 参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] David G. Lowe: Object Recognition from Local Scale-Invariant Features. ICCV 1999: 1150-1157"
   ]
  }
 ],
 "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.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
