Merge pull request #1811 from AdeelH/strtree

Bump pygeos, shapely, and geopandas versions

rastervision_core/rastervision/core/data/utils/geojson.py

@@ -192,7 +192,7 @@ def split_multi_geometries(geojson: dict) -> dict:
     def split_geom(geom: 'BaseGeometry') -> List['BaseGeometry']:
         # Split GeometryCollection into list of geoms.
         if geom.geom_type == 'GeometryCollection':
-            geoms = list(geom)
+            geoms = list(geom.geoms)
         else:
             geoms = [geom]
         # Split any MultiX to list of X.

rastervision_core/rastervision/core/evaluation/object_detection_evaluation.py

@@ -1,12 +1,12 @@
 from typing import TYPE_CHECKING, Dict, Tuple
+
 import numpy as np
-from shapely.strtree import STRtree
+import geopandas as gpd
 
 from rastervision.core.evaluation import (ClassificationEvaluation,
                                           ClassEvaluationItem)
 
 if TYPE_CHECKING:
-    from shapely.geometry import Polygon
     from rastervision.core.data import ObjectDetectionLabels
     from rastervision.core.data.class_config import ClassConfig
 
@@ -16,41 +16,63 @@ def compute_metrics(
         pred_labels: 'ObjectDetectionLabels',
         num_classes: int,
         iou_thresh: float = 0.5) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Compute per-class true positves, false positives, and false negatives.
+
+    Does the following:
+
+    1.  Spatially join ground truth (GT) boxes with predicted boxes.
+    2.  Compute intersection-overo-union (IoU) for each matched box-pair.
+    3.  Filter matches by ``iou_thresh``.
+    4.  For each GT box >1 matches, keep only the max-IoU one.
+    5.  For each pred box >1 matches, keep only the max-IoU one.
+    6.  For each class, c, compute:
+
+        a.  True positives (TP) := #matches where GT class ID == c and
+            pred class ID == c
+        b.  False positives := #preds where (class ID == c) minus TP
+        c.  False negatives := #GT where (class ID == c) minus TP
+
+    """
     gt_geoms = [b.to_shapely() for b in gt_labels.get_boxes()]
     gt_classes = gt_labels.get_class_ids()
     pred_geoms = [b.to_shapely() for b in pred_labels.get_boxes()]
     pred_classes = pred_labels.get_class_ids()
 
-    for pred_geom, class_id in zip(pred_geoms, pred_classes):
-        pred_geom.class_id = class_id
-    pred_tree = STRtree(pred_geoms)
-
-    def iou(a: 'Polygon', b: 'Polygon') -> float:
-        return a.intersection(b).area / a.union(b).area
-
-    def is_matched(geom) -> bool:
-        return hasattr(geom, 'iou_matched')
+    gt_df = gpd.GeoDataFrame(
+        dict(class_id=gt_classes, id=range(len(gt_geoms))), geometry=gt_geoms)
+    pred_df = gpd.GeoDataFrame(
+        dict(class_id=pred_classes, id=range(len(pred_geoms))),
+        geometry=pred_geoms)
+
+    gt_df.loc[:, '_geometry'] = gt_df.geometry
+    pred_df.loc[:, '_geometry'] = pred_df.geometry
+
+    match_df: gpd.GeoDataFrame = gt_df.sjoin(
+        pred_df,
+        how='inner',
+        predicate='intersects',
+        lsuffix='gt',
+        rsuffix='pred')
+
+    intersection = match_df['_geometry_gt'].intersection(
+        match_df['_geometry_pred'])
+    union = match_df['_geometry_gt'].union(match_df['_geometry_pred'])
+    match_df.loc[:, 'iou'] = (intersection.area / union.area)
+    match_df = match_df.loc[match_df['iou'] > iou_thresh]
+    match_df = match_df.sort_values('iou').drop_duplicates(
+        ['id_gt'], keep='last')
+    match_df = match_df.sort_values('iou').drop_duplicates(
+        ['id_pred'], keep='last')
 
     tp = np.zeros((num_classes, ))
     fp = np.zeros((num_classes, ))
     fn = np.zeros((num_classes, ))
 
-    for gt_geom, gt_class in zip(gt_geoms, gt_classes):
-        matches = [
-            g for g in pred_tree.query(gt_geom)
-            if (not is_matched(g)) and (g.class_id == gt_class)
-        ]
-        ious = np.array([iou(m, gt_geom) for m in matches])
-        if (ious > iou_thresh).any():
-            max_ind = np.argmax(ious)
-            matches[max_ind].iou_matched = True
-            tp[gt_class] += 1
-        else:
-            fn[gt_class] += 1
-
     for class_id in range(num_classes):
-        pred_not_matched = np.array([not is_matched(g) for g in pred_geoms])
-        fp[class_id] = np.sum(pred_not_matched[pred_classes == class_id])
+        tp[class_id] = sum((match_df['class_id_gt'] == class_id)
+                           & (match_df['class_id_pred'] == class_id))
+        fp[class_id] = sum(pred_df['class_id'] == class_id) - tp[class_id]
+        fn[class_id] = sum(gt_df['class_id'] == class_id) - tp[class_id]
 
     return tp, fp, fn
 

rastervision_core/requirements.txt

@@ -2,9 +2,9 @@ rastervision_pipeline==0.20.3-dev
 
 # These 3 should be updated together to ensure compatibility. Incompatibility
 # results in a warning about pygeos versions.
-pygeos==0.13
-shapely==1.8.4
-geopandas==0.12.0
+pygeos==0.14
+shapely==2.0.1
+geopandas==0.13.2
 
 numpy==1.25.0
 pillow==9.3.0

rastervision_pytorch_learner/rastervision/pytorch_learner/dataset/utils/aoi_sampler.py

@@ -1,7 +1,7 @@
 from typing import Sequence, Tuple
 
 import numpy as np
-from shapely.geometry import Polygon, LinearRing, MultiPolygon
+from shapely.geometry import Polygon, MultiPolygon
 from shapely.ops import unary_union
 from triangle import triangulate
 
@@ -21,7 +21,9 @@ def __init__(self, polygons: Sequence[Polygon]) -> None:
         merged_polygons = unary_union(polygons)
         if isinstance(merged_polygons, Polygon):
             merged_polygons = [merged_polygons]
-        self.polygons = MultiPolygon(merged_polygons)
+        elif isinstance(merged_polygons, MultiPolygon):
+            merged_polygons = list(merged_polygons.geoms)
+        self.polygons = merged_polygons
         self.triangulate(self.polygons)
 
     def triangulate(self, polygons) -> dict:
@@ -65,7 +67,7 @@ def triangulate_polygon(self, polygon: Polygon) -> dict:
         """Extract vertices and edges from the polygon (and its holes, if any)
         and pass them to the Triangle library for triangulation.
         """
-        vertices, edges = self.polygon_to_graph(polygon.exterior)
+        vertices, edges = self.polygon_to_graph(polygon)
 
         holes = polygon.interiors
         if not holes:
@@ -107,17 +109,17 @@ def triangulate_polygon(self, polygon: Polygon) -> dict:
         return out
 
     def polygon_to_graph(self,
-                         polygon: LinearRing) -> Tuple[np.ndarray, np.ndarray]:
-        """Given the exterior of a polygon, return its graph representation.
+                         polygon: Polygon) -> Tuple[np.ndarray, np.ndarray]:
+        """Given a polygon, return its graph representation.
 
         Args:
-            polygon (LinearRing): The exterior of a polygon.
+            polygon (Polygon): A polygon.
 
         Returns:
             Tuple[np.ndarray, np.ndarray]: An (N, 2) array of vertices and
             an (N, 2) array of indices to vertices representing edges.
         """
-        vertices = np.array(polygon)
+        vertices = np.array(polygon.exterior.coords)
         # Discard the last vertex - it is a duplicate of the first vertex and
         # duplicates cause problems for the Triangle library.
         vertices = vertices[:-1]

tests/pytorch_learner/dataset/test_aoi_sampler.py

@@ -10,10 +10,9 @@
 
 class TestAoiSampler(unittest.TestCase):
     def test_sampler(self, nsamples: int = 200):
-        """This test attempts to check if the points are distributed uniformly
-        within the AOI.
+        """Attempt to check if points are distributed uniformly within an AOI.
 
-        To do this, it performs the following steps:
+        To do this, we perform the following steps:
         - Create an AOI in the form of a plus-sign shape centered in a 6x6
         grid.
         - Create an AoiSampler for this AOI.
@@ -33,8 +32,7 @@ def test_sampler(self, nsamples: int = 200):
         Args:
             nsamples (int, optional): Number of points to sample. It is
                 important for the sample size to not be too large or the test
-                will become over-powered.
-                Defaults to 200.
+                will become over-powered. Defaults to 200.
         """
         np.random.seed(0)
 
@@ -57,7 +55,7 @@ def test_sampler(self, nsamples: int = 200):
         points = MultiPoint(aoi_sampler.sample(n=nsamples))
         # number of points in each block
         counts = np.array(
-            [len(block.intersection(points)) for block in blocks])
+            [len(block.intersection(points).geoms) for block in blocks.geoms])
 
         p_value = chisquare(counts).pvalue