Methodology for obtaining contact angles in rock sample images using image processing and polynomial fitting techniques

Estimating the contact angle in very complex rock pores presents some challenges to accurately identify the fluid–rock contact surface. This work presents a methodology to estimate the contact angle formed by the brine–rock and the brine–oil interfaces on processing high-resolution images provided by micro-CT scan. We focus the discussion on the limitations of the most popular computational techniques used to determine the contact angle and discuss how to select a practical way to evaluate it. The method consists of four steps: (1) processing the images to determine each fluid present in the image, (2) selection of the pixels that will be part of the contact interface of fluids and the contact point, (3) fitting polynomial equations for each interface and selection of the equation that gives the lowest error, (4) estimation of the contact angle based on the more appropriate polynomial equation. The contact angle is calculated based on the slope of the interfaces’ tangents at the contact point. Several types of approaches were tested to determine the contact interface and the contact point. In order to evaluate the applicability of our method, we use an analytically generated image and rock sample images. Potential errors between the angle obtained from the analytically generated image and the angle calculated from the method show the impact of the right selection of pixels during the image processing step. High sensitivity is also observed for the tangent values in the presence or absence of pixels from the rock sample analysis.


Introduction
The depletion of reservoirs evokes the necessity of applying enhanced oil recovery (EOR) methods, in order to prolong the field life cycle and increase the recovery factor as much as possible. As fewer oil discoveries are made, EOR becomes an option to maintain or increase the production of mature fields, playing a more significant role in the global supply of oil and gas. An important key to implement an EOR system is the determination of the rock wettability. Wettability, by definition, is the tendency of a fluid to adhere to a solid surface in the presence of another immiscible fluid (Zisman 1964) and can be quantified, at the pore scale, by the contact angle formed at the interface of the two fluids in balance with a solid surface. Small contact angles (0 ≤ θ ≪ 90°) represent a higher wettability, i.e., the fluid spread over a large area of the surface, while large contact angles (90° ≪ θ ≤ 180°) imply a lower wettability, i.e., the fluid forms a compact liquid droplet and minimizes its contact with the surface (Bracco and Holst 2013).
Calculations based on measured contact angle are fundamental in determining rock wettability characteristics, such as if a rock is oil-wet or water-wet. The pore-scale topology, interfacial tension, and contact angle, in turn, control the parameters for the multiphase flow in porous media such as capillary pressure, fluid saturation, and relative permeability (Gray and Miller 2011). Thus, the estimation of the contact angle helps to model multiphase flow in order to enhance oil and gas recovery (Pope and Braviere 1991) and geological CO 2 storage (Chalbaud et al. 2009;DePaolo and Cole 2013). One way to estimate the contact angle is through the use of high-resolution images of rock samples.
Recent advancements in high-resolution imaging techniques such as micro-CT scan provide an opportunity to describe heterogeneous characteristics such as mineral composition, geometry, and surface roughness (Morrow 1975;Wan et al. 2014), and in situ fluid distributions in pore structures under varying conditions (Andrew et al. 2013;Wildenschild and Sheppard 2013). Therefore, images need to go through image processing techniques in order to segment the image while preserving the edges (Andrew et al. 2014;Scanziani et al. 2017). Other studies start the analysis with already segmented images Klise et al. 2016), relying on the accuracy of previously performed image processing.
The most employed techniques to measure the contact angle are still performed using manual processes, such as using a goniometer after image processing (Andrew et al. 2014). However, it is hard for an inexperienced person to obtain the contact angle from the images, as it is difficult to define the interfaces and the contact point due to the images being blurry and the edges not being clear. Thus, image processing becomes an important tool during the process of obtaining a good segmentation, which in turn will provide a better estimate of the contact angle. Recently, several algorithms have been developed to automatically measure the contact angle Klise et al. 2016;Scanziani et al. 2017). One of the problems encountered in the automatic measurement is the correct selection of an interface between the two immiscible fluids and the solid phase. For a drop profile, Bateni et al. (2003) develop an automated polynomial fitting (APF) scheme to describe the interfaces. Scanziani et al. (2017) propose an approach that is based on the real-life insight that the fluid-fluid interface has a constant curvature. Alratrout et al. (2017) extract a surface mesh from a segmented voxelized image. Automated methods applied to estimate the contact angle through image processing are techniques relatively new and unexplored in the literature, and this work aims to develop further these methodologies. Thus, if such methods are applied correctly, they will provide a better understanding of the reservoir and benefit-enhanced oil recovery techniques.
The present work focuses on discussing the limitations of the most popular computational techniques used to determine the contact angle, and on how to select a practical way to evaluate it. Several types of approaches were tested in order to determine the pixels that belong to each interface of the contact point, as well as to analyze the range of the polynomial order employed. We use an analytically generated image as a reference to evaluate the applicability of our methodology, and estimation errors of the contact angle are obtained. The method and range of polynomial order with the best performance, i.e., with lowest overall error, are then chosen to be employed in slices of a carbonate reservoir sample generated by Alhammadi et al. (2017). The dataset is composed of high-resolution X-ray micro-tomography at subsurface conditions.
In the next sections, methods, materials, results, and conclusions are contemplated. The methods describe the algorithm employed to compute the contact angle. The materials show the images employed for the analyses. The results present and discuss the data. Finally, the principal findings are summarized in the conclusions.

Methods
In this section, the steps of the algorithm applied to calculate the contact angle are described. The steps are (1) image processing, (2) pixels selection, (3) polynomial fitting, and (4) estimation of the contact angle.

Image processing
In image processing, the coordinates of the contact point and the interfaces between brine-oil and brine-rock are identified using the open-source library scikit-image for the filters. Initially, Laplacian of Gaussian edge detection (Marr and Hildreth 1980) is employed to remove noise. Then, the image is segmented through multi-Otsu thresholding (Otsu 1979) into three categories: rock (solid phase), oil (light liquid phase), and brine (dense liquid phase). For edge detection, the Sobel filter (Sobel and Feldman 1968) is applied generating a binary image. In order to avoid removing any relevant pixel during the image processing, pixels with gradient magnitude greater than zero receive the value 1. In addition, pixels belonging to corners are removed to avoid possible problems during computation.
Since the binary image has only edges from the segmented image, its pixels can be defined as brine-oil interface, as brine-rock interface, or as three-phase contact points. The pixel selection for the interfaces or for the contact point is performed by firstly mapping only pixels with value 1 from the binary image, and then by comparing them with the segmented image. The mapping process is performed using a 3 × 3 matrix, in which each matrix element represents a pixel. The pixel with value 1 from the binary image is set in the center of the matrix. If only brine and oil are in the matrix, according to pixel information from the segmented image, the pixel that was set in the center of the matrix is defined as a pixel from the oil-brine interface; meanwhile, if the matrix has only brine and rock pixels, the center pixel is defined as belonging to the rock-brine interface. If oil, rock, and brine are in the same matrix, the corresponding pixel is defined as a contact point. Finally, for a pixel to be a contact point, its category in the segmented image must be a brine pixel. Figure 1 summarizes the classification procedure of pixels during the mapping process. Each image represents a 3x3 matrix, in which the cell with an "X" mark represents the pixel to be classified, and the colors black, light gray, and white represent the oil, the brine, and the rock, respectively. According to the procedure described above, in Fig. 1a, the pixel at the center of the matrix is classified as a brine-oil interface pixel, in Fig. 1b, the pixel is a brine-rock interface, and in Fig. 1c, it is classified as a contact point.

Pixel selection
Initially, the contact point is defined, from which the calculations are based. Starting with mapping the brine-rock interface, its pixels are labeled to ensure the connectivity of the elements with eight-connected neighborhood. Pixels possessing a different label than the contact point are removed from the pixel selection. Thereafter, in order to select a limited number of pixels, a region of interest centered at the contact point is defined according to a constant value. The same mapping process is performed for the brine-oil interface. The following steps present the differences among six methods chosen for the pixel selection.
From the segmented image, methods 1, 2, and 3 take into account only the brine pixels for the interfaces, while methods 4, 5, and 6 consider only oil and rock pixels for the brine/oil and brine/rock interfaces. The difference between 1 and 4, 2 and 5, and 3 and 6 lies only in the selection of the pixel category. Thus, from now on, the methods are referred to as 1/4, 2/5, and 3/6.
Since the applied threshold, which is performed after the Sobel filter, intends to make the edge thick, the methods are employed to remove the excess of pixels. The method 1/4 computes the contour of the interface, keeping only the pixels belonging to the contour and to their respective pixel category; the method 2/5 solely selects the pixels according to the category of the method; and in the method 3/6, a skeletonize algorithm is applied after selecting the pixel category. In some methods, the contact point is removed from the interfaces, and thus it must be re-added later. Table 1 summarizes these methods. Figure 2 shows the difference between methods for the same image. The light blue curve represents the brine-oil interface, the orange line characterizes the brine-rock interface, and the black dot is the contact point. Although the distinction of the pixel selection is subtle, as will be observed in the next section, the contact angles results can diverge greatly.

Polynomial fitting
The polynomial fitting is based on the approach presented by Bateni et al. (2003). The pixel coordinate from brine-oil and brine-rock interfaces, obtained in the pixel selection, is employed to create approximation curves used to describe each interface as a polynomial. The range of the polynomial fitting order is defined to vary from first to sixth order for each interface, storing only the polynomial coefficients that give the lowest error. Thus, the polynomial fitting has the flexibility to fit each interface to its best approximation curve, without the constraint of a fixed polynomial order. Equation 1 presents the computed error associated with the selection of the best polynomial order (Bateni et al. 2003): where S is the error, j is the jth pixel of the interface, P is the number of pixels, O is the polynomial order, Y polynomial is the y coordinate of the approximation curve and Y pixel is the y coordinate of the pixel obtained from the pixel selection.
Some issues were encountered during the polynomial fitting approach. If more than one pixel has the same value in the x-direction, by definition, the function cannot be defined. For this reason, in each set of interface coordinates, first the polynomials are fitted, and the polynomial coefficients that give the lowest error are stored. Then, the axes are flipped by 90 o , and the polynomials and the errors are computed again. If the error improves, i.e., less than the previous value, the new coefficients of the polynomial fitting are stored. This procedure allows a better adjustment of the approximation curve.
( Fig. 1 Mapping process for pixel classification

Contact angle
The derivative of a polynomial gives the slope of the curve. Thus, with the polynomial coefficients obtained from the polynomial fitting, the tangent of the interfaces at the same contact point can be computed. The contact angle θ is then calculated, with the slope components in x and y directions, through Eq. 2.
where ⃗ u is the tangent vector of the brine/rock interface, and ⃗ v is the tangent vector of the brine/oil interface.
For each image, results from the different methods are obtained by changing the maximum polynomial fitting order and the image DPI (dots per inch), in order to explore the influence of image resolution. The contact angle is presented as the mean value among the computed angles since more than one point can be detected in the surroundings due to the way how the contact points are detected.

Materials
The images employed for the analyses are presented in this section: analytically generated images and rock sample images. (

Analytically generated images
Initially, the presented methods are evaluated using analytically generated images, whose contact angle is known. The circle is drawn with a known function, in which it is possible to compute its tangent at any point over the circle, and the flat surface has the same tangent value for all cases. Thus, the analytical contact angle can be computed through the known tangents according to Eq. 2. Errors in the estimation of the contact angle are calculated based on the analytical and computed values. The analytically generated images are shown in Fig. 3, where the oil is portrayed as a semicircle (green), the rock is a flat surface passing through the circle (black), and the brine is the remainder of the image area (white). The five images are generated by varying the distance between the surface and the circle center. Alhammadi et al. (2017)'s dataset is composed of highresolution X-ray micro-tomography reservoir rock samples obtained after 20 pore volumes of waterflooding at subsurface conditions. Figure 4 shows the selected images in grayscale from the original image dataset, with the colors black, dark gray, and light gray corresponding to oil, brine, and rock, respectively.

Results and discussions
This section presents the results for the images generated analytically and the slices of rock sample from Alhammadi et al. (2017).

Analytically generated images
For each method, Fig. 5 shows the error percentage for five different contact angles according to the five images shown in Fig. 3. The abscissa corresponds to the contact angle value from the analytically generated image; the ordinate is the error percentage between the analytical values and the computed ones; and each color of the marker corresponds to the employed range of the polynomial fitting order.
Regardless of the method, contact angles lower than 40° are poorly estimated. The presence or absence of pixels affects smaller angles since small differences in high curvatures can cause significant differences in the slope computation of its fitting polynomial. Although the difference of error percentage between the angles is significant, Method 2 shows stable behavior regarding the polynomial order due to pixel selection. In general, the polynomial order ranging from 1 to 4 shows smaller errors, especially in Method 3, where it displays the lowest overall error. Thus, Methods 2 and 3, with the polynomial order ranging from 1 to 4, are selected for the upcoming analyses.
The influence of image resolution is also explored by fixing the number of selected pixels. For the same chosen angles, Fig. 5 presents the error percentage of Methods 2 and 3, where each marker color represents five different DPI values, ranging from 100 to 500 DPI. Some images show accuracy improvement when increasing the DPI value (Klise et al. 2016), but, overall, the best cases are the images with 100 and 200 DPI.
In order to describe the same curvature for the same interface, the image with higher DPI needs more pixels than the one with lower DPI, which is natural since the pixel density is higher. In addition, the disposition of interface pixels also changes, although visually imperceptible. Thus, in this methodology, increasing the DPI does not necessarily improve the results. Other factors, such as the image segmentation, also contribute toward the methodology performance, as shown in Fig. 6. Now, we select a criterion for choosing which method should be applied to the actual rock sample image. We consider that the best method will be the one displaying the lowest deviation regardless of the polynomial order. Since methods 2 and 3 possess a more stable behavior regarding the polynomial order, they are chosen to calculate the contact angle of a weakly water-wet sample from a carbonate reservoir. The results for the generated images show that the right selection of pixels, performed during the image processing step, and the polynomial fitting order used to describe the interfaces have a strong influence on observed errors between the angle measured from the generated image and the angle calculated using our methodology. The image resolution also affects the image processing results, although not intuitively, as shown in Fig. 5.

Rock sample images
Methods 2 and 3, with maximum polynomial order of four, are applied in three images of the weakly water-wet sample from the segmented dataset of Alhammadi et al. (2017). Method 2 results are shown in Fig. 7, and Method 3 in Fig. 8. Applying the same principle as in the generated images, the X-ray images are segmented into three defined grayscale colors, where the color black correspond to oil, gray as brine, and white represents rock. The blue curves represent the fitted polynomial used to describe the interfaces, and the red lines are the tangents of the blue curves obtained at the contact point. The contact angle is computed as the angle formed between the red lines in the gray color (brine).
Differences between the two methods in the fitted polynomial curves are subtle. However, in Figs. 7b and 8b, the tangent line from the brine-rock interface differ significantly. The contact angle results using Methods 2 and 3 are presented in Table 2.
Similar contact angle values are observed for Figs. 6 and 7 (a) and Figs. 6 and 7 (c), while Figs. 6 and 7 (b) have a variation of 9.93 o due to the difference in the tangent from the brine-rock interface. Although the blue curve seems to    have a similar curvature in both images, the tangent has high sensitivity to small variations in the polynomial due to the employed methods for pixel selection.

Conclusions
This work focused on discussing the difficulties of image processing techniques used to determine contact angles. Different methods were presented while varying the range of the polynomial fitting and the resolution of the image.
Several types of approaches were tested in order to obtain the interfaces and the contact point. In order to evaluate our methodology, analytically generated images with known contact angles were used to validate our results, to estimate the error of our methodology, and to verify the best methods. After this, the best methods (Methods 2 and 3) were tested further with X-ray images of a rock sample. In most cases using analytically generated images, we observed that a polynomial with a maximum order of four generated better results, and proper pixel selection can reduce the error significantly during the polynomial fitting. Increasing only the image resolution did not assure an effective improvement of the results, since employing the same image processing steps do not keep the same pixels consistency. In addition, results from the rock sample images showed that the tangents obtained from polynomial fitting were very sensitive to the presence or absence of pixels.
Some issues were encountered during the study. Performing image processing is a hard task when precision is required in order to generate a proper image segmentation, and correct pixel selection is mandatory for right determination of the interfaces and the contact point since small changes in the pixel disposal can affect significantly the final result. For future work, we intend to expand our study and analysis to three-dimensional images.