# Optimizing the check for overlapping features
---

### Overview
Currently checking for overlap by checking if each feature's bounding box overlaps with every other feature's bounding box. Although we are only checking unique combinations/pairs of bounding boxes, the algorithm is still $O(n^2)$ given that the number of combinations scales with

\begin{equation}
{}_n C_2 = \frac{n!}{2(n-2)!}
\end{equation}

Here we'll look at implementing a more efficient method, cell listing [1] ($O(n)$). Basically involves subdividing the spatial domain into cells (wider than the threshold distance) and then checking for overlap among only the features in each cell (and the neighboring cells).

[1] https://en.wikipedia.org/wiki/Cell_lists


```python
from itertools import combinations

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

%matplotlib inline
```

### Check if two boxes overlap
(This part is fine I think)


```python
def bboxes_overlap(bbox_1, bbox_2):
    """Determines if two bounding boxes overlap or coincide

    Parameters
    ----------
    bbox_1 : array-like (or 4-tuple)
        1st bounding box
        convention: (x_min, y_min, x_max, y_max)
    bbox_2 : array-like (or 4-tuple)
        2nd bounding box
        convention: (x_min, y_min, x_max, y_max)

    Returns
    -------
    overlap : bool
        True if bounding boxes overlap / coincide
        False otherwise

    References
    ----------
    [1] https://stackoverflow.com/a/20925869/5285918
    """
    # 2 tiles overlap iff their projections onto both x and y axis overlap
    # Overlap in 1D iff box1_max > box2_min AND box1_min < box2_max
    overlap = ((bbox_1[2] >= bbox_2[0]) & (bbox_1[0] <= bbox_2[2])) & \
              ((bbox_1[3] >= bbox_2[1]) & (bbox_1[1] <= bbox_2[3]))
    return overlap
```

### Create features and boxes


```python
# Set plane dimensions
Nx, Ny = (100, 100)

# Set bbox size
wx = 10
wy = 10

# Create a random distribution of features
N = 20
X = np.random.RandomState(15).uniform(0, Nx, size=N)
Y = np.random.RandomState(38).uniform(0, Ny, size=N)

# Create DataFrame of features and their bboxes
features = pd.DataFrame({'x': X, 'y': Y})
features['x_min'] = X - wx/2
features['y_min'] = Y - wy/2
features['x_max'] = X + wx/2
features['y_max'] = Y + wy/2

# Sample features
features.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x</th>
      <th>y</th>
      <th>x_min</th>
      <th>y_min</th>
      <th>x_max</th>
      <th>y_max</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>84.881770</td>
      <td>38.477312</td>
      <td>79.881770</td>
      <td>33.477312</td>
      <td>89.881770</td>
      <td>43.477312</td>
    </tr>
    <tr>
      <th>1</th>
      <td>17.889592</td>
      <td>85.970785</td>
      <td>12.889592</td>
      <td>80.970785</td>
      <td>22.889592</td>
      <td>90.970785</td>
    </tr>
    <tr>
      <th>2</th>
      <td>5.436321</td>
      <td>94.419964</td>
      <td>0.436321</td>
      <td>89.419964</td>
      <td>10.436321</td>
      <td>99.419964</td>
    </tr>
    <tr>
      <th>3</th>
      <td>36.153845</td>
      <td>70.282489</td>
      <td>31.153845</td>
      <td>65.282489</td>
      <td>41.153845</td>
      <td>75.282489</td>
    </tr>
    <tr>
      <th>4</th>
      <td>27.540093</td>
      <td>63.363410</td>
      <td>22.540093</td>
      <td>58.363410</td>
      <td>32.540093</td>
      <td>68.363410</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Make figure
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(features['x'], features['y'], 'o', color='C1', ms=3, lw=0)
for i, row in features.iterrows():
    p = Rectangle((features.loc[i, 'x_min'], features.loc[i, 'y_min']),
                  wx, wy, facecolor='none', lw=1, edgecolor='C0')
    ax.add_patch(p)
    ax.text(features.loc[i, 'x_min'], features.loc[i, 'y_max'],
            f'{i}', ha='left', va='top')
ax.grid(ls=':')
```

### Cell listing method


```python
# Keep track of overlapping features
overlapping = []
# Define cell parameters
cw = 2*wx  # cell width

# Loop through a grid in x, y to create cells
for x in np.arange(0, Nx, cw):
    for y in np.arange(0, Ny, cw):
        # Create cell
        cell = [x-cw, y-cw, x+2*cw, y+2*cw]
        # Get features in cell
        in_cell = features[((cell[0] < X) & (X < cell[2]) &\
                            (cell[1] < Y) & (Y < cell[3]))]
        # Combinations
        pairs = list(combinations(in_cell.reset_index().values, 2))

        # Loop through pairs of bboxes
        for (bbox_i, bbox_j) in pairs:
            if bboxes_overlap(bbox_i[-4:], bbox_j[-4:]):
                overlapping.append(bbox_i[0])
                overlapping.append(bbox_j[0])

overlapping = np.unique(overlapping)
overlapping
```




    array([ 1.,  3.,  4.,  7., 11., 12., 14., 15., 16., 19.])



### Plot overlapping features


```python
def plot_overlapping_features(features, overlapping, plot_feature_index=False):
    """Plot features showing where bboxes overlap"""
    # Make bboxes
    features['x_min'] = features['x'] - wx/2
    features['y_min'] = features['y'] - wy/2
    features['x_max'] = features['x'] + wx/2
    features['y_max'] = features['y'] + wy/2

    # Make figure
    fig, ax = plt.subplots(figsize=(5, 5))
    ax.plot(features['x'], features['y'], 'o', color='C1', ms=3, lw=0)
    for i, row in features.iterrows():
        if i not in overlapping:
            p = Rectangle((features.loc[i, 'x_min'], features.loc[i, 'y_min']),
                          wx, wy, facecolor='none', lw=1, edgecolor='C0')
        else:
            p = Rectangle((features.loc[i, 'x_min'], features.loc[i, 'y_min']),
                          wx, wy, facecolor='none', lw=1, edgecolor='C3')
        ax.add_patch(p)
        if plot_feature_index:
            ax.text(features.loc[i, 'x_min'], features.loc[i, 'y_max'],
                    f'{i}', ha='left', va='top')
    ax.grid(ls=':')
```


```python
plot_overlapping_features(features, overlapping, plot_feature_index=True)
```

## Timing comparison
---


```python
from tqdm.notebook import tqdm
from psf_extractor import remove_overlapping_features, remove_overlapping_features_bruteforce
```

### Bigger feature set


```python
# Set plane dimensions
Nx, Ny = (2000, 2000)

# Create a random distribution of features
N = 300
X = np.random.RandomState(15).uniform(0, Nx, size=N)
Y = np.random.RandomState(38).uniform(0, Ny, size=N)
features = pd.DataFrame({'x': X, 'y': Y})

# Create a bounding box for each bead
wx = 80
wy = 80
```

### Old method (brute-force)


```python
_, overlapping = remove_overlapping_features_bruteforce(features, wx=wx, wy=wy, return_indices=True)

out = f"Overlapping features... {len(overlapping)} / {N}"
print(out)
```


      0%|          | 0/44850 [00:00<?, ?it/s]


    Overlapping features... 250 / 300



```python
plot_overlapping_features(features, overlapping)
```

### New method (cell listing)


```python
_, overlapping = remove_overlapping_features(features, wx=wx, wy=wy, return_indices=True)

out = f"Overlapping features... {len(overlapping)} / {N}"
print(out)
```


      0%|          | 0/14 [00:00<?, ?it/s]


    Overlapping features... 250 / 300



```python
plot_overlapping_features(features, overlapping)
```

d254b10e0a79ae69c79a49782c93ec19a3ac0744
ipynb
Jupyter Notebook
notebooks/cell_listing.ipynb
hoogenboom-group/PSF-Extractor
68fc52cf75e75481744b13ff97332d3330a400c5
CC-BY-4.0
2021-12-01T09:47:47.000Z
2021-12-01T09:47:47.000Z
notebooks/cell_listing.ipynb
hoogenboom-group/PSF-Extractor
68fc52cf75e75481744b13ff97332d3330a400c5
CC-BY-4.0
2021-08-15T19:28:20.000Z
2021-12-17T13:55:11.000Z
notebooks/cell_listing.ipynb
hoogenboom-group/PSF-Extractor
68fc52cf75e75481744b13ff97332d3330a400c5
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES