{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 代码实现\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, \n",
    "                             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, \n",
    "                                     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*3\n",
    "    # 即[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) \\\n",
    "                    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) \\\n",
    "                    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,\n",
    "                        border, num_layers, g_image, corners, \n",
    "                        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空间中相邻的3层\n",
    "        img1, img2, img3 = dog_images[layer-1:layer+2]\n",
    "        # 对这3层进行切片，得到大小为3*3的3层\n",
    "        cube = np.array([img1[x-1:x+2, y-1:y+2], \n",
    "                         img2[x-1:x+2, y-1:y+2],\n",
    "                         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 \\\n",
    "            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 \\\n",
    "                        or y < border \\\n",
    "                        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",
    "    # 超出空间大小或者未不收敛，直接返回\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的黑塞矩阵\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 \\\n",
    "            (xy_hessian_trace ** 2) / xy_hessian_det < \\\n",
    "            ((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*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*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] \n",
    "           - cube[1, 0, 2] + cube[1, 0, 0]) / 4\n",
    "    dxs = (cube[2, 1, 2] - cube[2, 1, 0] \n",
    "           - cube[0, 1, 2] + cube[0, 1, 0]) / 4\n",
    "    dys = (cube[2, 2, 1] - cube[2, 0, 1] \n",
    "           - 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(\n",
    "                        -0.5 / ((size*1.5) ** 2) \n",
    "                        * (i ** 2 + j ** 2)\n",
    "                        )  \n",
    "                    histogram_index = int(\n",
    "                        round(orien * num_bins / 360.)\n",
    "                        )\n",
    "                    histogram[histogram_index % num_bins] += \\\n",
    "                                weight * value\n",
    "\n",
    "    # 对直方图进行高斯平滑\n",
    "    for n in range(num_bins):\n",
    "        smooth_histogram[n] = (6 * histogram[n] \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",
    "    # 开始寻找块状区域\n",
    "    threshold = 0.02\n",
    "    corners = []\n",
    "    scales = []\n",
    "    orientations = []\n",
    "    layers = []\n",
    "    \n",
    "    for layer, (image1, image2, image3) in \\\n",
    "        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], \n",
    "                                   image2[x-1:x+2, y-1:y+2], \n",
    "                                   image3[x-1:x+2, y-1:y+2], \n",
    "                                   threshold):\n",
    "                    # 如果是候选点，则进行进一步定位筛选，并返回其信息\n",
    "                    computeBlobAttribute(x, y, layer+1, dog_images, \n",
    "                                         sigma, threshold, border, \n",
    "                                         num_layers, g_images[layer], \n",
    "                                         corners, scales, orientations, \n",
    "                                         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, \n",
    "                        orientations, layers):\n",
    " \n",
    "    if len(corners) != len(scales) \\\n",
    "        or len(corners) != len(orientations):\n",
    "        raise ValueError(\n",
    "          '`corners`, `scales` and `orientations` \\\n",
    "          must all have the same length.')\n",
    "    \n",
    "    descriptors_list = []\n",
    "    \n",
    "    for pt, size, orien, layer in \\\n",
    "        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, \n",
    "            np.sqrt(g_image.shape[0]**2+g_image.shape[1]**2))))\n",
    "        \n",
    "        # 用列表依次存储块状区域内所有点的信息\n",
    "        i_index = []\n",
    "        j_index = []\n",
    "        value_list = []\n",
    "        orien_index = []\n",
    "        \n",
    "        # 三维数组存储16个窗口的8个方向，为防止计算时边界溢出，\n",
    "        # 在行、列的首尾各扩展一次\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)) \\\n",
    "                        + i * np.cos(np.deg2rad(orien))\n",
    "                j_rotate = j * np.cos(np.deg2rad(orien)) \\\n",
    "                        - 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 \\\n",
    "                    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 \\\n",
    "                        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] \\\n",
    "                            - g_image[i_inimg, j_inimg - 1]\n",
    "                        dy = g_image[i_inimg - 1, j_inimg] \\\n",
    "                            - g_image[i_inimg + 1, j_inimg]\n",
    "                        grad_value = np.sqrt(dx**2 + dy**2)\n",
    "                        # 计算梯度方向\n",
    "                        grad_orien = np.rad2deg(\n",
    "                            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 * (\n",
    "                            (i_rotate / win_s) ** 2 + \n",
    "                            (j_rotate / win_s) ** 2)\n",
    "                            )\n",
    "                        value_list.append(g_weight * grad_value)\n",
    "                        # 将梯度方向投影到8个方向\n",
    "                        # 这里的grad_orien是原图像中的梯度\n",
    "                        # 需要叠加上旋转之后的角度\n",
    "                        orien_index.append(\n",
    "                            (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个整数坐标都有一定的贡献\n",
    "    # （可以将该浮点坐标想象成一个正方体中的一点，它对正方体的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('sift1.png', cv2.IMREAD_COLOR)\n",
    "img2 = cv2.imread('sift2.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 = \\\n",
    "                                    detect_blobs(gray1)\n",
    "g_images2, corners2, scales2, orientations2, layers2 = \\\n",
    "                                    detect_blobs(gray2)\n",
    "\n",
    "# 获得两张图像的SIFT特征点\n",
    "descriptors1 = compute_descriptors(g_images1, corners1, scales1, \n",
    "                                   orientations1, layers1)\n",
    "descriptors2 = compute_descriptors(g_images2, corners2, scales2, \n",
    "                                   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": [
    "接下来，需要将寻找出的特征点进行匹配"
   ]
  },
  {
   "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] \\\n",
    "                                      - 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), \n",
    "                   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]), \n",
    "                   int(corners1[m[0]][1] + hres))\n",
    "            pt2 = (int(corners2[m[1]][0] + w1), \n",
    "                   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('sift1.png', cv2.IMREAD_GRAYSCALE)\n",
    "image2 = cv2.imread('sift2.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": [
    "可以很明显的看出，两幅图像中相同的特征点完美的匹配在了一起。再这之后，便可以进行图像拼接等一系列后续操作了。"
   ]
  }
 ],
 "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.0.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
