Rounding issue in local_thickness

[pmeal]

Rounding issue in local_thickness

[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.