{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第5章 边缘检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "边缘（Edge）是图像中局部亮度、纹理或颜色发生剧烈变化的位置。如图5-1所示，这些变化往往发生在物体与背景、物体与物体、区域与区域之间的边界上，边缘蕴含了图像中物体的形状信息和场景的结构信息，这些信息在视觉系统中扮演着至关重要的作用。印第安纳大学心理学教授Linda B. Smith提出形状信息相比其他信息在人眼识别物体中更为重要。例如图5-1中的右边的边缘图已经完全丢失了色彩信息，但是我们仍然可以识别出这些边缘所表示的语义——天鹅。David Marr的视觉理论也印证了这一点，根据他的理论，视觉信息的处理第一步是建立在从图像中提取的边缘等重要特征上的，也就是说边缘检测是视觉信息处理的基础。在本章中，我们将介绍边缘检测的数学模型以及几种经典边缘检测算子的实现。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-abff6d8eed86fa6999a0f5b95a287ea2_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图5-1 边缘检测示意</div>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 边缘检测的数学基础"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "简介中我们给出了边缘定性的定义，但是本节给出检测边缘需要定量的定义。给定一张数字图像$\\mathbf{I}$，边缘检测目标是计算出图像中每一点$[i,j]$处的边缘强度$\\mathbf{E}_s[i,j|I]$和边缘方向$\\mathbf{E}_{\\theta}[i,j|I]$。这就涉及到边缘的数学模型。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-0da47e95c674306a9f48bfcd5026e304_1440w.jpg\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图5-2 一维信号中的边缘</div>\n",
    "</center>\n",
    "\n",
    "图 5-2 显示了一个一维信号，大家可以把它看作二维图像像素强度沿一个轴向（$X$轴）的信号。根据简介中的定义，边缘所处的位置就是像素强度发生快速变化的位置。根据微积分的知识，我们知道函数的变化可以用函数的一阶导数来表示。也就是说，图像边缘就是图像一阶导数的极大值处。接下来，我们介绍图像一阶导数（微分）的计算。\n",
    "\n",
    "根据微积分的知识，一维函数$f(x)$的一阶微分定义为：\n",
    "\n",
    "$$\\frac{df(x)}{dx} = \\lim_{\\Delta x \\rightarrow 0} \\frac{f(x+\\Delta x)-f(x-\\Delta x)}{2\\Delta x}$$\n",
    "\n",
    "而我们知道，图像$\\mathbf{I}$是一个二维函数$\\mathbf{I}=f(x,y)$，所以我们将一维的一阶微分推广到二维：\n",
    "\n",
    "$$\\mathbf{I}_x=\\frac{\\partial f(x,y)}{\\partial x} = \\lim_{\\epsilon \\rightarrow 0} \\frac{f(x+\\epsilon, y)-f(x-\\epsilon,y)}{2\\epsilon}$$\n",
    "\n",
    "$$\\mathbf{I}_y=\\frac{\\partial f(x,y)}{\\partial y} = \\lim_{\\epsilon \\rightarrow 0} \\frac{f(x,y+\\epsilon)-f(x,y-\\epsilon)}{2\\epsilon}$$\n",
    "\n",
    "上式中的$[\\mathbf{I}_x,\\mathbf{I}_y]$通常被称为图像梯度（gradient），记为$\\nabla\\mathbf{I}$，即图像像素强度沿图像$X$轴和$Y$轴的一阶偏导数。由于数字图像是离散的，在数字图像中最小的位置距离单位为1个像素，即$\\epsilon=1$。那么在离散情况下，上式变为：\n",
    "\n",
    "$$ \\mathbf{I}_x[i,j] = \\frac{\\mathbf{I}[i+1, j]-\\mathbf{I}[i-1,j]}{2}$$\n",
    "\n",
    "$$ \\mathbf{I}_y[i,j] = \\frac{\\mathbf{I}[i,j+1]-\\mathbf{I}[i,j-1]}{2}$$\n",
    "\n",
    "所以，在离散的数字图像里，对图像的一阶微分是通过一阶差分实现的。我们将上式进行变形：\n",
    "\n",
    "$$ \\mathbf{I}_x[i,j] =\\frac{1}{2}\\sum^{1}_{m=-1}\\sum^{1}_{n=-1}\\mathbf{I}[i-m,j-n]\\mathbf{K}_x[m,n]$$\n",
    "\n",
    "$$ \\mathbf{I}_y[i,j]=\\frac{1}{2}\\sum^{1}_{m=-1}\\sum^{1}_{n=-1}\\mathbf{I}[i-m,j-n]\\mathbf{K}_y[m,n]$$\n",
    "\n",
    "其中，$\\mathbf{K}_x=\\begin{bmatrix}\n",
    "   0 & 0 & 0 \\\\\n",
    "   +1 & 0 & -1 \\\\\n",
    "   0 & 0 & 0\n",
    "  \\end{bmatrix}$和$\\mathbf{K}_y=\\begin{bmatrix}\n",
    "   0 & +1 & 0 \\\\\n",
    "   0 & 0 & 0 \\\\\n",
    "   0 & -1 & 0\n",
    "  \\end{bmatrix}$。这里我们令矩阵$\\mathbf{K}_x$中心位置的坐标为$[0,0]$，对$\\mathbf{K}_y[0,0]$同理。大家由上式的形式联想到什么没有？没错，上式正是卷积的形式，说明我们可以通过卷积计算图像的一阶差分，即\n",
    "  \n",
    "  $$\\mathbf{I}_x=\\mathbf{I}*\\mathbf{K}_x,~\\mathbf{I}_y=\\mathbf{I}*\\mathbf{K}_y$$\n",
    "  \n",
    "  从而得到图像每个像素点的梯度$\\nabla\\mathbf{I}[i,j]=[\\mathbf{I}_x[i,j],\\mathbf{I}_y[i,j]]$。图像$\\mathbf{I}$位置$[i,j]$处的边缘强度$\\mathbf{E}_s[i,j|\\mathbf{I}]$和边缘方向$\\mathbf{E}_{\\theta}[i,j|\\mathbf{I}]$可分别根据图像梯度$\\nabla\\mathbf{I}[i,j]$计算得到：\n",
    "  \n",
    "  $$\\mathbf{E}_s[i,j|\\mathbf{I}]=\\|\\nabla\\mathbf{I}[i,j]\\|=\\sqrt{\\mathbf{I}_x^2[i,j]+\\mathbf{I}_y^2[i,j]}$$\n",
    "  \n",
    "  $$\\mathbf{E}_{\\theta}[i,j|\\mathbf{I}]=\\tan^{-1}(\\frac{\\mathbf{I}_y[i,j]}{\\mathbf{I}_x[i,j]})+\\frac{\\pi}{2}$$\n",
    "\n",
    "\n",
    "注意，边缘方向和图像梯度方向相互垂直。由此可见，计算图像边缘的关键是计算图像的梯度，而图像梯度的是通过差分计算得到。在上面的计算中，图像差分通过以$\\mathbf{K}_x$和$\\mathbf{K}_y$为卷积核对图像进行卷积实现，这两个卷积核被称为差分算子，也是最简单的差分算子。边缘检测的主要工作之一就是设计更合理的差分算子从而更准确更鲁棒地计算图像梯度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3 边缘检测算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "边缘检测算法通常有三个步骤：\n",
    "\n",
    "1. 图像平滑：边缘检测算法主要是基于图像像素强度的微分，但微分对图像噪声很敏感，因为图像噪声也会导致图像像素强度地突变。因此我们首先对图像进行滤波平滑，以此来提高边缘检测算法对图像噪声的鲁棒性。\n",
    "\n",
    "2. 图像梯度计算：通过差分模板对平滑后的图像进行卷积，得到图像梯度；\n",
    "\n",
    "3. 图像边缘定位：根据梯度得到的边缘强度和方向定位图像中的边缘。\n",
    "\n",
    "由于图像平滑和差分运算都可以通过卷积完成，根据卷积的结合律，上述前两个步骤也可以合成一步完成。接下来，我们将介绍两种的边缘检测算法：Sobel边缘检测算法和Canny边缘检测算法。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3.1 Sobel边缘检测算法 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sobel边缘检测算法 [1] 是最经典的边缘检测算法之一，其核心是Sobel算子，该算子为两组$3{\\times}3$的卷积核，数学形式如下：\n",
    "\n",
    "$$\\mathbf{K}_x = \\begin{bmatrix}\n",
    "   +1 & 0 & -1 \\\\\n",
    "   +2 & 0 & -2 \\\\\n",
    "   +1 & 0 & -1\n",
    "  \\end{bmatrix},\\quad\n",
    "  \\mathbf{K}_y = \\begin{bmatrix}\n",
    "   +1 & +2 & +1 \\\\\n",
    "   0 & 0 & 0 \\\\\n",
    "   -1 & -2 & -1\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "我们选取沿$X$轴方向的Sobel算子，将其拆分为两个向量相乘：\n",
    "\n",
    "$$\n",
    "\\mathbf{K}_x = \\begin{bmatrix}\n",
    "   +1 & 0 & -1 \\\\\n",
    "   +2 & 0 & -2 \\\\\n",
    "   +1 & 0 & -1\n",
    "  \\end{bmatrix} = \\begin{bmatrix} \n",
    "  1 \\\\ 2 \\\\ 1\n",
    "  \\end{bmatrix}\\begin{bmatrix} \n",
    "  +1 & 0 & -1\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "前者类似于高斯滤波器（中间大两边小），目的是对图像中的噪声进行抑制；根据上一章节的介绍，后者是一个差分算子，作用于图像得到沿$X$轴方向的梯度，也就是垂直方向上的边缘。\n",
    "\n",
    "Sobel算子包含图像平滑与图像差分操作，一步完成上述前两个步骤。接下来，我们编程实现Sobel算子与图像的卷积，分别计算沿$X$轴和$Y$轴的图像梯度$\\mathbf{I}_x$和$\\mathbf{I}_y$。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:17:52.525249Z",
     "iopub.status.busy": "2022-09-20T09:17:52.524296Z",
     "iopub.status.idle": "2022-09-20T09:17:54.158634Z",
     "shell.execute_reply": "2022-09-20T09:17:54.157641Z",
     "shell.execute_reply.started": "2022-09-20T09:17:52.525122Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 788.5, 430.5, -0.5)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 295300 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# 输入图像\n",
    "img = cv2.imread('lena.jpeg')\n",
    "\n",
    "# x方向算子\n",
    "kx = np.array([\n",
    "    [1,0,-1],\n",
    "    [2,0,-2],\n",
    "    [1,0,-1]\n",
    "])\n",
    "\n",
    "# y方向Sobel算子\n",
    "ky = np.array([\n",
    "    [1, 2, 1],\n",
    "    [0, 0, 0],\n",
    "    [-1,-2,-1]\n",
    "])\n",
    "\n",
    "# 展示输入图像\n",
    "plt.imshow(img[:, :, ::-1])\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们接下来先检测$X$轴方向的梯度。利用已经学习过的卷积的知识，我们可以很轻松的编写出相应的代码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:17:54.161398Z",
     "iopub.status.busy": "2022-09-20T09:17:54.160807Z",
     "iopub.status.idle": "2022-09-20T09:17:54.593141Z",
     "shell.execute_reply": "2022-09-20T09:17:54.592392Z",
     "shell.execute_reply.started": "2022-09-20T09:17:54.161353Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 788.5, 430.5, -0.5)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 327396 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# x方向Sobel算子与图像进行卷积\n",
    "conv_x = cv2.filter2D(img, -1, kx)\n",
    "plt.imshow(conv_x[:, :, ::-1])\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以观察到，沿$X$轴方向的梯度$\\mathbf{I}_x$，也就是垂直方向上的边缘信息被有效检测出，如手臂的线条、帽子等。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:17:54.595076Z",
     "iopub.status.busy": "2022-09-20T09:17:54.594599Z",
     "iopub.status.idle": "2022-09-20T09:17:54.995872Z",
     "shell.execute_reply": "2022-09-20T09:17:54.995166Z",
     "shell.execute_reply.started": "2022-09-20T09:17:54.595039Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 788.5, 430.5, -0.5)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 324136 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# y方向Sobel算子与图像进行卷积\n",
    "conv_y = cv2.filter2D(img, -1, ky)\n",
    "plt.imshow(conv_y[:, :, ::-1])\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如上图所示，沿$Y$轴方向的梯度$\\mathbf{I}_y$，也就是水平方向的边缘信息也被很好的提取出，如眼睛的轮廓等。基于图像梯度$[\\mathbf{I}_x,\\mathbf{I}_y]$，代入章节5.2中的边缘强度和边缘方向公式，即可得到边缘检测的结果。\n",
    "\n",
    "在计算边缘强度时，为了提高效率我们通常使用不开平方的近似值：\n",
    "\n",
    "$$\\mathbf{E}_s[i,j|\\mathbf{I}]={\\|\\mathbf{I}_x[i,j] \\| + \\| \\mathbf{I}_y[i,j] \\|}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:17:54.997687Z",
     "iopub.status.busy": "2022-09-20T09:17:54.997235Z",
     "iopub.status.idle": "2022-09-20T09:17:55.202498Z",
     "shell.execute_reply": "2022-09-20T09:17:55.201514Z",
     "shell.execute_reply.started": "2022-09-20T09:17:54.997653Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 788.5, 430.5, -0.5)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 409364 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E = abs(conv_x) + abs(conv_y)\n",
    "plt.imshow(E[:, :, ::-1])\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3.2 Canny算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "虽然Sobel算子中加入了类似于高斯滤波的算子，使其对图像噪声有一定的抑制性，但是区分图像噪声和边缘本身是一个很困难的任务，特别是图像中常常会存在由杂乱纹理产生的虚假边缘。通过图像平滑能够一定程度上抑制噪声，但是，平滑也会使得图像变得模糊，降低边缘的强度，所以对图像的平滑程度应当是可以调节的，边缘检测的结果应该也是多尺度的（例如可以调节平滑程度）。接下来我们介绍一种多尺度边缘检测的算法：Canny边缘检测算法 [2]。\n",
    "\n",
    "Canny算法可能是计算机视觉中最广泛使用的边缘检测算法，它通常通过以下几步加以实现：\n",
    "\n",
    "1. 图像平滑：即抑制噪声；\n",
    "2. 图像梯度计算；\n",
    "3. 边缘强度非最大值抑制（non-maximum suppression，nms）；\n",
    "4. 基于滞后阈值化（hysteresis-based thresholding）的边缘定位。\n",
    "\n",
    "我们将按照其步骤分别进行详细介绍。\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3.2.1 图像平滑"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于边缘检测对图像噪声非常的敏感，因此在进行后续的操作之前，我们需要先对图像进行去噪。结合第3章中所学的知识，我们一般使用高斯滤波器对图像进行去噪。对图像平滑的程度取决于高斯核标准差$\\sigma$的大小。$\\sigma$越大，图像越平滑，但也越模糊；$\\sigma$越小，图像细节保留更多，但是噪声也越多。不同的$\\sigma$代表了不同的尺度，调节$\\sigma$的大小，在不同平滑程度的图像上进行边缘检测，得到的便是多尺度边缘检测结果。这里，我们先将高斯核的大小设置为5，$\\sigma$设置为1.5（一般高斯核的大小是$\\sigma$的3倍时滤波效果比较好）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:33:09.639798Z",
     "iopub.status.busy": "2022-09-20T09:33:09.638828Z",
     "iopub.status.idle": "2022-09-20T09:33:09.955227Z",
     "shell.execute_reply": "2022-09-20T09:33:09.954099Z",
     "shell.execute_reply.started": "2022-09-20T09:33:09.639752Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 251640 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 158428 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from utils import *\n",
    "# 设定高斯核的大小与方差\n",
    "kernel_size = 5\n",
    "sigma = 1.5\n",
    "\n",
    "# 导入一张图片\n",
    "img = Image.open('chong.png')\n",
    "\n",
    "# 将其转为灰度图像\n",
    "img = img.convert('L')\n",
    "img = np.array(img)\n",
    "\n",
    "# 构造高斯卷积核\n",
    "kernel = gaussian_kernel(kernel_size, sigma)\n",
    "\n",
    "# 与图像进行卷积\n",
    "smoothed = cv2.filter2D(img, -1, kernel)\n",
    "\n",
    "# 展示图像\n",
    "plot_image(img, 'Original image')\n",
    "plot_image(smoothed, 'Smoothed image')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3.2.2 图像梯度计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当我们获得了平滑后的图像之后，我们可以利用类似Sobel算子的方式来计算图像的水平方向梯度$\\mathbf{I}_x$和竖直方向梯度$\\mathbf{I}_y$；接着，我们便可以计算图像中每个像素梯度的大小$\\mathbf{E}_s[i,j|\\mathbf{I}]=\\|\\nabla\\mathbf{I}[i,j]\\|=\\sqrt{\\mathbf{I}_x^2[i,j]+\\mathbf{I}_y^2[i,j]}$，其中梯度方向是 $\\tan^{-1}(\\frac{\\mathbf{I}_y[i,j]}{\\mathbf{I}_x[i,j]})$，即垂直于边缘的方向，如图 5-3 所示。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic4.zhimg.com/80/v2-fe8bfe27580cf88ccaa26f17c4e9d503_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图5-3 梯度的方向</div>\n",
    "</center>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:34:23.260175Z",
     "iopub.status.busy": "2022-09-20T09:34:23.259063Z",
     "iopub.status.idle": "2022-09-20T09:34:23.564352Z",
     "shell.execute_reply": "2022-09-20T09:34:23.563493Z",
     "shell.execute_reply.started": "2022-09-20T09:34:23.260113Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 142552 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 164028 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 求x方向梯度\n",
    "def partial_x(img):\n",
    "    \n",
    "    # 获得图像的大小\n",
    "    Hi, Wi = img.shape  \n",
    "    out = np.zeros((Hi, Wi))\n",
    "    \n",
    "    # 这里我们将卷积核均值化\n",
    "    k = np.array([[0,0,0],[0.5,0,-0.5],[0,0,0]])\n",
    "    \n",
    "    # 对图像进行卷积\n",
    "    out = cv2.filter2D(img, -1, k)\n",
    "    \n",
    "    return out\n",
    "\n",
    "\n",
    "# 求y方向梯度\n",
    "def partial_y(img):\n",
    "    Hi, Wi = img.shape  \n",
    "    out = np.zeros((Hi, Wi))\n",
    "    k = np.array([[0,0.5,0],[0,0,0],[0,-0.5,0]])\n",
    "    out = cv2.filter2D(img, -1, k)\n",
    "    return out\n",
    "\n",
    "\n",
    "# 获得在两个方向上图像的梯度\n",
    "Gx = partial_x(smoothed)\n",
    "Gy = partial_y(smoothed)\n",
    "\n",
    "# 绘制图像梯度\n",
    "plot_image(Gx, 'Derivative in x direction')\n",
    "plot_image(Gy, 'Derivative in y direction')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:34:43.403784Z",
     "iopub.status.busy": "2022-09-20T09:34:43.403451Z",
     "iopub.status.idle": "2022-09-20T09:34:43.576826Z",
     "shell.execute_reply": "2022-09-20T09:34:43.575893Z",
     "shell.execute_reply.started": "2022-09-20T09:34:43.403745Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 270884 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 计算梯度以及方向\n",
    "def gradient(img):\n",
    "    G = np.zeros(img.shape)\n",
    "    theta = np.zeros(img.shape)\n",
    "    dx = partial_x(img)\n",
    "    dy = partial_y(img)\n",
    "    \n",
    "    # 获得图像梯度\n",
    "    G = np.sqrt(dx**2 + dy**2)\n",
    "    \n",
    "    # 获得梯度方向\n",
    "    theta = np.rad2deg(np.arctan2(dy, dx))\n",
    "    \n",
    "    # 将梯度方向的大小调整为0-360之间\n",
    "    theta %= 360\n",
    "\n",
    "    return G, theta\n",
    "\n",
    "G, theta = gradient(smoothed)\n",
    "\n",
    "# 展示利用图像梯度得到的结果\n",
    "plot_image(np.uint8(G), 'Gradient magnitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3.2.3  边缘强度非最大值抑制"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从上面的结果我们可以发现，直接从梯度得到的边缘线非常粗，显得很模糊，因此我们需要应用非最大值抑制将边缘细化。非最大值抑制沿着每个梯度方向找到边缘强度局部极大的像素点予以保留，去除掉该方向上边缘强度非局部极大的像素点。\n",
    "\n",
    "我们先将图像中每个像素的梯度方向量化为8个固定方向，即$\\{\\frac{k}{8}\\pi\\}_{k=0}^7$。接着，将当前像素的边缘强度与沿该像素的梯度方向（包括正反方向）的邻域内像素点的边缘强度进行比较。如果当前像素的边缘强度最大，则保留该像素为边缘中像素，反之则将其从边缘中移除。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-3ed4361a036a522611af90f887517f90_1440w.jpg\" width=200>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图5-4 非最大值抑制</div>\n",
    "</center>\n",
    "\n",
    "例如在图 5-4 中，我们考察当前点$[i,j]$，它的梯度方向如图 5-4 中箭头所示。我们分别沿着该梯度方向的正反方向，找到当前点$[i,j]$的两个邻域点$[i^{\\prime}, j^{\\prime}]$和$[i^{\\prime\\prime}, j^{\\prime\\prime}]$（通常选取当前点的8-邻域点）。如果当前点$[i,j]$的边缘强度在这三个点中最大，即$\\mathbf{E}_s[i,j|\\mathbf{I}]=\\max\\{\\mathbf{E}_s[i,j|\\mathbf{I}],\\mathbf{E}_s[i^{\\prime},j^{\\prime}|\\mathbf{I}],\\mathbf{E}_s[i^{\\prime\\prime},j^{\\prime\\prime}|\\mathbf{I}]\\}$，则保留当前像素在边缘中，反之将其从边缘中移除，并设$\\mathbf{E}_s[i,j|\\mathbf{I}]=0$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:35:19.921417Z",
     "iopub.status.busy": "2022-09-20T09:35:19.921099Z",
     "iopub.status.idle": "2022-09-20T09:35:20.814155Z",
     "shell.execute_reply": "2022-09-20T09:35:20.813262Z",
     "shell.execute_reply.started": "2022-09-20T09:35:19.921383Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 286088 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 非最大值抑制算法\n",
    "def non_maximum_suppression(G, theta):\n",
    "    \n",
    "    # 获得梯度图的大小\n",
    "    H, W = G.shape\n",
    "\n",
    "    # 将梯度方向投影到最近的45度角空间中\n",
    "    theta = np.floor((theta + 22.5) / 45) * 45\n",
    "    theta %= 360\n",
    "    \n",
    "    # 最后输出的图像梯度\n",
    "    out = G.copy()\n",
    "    \n",
    "    for i in range(1, H-1):\n",
    "        for j in range(1,W-1):\n",
    "            \n",
    "            # 当前像素点的角度大小，我们可以将其分为4类\n",
    "            angle = theta[i,j]\n",
    "            \n",
    "            if angle == 0 or angle == 180:\n",
    "                ma = max(G[i, j-1], G[i, j+1])\n",
    "                \n",
    "            elif angle == 45 or angle == 45 + 180:\n",
    "                ma = max(G[i-1, j-1], G[i+1, j+1])\n",
    "                \n",
    "            elif angle == 90 or angle == 90 + 180:\n",
    "                ma = max(G[i-1, j], G[i+1,j])\n",
    "                \n",
    "            elif angle == 135 or angle == 135 + 180:\n",
    "                ma = max(G[i-1, j+1], G[i+1, j-1])\n",
    "                \n",
    "            else:\n",
    "                print(angle)\n",
    "                raise\n",
    "            # ma是当前像素点相邻两个邻居点的像素梯度的最大值\n",
    "            \n",
    "            # 如果ma的值大于当前像素点的梯度值，则认为当前点非边缘，并将该点的梯度值设置为0\n",
    "            if ma > G[i,j]:\n",
    "                out[i,j]=0\n",
    "    return out\n",
    "\n",
    "\n",
    "nms = non_maximum_suppression(G, theta)\n",
    "plot_image(np.uint8(nms), 'Non-maximum suppressed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们能够很容易地从实验结果中看出，检测得到的图像边缘变得更加清晰了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3.2.4 基于滞后阈值化的边缘定位"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过非最大值抑制，我们将边缘细化成边缘线，接下来需要确定最终的边缘位置。我们需要对缘强度$\\mathbf{E}_s$进行阈值二值化，将其转化为非零即一的边缘图$\\mathbf{E}_b$，如果$\\mathbf{E}_b[i,j|\\mathbf{I}]$=1，那么像素点$[i,j]$是边缘点，否则不是边缘点。最简单的阈值化的方法是单一阈值法，即取一个单一阈值$T$对其二值化：\n",
    "\n",
    "$$\\mathbf{E}_b[i,j|\\mathbf{I}]= = \\begin{cases}\n",
    "1, & \\mathbf{E}_s[i,j|\\mathbf{I}]\\ge T \\\\\n",
    "0. & \\mathbf{E}_s[i,j|\\mathbf{I}] < T \\\\\n",
    "\\end{cases}$$\n",
    "\n",
    "但是这种单一阈值法很容易造成断裂的边缘。为了解决这个问题，Canny边缘检测算法使用双阈值，提出基于滞后阈值化的方法。这两个阈值一高一低，分别为$T_h$和$T_l$（$T_l<T_h$）。如果像素点$[i,j]$的边缘强度$\\mathbf{E}_s[i,j|\\mathbf{I}]$大于高阈值，则该点被标记为强边缘点，$\\mathbf{E}_b[i,j|\\mathbf{I}]=1$；而如果边缘强度小于低阈值，则该点为非边缘点，$\\mathbf{E}_b[i,j|\\mathbf{I}]=0$；如果边缘强度在两个阈值之间，则该点被标记为弱边缘点，需要进一步确定$\\mathbf{E}_b[i,j|\\mathbf{I}]$取值多少：\n",
    "\n",
    "$$\\mathbf{E}_b[i,j|\\mathbf{I}]= = \\begin{cases}\n",
    "1, & \\mathbf{E}_s[i,j|\\mathbf{I}]\\geq T_h \\\\\n",
    "0, & \\mathbf{E}_s[i,j|\\mathbf{I}] < T_l \\\\\n",
    "\\text{to be determined}. & T_l \\leq \\mathbf{E}_s[i,j|\\mathbf{I}] < T_h \\\\\n",
    "\\end{cases}$$\n",
    "\n",
    "这就是所谓的滞后阈值化。我们现在编程实现通过双阈值确定强弱边缘："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:36:20.026400Z",
     "iopub.status.busy": "2022-09-20T09:36:20.026093Z",
     "iopub.status.idle": "2022-09-20T09:36:21.063101Z",
     "shell.execute_reply": "2022-09-20T09:36:21.062198Z",
     "shell.execute_reply.started": "2022-09-20T09:36:20.026368Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 117084 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 213236 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 双阈值确定强弱边缘\n",
    "def double_thresholding(img, high, low):\n",
    "    # 初始化强弱边缘为bool矩阵\n",
    "    strong_edges = np.zeros(img.shape, dtype=np.bool_)\n",
    "    weak_edges = np.zeros(img.shape, dtype=np.bool_)\n",
    "    \n",
    "    # 获得输入图像大小\n",
    "    a,b = img.shape\n",
    "    for i in range(a):\n",
    "        for j in range(b):\n",
    "            # 大于Th阈值则为强边缘\n",
    "            if img[i,j] > high:\n",
    "                strong_edges[i,j] = 1\n",
    "                \n",
    "            # 小于Th阈值大于Tl阈值则为弱边缘\n",
    "            elif img[i,j] <= high and img[i,j] > low:\n",
    "                weak_edges[i, j] = 1\n",
    "\n",
    "    return strong_edges, weak_edges\n",
    "\n",
    "\n",
    "low_threshold = 3\n",
    "high_threshold = 6\n",
    "\n",
    "strong_edges, weak_edges = double_thresholding(nms, high_threshold, low_threshold)\n",
    "\n",
    "# 返回强弱边缘叠加的边缘图\n",
    "edges = strong_edges * 1.0 + weak_edges * 0.5\n",
    "\n",
    "plot_image(strong_edges, 'Strong Edges')\n",
    "plot_image(edges, 'Strong+Weak Edges')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "强边缘点是已确定的边缘点，而弱边缘点，当且仅当它们在已确定的边缘点的邻域内（如8-邻域）时，它们才被确定为边缘点。接下来，我们实现对弱边缘点的确定过程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:36:41.236410Z",
     "iopub.status.busy": "2022-09-20T09:36:41.235771Z",
     "iopub.status.idle": "2022-09-20T09:36:41.494605Z",
     "shell.execute_reply": "2022-09-20T09:36:41.493468Z",
     "shell.execute_reply.started": "2022-09-20T09:36:41.236364Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 166708 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 获得在(x，y)处点的邻居\n",
    "def get_neighbors(x, y, H, W):\n",
    "    neighbors = []\n",
    "    for i in (x - 1, x, x + 1):\n",
    "        for j in (y - 1, y, y + 1):\n",
    "            if i >= 0 and i < H and j >= 0 and j < W:\n",
    "                if (i == x and j == y):\n",
    "                    continue\n",
    "                neighbors.append((i, j))\n",
    "\n",
    "    return neighbors\n",
    "\n",
    "\n",
    "# 判断弱边缘点是否为边缘\n",
    "def link_edges(strong_edges, weak_edges):\n",
    "    # 获得图的大小\n",
    "    H, W = strong_edges.shape\n",
    "    \n",
    "    # 获得强边缘点的位置\n",
    "    indices = np.stack(np.nonzero(strong_edges)).T\n",
    "    \n",
    "    # 初始化最终的结果，为bool矩阵\n",
    "    edges = np.zeros((H, W), dtype=np.bool_)\n",
    "\n",
    "    weak_edges = np.copy(weak_edges)\n",
    "    edges = np.copy(strong_edges)\n",
    "    \n",
    "    # 将所有强边缘点的坐标组合成一个列表\n",
    "    q = [(i,j) for i in range(H) for j in range(W) if strong_edges[i,j]]\n",
    "    \n",
    "    while q:\n",
    "        \n",
    "        # pop()函数返回列表中的末尾元素，并将其从列表中删除\n",
    "        i, j = q.pop()\n",
    "        \n",
    "        # (a, b)是(i, j)的邻居\n",
    "        for a, b in get_neighbors(i, j, H, W):\n",
    "            \n",
    "            # 如果当前点是弱边缘点，由于其在强边缘点的邻域内，即可认为该点是边缘点\n",
    "            if weak_edges[a][b]:\n",
    "                # 为了避免对同一弱边缘点重复判断，便将其值设置为0\n",
    "                weak_edges[a][b] = 0\n",
    "                \n",
    "                # 在边缘点中添加该弱边缘点\n",
    "                edges[a][b] = 1\n",
    "                \n",
    "                # 由于该弱边缘点已经变成边缘点，因此需要将其加入q中，判断该点周围是否存在弱边缘点\n",
    "                q.append((a,b))\n",
    "                \n",
    "    return edges\n",
    "\n",
    "\n",
    "# 获得最终边缘图\n",
    "edges = link_edges(strong_edges, weak_edges)\n",
    "plot_image(edges, 'Final image')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "综上，我们将Canny边缘检测器封装成如下函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:37:16.571919Z",
     "iopub.status.busy": "2022-09-20T09:37:16.571615Z",
     "iopub.status.idle": "2022-09-20T09:37:18.504669Z",
     "shell.execute_reply": "2022-09-20T09:37:18.503735Z",
     "shell.execute_reply.started": "2022-09-20T09:37:16.571885Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 251640 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 167232 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def canny(img, kernel_size=5, sigma=1.5, high=6, low=3):\n",
    "    \n",
    "    # 获取高斯滤波器\n",
    "    kernel = gaussian_kernel(kernel_size, sigma)\n",
    "    # 将图像进行高斯滤波\n",
    "    smoothed = cv2.filter2D(img, -1, kernel)    \n",
    "    # 得到图像的梯度图\n",
    "    G, theta = gradient(smoothed)\n",
    "    # 对梯度图进行nms\n",
    "    nms = non_maximum_suppression(G, theta)\n",
    "    # 获得强弱边缘信息\n",
    "    strong_edges, weak_edges = double_thresholding(nms, high, low)\n",
    "    # 对弱边缘进行分类，得到最终的边缘图\n",
    "    edge = link_edges(strong_edges, weak_edges)\n",
    "\n",
    "    return edge\n",
    "\n",
    "\n",
    "plot_image(img, 'Original image')\n",
    "canny_img = canny(img)\n",
    "plot_image(canny_img, 'Edges of image')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们发现上面的图像中，这只蜥蜴的头部的边缘都被检测了出来，但是也包含了很多由细小纹理变化导致的边缘。如前所述，Canny算子可以通过调节高斯滤波的标准差以及核大小得到多尺度输出，从而使得边缘检测的结果对细小纹理变化不敏感。接下来，我们试着调整高斯核的大小与标准差，观察比较不同的平滑程度对边缘检测结果的影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-20T09:38:56.105174Z",
     "iopub.status.busy": "2022-09-20T09:38:56.104861Z",
     "iopub.status.idle": "2022-09-20T09:39:01.469897Z",
     "shell.execute_reply": "2022-09-20T09:39:01.468997Z",
     "shell.execute_reply.started": "2022-09-20T09:38:56.105139Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 129556 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 97852 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 76016 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "canny2 = canny(img, kernel_size=7, sigma=2)\n",
    "plot_image(canny2, 'Kernel Size: 7, Sigma: 2')\n",
    "\n",
    "canny2 = canny(img, kernel_size=9, sigma=2.5)\n",
    "plot_image(canny2, 'Kernel Size: 9, Sigma: 2.5')\n",
    "\n",
    "canny2 = canny(img, kernel_size=9, sigma=3)\n",
    "plot_image(canny2, 'Kernel Size: 9, Sigma: 3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们发现，当高斯核标准差变大时，边缘检测的结果会对细小纹理的变化更为鲁棒，但是重要的边缘也会出现丢失情况。因此，我们需要权衡图像平滑和降低噪声之间的权重。我们将上述函数全部封装进utils.py中。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.4 小结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一章中，我们介绍了边缘检测的数学模型，并学习了两种边缘检测算法——Sobel算子和Canny算法。这两种算法都是基于差分计算图像的梯度，边缘的强度和方向都可以通过梯度计算得到。Sobel算子虽然简单，但是检测结果容易受图像噪声的影响。Canny算子是应用最广泛的一种边缘检测算法，通过改变卷积核大小调节对图像的平滑程度，输出多尺度的边缘检测结果，并通过非最大值抑制和滞后阈值化方法分别提升边缘检测的位置准确性和对噪声的鲁棒性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.5 习题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 习题 5.1：本章中我们所学习使用的边缘检测方法都是基于图像一阶微分的，那么图像的二阶微分是否也可以实现边缘检测呢？如果可以，请实现基于二阶微分的边缘检测算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 引用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Sobel, Irwin. (2014). An Isotropic 3x3 Image Gradient Operator. Presentation at Stanford A.I. Project 1968.\n",
    "\n",
    "[2] Canny J. A computational approach to edge detection[J]. IEEE Transactions on pattern analysis and machine intelligence, 1986 (6): 679-698."
   ]
  }
 ],
 "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
