Remove shapely dependency (#204)

* Remove shapely dependency

* fixes

* remove print

* Add tests

Author: phil-scale

nucleus/metrics/filters.py

@@ -2,15 +2,18 @@
 
 from nucleus.prediction import PredictionList
 
-from .polygon_utils import BoxOrPolygonAnnotation, polygon_annotation_to_shape
+from .polygon_utils import (
+    BoxOrPolygonAnnotation,
+    polygon_annotation_to_geometry,
+)
 
 
 def polygon_area_filter(
     polygons: List[BoxOrPolygonAnnotation], min_area: float, max_area: float
 ) -> List[BoxOrPolygonAnnotation]:
     filter_fn = (
         lambda polygon: min_area
-        <= polygon_annotation_to_shape(polygon)
+        <= polygon_annotation_to_geometry(polygon).signed_area
         <= max_area
     )
     return list(filter(filter_fn, polygons))

nucleus/metrics/geometry.py

@@ -0,0 +1,184 @@
+from typing import List, Tuple, Union
+
+import numpy as np
+
+TOLERANCE = 1e-8
+
+
+class GeometryPoint:
+    def __init__(self, xy: Union[Tuple[float, float], np.ndarray]):
+        self.xy = np.array(xy)
+        self.x = xy[0]
+        self.y = xy[1]
+
+    def __repr__(self) -> str:
+        return f"GeometryPoint(xy=({self.xy[0]}, {self.xy[1]})"
+
+    def __add__(self, p: "GeometryPoint") -> "GeometryPoint":
+        return GeometryPoint(self.xy + p.xy)
+
+    def __sub__(self, p: "GeometryPoint") -> "GeometryPoint":
+        return GeometryPoint(self.xy - p.xy)
+
+    def __rmul__(self, scalar: float) -> "GeometryPoint":
+        return GeometryPoint(self.xy * scalar)
+
+    def __mul__(self, scalar: float) -> "GeometryPoint":
+        return GeometryPoint(self.xy * scalar)
+
+    # Operator @
+    def __matmul__(self, p: "GeometryPoint") -> float:
+        return self.xy @ p.xy
+
+    def length(self) -> float:
+        return np.linalg.norm(self.xy)
+
+    def cmp(self, p: "GeometryPoint") -> float:
+        return np.abs(self.xy - p.xy).sum()
+
+
+class GeometryPolygon:
+    def __init__(self, points: List[GeometryPoint]):
+        self.points = points
+
+        points_x = np.array([point.x for point in points])
+        points_y = np.array([point.y for point in points])
+        self.signed_area = (
+            points_x @ np.roll(points_y, 1) - points_x @ np.roll(points_y, -1)
+        ) / 2
+        self.area = np.abs(self.signed_area)
+
+    def __len__(self):
+        return len(self.points)
+
+    def __getitem__(self, idx):
+        return self.points[idx]
+
+    def __repr__(self) -> str:
+        points = ", ".join([str(point) for point in self.points])
+        return f"GeometryPolygon({points})"
+
+
+Segment = Tuple[GeometryPoint, GeometryPoint]
+
+
+# alpha * a1 + (1 - alpha) * a2 = beta * b1 + (1 - beta) * b2
+def segment_intersection(
+    segment1: Segment, segment2: Segment
+) -> Tuple[float, float, GeometryPoint]:
+    a1, a2 = segment1
+    b1, b2 = segment2
+    x2_x2 = b2.x - a2.x
+    y2_y2 = b2.y - a2.y
+    x1x2 = a1.x - a2.x
+    y1y2 = a1.y - a2.y
+    y1_y2_ = b1.y - b2.y
+    x1_x2_ = b1.x - b2.x
+
+    if np.abs(y1_y2_ * x1x2 - x1_x2_ * y1y2) < TOLERANCE:
+        beta = 1.0
+    else:
+        beta = (x2_x2 * y1y2 - y2_y2 * x1x2) / (y1_y2_ * x1x2 - x1_x2_ * y1y2)
+
+    if x1x2 == 0:
+        alpha = (y2_y2 + y1_y2_ * beta) / (y1y2 + TOLERANCE)
+    else:
+        alpha = (x2_x2 + x1_x2_ * beta) / (x1x2 + TOLERANCE)
+
+    return alpha, beta, alpha * a1 + (1 - alpha) * a2
+
+
+def convex_polygon_intersection_area(
+    polygon_a: GeometryPolygon, polygon_b: GeometryPolygon
+) -> float:
+    # pylint: disable=R0912
+    sa = polygon_a.signed_area
+    sb = polygon_b.signed_area
+    if sa * sb < 0:
+        sign = -1
+    else:
+        sign = 1
+    na = len(polygon_a)
+    nb = len(polygon_b)
+    ps = []  # point set
+    for i in range(na):
+        a1 = polygon_a[i - 1]
+        a2 = polygon_a[i]
+        flag = False
+        sum_s = 0
+        for j in range(nb):
+            b1 = polygon_b[j - 1]
+            b2 = polygon_b[j]
+            sum_s += np.abs(GeometryPolygon([a1, b1, b2]).signed_area)
+
+        if np.abs(np.abs(sum_s) - np.abs(sb)) < TOLERANCE:
+            flag = True
+
+        if flag:
+            ps.append(a1)
+        for j in range(nb):
+            b1 = polygon_b[j - 1]
+            b2 = polygon_b[j]
+            a, b, p = segment_intersection((a1, a2), (b1, b2))
+            if 0 < a < 1 and 0 < b < 1:
+                ps.append(p)
+
+    for i in range(nb):
+        a1 = polygon_b[i - 1]
+        a2 = polygon_b[i]
+        flag = False
+        sum_s = 0
+        for j in range(na):
+            b1 = polygon_a[j - 1]
+            b2 = polygon_a[j]
+            sum_s += np.abs(GeometryPolygon([a1, b1, b2]).signed_area)
+        if np.abs(np.abs(sum_s) - np.abs(sa)) < TOLERANCE:
+            flag = True
+        if flag:
+            ps.append(a1)
+
+    def unique(ar):
+        res = []
+        for i, _ in enumerate(ar):
+            if _.cmp(ar[i - 1]) > TOLERANCE:
+                res.append(_)
+
+        return res
+
+    ps = sorted(ps, key=lambda x: (x.x + TOLERANCE * x.y))
+    ps = unique(ps)
+
+    if len(ps) == 0:
+        return 0
+
+    tmp = ps[0]
+
+    res = []
+    res.append(tmp)
+    ps = sorted(
+        ps[1:],
+        key=lambda x: -(
+            (x - tmp) @ GeometryPoint((0, 1)) / (x - tmp).length()
+        ),
+    )
+    res.extend(ps)
+
+    return GeometryPolygon(res).signed_area * sign
+
+
+def polygon_intersection_area(
+    polygon_a: GeometryPolygon, polygon_b: GeometryPolygon
+) -> float:
+    na = len(polygon_a)
+    nb = len(polygon_b)
+    res = 0.0
+    for i in range(1, na - 1):
+        sa = [polygon_a[0], polygon_a[i], polygon_a[i + 1]]
+        for j in range(1, nb - 1):
+            sb = [polygon_b[0], polygon_b[j], polygon_b[j + 1]]
+            tmp = convex_polygon_intersection_area(
+                GeometryPolygon(sa), GeometryPolygon(sb)
+            )
+            res += tmp
+
+    return np.abs(res)

nucleus/metrics/polygon_utils.py

@@ -4,13 +4,13 @@
 
 import numpy as np
 from scipy.optimize import linear_sum_assignment
-from shapely.geometry import Polygon
 
 from nucleus.annotation import BoxAnnotation, PolygonAnnotation
 from nucleus.prediction import BoxPrediction, PolygonPrediction
 
 from .base import MetricResult
 from .errors import PolygonAnnotationTypeError
+from .geometry import GeometryPoint, GeometryPolygon, polygon_intersection_area
 
 BoxOrPolygonPrediction = TypeVar(
     "BoxOrPolygonPrediction", BoxPrediction, PolygonPrediction
@@ -20,31 +20,40 @@
 )
 
 
-def polygon_annotation_to_shape(
+def polygon_annotation_to_geometry(
     annotation: BoxOrPolygonAnnotation,
-) -> Polygon:
+) -> GeometryPolygon:
     if isinstance(annotation, BoxAnnotation):
         xmin = annotation.x - annotation.width / 2
         xmax = annotation.x + annotation.width / 2
         ymin = annotation.y - annotation.height / 2
         ymax = annotation.y + annotation.height / 2
-        return Polygon(
-            [(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)]
-        )
+        points = [
+            GeometryPoint((xmin, ymin)),
+            GeometryPoint((xmax, ymin)),
+            GeometryPoint((xmax, ymax)),
+            GeometryPoint((xmin, ymax)),
+        ]
+        return GeometryPolygon(points)
     elif isinstance(annotation, PolygonAnnotation):
-        return Polygon([(point.x, point.y) for point in annotation.vertices])
+        return GeometryPolygon(
+            [
+                GeometryPoint((point.x, point.y))
+                for point in annotation.vertices
+            ]
+        )
     else:
         raise PolygonAnnotationTypeError()
 
 
-def _iou(annotation: Polygon, prediction: Polygon) -> float:
-    intersection = annotation.intersection(prediction).area
+def _iou(annotation: GeometryPolygon, prediction: GeometryPolygon) -> float:
+    intersection = polygon_intersection_area(annotation, prediction)
     union = annotation.area + prediction.area - intersection
     return intersection / max(union, sys.float_info.epsilon)
 
 
 def _iou_matrix(
-    annotations: List[Polygon], predictions: List[Polygon]
+    annotations: List[GeometryPolygon], predictions: List[GeometryPolygon]
 ) -> np.ndarray:
     iou_matrix = np.empty((len(predictions), len(annotations)))
     for i, prediction in enumerate(predictions):
@@ -65,9 +74,13 @@ def _iou_assignments_for_same_reference_id(
         len(reference_ids) <= 1
     ), "Expected annotations and predictions to have same reference ID."
 
-    # Convert annotation and predictions to shapely.geometry.Polygon objects
-    polygon_annotations = list(map(polygon_annotation_to_shape, annotations))
-    polygon_predictions = list(map(polygon_annotation_to_shape, predictions))
+    # Convert annotation and predictions to GeometryPolygon objects
+    polygon_annotations = list(
+        map(polygon_annotation_to_geometry, annotations)
+    )
+    polygon_predictions = list(
+        map(polygon_annotation_to_geometry, predictions)
+    )
 
     # Compute IoU matrix and set IoU values below the threshold to 0.
     iou_matrix = _iou_matrix(polygon_annotations, polygon_predictions)

pyproject.toml

@@ -43,7 +43,6 @@ pydantic = "^1.8.2"
 isort = "^5.10.1"
 numpy = "^1.19.5"
 scipy = "^1.5.4"
-Shapely = "^1.8.0"
 
 [tool.poetry.dev-dependencies]
 poetry = "^1.1.5"

tests/metrics/test_geometry.py

@@ -0,0 +1,59 @@
+import numpy as np
+import pytest
+
+from nucleus.metrics.geometry import (
+    GeometryPoint,
+    GeometryPolygon,
+    convex_polygon_intersection_area,
+    polygon_intersection_area,
+    segment_intersection,
+)
+
+RECTANGLE1 = GeometryPolygon(
+    [
+        GeometryPoint((0, 0)),
+        GeometryPoint((100, 0)),
+        GeometryPoint((100, 100)),
+        GeometryPoint((0, 100)),
+    ]
+)
+
+RECTANGLE2 = GeometryPolygon(
+    [
+        GeometryPoint((50, 50)),
+        GeometryPoint((50, 150)),
+        GeometryPoint((150, 150)),
+        GeometryPoint((150, 50)),
+    ]
+)
+
+SEGMENT1 = (GeometryPoint((0, 0)), GeometryPoint((0, 50)))
+
+SEGMENT2 = (GeometryPoint((-25, 25)), GeometryPoint((25, 25)))
+
+
+def test_segment_intersection():
+    alpha, beta, intersection = segment_intersection(SEGMENT1, SEGMENT2)
+    assert alpha == pytest.approx(0.5)
+    assert beta == pytest.approx(0.5)
+    assert intersection.x == pytest.approx(0)
+    assert intersection.y == pytest.approx(25)
+
+
+def test_convex_polygon_intersection_area():
+    intersection = np.abs(
+        convex_polygon_intersection_area(RECTANGLE1, RECTANGLE2)
+    )
+    assert intersection == pytest.approx(2500)
+
+
+def test_polygon_intersection_area():
+    intersection = polygon_intersection_area(RECTANGLE1, RECTANGLE2)
+    assert intersection == pytest.approx(2500)
+
+
+def test_polygon_area():
+    assert RECTANGLE1.signed_area == -10000
+    assert RECTANGLE1.area == 10000
+    assert RECTANGLE2.signed_area == 10000
+    assert RECTANGLE2.area == 10000