Handle unsupported geometry types

sleap_roots/convhull.py

@@ -1,6 +1,7 @@
 """Convex hull fitting and derived trait calculation."""
 
 import numpy as np
+import logging
 from scipy.spatial import ConvexHull
 from scipy.spatial.distance import pdist
 from typing import Tuple, Optional, Union
@@ -30,13 +31,13 @@
 
     # Check for infinite values
     if np.isinf(pts).any():
-        print("Cannot compute convex hull: input contains infinite values.")
+        logging.info("Cannot compute convex hull: input contains infinite values.")
         return None
 
     # Ensure there are at least 3 unique non-collinear points
     unique_pts = np.unique(pts, axis=0)
     if len(unique_pts) < 3:
-        print("Cannot compute convex hull: not enough unique points.")
+        logging.info("Cannot compute convex hull: not enough unique points.")
         return None
 
     try:
@@ -384,16 +385,25 @@
 
     # Find the intersection between the hull perimeter and the extended line
     intersection = extended_line.intersection(hull_perimeter)
+    logging.debug(f"Intersection: {intersection}")
 
     # Compute the intersection points and add to lists
-    if not intersection.is_empty:
+    if intersection and not intersection.is_empty:
+        logging.debug("Intersection points found between the convex hull and the line.")
         intersect_points = extract_points_from_geometry(intersection)
-        above_line.extend(intersect_points)
-        below_line.extend(intersect_points)
+        if not intersect_points: # Ensure it's not an empty list
+            return np.nan, np.nan
+    else:
+        logging.debug("No intersection points found between the convex hull and the line.")
+        return np.nan, np.nan
 
+    # Add intersection points to the lists
+    above_line.extend(intersect_points)
+    below_line.extend(intersect_points)
+
     # Calculate areas using get_chull_area
-    area_above_line = get_chull_area(np.array(above_line)) if above_line else 0.0
-    area_below_line = get_chull_area(np.array(below_line)) if below_line else 0.0
+    area_above_line = get_chull_area(np.array(above_line)) if above_line else np.nan
+    area_below_line = get_chull_area(np.array(below_line)) if below_line else np.nan
 
     return area_above_line, area_below_line
 
@@ -475,19 +485,22 @@
 
     # Check for a valid or existing convex hull
     if hull is None or len(unique_pts) < 3:
+        logging.debug("Not enough unique points to compute convex hull intersections.")
         # Return two vectors of NaNs if not valid hull
         return (np.array([[np.nan, np.nan]]), np.array([[np.nan, np.nan]]))
 
     # Ensure rn_pts does not contain NaN values
     rn_pts_valid = rn_pts[~np.isnan(rn_pts).any(axis=1)]
     # Need at least two points to define a line
     if len(rn_pts_valid) < 2:
+        logging.debug("Not enough valid rn points to compute convex hull intersections.")
         return (np.array([[np.nan, np.nan]]), np.array([[np.nan, np.nan]]))
 
     # Ensuring r0_pts does not contain NaN values
     r0_pts_valid = r0_pts[~np.isnan(r0_pts).any(axis=1)]
     # Expect two vectors in the end
     if len(r0_pts_valid) < 2:
+        logging.debug("Not enough valid r0 points to compute convex hull intersections.")
         return (np.array([[np.nan, np.nan]]), np.array([[np.nan, np.nan]]))
 
     # Get the vertices of the convex hull
@@ -509,6 +522,7 @@
     leftmost_vector = np.array([[np.nan, np.nan]])
     rightmost_vector = np.array([[np.nan, np.nan]])
     if not is_leftmost_on_hull and not is_rightmost_on_hull:
+        logging.debug("Leftmost and rightmost r0 points are not on the convex hull.")
         # If leftmost and rightmost r0 points are not on the convex hull return NaNs
         return leftmost_vector, rightmost_vector
 
@@ -518,6 +532,7 @@
         rightmost_rn = rn_pts[np.argmax(rn_pts[:, 0])]
         m, b = get_line_equation_from_points(leftmost_rn, rightmost_rn)
     except Exception:
+        logging.debug("Could not find line equation between leftmost and rightmost rn points.")
         # If line equation cannot be found, return NaNs
         return leftmost_vector, rightmost_vector
 
@@ -545,11 +560,17 @@
 
     # Find the intersection between the hull perimeter and the extended line
     intersection = extended_line.intersection(hull_perimeter)
+    logging.debug(f"Intersection: {intersection}")
 
     # Get the intersection points
-    if not intersection.is_empty:
+    if intersection and not intersection.is_empty:
+        logging.debug(f"Intersection points found between the convex hull and the line: {intersection}.")
         intersect_points = extract_points_from_geometry(intersection)
+        if not intersect_points:  # Ensure it's not an empty list
+            logging.debug("No intersection points found after extraction.")
+            return leftmost_vector, rightmost_vector
     else:
+        logging.debug("No intersection points found between the convex hull and the line.")
         # Return two vectors of NaNs if there is no intersection
         return leftmost_vector, rightmost_vector
 

sleap_roots/points.py

@@ -1,6 +1,7 @@
 """Get traits related to the points."""
 
 import numpy as np
+import logging
 from matplotlib import pyplot as plt
 from matplotlib.lines import Line2D
 from shapely.geometry import Point, MultiPoint, LineString, GeometryCollection
@@ -31,6 +32,8 @@ def extract_points_from_geometry(geometry) -> List[np.ndarray]:
     >>> extract_points_from_geometry(geom_col)
     [array([1, 2]), array([1, 2]), array([3, 4]), array([0, 0]), array([1, 1]), array([2, 2])]
     """
+    logging.debug(f"Geometry type: {type(geometry).__name__}")
+    logging.debug(f"Geometry: {geometry}")
     if isinstance(geometry, Point):
         return [np.array([geometry.x, geometry.y])]
     elif isinstance(geometry, MultiPoint):
@@ -43,7 +46,9 @@ def extract_points_from_geometry(geometry) -> List[np.ndarray]:
             points.extend(extract_points_from_geometry(geom))
         return points
     else:
-        raise TypeError(f"Unsupported geometry type: {type(geometry).__name__}")
+        logging.info(f"Unsupported geometry type: {type(geometry).__name__}")
+        return []
+
 
 
 def get_count(pts: np.ndarray):

tests/test_points.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from shapely.geometry import Point, MultiPoint, LineString, GeometryCollection
+from shapely.geometry import Point, MultiPoint, LineString, GeometryCollection, MultiLineString
 from sleap_roots import Series
 from sleap_roots.lengths import get_max_length_pts
 from sleap_roots.points import (
@@ -779,6 +779,13 @@ def test_extract_from_linestring():
     results = extract_points_from_geometry(linestring)
     assert all(np.array_equal(result, exp) for result, exp in zip(results, expected))
 
+def test_extract_from_multilinestring():
+    multilinestring = MultiLineString([[(0, 0), (1, 1)], [(2, 2), (3, 3)]])
+    # Multilinestring is not supported, so it should return an empty list
+    expected = []
+    results = extract_points_from_geometry(multilinestring)
+    assert all(np.array_equal(result, exp) for result, exp in zip(results, expected))
+
 
 def test_extract_from_geometrycollection():
     geom_collection = GeometryCollection([Point(1, 2), LineString([(0, 0), (1, 1)])])
@@ -803,10 +810,35 @@ def test_extract_from_unsupported_type():
     with pytest.raises(NameError):
         extract_points_from_geometry(
             Polygon([(0, 0), (1, 1), (1, 0)])
-        )  # Polygon is unsupported
+        )  # Polygon is unsupported type and not imported from shapely.geometry
 
 
 def test_extract_from_empty_geometrycollection():
     empty_geom_collection = GeometryCollection()
     expected = []
     assert extract_points_from_geometry(empty_geom_collection) == expected
+
+
+@pytest.mark.parametrize(
+    "geometry, expected",
+    [
+
+        (MultiLineString([[(0, 0), (1, 1)], [(2, 2), (3, 3)]]), []),
+        (GeometryCollection([Point(1, 2), LineString([(0, 0), (1, 1)]), MultiLineString([[(0, 0), (1, 1)], [(2, 2), (3, 3)]])]), [np.array([1, 2]), np.array([0, 0]), np.array([1, 1])]), # GeometryCollection with MultiLineString
+    ],
+)
+def test_extract_from_multilinestring(geometry, expected):
+    results = extract_points_from_geometry(geometry)
+    assert all(np.array_equal(result, exp) for result, exp in zip(results, expected))
+
+
+@pytest.mark.parametrize(
+    "unexpected_input, expected_output",
+    [
+        ("24", []),
+        (None, []),
+        (5, []),
+    ],
+)
+def test_extract_from_unsupported_geometry(unexpected_input, expected_output):
+    assert extract_points_from_geometry(unexpected_input) == expected_output

tests/test_trait_pipelines.py

@@ -182,23 +182,22 @@ def test_younger_monocot_pipeline(rice_pipeline_output_folder):
 def test_older_monocot_pipeline(rice_10do_pipeline_output_folder):
     # Find slp paths in folder
     slp_paths = find_all_slp_paths(rice_10do_pipeline_output_folder)
-    assert len(slp_paths) == 1
+    assert len(slp_paths) == 2
     # Load series from slps
     rice_series_all = load_series_from_slps(
         slp_paths=slp_paths, h5s=False, csv_path=None
     )
-    assert len(rice_series_all) == 1
+    assert len(rice_series_all) == 2
     # Get first series
     rice_series = rice_series_all[0]
-    assert rice_series.series_name == "scan0K9E8BI"
 
     pipeline = OlderMonocotPipeline()
     all_traits = pipeline.compute_batch_traits(rice_series_all)
     rice_traits = pipeline.compute_plant_traits(rice_series)
 
     # Dataframe shape assertions
     assert rice_traits.shape == (72, 102)
-    assert all_traits.shape == (1, 901)
+    assert all_traits.shape == (2, 901)
 
     # Dataframe dtype assertions
     expected_rice_traits_dtypes = {
@@ -223,20 +222,20 @@ def test_older_monocot_pipeline(rice_10do_pipeline_output_folder):
 
     # Value range assertions for traits
     assert (
-        all_traits["crown_curve_indices_mean_median"] >= 0
+        all_traits["crown_curve_indices_mean_median"].dropna() >= 0
     ).all(), "crown_curve_indices_mean_median in all_traits contains negative values"
     assert (
-        all_traits["crown_curve_indices_mean_median"] <= 1
+        all_traits["crown_curve_indices_mean_median"].dropna() <= 1
     ).all(), (
         "crown_curve_indices_mean_median in all_traits contains values greater than 1"
     )
     assert (
-        (0 <= rice_traits["crown_angles_proximal_p95"])
-        & (rice_traits["crown_angles_proximal_p95"] <= 180)
+        (0 <= rice_traits["crown_angles_proximal_p95"].dropna())
+        & (rice_traits["crown_angles_proximal_p95"].dropna() <= 180)
     ).all(), "angle_column in rice_traits contains values out of range [0, 180]"
     assert (
-        (0 <= all_traits["crown_angles_proximal_median_p95"])
-        & (all_traits["crown_angles_proximal_median_p95"] <= 180)
+        (0 <= all_traits["crown_angles_proximal_median_p95"].dropna())
+        & (all_traits["crown_angles_proximal_median_p95"].dropna() <= 180)
     ).all(), "angle_column in all_traits contains values out of range [0, 180]"