{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第6章 角点检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果一个点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点，换言之，图像的梯度中的局部最大值所对应的像素点就是角点。在现实世界中，角点对应于物体的拐角，道路的十字路口、丁字路口等。角点是图像中的一个重要概念，在这一章中，我们将介绍如何对一张图像中的角点进行检测。\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;\">图6-1 角点样例</div>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2 Harris角点检测算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在正式介绍角点检测算法之前，我们先来学习一些预备知识。图 6-2 是图像中可能存在的三种不同的情况。为了更清楚的解释彼此的不同，我们在图像中使用一个小的滑动窗口。当滑动窗处于平坦区域时，沿着任何方向均没有明显变化，即滑动窗在横轴和纵轴方向梯度为0；当滑动窗处于边缘时，沿着边缘方向的梯度为0，在图 6-2 中即纵轴方向的梯度为0；当滑动窗处于角点位置时，无论哪一个方向的梯度都会发生明显变化。\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;\">图6-2 滑动窗可能位于的三个位置</div>\n",
    "</center>\n",
    "\n",
    "因此，根据上述讨论，我们可以按照以下方法来检测角点：使用一个固定大小的窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后像素灰度变化程度，如果存在任意方向上的滑动都有着较大灰度变化，那么我们可以认为该窗口中存在角点。其中，Harris算法便是一个代表性的角点检测算法。\n",
    "\n",
    "Harris算法主要分为三步：\n",
    "\n",
    "1. 当窗口同时向$X$和$Y$两个方向移动时，计算窗口内部的像素值变化量；\n",
    "2. 对于每个窗口，计算其对应的一个角点响应函数；\n",
    "3. 对该函数进行阈值处理，提取角点。\n",
    "\n",
    "我们将按照步骤分别详细介绍。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2.1 像素值变化量计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对一张图像$\\mathbf{I}$，考虑一个窗函数$\\mathbf{W}$沿着$(1,1)$方向移动$(u,v)$个单位，我们定义像素值变化量$E(u,v)$如下：\n",
    "\n",
    "$$E(u,v) = \\sum_{i, j}\\mathbf{W}(i,j)[\\mathbf{I}(i+u, j+v) - \\mathbf{I}(i,j)]^2$$\n",
    "\n",
    "其中，$\\mathbf{I}(i,j)$是当前像素的强度，$\\mathbf{I}(i+u, j+v)$是位移后像素的强度，$\\mathbf{W}(i,j)$表示窗函数，其中窗函数可以是平坦的，也可以是高斯的，是一个二维的滤波器，如图 6-3 所示。\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;\">图6-3 窗函数</div>\n",
    "</center>\n",
    "\n",
    "对于一个角点来说，$E(u,v)$会非常大。因此，我们可以最大化$E(u,v)$来检测图像中的角点。不过，使用上述公式计算会非常缓慢，因此，我们使用上述公式的一阶泰勒展开以获取其近似形式。\n",
    "\n",
    "对一张二维图像而言，其一阶泰勒展开可写作：\n",
    "\n",
    "$$\\mathbf{I}(i+u,j+v) \\approx \\mathbf{I}(i,j) + u\\mathbf{I}_x + v\\mathbf{I}_y$$\n",
    "\n",
    "其中，$\\mathbf{I}_x$是图像沿$X$轴的梯度，$\\mathbf{I}_y$是图像沿$Y$轴的梯度，即：\n",
    "\n",
    "$$\\mathbf{I}_x = \\frac{\\partial \\mathbf{I}(i,j)}{\\partial x}, \\mathbf{I}_y = \\frac{\\partial \\mathbf{I}(i,j)}{\\partial y}$$\n",
    "\n",
    "因此，我们可以将$E(u,v)$写成：\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "E(u,v)  &= \\sum_{i, j}\\mathbf{W}(i,j)[\\mathbf{I}(i, j) +u\\mathbf{I}_x + v\\mathbf{I}_y- \\mathbf{I}(i,j)]^2 \\\\\n",
    "& = \\sum_{i, j}\\mathbf{W}(i,j)(u\\mathbf{I}_x + v\\mathbf{I})^2 \\\\\n",
    "& = \\sum_{i, j}\\mathbf{W}(i,j)(u^2\\mathbf{I}_x^2+v^2\\mathbf{I}_y^2+2uv\\mathbf{I}_x\\mathbf{I}_y)\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "将$(u,v)$提取出，我们便得到最终的形式：\n",
    "\n",
    "$$E(u, v) = \\begin{bmatrix}u,v \\end{bmatrix}\\mathbf{M}\\begin{bmatrix}u\\\\v \\end{bmatrix}$$\n",
    "\n",
    "其中，矩阵$\\mathbf{M}$为：\n",
    "\n",
    "$$\\mathbf{M} = \\sum_{i,j}\\mathbf{W}(i,j)\\begin{bmatrix}\\mathbf{I}_x^2 & \\mathbf{I}_x\\mathbf{I}_y \\\\ \\mathbf{I}_x\\mathbf{I}_y & \\mathbf{I}_y^2\\end{bmatrix}$$\n",
    "\n",
    "不难发现$\\mathbf{M}$是一个对称矩阵。我们将$\\mathbf{M}$进一步处理为：\n",
    "\n",
    "$$\\mathbf{M} = \\mathbf{R}^{-1} \\begin{bmatrix} \\lambda_1 & 0 \\\\ 0 & \\lambda_2 \\end{bmatrix} \\mathbf{R} $$\n",
    "\n",
    "我们可以将$\\mathbf{R}$看作为旋转因子。经对角化处理后，我们可以将两个正交方向的变化分量提取出，即$\\lambda_1, \\lambda_2$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2.2 计算角点响应函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们计算$E(u,v)$的目的是要找到会引起像素值发生较大变化的窗口。根据6.2.1中推导的公式我们不难发现，像素值变化的大小取决于矩阵$\\mathbf{M}$。\n",
    "\n",
    "我们知道$\\mathbf{M} = \\mathbf{R}^{-1} \\begin{bmatrix} \\lambda_1 & 0 \\\\ 0 & \\lambda_2 \\end{bmatrix} \\mathbf{R} $。实际上，$\\mathbf{M}$可以想象成一个椭圆，其长度由特征值决定，方向由$\\mathbf{R}$决定，如图 6-4 所示。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"4.png\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图6-4 矩阵M示意</div>\n",
    "</center>\n",
    "\n",
    "椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系图 6-5 所示。共可分为三种情况：\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"5.png\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图6-5 椭圆特征值与图像的关系</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "1. 一个特征值大，另一个特征值小，对应图像中的边缘；\n",
    "2. 两个特征值都小，且近似相等，对应图像中的平面；\n",
    "3. 两个特征值都大，且近似相等，对应图像中的角点\n",
    "\n",
    "因此，我们可以通过计算$\\mathbf{M}$的两个特征值$\\lambda_1$和$\\lambda_2$的对图像角点进行检测。我们知道\n",
    "\n",
    "$${\\rm det}(\\mathbf{M}) = \\lambda_1\\lambda_2$$\n",
    "$${\\rm trace}(\\mathbf{M}) = \\lambda_1+\\lambda_2$$\n",
    "\n",
    "我们定义角点响应函数$\\theta$：\n",
    "\n",
    "$$\\theta = {\\rm det}(\\mathbf{M}) - \\alpha {\\rm trace}(\\mathbf{M})^2 = \\lambda_1\\lambda_2 - \\alpha(\\lambda_1+\\lambda_2)^2$$\n",
    "\n",
    "这里，我们计算矩阵的行列式和迹来代替计算矩阵的特征值。$\\alpha$是一个常数，通常取0.04-0.06。我们便可以将图 6-5 重制为图 6-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;\">图6-6 角点响应函数与图像的关系</div>\n",
    "</center>\n",
    "\n",
    "我们可以用角点响应函数来绝地个一个窗口是平面、边缘还是角点：\n",
    "\n",
    "1. 平面：该窗口在平坦区域上滑动，窗口内的灰度值基本不会发生变化，所以|$\\theta$|值非常小；\n",
    "2. 边缘：$\\theta$为负数，仅在水平或竖直方向有较大的变化量；\n",
    "3. 角点：$\\theta$很大，在水平、竖直两个方向上变化均较大的点。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2.3 提取角点"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以根据$\\theta$的值将该窗口窗口区域划分为平面、边缘或者角点。为了得到最优的角点，我们还可以使用非极大抑制。Harris角点检测的结果是一个矩阵，每一个位置表示该像素点对应的$\\theta$值。我们可以设定一个阈值，当$\\theta$的值大于该阈值时，我们便可以认为该窗口区域为角点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.3 算法实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们已经学习了Harris角点的相关知识，接下来，我们将编程实现Harris角点检测。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import *\n",
    "\n",
    "\n",
    "def responseFunc(M):\n",
    "    k = 0.04\n",
    "    det = np.linalg.det(M)\n",
    "    trace = np.trace(M)\n",
    "    R = det - k * trace **2\n",
    "    return R\n",
    "\n",
    "\n",
    "def harris_corners(src, NMS=False):\n",
    "    threshold = 0.01\n",
    "    h, w = src.shape[:2]\n",
    "    \n",
    "    # 将图像转为灰度图像\n",
    "    gray_image = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)\n",
    "    cornerPoint = np.zeros_like(gray_image,dtype=np.float32)\n",
    "    \n",
    "    # 计算图像沿横轴和纵轴的梯度\n",
    "    grad = np.zeros((h, w, 2), dtype=np.float32)\n",
    "    grad[:,:,0] = cv2.Sobel(gray_image,cv2.CV_16S,1,0)\n",
    "    grad[:,:,1] = cv2.Sobel(gray_image,cv2.CV_16S,0,1)\n",
    "\n",
    "    Ixx = grad[:,:,0] ** 2\n",
    "    Iyy = grad[:,:,1] ** 2\n",
    "    Ixy = grad[:,:,0] * grad[:,:,1]\n",
    "    \n",
    "    # 窗函数使用高斯核\n",
    "    Ixx = cv2.GaussianBlur(Ixx, (3, 3), sigmaX=2)\n",
    "    Iyy = cv2.GaussianBlur(Iyy, (3, 3), sigmaX=2)\n",
    "    Ixy = cv2.GaussianBlur(Ixy, (3, 3), sigmaX=2)\n",
    "\n",
    "    \n",
    "    for i in range(gray_image.shape[0]):\n",
    "        for j in range(gray_image.shape[1]):\n",
    "            struture_matrix = [[Ixx[i][j], Ixy[i][j]], [Ixy[i][j], Iyy[i][j]]]\n",
    "\n",
    "            # 计算角点响应函数\n",
    "            R = responseFunc(struture_matrix)\n",
    "            cornerPoint[i][j] = R\n",
    "\n",
    "    # 非最大抑制\n",
    "    corners = np.zeros_like(gray_image,dtype=np.float32)\n",
    "    counter = 0\n",
    "    maxValue = np.max(cornerPoint)\n",
    "    for i in range(cornerPoint.shape[0]):\n",
    "        for j in range(cornerPoint.shape[1]):\n",
    "            if NMS:\n",
    "                if cornerPoint[i][j] > threshold * maxValue and cornerPoint[i][j] == np.max(cornerPoint[max(0, i - 1):min(i + 1, h - 1), max(0, j - 1):min(j + 1, w - 1)] ):\n",
    "                    counter+=1\n",
    "                    corners[i][j] = 255\n",
    "            else:\n",
    "                if cornerPoint[i][j] > threshold * maxValue:\n",
    "                    counter+=1\n",
    "                    corners[i][j] = 255\n",
    "                    \n",
    "    # print('FeaturePoints=', counter)\n",
    "    return corners\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 84908 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 84532 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "img = cv2.imread('sudoku.png')\n",
    "\n",
    "plt.imshow(img[:, :, ::-1])\n",
    "plt.axis('off')\n",
    "plt.title('Original Input')\n",
    "plt.show()\n",
    "\n",
    "response = harris_corners(img)\n",
    "img[response == 255] = [0, 0, 255]\n",
    "plt.imshow(img[:, :, ::-1])\n",
    "plt.axis('off')\n",
    "plt.title('Harris Corner Response')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.4 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本章中，我们学习了角点检测的相关知识以及Harris角点检测算法。Harris角点检测是计算机视觉中常用的一个算法，随着人们研究的深入，Harris算法也在不断地被改进，如Shi-Tomasi算法等，感兴趣的同学可以自行查阅相关资料。接下来的章节中，我们将介绍计算机视觉中另一种常用的技术————特征匹配与检测。"
   ]
  }
 ],
 "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
}
