Introduction

A porous medium is essentially treated as a physical structure divided into a void space, referred as pores, and a solid matrix, named grains. In many applications of porous media, the main interest is sometimes related to the pores, particularly when the research focus is on fluid flow and mass transport, while in others situations the focus is on the grains, especially when the structural behavior under load is the relevant aspect of study. These occurrences are observed, for instance, in porous media like reservoir rocks, soil or trabecular bone.

Currently, most of the research on porous media microscale is based on imaging, particularly microCT (\(\mu\)CT) imaging of samples as shown in the following papers Sidorenko et al (2021), Sakugawa et al (2020), Irayani et al (2019), Andrä et al (2013a), Andrä et al (2013b). In addition, for simplicity, the 2D or 3D grayscale image is converted to binary, with the black pixels or voxels representing the pores and the white ones the grains, such that the color attribute of the pores is 0 and that of the grains is 255 (or 1). Notice that sometimes it is convenient to work with the complement by inverting it, turning white (255) to represent the pores and black (0) the grains.

When studying the fluid flow through a porous medium, the effective pores define the route through where the fluid may flow. In other words, the connected pores that have an inlet at one face and an outlet at another face are the ones that allow the fluid to move through the region or volume of interest (ROI or VOI). There are a couple of properties associated with the pores that have to be considered, such as the pore thickness, the pore connectivity and their degree of tortuosity, for instance. Additionally, it is also important to know what is the preferential direction of the pores in the ROI/VOI, what percentage of the pores are aligned horizontally, vertically or diagonally. The pore space main directions influence how the fluid will naturally or forced flow under the action of pressure or gravity, for instance.

The trabecular bone is a spongy-like structure, with the trabeculae as the solid matrix (grains) and the marrow filling the pore space. The predominant direction of the trabeculae is an important characteristic of the structure as it acts as load-bearing paths to spread out applied stress Roque and Alberich-Bayarri (2015). In soil physics, permeability and consolidation are the most variable soil properties. These parameters may vary with depth even in case of homogeneous soil layers. However, due to the anisotropy, the coefficient of consolidation and the coefficient of permeability in the horizontal direction are typically different from their values in the vertical direction Laskar and Pal (2018). The vertical and horizontal pore directions play a relevant role for water flow and nutrient transport to different ground depths as well as for the water retention curve.

In the field of carbonate reservoir rocks, porosity and permeability are quite important rock properties; they are, respectively, essentially responsible for the storage capacity and flow of hydrocarbons in the reservoir. Additionally, the pore connectivity, tortuosity and pore predominant directions are determinant for the fluid and electrical current to flow through the rock Tiab and Donaldson (2016).

The fluid flow in a porous medium depends on the pore space. On the other hand, the pore network influences the directional permeability, and thus, the predominant direction of the pore network constrains the fluid flow. For instance, air permeability is an important property of fabrics as it can be a decisive parameter for the functionality of the fabric. It has been shown that the pore direction in woven fabrics affect the air permeability Havlov (2014).

The special core analysis laboratory (SCAL) methods used to estimate porous medium quantities are costly and time-consuming, but more recently digital rock analysis Kalam (2012) has significantly minimized these downside effects. Digital rock physics (DRP) Berg et al (2017), Al-Marzouqi (2018) is a breakthrough technology to digital oilfield that merges high-resolution imaging, image processing techniques and numerical simulations on the pore scale based on porous medium core samples.

In this paper, it will be presented and discussed a technique to identify the predominant pore direction in a porous medium and through it do a partition of the pore space according to the pore orientation. The layout of the paper is as follows: First, Sect. 2 presents the theoretical framework with the 2D (2.1) and 3D (2.2) general equations defining the predominant pore directions and providing the partition of the pore space into vertical, horizontal and diagonal pores; then, Sect. 3 presents the results of the application of the techniques to 2D and 3D different case studies. Section 4 provides a brief discussion on the necessary calculations and computational costs to compute the 2D and 3D pore predominant directions. Finally, Sect. 5 concludes the paper.

Theoretical framework

For the purpose of the techniques that will be presented here, it will be considered 2D and 3D binary images of the porous medium. Let us consider initially a 2D binary image I of a porous medium and assume that I is an image of size \(n \times m\) pixels that a pixel \(p\in I\) is represented by a square and has coordinates \(p =(x, y)\) set at the center of the pixel, following the standard reference frame adopted in image processing, where the origin pixel (0, 0) is at the top left corner of the image. Here, it is adopted that in a binary image of a porous medium, the solid matrix (grains) is represented by black pixels with a 0 color attribute and the voids (pores) are white pixels with a 255 color attribute. Likewise, 3D binary images will have the same attributes as for the pixels, with the voxel coordinates set at the center of the voxel.

2D directional partition of pores

In this, it will be introduced the foundations for the identification of a 2D pore orientation in a porous medium. A pore will be classified according to their horizontal (H), vertical (V) or diagonal (D) predominant direction, and the pore (or grain) space will be partitioned into a subset of pores according to their predominant orientation. For that, given a single pore pixel, say \(p_0=(x_0,y_0) \in I\), the problem consists in how to identify what is its predominant direction. Figure 1 gives an illustration of the three major pore predominant directions defined in this study.

Fig. 1
figure 1

2D horizontal, vertical and diagonal directional partitions

Full size image

To determine which is the predominant direction of a pore, according to the partitions of the plane given in Fig. 1, the star length distribution principle Smit et al (1997) will be applied from a pore pixel \(p_0 = (x_0, y_0)\) with \(k=16\) rays (straight lines) being launched in every \(\theta = \pi /6\) interval, performing a total of \(N=192\) rays. According to the star length distribution principle, to identify the pore pixel predominant direction we look for the first grain (black pixel) to be intercepted by the ray. If no pixel is found, we take the pore pixel that has reached the boundary of the image \(I_{n\times m}\). Then, the Euclidean distance along the ray \(L(\theta _i), i =0, \ldots , 191\) is computed from the seed pixel \(p_0\) to the first grain one or to the boundary, as Fig. 2 illustrates.

Fig. 2
figure 2

Schematic illustration of the star length distribution on 2D pore space

Full size image

To obtain the 2D pore predominant direction, the following expressions are calculated for each partition:

Horizontal partition \({\mathcal {H}}\): \(0 \le \theta < \pi /6\), \(5\pi /6 \le \theta < 7\pi /6\) and \(11\pi /6 \le \theta < 2\pi\).

$$\begin{aligned} H(p_0) = \sum _{i=0}^{\frac{N}{12}-1} L(\frac{2\pi i}{N}) + \sum _{i={5N/12}}^{\frac{7N}{12}-1} L(\frac{2\pi i}{N}) + \sum _{i={11N/12}}^{N-1} L(\frac{2\pi i}{N}) \end{aligned}$$
(1)

Vertical partition \({\mathcal {V}}\): \(\pi /3 \le \theta < 2\pi /3\) and \(4\pi /3 \le \theta < 5\pi /3\).

$$\begin{aligned} V(p_0)= \sum _{i={N/6}}^{\frac{N}{3}-1} L(\frac{2\pi i}{N}) + \sum _{i={2N/3}}^{\frac{5N}{6}-1} L(\frac{2\pi i}{N}) \end{aligned}$$
(2)

Diagonal partition \({\mathcal {D}}\): \(\pi /6 \le \theta < \pi /3\), \(7\pi /6 \le \theta < 4\pi /3\) for \(D^+\) and \(2\pi /3 \le \theta < 5\pi /6\), \(5\pi /3 \le \theta < 11\pi /6\) for \(D^-\).

$$\begin{aligned} D(p_0) = \sum _{i={N/12}}^{\frac{N}{6}-1} L(\frac{2\pi i}{N}) + \sum _{i={7N/12}}^{\frac{2N}{3}-1} L(\frac{2\pi i}{N}) + \sum _{i={N/3}}^{\frac{5N}{12}-1} L(\frac{2\pi i}{N}) + \sum _{i={5N/6}}^{\frac{11N}{12}-1} L(\frac{2\pi i}{N}) \end{aligned}$$
(3)

where \(\theta _i = \frac{2\pi i}{N}\) corresponds to the angle of each ray launched from the seed pore pixel \(p_0\). The predominant direction of \(p_0\) is obtained from

$$\begin{aligned} Dir(p_0) = \max \{H(p_0), V(p_0), D(p_0)\}, \end{aligned}$$
(4)

that means the direction with the largest value defines the predominant pore direction (PPD) of \(p_0\).

Remind that \({\mathcal {H}}\), \({\mathcal {V}}\) and \({\mathcal {D}}\) are mathematical partitions such that they satisfy \({\mathcal {H}} \cap {\mathcal {V}} \cap {\mathcal {D}} = \emptyset\) and \({\mathcal {H}} \cup {\mathcal {V}} \cup {\mathcal {D}} = PS\), where PS is the pore space of the sample.

3D directional partition of the pores

The 3D pore direction partitions are a natural extension of the 2D case. Figure 3 illustrates the 3D directional partitions. Let us consider the usual image processing spherical coordinates and that from a pore voxel \(v_0\) a set of rays are launched according to \((r, \theta , \phi )\), where r are the rays departing from \(v_0\) to the nearest grain or to the boundary of the VOI, \(\theta\) is the polar angle, and \(\phi\) is the azimuthal angle. The Euclidean distance along a ray will be given by \(L(\theta _i, \phi _j)\), \(0 \le \theta _i \le \pi\) and \(0 \le \phi _j < 2\pi\).

Fig. 3
figure 3

3D horizontal (red), vertical (blue) and diagonal (green) directional partitions

Full size image

To identify the pore voxel predominant direction, let us consider that 16 rays are launched from a pore seed voxel \(v_0\) within each interval \(\theta = \pi /6\) and \(\phi = \pi /6\), leading to \(N_{\theta }= 24\) polar rays and \(N_{\phi }=48\) azimuthal rays, providing a total of \(N = 1152\) radial rays from each pore voxel. The Euclidean distance along a ray will be given by \(L(\theta _i, \phi _j)\), \(0 \le \theta _i \le \pi , i=1, \ldots , 24\) and \(0 \le \phi _j < 2\pi , j=1, \ldots , 48\). Notice that now the total number of rays is 6 times greater than for the 2D case, which causes an increase in the computational cost associated with the size of the problem, in other words, to the sample volume, its porosity and mathematical operations.

The 3D directional partitions are given by:

Horizontal partition \({\mathcal {H}}\): \(\pi /3 \le \theta < 2\pi /3\) and \(0 \le \phi < 2\pi\).

$$\begin{aligned} H(v_0)= \sum _{i=N_\theta /3}^{\frac{2N_\theta }{3}-1}\; \sum _{j=0}^{N_\phi - 1} L(\frac{\pi i}{N_\theta },\frac{2\pi j}{N_\phi }) \end{aligned}$$
(5)

Vertical partition \({\mathcal {V}}\): \(0 \le \theta < \pi /6\), \(5\pi /6 \le \theta \le \pi \;\) and \(0 \le \phi < 2\pi\).

$$\begin{aligned} V(v_0)= \sum _{i=0}^{\frac{N_\theta }{6}-1}\; \sum _{j=0}^{N_\phi -1} L(\frac{\pi i}{N_\theta },\frac{2\pi j}{N_\phi }) + \sum _{i=5N_\theta /6}^{N_\theta -1}\; \sum _{j=0}^{N_\phi -1} L(\frac{\pi i}{N_\theta },\frac{2\pi j}{N_\phi }) \end{aligned}$$
(6)

Diagonal partition \({\mathcal {D}}\): \(\pi /6 \le \theta < \pi /3\), \(2\pi /3 \le \theta < 5\pi /6\) and \(0 \le \phi < 2\pi\).

$$\begin{aligned} D(v_0)= \sum _{i=N_\theta /6}^{\frac{N_\theta }{3}-1}\; \sum _{j=0}^{N_\phi -1} L(\frac{\pi i}{N_\theta },\frac{2\pi j}{N_\phi }) + \sum _{i=2N_\theta /3}^{\frac{5N_\theta }{6}-1}\; \sum _{j=0}^{N_\phi -1} L(\frac{\pi i}{N_\theta },\frac{2\pi j}{N_\phi }) \end{aligned}$$
(7)

The predominant direction of a pore voxel \(v_0\) is obtained from

$$\begin{aligned} Dir(v_0) = \max \{H(v_0), V(v_0), D(v_0)\}; \end{aligned}$$
(8)

in other words, the direction with the largest value defines the predominant pore direction (PPD) of voxel \(v_0\).

PPD applications

In this section, application tests of the PPD identification and directional partitioning of the pores are presented. First, the PPDs for a 2D synthetic binary image are evaluated and then for a 2D binary image of a Berea sandstone rock sample. These are two simple cases just to illustrate the feasibility of the technique.

For the 3D case, three application tests are considered. First, a Berea sandstone and a Carbonate reservoir rock core samples are investigated, both were obtained from the MicroCT Images and Networks of Imperial College London database and, finally, an application considering the grain partitioning of a distal-radius trabecular bone structure. These are well-known type of porous medium structures; in particular, the rock core samples are typical of a sandstone and a carbonate reservoir.

The 2D and 3D algorithms for directional partitioning of a pore space were implemented in an in-house software using Phyton running on a Windows\(^\circledR\) OS. In addition, the computation of the number of pixels (2D) or voxels (3D) in the sample, the pore space porosity, the partition porosities, the effective pore networks and their porosities is also provided by the software.

2D synthetic and Berea image samples

As a simple and effective test for the PPD technique, a 2D synthetic binary image was created with 6 different geometrical pictures that include straight horizontal, vertical and diagonal lines, a curved line, two solid rectangles and two solid ellipses. Figure 4 shows, on the left, a 2D synthetic binary image and on the right the identification of its PPD partitions: horizontal in red, vertical in blue and diagonal in green. The straight lines show correctness of the 2D PPD. As the curved line is a digital curve, that is, a discrete line formed by pixels, the 2D PPD H, V and D partitions are present along the same line. The same occurs for the solid pictures. Remind that these are digital pictures and the 2D PPD is defined for an angular interval.

Fig. 4
figure 4

2D synthetic binary image (left) and its PPD partitions with H in red, V in blue and D in green (right)

Full size image

Figure 5 shows the application of the 2D PPD equations given in 2.1 for a 2D binary image of a Berea sandstone sample seen on the left and its PPD partitions on the right.

Fig. 5
figure 5

2D Berea Sandstone binary image (left) and the PPD partitions with H in red, V in blue, D in green (right)

Full size image

The 2D PPD technique can be applied for thin sections Saxena and Mavko (2016); Saxena et al (2017), and a 3D reconstruction of a sample can be made; nevertheless, this approach is not accurate as it does not take into account the rays in the z-direction as indicated in Sect. 2.2.

3D PPD applications

Let us consider the partitioning based on the 3D PPD Eqs. 5, 6, 7 and 8 applied to the Berea sandstone and Carbonate C1 reservoir rock core samples to show their pore predominant direction subspaces, and to the distal-radius trabecular bone to show its grain predominant direction subspaces.

Berea sandstone

The first sample is a Berea sandstone that has volume size 400 \(\times\) 400 \(\times\) 400 voxels and 5.345 \(\mu\)m resolution. Figure 6 shows the 3D full pore space in grayscale; red, blue and green are the horizontal, vertical and diagonal PPD subspaces, respectively, and finally the full pore space exhibiting the union of the three PPD subspace partitions. Just the horizontal, vertical and diagonal PPD partitions subspaces are displayed in Fig. 7.

Fig. 6
figure 6

3D Berea pore space (top left), its H (red), V (blue) and D (green) PPD partitions and full pore space exhibiting the union of the PPD partitions (bottom right)

Full size image
Fig. 7
figure 7

Individual view of Berea PPD partition subspaces, H (cyan), V (yellow) and D (magenta)

Full size image

Figure 8 shows the effective pore networks of the full sample and the H’, V’ and D’ PPD partitions obtained from the effective pore network.

Fig. 8
figure 8

3D effective pore networks for the Berea sandstone sample and the H’, V’ and D’ PPDs partitions from the effective pore network. Different colors indicate distinct partitions

Full size image

Notice that the H’, V’ and D’ PPDs obtained from the full effective pore network shown in Fig. 8 are not necessarily effective pores in the respective directions. Figure 9 shows the effective pore network from each one of these directional partitions. The vertical PPD partition shown in Fig. 8 does not have any effective pore network.

Fig. 9
figure 9

H (left) and D (right) effective directional pore networks for the Berea sandstone. Different colors indicate distinct effective pore networks

Full size image

Table 1 gives the Berea sandstone sample porosity \(\phi\), the PPD porosities \(\phi _{\{H,V,D\}}\), the effective pore space porosity \(\phi _e\), the PPD porosities \(\phi _{\{H', V', D'\}}\) obtained from the effective pore space and PPD porosities \(\phi _{\{eH,eV,eD\}}\) of the effective directional pore networks. Notice that 99% of the pore space are effective pores. From the H’, V’ and D’ partitions, only 67%, 0% and 1% of the pores form the fundamental directional effective networks.

Table 1 3D Berea sandstone sample porosity \(\phi\), PPD porosities \(\phi _{\{H,V,D\}}\), effective porosity \(\phi _e\), PPD porosities \(\phi _{\{H', V', D'\}}\) obtained from the effective pore space and PPD porosities \(\phi _{\{eH,eV,eD\}}\) of the effective directional pore networks
Full size table

Carbonate C1

The second application is for the Carbonate C1 sample from the database. Its volume size is 400 \(\times\) 400 \(\times\) 400 voxels and 2.85 \(\mu\)m resolution. Figure 10 shows the 3D full pore space, the horizontal, vertical and diagonal PPD pore subspaces and finally the full pore space exhibiting the union of the three pore subspace partitions.

Fig. 10
figure 10

3D Carbonate C1 pore space (top left), its H (red), V (blue) and D (green) PPD partitions and full pore space exhibiting the union of the PPD partitions (bottom right)

Full size image

Figure 11 shows the effective pore network of the full sample and the H’, V’ and D’ PPD obtained from the effective pore network. Figure 12 shows the effective pore networks from each one of these partitions.

Fig. 11
figure 11

3D effective pore networks for the carbonate C1 sample and the H’, V’ and D’ PPDs partitions from the effective pore networks. Different colors indicate distinct partitions

Full size image
Fig. 12
figure 12

From left to right, H, V and D effective directional pore networks for Carbonate C1. Different colors indicate distinct effective pore networks

Full size image

Table 2 gives the Carbonate C1 sample porosity \(\phi\), the PPD porosities \(\phi _{\{H,V,D\}}\), the effective pore space porosity \(\phi _e\), the PPD porosities \(\phi _{\{H', V', D'\}}\) obtained from the effective pore space and PPD porosities \(\phi _{\{eH,eV,eD\}}\) of the effective directional pore networks. Notice that 92% of the pore space are effective pores. From the H’, V’ and D’ partition, only 89%, 34% and 32% of the pores form the fundamental directional effective networks.

Table 2 3D Carbonate C1 sample porosity \(\phi\), PPD porosities \(\phi _{\{H,V,D\}}\), effective porosity \(\phi _e\), PPD porosities \(\phi _{\{H', V', D'\}}\) obtained from the effective pore space and PPD porosities \(\phi _{\{eH,eV,eD\}}\) of the effective directional pore networks
Full size table

Trabecular bone

Let us consider now a 3D distal-radius trabecular bone image obtained by \(\mu\)CT from an ex vivo sample, which has a volume of 257 \(\times\) 257 \(\times\) 240 voxels and resolution of 34 \(\mu\)m. This sample corresponds to an osteoporotic trabecular bone structure; additional details about this sample, named 269, can be found in Roque et al. (2013).

Figure 13 shows the 3D full grain space of the distal-radius trabecular bone, the horizontal, vertical and diagonal predominant grain direction (PGD) partitions and the full grain space exhibiting the union of the three partitions.

Fig. 13
figure 13

3D full grain space (top left) of the distal-radius trabecular bone sample, H (red), V (blue) and D (green) PGD and the full grain space exhibiting the union of the three subspace partitions

Full size image

Figure 14 shows the effective grain network of the sample and its H’, V’ and D’ PGD partitions. Figure 15 shows the effective fundamental directional grain network from each one of these partitions.

Fig. 14
figure 14

3D effective grain space of the trabecular bone sample and the H’, V’ and D’ PGPs from the effective grain space. Different colors indicate distinct effective partitions

Full size image
Fig. 15
figure 15

Effective fundamental directional TB networks. Different colors indicate distinct effective grain networks

Full size image

Table 3 gives the trabecular bone volume fraction (\(\varphi = \frac{BV}{TV}\), BV is bone volume and TV is total volume of the sample), the PGD TBVF \(\varphi _{\{H,V,D\}}\), the effective TBVF \(\varphi _e\), the PGD TBFV \(\varphi _{\{H', V', D'\}}\) obtained from the effective TBFV and PGD TBFV \(\varphi _{\{eH,eV,eD\}}\) of the effective directional trabecular networks. Notice that 98% of the trabecular bone space are effective trabeculae. From the H’, V’ and D’ partitions, only 10%, 59% and 2% of the trabeculae form the fundamental directional effective networks.

Table 3 3D distal-radius sample TBVF \(\varphi\), PGD TBFV \(\varphi _{\{H,V,D\}}\), effective TBVF \(\varphi _e\), PGD TBVF \(\varphi _{\{H', V', D'\}}\) obtained from the effective TB space and PGD TBVF \(\varphi _{\{eH,eV,eD\}}\) of the effective directional TB networks
Full size table

In this section, a set of two 2D and three 3D application tests of the PPD have shown how the technique can be used to identify the pore orientation on 2D or 3D microCT images, respectively, and partition the pore space of the sample (reservoir rocks for pore and trabecular bone for grain) into horizontal, vertical and diagonal directional subspaces. In addition, the effective pore networks were identified and their PPDs partitions were obtained with their porosities calculated and the results are summarized in Tables 1, 2 and 3.

Discussion

The predominant direction of a pore in a 2D or 3D space is dependent on the reference frame orientation and on the length of the rays. The image processing frame follows the right-hand rule, as such a pore direction changes due to frame rotations, and a change in the pore structure itself may change the length of the rays, which provoke orientation changes as well. A pore orientation is not affected by a simple change in the sample porosity as long as the pore structure where it belongs to is not modified. The PPD partitions (see Figs. 1 and 3) can be easily adjustable, leading to narrower or wider pore orientations.

The effective pore networks of a sample show the set of pores that are connected to each other from an inlet face to a distinct outlet face. They form the pore network that allows a fluid to flow throughout the sample. The H, V and D effective fundamental directional pore networks give an indication on how the pores get directionally connected to each other. This induces a natural direction choice for fluid flow under the action of a strength. From Table 1, it can be seen that the effective horizontal porosity of Berea sandstone is higher than the vertical and diagonal ones, which gives an indication that individually there will be more facility to horizontal flow than in the other two directions. Similarly, in the Carbonate C1 case (Table 2) the H effective fundamental directional porosity is higher than the other two directions, which indicates that a fluid would flow more naturally in the horizontal direction instead.

In the last case, if we look at the trabecular structure, the vertical effective fundamental directional trabecular bone volume fraction is higher than in the other two directions. According to previous studies Roque et al. (2013), this result was actually expected as the trabeculae are primarily aligned along the direction where they suffer more frequent stress; in such sample, it corresponds to the vertical direction (z-direction, distal direction) to provide more strength to the bone structure.

In Figs. 9 and 12, it can be seen that the effective fundamental directional pore network and in 15 the effective fundamental directional trabecular bone network were obtained from the respective H’, V’ and D’ partitions

The H’, V’ and D’ partitions shown in Figs. 8, 11 and 14 were obtained from the effective pore (or grain) network of the respective samples. However, in the process of identifying these partitions some of the effective pore (or grain) network loses connectivity between previous connected pores that had distinct predominant directions, turning to directional partitions where the pores (or grains) are no longer effective in the corresponding direction. That is why the effective fundamental directional pore (or grain) networks shown in Figs. 9, 11 and 15 have a much less porosity (or TBVF) as shown in Tables 1, 2 and 3.

Algorithm implementation and computational cost

The main difficulty to implement the algorithms to compute the 2D and 3D PPDs is to find out the first black (grain) pixel/voxel that is intercepted by the ray or, otherwise, the white (pore) pixel/voxel that are at the boundary of the region/volume of the sample.

Suppose a 2D image of size \( n\times m\) with \(P_{\scriptscriptstyle 2D}\) pore pixels. From each pore pixel, a total of 192 rays is launched (see Sect. 2.1), and thus, there will be \(192 \, \phi _{\scriptscriptstyle 2D}\, n\, m\) rays launched for the sample, where \(\phi _{\scriptscriptstyle 2D}\) is the 2D sample porosity, and for each one of the pores, there is the need to find out the first black pixel (grain) that is intercepted by the ray or the pixel at the boundary along the ray direction. For the example of 2D Berea sandstone given in Fig. 5, whose size is \(400 \times 400\) pixels and porosity \(\phi _{\scriptscriptstyle 2D} = 0.2012\), the cost is \(\approx 6.18 \times 10^6\) searches.

For the 3D image of size \( n\times m \times k\), with \(P_{\scriptscriptstyle 3D}\) pore voxels, a total of 1152 rays are launched from each pore voxel (see Sect. 2.2); therefore, there will be \(1152\, \phi _{\scriptscriptstyle 3D}\, n\, m\, k\) rays launched for the volume sample, where \(\phi _{\scriptscriptstyle 3D}\) is the 3D sample porosity, and for each one of the pores, there is a need to find out the first grain voxel that is intercepted by the ray or the voxel at the boundary along the ray direction.

If we just consider the total number of rays launched from each pore voxel as the size of the problem, the order of magnitude of the problem from 2D to 3D grows \(6k \frac{\phi _{\scriptscriptstyle 3D}}{\phi _{\scriptscriptstyle 2D}}\) times. The 2D search algorithm for the first black pixel intercepted by the ray is relatively simple; however, for the 3D case, it is more complex, which turns the computational cost much higher yet, demanding more computational resources for large 3D samples. The computational cost may impose a certain limitation for the PPD determination when the size of the samples and/or when the number of rays launched from the source pore pixel (2D) or pore voxel (3D) is increased.

The PPD and effective pore network identification algorithms for the 2D and 3D cases were implemented in an in-house software using Phyton. The applications presented in Sect. 3.2 were run in a PC with i7-8565U CPU, 1.80GHz, quad-core and 16Gb of RAM, running under Windows operating system. Table 4 gives the number of rays launched for each sample, the CPU time spent to run the 3D H, V, D PPD for the Berea sandstone and Carbonate C1 samples and the PGD for the distal-radius trabecular bone sample. It is also given the CPU time spent to find the effective pore (EF) or grain (EG) networks for these samples.

Table 4 Number of rays launched, CPU time spent by the computer to run the PPD/PGD and EF/EG for the samples
Full size table

Summary and conclusion

In summary, in this paper we have presented and discussed:

  • A technique and its computational implementation to identify the predominant pore (or grain) direction (PPD) as horizontal, vertical or diagonal, in a porous medium.

  • The partitioning of a 2D or 3D pore (or grain) space of a porous medium according to the PPD.

  • The algorithms for 2D and 3D PPD partitioning were implemented in Phyton, and a set of application tests were done: First, it was considered a single 2D synthetic image which included vertical, horizontal and diagonal digital straight lines, a curved line and solid pictures, and then a Berea sandstone sample. These have shown the accuracy of the technique. Secondly, a 3D Berea sandstone and Carbonate reservoir rock core image samples, that are available at the MicroCT Images and Networks of Imperial College London database, were considered as they are very typical of reservoir rocks. In addition, a 3D application was done considering the grain partitioning of an osteoporotic distal-radius trabecular bone structure.

  • The effective pore network of the 3D pore space was obtained, and the PPD and partition of the effective pore network were computed for each sample.

  • The 2D and 3D pore space porosities effective porosities, H, V and D PPD subspace porosities and their respective PPD effective pore network porosities were also obtained for the reservoir rocks and, similarly, for the grain volume fraction of the distal-radius trabecular bone, whose results can be found in Tables 1, 2 and 3.

The technique presented in this paper is novel providing an additional tool to improve the digital analysis of porous medium properties with respect to their pore orientations, as such it needs to be further explored. Currently, work is in progress to investigate how the PPD influences the pore network connectivity, its tortuosity and the permeability of reservoir rock core samples.