Handle unsupported geometry types (#91)

* add test data that breaks extract_points_from_geometry

* replace images filename in slp file

* update OlderMonocotPipeline notebook

* return empty list if unsupported geometry type

* update oldermonocot pipeline

* use logging statements for convex hull calculations

* modify tests for two plants

* check intersection for empty list, return nans

* add more test for unexpected inputs to extract_points_from_geometry

* include empty MultiLineString and remove duplicate test

* update OlderMonocot notebook

* add stunted 10 day-old rice predictions as fixture

* fix paths

* organize and add needed imports

* update OlderMonocot pipeline notebook

* black

* black

* test new fixture

* black

* black

* add more unexpected types

* use logging

* add test for inf values

* add more tests for multilinestring

* add tests for nex convex hull features using rice fixture

* black

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,20 +31,20 @@ def get_convhull(pts: np.ndarray) -> Optional[ConvexHull]:
 
     # 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:
         # Compute and return the convex hull
         return ConvexHull(unique_pts)
     except Exception as e:
-        print(f"Cannot compute convex hull: {e}")
+        logging.info(f"Cannot compute convex hull: {e}")
         return None
 
 
@@ -384,16 +385,27 @@ def get_chull_areas_via_intersection(
 
     # 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 +487,26 @@ def get_chull_intersection_vectors(
 
     # 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 +528,7 @@ def get_chull_intersection_vectors(
     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 +538,9 @@ def get_chull_intersection_vectors(
         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 +568,21 @@ def get_chull_intersection_vectors(
 
     # 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,8 @@ 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/fixtures/data.py

@@ -79,6 +79,12 @@ def rice_main_10do_slp():
     return "tests/data/rice_10do/0K9E8BI.crown.predictions.slp"
 
 
+@pytest.fixture
+def rice_10do_stunted_slp():
+    """Path to stunted root predictions for 10 day old rice."""
+    return "tests/data/rice_10do_pipeline_output/scan_7859150.model_221208_113552.multi_instance.n=574.root_crown.slp"
+
+
 @pytest.fixture
 def soy_folder():
     """Path to a folder with the predictions for 6 day old soy."""

tests/test_convhull.py

@@ -1,4 +1,9 @@
 from scipy.spatial import ConvexHull
+
+import numpy as np
+import pytest
+import logging
+
 from sleap_roots import Series
 from sleap_roots.convhull import (
     get_convhull,
@@ -10,11 +15,15 @@
     get_chull_division_areas,
     get_chull_areas_via_intersection,
     get_chull_intersection_vectors,
+    get_chull_intersection_vectors_left,
+    get_chull_intersection_vectors_right,
+    get_chull_area_via_intersection_above,
+    get_chull_area_via_intersection_below,
 )
 from sleap_roots.lengths import get_max_length_pts
 from sleap_roots.points import get_all_pts_array, get_nodes
-import numpy as np
-import pytest
+from sleap_roots.bases import get_bases
+from sleap_roots.angle import get_vector_angles_from_gravity
 
 
 @pytest.fixture
@@ -28,6 +37,14 @@ def valid_input():
     return rn_pts, pts, hull, (expected_area_above, expected_area_below)
 
 
+@pytest.fixture
+def invalid_input_with_inf_values():
+    pts = np.array(
+        [[[0, 0], [2, 2], [4, 0], [2, -2], [0, -4], [4, -4], [np.inf, np.inf]]]
+    )
+    return pts
+
+
 @pytest.fixture
 def invalid_pts_shape():
     rn_pts = np.array([[0, 0], [1, 1]])
@@ -120,6 +137,18 @@ def lateral_pts():
     )
 
 
+# test get_convhull with invalid input with inf values
+def test_infinite_values_logging(invalid_input_with_inf_values, caplog):
+    with caplog.at_level(logging.INFO):  # This message is INFO level
+        _ = get_convhull(invalid_input_with_inf_values)
+
+    # Check that the specific log message is in caplog
+    assert any(
+        "Cannot compute convex hull: input contains infinite values." in message
+        for message in caplog.messages
+    )
+
+
 # test get_convhull function using canola
 def test_get_convhull_canola(canola_h5, canola_primary_slp, canola_lateral_slp):
     # Set frame index to 0
@@ -188,17 +217,48 @@ def test_get_convhull_features_rice(rice_h5, rice_long_slp, rice_main_slp):
     )
     # Get the crown root from the series
     crown_pts = series.get_crown_points(frame_index)
+
+    # Get the r0 and r1 nodes from the crown root
+    r0_pts = get_bases(crown_pts)
+    r1_pts = get_nodes(crown_pts, 1)
+
     # Get the convex hull from the points
     convex_hull = get_convhull(crown_pts)
     perimeter = get_chull_perimeter(convex_hull)
     area = get_chull_area(convex_hull)
     max_width = get_chull_max_width(convex_hull)
     max_height = get_chull_max_height(convex_hull)
 
+    # Get the intersection vectors
+    left_vector = get_chull_intersection_vectors_left(
+        get_chull_intersection_vectors(r0_pts, r1_pts, crown_pts, convex_hull)
+    )
+    right_vector = get_chull_intersection_vectors_right(
+        get_chull_intersection_vectors(r0_pts, r1_pts, crown_pts, convex_hull)
+    )
+    # Get angles from gravity
+    left_angle = get_vector_angles_from_gravity(left_vector)
+    right_angle = get_vector_angles_from_gravity(right_vector)
+
+    # Get the intersection areas
+    area_above = get_chull_area_via_intersection_above(
+        get_chull_areas_via_intersection(r1_pts, crown_pts, convex_hull)
+    )
+    area_below = get_chull_area_via_intersection_below(
+        get_chull_areas_via_intersection(r1_pts, crown_pts, convex_hull)
+    )
+
+    assert left_vector.shape == (1, 2)
+    assert right_vector.shape == (1, 2)
+
     np.testing.assert_almost_equal(perimeter, 1450.6365795858003, decimal=3)
     np.testing.assert_almost_equal(area, 23722.883102604676, decimal=3)
     np.testing.assert_almost_equal(max_width, 64.341064453125, decimal=3)
     np.testing.assert_almost_equal(max_height, 715.6949920654297, decimal=3)
+    np.testing.assert_almost_equal(left_angle, 166.08852115310046, decimal=3)
+    np.testing.assert_almost_equal(right_angle, 0.04343543279020469, decimal=3)
+    np.testing.assert_almost_equal(area_above, 1903.040353098403, decimal=3)
+    np.testing.assert_almost_equal(area_below, 21819.842749506273, decimal=3)
 
 
 # test plant with 2 roots/instances with nan nodes
@@ -362,3 +422,51 @@ def test_no_convex_hull():
     assert np.isnan(
         right_vector
     ).all(), "Expected NaN vector for right_vector when hull is None"
+
+
+# Test with stunted 10 day-old rice sample that gave `MultiLineString` geometry type as r1 intersection with convex hull
+def test_stunted_rice_10do_multilinestring(rice_10do_stunted_slp):
+    # Load the series from the rice dataset
+    series = Series.load(
+        series_name="stunted_rice_test",
+        crown_path=rice_10do_stunted_slp,
+    )
+    # Get the crown roots from the series
+    crown_pts = series.get_crown_points(3)  # Frame index 3 was the problematic frame
+    # Get the convex hull from the points
+    convex_hull = get_convhull(crown_pts)
+    r0_pts = get_bases(crown_pts)
+    r1_pts = get_nodes(crown_pts, 1)
+    left_vector, right_vector = get_chull_intersection_vectors(
+        r0_pts, r1_pts, crown_pts, convex_hull
+    )
+    # Expected vectors are NaN since the intersection is a `MultiLineString` geometry type
+    # Check if both left_vector and right_vector are all NaNs
+    assert np.all(
+        np.isnan(left_vector)
+    ), "Left vector does not contain all NaNs as expected."
+    assert np.all(
+        np.isnan(right_vector)
+    ), "Right vector does not contain all NaNs as expected."
+
+
+# Similarly, `get_chull_intersection_areas` can be tested with the same stunted 10 day-old rice sample
+# to check if the areas are NaN when the intersection is a `MultiLineString` geometry type
+def test_stunted_rice_10do_multilinestring_areas(rice_10do_stunted_slp):
+    # Load the series from the rice dataset
+    series = Series.load(
+        series_name="stunted_rice_test",
+        crown_path=rice_10do_stunted_slp,
+    )
+    # Get the crown roots from the series
+    crown_pts = series.get_crown_points(3)  # Frame index 3 was the problematic frame
+    # Get the convex hull from the points
+    convex_hull = get_convhull(crown_pts)
+    r1_pts = get_nodes(crown_pts, 1)
+    area_above, area_below = get_chull_areas_via_intersection(
+        r1_pts, crown_pts, convex_hull
+    )
+    # Expected areas are NaN since the intersection is a `MultiLineString` geometry type
+    # Check if both area_above and area_below are NaNs
+    assert np.isnan(area_above), "Area above line does not contain NaN as expected."
+    assert np.isnan(area_below), "Area below line does not contain NaN as expected."