Rounding issue in local_thickness

[gourrier]

Hi,
I recently bumped into porespy which has exactly the bunch of features I was looking for.

Q0: (added at the end): why use scipy for array manipulations where numpy would be the logical package to use ? Sure this adds another import, requirements etc, but it would definitely be cleaner.

This includes the local_thickness function (I don't understand why on earth you want to remove it from the package ???).

One thing I find odd is the following line:
sizes = sp.unique(sp.around(dt, decimals=0))

---

Q1: This basically results in even numbers of pore radii (pixels). What is the point of doing this ?
I.e. the output from the former gives radii of 0, 1, 2, 3... => diameters of 0, 2, 4, 6... but diameters of 3, 5, 7... pixels should also be allowed ? What am I missing here ?

If the idea is to have integer numbers of pixels, then I think that the line should rather be:
sizes = sp.unique(np.round(dt*2.)/2.)
which gives half integers as radii an thus pixel numbers as diameters (Note, np is import numpy as np).

---

Q2: But then, this somehow results in a strange value distributions as, e.g. 2*2**.5, 3, 10**.5 (3 consecutive euclidian distances, i.e. (2,2), (3,0), (3,1)) would all result to values of 3.

So why not simply round to 2 digits to have more realistic distances and let the user decide what he wants to do with these values ? I.e.:
sizes = sp.unique(np.round(dt,decimals=2))
BTW: this gives a result which is closer to the local_thickness calculation of ImageJ (FIJI) and BoneJ (albeit with a factor 2 less since they return the diamter as a value of "thickness").

---

SUGGESTION: A simple effortless workaround would be to add a keyword: exact_pixels = False in the function and then replace the above line by:
if exact_pixels:
sizes = sp.unique(np.round(dt*2.)/2.)
else:
sizes = sp.unique(np.round(dt,decimals=2))

At least this is exactly what I will as a patch to this module if it doesn't change...

Cheers,
Aurélien

[jgostick]

Hi Aurelien
Thanks for your comments! I'll try to answer your questions:
(0) Scipy is a super set of numpy, so all array ops using scipy really just call numpy. Scipy has the very important ndimage module, which we use extensively. Our plan is to remove the code within the local thickness function, and have it call porosimetry behind the scenes. This way we only have to maintain one version of the code. Local thickness is just a porosimetry simulations with access limitations turned off. RE the sizes thing, we do this to reduce the number of points we perform the simulation on to save time. You are able to override this by sending in sizes yourself...but be prepared to wait!

Regarding Q1 and Q2, they sound like bugs, so I'll need more time to investigate. I'll get back to you soon.

[gourrier]

Thanks. I know about the scipy/numpy connectivity and I also use ndimage quite a bit, but it's probably the first time I see such heavy use of scipy for handling arrays this way (I generally see both imports and use of numpy for matrices). But this is just a lack of experience on my side... so don't worry about this part...

On the size idea, I understand better why you went for this option, now. The thing is that I'm currently running a script to measure porosity at the resolution limit, i.e. 1-3 pixels roughly. I know, it's really not optimal to run a distance map calcuation, but our voxels are (14 nm)^3 so this is actually quite good in terms of imaging... So in this specific and probably odd context, I need precision on the size distribution.

So all in all, maybe adding a simple keyword specifying how precise you want the calculation to be (and computing time you can afford) would fit both your everyday needs and my specific one.

[jgostick]

I just did a quick check with the following:

import scipy.ndimage as spim
import porespy as ps
im = ps.generators.blobs(shape=[500, 500])
lt = ps.filters.local_thickness(im=im, sizes=sp.unique(spim.distance_transform_edt(im)))

This works as you are looking for by specifying the exact sizes to use, BUT is does over 1000 loops! There are just too many sizes, so we need some way to do it by bins. The integer value of radius seems to make sense. Note that all the convolution take a structuring element of given radius, so doing diameter doesn't really work. The radius needs to be a whole number of voxels.

Anyway, can you try sending your desired sizes into the sizes argument and let us know if it works for you, and if this is a reasonable approach, as opposed to adding a flag to the function as you suggested? Using sizes essentially gives the user infinite control.

[jgostick]

Since this function works as we think it should, and @gourrier hasn't replied in a while, I'll assume the problems are all sorted and close this.