{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第7章 特征检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还记得上一章中我们学习过的角点吗？其实角点是图像中“特征点”的一种。所谓特征点，即能够表示图像特性或者局部特性的像素，这些像素叫称为特征点。通常，特征点是由关键点和描述子组成的。关键点是图像中含有“关键信息”的像素，这些点是一些十分突出的点，不会因光照条件的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点，也包括像素的位置、角度等信息。描述子是用来唯一描述关键点的一串数字，又称为描述符，与每个人的个人信息相似。今天，就让我们一起来学习如何检测特征点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2 SIFT算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SIFT的全称是 Scale Invariant Feature Transform（尺度不变特征变换），是由加拿大教授David G.Lowe在1999年在会议文章中提出，2004年发表在IJCV上， 是计算机视觉界近二十年来引用率最高的文章之一。SIFT的特征提取方面对计算机视觉近年来的发展影响深远，特别是几乎影响到了后续所有的角点提取和匹配算法。SIFT算法的主要步骤可包括为：\n",
    "\n",
    "1. 尺度空间极值检测：搜索所有尺度上的图像位置，通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。\n",
    "\n",
    "2. 关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。\n",
    "\n",
    "3. 方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。\n",
    "\n",
    "4. 关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化\n",
    "\n",
    "接下来，我们将详细介绍每一个步骤。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.1 尺度空间极值检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.1 尺度空间理论"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "尺度空间理论是指，在图像信息处理模型中引入一个被视为尺度的参数，通过连续改变尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等。\n",
    "\n",
    "尺度空间方法将传统的单尺度图像信息处理技术纳入尺度不断变化的动态分析框架中，更容易获取图像的本质特征。尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。\n",
    "\n",
    "尺度空间满足视觉不变性。该不变性的视觉解释如下：当我们用眼睛观察物体时，一方面当物体所处背景的光照条件变化时，视网膜感知图像的亮度水平和对比度是不同的，因此要求尺度空间算子对图像的分析不受图像的灰度水平和对比度变化的影响，即满足灰度不变性和对比度不变性；另一方面，相对于某一固定坐标系，当观察者和物体之间的相对位置变化时，视网膜所感知的图像的位置、大小、角度和形状是不同的，因此要求尺度空间算子对图像的分析和图像的位置、大小、角度以及仿射变换无关，即满足平移不变性、尺度不变性、欧几里德不变性以及仿射不变性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.2 尺度空间的表示"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一个图像的尺度空间$\\mathbf{L}(i,j,\\sigma)$定义为一个变化尺度的高斯函数$\\mathbf{G}(i,j,\\sigma)$与原图像$\\mathbf{I}(i,j)$的卷积，即：\n",
    "\n",
    "$$\\mathbf{L}(i,j,\\sigma) = \\mathbf{G}(i,j,\\sigma) * \\mathbf{I}(i,j)$$\n",
    "\n",
    "其中，$(i,j)$表示图像的像素位置，$\\sigma$是尺度空间因子，根据我们在滤波中学习的知识，我们知道$\\sigma$值越小，图像被平滑的越少，即相应的尺度就越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.3 高斯金字塔的构建"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "尺度空间可以使用高斯金字塔表示，高斯金字塔的构建分为两部分：\n",
    "\n",
    "1. 对图像做不同尺度的高斯模糊；\n",
    "2. 对图像做降采样(隔点采样)。\n",
    "\n",
    "图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：\n",
    "\n",
    "$$n=\\log_2\\{\\min(M,N)\\} - t, t \\in [0, \\log_2\\{\\min(M,N)\\})$$\n",
    "\n",
    "其中$M,N$为原图像的大小，$t$为塔顶图像的最小维数的对数值。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"1.png\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-1 高斯金字塔</div>\n",
    "</center>\n",
    "\n",
    "为了让尺度体现其连续性，高斯金字塔在简单降采样的基础上加上了高斯滤波。如图 7-1 所示，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，每组含有多张图像。另外，降采样时，高斯金字塔上一组图像的初始图像(底层图像)是由前一组图像的倒数第三张图像隔点采样得到的。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.4 高斯差分金字塔"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们知道尺度归一化的高斯拉普拉斯函数为$\\sigma^2 \\nabla^2\\mathbf{G}$\n",
    "\n",
    "高斯差分函数（Difference of Gaussian，简称DOG算子）$\\mathbf{D}(i,j,\\sigma)$，我们用差分近似代替微分，因为：\n",
    "\n",
    "$$\\frac{\\partial \\mathbf{G}}{\\partial \\sigma} = \\sigma \\nabla^2\\mathbf{G}$$\n",
    "\n",
    "所以，\n",
    "\n",
    "$$\\sigma \\nabla^2\\mathbf{G} = \\frac{\\partial \\mathbf{G}}{\\partial \\sigma} \\approx \\frac{\\mathbf{G}(i,j,k\\sigma) - \\mathbf{G}(i,j,\\sigma)}{k\\sigma - \\sigma} $$\n",
    "\n",
    "因此，我们可以得到：\n",
    "\n",
    "$$\\mathbf{D}(i,j,\\sigma)=\\mathbf{G}(i,j,k\\sigma) - \\mathbf{G}(i,j,\\sigma) \\approx {(k-1)\\sigma^2}\\nabla^2\\mathbf{G}$$\n",
    "\n",
    "其中k-1是个常数，并不影响极值点位置的求取。如图 7-2 所示，我们可以发现高斯拉普拉斯和高斯差分近似相同。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"2.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-2 高斯拉普拉斯和高斯差分的比较</div>\n",
    "</center>\n",
    "\n",
    "我们使用更高效的高斯差分算子代替拉普拉斯算子进行极值检测：\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf{D}(i,j,\\sigma) & = (\\mathbf{G}(i,j,k\\sigma) - \\mathbf{G}(i,j,\\sigma))* \\mathbf{I}(i,j) \\\\\n",
    "& = \\mathbf{L}(i,j,k\\sigma)-\\mathbf{L}(i,j,\\sigma)\\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "在实际计算时，我们使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如图 7-3 所示。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"3.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-3 高斯差分金字塔</div>\n",
    "</center>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.5 空间极值点检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关键点是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DOG相邻两层图像之间比较完成的。为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图3.4所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。 \n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"4.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-4 DOG空间极值检测</div>\n",
    "</center>\n",
    "\n",
    "由于要在相邻尺度进行比较，如图 7-3 右侧每组含4层的高斯差分金子塔，只能在中间两层中进行两个尺度的极值点检测，其它尺度则只能在不同组中进行。为了在每组中检测$S$个尺度的极值点，则DOG金字塔每组需$S+2$层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔每组需$S+3$层图像，实际计算时$S$在3到5之间。\n",
    "\n",
    "注意，局部极值点都是在同一个组当中进行的，所以肯定有这样的问题，某一组当中的第一个图和最后一个图层没有前一张图和下一张图，那该怎么计算？ 解决办法是，再用高斯模糊，在高斯金字塔多“模糊”出三张来凑数，所以在DOG中多出两张（同一组中第一层和最后一层不做极值点检测）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.1.6 构建尺度空间需要的参数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们定义$\\sigma(o,s)$为尺度空间坐标，其中$o$为组的索引，$s$为组内层的索引。在构建高斯金字塔时，组内每层的尺度坐标按如下公式计算：\n",
    "\n",
    "$$\\sigma(s) = \\sqrt{(k^s\\sigma_0)^2-(k^{s-1}\\sigma_0)^2}$$\n",
    "\n",
    "其中，$\\sigma_0$为初始尺度，可以取1.6，$s$为组内的层索引，不同组相同层的组内尺度坐标$\\sigma(s)$相同。组内下一层图像是由前一层图像按$\\sigma(s)$进行高斯模糊所得；$k$为组内总层数的倒数，即$k=2^{\\frac{1}{\\mathbf{S}}}$。由此，我们可以将高斯差分算子写成：\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf{D}(i,j,\\sigma) & = (\\mathbf{G}(i,j,\\sigma(s+1)) - \\mathbf{G}(i,j,\\sigma(s)))* \\mathbf{I}(i,j) \\\\\n",
    "& = \\mathbf{L}(i,j,\\sigma(s+1))-\\mathbf{L}(i,j,\\sigma(s))\\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.2 关键点定位"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以上方法检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DOG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.2.1 关键点的精确定位"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "离散空间的极值点并不是真正的极值点，图 7-5 显示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation）。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"5.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-5 离散空间与连续空间极值点的差别 </div>\n",
    "</center>\n",
    "\n",
    "为了提高关键点的稳定性，需要对尺度空间DOG函数进行曲线拟合。利用DOG函数在尺度空间的泰勒展开为：\n",
    "\n",
    "$$\\mathbf{D}(\\mathbf{X}) = \\mathbf{D}(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",
    "对应极值点，方程的值为：\n",
    "\n",
    "$$\\mathbf{D}(\\hat{\\mathbf{X}}) = \\mathbf{D}(0) +\\frac{1}{2}\\frac{\\partial\\mathbf{D}^\\mathsf{T}}{\\partial \\mathbf{X}}\\hat{\\mathbf{X}}$$\n",
    "\n",
    "其中，当$\\hat{\\mathbf{X}}$在任一维度上的偏移量大于0.5时（即$x,y$或$\\sigma$），这意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置，同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，另外，过小的点易受噪声的干扰而变得不稳定，所以将小于某个经验值（如0.03）的极值点删除。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2.2.2 消除边缘响应"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DOG算子会产生较强的边缘响应，因此我们需要剔除不稳定的边缘响应点。曲率可以通过一个2x2的Hessian矩阵$\\mathbf{H}$求出:\n",
    "\n",
    "$$\\mathbf{H} = \\begin{bmatrix}\n",
    "\\mathbf{D}_{xx} & \\mathbf{D}_{xy} \\\\\n",
    "\\mathbf{D}_{xy} &\\mathbf{D}_{yy}  \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "$\\mathbf{H}$的特征值为$\\alpha, \\beta$，${\\rm tr}(\\mathbf{H}) = \\mathbf{D}_{xx}+\\mathbf{D}_{yy}=\\alpha+\\beta$表示$\\mathbf{H}$对角线元素之和，${\\rm det}(\\mathbf{H}) = \\mathbf{D}_{xx}\\mathbf{D}_{yy} - (\\mathbf{D}_{xy})^2 = \\alpha\\beta$表示$\\mathbf{H}$的行列式。假设$\\alpha$是较大的特征值，令$\\alpha= r \\beta$，则：\n",
    "\n",
    "$$\\frac{{\\rm tr}(\\mathbf{H})^2}{{\\rm det}(\\mathbf{H})} = \\frac{(\\alpha+\\beta)^2}{\\alpha\\beta} = \\frac{(r+1)^2}{r}$$\n",
    "\n",
    "$\\frac{(r+1)^2}{r}$的值在两个特征值相等时最小，随着$r$的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测曲率是否在某域值r下，只需检测\n",
    "\n",
    "$$\\frac{{\\rm tr}(\\mathbf{H})^2}{{\\rm det}(\\mathbf{H})} <\\frac{(r+1)^2}{r}$$\n",
    "\n",
    "当上式成立时，保留关键点，反之剔除。r可以取值为10（参照原论文）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.3 关键点方向确定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了使描述符具有旋转不变性，我们需要利用图像的局部特征为给每一个关键点分配一个基准方向，可以使用图像梯度的方法求取局部结构的稳定方向。对于在DOG金字塔中检测出的关键点，采集其所在高斯金字塔图像$3\\sigma$邻域窗口内像素的梯度和方向分布特征。梯度的模值和方向如下：\n",
    "\n",
    "$$m(i,j)=\\sqrt{(\\mathbf{L}(i+1,j) - \\mathbf{L}(i-1,j))^2 +(\\mathbf{L}(i,j+1) - \\mathbf{L}(i,j-1))^2}$$\n",
    "$$\\theta(i,j) = \\tan^{-1}(\\frac{\\mathbf{L}(i,j+1)-\\mathbf{L}(i,j-1)}{\\mathbf{L}(i+1,j)-\\mathbf{L}(i-1,j)})$$\n",
    "\n",
    "在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度。如图 7-6 所示，直方图的峰值方向代表了关键点的主方向（为简化，图中只画了八个方向的直方图）。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"6.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-6 关键点方向直方图 </div>\n",
    "</center>\n",
    "\n",
    "方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。虽然仅有15％的关键点被赋予多个方向，但这可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值。\n",
    "\n",
    "至此，将检测出的含有位置、尺度和方向的关键点即是该图像的SIFT特征点。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.4 关键点描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。 \n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"7.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-7 把关键点周围区域旋转到主方向上 </div>\n",
    "</center>\n",
    "\n",
    "对每个关键点，让该关键点周围区域旋转到主方向上（如图 7-7 ），然后我们取以关键点为圆心半径为r的圆的外接正方形区域，将其分成4×4个子区域，其中子区域以像素为单位，如图 7-8 所示。在每个子区域内进行统计每个像素点在8个方向上的梯度幅值（8个方向为东，西，南，北，东南，东北，西南，西北）。梯度幅值可以经过高斯加权，也可以不加权。每个区域用长度为8的向量表示，将16个子区域的向量合并构成的长度128一维的向量就是该关键点的描述符。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"8.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图7-8 把关键点周围区域旋转到主方向上 </div>\n",
    "</center>\n",
    "\n",
    "至此，SIFT特征描述向量生成。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 编程实现\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在一节中，我们将通过编程实现SIFT算法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from scipy.ndimage.filters import maximum_filter\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "\n",
    "def generateGimage(image, sigma, num_layers = 8, k_stride = 2):\n",
    "    '''\n",
    "    先生成指定数量的高斯核，在对image进行高斯模糊，得到scale space\n",
    "    '''\n",
    "    sigma_res = np.sqrt(np.max([(sigma ** 2) - ((1) ** 2), 0.01]))\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",
    "    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",
    "    return np.array(gaussian_images)\n",
    "        \n",
    "def generateDoGSpace(gaussian_images):\n",
    "    \"\"\"\n",
    "    生成dog空间\n",
    "    \"\"\"\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",
    "\n",
    "def isLocalExtremum(l1, l2, l3, threshold):\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",
    "        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\n",
    "\n",
    "def compute1derivative(cube):\n",
    "    \"\"\"\n",
    "    计算得到一阶导数\n",
    "    \"\"\"\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",
    "def compute2derivative(cube):\n",
    "    \"\"\"\n",
    "    计算二阶导数\n",
    "    \"\"\"\n",
    "    center = cube[1, 1, 1]\n",
    "    dxx = cube[1, 1, 2] - 2 * center + cube[1, 1, 0]\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",
    "    dyy = cube[1, 2, 1] - 2 * center + cube[1, 0, 1]\n",
    "    dys = (cube[2, 2, 1] - cube[2, 0, 1] - cube[0, 2, 1] + cube[0, 0, 1]) / 4\n",
    "    dss = cube[2, 1, 1] - 2 * center + cube[0, 1, 1]\n",
    "    return np.array([[dxx, dxy, dxs], [dxy, dyy, dys],[dxs, dys, dss]])\n",
    "\n",
    "def computeOrien(pt, size, layer, g_image):\n",
    "    # 1.5倍的3sigma原则，决定圆的大小\n",
    "    radius = int(round(3 * size * 1.5))\n",
    "    image_shape = g_image.shape\n",
    "    num_bins = 36\n",
    "    histogram = np.zeros(num_bins)\n",
    "    smooth_histogram = np.zeros(num_bins)\n",
    "    orien_list = []\n",
    "\n",
    "    for i in range(-radius, radius + 1):\n",
    "        y = int(round(pt[1])) + i\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",
    "                if x > 0 and x < image_shape[1] - 1:\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",
    "                    value = np.sqrt(dx * dx + dy * dy)\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",
    "    for n in range(num_bins):\n",
    "        #高斯平滑\n",
    "        smooth_histogram[n] = (6 * histogram[n] + 4 * (histogram[n-1] + histogram[(n+1) % num_bins]) + 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)) if smooth_histogram[i] > smooth_histogram[i-1] and smooth_histogram[i] > smooth_histogram[(i+1)%num_bins])\n",
    "    # 大于最大值80%的，也要视为一个独立的blob\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\n",
    "\n",
    "def computeBlobAttribute(x, y, layer, dog_images, sigma, threshold, border, num_layers, g_image, corners, scales, orientations, layers):\n",
    "    gamma = 10\n",
    "    image_shape = dog_images[0].shape\n",
    "    out_flag = False\n",
    "    \n",
    "    '''\n",
    "    使用牛顿迭代法得到极值点\n",
    "    至多更新5次，如果未收敛，则认为该候选点不是极值点\n",
    "    我们选择的极值点的坐标不需要量化，应当尽量避免量化误差\n",
    "    '''\n",
    "    for iter_num in range(5):\n",
    "        img1, img2, img3 = dog_images[layer-1:layer+2]\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",
    "        # 分别得到一二阶导数\n",
    "        grad = compute1derivative(cube)\n",
    "        hessian = compute2derivative(cube)\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在space里\n",
    "        if x < border or x >= image_shape[0] - border or y < border or y >= image_shape[1] - border or layer < 1 or layer > num_layers - 2:\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",
    "        # 特征值都为正，且最大最小特征值的比值不超过10，此时该点一般不为边缘\n",
    "        if xy_hessian_det > 0 and (xy_hessian_trace ** 2) / xy_hessian_det < ((gamma + 1) ** 2) / gamma:\n",
    "            pt = ((y + update[0]), (x + update[1]))\n",
    "            size = sigma * (2 ** ((layer + update[2])))\n",
    "            # 计算keypoint的方向，由于一个blob可能根据方向变成多个blob，所以返回值是list\n",
    "            orien_list = computeOrien(pt, size, layer, g_image)\n",
    "            for tmp_orien in orien_list:\n",
    "                corners.append(pt)\n",
    "                scales.append(size)\n",
    "                layers.append(layer)\n",
    "                orientations.append(tmp_orien)\n",
    "    return \n",
    "\n",
    "def tirlinearInterpolation(i, j, value, orien, result_cube):\n",
    "    i_quant = int(np.floor(i))\n",
    "    j_quant = int(np.floor(j))\n",
    "    orien_quant = int(np.floor(orien)) % 8\n",
    "                \n",
    "    i_residual = i - i_quant\n",
    "    j_residual =  j - j_quant\n",
    "    orien_residual = (orien - orien_quant) % 8\n",
    "    \n",
    "    c1 = i_residual * value\n",
    "    c0 = (1 - 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",
    "    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",
    "    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 \n",
    "\n",
    "def detect_blobs(image):\n",
    "    \"\"\"Laplacian blob detector.\n",
    "    \n",
    "    Args:\n",
    "    - image (2D float64 array): A grayscale image.\n",
    "    \n",
    "    Returns:\n",
    "    - corners (list of 2-tuples): A list of 2-tuples representing the locations\n",
    "        of detected blobs. Each tuple contains the (x, y) coordinates of a\n",
    "        pixel, which can be indexed by image[y, x].\n",
    "    - scales (list of floats): A list of floats representing the scales of\n",
    "        detected blobs. Has the same length as `corners`.\n",
    "    - orientations (list of floats): A list of floats representing the dominant\n",
    "        orientation of the blobs.\n",
    "    \"\"\"\n",
    "    sigma = 1.6\n",
    "    num_layers = 4\n",
    "    border = 5\n",
    "    k_stride = 1\n",
    "    g_images = generateGimage(image, sigma, num_layers, k_stride)\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], image3[x-1:x+2, y-1:y+2], threshold):\n",
    "                    # 如果是备选点，就计算其各项属性\n",
    "                    computeBlobAttribute(x, y, layer+1, dog_images, sigma, threshold, border, num_layers, g_images[layer], corners, scales, orientations, layers)\n",
    "                    \n",
    "    return g_images, corners, scales, orientations, layers  \n",
    "\n",
    "\n",
    "def compute_descriptors(g_images, corners, scales, orientations, layers):\n",
    "    \"\"\"Compute descriptors for corners at specified scales.\n",
    "    \n",
    "    Args:\n",
    "    - image (2d float64 array): A grayscale image.\n",
    "    - corners (list of 2-tuples): A list of (x, y) coordinates.\n",
    "    - scales (list of floats): A list of scales corresponding to the corners.\n",
    "        Must have the same length as `corners`.\n",
    "    - orientations (list of floats): A list of floats representing the dominant\n",
    "        orientation of the blobs.\n",
    "    \n",
    "    Returns:\n",
    "    - descriptors (list of 1d array): A list of desciptors for each corner.\n",
    "        Each element is an 1d array of length 128.\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",
    "        g_image = g_images[layer]\n",
    "        # 读取blob的各项信息\n",
    "        x, y = np.round(np.array(pt)).astype(np.int32)\n",
    "        orien = 360 - orien\n",
    "        '''\n",
    "        根据论文计算各个窗口的大小\n",
    "        小窗口边长为1.5sigma\n",
    "        '''\n",
    "        win_s = 1.5 * 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",
    "        # 用list依次存储方框内所有点的信息\n",
    "        i_index = []\n",
    "        j_index = []\n",
    "        value_list = []\n",
    "        orien_index = []\n",
    "        # 三维数组存储16个窗口的8个方向，多出来的2是作为border，用来计算三线性插值的\n",
    "        result_cube = np.zeros((4 + 2, 4 + 2, 8)) \n",
    "        \n",
    "        for i in range(-win_l, win_l+1):\n",
    "            for j in range(-win_l, win_l+1):\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",
    "                tmp_i = (i_rotate / win_s) + 2 - 0.5\n",
    "                tmp_j = (j_rotate / win_s) + 2 - 0.5\n",
    "\n",
    "                if tmp_i > -1 and tmp_j > -1 and tmp_i < 4 and tmp_j < 4:\n",
    "                    i_inimg = int(round(y + i))\n",
    "                    j_inimg = int(round(x + j))\n",
    "                    if i_inimg > 0 and j_inimg > 0 and i_inimg < g_image.shape[0]-1 and j_inimg < g_image.shape[1]-1:\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",
    "                        grad_orien = np.rad2deg(np.arctan2(dy, dx)) % 360\n",
    "                        i_index.append(tmp_i)\n",
    "                        j_index.append(tmp_j)\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",
    "                        orien_index.append((grad_orien - orien) * 8 / 360)\n",
    "        \n",
    "        # 进行三线性插值\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",
    "        # 归一化与截断\n",
    "        descriptor = result_cube[1:-1, 1:-1, :].flatten()  \n",
    "        l2norm = np.linalg.norm(descriptor)\n",
    "        threshold = l2norm * 0.2\n",
    "        descriptor[descriptor > threshold] = threshold\n",
    "        descriptor /= l2norm\n",
    "        descriptors_list.append(descriptor)\n",
    "    return descriptors_list\n",
    "\n",
    "\n",
    "def draw(img, 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\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread('graf.png', cv2.IMREAD_COLOR)\n",
    "gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) / 255.0\n",
    "g_images, corners, scales, orientations, layers = detect_blobs(gray)\n",
    "descriptors = compute_descriptors(g_images, corners, scales, orientations, layers)\n",
    "img1 = draw(img, corners)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 204716 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 204040 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Original Input')\n",
    "plt.axis('off')\n",
    "plt.imshow(img[:, :, ::-1])\n",
    "plt.show()\n",
    "\n",
    "plt.title('Output')\n",
    "plt.axis('off')\n",
    "plt.imshow(img1[:, :, ::-1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一章中，我们学习了SIFT算法的具体原理，并动手实现了SIFT。SIFT具有多种优点，它稳定性高，对旋转、尺度缩放、亮度变化等能够保持不变；独特性好，信息量丰富。但是，SIFT仍然存在一些缺点，比如实时性不高等。接下来，我们将介绍SIFT的一个实际应用————图像拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
