New submodule for feret calculations

[ns-rse]

New submodule measures.feret with functions for calculation of min and max feret and the associated co-ordinates.

Utilises code from Skan developer see gist and as suggested it adds tests. In doing so I found sorting points prior to calculation of upper and lower convex hulls missed some points.

I'm also not sure about the rotating_calipers() function as it doesn't seem to return all pairs formed across the points from the upper and lower convex hulls and so have included but not used the all_pairs() function which does, although if used in place this doesn't return the minimum feret correctly, rather just the smallest distance.

Currently some of the tests for the curved_line fail too.

The intention is to use this without instantiating the GrainStats class to be used in profiles (#748) and could serve as a concise stand-alone set of functions outside of GrainStats class which currently has static methods for the calculations (#750).

Benchmarking

In light of #750 and re-implementing feret calculations as a new sub-module I wanted to see how this compares in terms of performance to those in GrainStats() class so have run some very basic benchmarking. Note this is not the full pipeline of taking labelled images and finding the outlines/edge points which are required as inputs to the calculations, its purely on the calculation of min-max feret from the edge points.

import timeit

import numpy as np
from skimage import draw

from topostats.measure import feret
from topostats.grainstats import GrainStats

holo_circle = np.zeros((14, 14), dtype=np.uint8)
rr, cc = draw.circle_perimeter(6, 6, 5)
holo_circle[rr, cc] = 1

holo_circle_edge_points = np.argwhere(holo_circle == 1)


%timeit feret.min_max_feret(holo_circle_edge_points)

83 µs ± 686 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

%timeit GrainStats.get_max_min_ferets(holo_circle_edge_points)

1.06 ms ± 10.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

So this new implementation appears faster.

[codecov]

Codecov Report

Attention: Patch coverage is 99.37107% with 1 lines in your changes are missing coverage. Please review.

Project coverage is 85.51%. Comparing base (4ae1b70) to head (4e3db2e).
Report is 49 commits behind head on main.

Files	Patch %	Lines
topostats/measure/feret.py	99.37%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #755      +/-   ##
==========================================
+ Coverage   84.73%   85.51%   +0.77%     
==========================================
  Files          21       22       +1     
  Lines        3196     3368     +172     
==========================================
+ Hits         2708     2880     +172     
  Misses        488      488              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ns-rse]

Please hold back on this, I've been working on simplifying the tests and think there is a problem with the Rotating Caliper algorithm so will be checking that next week.

[ns-rse]

The inverse gradient method of calculating the triangle height formed when determining the minimum feret results in ZeroDivisionError when the base of the triangle is horizontal so switching methods (thanks for the comprehensive suggestions of both @SylviaWhittle).

In order to check everything I wrote a plotting function (see previous commit) and whilst many of the minimum feret distances are correct there are still instances where the triangle height is calculated outside of the convex hull which is baffling. For the most part not an issue but in the holo_ellipse_angled this is an issue.

Figures below show the problem.

The same Graham Scan algorithm is used for the rotating caliper pairs as in GrainStats so unsure if this is an issue there and can't check as co-ordinates aren't currently returned (possibly bounding box may be informative as I think I saw that was calculated for the feret distances, but inclined to solve this in current code base).

Slightly confusing plots but they show the lower and upper hulls in red/green (apologies if anyone is colourblind) with the rotating calipers in yellow. The Triangle heights are in blue and it is these that are wrong for some, in most cases this isn't an issue as the min feret (in magenta as is maximum feret) are elsewhere)....

BUT that won't always be the case as the following shows, the minimum feret is determined to be a line outside of the convex hulls and taking the trace here would give an incorrect set of heights.

These tests may seem/feel like over-simplifications but if the code isn't correct for these then we have no guarantee it will be correct for more complex shapes.

I feel like I'm close to solving this and its important to me that its correct as I would hate to have results published that contained a corrigendum.

Once I've fixed these I will reinstate the tests for filled objects and images with multiple masks (the later isn't essential but it may be useful to users using the functions interactively and will "just work (TM)" once everything else is correct as its just a wrapper).

[ns-rse]

Started investigating this and have worked out where the problem arises and have some thoughts on how to resolve it.

Taking this small ellipse...

The rotating_calipers() function finds the first caliper between [5,2] and [1,2] that is the horizontal line across the middle of the ellipse.

The next point on the adjacent hull is [4,1] thus the triangle is formed by [[5,2],[4,1],[1,2] which is a scalene triangle

Extend the base and we have...

The height of the triangle, which is from the apex to the perpendicular intersection with the base is outside of the original triangle...

But this doesn't matter its not the minimum feret

In this example its not but as noted above this tilted ellipse has a similar error and the minimum feret (short pink line) is outside of the shape.

Checking this we have...

Caliper/Triangle 1

Problematic A scalene triangle where the co-ordinates of the height (which we are using as the feret distance) lie outside of the convex hull. This isn't so much of a problem as caliper/triangle 4 as the height is > than the minimum feret so it would never be chosen but we should still be excluding this.

Caliper/Triangle 2

Problematic A scalene triangle where the co-ordinates of the height (which we are using as the feret distance) lie outside of the convex hull. This isn't so much of a problem as caliper/triangle 4 as the height is > than the minimum feret so it would never be chosen but we should still be excluding this.

Caliper/Triangle 3

NB This is a right-angled triangle so the height is within the bounds of the triangle and what we want.

Caliper/Triangle 4

Offender this is the triangle who's height is equal to the minimum feret but lies outside of the convex hull.

Caliper/Triangle 5

Caliper/Triangle 6

NB This is a right-angled triangle so the height is within the bounds of the triangle and what we want, you can't see it because its the black side of the triangle and I only wrote a short hacky function to draw these plots.

Ok, this is doing what we have asked which is find the height of the triangle formed by the calliper pairs and the next point on the convex hull, but its not what we want to know and its actually wrong if the triangle formed is scalene.

Why are the minimum feret co-ordinates outside of the triangle and not the height of caliper/triangle 2?

The rotating_calipers() function returns a generator, for each pass round the rotating calipers it yields the min_feret, the caliper co-ordinates and the co-ordinates of the min_feret (the apex to the point perpendicular to the base of the triangle, which as we have seen can be outside of scalene triangles).

The min_ferets and min_feret_coords are then zip() up and passed to the min() function which returns the minimum. It does so based on the first element, if any of the first elements (i.e. min_feret) match then they are sorted by the min_feret_coords.

This list looks like this...

[
    (5.0, [[6.0, 2.0], [1, 2]]),
    (5.0, [[1.0, 5.0], [6, 5]]),
    (2.82842712474619, [[3.0000000000000004, 2.0000000000000004], [1, 4]]),
    (2.82842712474619, [[0.0, 3.0], [2, 1]]),      
    (7.071067811865475, [[6.661338147750939e-16, 2.999999999999999], [5, 8]]), 
    (7.071067811865475, [[6.0, 7.0], [1, 2]])
]

The third and fourth elements have the same minimum feret distance (and you can see that the min_feret in the coloured plots in pink is indeed the same distance as the caliper/triangle 3 height) but because the co-ordinates of the feret line in element 4 of this list are smaller it is selected as the minimum feret.

What to do?

The Graham Scan algorithm is the basis for this implementation (as it is in the current grainstats.GratinStats.get_max_min_feret() function, note the calls to GrainStats.get_triangle_height() which is taking the caliper and the next point on the convex hull, so the problem already exists in the code base and it is possible the minimum feret distances reported are wrong).

If algorithms don't work on small examples we can't expect them to generalise well and I would never have come across this issue if I hadn't created the various small simple shapes and tested them!

I think a sensible approach here might be to have a check made before finding the minimum feret that checks if the co-ordinates of the triangle height are on the convex hull. If either of these points are outside of the shape we can drop that feret distance/co-ordinates from the list prior to finding the minimum which means we should only ever get feret distances that are within the convex hull and we therefore get not just the correct distance but the correct co-ordinates through which we can take a profile (which is what has stimulated this work in the first place!).

I know you're incredibly busy @SylviaWhittle but as you tussled with this problem in the past your thoughts would be appreciated as would any thoughts from @MaxGamill-Sheffield or @llwiggins .

Possible solutions

Various solutions for finding whether points are within the convex hull are suggested in this thread along with timings, some parallel, some using numba for Just in Time compilation to speed things up, others using existing packages such as matplotlib which we already have as a package dependency or shapely which we don't currently use but may have utility beyond the current problem (constrained to 2-D images though so maybe not as lots comes down to the height erosion maybe different from that which is currently used and may improve on some of the current methods 🤷 ).

[ns-rse]

Notes on solutions for excluding triangle heights from scalene triangles that are outside of the convex hull

Matplotlib Path

Gave this method a whirl but as noted in the documentation.

The current algorithm has some limitations:

• The result is undefined for points exactly at the boundary (i.e. at the path shifted by radius/2).
• The result is undefined if there is no enclosed area, i.e. all vertices are on a straight line.
• If bounding lines start to cross each other due to radius shift, the result is not guaranteed to be correct.

Had a play with the value of radius but without success which excludes this method as many of the points we will be assessing will be on the boundary.

Did like the fact you could assess a vector of co-ordinates rather than having to loop over them though.

is_inside_sm()

This method looks promising as it works for points on the edge of the polygon (where Matplotlib failed) but being wary of copy-pasta-ing code would want to have robust tests in place (the above link to StackOverflow suggests it is fast and scales well though).

Shapely

Shapely looks like a very useful package. There are a number of functions for determining whether a point is within a polygon.

I've gone with this package and had to make another check because if the line was identical to the edge which may arise then Shapely contains() returns False and the line (a possible minimum feret distance) is marked as being outside which is not what we want.

As an example the angled ellipse previously identified the minimum feret as being outside of the the hull (the shorter pink line)...

That line is now excluded and shows as blue in the following as it is one of the triangle heights, but the minimum feret (which is the same distance) is now correctly labelled as the short pink line.

Aside

Also the shapely.orientated_envelope() may be a faster solution to the overall problem returning the bounding rectangle of an object with the minimum area (i.e. defined by the min/max feret distances). Think I saw similar calculations in GrainStats class perhaps?

[ns-rse]

This is ready for review.

Once approved and merged I shall proceed with extracting the profiles along the min/max feret lines (the main feature from #748) and addressing #798 as I think the code implemented in GrainStats suffered from the above mentioned problem that scalene triangles causes when calculating minimum feret distances.

Test coverage is low on this PR as I've not tested the topostats.measure.feret.plot_feret() function I wrote to aid investigation. Can address this but its really for debugging and would bloat the repository adding figures for each scenario (although could easily keep them as small as possible in order to minimise that side effect).

[ns-rse]

Update - Decided to try the ray tracer algorithm and as we want to include points on the convex hull have used >= and <= but have a slight problem with tests failing for points on the left and top edge of shapes as the odd value gets flipped twice. I think since these are convex hulls rather than concave we might be able to to use the number of times crosses are made being >= 1 as indicative of a point being on or within but will check next week.

[ns-rse]

After discussion with @SylviaWhittle and @llwiggins it was decided that because of the issue of minimum feret coordinates sometimes falling outside of the convex hull that we should not use it as a consistent line through the grain/image.

Instead we will use the maximum feret and a line perpendicular to the mid-point of the maximum feret.

As such all attempts to work out if minimum feret coordinates/lines were inside or on the edge of a polygon have now been removed and the function returns a dictionary of...

• min_feret
• min_feret_coords
• max_feret
• max_feret_coords

The dictionary should be conducive to being saved in the HDF5 output (perhaps a separate issue to do so).

@SylviaWhittle : Would it be useful to have an additional element in the returned dictionary for the triangle co-ordinates from which the minimum feret (triangle height) co-ordinates are derived?

[SylviaWhittle]

I ran this locally using some test shapes that I generated and it works perfectly. Returning the coordinates of the min and max ferets is really useful, thank you.

The code reads beautifully and is very comprehendible. I also like the several options that you've added for running it.

Thank you for your patience and sorry for not reviewing it last week as requested.