Crofton perimeter and Euler number

[yg42]

Description

This PR consists in 2 changes in regionprops functions

1. A new Euler number function is proposed, that uses convolution to get configurations codes, followed by a scalar product. This algorithm is in O(m), with m the total number of pixels in the image. It works for 2D and 3D images. Notice that I propose this method because it seemed that there was an error in the previous version (I compared with the matlab version). I got the coefficients from [1] for 2D images and [2] for 3D images.
   [1] S. Rivollier. Analyse d’image geometrique et morphometrique par diagrammes de forme et voisinages adaptatifs generaux. PhD thesis, 2010. Ecole Nationale Superieure des Mines de Saint-Etienne. https://tel.archives-ouvertes.fr/tel-00560838
   [2] Ohser J., Nagel W., Schladitz K. (2002) The Euler Number of Discretized Sets - On the Choice of Adjacency in Homogeneous Lattices. In: Mecke K., Stoyan D. (eds) Morphology of Condensed Matter. Lecture Notes in Physics, vol 600. Springer, Berlin, Heidelberg

2. A new perimeter method, based on the Crofton formula [3], is proposed. The same algorithm as for Euler number is used. It works only for 2D images.
   I added an example to illustrate the approximation of the 2 different perimeter methods with 4 and 8 connectivity. The perimeter of a rotated square is represented as a function of the rotation angle.
   [3] https://en.wikipedia.org/wiki/Crofton_formula

Note: There is an interesting link with #3812 and #3797, because minkowski functionals (volume, surface are, integral of mean curvature and Euler number in 3D) can be evaluated by the same method proposed in the previous references (see [4]). You can also refer to
Geometric measure in 2D and 3D for matlab, by D. Legland where a similar method is implemented using a marching cube approach.

[4] Ohser, J., & Schladitz, K. (2009). 3D images of materials structures: processing and analysis. John Wiley & Sons.

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.
• Consider backporting the PR with @meeseeksdev backport to v0.14.x

[emmanuelle]

thank you @yg42 ! I'll review your PR in the next couple of days, I'm very excited about having these features for 3D images. Maybe @alexdesiqueira is also an appropriate reviewer.

[emmanuelle]

So I started playing a bit adapting your example and it's possible that there is a bug somewhere because the result below looks fishy to me:

from skimage.measure import perimeter
from skimage.measure import crofton_perimeter
import numpy as np
from skimage import draw

img = np.zeros((100, 100))
rr, cc = draw.circle_perimeter(50, 50, 20)
img[rr, cc] = 1

print(perimeter(img, 4))
print(perimeter(img, 8))

print(crofton_perimeter(img, 2))
print(crofton_perimeter(img, 4))
----
131.88225099390854
131.88225099390854
251.32741228718345
223.40713078307576

I would not expect the Crofton perimeter to be so different from perimeter... Did I miss something? I'll keep on investigating.

[yg42]

@emmanuelle
My implementation of Crofton perimeter supposes that you have filled objects. If there is a hole, the inner contour is also added in the perimeter value.

Perimeter function only counts the outer contour, is this an expected behavior? I would say no. If I test on a disk I get the same result. My proposition of Crofton perimeter gives approximately the same perimeter for a disk.

from skimage.measure import perimeter
from skimage.measure import crofton_perimeter
import numpy as np
from skimage import draw

img = np.zeros((100, 100))
rr, cc = draw.circle_perimeter(50, 50, 20)
img[rr, cc] = 1
print(perimeter(img, 4))

rr,cc = draw.circle(50, 50, 20);
img[rr, cc] = 1
print(perimeter(img, 4))
crofton_perimeter(img, 4)
----
131.88225099390854
131.88225099390854
127.71373126734748

[emmanuelle]

Oh yes my example was actually very convoluted :-), sorry about that. I was trying to find an example for which crofton perimeter would get more accurate results than the perimeter function with 4 neighbors, because in the original gallery example you provided, they look really close.

[emmanuelle]

Would it be possible to explain more the pros and cons of the two algorithms?

[emmanuelle]

all values strictly greater than zero

[emmanuelle]

is it possible to find a more self-explaining name than F?

[yg42]

I changed it to config, as this describes a configuration neighbourhood.

[yg42]

Oh yes my example was actually very convoluted :-), sorry about that. I was trying to find an example for which crofton perimeter would get more accurate results than the perimeter function with 4 neighbors, because in the original gallery example you provided, they look really close.

There is no big difference, which is good!

As you know, the perimeter evaluation is really tricky and there is no absolute method to measure it. You can only have approximations. Crofton perimeter, as well as previously present perimeter function in skimage are one of them. OpenCV has a different approach: it approximates the contour with polygons and evaluates the arc length. In fact, the question of evaluating the perimeter is linked to finding the original boundary of the discretized object. Maybe I will work on this last approach and try to include it in skimage.

In 3D, this is the same problem as surface area evaluation, I know you mentioned it in feature requests.

This is not easy to evaluate the efficiency of a perimeter function, and I wouldnt dare saying that one method is better than the other.

[emmanuelle]

2D or 3D

[yg42]

ok

[emmanuelle]

I would remove "see variable config" since we don't expect users reading the docstring to read the code as well. However you can move this to a comment within the function.

[yg42]

ok

[emmanuelle]

@yg42 regarding perimeter evaluation based on polygons, we already have the building blocks for this with marching cubes in 3d and measure.find_contours in 2D, but indeed we're not using it in regionprops.

[emmanuelle]

Out of curiosity, why did you want to have Crofton perimeter in scikit-image? The Euler number in 3D was clearly missing (all the more if the 2D version is buggy), but I have more mixed feelings about adding another perimeter function.

[emmanuelle]

Thank you very much @yg42. I left a couple of minor comments, besides we also need

• to decide whether the Crofton perimeter is needed in scikit-image or not (@skimage/core ?)
• if yes, to flesh out a bit the introduction paragraph of the perimeter gallery example to explain the different between perimeter functions
• and to add a gallery example with the Euler characteristic.

I know it's a lot of work, please bear with us!

[yg42]

@emmanuelle
Thank you for your comments, I really appreciate working on this proposition. I will work on this until the PR is completed.

About the Crofton perimeter interest, I wanted to make a presentation to my colleagues who work in the field of chemical engineering. As most of the users of skimage (I suppose), they are not aware that a perimeter in discrete objects is only an approximation, and I wanted to illustrate this by simple examples (I added the perimeter evaluation of a rotated square in the examples). For this purpose, I coded the Crofton perimeter. As this code is really similar to the Euler characteristic evaluation, I detected some problems on it and decided to submit these two methods.

In my opinion, Crofton perimeter has to be in the measure submodule. Does it have to be in regionprops? It can create confusion for general users, but it is an interesting property, mathematically elegant, and practically quite precise for high sampling and general objects.

I will soon add explanations and discussion in the gallery example, even if the Crofton perimeter is not included, it can be interesting to warn people about this approximation.

I will propose an example illustrating the Euler number, but I have to think of it.
yann

[emmanuelle]

Thank you @yg42 ! One thing I forgot: since you're adding new features and the behavior of euler_number changed, this should be mentioned in the release notes doc/release/release_dev.rst

[yg42]

ok, I mentioned it as a bug and as a new feature (the function is now a method and has neighbourhood parameters).
I also modified the description in regionprops in order to describe the behavior for 3D images as well as for label images.

I will wait for a decision about the crofton perimeter.

[yg42]

I did my best to rebase on trunk. Just tell me if I still have things to do. :)

[yg42]

@emmanuelle
What should I do now ?

[alexdesiqueira]

@jni @hmaarrfk I get afraid when I see lots of old PRs within one, like here. Should I?
Other than that, LGTM.

[yg42]

@alexdesiqueira @jni @emmanuelle
Is there something to do ?
yann

[jni]

@yg42 apologies, I have had an extremely busy past few weeks. I will try to review this properly next week, if someone doesn't beat me to it!

[soupault]

I think this is a very exciting PR, but I need a bit more time to read through the theoretical aspects.
The rebase did go wrong, it seems, let's see what we can do.

[yg42]

Hi all,
Its been a long time since this discussion started, I am not sure what needs to be done. I am kind of busy these days, but if needed, I could relaunch this PR.
all the best
yann