Possible changes/addition to the label surface function

[DMaxJ]

While working with extracted networks, it becomes essential to find the surface pores and add labels to them. The present function inside topotools allow that functionality, but it has some limitations. 1) The function

label_faces(network, tol=0.0, label='surface')

takes the tolerance parameter as a fraction of the xspan, yspan or zspan. Thus, if the dimensions of the domain are not same in all directions it will label different fractions of pores (e.g. for a network having dimensions 800x400x150, the labeled portion will be different in each direction ) 2) Inside this function it calls another function (if required)

find_surface_pores(network, label=label)

which is not very easy to control. When used with default parameters it ends up producing much deviation when compared to networks extracted using the maximal ball algorithm.

Now, I have two suggestions: Suggestion 1

def custom_label_faces(network, tol=0.0, label='surface'):
    label = label.split('.', 1)[-1]
    if 'pore.'+label not in network.labels():
        find_surface_pores(network, label=label)
    Psurf = network['pore.'+label]
    crds = network['pore.coords']
    xmin, xmax = sp.amin(crds[:, 0]), sp.amax(crds[:, 0])
    xspan = xmax - xmin
    ymin, ymax = sp.amin(crds[:, 1]), sp.amax(crds[:, 1])
    yspan = ymax - ymin
    zmin, zmax = sp.amin(crds[:, 2]), sp.amax(crds[:, 2])
    zspan = zmax - zmin
    dims = dimensionality(network)
    if dims[0]:
        network['pore.front'] = (crds[:, 0] <= (xmin + tol)) * Psurf
        network['pore.back'] = (crds[:, 0] >= (xmax - tol)) * Psurf
    if dims[1]:
        network['pore.left'] = (crds[:, 1] <= (ymin + tol)) * Psurf
        network['pore.right'] = (crds[:, 1] >= (ymax - tol)) * Psurf
    if dims[2]:
        network['pore.top'] = (crds[:, 2] >= (zmax - tol)) * Psurf
        network['pore.bottom'] = (crds[:, 2] <= (zmin + tol)) * Psurf

Here the tolerance is used as a fixed values. This can be chosen depending on the average pore diameter of the network. However, this method still suffers if the surface pores are not found properly.

Suggestion 2

def custom_label_faces2(network, tol=0.0):
    crds = network['pore.coords']
    xmin, xmax = sp.amin(crds[:, 0]), sp.amax(crds[:, 0])
    xspan = xmax - xmin
    ymin, ymax = sp.amin(crds[:, 1]), sp.amax(crds[:, 1])
    yspan = ymax - ymin
    zmin, zmax = sp.amin(crds[:, 2]), sp.amax(crds[:, 2])
    zspan = zmax - zmin
    dims = dimensionality(network)
    if dims[0]:
        network['pore.front'] = crds[:, 0] <= (xmin + tol)
        network['pore.back'] = crds[:, 0] >= (xmax - tol)
    if dims[1]:
        network['pore.left'] = crds[:, 1] <= (ymin + tol)
        network['pore.right'] = crds[:, 1] >= (ymax - tol)
    if dims[2]:
        network['pore.top'] = crds[:, 2] >= (zmax - tol)
        network['pore.bottom'] = crds[:, 2] <= (zmin + tol)

This method takes out the find_surface_pores(network, label=label) function and allows to label pores within a certain fixed range from the extremities of the network. When a tolerance value with 2 to 2.5 times the average or maximum pore radius (the factor is variable) is used, this produces surface pores similar to maximal ball algorithm.

I know, this is not exactly a proper method of finding surface pores, but it somehow works. If such a functionality can be added to openpnm, it will be a great help

[jgostick]

Thanks for the thoughtful suggestion. I really appreciate that you put some sample code down.

The issue is that there is NO universal/rigorous way to find surface pores. The 'proximity to extremity' approach works on networks with a narrow pore size distribution, but on wide psds would find many small pores on surfaces, or if the threshold were too small would find no big pores. The way the functions currently work in openpnm is to use a Delaunay tessellation of the pores with some fictitious points outside the domain. This creates a sort of 'line of sight' connection, which lets us identify surface pores as those who formed a connection with the fictitious points. You are able to add more points (called 'markers') when calling the function, to help find more surface pores. This is the only feasible, general and reasonably rigorous way to do it that I've found. It could be done with image processing as well, but this would make it infeasible due to memory requirements.

[DMaxJ]

Thanks. I understood your point. This was just my kind of quick fix which worked at least for the types of materials that I deal with.

[jgostick]

I'm going to close this. Finding pores based on thresholding the coordinates is a fine approach in many cases, but is not really general enough.