3.0-tl-feature-matching.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3: 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. Still, if it does not fit your idea you may ignore it.\n",
    "\n",
    "\n",
    "## 1. 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 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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 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",
    "### 2.1. Key Point Detection using Harris Corner Detector\n",
    "\n",
    "Complete the implementation below. 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.\n"
   ]
  },
  {
   "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",
    "# use this function to turn the original image into grayscale\n",
    "def to_grayscale(img):\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",
    "\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",
    "    \n",
    "    # 1. convert to grayscale\n",
    "\n",
    "    # 2. compute image gradients in x and y direction (smooth the image before doing that)\n",
    "\n",
    "    # 3. compute structure tensor\n",
    "\n",
    "    # 4. compute harris score\n",
    "\n",
    "    # 5. perform non-maximum surpression\n",
    "    \n",
    "    # 6. compute angles"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.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)"
   ]
  },
  {
   "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",
    "    # add padding such that also patches at the borders can be extracted\n",
    "    pad_size =\n",
    "    padded_img =\n",
    "    \n",
    "    # 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": {
    "scrolled": false
   },
   "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)"
   ]
  },
  {
   "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",
    "### 2.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 weigh 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",
    "\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\n",
    "    orientations = \n",
    "    magnitudes = \n",
    "    \n",
    "    # 5. magnitude weighting\n",
    "    \n",
    "    # 6. create features\n",
    "    features = []\n",
    "    for k in range(len(keypoints)):\n",
    "        # compute histograms\n",
    "    \n",
    "        # stack histograms together\n",
    "        feat_vector = ...\n",
    "        \n",
    "        # normalize feature vector\n",
    "    \n",
    "    # 7. return features\n",
    "    "
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 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",
    "D = \n",
    "\n",
    "plt.imshow(D)\n",
    "plt.show()"
   ]
  },
  {
   "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 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",
    "    \n",
    "    # 1. compute distances\n",
    "    \n",
    "    # 2. retrieve best matches for key points from view 1 and 2\n",
    "    \n",
    "    # 3. perform ratio checks\n",
    "    \n",
    "    # 4. remove duplicates\n",
    "    \n",
    "    # 5. return all matches between both views and their distances"
   ]
  },
  {
   "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 = \n",
    "\n",
    "# choose k best matches\n",
    "k = \n",
    "top_k_matches =\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",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}