diff --git a/tests/conftest.py b/tests/conftest.py index 540470303cf..d535c6f797d 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -4,10 +4,13 @@ from pathlib import Path import numpy as np +import numpy.typing as npt import pandas as pd import pySPM import pytest import yaml +from skimage import draw, filters +from skimage.morphology import skeletonize import topostats from topostats.filters import Filters @@ -164,32 +167,32 @@ def plotting_config_with_plot_dict(default_config: dict) -> dict: @pytest.fixture() -def image_random() -> np.ndarray: +def image_random() -> npt.NDArray: """Random image as NumPy array.""" rng = np.random.default_rng(seed=1000) return rng.random((1024, 1024)) @pytest.fixture() -def small_array() -> np.ndarray: +def small_array() -> npt.NDArray: """Small (10x10) image array for testing.""" return RNG.random(SMALL_ARRAY_SIZE) @pytest.fixture() -def small_mask() -> np.ndarray: +def small_mask() -> npt.NDArray: """Small (10x10) mask array for testing.""" return RNG.uniform(low=0, high=1, size=SMALL_ARRAY_SIZE) > 0.5 @pytest.fixture() -def synthetic_scars_image() -> np.ndarray: +def synthetic_scars_image() -> npt.NDArray: """Small synthetic image for testing scar removal.""" return np.load(RESOURCES / "test_scars_synthetic_scar_image.npy") @pytest.fixture() -def synthetic_marked_scars() -> np.ndarray: +def synthetic_marked_scars() -> npt.NDArray: """Small synthetic boolean array of marked scar coordinates corresponding to synthetic_scars_image.""" return np.load(RESOURCES / "test_scars_synthetic_mark_scars.npy") @@ -777,7 +780,7 @@ def minicircle_all_statistics() -> pd.DataFrame: # Skeletonizing Fixtures @pytest.fixture() -def skeletonize_circular() -> np.ndarray: +def skeletonize_circular() -> npt.NDArray: """A circular molecule for testing skeletonizing.""" return np.array( [ @@ -807,13 +810,13 @@ def skeletonize_circular() -> np.ndarray: @pytest.fixture() -def skeletonize_circular_bool_int(skeletonize_circular: np.ndarray) -> np.ndarray: +def skeletonize_circular_bool_int(skeletonize_circular: np.ndarray) -> npt.NDArray: """A circular molecule for testing skeletonizing as a boolean integer array.""" return np.array(skeletonize_circular, dtype="bool").astype(int) @pytest.fixture() -def skeletonize_linear() -> np.ndarray: +def skeletonize_linear() -> npt.NDArray: """A linear molecule for testing skeletonizing.""" return np.array( [ @@ -846,9 +849,116 @@ def skeletonize_linear() -> np.ndarray: @pytest.fixture() -def skeletonize_linear_bool_int(skeletonize_linear) -> np.ndarray: +def skeletonize_linear_bool_int(skeletonize_linear) -> npt.NDArray: """A linear molecule for testing skeletonizing as a boolean integer array.""" return np.array(skeletonize_linear, dtype="bool").astype(int) -# Curvature Fixtures +# Pruning and Height profile fixtures +# +# Skeletons are generated by... +# +# 1. Generate random boolean images using scikit-image. +# 2. Skeletonize these shapes (gives boolean skeletons), these are our targets +# 3. Scale the skeletons by a factor (100) +# 4. Apply Gaussian filter to blur the heights and give an example original im +# for. + + +def _generate_heights(skeleton: npt.NDArray, scale: float = 100, sigma: float = 5.0, cval: float = 20.0) -> npt.NDArray: + """Generate heights from skeletons by scaling image and applying Gaussian blurring. + + Uses scikit-image 'skimage.filters.gaussian()' to generate heights from skeletons. + + Parameters + ---------- + skeleton : npt.NDArray + Binary array of skeleton. + scale : float + Factor to scale heights by. Boolean arrays are 0/1 and so the factor will be the height of the skeleton ridge. + sigma : float + Standard deviation for Gaussian kernel passed to `skimage.filters.gaussian()'. + cval : float + Value to fill past edges of input, passed to `skimage.filters.gaussian()'. + + Returns + ------- + npt.NDArray + Array with heights of image based on skeleton which will be the backbone and target. + """ + return filters.gaussian(skeleton * scale, sigma=sigma, cval=cval) + + +def _generate_random_skeleton(**extra_kwargs): + """Generate random skeletons and heights using skimage.draw's random_shapes().""" + kwargs = { + "image_shape": (128, 128), + "max_shapes": 20, + "channel_axis": None, + "shape": None, + "allow_overlap": True, + } + heights = {"scale": 100, "sigma": 5.0, "cval": 20.0} + random_image, _ = draw.random_shapes(**kwargs, **extra_kwargs) + mask = random_image != 255 + skeleton = skeletonize(mask) + return {"img": _generate_heights(skeleton, **heights), "skeleton": skeleton} + + +@pytest.fixture() +# def pruning_skeleton_loop1(heights=heights) -> dict: +def skeleton_loop1() -> dict: + """Skeleton with loop to be retained and side-branches.""" + return _generate_random_skeleton(rng=1, min_size=20) + + +@pytest.fixture() +def skeleton_loop2() -> dict: + """Skeleton with loop to be retained and side-branches.""" + return _generate_random_skeleton(rng=165103, min_size=60) + + +@pytest.fixture() +def skeleton_linear1() -> dict: + """Linear skeleton with lots of large side-branches, some forked.""" + return _generate_random_skeleton(rng=13588686514, min_size=20) + + +@pytest.fixture() +def skeleton_linear2() -> dict: + """Linear Skeleton with simple fork at one end.""" + return _generate_random_skeleton(rng=21, min_size=20) + + +@pytest.fixture() +def skeleton_linear3() -> dict: + """Linear Skeletons (i.e. multiple) with branches.""" + return _generate_random_skeleton(rng=894632511, min_size=20) + + +# Helper functions for visualising skeletons and heights +# +# def pruned_plot(gen_shape: dict) -> None: +# """Plot the original skeleton, its derived height and the pruned skeleton.""" +# img_skeleton = gen_shape() +# pruned = topostatsPrune( +# img_skeleton["heights"], +# img_skeleton["skeleton"], +# max_length=-1, +# height_threshold=90, +# method_values="min", +# method_outlier="abs", +# ) +# pruned_skeleton = pruned._prune_by_length(pruned.skeleton, pruned.max_length) +# fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2) +# ax1.imshow(img_skeleton["skeleton"]) +# ax2.imshow(img_skeleton["heights"]) +# ax3.imshow(pruned_skeleton) +# plt.show() + + +# pruned_plot(skeleton_loop1) +# pruned_plot(skeleton_loop2) +# pruned_plot(skeleton_linear1) +# pruned_plot(skeleton_linear2) +# pruned_plot(skeleton_linear3) diff --git a/tests/measure/test_height_profiles.py b/tests/measure/test_height_profiles.py new file mode 100644 index 00000000000..13b366ca9b4 --- /dev/null +++ b/tests/measure/test_height_profiles.py @@ -0,0 +1,449 @@ +"""Tests for height profiles.""" + +from __future__ import annotations + +import numpy as np +import numpy.typing as npt +import pytest + +from topostats.measure import height_profiles + +# pylint: disable=too-many-lines + + +@pytest.mark.parametrize( + ("img", "skeleton", "interpolate_conf", "target"), + [ + pytest.param( + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + {"method": "linear"}, + np.asarray([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]), + id="Basic circle, interpolation : linear", + ), + pytest.param( + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + {"method": "linear"}, + np.asarray([6, 7, 8, 9, 10, 9, 8, 7, 6]), + id="Diagonal line, interpolation : linear", + ), + pytest.param( + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + {"method": "cubic"}, + np.asarray( + [ + 6.00003876621563, + 7.000004717194621, + 7.9999928244616285, + 8.99999298402839, + 9.999970980789172, + 8.999992984028388, + 7.999992824461623, + 7.000004717194617, + 6.0000387662156305, + ] + ), + id="Diagonal line, interpolation : cubic", + ), + pytest.param( + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ), + {"method": "quintic"}, + np.asarray( + [ + 5.999989640371591, + 7.000016244559416, + 8.00000614600215, + 8.999994848347821, + 9.999994179877497, + 8.999996962477201, + 8.000007380709196, + 7.000012865886609, + 5.999980338747783, + ] + ), + id="Diagonal line, interpolation : quintic", + ), + ], +) +def test_interpolate_height_profile_basic( + img: npt.NDArray, skeleton: npt.NDArray, interpolate_conf: dict, target: npt.NDArray +) -> None: + """Test interpolate_height_profile() with basic shapes and heights.""" + np.testing.assert_array_almost_equal( + height_profiles.interpolate_height_profile(img, skeleton, **interpolate_conf), target, decimal=18 + ) + + +@pytest.mark.parametrize( + ("img", "interpolate_conf", "target"), + [ + pytest.param( + "skeleton_loop1", + {"method": "linear"}, + np.asarray( + [ + 4.67065435e-19, + 5.27325236e-19, + 5.73001876e-19, + 6.00077154e-19, + 6.02620717e-19, + 5.84638629e-19, + 5.49900643e-19, + 5.01019148e-19, + 4.44454025e-19, + 3.84704618e-19, + 3.24727212e-19, + 2.71579881e-19, + 2.24030707e-19, + 1.82928477e-19, + 1.52009449e-19, + 1.29002796e-19, + 1.12210549e-19, + 1.02377000e-19, + 9.88337895e-20, + 9.86897569e-20, + 1.01606876e-19, + 1.08936285e-19, + 1.18189308e-19, + 1.29587354e-19, + 1.45118712e-19, + 1.64343073e-19, + 1.87514816e-19, + 2.16107520e-19, + 2.50541210e-19, + 2.89309484e-19, + 3.31579050e-19, + 3.77100234e-19, + 4.22750106e-19, + 4.67427278e-19, + 5.12015950e-19, + 5.56562119e-19, + 6.02714720e-19, + 6.52632142e-19, + 7.06722254e-19, + 7.64764648e-19, + 8.24194642e-19, + 8.77655946e-19, + 9.22918373e-19, + 9.55310478e-19, + 9.68676402e-19, + 9.65855291e-19, + 9.48983887e-19, + 9.21105154e-19, + 8.88051609e-19, + 8.54278066e-19, + 8.23290173e-19, + 7.97837145e-19, + 7.78036698e-19, + 7.64027816e-19, + 7.54952213e-19, + 7.49473349e-19, + 7.46876370e-19, + 7.45973399e-19, + 7.45904036e-19, + 7.46705469e-19, + 7.48013031e-19, + 7.49133335e-19, + 7.50512257e-19, + 7.51909422e-19, + 7.52391803e-19, + 7.51626216e-19, + 7.48594362e-19, + 7.41889223e-19, + 7.30464673e-19, + 7.12964755e-19, + 6.88595014e-19, + 6.58882186e-19, + 6.22972271e-19, + 5.81557780e-19, + 5.40406540e-19, + 4.97897592e-19, + 4.55159488e-19, + 4.16411088e-19, + 3.81464749e-19, + 3.49043423e-19, + 3.21295477e-19, + 3.00311339e-19, + 2.84251853e-19, + 2.75044356e-19, + 2.77171887e-19, + 2.90484914e-19, + 3.17730359e-19, + 3.62941157e-19, + 4.27182722e-19, + 5.09548575e-19, + 6.08585652e-19, + 7.20418643e-19, + 8.35361178e-19, + 9.44196778e-19, + 1.03821903e-18, + 1.10607841e-18, + 1.14204135e-18, + 1.14240973e-18, + 1.10717643e-18, + 1.04225817e-18, + 9.53125935e-19, + 8.44475011e-19, + 7.27831261e-19, + 6.10118898e-19, + ] + ), + id="Skeleton loop 1", + # marks=pytest.mark.skip(), + ), + pytest.param( + "skeleton_loop2", + {"method": "linear"}, + np.asarray([4.40393772e-19, 2.04340231e-19, 4.34636982e-19]), + id="Skeleton loop 2", + marks=pytest.mark.skip(reason="Looks weird, expect more than three points."), + ), + pytest.param( + "skeleton_linear3", + {"method": "linear"}, + np.asarray( + [ + 4.60623136e-19, + 5.10302665e-19, + 5.37095486e-19, + 5.37165591e-19, + 5.10075873e-19, + 4.58990375e-19, + 3.92579809e-19, + 3.19427408e-19, + 2.46754846e-19, + 1.80615293e-19, + 1.27149020e-19, + 8.58026832e-20, + 5.56100490e-20, + 3.52552496e-20, + 2.32366177e-20, + 1.81760266e-20, + 1.76398399e-20, + 2.07214078e-20, + 2.66782218e-20, + 3.59301348e-20, + 4.79535313e-20, + 6.15538262e-20, + 7.59615427e-20, + 9.06585187e-20, + 1.07280567e-19, + 1.23941969e-19, + 1.41789438e-19, + 1.63283229e-19, + 1.93542424e-19, + 2.38155898e-19, + 3.00928079e-19, + 3.86080676e-19, + 4.95454001e-19, + 6.28051594e-19, + 7.73088313e-19, + 9.17635389e-19, + 1.04659862e-18, + 1.14401194e-18, + 1.19613838e-18, + 1.20151215e-18, + 1.16362278e-18, + 1.09164268e-18, + 9.98631761e-19, + 8.99456795e-19, + 8.04297034e-19, + 7.19380362e-19, + 6.47489041e-19, + 5.89873525e-19, + 5.40925921e-19, + 4.97410584e-19, + 4.56471131e-19, + 4.16658215e-19, + 3.76537601e-19, + 3.35352015e-19, + 2.93957893e-19, + 2.53511995e-19, + 2.15748369e-19, + 1.81398395e-19, + 1.51714946e-19, + 1.27775160e-19, + 1.10442752e-19, + 1.00520565e-19, + 9.82528292e-20, + 1.04083463e-19, + 1.18546517e-19, + 1.42786089e-19, + 1.77983378e-19, + 2.25405257e-19, + 2.86182615e-19, + 3.60671679e-19, + 4.47865755e-19, + 5.43429420e-19, + 6.41590773e-19, + 7.35104148e-19, + 8.16529727e-19, + 8.79666645e-19, + 9.21252909e-19, + 9.41628849e-19, + 9.44526855e-19, + 9.34636518e-19, + 9.18033824e-19, + 8.99737937e-19, + 8.82842133e-19, + 8.66510367e-19, + 8.47388128e-19, + 8.27092889e-19, + 8.05502825e-19, + 7.83381220e-19, + 7.59121805e-19, + 7.38876243e-19, + 7.26220278e-19, + 7.22418772e-19, + 7.26808537e-19, + 7.36069916e-19, + 7.48213880e-19, + 7.58164506e-19, + 7.60680300e-19, + 7.48834881e-19, + 7.18989781e-19, + 6.71276589e-19, + 6.06946427e-19, + 5.29400316e-19, + ] + ), + id="Linear skeleton 3, row 81", + ), + ], +) +def test_interpolate_height_profile_images(img: dict, interpolate_conf: dict, target: npt.NDArray, request) -> None: + """Test interpolate_height_profile() with more realistic images.""" + _img = request.getfixturevalue(img) + interpolated_heights = height_profiles.interpolate_height_profile(_img["img"], _img["skeleton"], **interpolate_conf) + np.testing.assert_array_almost_equal(interpolated_heights, target, decimal=18) diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Four arrays of different length (one even, two odd).png b/tests/resources/img/height_profiles/test_plot_height_profiles_Four arrays of different length (one even, two odd).png new file mode 100644 index 00000000000..6ce21cdbbc1 Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Four arrays of different length (one even, two odd).png differ diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Single height profile.png b/tests/resources/img/height_profiles/test_plot_height_profiles_Single height profile.png new file mode 100644 index 00000000000..04a640ba9df Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Single height profile.png differ diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Three arrays of different length (one even, one odd).png b/tests/resources/img/height_profiles/test_plot_height_profiles_Three arrays of different length (one even, one odd).png new file mode 100644 index 00000000000..2341b548d6a Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Three arrays of different length (one even, one odd).png differ diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is even).png b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is even).png new file mode 100644 index 00000000000..93be6eedd18 Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is even).png differ diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is odd).png b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is odd).png new file mode 100644 index 00000000000..31176eba02e Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of different length (diff in length is odd).png differ diff --git a/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of same length.png b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of same length.png new file mode 100644 index 00000000000..93be6eedd18 Binary files /dev/null and b/tests/resources/img/height_profiles/test_plot_height_profiles_Two arrays of same length.png differ diff --git a/tests/test_plotting.py b/tests/test_plotting.py index 7ce3f91c6d1..449f6b9e331 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -1,8 +1,12 @@ """Tests for the plotting module.""" +from __future__ import annotations + from importlib import resources from pathlib import Path +import numpy as np +import numpy.typing as npt import pandas as pd import pytest import yaml @@ -10,7 +14,7 @@ import topostats from topostats.entry_point import entry_point -from topostats.plotting import TopoSum, toposum +from topostats.plotting import TopoSum, _pad_array, plot_height_profiles, toposum # pylint: disable=protected-access @@ -192,3 +196,82 @@ def test_plot_violin_multiple_directories(toposum_object_multiple_directories: T """Test plotting Kernel Density Estimate and Histogram for area with multiple images.""" fig, _ = toposum_object_multiple_directories.sns_violinplot() return fig + + +@pytest.mark.mpl_image_compare(baseline_dir="resources/img/height_profiles/") +@pytest.mark.parametrize( + ("height_profile"), + [ + pytest.param(np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), id="Single height profile"), + pytest.param( + [ + np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), + np.asarray([0, 0, 0, 2, 4, 5, 5, 5, 4, 2, 0, 0, 0]), + ], + id="Two arrays of same length", + ), + pytest.param( + [ + np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), + np.asarray([0, 0, 2, 4, 5, 5, 5, 4, 2, 0, 0]), + ], + id="Two arrays of different length (diff in length is even)", + ), + pytest.param( + [ + np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), + np.asarray([0, 0, 2, 4, 5, 5, 5, 4, 2, 0, 0, 0]), + ], + id="Two arrays of different length (diff in length is odd)", + ), + pytest.param( + [ + np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), + np.asarray([0, 0, 2, 4, 5, 5, 5, 4, 2, 0, 0]), + np.asarray([0, 0, 1, 5, 6, 7, 6, 5, 1, 0, 0, 0]), + ], + id="Three arrays of different length (one even, one odd)", + ), + pytest.param( + [ + np.asarray([0, 0, 0, 2, 3, 4, 4, 4, 3, 2, 0, 0, 0]), + np.asarray([0, 0, 2, 4, 5, 5, 5, 4, 2, 0, 0]), + np.asarray([0, 0, 1, 5, 6, 7, 6, 5, 1, 0, 0, 0]), + np.asarray([0, 0, 1, 4, 1, 0, 0]), + ], + id="Four arrays of different length (one even, two odd)", + ), + ], +) +def test_plot_height_profiles(height_profile: list | npt.NDArray) -> None: + """Test plotting of height profiles.""" + fig, _ = plot_height_profiles(height_profile) + return fig + + +@pytest.mark.parametrize( + ("arrays", "max_array_length", "target"), + [ + pytest.param( + np.asarray([1, 1, 1]), + 3, + np.asarray([1, 1, 1]), + id="Array length 3, max array length 3", + ), + pytest.param( + np.asarray([1, 1, 1]), + 5, + np.asarray([0, 1, 1, 1, 0]), + id="Array length 3, max array length 5", + ), + pytest.param( + np.asarray([1, 1, 1]), + 4, + np.asarray([1, 1, 1, 0]), + id="Array length 3, max array length 4", + ), + ], +) +def test_pad_array(arrays: list, max_array_length: int, target: list) -> None: + """Test padding of arrays.""" + np.testing.assert_array_equal(_pad_array(arrays, max_array_length), target) diff --git a/topostats/measure/height_profiles.py b/topostats/measure/height_profiles.py new file mode 100644 index 00000000000..a518a2cf114 --- /dev/null +++ b/topostats/measure/height_profiles.py @@ -0,0 +1,57 @@ +"""Derive height profiles across the images.""" + +from __future__ import annotations + +import logging +import warnings + +import numpy as np +import numpy.typing as npt +from scipy import interpolate + +from topostats.logs.logs import LOGGER_NAME +from topostats.measure import feret + +LOGGER = logging.getLogger(LOGGER_NAME) + +# Handle warnings as exceptions (encountered when gradient of base triangle is zero) +warnings.filterwarnings("error") + + +def interpolate_height_profile(img: npt.NDArray, skeleton: npt.NDArray, **kwargs) -> npt.NDArray: + """ + Interpolate heights along the maximum feret. + + Interpolates the height along the line of the maximum feret using SciPy `scipy.interpolateRegularGridInterpolator + https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html`. Arguments can + be passed using 'kwargs'. + + Parameters + ---------- + img : npt.NDArray + Original image with heights. + skeleton : npt.NDArray + Binary skeleton. + **kwargs : dict + Keyword arguments passed on to scipy.interpolate.RegularGridInterpolator(). + + Returns + ------- + npt.NDArray + Interpolated heights between the calculated feret co-cordinates. + """ + # Extract feret coordinates + feret_stats = feret.get_feret_from_mask(skeleton) + x_coords = np.sort(feret_stats["max_feret_coords"][:, 0]) + y_coords = np.sort(feret_stats["max_feret_coords"][:, 1]) + # Evenly spaced points on x-axis + x_points = np.linspace(np.min(x_coords), np.max(x_coords), np.max(x_coords) - np.min(x_coords) + 1) + # Corresponding y values for each x + y_points = np.interp(x_points, x_coords, y_coords) + # Combine into array + xy_points = np.vstack((x_points, y_points)).T + # Interpolate + interp = interpolate.RegularGridInterpolator( + points=(np.arange(img.shape[0]), np.arange(img.shape[1])), values=img, **kwargs + ) + return interp(xy_points) diff --git a/topostats/plotting.py b/topostats/plotting.py index 6bf51b7116c..1e702bae4f6 100644 --- a/topostats/plotting.py +++ b/topostats/plotting.py @@ -10,9 +10,10 @@ import sys import yaml import matplotlib.pyplot as plt +import numpy as np +import numpy.typing as npt import pandas as pd import seaborn as sns -import numpy as np from topostats.io import read_yaml, save_pkl, write_yaml, convert_basename_to_relative_paths from topostats.logs.logs import LOGGER_NAME @@ -473,3 +474,64 @@ def run_toposum(args=None) -> None: # Plot statistics toposum(config) + + +def plot_height_profiles(height_profiles: list | npt.NDArray) -> tuple: + """ + Plot height profiles. + + Parameters + ---------- + height_profiles : npt.NDArray + Single height profile (1-D numpy array of heights) or array of height profiles. If the later the profiles plot + will be overlaid. + + Returns + ------- + tuple + Matplotlib.pyplot figure object and Matplotlib.pyplot axes object. + """ + # If we have only one profile put it in a list so we can process + if not isinstance(height_profiles, list): + height_profiles = [height_profiles] + # We have 1-D arrays of different sizes, we need to know the maximum length and then pad shorter ones so that the + # profiles will roughly align in the middle + max_array_length = max(map(len, height_profiles)) + + # Pad shorter height profiles to the length of the longest + padded_height_profiles = [_pad_array(profile, max_array_length) for profile in height_profiles] + fig, ax = plt.subplots(1, 1) + max_y = 0 + for height_profile in padded_height_profiles: + ax.plot(np.arange(max_array_length), height_profile) + max_y = max_y if max(height_profile) < max_y else max(height_profile) + 1 + ax.margins(0.01, 0.1) + return fig, ax + + +def _pad_array(profile: npt.NDArray, max_array_length: int) -> npt.NDArray: + """ + Pad array so that it matches the largest profile and plots are somewhat aligned. + + Centering is done based on the mid-point of longest grain and heights of zero ('0.0') are used in padding. + + Parameters + ---------- + profile : npt.NDArray + 1-D Height profile. + max_array_length : int + The longest height profile across a range of detected grains. + + Returns + ------- + npt.NDArray + Array padded to the same length as max_array_length. + """ + profile_length = profile.shape[0] + array_diff = max_array_length - profile_length + pad = array_diff // 2 + if array_diff % 2 == 0: + left, right = pad, pad + else: + left, right = pad, (pad + 1) + return np.pad(profile, (left, right))