Merge pull request #859 from AFM-SPM/ns-rse/748-height-profiles

feature(measure): Adds the height_profiles sub-module

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,114 @@ 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 image with heights
+
+
+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 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)

tests/measure/conftest.py

@@ -0,0 +1,167 @@
+"""Fixtures for testing sub-modules of the measure module."""
+
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+import pytest
+from skimage import draw
+
+# pylint: disable=redefined-outer-name
+
+
+@pytest.fixture()
+def tiny_circle() -> npt.NDArray:
+    """Tiny circle."""
+    tiny_circle = np.zeros((3, 3), dtype=np.uint8)
+    rr, cc = draw.circle_perimeter(1, 1, 1)
+    tiny_circle[rr, cc] = 1
+    return tiny_circle
+
+
+@pytest.fixture()
+def small_circle() -> npt.NDArray:
+    """Small circle."""
+    small_circle = np.zeros((5, 5), dtype=np.uint8)
+    rr, cc = draw.circle_perimeter(2, 2, 2)
+    small_circle[rr, cc] = 1
+    return small_circle
+
+
+@pytest.fixture()
+def tiny_quadrilateral() -> npt.NDArray:
+    """Tiny quadrilateral."""
+    return np.asarray(
+        [
+            [0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0],
+            [0, 1, 0, 0, 1, 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],
+        ],
+        dtype=np.uint8,
+    )
+
+
+@pytest.fixture()
+def tiny_square() -> npt.NDArray:
+    """Tiny square."""
+    return np.asarray([[0, 0, 0, 0], [0, 1, 1, 0], [0, 1, 1, 0], [0, 0, 0, 0]], dtype=np.uint8)
+
+
+@pytest.fixture()
+def tiny_triangle() -> npt.NDArray:
+    """Tiny triangle."""
+    return np.asarray([[0, 0, 0, 0], [0, 1, 1, 0], [0, 1, 0, 0], [0, 0, 0, 0]], dtype=np.uint8)
+
+
+@pytest.fixture()
+def arbitrary_triangle() -> npt.NDArray:
+    """Arbitrary triangle."""
+    return np.array([[1.18727719, 2.96140198], [4.14262995, 3.17983179], [6.92472119, 0.64147496]])
+
+
+@pytest.fixture()
+def tiny_rectangle() -> npt.NDArray:
+    """Tiny rectangle."""
+    return np.asarray([[0, 0, 0, 0], [0, 1, 1, 0], [0, 1, 1, 0], [0, 1, 1, 0], [0, 0, 0, 0]], dtype=np.uint8)
+
+
+@pytest.fixture()
+def tiny_ellipse() -> npt.NDArray:
+    """Tiny ellipse."""
+    return np.asarray(
+        [
+            [0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0],
+            [0, 1, 0, 1, 0],
+            [0, 1, 0, 1, 0],
+            [0, 1, 0, 1, 0],
+            [0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0],
+        ],
+        dtype=np.uint8,
+    )
+
+
+@pytest.fixture()
+def holo_circle() -> npt.NDArray:
+    """Holo circle."""
+    holo_circle = np.zeros((7, 7), dtype=np.uint8)
+    rr, cc = draw.circle_perimeter(3, 3, 2)
+    holo_circle[rr, cc] = 1
+    return holo_circle
+
+
+@pytest.fixture()
+def holo_ellipse_vertical() -> npt.NDArray:
+    """Holo ellipse (vertical)."""
+    holo_ellipse_vertical = np.zeros((11, 9), dtype=np.uint8)
+    rr, cc = draw.ellipse_perimeter(5, 4, 4, 3)
+    holo_ellipse_vertical[rr, cc] = 1
+    return holo_ellipse_vertical
+
+
+@pytest.fixture()
+def holo_ellipse_horizontal() -> npt.NDArray:
+    """Holo ellipse (horizontal)."""
+    holo_ellipse_horizontal = np.zeros((9, 11), dtype=np.uint8)
+    rr, cc = draw.ellipse_perimeter(4, 5, 3, 4)
+    holo_ellipse_horizontal[rr, cc] = 1
+    return holo_ellipse_horizontal
+
+
+@pytest.fixture()
+def holo_ellipse_angled() -> npt.NDArray:
+    """Holo ellipse (angled)."""
+    holo_ellipse_angled = np.zeros((8, 10), dtype=np.uint8)
+    rr, cc = draw.ellipse_perimeter(4, 5, 1, 3, orientation=np.deg2rad(30))
+    holo_ellipse_angled[rr, cc] = 1
+    return holo_ellipse_angled
+
+
+@pytest.fixture()
+def curved_line() -> npt.NDArray:
+    """Curved line."""
+    curved_line = np.zeros((10, 10), dtype=np.uint8)
+    rr, cc = draw.bezier_curve(1, 5, 5, -2, 8, 8, 2)
+    curved_line[rr, cc] = 1
+    return curved_line
+
+
+@pytest.fixture()
+def filled_circle() -> npt.NDArray:
+    """Circle."""
+    filled_circle = np.zeros((9, 9), dtype=np.uint8)
+    rr, cc = draw.disk((4, 4), 4)
+    filled_circle[rr, cc] = 1
+    return filled_circle
+
+
+@pytest.fixture()
+def filled_ellipse_vertical() -> npt.NDArray:
+    """Ellipse (vertical)."""
+    filled_ellipse_vertical = np.zeros((9, 7), dtype=np.uint8)
+    rr, cc = draw.ellipse(4, 3, 4, 3)
+    filled_ellipse_vertical[rr, cc] = 1
+    return filled_ellipse_vertical
+
+
+@pytest.fixture()
+def filled_ellipse_horizontal() -> npt.NDArray:
+    """Ellipse (horizontal)."""
+    filled_ellipse_horizontal = np.zeros((7, 9), dtype=np.uint8)
+    rr, cc = draw.ellipse(3, 4, 3, 4)
+    filled_ellipse_horizontal[rr, cc] = 1
+    return filled_ellipse_horizontal
+
+
+@pytest.fixture()
+def filled_ellipse_angled() -> npt.NDArray:
+    """Ellipse (angled)."""
+    filled_ellipse_angled = np.zeros((9, 11), dtype=np.uint8)
+    rr, cc = draw.ellipse(4, 5, 3, 5, rotation=np.deg2rad(30))
+    filled_ellipse_angled[rr, cc] = 1
+    return filled_ellipse_angled