PMEAL / porespy

A set of tools for characterizing and analying 3D images of porous materials
https://porespy.org
MIT License
309 stars 100 forks source link

Method for increasing max_iters for saddle point removal with SNOW #276

Closed 1minus1 closed 3 years ago

1minus1 commented 4 years ago

I am working with some large images (1000^3), and I am getting warnings that:

Maximum number of iterations reached, consider running again with a larger value of max_iters

This is happening in SNOW's call to trim_saddle_points, which has max_iters set to 200. There is no interface for increasing max_iters to complete saddle point removal, and the output network ends up giving errors during a subsequent invasion percolation, where networks made from smaller images work just fine (so it might be somehow malformed?).

I would prefer for the algorithm to take as many iterations as it needs, or at least have an option to tell it to do so.

jgostick commented 4 years ago

It should not need that many iterations...it's probably a sign that something is not quite right with the image, or that the function is not handling something correctly. I have been meaning to dig into that function anyway, to handle this issue better and to also work faster (i.e. with numba).

However, if you want to experiment with the parameter, you can call the individual steps of the snow procedure directly, as outlined in this example

1minus1 commented 4 years ago

It happens with fairly standard structures. Here's an example script. I tested this in the feature_skip_snow_steps branch, but it should give the same result in dev if you remove the last two arguments to SNOW. It complains about not enough iters for the 1000^3 structure, but not 100^3. For what it's worth, I think these images are fine. When network extraction completes, I have been able to replicate mercury porosimetry results from porespy's pixel method with openpnm's network method for the same source image.

import porespy as ps

dims = [[100, 100, 100], [1000, 1000, 1000]]
fiber_ns = [41, 4362]

for dim, fiber_n in zip(dims, fiber_ns):
    im = ps.generators.cylinders(dim, 5, fiber_n, 0, 90)
    print(ps.metrics.porosity(im))
    snownw = ps.networks.snow(im,
                              voxel_size=1,
                              boundary_faces=['left', 'right', 'front', 'back', 'top', 'bottom'],
                              marching_cubes_area=False,
                              skip_trim_saddle=False,
                              skip_trim_nearby=False)

Output for the 1000^3:

0.810147299
------------------------------------------------------------
Beginning SNOW Algorithm
Converting supplied image (im) to boolean
Peforming Distance Transform
Applying Gaussian blur with sigma = 0.4
Initial number of peaks:  46089
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Maximum number of iterations reached, consider running again with a larger value of max_iters
Peaks after trimming saddle points:  34156
Peaks after trimming nearby peaks:  31748
1minus1 commented 4 years ago

I created a branch where max_iters was set to 50000, and still was getting this result.

jgostick commented 4 years ago

Ya, I suspect there is a condition other than the intended stop conditions that is screwing things up. Can you perhaps try to pad the marker image with 0's on all sides? I'm wondering if the makers are touching the edge of the image and that is the problem? I probably can't dig into this for another week or two though, but once the summer semester is over, I can get on it.

1minus1 commented 4 years ago

There is the same effect when the source image is padded with zeros.

import porespy as ps
import numpy as np
import imageio

dims = [[100, 100, 100], [1000, 1000, 1000]]
fiber_ns = [41, 4362]

for dim, fiber_n in zip(dims, fiber_ns):
    im = ps.generators.cylinders(dim, 5, fiber_n, 0, 90)
    # imageio.mimsave(str(dim[0])+'nopad.tif', np.array(im, dtype=np.int32))
    im = np.pad(im, ((1, 1), (1, 1), (1, 1)), 'constant')
    # imageio.mimsave(str(dim[0])+'pad.tif', np.array(im, dtype=np.int32))
    print(ps.metrics.porosity(im))
    snownw = ps.networks.snow(im,
                              voxel_size=1,
                              boundary_faces=['left', 'right', 'front', 'back', 'top', 'bottom'],
                              marching_cubes_area=False,
                              skip_trim_saddle=False,
                              skip_trim_nearby=False)
jgostick commented 3 years ago

this is not necessary since you can just call each function separately. i will still look into the function to see if i can figure out why it happens in the first place.