trjhtr

Background:

Motivation behind writing the following code is originated in the area of computer vision. More specifically – image rectification. In order to obtain rectified images, one has to find a set of matching features/keypoints on both images beforehand. My code is supposed to operate on these points.

When calculating rectifications for a series of aerial images, I was obtaining poor results. Around 80% of rectified images were squeezed and tilted at a large angle. For example:

My professor suggested that this problem could be solved by choosing keypoints that lay far apart from each other.

When looking for an algorithm that could select the most scattered points, I came across two posts on Stack Exchange:

• 
  Farthest point algorithm in Python¹


• 
  Selecting most scattered points from a set of points

In the second post I found a link leading to the following article: Efficiently selecting spatially distributed keypoints
for visual tracking. There they describe an algorithm that takes into account weights of the points. Since it's possible to sort the keypoints by a score of similarity, this looked like a way to go. I implemented their Suppression via Disk Covering (SDC) algorithm and got ~70% good quality rectified images comparing to previous ~20%. Here is an example of a good one:

Idea of the algorithm is pictured on the following image taken from the article:

Relevant steps (citing the article):

1) Input image for which a set of, for example, k = 20 strong, well-distributed keypoints is sought.

2) Keypoints found by detector.

3) The first (strongest) keypoint is selected and all cells within the approximated radius r are covered.

4) Strongest uncovered point is selected, surrounding cells covered.

5) Finally, with this radius, five keypoints have been selected. This is below the desired k, so a new iteration (6)) is started with a smaller r.

7) More than k points are selected and still uncovered points left:
r is too small, the iteration can be aborted and the next iteration (8)) started.

9) Finally, with this r , exactly k keypoints have been selected and are returned as result

What I want reviewed:

• Code organization, separating logic to functions, DRY, style, choice of names, etc.


• Performance. (It's possible that in the future this code will run in real time. Right now it's not critical, so I didn't do any profiling. But by taking a look on the current state of the code, I don't see how it could be sped up.)


• Missed bugs, edge-cases.


• Alternative algorithms.


• Is this an XY problem? I am new to computer vision and I'm afraid that I could take a wrong way to tackle this problem. I am thinking to ask a separate question about tilted rectified images on Signal Processing Stack Exchange as well.
  Edit: This is done now: Skewed rectified aerial images
  

Code:

suppression_via_disk_covering.py

from functools import partial
from typing import Tuple

import numpy as np


def select(points: np.ndarray,
 *,
 image_shape: Tuple[int, int],
 count: int,
 count_delta: int = 1,
 radius: int = 10,
 radius_delta: int = 2,
 max_iterations_count: int = 15,
 min_cell_size: int = 2,
 max_cell_size: int = 100) -> np.ndarray:
 """
 Selects points by a Suppression via Disk Covering algorithm.
 For more details see:
 http://lucafoschini.com/papers/Efficiently_ICIP11.pdf
 :param points: original set that should be ordered by distance
 :param image_shape: shape of an image
 :param count: number of output points
 :param count_delta: let k = `count` and Δk = `count_delta`,
 if number of found points is within [k; k + Δk],
 return top-k points
 :param radius: initial radius of area where points will be removed
 :param radius_delta: determines width of cells
 :param max_iterations_count: prevents infinite loop
 :param min_cell_size:
 :param max_cell_size:
 :return: mask array with selected strong scattered keypoints
 """
 if len(points) < count:
 raise ValueError('Not enough points to select.')

 grid_resolution = radius_delta * radius / np.sqrt(2)

 max_count = count + count_delta

 points_mask = partial(selected_points_mask,
 points,
 image_shape=image_shape,
 count=max_count,
 radius=radius)

 for _ in range(max_iterations_count):
 result_mask = points_mask(grid_resolution=grid_resolution)

 selected_points_count = result_mask.sum()

 if selected_points_count == count:
 return result_mask

 if count < selected_points_count <= max_count:
 return erase_extra_points(result_mask,
 count=count)

 if selected_points_count < count:
 max_cell_size = grid_resolution
 grid_resolution -= (grid_resolution - min_cell_size) / 2
 else:
 min_cell_size = grid_resolution
 grid_resolution += (max_cell_size - grid_resolution) / 2

 raise ValueError('Number of iterations exceeded.')


def selected_points_mask(points: np.ndarray,
 *,
 grid_resolution: float,
 image_shape: Tuple[int, int],
 count: int,
 radius: int) -> np.ndarray:
 """
 Calculates boolean mask corresponding to array of input points.
 True values are for those points
 that will be selected as scattered enough from each other.
 In case if there were too many points found,
 the mask still will be returned.
 :param points: input array
 :param grid_resolution: size of a cell in a grid
 :param image_shape:
 :param count: number of points to select
 :param radius: as number of cells where points won't be selected
 :return: boolean array with True values for selected points
 """
 points_grid_indices = (points // grid_resolution).astype(int)

 grid_shape = (int(image_shape[0] // grid_resolution) + 1,
 int(image_shape[1] // grid_resolution) + 1)
 grid = np.full(shape=grid_shape,
 fill_value=False)

 result_mask = np.full(shape=points.shape[0],
 fill_value=False)

 for index, point_grid_index in enumerate(points_grid_indices):
 if grid[tuple(point_grid_index)]:
 continue

 result_mask[index] = True

 if result_mask.sum() > count:
 break

 mask = circular_mask(grid.shape,
 center=point_grid_index,
 radius=radius)
 grid[mask] = True

 return result_mask


def circular_mask(array_shape: Tuple[int, int],
 *,
 center: Tuple[int, int],
 radius: int) -> np.ndarray:
 """
 Returns 2d array with applied a disc shaped mask over it.
 For more details see:
 https://stackoverflow.com/questions/8647024/how-to-apply-a-disc-shaped-mask-to-a-numpy-array
 https://stackoverflow.com/questions/44865023/circular-masking-an-image-in-python-using-numpy-arrays
 :param array_shape: shape of original image
 :param center: center of the disc
 :param radius: radius of the disc
 :return: boolean array with applied circular mask
 """
 y, x = np.ogrid[-center[0]:array_shape[0] - center[0],
 -center[1]:array_shape[1] - center[1]]
 return x * x + y * y <= radius * radius


def erase_extra_points(array: np.ndarray,
 *,
 count: int) -> np.ndarray:
 """
 Let n = `count`, sets to False all elements
 after the n-th occurrence of a True element.
 :param array: input boolean array
 :param count: number of True elements to remain
 :return:
 """
 array = array.copy()
 last_true_index_to_remain = np.where(array)[0][count]
 array[last_true_index_to_remain:] = False
 return array

Examples of usage:

Simple example without opencv, ignoring ordering of points by strength:

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

import suppression_via_disc_covering as sdc

image_shape = (1500, 2000)
points = np.random.uniform(low=(0, 0), 
 high=image_shape, 
 size=(100, 2))
mask = sdc.select(points,
 image_shape=image_shape,
 count=8)
plt.scatter(points[:, 0],
 points[:, 1])
plt.scatter(points[mask, 0],
 points[mask, 1],
 color='r')

Output:

Example with opencv (Shi-Tomasi corner detector):

Taking the following image:
duck.jpg

import cv2

import suppression_via_disc_covering as sdc

image = cv2.imread('duck.jpg')
image = cv2.cvtColor(image,
 cv2.COLOR_BGR2GRAY)

corners = cv2.goodFeaturesToTrack(image,
 maxCorners=1000,
 qualityLevel=0.01,
 minDistance=10)
corners = corners.reshape(-1, 2)

for corner in corners:
 cv2.circle(image,
 center=tuple(corner),
 radius=3,
 color=0,
 thickness=-1)

mask = sdc.select(corners,
 image_shape=image.shape[::-1],
 count=25)

good_corners = corners[mask, :]

image = cv2.cvtColor(image,
 cv2.COLOR_GRAY2RGB)

for corner in good_corners:
 cv2.circle(image,
 center=tuple(corner),
 radius=3,
 color=(0, 0, 255),
 thickness=-1)

cv2.imshow('image',
 image)
cv2.waitKey(0)

Output:

P.S.: I'm not including an example with pairs of aerial images, since there is much more code there that should be reviewed separately. Moreover, I still consider that 70% of good results is not good enough and hence this example would be considered as "doesn't work as intended".

¹I tried this algorithm as well and quality of the results improved from ~20% of good images to ~45%. This is not good enough.