diff --git a/Cargo.toml b/Cargo.toml
index 1512793..e0c427f 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "neighborhood_analysis"
-version = "0.1.4"
+version = "0.1.5"
 authors = ["Mr-Milk <zym.zym1220@gmail.com>"]
 edition = "2018"
 license = "MIT"
@@ -28,6 +28,7 @@ rand = "0.7.3"
 kdbush = "0.2.0"
 rayon = "1.4.0"
 rstar = "0.8.2"
+spade = "1.8.2"
 
 [profile.dev]
 opt-level = 3
diff --git a/README.md b/README.md
index 6d7fe5c..9ec313e 100644
--- a/README.md
+++ b/README.md
@@ -78,6 +78,17 @@ z_score = comb_bootstrap(X, Y, neighbors, ignore_self=True)
 
 ```python
 
+def get_bbox(points_collections):
+    """A utility function to return minimum bounding box list of polygons
+    
+        Args:
+            points_collections: List[List[(float, float)]]; List of 2d points collections
+        
+        Return:
+            A dictionary of the index of every points, with the index of its neighbors
+
+    """
+
 def get_point_neighbors(points, r):
     """A utility function to search for neighbors
     
diff --git a/neighborhood_analysis/__init__.py b/neighborhood_analysis/__init__.py
index 56deacc..80dde97 100644
--- a/neighborhood_analysis/__init__.py
+++ b/neighborhood_analysis/__init__.py
@@ -1 +1 @@
-from .neighborhood_analysis import CellCombs, get_point_neighbors, get_bbox_neighbors, comb_bootstrap
+from .neighborhood_analysis import CellCombs, get_bbox, get_point_neighbors, get_bbox_neighbors, comb_bootstrap
diff --git a/src/lib.rs b/src/lib.rs
index c3e6afe..3e94019 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -6,6 +6,7 @@ use rand::seq::SliceRandom;
 use itertools::Itertools;
 use kdbush::KDBush;
 use rstar::{RTree, RTreeObject, AABB};
+use spade::{BoundingRect};
 use std::collections::HashMap;
 use counter::Counter;
 use rayon::prelude::*;
@@ -18,12 +19,38 @@ use pyo3::wrap_pyfunction;
 #[pymodule]
 fn neighborhood_analysis(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_class::<CellCombs>()?;
+    m.add_wrapped(wrap_pyfunction!(get_bbox))?;
     m.add_wrapped(wrap_pyfunction!(get_point_neighbors))?;
     m.add_wrapped(wrap_pyfunction!(get_bbox_neighbors))?;
     m.add_wrapped(wrap_pyfunction!(comb_bootstrap))?;
     Ok(())
 }
 
+
+/// A utility function to return minimum bounding box list of polygons
+///
+/// Args:
+///     points_collections: List[List[(float, float)]]; List of 2d points collections
+///
+/// Return:
+///     A dictionary of the index of every points, with the index of its neighbors
+///
+#[pyfunction]
+fn get_bbox(points_collections: Vec<Vec<(f64, f64)>>)
+-> Vec<(f64, f64, f64, f64)> {
+
+    let bbox: Vec<(f64, f64, f64, f64)> = points_collections.par_iter().map(|p| {
+        let points: Vec<[f64;2]> = p.iter().map(|ps| [ps.0, ps.1]).collect();
+        let rect = BoundingRect::from_points(points);
+        let lower: [f64;2] = rect.lower();
+        let upper: [f64;2] = rect.upper();
+        (lower[0], lower[1], upper[0], upper[1])
+    }).collect();
+
+    bbox
+}
+
+
 /// A utility function to search for point neighbors using kd-tree
 ///
 /// Args:
diff --git a/tests.py b/tests.py
index 29ed703..78d219a 100644
--- a/tests.py
+++ b/tests.py
@@ -1,6 +1,6 @@
 import numpy as np
 import neighborhood_analysis as na
-from neighborhood_analysis import CellCombs, get_point_neighbors, get_bbox_neighbors,comb_bootstrap
+from neighborhood_analysis import CellCombs, get_bbox, get_point_neighbors, get_bbox_neighbors,comb_bootstrap
 
 from time import time
 
@@ -9,13 +9,17 @@
 corr_types = np.random.choice(types, 10000)
 points = [(x, y) for (x, y) in points]
 
-bbox = []
-for _ in range(100):
-    ix1, ix2 = np.random.choice(range(len(points)), 2)
-    if points[ix1][0] < points[ix2][0]:
-        bbox.append((*points[ix1], *points[ix2]))
-    else:
-        bbox.append((*points[ix2], *points[ix1]))
+polygons = []
+for _ in range(10000):
+    ixs = np.random.choice(range(len(points)), 5)
+    polygon = []
+    for x in ixs:
+        polygon.append(points[x])
+
+start = time()
+bbox = get_bbox(polygons)
+end = time()
+print(f"Get bbox used {(end-start):.2f}s")
 
 start = time()
 neighbors = get_bbox_neighbors(bbox, 2)