feature(measure): Adds the height_profiles sub-module

[ns-rse]

Further work on height profiles this PR adds the topostats.measure.height_profiles sub-module which for a given pair of feret coordinates...

1. Determines the orientation of of the feret co-ordinates relative to horizontal.
2. Rotates the image so that the feret is horizontal.
3. Recalculates the co-ordinates of the feret after rotation.
4. Extracts the height profile along the feret.

Includes a function in topostats.plotting.plot_height_profiles() which produces a line-plot of multiple height profiles.

Test shapes defined tests/measure/test_feret.py moved to fixtures in tests/measure/conftest.py to avoid duplication and tests for all functions are included. Because these fixtures are then used in @pytest.mark.parameterize() it has necessitated using the request fixture and its .getfixturevalue() method.

More on this can be read in this blog
post.

ToDo...

Still a fair few steps to integrate this into the processing.

• Add configuration options to topostats/default_config.yaml of whether to calculate height_profiles.
• Update GrainStats() to calculate the height profile for the image being processed if required based on configuration.
• Return the height_profile (1-D Numpy array) as part of the tuple GrainStats() returns.
• Collect height_profile across grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in
  plot_height_profiles()).
• Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc.
• Currently only configured to give height profiles across maximum feret, would need additional work and consideration for profiles perpendicular to the mid-point of the maximum feret (but code is structured such that the extraction of this would be relatively easy, but would require additional work on plotting and returning data).

Related : #748 #755

[ns-rse]

@llwiggins has identified a scenario where behaviour is not as expected.

import numpy as np
from topostats.measure import height_profiles

circle = np.array([
    [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]
])

binary_mask = circle

height_profiles.extract_feret_profiles(circle, binary_mask)

This fails the following test...

import pytest
import numpy.typing as npt

@pytest.mark.parametrize(
    ("img", "target"),
    [
        pytest.param(
            {
                "img": np.asarray(
                    [
                        [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],
                    ]
                ),
                "skeleton": np.asarray(
                    [
                        [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],
                    ]
                ),
            },
            np.asarray([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]),
        ),
    ],
)
def test_extract_feret_profiles2(img: npt.NDArray, target: npt.NDArray, request) -> None:
    """Test extract_feret_profiles()."""
    # _img = request.getfixturevalue(img)
    # Convert boolean to 0/1 and extract profiles
    feret_profiles = height_profiles.extract_feret_profiles(img=img["img"], skeleton=np.where(img["skeleton"], 1, 0))
    # This fails...
    # feret_profiles = height_profiles.extract_feret_profiles(img=_img["img"], skeleton=_img["skeleton"])
    np.testing.assert_array_almost_equal(feret_profiles, target, decimal=21)

E           AssertionError: 
E           Arrays are not almost equal to 21 decimals
E           
E           Mismatched elements: 10 / 11 (90.9%)
E           Max absolute difference: 1
E           Max relative difference: 1.
E            x: array([0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0])
E            y: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

Why does this fail?

Whilst the image doesn't need rotating this fails because the maximum feret co-ordinates after rotating (which wasn't in effect performed) are a vertical line rather than the expected horizontal line...

rotated_feret_stats={'max_feret': 10.0, 
                     'min_feret': 8.48528137423857, 
                     'max_feret_coords': array([[10,  5],
                                                [ 0,  5]]),
                     'min_feret_coords': array([[1., 3.],
                                                [7., 9.]])}

And the profile() function is very basic, taking an image and the row from which the profile is to be returned...

def profile(img: npt.NDArray, row: int) -> npt.NDArray:
    """
    ...
    """
    try:
        return img[row]
    except IndexError as e:
        raise IndexError("The slice (row) is outside the image boundary.") from e

And so because the row being used is rotated_feret_stats["max_feret_coords"][1, 0] a slice is taken across the first row which is why we see the returned value.

[ns-rse]

@llwiggins I've revised how the profiles are obtained and ditched the rotation method instead using scipy.interpolate.RegularGridInterpolator.

As a consequence the topostats.measure.heights_profile sub-module is considerably smaller and consists of a single function interpolate_height_profiles() which performs interpolation of heights directly between two points by first getting evenly spaced points on the x-axis between the maximum feret co-ordinates and then interpolating the corresponding y-coordinates. These x and y co-ordinates then have the height interpolated using RegularGridInterpolator.

Arguments can be passed via **kwargs to RegularGridInterpolator() and the effect of modifying the method is tested. I've included your example of a basic circle which caused problems previously and those tests pass now (whenever a scenario that causes failure is encountered its a good candidate for a test so thank you).

I've found for some reason that interpolation (using the method = linear) on M$-Win gives different results so have marked those tests to be conditionally skipped if running on a Windows system.

[SylviaWhittle]

Been reading the PR and looks good to me. Much prefer the interpolation method, so thank you. Not given it a thorough enough read & tinker with to approve yet, but I like the method and implementation! I'll test it out soon (once my report is done at the latest). Apologies for not being active much on TopoStats recently.

[llwiggins]

Thank you @ns-rse for revising the original method for obtaining profiles - I really like the use of scipy's RegularGridInterpolator() here. I've tested the code out with various dummy examples, plotting graphs throughout the workflow to assess different stages, and all seems to be working as expected. Happy for this work to be merged. 😄