Algorithm for post-processing of tomography images to calculate the dimension-geometric features of porous structures

An algorithm for post-processing of the grayscale 3D computed tomography (CT) images of porous structures with the automatic selection of filtering parameters is proposed. The determination of parameters is carried out on a representative part of the image under analysis. A criterion for the search for optimal filtering parameters based on the count of "levitating stone" voxels is described. The stages of CT image filtering and its binarization are performed sequen-tially. Bilateral and anisotropic diffuse filtering is implemented; the Otsu method for unbalanced classes is chosen for binarization. Verification of the proposed algorithm was carried out on model data. To create model porous structures, we used our image generator, which implements the function of anisotropic porous structures generation. Results of the post-processing of real CT images containing noise and reconstruction artifacts by the proposed method are discussed.


Introduction
Porous structures, i.e., structures containing pores of different volumes and shapes in the internal volume, are of great interest in the modern world. Studying the internal structure of such objects allows us not only to better understand the processes occurring in geological rocks [1,2], snow layers [3], ceramics [4], polymer structures [5], and in many other similar materials but also to synthesize porous structures with the given properties. The latter are used in industry as filters [6], in medicine as permanent or bioresorbable implants that ensure the budding of bone or vascular tissues [7], etc.
Knowledge of the spatial (3D) structure of a sample allows us to calculate such essential characteristics as its porosity (the ratio of pore volume to total volume), specific surface area (porous surface area per unit volume), permeability (the ability to filter fluids in the presence of a pressure drop) and pore size distribution, as well as to simulate the propagation of fluids in the sample (moving liquids and gases in geological samples, melting of snow, cleaning fluids with ceramic filters, etc.). Since the 1980s, the method of computed tomography (CT), as an innovative method, began to be used to visualize the internal volume of core samples to study the oil saturation of reservoirs [8].
This method allows us to analyze the dimensional and geometric features of microporous structures, based on which the characteristics mentioned above are calculated, regardless of their origin, pore location, and structure regularity. Unlike microscopic methods (electronic, Xray, or visible range), X-ray tomography does not require the destruction or special preparation of samples. A sample is probed with a wide beam of X-rays at different angles. The recorded projections are used to reconstruct the spatial distribution of the X-ray attenuation coefficient, which describes the morphological structure of the object. However, even if the studied porous object has a homogeneous composition, and the morphological structure is represented by empty pores and a homogeneous matrix material (a matrix is the basis of the porous structure, the material inside which pores exist), the reconstructed image will not be binary, but in gradations gray. This is explained by the statistical nature of the radiation source, the uneven response of the cells, and the limited spatial resolution of the recording equipment and artifacts of the reconstruction itself. To calculate the characteristics of a porous structure, such as porosity or surface area, as well as to simulate the hydrodynamic processes of the fluid's behavior in porous media, the segmentation of a tomography image, or its binarization if the matrix is homogeneous, is required [9 -12]. The question of choosing a segmentation / binarization method is not straightforward [13,14]; methods that work perfectly on one type of image, for example, documents [15], cannot be transferred entirely to CT images. Small errors in the segmentation of images of text documents, for example, in determining a threshold, in general, do not falsify the required information, while for porous structures, this can lead to significant errors in determining the porosity and specific surface area. For example, when, in a real sample, the pores are connected to each other. After such binarization, they can be separated, and vice versa, which will lead to errors in calculating permeability and simulating hydrodynamic processes. Thus, the quality of the segmentation of the original image is critical for studying porous structures.
Improving the binarization accuracy of a CT image can be achieved by performing its preliminary filtering from high-frequency noise while preserving the boundaries, since filtering methods that do not take into account the semantic structure of the image, for example, the Gaussian filter, can lead to strong distortions of the structure [16]. To filter images of porous structures, the following methods are usually used: 1. Median filtering [17 -19]. 2. Anisotropic diffusion (anisotropic filtration) [20 -23].
3. Bilateral filtration [24 -26]. 4. Guided filtration [27,28]. The listed methods are multiparametric; therefore, the question is to choose the optimal values of the parameters for each of the methods. The quality of the resulting image depends on this. For example, if parameters # 1 are used in the anisotropic filtering method (the number of iterations (niter) = 1; the conductivity coefficient (kappa) = 5; the diffusion rate (gamma) = 0.1), then the resulting image ( Fig. 1b) will slightly differ from the original image (Fig. 1a); consequently, all problems associated with noise during reconstruction will remain. If parameters # 2 are used (niter = 20; kappa = 120; gamma = 0.25), then we will obtain a very blurry image, which will lead to a significant distortion of the object boundaries (Fig. 1c). Visual analysis of the images leads to the conclusion that random selection of parameters results in random quality. We failed to find algorithms for finding the optimal parameters. Perhaps, this is because, in order to solve the problem, it is necessary to construct a criterion by which the filtering quality will be fulfilled, and in the absence of knowledge about the real morphological structure of the studied object, the choice of such a criterion is far from obvious.
There is noise on the grayscale CT image, and, as a result, artifacts (distortions of the real structure of the object) appear during the binarization process. In the previous publication [29], such artifacts were called "levitating stones" -these are conglomerates of voxels classified as a matrix but entirely surrounded by voxels of pores. Artifacts in the form of phantom closed pores also appear along with them. Examples of such artifacts are shown in Fig. 2. The appearance of artifacts after binarization is not a problem exclusively for the analysis of CT images but also occurs when processing other types of images [30,31].
In our work, we propose using the estimate of the voxels number of such "levitating stones" as a criterion for finding the optimal values of the filtration parameters since, unlike closed pores, the existence of "levitating stones" is not physically possible and, therefore, can be an indicator of the image processing quality. A separate task in processing images of porous structures is the search for a representative volume (RV). As a rule, real reconstructed images have sizes of 109 voxels or more, and multiple processing of such an image in full volume in the process of finding the optimal values of the filter parameters can be very laborious. To reduce the running time of the algorithm, it is proposed choosing an RV that has a smaller size than the original image, on which the parameters will be determined.
Usually, some controlled value is used as a parameter for determining RV. As a rule, when analyzing porous structures, the value of the porosity or specific surface is used [32,33]. In this case, the minimum size of an image fragment at which the controlled value with some error is equal to the same value calculated over the entire image is taken as the RV. This is not suitable for us, because, firstly, to calculate such a value, the image should be first filtered and binarized, and for this, the filtering parameters that we want to determine should be selected, which leads us to a vicious circle. Secondly, we are only interested in the presence of a material / pore boundary in the RV and a significant (representative) noise value and this may require a much smaller volume than in the case of determining the RV using the standard approach described above.
We take the minimum image size as the RV, which contains sufficient information about the structure of a porous sample and artifacts appearing in the process of binarization [34].
In this paper, we propose an algorithm for finding the optimal values of the filtering parameters based on the criterion of the number of "levitating stone" voxels, which includes the following steps (a detailed description of the steps is given in the corresponding sections of the paper): 1. Binarization of the initial tomography 3D image (see "Binarization"). 2. Search for the RV, on the obtained binary image (see "Method of searching for a representative volume"). 3. Determination of the optimal values of the filtration parameters by the simulated annealing (see "Search for the optimal values of the filtration parameters"). At this stage, the following are iteratively performed:  RV filtering (see "Filtration. Description of filters and their parameters");  Counting the number of "levitating stones" (see "Description of the method for detecting "levitating stones" and closed pores"). 4. Filtration and binarization of the tomography 3D image using the found parameter values. To verify the quality of the algorithm operation, we used software-generated porous structures (phantoms) with a subsequent imitation of a real tomography experiment to obtain noisy 3D images of porous samples. A description of the phantom generator is given in the section "Algorithm for generating model porous structures". The procedure for simulating the experiment is presented in the sections "Noise on the CT image due to an imperfect source and detector" and "Simulating a tomography experiment". A comparison of the image obtained as a result of the algorithm operation with the initial phantom was carried out according to the Jaccard coefficient (see "Criterion for evaluating the quality of the algorithm operation").
The results of testing the algorithm using two types of filters (anisotropic and bilateral) are presented in the section "Testing the algorithm on model CT images". Using the median filter did not show high results. The guided filtering is difficult to implement.
We also applied the proposed algorithm for processing CT images of real ceramic porous structures fabricated by self-propagating high-temperature synthesis [35] (see "Processing of the CT image of a cermet membrane").

Algorithm for generating model porous structures
To create model porous structures, we have implemented a phantom generator with the ability to set anisot-ropy (different pore sizes along three mutually orthogonal directions). Its distinguishing feature from existing generators [36 -38] is the ability to set any angle of rotation and inclination of the anisotropy axes relative to the axes of the 3D image. Geological porous materials often demonstrate a certain degree of anisotropy due to rock formation processes [38]. Anisotropy is also observed when studying the porous structures of snow [39]. Thus, when creating model phantoms of porous samples, it is necessary to be able to set anisotropy, and in the general case, anisotropy of the sample can be multicomponent, internally non-orthogonal, and non-coplanar to the axes of the phantom voxel array. A requirement for the possibility of specifying anisotropy in three-dimensional space with arbitrarily chosen directions of anisotropy vectors follows from here. In our case, it has been realized by creating a separate kernel of a three-dimensional Gaussian filter, properly oriented in space, for subsequent convolution with an undistorted image. To generate undistorted porous structures, a modified algorithm of porespy.generators.blobs from the PoreSpy package has been used [36].
The following sequence of steps has been implemented. 1. A three-dimensional array of the required size is created, filled with uniformly distributed random numbers from the interval [0, 1). 2. A kernel of a Gaussian filter is created: the specified scale ratios of the array size to the pore size and material size determine the value of the rootmean-square deviation σ of the Gaussian along each of the three axes, and the direction of the anisotropy axis determines the tilt angle of the vertical axis of the filter. The filter size for the data used in this work was ± 4 root-mean-square deviations. 3. The convolution of the original image, created using a random number generator, is conducted with the resulting Gaussian kernel. 4. To optimize the computation of the convolution, a check is made for the orthogonality of the anisotropy axis to the axes of the array and, when the condition of separability and symmetry of the kernel is met, the three-dimensional convolution was replaced by a combination of three onedimensional convolutions. 5. The resulting 3D image is subjected to histogram equalization. 6. Binarization of the 3D image is carried out according to the threshold equal to the value of the desired porosity (the ratio of the pore volume to the total volume) of the phantom. For further simulating the tomography experiment, the absorption coefficient of the material was taken equal to 0.25 mm -1 , which corresponds to the absorption of Al 2 O 3 at a probe radiation wavelength of 33 keV.
Noise on the CT image due to an imperfect source and detector In the method of X-ray computed tomography, the final image presented as a result of the operation of the to-mography method is not the image recorded by the detector, but the result of applying a multi-stage mathematical apparatus to projections measured at different angles, where only the final step is the reconstruction itself. Due to the imperfection of both the radiation source and the detector, the recorded projection images contain noise, which is converted into noise in the final 3D image during the reconstruction process.
The attenuation of X-ray radiation as it passes through an absorbing medium is described by the Beer-Lambert-Bouguer law. If an infinitely thin monochromatic ray with a unit intensity passes in the direction given by the straight line L through the medium described by the distribution of the linear attenuation coefficients μ, then the output intensity value A can be calculated according to the expression: where ρ is the vector that determines the position of the detector pixel relative to its center, φ is the projection angle at which the object is probed, and dl is the segment of the straight line L. The process of generating X-ray quanta by a laboratory source is well described by the Poisson process [41]. The value in the pixel of the image submitted for reconstruction can be calculated as follows: where S is the expected value of the brightness in a pixel, I 0 is the intensity of the source determined by the anode material and the source mode used (voltage, current), P(I 0 ) is the realization of the Poisson process, n is the number of frames made (tomography imaging frames), and i is the index number of shooting frame. This expression assumes that the response for each pixel is generated independently. In practice, this cannot be realized due to the diffraction of radiation on the sample and the diffusion of charge in the cells of the CCD matrix. To take into account the influence of these effects in the simulation, the operation of the convolution of the calculated projection image with a two-dimensional Gaussian filter with a size of 3 × 3 pixels is introduced: Here, G is the Gaussian function. Sigma is selected in such a way that 60 % of its original signal remains in the central pixel, and the remaining 40 % are distributed over neighboring pixels [42]. This mathematical model of signal formation on the image is used to calculate projection images from the generated anisotropic phantoms.

Simulating a tomography experiment
Obtaining porous objects for adjusting the operation of the algorithm is carried out by the numerical simulation of a tomography experiment. From the generated phantom of the porous object (where the absorption coefficient of the simulated material is set for each voxel of the matrix), tomography projections are calculated using the Radon transform. The resulting projections are modified taking into account the parameters of the intensity of the radiation source in accordance with the above mathematical model of image formation. The resulting noisy projection images are given for reconstruction, which is performed by the filtered back projection algorithm (FBP). This is described in more detail in [34]. As a result, two 3D images are obtained -the original phantom and the corresponding reconstructed noisy image. The presence of a phantom allows us to evaluate the quality of the results obtained after filtering / binarization with the original image (the Jaccard coefficient is used) and draw conclusions about the correctness of the algorithm operation. It is also possible to trace the "life history" of "levitating stones" from these images -some of them are initially matrix voxels and "split off" in the process of binarization. The existence of such "split stones" indicates the impossibility of naive processing of the obtained images by removing all "levitating stones" after binarization.

Binarization
The tomography reconstruction process produces a grayscale image at the output. To determine the interested parameters of the porous structure, its segmentation or binarization in the case of a homogeneous matrix is necessary. In this work, the unbalanced Otsu method is used for binarization [43], which is a global binarization method that divides the histogram of image brightness into two classes, using the assumption that these two classes are realizations of normal distributions with different weights.

Method of searching for a representative volume
In order for the result of our algorithm operation with the optimal values of the filtering parameters obtained on the image fragment to coincide with the result of the algorithm operation, for which the parameters are determined over the entire image, the fragment should be representative. In other words, it should contain a sufficient number of elements that are determinative in the choice of filter parameters. Since the original object is twocomponent (material / pores), we need the presence of a border on the selected fragment, and we also need a sufficient number of "levitating stones", as their number is a criterion for determining the optimal parameters.
A representative element (RE) as a fragment of an image with a size equal to RV is determined. The search for the RE is carried out as follows. Rows of voxels are randomly se-lected from the image along a given direction (say Z), which, in turn, are divided into fragments of different lengths. Let us introduce the repeatability function R (ϵ, L f ), which shows how often the selection contains fragments of length L f that are similar to each other, the degree of similarity ϵ ranges from 0 to 1. The value of RV along the chosen direction depends on the accuracy parameter ϵ.
The dependence of R (ϵ, L f ) on L f for porous, not strictly periodic structures is a decreasing function that reaches a plateau. The position of the plateau level varies depending on the periodicity of the structure. For clarity, the dependence of R (ϵ, L f ) on L f can be normalized from 0 to 1, as shown in Fig. 3 (||R(ϵ, L f )|| is the normalized value of the function R(ϵ, L f )). It can then be argued that the curve reaches a plateau at the shortest length min(L f ), where the deviation of the normalized value of the function from its global minimum is no more than 1 %.

Fig. 3. Repeatability function dependence from length
The min(L f ) value, where the graph of the normalized function R (ϵ, L f ) reaches a plateau, determines the representative elementary length (EL) along the selected direction. It follows from this that any fragment of a shorter length than the EL has a relatively high repeatability. The larger ϵ, that is, the more accurate the match must be for the selected fragment with the length of L f , the smaller the obtained EL will be, and vice versa, the smaller ϵ is (the lower the accuracy requirement is), the larger EL will be. The value ϵ = 0.5 is used. Similarly, we calculate the EL values along the mutually perpendicular directions X (ELX) and Y (ELY), orthogonal to Z (ELZ). The size d of the RV is expressed as: After that, RE should be selected. To do this, a fragment of the image with a size equal to RV is cut, taking into account the fact that the specific number of "levitating stones" of the selected fragment N Vr should be close to the specific amount averaged over the entire image N V : The requirement for RE to contain at least one closed pore voxel is guaranteed by the presence of boundaries in RE, which is required for the correct selection of the optimal values of the filter parameters. An example of RE for a reconstructed tomography image is shown in Fig. 4.

Filtration. Description of filters and their parameters
The appearance of the mentioned above binarization artifacts in the form of "levitating stones" and phantom closed pores is the result of an incorrect classification of the pixels of the reconstructed image u with a high noise level. To eliminate them, it is necessary to apply noise reduction methods to the reconstructed images. In this work, we used two noise reduction methods with saving boundaries -anisotropic and bilateral filtration.
The principles of filter operation will be described further. We start with a description of the anisotropic filter [44]. It saves boundaries in areas of the u image with a large gradient value. Such areas usually correspond to the boundaries of objects in the image; in our case, these are the boundaries between the matrix and the pores. The principle of the filter operation is based on representing the process of blurring an image as a diffusion process, defined by the following expression: where u N is the reconstructed image u at the N filtration iteration, λ is the diffusion rate coefficient, NSEW are four directions, and g (|u|) is the diffusion coefficient. The diffusion coefficient [44] is as follows: where K is the user-defined conductivity coefficient, determined empirically. This filter is iterative, which can be attributed to its disadvantages, as it slows down the calculations. Before describing the bilateral filter, a generalized form of the filtration process b i =∑ j W ij a j is introduced. Here, a and b are input and output images, respectively, i and j are pixel indices, and W is the filter kernel. The bilateral filter [45] uses not only information about the spatial proximity of pixels in a given window, like a classical Gaussian filter, but also about the proximity of pixel brightness values. The filter kernel is defined as follows: where k is the normalization parameter, p i and p j are vectors specifying the coordinates of i and j pixels of the u image, respectively, u i and u j are brightness values of the i and j pixels, and the parameters σ s and σ c determine the spatial proximity and proximity of the brightness values, respectively.

Search for the optimal values of the filtration parameters
Both filters considered above have three user-defined parameters, and the question is of which parameter values should be chosen for processing specific images. As a criterion for choosing the optimal values of the filter parameters, we propose to use the number of "levitating stone" voxels, i.e. the optimal values are those at which the number of "levitating stones" is minimal.
As shown in [46], the function of the dependence of the number of "levitating stones" on the values of the parameters of the applied filter has local minima, while the method of enumeration of values either has the probability of missing the minimum due to too large a step, or its application is redundant (computationally expensive) if the step is small. Numerical methods are usually used to solve the minimization problem, but there is no analytical expression for the relationship between the parameters and the number of "levitating stones". In this case, the socalled "black box" optimization methods are used [47]. In this work, we use the simulated annealing method [48], as it is a well-proven method with respect to problems similar to ours.
The principle of operation of the simulated annealing method consists in a step-by-step descent to the minimum of the value of the objective function, while the method has the ability to "exit" from local minima due to the rise to a greater value. To apply the method, it is necessary to determine the type of optimized function, the function of choosing new parameters, and the temperature function, according to which the velocity of descent and the probability of transition to a higher value are set. The following sequence of steps is used: 1. An initial triplet of parameters is selected from the given intervals of parameters, randomly, according to a uniform distribution, and, using the algorithm proposed by the authors, described below, the number of "levitating stones" obtained after filtering is calculated; 2. A random set of parameters is determined according to the Gaussian distribution with the average values selected in the previous step, and the number of "levitating stones" for it is calculated;  If the received value has decreased, then this set of parameters is considered the current minimum.  If the value has increased, then the probability of transition to that "higher" position is calculated. This probability depends on the given value of the temperature function and the distance of the current set of filter parameters in the parameter space. The lower the temperature and the more the considered set of parameters differs from the previous one, the lower the probability of transition. If the transition is made, then this parameter set becomes the current minimum; otherwise, this parameter set is discarded. This step is needed in order to avoid a descent to a local minimum. 3. The temperature decreases in accordance with the selected law (we used a linear one), and Step 2 is repeated until the temperature is higher than the specified minimum.
To select the range of values of the anisotropic filter parameters, an analysis of the literature on the topic [49 -54] was carried out. The conductivity coefficient K is the threshold by which it is determined whether a pixel belongs to the boundary, so its values can vary within the image brightness, for eight-bit images -from 0 to 255. When the diffusion rate coefficient λ > 0.25, the filter is unstable, and when λ < 0.01, the filtering effect is not noticeable [49]. The analysis also showed that, even at small λ values, it makes no sense to make the number of iterations more than 30 due to excessive blurring of the image [49 -53]. The choice of the range of parameter values for the bilateral filter was carried out according to [54]. Selected ranges: a window size w from 3 to 7; σ s from 0.4 to 1.0; σ c from 0.4 to 2.0.

Description of the method for detecting "levitating stones" and closed pores
For each voxel of the image  with integer coordinates (x, y, z) in three-dimensional space, there is a set N 6 () of 6 voxels with adjacent faces, the coordinates of which are presented as follows: voxels with adjacent vertexes: voxels with adjacent edges: as well as a set of 26 voxels-neighbors N 26 (), which have both common edges and common vertexes, and whose coordinates are specified as: 26 6 N ( ) N ( ) N ( ) N ( ).
According to [54], in the binary image, two voxels α and β are m-adjacent (m = 6, 26) if the β element is included in the set N m (), and and β have the same brightness. A discrete m-path from voxel A to voxel B, which has the same brightness as A, is a set of m-adjacent voxels, which has elements m-adjacent to voxels A and B. Two voxels of the same image brightness will be called m-connected if there is an m-connected discrete path for them; an m-connected set of voxels R is defined as a set of voxels of the same brightness if any two voxels from this set are m-connected.
Based on the above definitions, the matrix of the porous structure is m-connected voxels of brightness 1, which are m-connected to the edge of the image. Likewise, the open pore system will be an m-connected set of voxels of brightness 0, which is m-connected with the set of pore voxels at the edge of the image.
In our algorithm, a "levitating stone" is defined as a voxel or a group of voxels that are not m-connected to the matrix of the porous structure. In turn, a "closed pore" is defined as a voxel or a group of voxels that cannot be called m-connected to the system of open pores. In this paper, we have chosen 6-connectivity to define "levitating stones". We have obtained stable results after simulating a tomography experiment, according to which the following conclusion has been made with confidence: we find more "levitating stones" on phantoms with a high level of noise than on phantoms with a low level of noise. When 26-connectivity is used, the result of determining the number of "levitating stones" is unstable and inconsistent; a phantom with a low noise level may have more "levitating stones" than a phantom with a high noise level. This is due to the fact that, after binarization of a grayscale image with a high noise level, a large number of voxels, falsely referred to as the "matrix" class and located deep in the pore volume, are connected through the vertexes with their neighbors, and through them with the matrix of the porous structure. Thus, these voxels are mistakenly classified as "matrix" rather than "levitating stones".

Criterion for evaluating the quality of the algorithm operation
To compare the image obtained as a result of the proposed algorithm with the original phantom, the Jaccard coefficient, defined as the ratio of the intersection of two sets (for example, A and B) to their union, is used: In the case of complete coincidence in two binary images, the Jaccard coefficient is 1. If the images do not have the same areas, then the coefficient will be 0. The intermediate results, by their proximity to 0 or 1, characterize the degree of similarity between two images.
The Jaccard coefficient, chosen as a measure of the similarity of the images before and after the filtering stage, correlates well with the number of "levitating stone" voxels (Fig. 5). Each point in the figure is calculated with a random set of bilateral filter parameters taken from the specified range in the section "Search for the optimal values of the filtration parameters".

Testing the algorithm on model CT images
To check the operation of the proposed algorithm, the generated phantoms and the corresponding model reconstructed 3D images are used. The size of each image is 1000 × 1000 voxels in the horizontal section and 500 voxels in height. The porosity of the sample is set equal to 30 %. The simulation of the experiment is carried out using the following parameters: 360 projections with a step of 0.5 degrees and an intensity of the radiation source I 0 equal to 1000, 1500, and 3000 photons. The reconstructed images for 1000, 1500, and 3000 photons are named Sample 1, Sample 2, and Sample 3, respectively. Examples of the obtained two-dimensional sections of three-dimensional images are shown in Figure 6. The pores are represented by black pixels. For the chosen parameter of the accuracy of determining RV ϵ = 0.5, the dimensions of RE of our phantom objects for the above parameters are 353, 453, and 703 voxels, respectively. The images are shown in Fig. 6. First, the reconstructed image is binarized using the method described above. The obtained binary image is used to calculate the number of "levitating stones" and determine the Jaccard similarity coefficient with the original generated image (phantom). The number of voxels with "levitating stones" is 2.3 × 10 7 , 2.1 × 10 7 , and 1.6 × 10 7 (hereinafter, the triplets of obtained values for sample1, sample2, and sample3 are given, respectively). At the same time, the number of "stones" that originally belonged to the matrix, but are "split off" after the simulation of the experiment and subsequent binarization is 185341, 86018, and 17640, respectively. The Jaccard coefficient for the obtained binary images is 0.6037, 0.6592, and 0.7702. The results obtained after the work of the proposed algorithm are further compared with these indicators.
First of all, the idea of direct removal of all voxels of "levitating stones" is verified -after such a procedure, the values of the Jaccard coefficient are 0.6038, 0.6595, and 0.7722, which practically does not differ from the initial values and suggests that it is necessary to apply more sophisticated approaches to image processing.
Since the annealing method does not guarantee the stability of the result, five runs of the algorithm were carried out on each model sample, and the obtained results were averaged. The bilateral filter and the anisotropic diffusion filter were also tested separately. Fig. 7 shows images of fragments of the original phantom and the corresponding fragments of images of different stages of the algorithm operation.
As a result of the algorithm operation, the average number of "levitating stones" in the case of the bilateral filter decreased on average to 7.3 × 10 4 , 3.6 × 10 4 , and 280, and in the case of the anisotropic diffusion filter, to 275, 152, and 206 (Fig. 8a). The number of "split off stones" was 121, 37, and 12 for the bilateral and 195, 33, and 3 for the anisotropic filter -a decrease in this value indirectly indicates that the resulting matrix structure is closer to the initial one than it was before filtration. The Jaccard coefficient of the processed images averaged 0.9774, 0.9812, and 0.9863 for the bilateral filter and 0.9806, 0.9828, and 0.9858 for anisotropic diffusion (Fig. 8b). As shown, the anisotropic diffusion filter gives results comparable to the bilateral filter, but in terms of velocity, it is inferior to the bilateral filter -it takes on average twice as much time to process the image, which is most likely due to the iterative essence of the anisotropic filtering method. The obtained verification results allow us to conclude that the proposed algorithm gives stable results and can be used to process tomography images of real porous structures.

Processing of the CT image of a cermet membrane
The proposed method was tested for post-processing of a CT image of a cermet membrane made of Al2O3 [29]. The measurements were carried out on a microtomograph of the Federal Research Center "Crystallography and Photonics" RAS. The image was reconstructed using our own software [55,56].
The porosity value obtained by hydrostatic weighing is 30.5 %. The specific surface area measured by the BET method (Brunauer-Emmett-Teller) was 0.2 m 2 /g. The pores do not exhibit pronounced anisotropy. The structure contains single large pores, which creates porosity inhomogeneity in the sample.
A fragment of the reconstructed image (representing a rectangular parallelepiped with dimensions of 500 × 500 × 700 voxels) was binarized and filtered, followed by binarization. The filter parameter values were determined according to our proposed method. The size of the RE was 97 3 . As a result of processing with an anisotropic filter, the number of "levitating stones" decreased from 1.5 × 10 7 to 1.7 × 10 5 ; with a bilateral filter, the number decreased to -6.8 × 10 5 . The porosity and specific surface area values measured from the obtained image were 30.3 % and 0.207 m 2 / g, respectively, which is in good agreement with the values measured by standard methods of materials science. Fig. 9 shows the corresponding images. When processing a real porous structure, the need to distinguish between the concepts of RV and RE is especially pronounced. If the distribution of pores and material can be considered uniform, a fragment of the RV can be "cut out" from any place, and the RE can be obtained. If the distribution cannot be considered uniform, then the RE cannot be chosen arbitrarily. We have presented in the paper a criterion for its correct choice. Figure 9 shows that there are several large voids. An RE containing only voxels of pores from these voids can lead to an error in the search for optimal values of the filtration parameters. Since there are many triplets of parameter values and not all of them are optimal, some of them lead to complete image blurring. Therefore, the algorithm proposed in the paper, in which the work is carried out not only with the dimensional parameters of the RV but also with a correctly selected RE in the volume, is expedient.

Conclusion
This paper presents a new algorithm for postprocessing 3D CT images of porous structures. To optimize the time required to calculate the optimal values of the image filtration parameters, their calculation was carried out on a representative part of the analyzed image. The stages of CT image filtering and its binarization were performed sequentially. The annealing method was used to find the optimal parameters at the filtration stage. Bilateral and anisotropic filtering was implemented, and the Otsu method for unbalanced classes was chosen for binarization. The verification of the proposed method was carried out on model data. To obtain model images, a generator of porous structures was implemented. The Jaccard coefficient was used to assess the proximity of the model and processed images.
The results of post-processing by the proposed CT method of images of a cermet membrane, measured on a microtomography setup of the Federal Research Center "Crystallography and Photonics" RAS, showed that the porosity value calculated from the post-processed image agrees with the value obtained by measuring the porosity by the alternative method.