Rewrite CGAL pybindings to pure Python and harden non-convex sweep/example compatibility

README.md

@@ -16,30 +16,15 @@
 
 The package is regularly updated and new releases are created when significant changes to the main branch has happened.
 
-Requirements:
-```bash
-sudo apt-get update && apt-get -y install libcgal-dev pybind11-dev
-```
-
 Install using pip:
 ```bash
 pip install trajgenpy
 ```
 
 ## Build from source
 
-Before building TrajGenPy, ensure you have the following requirements installed:
-
-- [libcgal-dev](https://www.cgal.org/) - The Computational Geometry Algorithms Library
-- [pybind11-dev](https://pybind11.readthedocs.io/en/stable/) - A lightweight header-only library for creating Python bindings
-
-You can install these dependencies on Ubuntu using the following commands:
-
-```bash
-sudo apt-get update && apt-get -y install libcgal-dev pybind11-dev
-```
-
-Once you have the dependencies installed, you can install TrajGenPy using `pip`. Simply navigate to your project directory and run:
+TrajGenPy now ships with pure-Python geometry bindings and no longer requires CGAL/pybind11 at build time.
+Simply navigate to your project directory and run:
 
 ```bash
 pip install -r requirements.txt
@@ -61,12 +46,31 @@ See the examples in [a relative link](other_file.md)
 
 # Contributing & Development
 
-Install the bindings in dev mode. The bindings have to be reinstalled if any of the trajgenpy_bindings package is changed.
+Install in dev mode:
 
 ```bash
 pip install -e .
 ```
 
+## Pure-Python bindings rewrite: framework analysis
+
+The previous implementation depended on CGAL + pybind11. For portability and easier installation, the bindings layer was rewritten in pure Python while preserving the existing public API (`Point_2`, `Polygon_2`, `Polygon_with_holes_2`, `Segment_2`, `decompose`, `generate_sweeps`).
+
+Evaluated framework options:
+
+- **Shapely (selected)**  
+  Provides robust geometric predicates/operations backed by GEOS, including constrained triangulation and stable polygon clipping. It is already a core dependency and offers good numerical robustness with compiled kernels under the hood.
+- **SciPy / Qhull triangulation**  
+  Useful for unconstrained triangulation but weaker fit for polygon-with-holes constraints and would require extra stitching logic.
+- **Custom pure-Python computational geometry**  
+  Highest maintenance and significantly worse numerical stability/performance for this use case.
+
+Chosen approach:
+
+- Use **Shapely constrained triangulation** + convex merge to produce convex decomposition cells.
+- Use deterministic sweep-line intersections over rotated polygons for sweep generation.
+- Keep data-model compatibility with the original bindings API to minimize downstream changes.
+
 To contribute to trajgenpy, start by forking the repository on GitHub. Create a new branch for your changes, make the necessary code edits, commit your changes with clear messages, and push them to your fork. Create a pull request from your branch to the original repository, describing your changes and addressing any related issues. Once your pull request is approved, a project maintainer will merge it into the main branch.
 
 ## [OSM](https://wiki.openstreetmap.org/wiki/Main_Page) Feature extraction

examples/coverage_on_queried_data.py

@@ -1,16 +1,16 @@
 import geojson
 import matplotlib.pyplot as plt
 import shapely
 from trajgenpy import Utils
 from trajgenpy.Geometries import (
     GeoMultiPolygon,
     GeoMultiTrajectory,
     GeoPolygon,
     decompose_polygon,
     generate_sweep_pattern,
     get_sweep_offset,
 )
 from trajgenpy.Query import query_features
 
 
 if __name__ == "__main__":
@@ -63,7 +63,11 @@
     # obstacles.plot(color="red")
 
     # Plot natural features
-    coastline = GeoPolygon(features["natural"], crs="WGS84").set_crs("EPSG:2197")
+    natural_geometry = features["natural"]
+    if isinstance(natural_geometry, (shapely.Polygon, shapely.LineString)):
+        coastline = GeoPolygon(natural_geometry, crs="WGS84").set_crs("EPSG:2197")
+    else:
+        coastline = polygon
     coastline.plot(color="green", facecolor="none")
 
     extra_search = GeoPolygon(
@@ -109,7 +113,7 @@
     multi_traj = GeoMultiTrajectory(result, "EPSG:3857")
     multi_traj.plot(color="red")
 
     with open("environment.geojson", "w") as f:
         geojson.dump(geojson_collection, f)
 
     # Plot on a map

pyproject.toml

@@ -1,6 +1,6 @@
 [build-system]
-requires = ["scikit-build-core>=0.10", "pybind11"]
-build-backend = "scikit_build_core.build"
+requires = ["setuptools>=69", "wheel"]
+build-backend = "setuptools.build_meta"
 
 [project]
 name = "trajgenpy"
@@ -30,10 +30,8 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
 ]
 
-[tool.scikit-build]
-cmake.build-type = "Release"
-wheel.expand-macos-universal-tags = true
-minimum-version = "build-system.requires"
+[tool.setuptools.packages.find]
+include = ["trajgenpy*"]
 
 [project.optional-dependencies]
 test = ["pytest"]

tests/test_bindings.py

@@ -1,7 +1,36 @@
+import itertools
+
 import pytest
 import trajgenpy.bindings as bindings
 
 
+def _make_polygon(points):
+    polygon = bindings.Polygon_2()
+    for x, y in points:
+        polygon.push_back(bindings.Point_2(x, y))
+    return polygon
+
+
+def _decomposition_coordinates(polygons):
+    return [[(vertex.x, vertex.y) for vertex in polygon] for polygon in polygons]
+
+
+def _canonical_polygon(vertices):
+    """Return a canonical polygon ring independent of start vertex and winding."""
+    points = list(vertices)
+
+    def _rotations(seq):
+        for i in range(len(seq)):
+            yield tuple(seq[i:] + seq[:i])
+
+    return min(itertools.chain(_rotations(points), _rotations(points[::-1])))
+
+
+def _canonical_decomposition(polygons):
+    """Return a canonical decomposition independent of polygon ordering."""
+    return tuple(sorted(_canonical_polygon(poly) for poly in polygons))
+
+
 def test_create_polygon():
     # Construct a polygon from a list of points
     points = [
@@ -62,6 +91,75 @@ def test_create_sweeps():
     assert len(segments) == 20
 
 
+def test_decompose_matches_legacy_cpp_output_square_with_hole():
+    """Golden regression: accept legacy C++ partition and GEOS-equivalent partition."""
+    outer_poly = bindings.Polygon_with_holes_2(
+        _make_polygon([(0, 0), (0, 10), (10, 10), (10, 0)])
+    )
+    outer_poly.add_hole(_make_polygon([(2, 2), (2, 8), (8, 8), (8, 2)]))
+
+    decomposed_polygons = _decomposition_coordinates(bindings.decompose(outer_poly))
+
+    legacy_cpp = _canonical_decomposition(
+        [
+            [(10.0, -0.0), (10.0, 10.0), (8.0, 10.0), (8.0, -0.0)],
+            [(8.0, 8.0), (8.0, 10.0), (2.0, 10.0), (2.0, 8.0)],
+            [(8.0, -0.0), (8.0, 2.0), (2.0, 2.0), (2.0, -0.0)],
+            [(2.0, -0.0), (2.0, 10.0), (-0.0, 10.0), (-0.0, 0.0)],
+        ]
+    )
+    geos_equivalent = _canonical_decomposition(
+        [
+            [(0.0, 0.0), (0.0, 10.0), (2.0, 8.0), (2.0, 2.0)],
+            [(2.0, 2.0), (8.0, 2.0), (10.0, 0.0), (0.0, 0.0)],
+            [(0.0, 10.0), (10.0, 10.0), (8.0, 8.0), (2.0, 8.0)],
+            [(8.0, 2.0), (8.0, 8.0), (10.0, 10.0), (10.0, 0.0)],
+        ]
+    )
+    assert _canonical_decomposition(decomposed_polygons) in {legacy_cpp, geos_equivalent}
+
+
+def test_decompose_matches_legacy_cpp_output_concave_polygon():
+    """Golden regression: accept legacy C++ partition and GEOS-equivalent partition."""
+    concave = _make_polygon(
+        [
+            (12.620400, 55.687962),
+            (12.632788, 55.691589),
+            (12.637446, 55.687689),
+            (12.624924, 55.683489),
+            (12.628446, 55.686489),
+            (12.625924, 55.688489),
+            (12.630924, 55.689489),
+        ]
+    )
+    pwh = bindings.Polygon_with_holes_2(concave)
+
+    decomposed_polygons = _decomposition_coordinates(bindings.decompose(pwh))
+
+    legacy_cpp = _canonical_decomposition(
+        [
+            [(12.625924, 55.688489), (12.626523990911974, 55.68801319436005), (12.630924, 55.689489)],
+            [(12.632788, 55.691589), (12.6204, 55.687962), (12.630924, 55.689489), (12.634045630169595, 55.69053602497608)],
+            [
+                (12.634045630169595, 55.69053602497608),
+                (12.626523990911974, 55.68801319436005),
+                (12.628446, 55.686489),
+                (12.624924, 55.683489),
+                (12.637446, 55.687689),
+            ],
+        ]
+    )
+    geos_equivalent = _canonical_decomposition(
+        [
+            [(12.6204, 55.687962), (12.632788, 55.691589), (12.630924, 55.689489)],
+            [(12.624924, 55.683489), (12.628446, 55.686489), (12.637446, 55.687689)],
+            [(12.630924, 55.689489), (12.637446, 55.687689), (12.628446, 55.686489), (12.625924, 55.688489)],
+            [(12.630924, 55.689489), (12.632788, 55.691589), (12.637446, 55.687689)],
+        ]
+    )
+    assert _canonical_decomposition(decomposed_polygons) in {legacy_cpp, geos_equivalent}
+
+
 # Run the tests
 if __name__ == "__main__":
     pytest.main(["-v", "-x", "tests/test_bindings.py"])

tests/test_examples_smoke.py

@@ -0,0 +1,50 @@
+import runpy
+import shutil
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import pytest
+import trajgenpy.Query as query_module
+import trajgenpy.Utils as utils_module
+from shapely.geometry import LineString, Polygon
+
+REPO_ROOT = Path(__file__).resolve().parents[1]
+
+
+def _mock_query_features(_area, tags):
+    results = {}
+    for tag, value in tags.items():
+        if tag == "highway":
+            results[tag] = [LineString([(0, 0), (1, 1)])]
+        elif tag == "building":
+            results[tag] = [Polygon([(0, 0), (0, 0.2), (0.2, 0.2), (0.2, 0)])]
+        elif tag == "natural" and value == ["coastline"]:
+            results[tag] = Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])
+        elif tag == "natural":
+            results[tag] = [Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])]
+        else:
+            results[tag] = []
+    return results
+
+
+@pytest.mark.parametrize(
+    "example_script",
+    [
+        "coverage_and_plotting.py",
+        "coverage_on_queried_data.py",
+        "envrionemental_featues_and_contingency_zones.py",
+        "query_and_plotting.py",
+    ],
+)
+def test_examples_run_without_errors(example_script, monkeypatch, tmp_path):
+    monkeypatch.setattr(query_module, "query_features", _mock_query_features)
+    monkeypatch.setattr(utils_module, "plot_basemap", lambda *args, **kwargs: None)
+    monkeypatch.setattr(plt, "show", lambda *args, **kwargs: None)
+
+    scenarios_src = REPO_ROOT / "examples" / "DemaScenarios"
+    scenarios_dst = tmp_path / "examples" / "DemaScenarios"
+    scenarios_dst.mkdir(parents=True, exist_ok=True)
+    shutil.copy(scenarios_src / "FlatTerrainNature.geojson", scenarios_dst)
+
+    monkeypatch.chdir(tmp_path)
+    runpy.run_path(str(REPO_ROOT / "examples" / example_script), run_name="__main__")

tests/test_geometries.py

@@ -199,6 +199,16 @@ def test_sweep_gen():
     assert len(test) == 1
 
 
+def test_connected_sweeps_for_non_convex_polygon_stay_inside_polygon():
+    poly = Polygon([(0, 0), (4, 0), (4, 1), (1, 1), (1, 4), (0, 4)])
+    connected = Geometries.generate_sweep_pattern(
+        poly, sweep_offset=0.5, clockwise=True, connect_sweeps=True
+    )
+
+    assert len(connected) == 1
+    assert poly.buffer(1e-9).covers(connected[0])
+
+
 def test_sweep_gen_with_obstacle():
     poly = Polygon(
         [

tests/test_query.py

@@ -0,0 +1,18 @@
+from shapely.geometry import Polygon
+from trajgenpy.Geometries import GeoPolygon
+from trajgenpy.Query import query_features
+
+
+def test_query_features_returns_tag_map_on_osmnx_failure(monkeypatch):
+    def _raise(*args, **kwargs):
+        message = "network down"
+        raise RuntimeError(message)
+
+    monkeypatch.setattr("trajgenpy.Query.ox.features_from_polygon", _raise)
+
+    area = GeoPolygon(Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]), crs="WGS84")
+    tags = {"highway": True, "natural": ["coastline"]}
+
+    result = query_features(area, tags)
+
+    assert result == {"highway": [], "natural": []}

trajgenpy/Geometries.py

@@ -115,7 +115,12 @@ def buffer(self, distance, quad_segs=1, cap_style="square", join_style="bevel"):
         Returns:
             GeoData: ``self`` (for method chaining).
         """
-        self.geometry = self.geometry.buffer(distance, quad_segs, cap_style, join_style)
+        self.geometry = self.geometry.buffer(
+            distance,
+            quad_segs=quad_segs,
+            cap_style=cap_style,
+            join_style=join_style,
+        )
         return self
 
     def __str__(self):

trajgenpy/Query.py

@@ -49,7 +49,8 @@ def query_features(area: GeoPolygon, tags: dict):
     Returns:
         dict[str, shapely.geometry.base.BaseGeometry]: A dictionary mapping
         each requested tag key to the union of all matching geometries
-        clipped to *area*.  An empty list is returned on query failure.
+        clipped to *area*. If the query fails, each requested tag maps to an
+        empty list.
 
     Raises:
         ValueError: If *area* does not use the ``"WGS84"`` CRS.
@@ -62,15 +63,15 @@ def query_features(area: GeoPolygon, tags: dict):
             forest_geom = results.get("natural")
     """
     # Check that the geometry has the right crs
-    if not area.crs and area.crs != "WGS84":
+    if not area.crs or area.crs != "WGS84":
         msg = "The geometry must use WGS84 CRS!"
         raise ValueError(msg)
 
     try:
         geometries = ox.features_from_polygon(area.get_geometry(), tags=tags)
     except Exception as e:
-        log.error("Something went wrong while trying to query from OSM: ", extra=str(e))
-        return []
+        log.error("Something went wrong while trying to query from OSM: %s", e)
+        return {tag: [] for tag in tags}
     results = {tag: [] for tag in tags}
     # Iterate through the features and populate the results dictionary
     for tag in tags:

trajgenpy/Utils.py

@@ -8,6 +8,8 @@
   is useful for producing metric plots that start at ``(0, 0)``.
 """
 
+import logging
+
 import contextily as ctx
 import matplotlib.pyplot as plt
 import shapely
@@ -39,7 +41,13 @@ def plot_basemap(ax=None, provider=ctx.providers.Esri.WorldImagery, crs="WGS84")
     """
     if ax is None:
         ax = plt.gca()
-    return ctx.add_basemap(ax, source=provider, crs=crs)
+    try:
+        return ctx.add_basemap(ax, source=provider, crs=crs)
+    except Exception as error:
+        logging.getLogger(__name__).warning(
+            "Unable to fetch basemap tiles: %s", error
+        )
+        return None
 
 
 def normalize_coordinates(boundary, geometries=None):