# Throat Finding Algorithms based on Throat Types

The three-dimensional geometry and connectivity of pore space determines the flow of single-phase incompressible flow. Herein I report on new throat finding algorithms that contribute to finding the exact flow-relevant geometrical properties of the void space, including high porosity samples of X2B images, three-dimensional synchrotron X-ray computed microtomographic images, and amounting to over 20 axis that comes from the 3DMA-Rock software package. To find accurate throats, we classify three major throat types: mostly planar and simply connected type, non-planar and simply connected type, and non-planar and non-simply connected type. For each type, we make at least one algorithm to find the throats. Here I introduce an example that has a non-planar and simply connected throat, and my solution indicated by one of my algorithms. My five algorithms each calculate the throat for each path. It selects one of them, which has the smallest inner area. New algorithms find accurate throats at least 98 samples (over 20 image. The new calculation uses three mathematical concepts: i) differentiability, ii) implicit function theorem, iii) line integral. The result can convert the discrete boundary of the XMCT image to the real boundary. When the real boundary has an arc shape, the new calculation has less than 1

## Authors

• 1 publication
• ### Genus Computing for 3D digital objects: algorithm and implementation

This paper deals with computing topological invariants such as connected...
12/25/2009 ∙ by Li Chen, et al. ∙ 0

• ### Finding Tutte paths in linear time

It is well-known that every 2-connected planar graph has a Tutte path, i...
12/11/2018 ∙ by Therese Biedl, et al. ∙ 0

• ### Practical I/O-Efficient Multiway Separators

We revisit the fundamental problem of I/O-efficiently computing r-way se...
07/06/2021 ∙ by Svend C. Svendsen, et al. ∙ 0

• ### The Weisfeiler-Leman Dimension of Planar Graphs is at most 3

We prove that the Weisfeiler-Leman (WL) dimension of the class of all fi...
08/24/2017 ∙ by Sandra Kiefer, et al. ∙ 0

• ### Ray-marching Thurston geometries

We describe algorithms that produce accurate real-time interactive in-sp...
10/29/2020 ∙ by Rémi Coulon, et al. ∙ 0

• ### On multiple connectedness of regions visible due to multiple diffuse reflections

It is known that the region V(s) of a simple polygon P, directly visible...
06/02/2003 ∙ by Sudebkumar Prasant Pal, et al. ∙ 0

• ### The Shape of an Image: A Study of Mapper on Images

We study the topological construction called Mapper in the context of si...
10/24/2017 ∙ by Alejandro Robles, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Understanding the precise structure of pore space is integral to fully understanding flow through porous media. This understanding is essential in research on single-phase incompressible flow in situations where another dominant flow factor is the fluid-wall surface tension [1]. Determining the relationship of pore structure and bulk flow properties is complicated since the structure of pore space is so very random and awkward because of the challenge of precisely measuring pore space.

A throat is the smallest corss-section area that corresponds a branch-branch medial axis path. For non-crossed throats, their outer perimeter voxels have to exist on the boundary grain voxels. The 3DMA-Rock software package [2] has three major throat-finding algorithms; (1) the wedge-based algorithm [3],(2) the Dijkstra-based shortest length algorithm [1], and (3) the planar dilation algorithm [4]. The wedge-based algorithm yields acceptable measures in only low porosity samples. It uses wedges to find the nearest boundary grain voxels and connects voxels in a way that allows the determination of border voxels which allow the connection of each wedge. For high porosity samples (over 30), It is very diffcult to find the connedcted throat-perimeter paths on the boundary grain voxels. The primary deficiency of the Dijkstra-based shortest length algorithm is that the shortest path perimeter does not reveal the smallest cross area (Figure 1). The main deficiency of the planar dilation algorithm is in finding approximated throats. In general, the throat shape is a non-planar type. The throat determined by the employment of the planar dilation algorithm is capable of spreading pore space without any restriction.

The goal in this paper is to find the accurate throat area (Figure 2 right panel) and calculate the circumference based on known mathematical formulae. Here, I present new throat-finding algorithms using the concepts in vector space and spherical coordinate system (decribed in section 2). In order to find an accurate throat, I need to classify throat shapes into three main types; (1) the simply-connected planar type, (2) the simply-connected non-planar type, and (3) the non-simply-connected non-planar type. For each thorat shape, I construct five algorithms; two of them for the simply-connected planar type, two of them for the simply-connected non-planar type, and one of them for the non-simply-connected non-planar type. The first algorithm is designed to reduce the total calculation time. It indicates the simply-connected planar throat perimeter, and only uses the visible path (the member voxels of the path are visible from the MA path, invisible path is not) from the pertinent MA path. Usually the perimeter (the solution of the first algorithm) is not a throat, but we can reduce the total calculated CPU time (involving 4 other of our algorithms) using this possible candidate perimeter. The second algorithm indicates the simply-connected planar throat perimeter. This second algorithm uses visible paths and invisible paths (from the medial axis voxel) that together make an enclosed loop. All 26-connected perimeter voxel paths are on the plane. The third algorithm is appropriate for determining the simply-connected non-planar throat perimeter. The perimeter voxels, in cases where we employ this third algorithm, consist of visible voxels and invisible voxels from the pertinent medial axis path. The visible voxels are on the plane, and the third algorithm connects discontinuous voxels using non-planar voxels (generally, the voxel set is not on the same plane of visible voxel set). The fourth algorithm is designed to determine the largely undulate throat perimeter (See Fig. 2). This algorithm makes four vertical and horizontal wedges in each neighboring direction (based on the normal plane to the medial axis path). For each wedge, the third algorithm finds the nearest voxels to the MA voxel, and connects eight perimeter voxels using Dijkstra’s algorithm

[5]. The fifth algorithm is designed to determine non-simply non-planar throat. Also, this algorithm is useful when a small number of visible voxels diverge from the perimeter. After finding the visible continuous voxel set, the fifth algorithm tries to connect them. When this algorithm fails to connect voxels, we delete the voxels from the member of the candidate perimeter voxel set; we keep the path with the shortest distance from the MA path when there are two or more paths. To deal with discontinuous voxels in the perimeter voxel set, this algorithm uses Dijkstra’s algorithm. When Dijkstra’s algorithm finds several paths, our algorithm uses the innermost path from the medial axis voxel that results in affording the smallest throat area. These algorithms can detect more than 98 throats over higher than 29 porosity samples.

## 2 Throat-Finding Algorithms

### 2.1 Simply Connected Planar Throat-Finding algorithms

The first algorithm is to find throat on a plane, called planar throat type. In this algorithm, both medial axis modiification [2] and medial axis extraction [6], which computed within the void space of a segmented [7] 3D XCMT image, are used. For constructing the first throat-finding algorithm, I need to check all possible perimeters of a throat from the first to the last voxel in the path of the medial axis since the area of throat does not show continuously change. The area of a throat can be varied by the position and direction of the medial axis voxels in three dimensional space. Here, I present two simply connected planar throat-finding algorithms. The first algorithm can be used when the boundary grain voxel set is visible from the medial axis path. The second algorithm is used when some voxels in the boundary grain voxel set is not visible from the medial axis path.

#### 2.1.1 The first Simply Connected Planar Throat-Finding algorithm

The first algorithm is the following steps. Here, I define
(1) : -th voxel of the medial axis path,
(2) : unit tangent vector of the medial axis path,
(3) : unit directional vector where is a degree of polar angle from and is a degree of azimuth angle in spherical coordinate system
(4) : a plane that has the normal vector at the center of
(5) : unit directional vector on the plane

Step 1 : Compute the tangent vector at the center of voxel in the path

Step 2 : Make unit directional vector set of with zenith vector

Step 3 : Make a plane with a point at the center of and a normal verctor

Step 4 : Make 360 of unit directional vector set on

Step 5 : Find the boundary grain voxel set using (Finding method will be discussed in below)

Step 6 : Check whether the connectivity of the boundary grain voxel set from Step 5 has 26 connectivity or not

Step 7 : If the connectivity of the boundary grain voxel set shows 26 connectivity then extend it to 6 connected closed loop voxel set. Otherwise, I use different throat-finding algorithms (these algorithms will be discussed in section 2.1.2 and section 2.2.1)

Step 8 : Calculate the area of the closed loop voxel set using 26 connected perimeter voxel set from Step 6.

Step 9 : Repeat from Step 1 to 7 for each to get local minimum throat

The main idea of this algorithm is to use the visible boundary grain voxels from the medial axis voxel on a given plane. All voxels originate as specific points within each coordinate. For the medial axis voxel , , where is the last voxel number in the medial axis path, we choose 5 consecutive medial axis voxels (the four proximal – at right and left, two at each side) to find the zenith direction , which is the solution of the least square problem of the 5 involved voxels. If the medial voxel position is located within the two outermost voxels on each side, I choose the 5 outermost in the medial axis path that includes the medial axis voxel (Figure 4 B). If the total number of voxels in the path is less than 5, I use them to calculate the directional vector solution of the least square problem. Using the center of = and = , I can make the sequence (Figure 4 C), defined by the polar and azimuth angles, and , =1,2,, 45, = 1,2,, 45, including . The center of and its normal vector create the plane. I now need to make two unit vector sets: and where = (Figure 4 D). The ray set on the plane is normalized by the orthogonal projection of the unit vector set . The unit ray set has an angle difference of , and is on the coordinate plane that consists of the two coordinate axes that have two small absolute components of . The set uses the formula of one of the three different sets:

 {(cos(θo),sin(θo),0)},{(cos(θo),0,sin(θo))},{(0,cos(θo),sin(θo))}

where = .

I choose the set when and . The two others are selected similarly. The directional vector = starts from the center of and spreads until it touches the respective boundary grain voxel (Figure 5). Each ray yields three coefficient sets: , , and , = , and = 1,2,3 that satisfy = , and , when is not zero.If is zero, then the ray doesn’t go in the β coordinate direction. Those coefficients represent the boundary points that touch the voxels on the ray. Using the mixed set of , , and with increasing order, we can determine the voxel order that the pertinent ray touches. The order of the new set, at each point, indicates the direction of the voxel movement from If two coefficients among the mixed set are the same, then the ray touches the edge of the next voxel, with the related coordinates indicating the coefficients from the current position. If three coefficients are the same, then the ray touches the vertex of the next voxel, with the related coordinates indicating the coefficients from the current position. When the ray touches an edge or a vertex, we only check the next position to see whether the voxel is grain or void on the ray, because the void space has 26- connectivity. When we construct a throat barrier with 6-connectivity using a ray, we add all touching void voxels. If a ray touches an edge and all contiguous voxels are void, then the three void voxels will be counted among the throat barrier voxels. If a ray touches a vertex and all contiguous voxels are void, then the seven void voxels will be counted among the throat barrier voxels. The ray set is used to search for the boundary grain voxels that are the candidates for the perimeter voxel set on (Figure 4 E). Usually finded grain voxel set has 26-connectivity, so we change the connectivity to the grain’s (6-connectivity). The set is comprised of blue voxels to the 6-connected set by adding the neighboring diagonally positioned boundary grain voxels on the plane. Sometimes an added voxel is not on . The candidate perimeter voxel set (Figure 4 F). The set is comprised of green Voxels) consists of piecewise 6-connected parts. When the 6-connected perimeter voxel set makes a closed loop, we calculate the area with the 26-connected voxel set. After we find all possible perimeters for all of the plane , we compare all of them to find a possible candidate throat perimeter.

#### 2.1.2 The Second Simply Connected Planar Throat-Finding algorithm

The second planar throat-finding algorithm follows the same process as the first, which constructs the 6-connected perimeter voxel set. This algorithm calculates the triangular area using the centers of the finded 26-connected boundary grain voxel set and . This algorithm allows the progression to the next step, when the summation of the triangular areas is smaller than the minimum area of the first planar throat-finding algorithm (Figure 6 B). This algorithm finds invisible boundary grain voxels (for the unseen, distended perimeter) using Dijkstra’s algorithm. The candidate perimeter voxel set (Figure The set is comprised of green Voxels) consists of piecewise 6-connected parts. For each discontinuous path (Figure 6 C Two orange voxels), we use Dijkstra’s algorithm to connect the end path voxels to the boundary grain voxels. This connection can be understood conceptually in the following way: The boundary grain voxels are on the plane . The basic criteria used to determine that a grain voxel is on is * , where , , and are 0 or 1, and at least one 1 is included. The region of interest concerning Dijkstra’s algorithm is conceptually within a cube of variable side length. The side length changes from to , where is the side length of a minimum inclusion cube sharing the same center (Figure 6 C is the side length of the rose square). For each discontinuity, if there is no voxel connection in the cube, we assume it is a non-perimeter and we try to apply other algorithms. When we have found the perimeter voxel set that has 6-conectivity making a single and efficient encirclement (Figure 6 D Green voxels), we calculate the area with the 26-connected perimeter voxel set to find the minimum surface area using 3DMA. The 26-connected perimeter voxel set is selected within the 6-connected perimeter voxel set. After we find all possible perimeters for all of the plane , we compare all of them to find a possible candidate throat perimeter. Our algorithm has only one rounding number, and applies Dijkstra’s algorithm to the small discontinuous consecutive voxels.

### 2.2 Simply Connected Non-Planar Throat-Finding algorithms

I made two throat-finding algorithms for the undulating throat type. The first undulating throat-finding algorithm is designed to find simply connected throats. Sometimes, this algorithm cannot find the smallest throat area, but it can indicate the throat perimeter when the other three throat-finding algorithms fail to find the perimeter.

#### 2.2.1 The first Simply Connected Non-Planar Throat-Finding algorithm

The first non-planar throat-finding algorithm is the same as the planar throat-finding algorithm, except concerning the region of interest in applying Dijkstra’s algorithm in dealing with the nearest boundary grain voxel from medial axis voxel for each step. For this algorithm, we select as the possible candidate the perimeter of a (conceptual) cube where the side length is the minimum cube that envelops two voxels with the same center in 3D (Figure 7 A). Now, 3DMA calculates the area with the 26-connected perimeter (Figure 7 B), derived from the 6-connected perimeter. After we find all possible perimeters for all of the plane , we compare all of them to find a possible candidate throat perimeter.

#### 2.2.2 The second Simply Connected Non-Planar Throat-Finding algorithm

The second non-planar algorithm can be used to analyze largely undulate type throats. We allow that the maximum angle of the fluctuation of this algorithm is from the base. Our assumption is that some perimeter voxels are situated nearest to an MA voxel. If the throat undulates a lot and satisfies our assumptions, then this algorithm can deduce its perimeter. Before constructing the plane , we use the same procedure as the plane of the planar algorithm. We draw the 91 unit ray set , θ=0,1,⋯,90, on the plane , drawn from the projections and normalization of the 91 unit ray set , = , on the coordinate axis plane (Figure 8 A). The set is selected among the three different sets:

 {(cos(θo),sin(θo),0)},{(cos(θo),0,sin(θo))},{(0,cos(θo),sin(θo))}

where = .. This 90-degree angle makes the first wedge. The selected order of the set is the same order of the maximum component value of . The set on is generated by the orthogonal projection of onto the plane. Using the two related ray sets and , we find the nearest grain boundary voxel and its angle ( = ) from . The voxel is on and in the first wedge with the angle. The next step is making three horizontal wedges on (Figure 8 B). Their center angles are , , and and their respective wedge angles are each . The three 61-unit ray sets are selected among (Figure 8 B). The selected angle sets are below:

 {θ2|θ2=ϵ+60o,ϵ+61o,⋯,ϵ+120o}, {θ3|θ3=ϵ+150o,ϵ+151o,⋯,ϵ+210o} {θ4|θ4=ϵ+240o,ϵ+241o,⋯,ϵ+300o}

For each horizontal wedge, we find the three nearest boundary grain voxels , , and their angles (=2,3,4) from on using 61 rays (Figure 8 B). We add four vertical wedges with center angles ,=1,2,3,4, where the angles are bisected between each angle of the set , = 1,2,3,4. Each vertical wedge’s azimuth angle varies from to from the new zenith . We make a 91-ray set in the vertical wedge, where each unit vector’s azimuth angle varies from to (Figure 8 C). The first vertical unit ray set is , where and =. The other three vertical unit ray sets are chosen similarly. Using four 91 ray sets, we find the four nearest boundary grain voxels , , , and to . There are eight candidate perimeter voxels. The final step is to connect those 8 voxels using Dijkstra’s algorithm in dealing with the nearest boundary grain voxel from the medial axis voxel for each step (Figure 8 D). If there is a complete 6-connected perimeter, 3DMA calculates the area with a 26-connected perimeter(Figure 8 D).

### 2.3 Non-Simply Connected Non-Planar Throat-Finding algorithm

The MA algorithm only preserves the topological structure of the given void space, so a throat’s type cannot be simply connected sometimes. For example, when there is a large peak extending from the bottom. If the path length is shorter than the peak, and the peak contains the projection of the path and its cluster voxels, then the throat has to be considered a non-simply connected type. This algorithm finds some other throat types, and is excellently suited to do so. Also, this algorithm can be used when perimeter parts indicate other pores’ boundaries. At this time, the wrong path parts can be nixed. When the throat type is the same as the previous type, the fourth algorithm can find similar perimeters. For the non-simply-connected non-planar throat type, we use the same procedure to find the plane and the 6-connected candidate perimeter set as the second planar algorithm. The set The 6-connected voxel set is comprised of blue voxels). The set consists of piecewise continuous paths (Figure 9 C). This algorithm starts with a continuous path among the set that constitutes the shortest voxel distance in the set from . Also, the starting voxel can be selected as the longest path or the largest angle. The algorithm connects each piece using Dijkstra’s algorithm in dealing with the nearest boundary grain voxel from ma voxel ν for each step. When this, our forth algorithm, fails to connect two pieces (Figure 9 D), the algorithm deletes a part of the set (Figure 9 D brown voxels). By deleting the discontinuous path, we connect to the next path (Figure 9 E). This algorithm’s candidate perimeter consists of the compositions of the piecewise continuous parts of ( is a visible perimeter voxel set, and is 6-connected voxel set that comes from ) and their connection (Figure 9 F). Sometimes this algorithm deletes lots of continuous paths of , so it requires checking the rounding number. The rounding number is defined as the sum of the angles of the projected perimeter voxels on . The perimeter has to have one rounding number. Now, the algorithm changes the perimeter to 26-connetivity, and 3DMA calculates the outer area. If the set - i=1,2,⋯, the total number of the perimeter voxels - is of the outer perimeter voxels of this type, the subtraction of the inner grain voxels’ area is necessary to calculate the throat area. We make a unit ray set , i=1,2,⋯, the total number of the perimeter voxels. A unit ray  ̵̂ indicates the respective position of the perimeter voxel relative to the MA voxel. Each  ̵̂ spreads to the perimeter voxel from . Our algorithm works when \$̂ go through the grain voxels (Figure 9 G), our algorithm saves the initial and final positions. We subtract the quadrilateral area from the outer area (Figure 9 H).

## 3 Point Set algorithm

The perimeter is used for the simulations of the drainage processes. The entry condition equation involves the area and perimeter. Digitized images are constituted of discrete values for each voxel. After segmentation, we only know from the data whether the voxel is grain or not. The segmented images consist of cubes, so we have to think about the difference between our calculations from the segmented images and the real length and area. Construction of the polygon (which is the boundary of the object in the CT image) by linear connection of the mid-points yields the relative error of its length for any random boundary. The reason for the relative error is that the mid-point connection contains only horizontal, vertical, and diagonal lines. A line integral needs the slope between two points =(,), =(+,+), so that the calculation expresses the real boundary length. The mid-point calculation of the length in CT images may yield a relative error of more than 8. Thus we needed to find a new method to calculate length with a new point set in the CT images. We use the 6-connected perimeter voxel set to calculate the exact length, and we use the 26-connected perimeter voxel set to calculate the area in 3D space. The mathematical concepts of the point set consist of: i) Differentiability, ii) Implicit Function Theorem, iii) Line Integral. For non-differentiable parts of the perimeter, we follow the inner part of the perimeter. When we apply a line integral to a CT image, we require only one mid-point on each instead of multiple points. (The extant mid-point method involves several horizontal and vertical points on each .) For each , the length is a multiple of the voxel size and all y values are the same. The mid-point locates the adjacent face of the perimeter voxel, and the point used is the middle of within part of the perimeter voxel set. We only use the perimeter voxel part that faces the 6-connected barrier voxel set in 3D space. To easily understand our algorithm, we here give an example for a line integral in a CT image in Figure 13 A and B. We have to know which parts of the perimeter voxel set are related to the mathematical concepts of differentiability, line integral, and implicit function theorem. Using these math concepts, we select the point set on the boundary between perimeter voxels and throat barrier voxels, because the perimeter voxels are located outside of the throat. We will divide the perimeter voxel set into several parts to apply these three concepts. The first step of the new algorithm distinguishes which perimeter parts are related to the differential and non-differential functions. We start with these. For the non-differential parts (U, T, and L shapes), we select points that identify the perimeter section shape. Those shapes consist of the combination of continuous horizontal and vertical voxels, and each part of the (respective) continuous vertical and horizontal part has at least 3 voxels. The first point is on the inner-facing part of the initial voxel, in the non-differential part in 2D space. From the second point on, points are selected: at each right angle in the non-differential perimeter part. Let an be a perimeter array and be a point set, and let be a subset of . U and T shaped voxels consist of horizontal-vertical-horizontal or vertical-horizontal-vertical voxel parts, and each part has the same x or y value in each of at least three voxels. To classify the U and T shaped voxels, a possible condition is 2, 1, and 2. The point (1 i 4) indicates vertex points (Figure 10 A and B). If the distances between the created points from each two sequential voxels are greater than or equal to 2 voxels, then the part satisfies the U and T shape conditions and four points are added into the point set. There are eight possible basic U and T shapes. The L shaped perimeter voxel part consists of horizontal-vertical or vertical-horizontal voxel parts, and each part has the same x or y value in each of at least three voxels. To classify the L shaped perimeter voxel part, one possible condition is 2 and 2. This process is similar to the process of finding the U and T shapes. Three points , , and are our vertex point. If the distances between the created points from each two sequential voxels are greater than or equal to 2 voxels, then the part satisfies the L shape condition and three points are added into the points set (Figure 10 C). The second step of the new algorithm distinguishes the perimeter part related to the implicit function theorem that will give us the C shape. The distinguishing condition of the C shape is that two edge voxels of the vertically altered perimeter part changes the y-coordinate direction to connect in the same direction in 2D space. Some of the U shapes are constituents of the C shape group. When a perimeter part satisfies the U and C shape conditions together, the U shape condition is given priority; we use that first. A point that is within the point set of a C shape voxel set related to implicit function theorem is selected in the middle of voxels that are vertical. The last classifying step is to select points related to the mathematical concept of the line integral. The line integral is understood as a linear connection of the points on the curve. For each ∆x that continues with the same y value in the perimeter voxel set, we select the mid-point – an 8-connected voxel set in 2D (26-connected voxels in 3D space). Each line integral part starts and ends with either a U, T, L, or C shape. When our understanding of the boundary involves the implicit function theorem, the boundary point uses a point in the C shape. When our understanding of the boundary involves the non-differential part, then the boundary point uses a point in any of the U, T, or L shapes.

## 4 Conclusions

### 4.1 Throat-Finding Algorithms

To find accurate throat, we calculate all possible areas with five algorithms. The results can classify the throat types. In natural samples, non-simply connected type is rare. We test our new algorithms with a hypothetical and then with real high porosity samples of more than 20. When two cylinders cross each other with a small cross area, the throat is saddle-shaped (Figure 11). Our algorithms suggest two possible throat candidates: a planar throat and a non-planar throat (Figure 11 A, B, and C). The throat we analyze is the smallest one (Figure 11 D), the actual. The second test sample is sphere pack with 9*9*9 spheres with a rectangular bounday. My algorithm can find all throats (Figure 12). Also, we test 8 of the S3 samples and 4 sections of the Handford wet inlet t241 samples. The results for these twelve tests are shown in Table 1. For cubic samples (512*512*512

), the total CPU time consumed for the throat finding algorithms and probability density function of throat area, pore volume, and coordination number is 4 to 10 hours. Our new algorithms include a crossed-throat algorithm [4]. The Pdf of throat area, pore volume, and coordination number is shown in Figure 13.

Table 1 : Throat type and success ratio of the new algorithm

### 4.2 Point Set Algorithm

To select a point set in 3D space, we project the real planar perimeter voxel set (in 3D space) onto the axial plane (2D) consisting of two coordinates, the lowest two absolute values among three components of the normal vector of the plane. To test our new algorithm, we apply our algorithm to lines with different slopes and circles with different radii. For the linear boundary, we draw the linear boundary , where =1,2,⋯,89, and the domain is from 0 to [50,000 ], where [x] is a maximum integer less than or equal to x. We make digitized images of void and grain voxel areas. The grain voxels are selected when we have more than half of the voxel area below the linear boundary. The mid-point algorithm’s maximum relative error for the linear boundary is 8.23, and the new algorithm’s relative error is 1.29 (Figure 14 C). The circular boundary has two perimeter parts that relate to the C shape, and the boundary doesn’t have any U, T, or L shapes. We can apply two line integrals which indicate the upper boundary and the lower boundary. The mid-point algorithm’s relative error oscillates around 5.5, but the new algorithm has less than a 1 relative error (Figure 14 D). For non-planar throat type, triangular throat area incluiding a medial axis voxel can over-calculate its area (Figure 15 The orange ellipse area is my interested area. B, C, and D are the same image with triangular throat calculation including a medial axis voxel. E is my preferred solution.)

Figures

Figure 1 : It shows that the short path does not always indicate the small area. The area (light blue color) of a throat in the right panel is smaller than the area in the left panel but the distance of perimeter of boundary voxel (green color) is opposite. (Brown color : grain voxel, green color : perimeter and boundary grain voxel and light blue color : throat barrier voxels)

Figure 2 : Example of Two cross cylinders: left figure is an example of non-planer throat type and the right figure shows the real throat (blue color) of left fiugre.

Figure 3 : Test example for Dijkstra’s algorithm and possible paths. Blue voxels are possible path grain voxels. Green color indicates the initial and terminal grain voxels. Left panel shows all possible paths using Dijkstra’s algorithms and numbers indicate step number. Middle panel shows an example of errant approach using Dijkstra’s algorithm with smallest tirangular area and Right panel illustrates real boudary grain path using Dijkstra’s algorithm with shortest distance for each step voxel.

 A B C D E F

Figure 4 : The first algorithm processing diagram. A : vertical subsection slide of 3D image (dark red, green and white colors indicate the grain, the medial axis, and void space, respectively). B : finding the tangent vector at the center of using five consecutive voxels. C : finding unit directional vector . This vector is an unit directional vector in spherical coordinate system whose zenith directions is in the same direction of . D : constructing a plane, with both a point (the center of ) and a normal vector and an unit directional vector set {} on the plane . E : finding the boundary grain voxel set using the unit directional vector set {}. F : changing the connectivity from 26 to 6 connectivity of the founded boundary grain voxel set from E.

Figure 5 : A unit ray spread in the digitized image. Yellow colored vector is in the direction of the unit vector and touches the boundary grain voxel.

 A B C D

Figure 6 : The second algorithm processing diagram. A : detects the visible boundary grain voxels (blue colored voxels). B : Finding triangular area using both two centers of boudary grain voxels and the center of medial voxel (red colored voxel in the middle). C : The minimum (pink rectangle) cube encompasses two non-connected (orange color) voxels. The side length of the cube is . The interested region is between + 4 and +8 to find accurate throats. D : Final boundary grain voxels using the second algorithm.

A B

Figure 7 : Detecting the boundary grain path in space using Dijkstra’s algorithm. A : Dijkstra’s algorithm is applied to interested region cube in space. B : detected boundary grain path in space.

 A B C D

Figure 8 : The processing of non-planar algorithm. A : detection of the nearest boundary grain voxel using where is between and . B : Detecting another nearest boundary grain voxel with different . C : make four verticle wedge (green ray) and dectect the boudary grain voxels not on the plane . D : Detected the boundary grain voxels in space.

 A B C D E F G H

Figure 9 : Non planar algorithm processing diagram. A: visible voxel sets (blue) from the medial axis on a plane. B : calculating the sum of triangular area using both the visible voxel sets and a medial voxel (). C : connect to disjoint consecutive voxel using Dijkstra’s algorithm. D : nixed voxel set (brown voxels). E : New Region of interest. F : Outer 26 connected perimeter voxel set. G : Find 4 points on the boundary of inner grain voxels. H : The inner grain boundary and the outer perimeter voxel set.

 A B C

Figure 10 : The method how to find points with respect to non-differential section. A : Perimeter voxels and its points set in the T shape. B : Perimeter voxels and its points set in the U shape. C : Perimeter voxels and its points set in the L shape.

 A B C D

Figure 11: Planar and non-planar throat types. Presented algorithms in this paper detect both planar and non-planar throat (red and light blue color in A and B), then I choose the ideal throat by taking the minimum of both detected planar and non-planar throats. Red color : planar throat barrier voxels. light blue color : non-planar throat barrier voxel set, purple color : final throat barrier voxel set

A B
C D
E

Figure 12 : Sphere pack sample with 5 representative throats. Blue color area indicate one of five representative throats. This sample has 3 layer spheres and each layer has 9 spheres.

A B
C

Figure 13 : Throat area (A), Pore volume distribution (B), and Coordination number (C) computed from the XCMT image of the Handford wet inlet t241 samples using 5 algorithms

 A B C D

Figure 14 : A,B : Example of calculating the exact linear distance (A) and circumferences (B) using the point set algorithm. C,D : Relative error in the digitized image of the real boundary (a) line with different angel. (b) circle with different radii.

A B
C D
E

Figure 15 : Counter example to find throat area using the sum of triangular area (in 3DMA thorat area calculation algorithm). Existing algorithm can not calculate the exact throat area of this surface. A : blue colored curve shows a perimeter of the surface. 3DMA throat algorithm cannot calculate the exact throat in orange colored region. C, D : different point of view of B. E : desired throat surface.

## References

• [1] W. B. Lindquist, A. Venkatarangan, J. Dunsmuir, T. f. Wong. Pore and throat size distributions measured from synchrotron X-ray tomographic images of Fontainebleau sandstones. J Geophys Res 2005;105:21508-28.
• [2] W. B. Lindquist. 3DMA-Rock General Users Manual, Department of Applied Mathematics and Statistics State University of New York at Sony Brook, 1999.
• [3] H. Shin, W. B. Lindquist, D. L. Sahagian, S. R. Song. Analysis of the vesicular structure of basalts. Comput Geosci 2005; 31:473-87.
• [4] J.-W. Kim, D. Kim, W. B. Lindquist. A Re-examination of Throats. Water resources research 2013;7615-7626.
• [5] E. W. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik 1959; 1, 269-271.
• [6] T. C. Lee, R. L. Kashyap, C. N. Chu. Building skeleton models via 3-D medial surface / axis thining algorithms. CVGIP: Graph Models Image Proc 1994;56:462-78.
• [7] W. Oh, W. B. Lindquist. Image thresholding by indicator kriging. IEEE Trans Pattern Anal Mach Intell 1999;21:590-602.