# Efficiently selecting spatially distributed weighted points

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:
To improve the results my professor suggested selecting only the top X% of the keypoints, as there are probably many bad quality keypoints (keypoints have weights or similarity scores by which we could sort them). However, that selection could introduce a strong bias, as, probably, one side of the image has much better keypoints than the other side. So, we would want them to be distributed evenly on the image.

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

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 asked a separate question about tilted rectified images: Skewed rectified aerial images

Code:

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,
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],
: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.')

max_count = count + count_delta

points,
image_shape=image_shape,
count=max_count,

for _ in range(max_iterations_count):

if selected_points_count == count:

if count < selected_points_count <= max_count:
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.')

*,
grid_resolution: float,
image_shape: Tuple[int, int],
count: int,
"""
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)

fill_value=False)

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

break

center=point_grid_index,

*,
center: Tuple[int, int],
"""
Returns 2d array with applied a disc shaped mask over it.
For more details see:
:param array_shape: shape of original image
:param center: center 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 , 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))
image_shape=image_shape,
count=8)
plt.scatter(points[:, 0],
points[:, 1])
color='r')


Output:

Example with (Shi-Tomasi corner detector):

Taking the following image:
duck.jpg

import cv2

import suppression_via_disc_covering as sdc

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),
color=0,
thickness=-1)

image_shape=image.shape[::-1],
count=25)

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

for corner in good_corners:
cv2.circle(image,
center=tuple(corner),
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".

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

• Well written question and clean, professional code, no suggestions re. code style. Have you looked into SIFT/SURF image features, if using them for point detection might yield better, more resilient keypoints? Or using K-means (or another kind of) clustering to select spatially distributed "centroid" keypoints? May 15 '18 at 14:51
• Thanks, @scnerd! I did not consider using SIFT/SURF because I was under impression that ORB yields better results after reading tutorial. But quick googling led me to the following article, and there I can see that ORB wins only by time. I will try them out. Regarding clustering, I got the same idea from a data scientist friend. Sounds like something worth investigating. May 15 '18 at 15:25
• I have recently published a paper that tackles the problem of homogeneous keypoint distribution on the image. C++, Python, and Matlab interfaces are provided in this repository. Moreover, comparison to SDC algorithm is in the paper. Oct 19 '18 at 7:12
• @AlexBailo Thanks for sharing! I took a look at Python code and I see that it could use some refactoring. Would you like to post it here so me or other people review it? :) Oct 19 '18 at 9:17
• Your keypoints have coordinates and some quality metric, let's call those x_i, y_i and Q_i. Try to fit a function f(x, y) = ax + by + c to the keypoints' quality metric. Then create a corrected quality metric Q'_i = Q_i / f(x_i, y_i). Select the best keypoints from the correct quality metric to do your visual tracking. If it works, then it's a faster algorithm than farthest points, and possibly it selects a better set of keypoints (it doesn't constrain the selection based on position relative to nearby keypoints). Oct 20 '18 at 9:08

I've read this question and the associated code a couple of times because I really wanted to review it. First of all I'd say that there's not much to review because it's really well written (and it's a complex subject).

Usage of functools.partial

This may be opiniated, but I don't think you should be using partial here.

I think you used it to make the code cleaner, but all in all, what's the difference between these two pieces of code :

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

for _ in range(max_iterations_count):


and

for _ in range(max_iterations_count):
image_shape=image_shape,
count=max_count,
grid_resolution=grid_resolution)


There are two other reasons I don't think you should use it :

1. You only re-use the partial function once

grid_shape = (int(image_shape[0] // grid_resolution) + 1,
int(image_shape[1] // grid_resolution) + 1)


You don't need to use the int conversion here, integer division // will return an int anyways.

That's a small performance improvement, but you compute radius*radius often while it could be computed once. You could create a radius_pow_2 variable and pass this to your circular_mask function.

In the select function, I'd be inclined to rename count to k. This might not be a popular decision, but you use k everywhere when you explain the algorithm, so it's very clear what it's supposed to do. It's also a pretty popular parameter (think K-Means Clustering or K Nearest Neighbours).

I also think you should revisit the documentation for this parameter : :param max_iterations_count: prevents infinite loop. The idea, if I understood correctly, isn't to prevent an infinite loop, but to set a "time limit" where you accept that the algorithm isn't finding a reasonable solution and this difference is pretty important.

This parameter do could also use some love : :param radius: as number of cells where points won't be selected;. It's not very clear what it means (even though it's clear what it does in your post, but the documentation should be clear, otherwise why have it.)

You throw ValueError, first I don't think that's the right exception (I also think the slim choice of exceptions we can throw in Python is... way too slim.)

Second I think it could be more detailed as to why the algorithm didn't find a solution. Was K too large? Was the initial radius too big/small? I'm pretty sure that by analyzing the responses your algorithm gave while iterating, you could give a little more "meat" to your exception message. While that might not be a good solution for a real time system, it could be an interesting addition for debugging.

If I understood correctly, in the erase_extra_points method, you basically delete every points after the k-est one. This as a consequence that the points near the bottom of your image would be deleted (again, if I understood correctly), without concern towards the importance of the said points. Even if I'm mistaken in my previous sentence, the idea is that deleting points with such a "simple" algorithm could hurt your performance.

Finally, I'm no CV expert either, but if there's one thing I've learned is that if you have a performance of 45% and you need it to be much higher, it might be wise to get more creative and think of other solutions.

• Thank you for the review! You hit some good points. A comment on usage of partial: I used it to show that all the arguments of selected_points_mask will be fixed except the grid_resolution which will be changing in the loop. While this was clear to me, I actually never saw anyone else using partial in a similar manner. So yes, it will be better to remove it following the principle of least astonishment. Regarding erase_extra_points : looks like you are right and this function should be redesigned. Good catch! Aug 5 '19 at 9:42
• About max_iterations_count: IIRC, the algorithm indeed hit infinite loops for some input data. So, providing more info in that ValueError message is important to explain what really happened. And a comment on your final remark: In our research we considered several ways of counting trees on aerial images. One of them was using stereovision which proved to be slow and did not always give good results. The alternative using color-based segmentation was a more appropriate choice for this task. As a reference, here is a link to our work: mdpi.com/2072-4292/11/3/316 Aug 5 '19 at 10:05
• @Georgy This is some very cool work I'll take a look at it Aug 5 '19 at 13:39