Revision d9e283f480a2592c4617102783c60d297754188b (click the page title to view the current version)

Corner Detection

Changes from d9e283f480a2592c4617102783c60d297754188b to f7bdeeefccdadbcf2754a6d5b466362d57ba9f8f

---
title: Corner Detection
categories: session
---


**Reading** Ma 4.3 and 4.A

**Warning** The textbook starts Chapter 4 by discussing *tracking*,
which means that motion as a function of time is considered as well
as the image as a funciton of spatial co-ordinates.
This is a lot of concepts and quantities to process at the same time.
We will instead start by discussing features in a still image.
When we have a good idea of what features are and how they behave,
we shall introduce motion.

**Briefing** [Corner Lecture]()

# Exercises

## Learning Objectives

1.  What makes a feature in visual terms?
2.  What makes a feature in mathematical terms?
3.  How do we differentiate a sampled signal?
4.  How does the Harris corner detector work?

## Setup

We will use scipy for a small part of the exercises,
if you haven't already, run 'pip install scipy'.

We need an image to work with, you can either load an image from disk or
capture a new image from the webcam with the below code.

```python
# Imports
import cv2 as cv
import numpy as np

# Capture a single frame
vid = cv.VideoCapture(0)
image = None

# Capture frames until we click the space button, use that image
while True:
    _, frame = vid.read()

    # Display the resulting frame
    cv.imshow("frame", frame)

    # Click space to capture frame
    k = cv.waitKey(1)
    if k % 256 == 32:
        # SPACE pressed
        image = frame
        break

# After the loop release the cap object
vid.release()
# Destroy all the windows
cv.destroyAllWindows()
```

As we will be working with a lot of different images in this exercise,
it is recommended to save it to disk, we can do that with

```python
# Imports
from pathlib import Path

# First we create a path to an images directory
p = Path.cwd() / "images" # <--- current working directory + /images
if not p.is_dir():
    p.mkdir()

# Then we save the image to the directory with name "frame.jpg"
cv.imwrite(str(p / "frame.jpg"), image)
```

Now we convert the image to gray-scale
We will be using opencv for today's exercises, you can install it with `pip install opencv-python`.
We might also use scipy, which can be downloaded the same way (`pip install scipy`)

```python
# Convert frame to grayscale
image_gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
We will also be working on a image, e.g. the valve image from https://upload.wikimedia.org/wikipedia/commons/f/f0/Valve_original_%281%29.PNG .
Load the image and convert it to grayscale with:

# Save gray to images
cv.imwrite(str((p / "gray.jpg")), image_gray)
```
  ```python
  import cv2 as cv

## Exercise 1, Sobel/Derivative filter
  # Load image, replace "path" with the image path
  img = cv.imread("path", 1)
  # Convert img to BGR, then to grayscale
  img = cv.cvtColor(img, cv.COLOR_RGB2BGR)
  img_gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
  ```

The first exercise is to implement a Sobel-filter.
Recall from the theory that we need to implement two 3x3 kernels to convolve with the original image.
## Exercise 1

This can be done using scipy.signal.convolve2d ( https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html )
Learning goal: 1D derivatives and 1D convolutions

(Note: For a larger challenge you can also try implementing your own algorithm for the convolve function using numpy)
**Part 1**<br>
Extract one row from the grayscale image, and visualize it (e.g. with matplotlib.pyplot).<br>
Does the values correspond to what you would expect from the row?

Code answers from here on will be collapsed, we recommend that you try to
implement them yourself before reading an answer.
**Part 2**<br>
Convolve the row with a $[1/2,-1/2]$ kernel, using numpy, and visualize.<br>
Does the values make sense?

<details>
  <summary>Hint (Click to expand)</summary>
  <summary>Hint 1 (Click to expand)</summary>
  </br>
  Use `scipy.signal.convolve2d(<image>, <filter>, boundary='symm', mode='same')`
</details>
  Create a new 1D array with `np.zeros(<shape>)` and iterate over the row and the kernel.
  Remember that the resulting array should be smaller than the original.
  </details>
&nbsp;

<details>
  <summary>Solution (Click to expand)</summary>

  <summary>Hint 2 (Click to expand)</summary>
  </br>
  If you are not able to get a result using numpy, use scipy.signal (can also be used to compare your result):
  ```python
  from scipy import signal

  sobel_x = np.array([[-1, 0, 1],
                      [-2, 0, 2],
                      [-1, 0, 1]])
  sobel_y = np.array([[-1, -2, -1],
                      [0, 0, 0],
                      [1, 2, 1]])

  grad_x = signal.convolve2d(image_gray, sobel_x, boundary='symm', mode='same')
  grad_y = signal.convolve2d(image_gray, sobel_y, boundary='symm', mode='same')
  row_d = signal.convolve(row, kernel)
  # or
  # row_d = signal.correlate(row, cv.flip(kernel, -1))
  ```
</details>
&nbsp;

You should then show the images using cv.imshow or save using cv.imwrite, as we did earlier.
Discuss the results.

You can compare the results of your implementation with the built in function
`cv.Sobel`.

```python
grad_x_cv = cv.Sobel(image_gray, cv.CV_32F, 1, 0, ksize=3)  # gradient along x axis,
grad_y_cv = cv.Sobel(image_gray, cv.CV_32F, 0, 1, ksize=3)  # gradient along y axis,
```

Using `cv.Sobel`, adjust the kernel-size to see how that changes the result.

Compute the magnitude and orientation of the derivatives using
```python
magnitude = np.sqrt((grad_x ** 2) + (grad_y ** 2))
orientation = np.arctan2(grad_y, grad_x) * (180 / np.pi) % 180
```

Show/Save the images, and discuss.

## Exercise 2, Eigenvalues and eigenvectors

TODO: ?
<details>
  <summary>Hint 3 (Click to expand)</summary>
  </br>
  If you use cross-correlation instead of convolution, flip the kernel.
</details>
&nbsp;

As this is not a trivial step, we will calculate the EigenVals and Vecs using
`cv.cornerEigenValsAndVecs`, note that this function also applies a sobel-filter
so we will call this function using the grayscale image.
**(Optional, if you have time) Part 3**<br>
Repeat part 2 with an image column (instead of row) and the transpose of the kernel.

```python
bsize = 3  # the size of neighbourhood considered for corner detection
ksize = 3  # Aperture parameter of the Sobel derivative used.
eigenv = cv.cornerEigenValsAndVecs(image_gray, bsize, ksize)
```
## Exercise 2

The cornerEigenValsAndVecs algorithm will
Learning goals: Blur filters

1. Calculate the derivatives $dI/dx$ and $dI/dy$ using the sobel filter (as we did in exercise 1)
2. For each pixel $p$, take a neighborhood $S(p)$ of blockSize `bsize*bsize`,
Calculate the covariation matrix $M$ of the derivatives over the neighborhood as:
**Part 1**<br>

$$M = \begin{bmatrix}
        \sum_{S(p)}(dI/dx)^2 & \sum_{S(p)}(dI/dx)(dI/dy) \\
        \sum_{S(p)}(dI/dx)(dI/dy) & \sum_{S(p)}(dI/dy)^2
$$1/16 \begin{bmatrix}
        1 & 2 & 1 \\
        2 & 4 & 2 \\
        1 & 2 & 1 \\
      \end{bmatrix}
$$

After this the eigenvectors and eigenvalues of $M$ are calculated and returned.

This should result in a $h*w*6$ array with values ($\lambda_1,\lambda_2,x_1,y_1,x_2,y_2$) where

* $\lambda_1,\lambda_2$ are the non-sorted eigenvalues of $M$
* $x_1,y_1$ are the eigenvectors corresponding to $\lambda_1$
* $x_2,y_2$ are the eigenvectors corresponding to $\lambda_2$

Analyze? TODO: How can students learn from this exercise?

## Exercise 3, Harris-criterion

Recall from the theory that we can calculate the harris-criterion with:
$$\lambda_1 * \lambda_2 - k * (\lambda_1 + \lambda_2)^2$$

Create a new image of the same size as image_gray, and calculate the
harris criterion for all pixels.

<details>
  <summary>Solution (Click to expand)</summary>
Works as an approximation of a $3x3$ gaussian blur filter.<br>
Using `scipy.signal.convolve2d`, `scipy.signal.correlate2d` or `cv.filter2d`, apply the filter
to the entire grayscale image and either show the image or write it to file with `cv.imwrite`.
&nbsp;
How does the filter affect the image?

  A simple (naive) solution, is to create an empty image and
   iterate over all pixels like this:
**Part 2**<br>
Compare the above result with the result from `cv.GaussianBlur(img_gray, (3, 3), -1)`

   (Note that this is a very computationally expensive solution,
   and takes a long time to complete)
**(Optional) Part 3**<br>
Using the function from part 2, test out different kernel sizes and compare the difference.

  ```python
  # Create an empty image of size image_gray.shape
  c_x = np.zeros(image_gray.shape, dtype=np.float32)
## Exercise 3

  # Define k, e.g. 0.04
  k = 0.04  # Harris free parameter
Learning goals: 2D Derivatives

  # Iterate over image and set values
  for i in range(c_x.shape[0]):
    for j in range(c_x.shape[1]):
      lambda_1 = eigenv[i, j, 0]
      lambda_2 = eigenv[i, j, 1]
      c_x[i, j] = lambda_1 * lambda_2 - k * np.power((lambda_1 + lambda_2), 2)
  ```
**Part 1**<br>
Apply the sobel operator
$$G_x = 1/8 \begin{bmatrix}
        1 & 0 & -1 \\
        2 & 0 & -2 \\
        1 & 0 & -1 \\
      \end{bmatrix}
$$

  Another, faster solution:
as you did in exercise 2.1.<br>
How does this relate to the 1D derivative you did in exercise 1.2?

  ```python
  # Define k, e.g. 0.04
  k = 0.04  # Harris free parameter
**Part 2**<br>
Repeat part 1 with $G_y$
$$G_y = 1/8 \begin{bmatrix}
        1 & 2 & 1 \\
        0 & 0 & 0 \\
        -1 & -2 & -1 \\
      \end{bmatrix}
$$

  # Extract eigen-values
  lambdas_1 = eigenv[:,:,0]
  lambdas_2 = eigenv[:,:,1]
**Part 3**<br>
Compute and visualize the gradient magnitude of the image, using the results from part 1 and 2.

  # Compute harris-criterion
  c_x = lambdas_1 * lambdas_2 - k * np.power((lambdas_1 + lambdas_2), 2)
  ```
**(Optional, only if you have time) Part 4**<br>

</details>
&nbsp;
Using the same method as in exercise 1.2 and 1.3, compute the gradient of all
rows and columns, and compute the magnitude (as in exercise 3.3), compare with the
magnitude from 3.3.

Take a look at some of the values you calculated above and compare them to the
original image (specifically where you expect the algorithm to detect corners).
You can slice the image with the `c_x[:,:]` command, the syntax is
[starting_row(including):ending_row(not including),starting_column(including):ending_column(not including)]
, e.g. if you want to look at the top left corner around ~10% from the edge on a
480x480 image, you can use `c_x[40:60, 40:60]` or for the same but the bottom right
corner, use e.g. `c_x[-60:-40,-60:-40]`.
## Exercise 4

Learning goals: Introduction to harris corner detector

## Exercise 4, Simple corner-detector
**Part 1**<br>

For the last exercise for corner-detection, you should draw a circle around all the
corners you found on the original image/frame.
Consider the grayscale image we have been working with, and the gradient magnitude from 3.3.<br>
Where do you expect to find corners?

Start by picking a threshold T based on the information you found in the last exercise
(or by using `T = min_cx + (max_cx - min_cx) * 50 / 100)` ). Now draw a circle around
all pixels exceeding this value, using `cv.circle`.
Apply opencv's built in harris-detector.<br>
E.g. `cv.cornerHarris(img_gray, block_size, kernel_size, k)` with `block_size = 2`, `kernel_size = 5` and `k = 0.06`.<br>
Here, block_size is the size of neighbourhood considered for corner detection, kernel_size is the size of the sobel derivative kernel, while k is the harris free parameter.

Adjust the threshold T and try with different values to see how this changes the detection.
Make a copy of the original image (with colors) and make circles around any corners found by the harris-detector.<br>
Example code for drawing circles is added below.

<details>
  <summary>Solution (Click to expand)</summary>

  A simple solution, is to iterate over all pixels and check compare the value against
  the threshold.

  ```python
  min_cx = np.min(c_x)
  max_cx = np.max(c_x)
  T = min_cx + (max_cx - min_cx) * 50 / 100)
  f_circles = frame

  # Drawing a circle around corners (naive, slow)
  for i in range(f_circles.shape[0]):
    for j in range(f_circles.shape[1]):
      if c_x[i, j] > T:
        cv.circle(f_circles, (j, i), 5, (0, 0, 255), 2)
  ```

  Another solution is to create a list of all pixels above the threshold before
  iterating over them:

  <summary>Code</summary>
  </br>
  ```python
  min_cx = np.min(c_x)
  max_cx = np.max(c_x)
  T = min_cx + (max_cx - min_cx) * 50 / 100)
  f_circles = frame
  cx = cv.cornerHarris(img_gray, bsize, ksize, k)

  mask = c_x > T
  corners = np.argwhere(mask).tolist()
  T = 0.1 # Threshold
  c_image = img

  for corner in corners:
    cv.circle(f_circles, tuple(corner), 5, (0, 0, 255), 2)
  for i in range(c_image.shape[0]):
      for j in range(c_image.shape[1]):
          if c_x[i, j] > T:
              cv.circle(c_image, (j, i), 2, (0, 0, 255), 2)
  ```

</details>
&nbsp;

Now do the same to an the original image based on the result from `cv.cornerHarris` and compare with your own result.

Using `cv.cornerHarris`, adjust block_size, aperture_size and the harris parameter (k),
to see how they affect the result.

## Exercise 5, SIFT

## Exercise 6, Feature matching

Save/visualize the result, how does it compare with your expectations?

# Additional Reading
**Part 2**<br>
Adjust the threshold `T` when drawing circles, what does this do?

+ OpenCV/Python Tutorial
    - Background: [Understanding Features](https://docs.opencv.org/master/df/d54/tutorial_py_features_meaning.html)
    - [Harris Corner Detection](https://docs.opencv.org/master/dc/d0d/tutorial_py_features_harris.html)
    - Overview
      [Feature Detection and Description](https://docs.opencv.org/master/db/d27/tutorial_py_table_of_contents_feature2d.html)
+ Ma 2004 Chapter 11.1-2.
**(Optional) Part 3**<br>
Adjust the kernel_size (must be positive and odd), block_size and/or k, and observe how they change the result.