Repository: talmolab/sleap-roots (language: python)

Issue:
Trait: Network Length Distribution (Ldist)
**GiA Trait**: Network Length Distribution (Ldist)
**Units**: n/n
**Description**: The fraction of network pixels found in the lower 2/3 of the network. The lower 2/3 of the network is defined based on the network depth.
**Analogous SLEAP Trait**: fraction of total lengths

**Implementation**:
```py
def get_network_length_fraction(points: np.ndarray, fraction: float = 2/3) -> float:
    """Return fraction of root network length that is in the lower fraction of the plant."""
```

This could be achieved by:

1. Find the bounding box of all the plant points --> `bbox = (left_x, top_y, width, height)`.
2. Find the bounding box of the lower fraction:
    ```py
    lower_height = bbox[3] * fraction
    lower_bbox = (bbox[0], bbox[1] - (bbox[3] - lower_height), bbox[2], lower_height)
    ```
3. Intersect `lower_bbox` with the points (tip: use [shapely](https://shapely.readthedocs.io/en/stable/manual.html#collections))
4. For every segment, compute length, sum the total --> `lower_length`
5. Subtract `lower_length` from total length


Discussion / hints:

Fix patch:
diff --git a/sleap_roots/__init__.py b/sleap_roots/__init__.py
index 08cb5ee..c3891bc 100644
--- a/sleap_roots/__init__.py
+++ b/sleap_roots/__init__.py
@@ -4,6 +4,7 @@
 import sleap_roots.bases
 import sleap_roots.convhull
 import sleap_roots.ellipse
+import sleap_roots.networklength
 import sleap_roots.scanline
 import sleap_roots.series
 from sleap_roots.series import Series
diff --git a/sleap_roots/networklength.py b/sleap_roots/networklength.py
new file mode 100644
index 0000000..e202872
--- /dev/null
+++ b/sleap_roots/networklength.py
@@ -0,0 +1,162 @@
+"""Fraction of root network length in the lower fraction of the plant."""
+
+import numpy as np
+from shapely import LineString, Polygon
+from sleap_roots.bases import get_root_lengths
+from sleap_roots.convhull import get_convhull_features
+from typing import Tuple
+
+
+def get_bbox(pts: np.ndarray) -> Tuple[float, float, float, float]:
+    """Return the bounding box of all landmarks.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+
+    Returns:
+        Tuple of four parameters in bounding box:
+        left_x, the x axis value of left side
+        top_y, the y axis value of top side
+        width, the width of the bounding box
+        height, the height of bounding box.
+    """
+    # reshape to (# instance, 2) and filter out NaNs.
+    pts2 = pts.reshape(-1, 2)
+    pts2 = pts2[~(np.isnan(pts2).any(axis=-1))]
+
+    # get the bounding box
+    left_x, top_y = np.min(pts2[:, 0]), np.min(pts2[:, 1])
+    width, height = np.max(pts2[:, 0]) - np.min(pts2[:, 0]), np.max(
+        pts2[:, 1]
+    ) - np.min(pts2[:, 1])
+    bbox = (left_x, top_y, width, height)
+    return bbox
+
+
+def get_network_width_depth_ratio(pts: np.ndarray) -> float:
+    """Return width to depth ratio of bounding box for root network.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+
+    Returns:
+        float of bounding box width to depth ratio of root network.
+    """
+    # get the bounding box
+    bbox = get_bbox(pts)
+    width, height = bbox[2], bbox[3]
+    if width > 0 and height > 0:
+        ratio = width / height
+        return ratio
+    else:
+        return np.nan
+
+
+def get_network_solidity(pts: np.ndarray) -> float:
+    """Return the total network area divided by the network convex area.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+
+    Returns:
+        float of the total network area divided by the network convex area.
+    """
+    # get the bounding box area
+    bbox = get_bbox(pts)
+    width, height = bbox[2], bbox[3]
+    bbox_area = width * height
+
+    # get the convex hull area
+    convhull_features = get_convhull_features(pts)
+    conv_area = convhull_features[1]
+
+    if bbox_area > 0 and conv_area > 0:
+        ratio = bbox_area / conv_area
+        return ratio
+    else:
+        return np.nan
+
+
+def get_network_distribution(pts: np.ndarray, fraction: float = 2 / 3) -> float:
+    """Return the root length in the lower fraction of the plant.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+        fraction: the network length found in the lower fration value of the network.
+
+    Returns:
+        float of the root network length in the lower fraction of the plant.
+    """
+    # get the bounding box
+    bbox = get_bbox(pts)
+    left_x, top_y, width, height = bbox[0], bbox[1], bbox[2], bbox[3]
+
+    # get the bounding box of the lower fraction
+    lower_height = bbox[3] * fraction
+    lower_bbox = (bbox[0], bbox[1] + (bbox[3] - lower_height), bbox[2], lower_height)
+
+    # convert bounding box to polygon
+    polygon = Polygon(
+        [
+            [bbox[0], bbox[1] + (bbox[3] - lower_height)],
+            [bbox[0], bbox[1] + height],
+            [bbox[0] + width, bbox[1] + height],
+            [bbox[0] + width, bbox[1] + (bbox[3] - lower_height)],
+        ]
+    )
+
+    # filter out the nan nodes
+    points = list(pts)
+    pts_nnan = []
+    for j in range(len(points)):
+        pts_j = points[j][~np.isnan(points[j]).any(axis=1)]
+        pts_nnan.append(pts_j)
+
+    # get length of lines within the lower bounding box
+    root_length = 0
+    for j in range(len(points)):
+        # filter out nan nodes
+        pts_j = points[j][~np.isnan(points[j]).any(axis=1)]
+        linestring = LineString(points[j])
+        if pts_j.shape[0] > 1:
+            if linestring.intersection(polygon):
+                intersection = linestring.intersection(polygon)
+                root_length += (
+                    intersection.length if ~np.isnan(intersection.length) else 0
+                )
+
+    return root_length
+
+
+def get_network_length(pts: np.ndarray) -> float:
+    """Return all primary or lateral root length one frame.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+
+    Returns:
+        float of primary or lateral root network length.
+    """
+    if np.sum(get_root_lengths(pts)) > 0:
+        length = np.sum(get_root_lengths(pts))
+        return length
+    else:
+        return 0
+
+
+def get_network_distribution_ratio(pts: np.ndarray, fraction: float = 2 / 3) -> float:
+    """Return ratio of the root length in the lower fraction over all root length.
+
+    Args:
+        pts: Root landmarks as array of shape (..., 2).
+        fraction: the network length found in the lower fration value of the network.
+
+    Returns:
+        float of ratio of the root network length in the lower fraction of the plant
+        over all root length.
+    """
+    if np.sum(get_root_lengths(pts)) > 0:
+        ratio = get_network_distribution(pts, fraction) / np.sum(get_root_lengths(pts))
+        return ratio
+    else:
+        return np.nan


Test patch:
diff --git a/tests/test_networklength.py b/tests/test_networklength.py
new file mode 100644
index 0000000..8b82555
--- /dev/null
+++ b/tests/test_networklength.py
@@ -0,0 +1,190 @@
+import pytest
+import numpy as np
+from sleap_roots import Series
+from sleap_roots.networklength import get_bbox
+from sleap_roots.networklength import get_network_distribution
+from sleap_roots.networklength import get_network_distribution_ratio
+from sleap_roots.networklength import get_network_length
+from sleap_roots.networklength import get_network_solidity
+from sleap_roots.networklength import get_network_width_depth_ratio
+
+
+@pytest.fixture
+def pts_nan3():
+    return np.array(
+        [
+            [
+                [852.17755127, 216.95648193],
+                [np.nan, 472.83520508],
+                [844.45300293, np.nan],
+            ]
+        ]
+    )
+
+
+def test_get_bbox(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    bbox = get_bbox(pts)
+    np.testing.assert_almost_equal(
+        bbox, [1016.7844238, 144.4191589, 192.1080322, 876.5622253], decimal=7
+    )
+
+
+def test_get_bbox_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    bbox = get_bbox(pts)
+    np.testing.assert_almost_equal(
+        bbox, [796.2611694, 248.6078033, 64.3410645, 715.6949921], decimal=7
+    )
+
+
+def test_get_bbox_nan(pts_nan3):
+    pts = pts_nan3
+    bbox = get_bbox(pts)
+    np.testing.assert_almost_equal(
+        bbox, [852.1775513, 216.9564819, 0.0, 0.0], decimal=7
+    )
+
+
+def test_get_network_width_depth_ratio(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    ratio = get_network_width_depth_ratio(pts)
+    np.testing.assert_almost_equal(ratio, 0.2191607471467916, decimal=7)
+
+
+def test_get_network_width_depth_ratio_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    ratio = get_network_width_depth_ratio(pts)
+    np.testing.assert_almost_equal(ratio, 0.0899001182996162, decimal=7)
+
+
+def test_get_network_width_depth_ratio_nan(pts_nan3):
+    pts = pts_nan3
+    ratio = get_network_width_depth_ratio(pts)
+    np.testing.assert_almost_equal(ratio, np.nan, decimal=7)
+
+
+def test_get_network_solidity(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    ratio = get_network_solidity(pts)
+    np.testing.assert_almost_equal(ratio, 2.0351966283618834, decimal=7)
+
+
+def test_get_network_solidity_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    ratio = get_network_solidity(pts)
+    np.testing.assert_almost_equal(ratio, 1.9411037610434734, decimal=7)
+
+
+def test_get_network_solidity_nan(pts_nan3):
+    pts = pts_nan3
+    ratio = get_network_solidity(pts)
+    np.testing.assert_almost_equal(ratio, np.nan, decimal=7)
+
+
+def test_get_network_distribution(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    fraction = 2 / 3
+    root_length = get_network_distribution(pts, fraction)
+    np.testing.assert_almost_equal(root_length, 589.4322131363684, decimal=7)
+
+
+def test_get_network_distribution_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    fraction = 2 / 3
+    root_length = get_network_distribution(pts, fraction)
+    np.testing.assert_almost_equal(root_length, 477.77168597561507, decimal=7)
+
+
+def test_get_network_distribution_nan(pts_nan3):
+    pts = pts_nan3
+    fraction = 2 / 3
+    root_length = get_network_distribution(pts, fraction)
+    np.testing.assert_almost_equal(root_length, 0, decimal=7)
+
+
+def test_get_network_length(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    length = get_network_length(pts)
+    np.testing.assert_almost_equal(length, 971.0504174567843, decimal=7)
+
+
+def test_get_network_length_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    length = get_network_length(pts)
+    np.testing.assert_almost_equal(length, 798.5726441151357, decimal=7)
+
+
+def test_get_network_length_nan(pts_nan3):
+    pts = pts_nan3
+    length = get_network_length(pts)
+    np.testing.assert_almost_equal(length, 0, decimal=7)
+
+
+def test_get_network_distribution_ratio(canola_h5):
+    series = Series.load(
+        canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    fraction = 2 / 3
+    ratio = get_network_distribution_ratio(pts, fraction)
+    np.testing.assert_almost_equal(ratio, 0.6070047, decimal=7)
+
+
+def test_get_network_distribution_ratio_rice(rice_h5):
+    series = Series.load(
+        rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
+    )
+    primary, lateral = series[0]
+    pts = primary.numpy()
+    fraction = 2 / 3
+    ratio = get_network_distribution_ratio(pts, fraction)
+    np.testing.assert_almost_equal(ratio, 0.5982820592421038, decimal=7)
+
+
+def test_get_network_distribution_ratio_nan(pts_nan3):
+    pts = pts_nan3
+    fraction = 2 / 3
+    ratio = get_network_distribution_ratio(pts, fraction)
+    np.testing.assert_almost_equal(ratio, np.nan, decimal=7)