Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add convex hull function #13

Merged
merged 4 commits into from
Jan 19, 2023
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions sleap_roots/convhull.py
Original file line number Diff line number Diff line change
@@ -1 +1,56 @@
"""Convex hull fitting and derived trait calculation."""

import numpy as np
from scipy.spatial import ConvexHull, convex_hull_plot_2d
from scipy.spatial.distance import pdist
from typing import Tuple


def get_convhull(
Copy link
Contributor

@talmo talmo Jan 10, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is maybe a semantic preference, but I'd call this get_convhull_features.

The convex hull is technically defined by a set of points, so that's what I'd expect to get back from calling a function named get_convhull.

A better organization here might be to split this into two functions, where one just gets the convhull points:

def get_convhull(pts: np.ndarray) -> Optional[ConvexHull]:
    pts = pts.reshape(-1, 2)
    pts = pts[~(np.isnan(pts).any(axis=-1))]

    if len(pts) <= 2:
        return None

    # Get Convex Hull
    hull = ConvexHull(pts)

    return hull

And the second one which takes the ConvexHull instance and gives you back features. I appreciate it's more cumbersome to have to call two functions, so you could add a convenience checker to compute the convhull if needed:

def get_convhull_features(pts: Union[np.ndarray, ConvexHull]):
    hull = pts if type(pts) == ConvexHull else get_convhull(pts)

    if hull is None:
        return np.full((7,), np.nan)
    
    # perimeter
    perimeter = hull.area
    ...

pts: np.ndarray,
) -> Tuple[float, float, float, float, float, float, float]:
"""Get the convex hull for the points per frame.

Args:
pts: Root landmarks as array of shape (..., 2).

Returns:
A tuple of (perimeters, areas, longest_dists, shortest_dists, median_dists,
max_widths, max_heights) containing perimeter, area, longest distance between
vertices, smallest distance between vertices, median distance between vertices,
maximum width and maximum height of the convex hull.

If the convex hull fitting fails, NaNs are returned.
"""
pts = pts.reshape(-1, 2)
pts = pts[~(np.isnan(pts).any(axis=-1))]

if len(pts) <= 2:
return np.full((7,), np.nan)

# Get Convex Hull
hull = ConvexHull(pts)
# perimeter
perimeter = hull.area
# area
area = hull.volume
# longest distance between vertices
longest_dist = np.nanmax(pdist(hull.points[hull.vertices], "euclidean"))
# smallest distance between vertices
shortest_dist = np.nanmin(pdist(hull.points[hull.vertices], "euclidean"))
# median distance between vertices
median_dist = np.nanmedian(pdist(hull.points[hull.vertices], "euclidean"))
# max 'width'
max_width = np.nanmax(pts[:, 0]) - np.nanmin(pts[:, 0])
# max 'height'
max_height = np.nanmax(pts[:, 1]) - np.nanmin(pts[:, 1])

return (
perimeter,
area,
longest_dist,
shortest_dist,
median_dist,
max_width,
max_height,
)
143 changes: 143 additions & 0 deletions tests/test_convhull.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
from sleap_roots import Series
from sleap_roots.convhull import get_convhull
import numpy as np
import pytest


@pytest.fixture
def pts_nan31_5node():
return np.array(
[
[
[852.17755127, 216.95648193],
[np.nan, np.nan],
[np.nan, np.nan],
[np.nan, np.nan],
[816.71142578, 808.12585449],
],
[
[852.17755127, 216.95648193],
[np.nan, np.nan],
[837.03405762, 588.5123291],
[828.87963867, 692.72009277],
[816.71142578, 808.12585449],
],
]
)


@pytest.fixture
def pts_nan_5node():
return np.array(
[
[
[852.17755127, 216.95648193],
[np.nan, np.nan],
[np.nan, np.nan],
[np.nan, np.nan],
[816.71142578, 808.12585449],
],
]
)


# test canola model
def test_get_convhull_canola(canola_h5):
series = Series.load(
canola_h5, primary_name="primary_multi_day", lateral_name="lateral_3_nodes"
)
primary, lateral = series[0]

primary_points = primary.numpy().reshape(-1, 2)
lateral_points = lateral.numpy().reshape(-1, 2)
convex_hull_points = np.concatenate((primary_points, lateral_points), axis=0)

(
perimeters,
areas,
longest_dists,
shortest_dists,
median_dists,
max_widths,
max_heights,
) = get_convhull(convex_hull_points)

np.testing.assert_almost_equal(perimeters, 1910.0476127930017, decimal=3)
np.testing.assert_almost_equal(areas, 93255.32153574759, decimal=3)
np.testing.assert_almost_equal(longest_dists, 884.6450178192455, decimal=3)
np.testing.assert_almost_equal(shortest_dists, 185.15460001685398, decimal=3)
np.testing.assert_almost_equal(median_dists, 404.77175083902506, decimal=3)
np.testing.assert_almost_equal(max_widths, 211.279296875, decimal=3)
np.testing.assert_almost_equal(max_heights, 876.5622253417969, decimal=3)


# test rice model
def test_get_convhull_rice(rice_h5):
series = Series.load(
rice_h5, primary_name="main_3do_6nodes", lateral_name="longest_3do_6nodes"
)
primary, lateral = series[0]

primary_points = primary.numpy().reshape(-1, 2)
lateral_points = lateral.numpy().reshape(-1, 2)
convex_hull_points = np.concatenate((primary_points, lateral_points), axis=0)

(
perimeters,
areas,
longest_dists,
shortest_dists,
median_dists,
max_widths,
max_heights,
) = get_convhull(convex_hull_points)

np.testing.assert_almost_equal(perimeters, 1458.8585933576614, decimal=3)
np.testing.assert_almost_equal(areas, 23878.72090798154, decimal=3)
np.testing.assert_almost_equal(longest_dists, 720.0494276295367, decimal=3)
np.testing.assert_almost_equal(shortest_dists, 139.99063366108192, decimal=3)
np.testing.assert_almost_equal(median_dists, 480.1270688506847, decimal=3)
np.testing.assert_almost_equal(max_widths, 64.4229736328125, decimal=3)
np.testing.assert_almost_equal(max_heights, 720.0375061035156, decimal=3)


# test plant with 2 roots/instances with nan nodes
def test_get_convhull_nan(pts_nan31_5node):
(
perimeters,
areas,
longest_dists,
shortest_dists,
median_dists,
max_widths,
max_heights,
) = get_convhull(pts_nan31_5node)

np.testing.assert_almost_equal(perimeters, 1184.6684128638494, decimal=3)
np.testing.assert_almost_equal(areas, 2276.1159928281368, decimal=3)
np.testing.assert_almost_equal(longest_dists, 592.2322796929061, decimal=3)
np.testing.assert_almost_equal(shortest_dists, 104.52632471064256, decimal=3)
np.testing.assert_almost_equal(median_dists, 296.20807552792064, decimal=3)
np.testing.assert_almost_equal(max_widths, 35.46612548999997, decimal=3)
np.testing.assert_almost_equal(max_heights, 591.16937256, decimal=3)


# test plant with 1 root/instance with only 2 non-nan nodes
def test_get_convhull_nanall(pts_nan_5node):
(
perimeters,
areas,
longest_dists,
shortest_dists,
median_dists,
max_widths,
max_heights,
) = get_convhull(pts_nan_5node)

np.testing.assert_almost_equal(perimeters, np.nan, decimal=3)
np.testing.assert_almost_equal(areas, np.nan, decimal=3)
np.testing.assert_almost_equal(longest_dists, np.nan, decimal=3)
np.testing.assert_almost_equal(shortest_dists, np.nan, decimal=3)
np.testing.assert_almost_equal(median_dists, np.nan, decimal=3)
np.testing.assert_almost_equal(max_widths, np.nan, decimal=3)
np.testing.assert_almost_equal(max_heights, np.nan, decimal=3)