{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 代码实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一节中，我们将编程实现图像拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:14:27.878595Z",
     "iopub.status.busy": "2022-09-16T03:14:27.878175Z",
     "iopub.status.idle": "2022-09-16T03:14:27.904993Z",
     "shell.execute_reply": "2022-09-16T03:14:27.903759Z",
     "shell.execute_reply.started": "2022-09-16T03:14:27.878559Z"
    }
   },
   "outputs": [],
   "source": [
    "from utils import *\n",
    "\n",
    "# 利用RANSAC算法剔除无效点\n",
    "def compute_affine_xform(corners1, corners2, matches):\n",
    "    \n",
    "    # 设置最大迭代次数\n",
    "    iteration = 50\n",
    "    M_list = []\n",
    "    \n",
    "    # 存放每一次模型的正常样本个数\n",
    "    inlier_num_list = []\n",
    "    for _ in range(iteration):\n",
    "        \n",
    "        # 随机选择4个特征点对来拟合模型\n",
    "        sample_index = random.sample(range(len(matches)), 4)\n",
    "        x1_s, y1_s = corners1[matches[sample_index[0]][0]]\n",
    "        x1_t, y1_t = corners2[matches[sample_index[0]][1]]\n",
    "        x2_s, y2_s = corners1[matches[sample_index[1]][0]]\n",
    "        x2_t, y2_t = corners2[matches[sample_index[1]][1]]\n",
    "        x3_s, y3_s = corners1[matches[sample_index[2]][0]]\n",
    "        x3_t, y3_t = corners2[matches[sample_index[2]][1]]\n",
    "        x4_s, y4_s = corners1[matches[sample_index[3]][0]]\n",
    "        x4_t, y4_t = corners2[matches[sample_index[3]][1]]\n",
    "        \n",
    "        # 编写矩阵A\n",
    "        A = np.array([[x1_s, y1_s, 1, 0, 0, 0, \n",
    "                       -x1_t*x1_s, -x1_t*y1_s, -x1_t],\n",
    "                      [0, 0, 0, x1_s, y1_s, 1, \n",
    "                       -y1_t*x1_s, -y1_t*y1_s, -y1_t],\n",
    "                      [x2_s, y2_s, 1, 0, 0, 0, \n",
    "                       -x2_t*x2_s, -x2_t*y2_s, -x2_t],\n",
    "                      [0, 0, 0, x2_s, y2_s, 1, \n",
    "                       -y2_t*x2_s, -y2_t*y2_s, -y2_t],\n",
    "                      [x3_s, y3_s, 1, 0, 0, 0, \n",
    "                       -x3_t*x3_s, -x3_t*y3_s, -x3_t],\n",
    "                      [0, 0, 0, x3_s, y3_s, 1, \n",
    "                       -y3_t*x3_s, -y3_t*y3_s, -y3_t],\n",
    "                      [x4_s, y4_s, 1, 0, 0, 0, \n",
    "                       -x4_t*x4_s, -x4_t*y4_s, -x4_t],\n",
    "                      [0, 0, 0, x4_s, y4_s, 1, \n",
    "                       -y4_t*x4_s, -y4_t*y4_s, -y4_t]\n",
    "                      ])\n",
    "        \n",
    "        # 求解A的特征向量\n",
    "        _,_, v = np.linalg.svd(A)\n",
    "        \n",
    "        # 取最小特征值对应的特征向量作为最终的结果\n",
    "        M = np.reshape(v[-1], (3, 3))\n",
    "        \n",
    "        inlier_num = 0\n",
    "        # 统计正常样本个数\n",
    "        for (index1, index2) in matches:\n",
    "            # coord是齐次坐标系的坐标\n",
    "            coord1 = [corners1[index1][0], corners1[index1][1], 1]\n",
    "            coord2 = [corners2[index2][0], corners2[index2][1], 1]\n",
    "            # 计算将coord1齐次变换之后的坐标\n",
    "            mapcoor = np.dot(M, coord1)\n",
    "            # 将齐次坐标系中的w置为1\n",
    "            mapcoor = mapcoor / mapcoor[-1]\n",
    "            if np.linalg.norm(coord2 - mapcoor) < 5:\n",
    "                inlier_num += 1\n",
    "                \n",
    "        # 将正常样本个数和对应的变换矩阵M记录在列表里\n",
    "        M_list.append(M)\n",
    "        inlier_num_list.append(inlier_num)\n",
    "    \n",
    "    # 获取列表中正常样本值最大的元素的下标\n",
    "    best_index = np.argmax(inlier_num_list)\n",
    "    # 获取对应的变换矩阵\n",
    "    xform = M_list[best_index].astype(np.float64)\n",
    "    \n",
    "    # 统计属于模型异常样本的特征点对\n",
    "    outlier_labels = []\n",
    "    for (index1, index2) in matches:\n",
    "        coord1 = [corners1[index1][0], corners1[index1][1], 1]\n",
    "        coord2 = [corners2[index2][0], corners2[index2][1], 1]\n",
    "        mapcoor = np.dot(xform, coord1)\n",
    "        mapcoor = mapcoor/mapcoor[-1]\n",
    "        if np.linalg.norm(coord2 - mapcoor) < 12:\n",
    "            outlier_labels.append(1)\n",
    "        else:\n",
    "            outlier_labels.append(0)\n",
    "    \n",
    "    return xform, outlier_labels  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了在绘制时能够体现出异常点对，我们改写第7章中的draw_matches()函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_matches(image1, image2, corners1, corners2, \n",
    "                 matches, outliers=None):\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",
    "            if outliers:\n",
    "                cv2.line(match_image, pt1, pt2, (255, 0, 0))\n",
    "            else:\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",
    "            if outliers:\n",
    "                cv2.line(match_image, pt1, pt2, (255, 0, 0))\n",
    "            else:\n",
    "                cv2.line(match_image, pt1, pt2, (0, 0, 255))\n",
    "        \n",
    "    return match_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之后，我们编写图像变换与缝合的相关函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 图像变换与缝合\n",
    "def stitch_images(image1, image2, xform):\n",
    "    # 将图像按照变换矩阵的方向进行变换\n",
    "    new_image = cv2.warpPerspective(image1, xform, \n",
    "                        (image2.shape[1],image2.shape[0]))\n",
    "    image = new_image.copy()\n",
    "    \n",
    "    # 获取变换后的图像位置坐标\n",
    "    h, w, _ = image1.shape\n",
    "    pts = np.float32([[0, 0],\n",
    "                      [0, h - 1],\n",
    "                      [w - 1, h - 1],\n",
    "                      [w - 1, 0]]).reshape(-1, 1, 2)\n",
    "    dst = cv2.perspectiveTransform(pts, xform)\n",
    "    \n",
    "    # 画出缝合线\n",
    "    image = cv2.polylines(image2, [np.int32(dst)], True, \n",
    "                          255, 1, cv2.LINE_AA)\n",
    "    \n",
    "    # 将重叠部分的图像求平均\n",
    "    image = np.where(new_image>0, 0.5*new_image+0.5*image, image)\n",
    "    return image.astype(np.uint8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "导入两张图像$\\boldsymbol{A}$和$\\boldsymbol{B}$，我们先计算图像的变换矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:14:27.906774Z",
     "iopub.status.busy": "2022-09-16T03:14:27.906345Z",
     "iopub.status.idle": "2022-09-16T03:15:33.710850Z",
     "shell.execute_reply": "2022-09-16T03:15:33.709748Z",
     "shell.execute_reply.started": "2022-09-16T03:14:27.906730Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 168644 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 读取图像\n",
    "img1 = cv2.imread('stitch1.jpg', cv2.IMREAD_COLOR)\n",
    "img2 = cv2.imread('stitch2.jpg', cv2.IMREAD_COLOR)\n",
    "gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) / 255.0\n",
    "gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY) / 255.0\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",
    "# 计算两张图像的特征点\n",
    "descriptors1 = compute_descriptors(g_images1, corners1, scales1, \n",
    "                                   orientations1, layers1)\n",
    "descriptors2 = compute_descriptors(g_images2, corners2, scales2, \n",
    "                                   orientations2, layers2)\n",
    "\n",
    "# 匹配两张图像中的特征点\n",
    "matches = match_descriptors(descriptors1, descriptors2)\n",
    "\n",
    "image1 = cv2.imread('stitch1.jpg', cv2.IMREAD_GRAYSCALE)\n",
    "image2 = cv2.imread('stitch2.jpg', cv2.IMREAD_GRAYSCALE)\n",
    "\n",
    "# xform为变换矩阵，outlier_labels为模型的异常样本\n",
    "xform, outlier_labels = compute_affine_xform(corners1, corners2, \n",
    "                                             matches) \n",
    "# 展示两张图像特征点的匹配\n",
    "match_image = draw_matches(image1, image2, corners1, corners2, \n",
    "                           matches, outliers=outlier_labels)\n",
    "\n",
    "plt.imshow(match_image)\n",
    "plt.axis('off')\n",
    "plt.plot()    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之后，我们将其中一张图像进行变换，并与另一张图像拼接。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-09-16T03:15:33.713310Z",
     "iopub.status.busy": "2022-09-16T03:15:33.712886Z",
     "iopub.status.idle": "2022-09-16T03:15:33.951886Z",
     "shell.execute_reply": "2022-09-16T03:15:33.950889Z",
     "shell.execute_reply.started": "2022-09-16T03:15:33.713272Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 207692 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 假设 img1 和 img2 是 BGR 格式的图像，将它们转换为 RGB\n",
    "img1_rgb = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB)\n",
    "img2_rgb = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)\n",
    "\n",
    "# 使用转换后的 RGB 图像拼接\n",
    "image = stitch_images(img1_rgb, img2_rgb, xform)\n",
    "\n",
    "# 确保图像是 uint8 格式，然后显示\n",
    "plt.imshow(image.astype(np.uint8))\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "test",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
