4.0-tl-feature-matching.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 5: Feature Extraction and Matching\n",
    "\n",
    "In this exercise, you will learn:\n",
    "\n",
    "- How to find key points in images and describe their features\n",
    "- How to match key points between two views of the same scene\n",
    "\n",
    "For some functions we already provide some draft implementation that you just need to complete. This is supposed to help you identifying the next steps.\n",
    "\n",
    "\n",
    "## 0. Setup\n",
    "\n",
    "Load all libraries and both views that we will work with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.typing as npt\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load images and resize\n",
    "view1 = np.array(Image.open(\"data/exercise3/view1.png\")) / 255\n",
    "view2 = np.array(Image.open(\"data/exercise3/view2.png\")) / 255\n",
    "\n",
    "# Check resolution\n",
    "print(f\"View 1 resolution: {view1.shape[0]}x{view1.shape[1]}\")\n",
    "print(f\"View 2 resolution: {view2.shape[0]}x{view2.shape[1]}\")\n",
    "\n",
    "# Show both views\n",
    "_, axes = plt.subplots(1,2, figsize=(12,6))\n",
    "axes[0].imshow(view1)\n",
    "axes[1].imshow(view2)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Key Point Detection and Feature Extraction\n",
    "\n",
    "You might want to look again into the concepts explained in slide decks \"04_sparse_geometry_1_appearance_matching.pdf\" and \"03_image_processing_2.pdf\n",
    "\n",
    "In this part of the exercise, you will learn how to build your own key point detector and feature descriptor. Generally, you should rather use the implementation of traditional and well-known feature descriptors such as [SIFT](https://link.springer.com/article/10.1023/B:VISI.0000029664.99615.94), [ORB](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6126544), etc. But here, we will learn how to come up with our own method from scratch.\n",
    "\n",
    "As a reminder: Feature extraction consists of two step. First, there is detection/localization of key points. Second, computing a feature vector for each key point to uniquely describe it.\n",
    "\n",
    "### 1.1. Key Point Detection using Harris Corner Detector\n",
    "\n",
    "Complete the implementation below. First you have to implement some helper functions (```smooth_image```, ```image_gradient```, ```compute_structure_tensor_Q```, ```compute_harris_score```, ```non_maximum_suppression```, ```compute_angles```). Finally you detect keypoints by implementing the function ```detect_corner_keypoints``` (and call your implemented helper functions from there). It is supposed to be a function that takes an image and some hyperparameters as inputs and returns a list of coordinates describing the location of detected key points. Furthermore, the orientation of the key point shall be returned. For key point detection you should apply the Harris detector to the image as dicussed in the lecture. For key point orientation, you should return the angle of the local gradient at each detected key point. The angle is important so that corresponding key points from differently rotated view points can be aligned. Hint: Use np.arctan2(x,y) for angle computation to get full 360 degree orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this time you may utilize library implementations of\n",
    "# the image processing filters you learned about\n",
    "from scipy.ndimage import sobel, gaussian_filter, maximum_filter, convolve\n",
    "\n",
    "def to_grayscale(img: npt.ArrayLike) -> np.ndarray:\n",
    "    channel_weights = [0.299, 0.587, 0.114]\n",
    "    weighted_img = img * np.reshape(channel_weights, (1,1,3))\n",
    "    return np.sum(weighted_img, axis=2)\n",
    "\n",
    "def smooth_image(image: npt.ArrayLike, sigma: int) -> np.ndarray:\n",
    "    \"\"\"Smooth the image with a gaussian filter\"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def image_gradient(image: npt.ArrayLike, axis: int) -> np.ndarray:\n",
    "    \"\"\"Compute the image gradients with the sobel filter along the specified axis.\n",
    "\n",
    "    :param image: grayscale image\n",
    "    :param axis: specified axis\n",
    "    :return array with the image gradient along specified axis\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def compute_structure_tensor_Q(Ix: npt.ArrayLike, Iy: npt.ArrayLike, window_size: int) -> np.ndarray:\n",
    "    \"\"\"Compute the structure tensor Q.\n",
    "\n",
    "    :param Ix: horizontal image gradient\n",
    "    :param Iy: vertical image gradient\n",
    "    :param window_size: window size for calculation\n",
    "    :return array with structure tensor Q\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def compute_harris_score(Q: npt.ArrayLike, k: float) -> np.ndarray:\n",
    "    \"\"\"Compute the harris score.\n",
    "\n",
    "    :param Q: structure tensor\n",
    "    :param k: hyperparameter\n",
    "    :return array with harris scores\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def non_maximum_suppression(H: npt.ArrayLike, threshold: int) -> np.ndarray:\n",
    "    \"\"\"Non-maximum surpression.\n",
    "\n",
    "    :param H: Harris scores\n",
    "    :param threshold: H > threshold -> pixel is corner candidate\n",
    "    :return array keypoints as boolean mask\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def compute_angles(Ix: npt.ArrayLike, Iy: npt.ArrayLike, keypoint_mask: npt.ArrayLike) -> np.ndarray:\n",
    "    \"\"\"Compute the angles of the keypoints.\n",
    "\n",
    "    :param Ix: horizontal image gradient\n",
    "    :param Iy: vertical image gradient\n",
    "    :param keypoint_mask: array with the detected keypoints\n",
    "    :return array containing the angle for each keypoint\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "    \n",
    "def detect_corner_keypoints(img, window_size=5, k=0.05, threshold=10):\n",
    "    \"\"\"\n",
    "    This function takes an input image and computes harris corner locations and orientations\n",
    "    :param img: the image\n",
    "    :param window_size: size of the window in which we look for corners\n",
    "    :param k: factor to compute the harris score, usually between 0.04 and 0.06\n",
    "    :param threshold: Harris scores below this threshold are not considered\n",
    "    :return: a tuple where the first element is a 2D numpy array of keypoint coordinate pairs\n",
    "             and the second element is a numpy array of local gradient angle\n",
    "    \"\"\"\n",
    "    # 1. convert to grayscale (use the already implemented function above)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 2. smooth the image (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 3. compute image gradients in x and y direction (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 4. compute structure tensor Q (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 5. compute harris score (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 6. perform non-maximum surpression (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 7. compute angles (implement the function above and call it)\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 8. return list containing pairs of coordinates\n",
    "    # TODO: implement it"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute the key points and visualize the corners you detected..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_keypoints(img, points):\n",
    "    plt.imshow(img)\n",
    "    plt.scatter(points[:,0], points[:,1], marker='x', color='r')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "keypoints = {}\n",
    "for i, target_image in  enumerate([view1, view2]):\n",
    "    points, angles = detect_corner_keypoints(target_image, window_size=5, threshold=10)\n",
    "    keypoints[f\"view{i+1}\"] = (points, angles)\n",
    "    plot_keypoints(target_image, points)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you did everything correctly you should already be able to spot some key points that correspond to the same real point in both views. Sofar the key points are robust in terms of rotation and translation. To become more robust against changes in scale, we would need to further extend our implementation. However this not in the scope of the exercise. But here is a sketch of the idea:\n",
    "\n",
    "To achieve scale invariance we will detect corner key points at multiple scales of the input image. First, we will create a Gaussian pyramid of the image, i.e, we will blur the image and subsample by a factor of 2. We continue doing so until the image reaches a limit, say 32x32 pixels... Then, we will compute key points for each scale image and determine their location with respect to the original resolution. Afterwards, we supress duplicates."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Select patches around key points\n",
    " \n",
    "Right now, we are aware of the keypoints in the image. However, we would like to have a unique description of every keypoint to rediscover it also in other views of the same scene. Therefore, we need to take a look at the proximity of the key point and encode the information we find there. We choose a patch of size 16x16 around each key point to extract information and transform them into a descriptor. The patches are rotated according to the key point orientation.\n",
    "\n",
    "First, use the following plotting function to visulize the patches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.patches as patches\n",
    "from matplotlib.transforms import Affine2D\n",
    "\n",
    "def plot_patches(img, keypoints, angles):\n",
    "    \"\"\"\n",
    "    :param img: image\n",
    "    :param keypoints: Numpy array containig the keypoints\n",
    "    :param rotations: Numpy array of length=len(keypoints) containing\n",
    "                      the patch rotation\n",
    "    \"\"\"\n",
    "    scales = np.ones(len(keypoints))*16    \n",
    "    ax = plt.gca()\n",
    "    ax.imshow(img)\n",
    "    ax.scatter(points[:,0], points[:,1], marker='.', alpha=0.7, color='r')\n",
    "    for kp, angle, length in zip(keypoints, angles, scales):\n",
    "        rect = patches.Rectangle(kp - length / 2, length, length, linewidth=1,\n",
    "                                 edgecolor='r', facecolor='none')\n",
    "        transform = Affine2D().rotate_deg_around(*kp, angle) + ax.transData\n",
    "        rect.set_transform(transform)\n",
    "        ax.add_patch(rect)\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, target_image in  enumerate([view1, view2]):\n",
    "    points, angles = keypoints[f\"view{i+1}\"]\n",
    "    plot_patches(target_image, points, angles)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, complete the function that extracts the patches from a given view and returns them as 16x16 grayscale numpy arrays. Obviously, the array of patches shall have the same ordering as the key points. You may use `skimage.transform.rotate` to rotate the image before extracting the patches. Also, cropping the image before rotating will speed up your runtime.\n",
    "\n",
    "Hint: You can pad the image using `np.pad` if the patches are extracted add the borders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import rotate\n",
    "\n",
    "def extract_patches(img, keypoints, angles, patch_size=16):\n",
    "    \"\"\"\n",
    "    This function extracts oriented patches around the detected key points and returns\n",
    "    them as grayscale images\n",
    "    :param img: the input image\n",
    "    :param keypoints: the extracted keypoints\n",
    "    :param angles: the orientation of the keypoints\n",
    "    :param patch_size: the pixel length of each patch in x,y directions\n",
    "    :return: a 3D Numpy array containing all grayscale patches. The first dimension\n",
    "             is the number of key points/patches. The second and third is 'patch_size'.\n",
    "             \n",
    "    \"\"\"\n",
    "    # first convert to grayscale\n",
    "    if len(img.shape) == 3:\n",
    "        monochrome_img = to_grayscale(img)\n",
    "    else:\n",
    "        monochrome_img = img\n",
    "    \n",
    "    # TODO: add padding such that also patches at the borders can be extracted\n",
    "    pad_size =\n",
    "    padded_img =\n",
    "    \n",
    "    # TODO: extract patches\n",
    "    patches = []\n",
    "    for (x,y), angle in zip(keypoints, angles):\n",
    "        ...\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patches_view1 = extract_patches(view1, *keypoints[\"view1\"])\n",
    "patches_view2 = extract_patches(view2, *keypoints[\"view2\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_patch_crops(patches):\n",
    "    # define image grind\n",
    "    cols = 10\n",
    "    rows = int(np.ceil(len(patches) / cols))\n",
    "    \n",
    "    # create subplots\n",
    "    _, axes = plt.subplots(rows, cols, figsize=(16, rows * 3 ))\n",
    "    axes = axes.flatten()\n",
    "    for patch, ax in zip(patches,axes):\n",
    "        ax.imshow(patch)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_patch_crops(patches_view1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_patch_crops(patches_view2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you analyze the extracted patches from both views you should notice that some of them already look quite similar to patches from the other view. Now, we will construct a feature descriptor from them.\n",
    "\n",
    "### 1.3. Computing the Feature Descriptor\n",
    "\n",
    "There are different options to utilize the patch information and transform them into a feature representation. An intuitive approach is to compare the images by their color or color histograms. Yet, this reduces the robustness of our method to changes in lighting. Instead, we use image gradients as motivated in the [SIFT](https://link.springer.com/article/10.1023/B:VISI.0000029664.99615.94) paper and in the lecture.\n",
    "\n",
    "To implement this descriptor you may reuse most of your code from the `extract_patches` method in the previous section. You again extract patches but this time not from the grayscale image but from the gradient images, i. e., the gradient image in $x$ and $y$ direction. Afterwards, you compute the gradient length and orientation for every pixel in the patch. The gradient magnitudes are weighted by a gaussian kernel of the same size as the patch and $\\sigma=0.5 \\cdot patch size$. Finally, the patch is separated into a 4x4 grid where each cell consists of an 4x4 pixels subpatch. For each of these subpatches an 8-bin histogram is computed over the orientations of the gradients in that cell. For example, all gradients with orientation between 0 - 45 degrees are counted in the first histogram bin. All gradients with orientation between 45-90 in the second and so on. Instead of increasing the counter by one for each gradient, we increase it by the weighted magnitude. In total, we have eight counts per cell and 16 cells in total which results in a 128-dimenstional feature vector which we will use to describe the keypoint. In the end, this feature vector shall be nomalized to unit length.\n",
    "\n",
    "Hint: For histogram binning you can use `np.histogram` which already provides you with weighted bin counting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "\n",
    "# you may use this function to get the gaussian weights to weight the gradient magnitudes\n",
    "def get_gaussian_weights(std, size):\n",
    "    mvn = multivariate_normal(mean=[0,0], cov=[[std,0.],[0., std]])\n",
    "    xs = np.arange((1 - size) / 2, (size + 1) / 2)\n",
    "    xs, ys = np.meshgrid(xs,xs)\n",
    "    pos = np.dstack((xs, ys))\n",
    "    return mvn.pdf(pos)\n",
    "\n",
    "def get_gradient_orientations(grad_patches_x: npt.ArrayLike, grad_patches_y: npt.ArrayLike) -> np.ndarray:\n",
    "    \"\"\"Compute the gradient orientation.\n",
    "\n",
    "    :param grad_patches_x: gradient patch in x direction [number of keypoints, patch_size, patch_size]\n",
    "    :param grad_patches_y: gradient patch in y direction [number of keypoints, patch_size, patch_size]\n",
    "    :return array with the gradient orientations [number of keypoints, patch_size, patch_size]\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def get_gradient_magnitudes(grad_patches_x: npt.ArrayLike, grad_patches_y: npt.ArrayLike) -> np.ndarray:\n",
    "    \"\"\"Compute the gradient magnitudes.\n",
    "\n",
    "    :param grad_patches_x: gradient patch in x direction [number of keypoints, patch_size, patch_size]\n",
    "    :param grad_patches_y: gradient patch in y direction [number of keypoints, patch_size, patch_size]\n",
    "    :return array with the gradient magnitudes [number of keypoints, patch_size, patch_size]\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def add_magnitude_weighting(magnitudes: npt.ArrayLike, patch_size: int):\n",
    "    \"\"\"Add weighting for the magnitudes (with get_gaussian_weights() function).\n",
    "\n",
    "    :param magnitudes: gradient magnitudes [number of keypoints, patch_size, patch_size]\n",
    "    :param patch_size: patch size\n",
    "    :return array with weighted magnitudes [number of keypoints, patch_size, patch_size]\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def create_features(keypoints: npt.ArrayLike, magnitudes: npt.ArrayLike, orientations: npt.ArrayLike, patch_size: int, cell_size: int, histogram_bins: int) -> np.ndarray:\n",
    "    \"\"\"Create feature vectors.\n",
    "\n",
    "    :param keypoints: keypoint coordinates\n",
    "    :param magnitudes: magnitudes\n",
    "    :param orientations: orientations\n",
    "    :param patch_size: patch size\n",
    "    :param cell_size: cell size\n",
    "    :param histogram_bins: histogram_bins\n",
    "    :return array with feature vectors [number of keypoints, feature vector size]\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def create_descriptors(img, keypoints, angles, patch_size=16, cell_size=4, histogram_bins=8):\n",
    "    \"\"\"\n",
    "    This function creates descriptors from oriented patches around each key point\n",
    "    :param img: the input image\n",
    "    :param keypoints: the extracted keypoints\n",
    "    :param angles: the orientation of the keypoints\n",
    "    :param patch_size: the pixel length of each patch in x,y directions\n",
    "    :param cell_size: the size of each gradient histogram cell\n",
    "    :param histogram_bins: the number of bins per histogram\n",
    "    :return: a 2D Numpy array containing all feature descriptors. The first dimension\n",
    "             is the number of key points/patches. The second and third is 'patch_size'.\n",
    "             \n",
    "    \"\"\"\n",
    "    \n",
    "    assert patch_size % cell_size == 0, \"patch_size must be evenly divisible by cell_size\"\n",
    "    \n",
    "    # 1. first convert to grayscale\n",
    "    monochrome_img = to_grayscale(img)\n",
    "\n",
    "    # 2. compute image gradients in x and y direction\n",
    "    smoothed_img = gaussian_filter(monochrome_img, sigma=1)\n",
    "    Ix, Iy = sobel(smoothed_img, axis=1), sobel(smoothed_img, axis=0)\n",
    "    \n",
    "    # 3. compute gradient patches\n",
    "    grad_patches_x = extract_patches(Ix, keypoints, angles, patch_size)\n",
    "    grad_patches_y = extract_patches(Iy, keypoints, angles, patch_size)\n",
    "    \n",
    "    # 4. compute gradient orientation and magnitude for each pixel in each patch (implement the functions above and call them)\n",
    "    orientations = # TODO: implement and call get_gradient_orientations\n",
    "    magnitudes = # TODO: implement and call get_gradient_magnitudes\n",
    "    \n",
    "    # 5. magnitude weighting (implement the function above and call it)\n",
    "    # TODO: implement and call add_magnitude_weighting\n",
    "    \n",
    "    # 6. create features (implement the function above and call it)\n",
    "    # TODO: implement and call create_features\n",
    "    \n",
    "    # 7. return features\n",
    "    # TODO: implement it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = {}\n",
    "features[\"view1\"] = create_descriptors(view1, *keypoints['view1'], patch_size=16, cell_size=4)\n",
    "features[\"view2\"] = create_descriptors(view2, *keypoints['view2'], patch_size=16, cell_size=4)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Feature Matching\n",
    "\n",
    "Now that we have found image key points and assigned them a (hopefully) unique descriptor, it is time to find corresponding points in both images. There are sophisticated matching algorithms which avoid checking every key point from one view with the one from another view. However, for the sake of this exercise comparing all key points is our way to go. As we encoded the appearance of each key point by a feature vector, we can easily compute its sum of squared differences to all other feature vectors.\n",
    "\n",
    "Hence, your first task is to compute the sum of squared differences between all key points in both views.\n",
    "The result should be a matrix $D$ where $d_{ij} = \\| f_i-f_j\\|_2^2$ and $f_i, f_j$ are the feature vectors from key points i and j. Note that key points i come from view 1 and key points j from view 2.\n",
    "\n",
    "$D$ can be displayed by matplotlib which shows you how close the points are to others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_view1 = features[\"view1\"]\n",
    "features_view2 = features[\"view2\"]\n",
    "\n",
    "def get_distances(features_view1: npt.ArrayLike, features_view2: npt.ArrayLike) -> np.ndarray:\n",
    "    \"\"\"Compute the distance matrix between the feature vectors.\n",
    "\n",
    "    :param features_view1: feature vectors of view 1 [number of keypoints view 1, feature vector size]\n",
    "    :param features_view2: feature vectors of view 2 [number of keypoints view 2, feature vector size]\n",
    "    :return array with distance matrix [number of keypoints 1, number of keypoints view 2]\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "D = # TODO: implement and call get_distances\n",
    "\n",
    "plt.imshow(D)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there is a lot of ambiguity for some points. Hence, it may not be enough to choose the best pair of points but only pairs where the second best match is significantly worse. A simple method to do so is a ratio test. We check if the ratio between the shortest distance $d_1$ of the best match and the distance $d_2$ of the second best match is smaller than some threshold $T$. If so, we return the matched pair otherwise no match is returned.\n",
    "\n",
    "Using this as a check, you can implement a simple matching function. You can verify the outputs by using the visualization function provided below. The function you are about to implement should also return the distance of the points in each match. You may sort the list of matches to see the best ones on top.\n",
    "\n",
    "Hint: Only because key point P1 from view1 matches best with P2 from view2 doesn't mean P2 matches best with P1. The matching relation is not symmetric. So, key points might be used in different matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import ConnectionPatch\n",
    "\n",
    "# use this function to visualize the matches\n",
    "def plot_matches(img1, img2, keypoints1, keypoints2, matches):\n",
    "    fig, axes = plt.subplots(1,2,figsize=(16,7))\n",
    "    \n",
    "    # draw images\n",
    "    axes[0].imshow(img1)\n",
    "    axes[1].imshow(img2)\n",
    "    \n",
    "    # draw matches\n",
    "    for index_1, index_2 in matches:\n",
    "        kp1, kp2 = keypoints1[index_1], keypoints2[index_2]\n",
    "        con = ConnectionPatch(xyA=kp1, coordsA=axes[0].transData,\n",
    "                              xyB=kp2, coordsB=axes[1].transData, color='r')\n",
    "        fig.add_artist(con)\n",
    "        axes[0].plot(*kp1, color='r', marker='x')\n",
    "        axes[1].plot(*kp2, color='r', marker='x')\n",
    "    \n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_matches(distances: npt.ArrayLike, threshold: float) -> np.ndarray:\n",
    "    \"\"\"Get matching keypoints.\n",
    "\n",
    "    :param distances: the keypoint feature vector distance matrix [number of keypoints view 1, feature vector size]\n",
    "    :param threshold: distance < threshold -> match candidate \n",
    "    :return array with matches, array with distances of the matches\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def remove_duplicate_matches(matches: npt.ArrayLike, match_distances: npt.ArrayLike) -> np.ndarray:\n",
    "    \"\"\"Remove duplicate matches.\n",
    "\n",
    "    :param matches: matched keypoints\n",
    "    :param match_distances: distances between matched keypoints\n",
    "    :return array with distinct matches, array with corresponding distances of the distinct matches\n",
    "    \"\"\"\n",
    "    # TODO: implement it\n",
    "\n",
    "def match_keypoints(features_view1, features_view2, threshold=0.7):\n",
    "    \"\"\"\n",
    "    :param features_view1: a 2D numpy array containing the feature vectors for each keypoint in view 1\n",
    "    :param features_view2: a 2D numpy array containing the feature vectors for each keypoint in view 2\n",
    "    :param threshold: the ratio threshold\n",
    "    :return: Two arrays are returned. First, a 2D numpy array where each row \n",
    "             consists of two indices forming a match. The first index corresponds \n",
    "             to to the row number in features_view1 and the second to \n",
    "             the row number in features_view2. The second array is the distance between\n",
    "             the points of a match.\n",
    "    \"\"\"\n",
    "    # 1. compute distances\n",
    "    # TODO: implement it\n",
    "    \n",
    "    # 2. retrieve best matches for key points from view 1 and 2 and perform ratio checks \n",
    "    # TODO: implement the function get_matches above and call it\n",
    "    \n",
    "    # 3. remove duplicates\n",
    "    # TODO: implement the function remove_duplicate_matches above and call it\n",
    "    \n",
    "    # 4. return all matches between both views and their distances\n",
    "    # TODO: implement it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matches, match_distances = match_keypoints(features[\"view1\"], features[\"view2\"])\n",
    "\n",
    "# sort matches by distance\n",
    "sorted_matches = # TODO: implement it\n",
    "\n",
    "# choose k best matches\n",
    "k = # TODO: implement it\n",
    "top_k_matches = # TODO: implement it\n",
    "\n",
    "# visualize\n",
    "plot_matches(view1, view2, keypoints[\"view1\"][0], keypoints[\"view2\"][0], top_k_matches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}