# Spatial Join Pandas Dataframes of Bounding Boxes (cross match)

Problem Statement Imagine there are two almost identical images with annotations (bounding boxes for certain objects), where one is so-called golden image (template) containing all must-have objects (ground truth), and the other one is the query image (input) let's say with lesser objects.

Given that we already have Pandas dataframe represnetation of objects with bounding box infos for each of the images like, how one can perform a spatial left join on template between bounding boxes (objects) so that we can easily identify the missing objects in the input?

Example. Let's say the template looks like:

and the corresponding template dataframe:

          name  xmin  ymin  xmax  ymax
0     big fish   251   504   485   654
1   small fish   583   572   748   660
2     big fish  1080   484  1236   597
3     big fish   574   122  1076   505
4     big fish  1351   187  1583   369
5   small fish   369    31   506   115
6   small fish  1081   148  1111   190
7   small fish   684   505   732   535
8   small fish   939   521   992   570
9   small fish   417   661   497   705
10  small fish   743   598   792   642
11  small fish   667   657   708   691


And the input image looks like:

and the corresponding input dataframe:

         name  xmin  ymin  xmax  ymax
0  small fish   342    16   478   101
1    big fish   221   490   459   646
2  small fish   579   564   723   641
3    big fish  1342   161  1558   337
4    big fish   557   102  1045   492
5  small fish  1049   132  1087   176
6  small fish   389   652   484   694
7  small fish   914   514   964   556
8  small fish   639   640   688   676


Expected Result In this example, there fishes are not present on the input image (missing), and I would lile to extract and identify that information by cross matching objects between the template and input dataframes. Ideally, I seek to have a subset of template dataframe only containing the missing objects in the input image:

          name  xmin  ymin  xmax  ymax
0     big fish  1080   484  1236   597
1   small fish   684   505   732   535
2   small fish   743   598   792   642


and the overlay on the input image would look like:

Attempts I have tried to use geopandas.GeoDataFrame.sjoin, borrowing functionalities for Points, Polygons for merging the dataframes. When bounding boxes are reasonably distanced from one and other, it would work as expected. However, when are are in proximity of one another, and often even overlaps, then geopandas sjoin and any other merging functionalities wouldn't work.

I have also tried to use distances (centers of bounding boxes) together with IoU (Intersection over Union) to cross match these geometries, but it wouldn't inheritly know how far it should look to cross match, and defining the threshold is not ideal, because we simply wouldn't know how many objects to expect to include or exclude (unless it is harded in the logic and definitely hard to maintain).

Question: Is there a better, smarter and efficient way to accomplish this?

P.S. (Materials): In order to make things easy to contribute and access these images, xmls, and dataframes, I have put everything in a pulbic Github repo, containing also a Notebook with some functions and steps using geopandas.GeoDataFrame.sjoin!

• Very well laid out question: you even did a git repo for the code! Feb 25 at 21:35
• It looks like it could be a classic tracking problem instead of what you are describing "there are two almost identical images with annotations". Do you mean you are looking at two consecutive frames in a movie or do you actually have one reference frame and all you care about is comparing fishes in other frames with that golden frame? Asking because this seems like a tracking problem but you framed as something slightly different. Feb 25 at 23:37
• @Justtrying The latter, "you actually have one reference frame and all you care about is comparing fishes in other frames with that golden frame"! This is not a tracking problem. I would happily adopt apporaches available in solving a tracking problem though, if applicable. Feb 26 at 21:54
• IoU is not going to be useful if there is no overlap between bounding boxes. So if the other frame is taken a few seconds after the golden frame, that will not work. Similarly, if you derive a criterion that works for two frames 1/30 second apart, it might not work consistently with frames 20/30 seconds apart. In any case, the Munkres algorithm I suggest in my answer works for the assignment. Feb 26 at 23:08
• The cost matrix in the Munkre algorithm is up to you to define. Bounding boxes IoU and centroid distances assume there was little movement between frames. Other ways to define the cost matrix would be: differences between the areas of the fishes could work. something like abs(a_i - a_j) / (a_i +a_j) which doesn't rely on bounding boxes overlapping. You could also use feature detectors and come up with a cost for each potential match. Whatever solution you decide to use will be heavily dependent on the nature of the problem, I still feel like the framing is too vague to give you a great answer Feb 26 at 23:13

Edited after comments in other answers: Generally speaking, I would reframe this as a linear sum assignment problem.

This can be solved using a modified version of Munkres algorithm allowing a cost for non-assignment. wich time complexity is pretty bad ($$O(n^3)$$), but will work for a dozen fishes. For, reference, the Matlab version allows you to handle tracks that end and start across frames, i.e fishes that disappear and appear between frames.

To use the Munkres algorithm, you need to define a cost matrix, with $$N_{tracks}$$ rows (first frame) and $$N_{detections}$$ columns (second frame).

The Munkres algorithm will minimize the global assignment cost.

### Case 1: significant overlap of bounding boxes across frames (tracking problem):

For the track $$i$$ in the first frame and detection $$j$$ in the next one, you can define the cost as $$IoU(i, j)$$ which is the intersection over union of the two bounding boxes for track $$i$$ and detection $$j$$. You could also consider using the distance between the centroid of the bounding boxes $$d(i, j)$$ or a combination of the two with a total cost such as $$C(i, j) = IoU(i, j) + \alpha \times d(i, j)$$ with $$\alpha$$ a parameter to determine to tune the respective weight of each cost in the full cost matrix. If you are only using bounding boxes the IoU is pretty easy to compute.

### Case 2: no significant overlap of bounding boxes across frames (detection problem):

In that case, you cannot rely on positional information. But hopefully, the fish shapes remain largely unchanged. So you can build descriptors/features, for instance:

• bounding box area: measure the number of pixels in the bounding box (this assumes the fish also didn't change orientation drastically, as fishes are pretty flat so the area from the side will be very different from the area from the front. You could consider using the longest side of the bounding box to mitigate this
• color composition: create a binned RGB histogram from all the pixels in the fish (ideally you would have access to a finer segmentation than just a bounding box to make it less sensitive to the background color)

You could also use feature descriptors such as SIFT, AKAZE, etc...

But it all comes down to the same two steps:

• find a good way to compare any pair of objects across frames
• make an optimal decision about how to match them across frames and how to decide which are missing

The second part will always be a linear sum assignment problem.

So the only thing now is that the scipy version doesn't offer the option to specify the unassignedTrackCost or unassignedDetectionCost like the Matlab version does. And this is actually what will allow you to handle fishes appearing or disappearing the way the Matlab version does. So you will need to modify it. Looking at the picture below, you now have the costMatrix and you need to build the bigger matrix to be able to handle the cases when fish appear or disappear.

Once you have managed to create the full cost matrix you can solve it using linear_sum_assignment and then find the tracks (resp. detections) that were assigned to dummy detections (resp. tracks).

## Implementation

### Getting the Cost matrix (distance only)

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.spatial.distance import cdist
from scipy.optimize import linear_sum_assignment

def transform(df):
df["centroid_x"] = (df["xmax"] + df["xmin"]) / 2
df["centroid_y"] = (df["ymax"] + df["ymin"]) / 2
return df

df0 = transform(df0)
df1 =  transform(df1)

distance_cost_matrix = cdist(df0[["centroid_x", "centroid_y"]], df1[["centroid_x", "centroid_y"]])
cost_mat = np.log(distance_cost_matrix) #np.log(np.multiply(area_cost_matrix, distance_cost_matrix))


### Modified Munkres algorithm

def pseudo_inf(cost_mat, inf_func):
pseudo_inf_val = inf_func(cost_mat[cost_mat != np.inf])
pseudo_cost_mat = cost_mat.copy()
pseudo_cost_mat[pseudo_cost_mat == np.inf] = pseudo_inf_val
return pseudo_cost_mat, pseudo_inf_val

def get_costs(cost_mat, row_ind, col_ind):
costs = [cost_mat[i, j] for i, j in zip(row_ind, col_ind)]
return costs

def assign_detections_to_tracks(
cost_mat,
cost_of_non_assignment=None
):
# in case there are infinite value, replace them by some pseudo infinite
# values
inf_func =  lambda x: np.max(x) * 2
pseudo_cost_mat, pseudo_inf_val = pseudo_inf(cost_mat, inf_func)
assigned_rows = []
unassigned_rows = []
assigned_cols = []
unassigned_cols = []
full_cost_mat = None
# basic case, handled by linear_sum_assignment directly
if cost_of_non_assignment is None:
assigned_rows, assigned_cols = linear_sum_assignment(pseudo_cost_mat)
assignment_costs = get_costs(cost_mat, assigned_rows, assigned_cols)
# if one cost of non assignment is provided, use it
else:
# build the pseudo-array
top_right_corner = np.full((cost_mat.shape[0], cost_mat.shape[0]), pseudo_inf_val)
np.fill_diagonal(top_right_corner, cost_of_non_assignment)
bottom_left_corner = np.full((cost_mat.shape[1], cost_mat.shape[1]), pseudo_inf_val)
np.fill_diagonal(bottom_left_corner, cost_of_non_assignment)
top = np.concatenate((cost_mat, top_right_corner), axis=1)
zero_corner = np.full(cost_mat.T.shape, 0)
# zero_corner = np.full(cost.shape,cost_of_non_assignment)
bottom = np.concatenate((bottom_left_corner, zero_corner), axis=1)
full_cost_mat = np.concatenate((top, bottom), axis=0)
# apply linear assignment to pseudo array
row_idxs, col_idxs = linear_sum_assignment(full_cost_mat)
# get costs
for row_idx, col_idx in zip(row_idxs, col_idxs):
if row_idx < cost_mat.shape[0] and col_idx < cost_mat.shape[1]:
assigned_rows.append(row_idx)
assigned_cols.append(col_idx)
elif row_idx < cost_mat.shape[0] and col_idx >= cost_mat.shape[1]:
unassigned_rows.append(row_idx)
elif col_idx < cost_mat.shape[1] and row_idx >= cost_mat.shape[0]:
unassigned_cols.append(col_idx)
#
full_costs = get_costs(full_cost_mat, row_idxs, col_idxs)
assignment_costs = get_costs(
full_cost_mat, assigned_rows, assigned_cols)
return assigned_rows, assigned_cols, unassigned_rows, unassigned_cols, full_cost_mat, assignment_costs


### Peforming the assignment:

cost_of_non_assignment was tuned looking at a histogram of cost_mat.ravel()

(
assigned_rows,
assigned_cols,
unassigned_rows,
unassigned_cols,
full_costs,
assignment_cost
) = assign_detections_to_tracks(cost_mat, cost_of_non_assignment=10)
unassigned_rows, unassigned_cols


### Results

df0.iloc[unassigned_rows, :]

          name  xmin  ymin  xmax  ymax centroid_x centroid_y
0     big fish  1080   484  1236   597   1158.0     540.5
1   small fish   684   505   732   535    708.0     520.0
2   small fish   743   598   792   642    767.5     620.0


it works!

• Wouldn't be able to show case your suggestion, especially Case 2 using linear sum assignment problem? Is there is a straightforward implementation for this algorithm? Thanks Feb 27 at 21:34
• Yes I edited again. I wrote that pretty quickly so it might not be perfect but you should have all you need. The threshold I picked for cost_of_non_assignment is by no mean universal it just worked (and was fairly robust) for the example you shared Feb 28 at 2:27
• Also I don't handle the two separate cases where you can specify both cost_of_unassigned_track and cost_of_unassigned_detection. Feb 28 at 2:30
• Wow, very detailed and impressive response, thanks for the contribution. Please give me a day to go over your response and will get back to you shortly. From now it looks like you are getting the bounty!! Feb 28 at 9:06
• I have tested shortly as promised, and it works for a number of tests I have carried out. Many thanks again for the contribution, already accepted it as a solution. Mar 1 at 11:01

One option is a brute-force approach.

Write a custom function that compares a single input box to a single template box and determines if they meet the criteria or not.

Then loop through the template dataframe row-by-row. Compare the box for the current template row to the box in every row in the input dataframe. Use the function to check whether the two boxes meet the criteria or not. If yes, remove the row from the template dataframe. The remaining rows in the template dataframe are the final result.

This algorithm has O(nm) runtime complexity (n is the length of the template dataframe and m is the length of the input dataframe).

• This greedy algorithm is dependent on the order in which you would process the row which is likely undesirable. Unless you are very lucky and are able to find a criterion that would never give you any false positive/ false matches, there will likely be situations where a template box might match two boxes in the query image. It might work for the specific images shared here but for fishes swimming at different depths and in any direction, you will likely end up in facing the situation I discussed. This is a common issue with tracking multiple object. Feb 26 at 1:03
• @Justtrying If the criteria is a real number measure and not just a boolean then the order can be determined as sort order, using trees to only compute overlapping boxes. OP mentioned various real valued measures in the question and more are available so this should not be a concern. Feb 26 at 1:55
• Well maybe but this solution also doesn't provide a global optimum and could lead to incorrect detections. @Brian Spiering says "This algorithm assures correctness" but I don't think that's the case. Feb 26 at 4:08
• Great points! I have revised my answer. Feb 26 at 13:32
• That was my first solution for sometime now, e.g. brute-forcing via Intersection over Union (IoU), or even distanced-based cross-matching. As you said it yourself, first it has O(nm) runtime complexity; not desirable, ans imost importantly is not reliableIt fails if the criteria is not met e.g. bounding box being small or images (and cosequently boudning boxes) have different sizes (translation, rotation, zoom etc.) and a proper image registeration is needed, which is not also easy and cheap! That is why it is not what I am looking for, thansk though! Feb 26 at 22:07

TL;DR: I would focus on a good measure first (described below, but IntOvUn is fine too) and then use KD-Tree in Scipy to speed up the computation if needed: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html

I can't imagine being able to rely on indexes built into Pandas or some other table system or database that would have the flexibility to handle your problem.

IntersectOverUnion seems like a decent measure for bounding boxes. But I think the percent of a given input box that overlaps with the template box is an even better measure since one really only cares about 1) what is the closest template box 2) how wrong the input boxes are given the closest matching template box. Or maybe some weighted combination of the above and IntOvUnion. But that should give you an assignment of each input box to the most overlapping template box, if there is one.

You can just double for loop if your data set is small enough for that not to be an issue. Given that the squared factor is within each image as opposed to over the whole dataset this really shouldn't be too much of a slowdown. Nonetheless, there is are some data structures especially made to see if intervals, rectangles, rectangloids and their higher dimensional analogs overlap in O(log(n)) time. Interval trees are probably the most appropriate for this narrow task. But there are no implementations of multidimensional interval trees, even 2D ones.

So instead we settle for SciPy's KD-Tree. You have to make the careful translation between points and intervals and back. But a KD-Tree can give you log(n) lookups.

I hope that helps.

• v nice answer. i'm looking now at both KD-tree and the polygon overlap now Feb 26 at 1:39
• I am aware of Shapely. In the git repo I have attached, I have tried spatial matching between geometrical polygons, but for difficult problems i.e. images with a mix of small, large and tight bounding boxes, it fails! See the Notebook please. Feb 26 at 21:59
• @TwinPenguins I just ran the notebook locally. Let's ignore Shapely for now. Any "failures" seem to be due to a high density of bounding boxes that causes almost every box in the golden data to be partially overlapping with some box in the input data even though the box in the input data isn't even trying to cover the same fish. So if you look for just a boolean does_overlaps you'll get near 100% overlaps even though most of those overlaps aren't meaningful. You need a non-boolean, soft measure and a threshold to avoid what are effectively false detections of comparable bounding boxes. Feb 26 at 23:10
• @TwinPenguins I would advise to just avoid the optimization and complexity with geopandas sjoins, Shapely, or KD-Trees. Your data is small. Do a simple double for loop with a soft measure and a threshold and get that to behave how you want before deciding things about indexing and data structures because the latter two will be determined by what you end up happy with. Feb 26 at 23:13