Introduction

Geometrical and topological representations of the internal structures of samples
used in Bio-CAD (bones, biomaterials, rocks and other material samples) are
necessary to evaluate their physical properties. Such samples have two disjoint
spaces, the pore space and solid space; thus, they can be represented with binary
models. The pore space, in turn, is made up of a collection of pores that can be
intuitively defined as local openings interconnected by narrow apertures called
throats that limit the access to a larger pore (^{Garboczi, Bentz & Martys, 1999}) (Figure 1).

Usually, porous datasets are considered adequate for analysis if they have their
holes homogeneously distributed along the sample, and the pore size is notably
smaller than the dataset size (^{Garboczi et
al., 1999}); that is, they contain complete holes and the
object resolution is high enough to clearly define their shape.

The study of porous materials’ properties is of great utility in several disciplines.
In medicine, they are used to evaluate the degree of osteoporosis and the adequacy
of synthetic biomaterial implants for bone regeneration, among other applications.
Bone regeneration occurs in the cavities of the implants where blood can flow.
Biomaterial implants can be designed as tissue scaffolds, that is, extracellular
matrices onto which cells can attach and then grow and form new tissues (^{Giannitelli, Accoto, Trombetta & Rainer,
2014}; ^{León-Mancilla, Araiza-Téllez,
Flores-Flores & Piña-Barba, 2016}). In geology, the porous morphology
and pore sizes are related to the mineral fabric (^{Klaver, Desbois, Urai & Littke, 2012}) and to its oil-bearing and
hydrological properties (^{Silin, Jin & Patzek,
2004}).

Likewise, brine inclusions in sea ice can be formalized in terms of their morphology
and pore connectivity (^{Golden et
al., 2007}). Silica sands need to have a high sphericity, in that
where the rounder and more spherical the particle is, the more resistant that
particle is to crushing or fragmenting (^{Schroth,
Istok, Ahearn & Selker, 1996}). In engineering, the durability of
cementitious materials is associated with certain mechanical and transport
properties that can be evaluated based on the properties of the pore space (^{Stroeven & Guo, 2006}).

There are several methods for describing porous materials. Some of them have been traditionally computed using 2D methodologies, and their results extended to 3D by means of stereologic techniques, while other parameters can be computed directly from the volume dataset.

Background

Pore-size distribution

Pore-size distribution is usually computed with mercury intrusion porosimetry (MIP), an in-vitro experiment based on the capillary law governing liquid penetration into small pores. In this technique, mercury is intruded into the sample at increasing pressures, causing the fluid to flow through smaller apertures.

Virtual porosimetry is an in-silico experimentation for general porous
materials that simulates MIP at incremental pressures (^{Garboczi et al., 1999}). All separated
components of the pore space filled at a given pressure are considered pores
and labeled with the diameter corresponding to the applied pressure. This is
a flood-fill methodology that uses a previously computed 2D (surface)
skeleton as a guide to simulate mercury intrusion from the entry points into
the pore space. Moreover, the skeleton is labeled at each point with the
corresponding distance (the maximum radius) used to allow or stop the
mercury intrusion simulation, i.e., when this radius is smaller than the
radius for the current pressure, the simulated intrusion at this pressure
stops. Regions corresponding to these smaller radii are the throats. This
method is applied to biomaterial samples (^{Vergés, Ayala, Grau & Tost, 2008b}) and soil samples (^{Delerue, Lomov, Parnas, Verpoest & Wevers,
2003}). The latter uses the obtained pore network to compute the
permeability. Several methods define the initial set of entry points as
those pore voxels which are connected to exterior; however, it can be
predefined to allow more freedom to the user in the simulation.

There are other methods that do not require skeleton computation to simulate
MIP (^{Rodríguez, Cruz, Vergés & Ayala,
2011}); instead, they use an iterative process that considers the
diameters corresponding to pressures. At each iteration, geometric tests
detect throats for the corresponding diameter, and a connected-component
labeling process collects the region invaded by the mercury.

There are other approaches to compute the pore-size histogram and the pore
graph. Some are related to the shape-analysis discipline and use a 2D
skeleton as a tool to devise the shape and size of pores. These approaches
are heuristic methods that cover the pore space with overlapping spheres, so
that the pores are computed as the unions and differences of maximal spheres
centered at skeleton points. These methods can be applied to sand samples
(^{Silin et al.,
2004}), bone scaffolds (^{Cleynenbreugel, 2005}), and biomaterial samples (^{Vergés et al.,
2008a}).

Other approaches based on the 1D (curve)-skeleton computation detect throats
as the absolute minima of the skeleton. Thus, pores are defined as the
regions limited by throats and solid space (^{Lindquist & Venkatarangan, 1999}). A graph is obtained
directly from the 1D skeleton, in which nodes correspond to pores and edges
to throats (^{Liang, Ioannidis & Chatzis,
2000}). Alternatively, the minimal cost paths connecting boundary
points can be computed instead of skeletons to allow methods based on
porosimetry, as well as sphere positioning, to be applied (^{Schena & Favretto, 2007}).

Another technique for computing pore-size histograms is granulometry. This
methodology consists in the application of successive morphological openings
and does not require the computation of a skeleton (^{Cnudde, Cwirzen, Maddchaele & Jacobs, 2009}; ^{Hilpert & Miller, 2001}; ^{Schulz, Becker, Wiegmann, Mukherjee &
Wang, 2007}; ^{Vogel, Tölke, Schulz,
Krafczyk & Roth, 2005}). A specific set of discrete spheres
with diameters *D* such that
*S*(*D*) ⊂
*S*(*D*+1) can be used, thereby ensuring
that a larger ball will never reach a cavity in which a smaller ball cannot
enter (^{Hilpert & Miller, 2001}).
Physical MIP has been compared with a granulometry-based method (^{Cnudde et al.,
2009}).

Throats can also be detected based on the negative Gaussian curvature of the
surface. This method is used to decompose the solid space into single
particles and to compute mechanical properties such as grain size and
coordination number (number of contacts per grain) (^{Theile & Schneebeli, 2011}).

Representation Models

In most of the reported literature, the operations to study binary volume datasets are performed directly on the classical voxel model. However, in the field of volume analysis and visualization, several alternative models have been devised for specific purposes.

Hierarchical structures such as octrees and kd-trees have been used for Boolean
operations (^{Samet, 1990}),
connected-component labeling (CCL) (^{Dillencourt,
Samet & Tamminen, 1992}), and thinning (^{Quadros, Shimada & Owen, 2004}; ^{Wong, Shih & Su, 2006}). Octrees are used as a means of
compacting regions and getting rid of the large amount of empty space in the
extraction of isosurfaces (^{Wilhems & Gelder,
1992}). Their hierarchy is suitable for multi-resolution when dealing
with very large data sets (^{Andújar, Brunet &
Ayala, 2002}; ^{LaMar, Hamann & Joy,
1999}), as well as to simplify isosurfaces (^{Vanderhyde & Szymczak, 2008}). Kd-trees have been used
to extract two-manifold isosurfaces (^{Greß &
Klein, 2004}).

There are models that store surface voxels, thereby gaining storage and
computational efficiency. The semi-boundary representation affords direct access
to surface voxels and performs fast visualization and manipulation operations
(^{Grevera, Udupa & Odhner, 2000}).
Certain methods of erosion, dilation and CCL use this representation (^{Thurfjell, Bengtsson & Nordin, 1995}).
The slice-based binary shell representation stores only surface voxels and is
used to render binary volumes (^{Kim, Seo &
Shin, 2001}).

In this paper, the Compact Union of Disjoint Boxes (CUDB) decomposition model is used. With this model the geometry is partitioned in a non-hierarchical, sweep-based way. An algorithm has been devised to detect the throats based on this model.

OPP and the CUDB model

A binary voxel model represents an object as the union of its foreground voxels,
and its continuous analog is an orthogonal pseudo-polyhedra (OPP) (^{Lachaud & Montanvert, 2000}). OPP have
been used in 2D to represent the extracted polygons from numerical control data
(^{Park & Choi, 2001}). Some 3D
applications of OPP are: general computer graphics applications, such as
geometric transformations and Boolean operations (^{Bournez, Maler & Pnueli, 1999}; ^{Esperança & Samet, 1998}); skeleton computation (instead of
iterative peeling techniques) (^{Eppstein &
Mumford, 2010}; ^{Martı́nez, Pla &
Vigo, 2013}); orthogonal hull computation (^{Biedl & Genç, 2011}; ^{Biswas, Bhowmick, Sarkar & Bhattacharya, 2012}); boundary
extraction (^{Vigo, Pla, Ayala & Martínez,
2012}); and model simplification (^{Cruz-Matías & Ayala, 2014}). OPP have been also used in theory of
hybrid systems to model the solutions of reachable states (^{Bournez et al., 1999}; ^{Dang & Maler, 1998}).

The Compact Union of Disjoint Boxes (CUDB) (^{Cruz-Matías & Ayala, 2017}) is a special kind of spatial
partitioning representation, where an OPP is decomposed in a list of disjoint
boxes. These models can be obtained from the voxel model and other alternative
models. The conversions algorithms have been published (^{Cruz-Matías & Ayala, 2017}).

Let *P* be an OPP and *Π*
_{
c
} a plane whose normal is parallel (without loss of generality) to the X
axis, intersecting it at *x = c*, where *c = C*
_{1},...,*C*
_{
n
} . More formally, *δ* is an arbitrarily small quantity. Then, *P* and *C*
_{
s
} such that *C*
_{
i
}
*< C*
_{
s
}
*< C*
_{(i+1)}, is called a section of *P*. Figure 2 shows an OPP with its cuts and
sections perpendicular to the X axis, since work is done with bounded regions,
*S*
_{0} (*P*)=∅ and *S*
_{
n
} (*P*)=∅, where *n* is the total number of
cuts along a given coordinate axis.

An OPP can be represented with a sequence of orthogonal prisms represented by their section. Moreover, if the same reasoning is applied to the representative section of each prism, an OPP can be represented as a sequence of boxes. CUDB represents an OPP with such an ordered sequence of boxes in a compact way, as many boxes generated by the aforementioned split process are merged into one in several parts of the model. CUDB is axis-aligned like octrees and bintrees, but the partition is done along the object geometry as in binary space partitioning (BSP). Depending on the order of the axes, chosen to split the data, a 3D object can be decomposed into six different sets of boxes: XYZ, XZY, YXZ, YZX, ZXY, ZYX, and the set will be ordered according to the chosen configuration. Figure 3 illustrates two possible decompositions of the model in Figure 2 (left). In 2D, an object can be decomposed into two different sets: XY and YX.

In the CUDB model, the adjacency information of the boxes is stored. Each box has
neighboring boxes in only two orthogonal directions: A- and B-direction, and for
each one there are two opposite senses; so, four arrays of pointers to the
neighboring boxes (two for each direction) are enough to preserve the adjacency
information that is required for future operations. These arrays are defined as
A-backward neighbors (ABN), A-forward neighbors (AFN), B-backward neighbors
(BBN) and B-forward neighbors (BFN). For more details of this model see (^{Cruz-Matías & Ayala, 2017}).

Methodology

Geometric Pore Space Partition

Unlike traditional approaches, which require the skeleton, the method adopted in this study uses the 2D CUDB-encoded pore space as input. It consists of two different sorted sets of boxes for 2D samples, in which a box represents a cluster of pixels.

The basic steps of our algorithm are:

Remove Unreachable Regions

As in real porosimetry, only those components connected to the exterior (reachable components) are able to be filled, so totally interior components (cavities) will not be considered.

In order to remove the unreachable regions, a CUDB-based connected component labeling is applied, where those boxes that are connected to the exterior are labeled with the same label; boxes with a different label are removed. Figure 4 illustrates the removal process with a 2D example.

2D Throats Detection

In the method used in this study, throats can be orthogonal or oblique and all of them are shaped as lines (Figure 5). In order to detect orthogonal and oblique throats, the pore space must be exhaustively scanned in the two 2D main directions.

Orthogonal Throats Detection

Boxes in 2D-CUDB have neighboring boxes only in the A-direction; therefore, orthogonal throats can exist in any of these directions. As the main axis that splits the model is A, orthogonal throats are detected in this direction. Thus, the two AB-sorted CUDB encodings of the pore space are required to detect all orthogonal throats.

Let *β*
_{
i
} and *β*
_{
j
} be two A-adjacent boxes and ^{1} of
their projections, respectively, over a segment perpendicular to the
A-coordinate. There is an orthogonal throat between *β*
_{
i
} and *β*
_{
j
} if

Note that the orthogonal throat is the intersection of two adjacent boxes, as shown in Figure 5 (left) and Figure 5 (middle). Algorithm 1 shows the steps of the orthogonal throat detection process for a single AB-sorting (Figure 6).

Oblique Throats Detection

Note that an oblique 2D throat is an oblique segment. The next process generates the oblique throats for an AB-sorted CUDB which, due to the characteristics of the CUDB model, are the same generated for the corresponding BA-sorted CUDB.

Let *β*
_{
i-1
} , *β*
_{
i
} and *β*
_{
i-1
} be three consecutive A-adjacent boxes in an AB-sorted CUDB, and let *β*
_{
i
} projection over a segment perpendicular to the A-coordinate. Then, an
oblique throat exists between *β*
_{
i-1
} and *β*
_{
i+1
} if the next conditions are satisfied (Figure 7 (a) illustrates these conditions):

However, oblique throats are not restricted to three consecutive boxes. Depending
on the geometry of the object, an oblique throat can occur between two boxes,
*β*
_{
S
} and *β*
_{
T
} , which have more than one intermediate box between them, as shown in
Figure 6 (b) . In these cases, all the
previous conditions must be satisfied for *β*
_{
S
} and *β*
_{
T
} instead of *β*
_{
i-1
} and *β*
_{
i+1
} , respectively. Therefore, these conditions are rewritten as:

Condition 1 means that *β*
_{
S
} and *β*
_{
T
} . In addition, Figure 2 shows why it
is necessary to consider open sets for the above conditions. Otherwise, the
aforementioned conditions would be satisfied, and the algorithm would detect
oblique throats where none exists.

The previous conditions are used to detect an oblique throat between
*β*
_{
S
} and *β*
_{
T
} . Although a rectangle is obtained (see the rectangle highlighted in
orange color in Figures 7a and 7b), the strategy to compute the throat is by
considering the diagonal included in the mentioned rectangle (see the inclined
line highlighted in yellow color in Figures
7a and 7b). This strategy allows
to separate the pore space in two regions in the current scan direction.

A box *β*
_{
S
} can generate several oblique throats with its consecutive boxes in the
A-direction (Figure 8a). Therefore, in
order to detect all the possible oblique throats, all consecutive boxes of
*β*
_{
S
} must be scanned. The method used here searches for any configuration that
satisfies the oblique throat conditions.

Conditions 1 and 3 are required conditions to stop the search process, while the
others are not. The search stops in case that Condition 1 is not satisfied, as
every subsequent box *β*
_{
T
} will not satisfy Condition 4 in the following iterations (Figure 7). When Condition 1 is not satisfied,
every subsequent box *β*
_{
T
} will not either satisfy Condition 2 or 4 in the following iterations
(Figure 7). Moreover, the search must
be also stopped when the distance between *β*
_{
S
}
*.V*1 and *β*
_{
T
}
*.V*1 in the A-coordinate is greater than a certain distance
*D*, this distance has been defined as the smaller dimension
of the model. However, when conditions 2 and 4 are not satisfied, there may
still be subsequent boxes that define throats that satisfy all conditions (Figures 8d and 8e).

Algorithms 2 and 3 show the steps of the oblique throat detection using a recursive strategy for a single AB-sorting (Figure 9 and 10). Figure 11 shows an example of throats computation. In this Figure, black and white represent the pore and solid spaces respectively.

Partitioning of the Pore Space

An accurate identification of pore throats is important in the partitioning of the pore space. It is necessary to be careful not to generate too many throats which can produce a large number of small pores, which would correspond primarily to the roughness of the pore wall interface.

From the point of view of the utility of the resulting pore and throats size
distributions in network simulations, it is desirable to partition the pore
space into nodal pores (pores that are also topological nodes) separated by
throat cross-sectional surfaces (^{Ioannidis,
Chatzis & Kwiecien, 1999}). This is possible if pore necks are
identified as the minimum distance that connect the solid space. In methods like
those by ^{Liang et al.
(2000)} and ^{Lindquist &
Venkatarangan (1999)}, throats are identified through a search for
minima in the hydraulic radius of individual pore space, where each branch of
the skeleton (link) is considered as a separate pore.

Although the method used follows the previous idea, the authors of the study do not rely on prior computation of the skeleton; instead, all the possible throats generated with the strategy described below are computed, then, such throats are analyzed in order to obtain those that better split the pore space.

The strategy to split the pore space is to determine the shortest throats that join the solid regions. Therefore, the first step consists in sorting the detected throats by length. To obtain the solid regions, a CCL process is applied to the CUDB-encoded solid space. For instance, Figure 11 shows that this sample has 10 solid regions.

Once the throats have been sorted and the solid regions have been detected, it is
necessary to determine the final throats that will split the pore space. Let
*S*
_{
A
} and *S*
_{
B
} be two solid regions, each of the ordered throats is analyzed to obtain
the pair (*S*
_{
A
}
*, S*
_{
B
} ) joined by this throat. This detection is done by evaluating each box of
each CUDB-encoded solid region. However, this process can be speeded up by using
the bounding box to do previous discards.

Figure 12 (left) shows the resulting
throats after the previous strategy. Note in this sample that although there are
no throats joining the same solid regions, several intersected throats exist.
Some of these throats are longer than some paths formed by other throats. This
problem can be solved by applying a Dijkstra's algorithm (^{Dijkstra, 1959}). Every time a throat is found for the pair
(*S*
_{
A
}
*, S*
_{
B
} ), it is examined whether there is a shorter path formed with previously
calculated throats; if this is the case, the throat is omitted. Implementation
of these criteria is shown in Figure 12
(middle).

Once the final throats have been detected, the object is subdivided along each one of them to partition the pore space into its constituent pores. To achieve this subdivision, a unit-width orthogonal polygon 2D is built for each throat. This OPP is obtained by the 2D scan-conversion of the segment which represent the throat. Therefore, the difference between the OPP-object and the OPP-throat produces the object subdivided along the throat. Figure 13 illustrates the subdivision of the object along two oblique throats.

Once the pores have been correctly separated, their area is straightforward computed using the CUDB-encoded pore.

Results and discussion

The method presented here, like those based on prior computation of the skeleton, is a geometric approach that computes the expected solution.

Figure 14 shows a well-defined structure for which it is easy to determine the theoretical solution. This model has a size of 450 x 450 pixels and contains squares of 85, 58 and 30 pixels per side connected with small lines (throats) of size 2 and 4 pixels. Note in this figure that the method is invariant to rotation to detect the throats and to divide the pore space correctly. For more clarity, the pores after the subdivision along the throats have been colored.

The adopted approach has been compared with a skeleton-based approach. The model used
to explain the partitioning process is the same used in ^{Liang et al. (2000)}. In Figure 9 (middle and right) both methods get a quite similar
partitioning of the pore space.

Also, the effects of resolution on the method used here have been tested. Resizing the image does not affect the topology of the detected throats. Although the bigger the image, the more throats are generated, the method to detect the final throats allows to obtain an almost identical pore space partition. Figure 15 shows the minor effects with increasing resolution. Of course, the accuracy of measurements is expected to increase as resolution increases.

Finally, sections of several real porous samples have been tested, as shown and described in Table 1 below. Figures 16 to 19 depict the corresponding partitioned porous space and its pore-size distribution.

Spheres: A synthetic model with spheres of a variable size.

PLA: A biomaterial sample consisting of polylactic acid (PLA) with a 16 µm

^{3}voxel resolution.Soygel: Sample of hydroxyapatite with a soy-based foaming agent for bone prosthetics with a 16 µm3 voxel resolution.

Stone: A stone sample consisting of sedimentary rock from a Lybian oil-bearing unit with a 4.4 µm

^{3}voxel resolution.

Dataset | Size | Time (ms) |

Spheres | 450 x 540 | 96 |

PLA | 314 x 237 | 56 |

Soygel | 251 x 251 | 843 |

Stone | 271 x 179 | 25 |

Source: Authors’ own elaboration.

All the real samples have been scanned by Trabeculae® and segmented by thresholding and applying noise filtering. The corresponding programs have been written in C++ and tested on a PC Intel®Core i7-4600M CPU@2.90GHz with 7.6 GB RAM and running Linux.

Based on the analysis of the figures and the histograms, it can be concluded that the obtained pore size distribution is accurate enough. A more detailed discussion of the comparison of physical and in-silico experimentation is beyond the scope of this paper.

Conclusions

A new method to characterize the pore space has been proposed; it does not require prior computation of the skeleton. As this prior preprocessing is very time-consuming, this approach achieves a noticeable reduction in time. The key feature of the presented method is the process to determine the throats that make it possible to separate regions corresponding to different pores.

The in-silico approaches are based on the geometry of the analyzed samples, and the results obtained are the expected theoretical results, as shown in the results section.

As a future work, several topics must still be studied. At present, the authors are starting to study how to extend this method directly to 3D models and will make contact with biomaterials and geology research teams that provide new problems to study and other kinds of real samples. Further, new methods will be developed for other parameters, such as connectivity and anisotropy. Moreover, the authors of this research are also interested in applying time-varying techniques based on a CUDB model.