Abstract
Thresholding is an efficient step to extract quantitative information since the potential artefacts are often introduced by the point-spread effect of tomographic imaging. The thresholding value was previously selected only relying on engineering experience or histogram of tomographic image, which often presents a great challenge to determine an accurate thresholding value for various applications. As the tomographic image features often do not provide sufficient information to choose the best thresholding value, the information implicit in the measured boundary data is introduced into the thresholding process in this paper. A projection error, the relative difference between the computed boundary data of current segmentation and the measured boundary data, is proposed as a quantitative measure of such image segmentation quality. Then, a new multistep image segmentation process, called size projection algorithm (SPA), is proposed to automatically determine an optimal thresholding value by minimising the projection error. Results of simulation and experiment demonstrate the improved performance of the SPA-based tomographic image segmentation. An application of size projection algorithm for gas-water two-phase flow visualisation is also reported in this paper.
1. Introduction
Tomography is a highly promising technique for providing 2D/3D internal configurations of physical objects utilising the feature of “seeing through” [1]. Particularly, electrical impedance tomography (EIT) has been successfully applied in many industrial [2] and medical applications [3, 4] by offering a relative fast, cost-effective, nonintrusive, and visualising solution. For a given “object,” the measured boundary data of the object are acquired along with a certain sensing strategy and boundary electrodes, and then an image of the object (i.e., tomographic image) is reconstructed from the set of measured boundary data [5]. Here, in EIT technique [1, 6], the measured boundary data are electrical voltage signals, and the tomographic image is represented by conductivity difference in grey-value levels. In this paper, we focus on the tomographic reconstruction of single-material/tissue object, such as bubble/oil/solid in water [7] and lung in human thorax [3, 6, 8].
The tomographic reconstructions, usually continuous grey-level or colour palette-based images, need be segmented for extracting quantitative information of objects, such as the object size or shape. Although many image segmentation algorithms exist [9], such image segmentation is usually performed to categorise pixels into two different classes using a thresholding value. Commonly, a fixed global value based on engineering experience or simulation [10–12] is simple and direct for tomographic image thresholding, but its selection is often somewhat arbitrary. For example, Wang et al. utilised an empirical value to estimate the central air-core size in a working hydrocyclone [10] and to image large bubbles in flowing water [11]. Afterwards, many histogram-based methods [13–15] were employed to automatically select the “optimal” thresholding value corresponding to various optimisation measures. For example, a typical histogram-based method (Otsu’s method) was used to remove the background from the target regions [16, 17] or to extract the region of interest (ROI) [18].
To the authors’ knowledge, the previous thresholding value selection methods only trust the information of tomographic reconstruction, while abandoning the information implicit in the measured boundary data. However, these tomographic images suffer from serious reconstruction artefacts. For the mathematical model of EIT [19], a linear approximation was introduced to solve the nonlinear inverse problem, based on the condition of an infinitesimal change in conductivity. But the precondition is often not satisfied in the actual reconstruction process, which introduces a certain artefact in the tomographic reconstruction. Moreover, the ill-posed nature of inverse problem [6] results in worse reconstruction quality of tomographic images. In order to reduce the impacts of these artefacts in the selection of thresholding value, it is considered that the measured boundary data are introduced into the image thresholding process. By computing the boundary data with respect to the segmented image using a forward solver, the relative difference between the computed boundary data of current segmentation and the measured boundary data, called projection error, can be calculated. This yields a quantitative measure of such image segmentation quality, that is, the relative difference between the segmented image and actual image of objects. From this point, a new multistep image segmentation process, called size projection algorithm (SPA), is proposed for automatically determining an optimal thresholding value by minimising the projection error.
The rest of this paper is organized as follows. In Section 2, the image segmentation problem for EIT tomographic reconstruction and our SPA-based image thresholding method are described. Simulation has been conducted to evaluate the performance of the proposed image thresholding method in Section 3, and the proposed method is employed in phantom experiment and an application of gas-water flow in Section 4. Section 5 concludes the paper.
2. Methods
2.1. Tomographic Reconstruction
In EIT technique, according to the measured boundary data (i.e., measured boundary voltages), we want to obtain the actual distribution of materials (i.e., actual object distribution) which is often unknown in practice. Commonly, various tomographic algorithms [20] are employed to reconstruct a grey-value image (i.e., computed object distribution) which is often a 2D grey-level or colour palette-based image. The 2D image was demonstrated and widely used as an approximate representation of the cross-sectional configuration of the actual material distribution even if it is not exactly the same. As the tomographic image contains potential reconstruction artefacts, it only reveals qualitative information (i.e., location and shape of object), as shown in Figures 1(a) and 1(b).
(a)
(b)
(c)
The specific motivation of this work is to segment the tomographic image for extracting quantitative information, such as the size and boundary of objects. In what follows, the tomographic image (i.e., ) to be segmented is a tomographic reconstruction of a single-material object which can be expressed by (i.e., {}). Here, 0 and 1 are the absence and existence of object in the background field. From the domain of measurement, the measured boundary data with respect to are collected by a typical EIT system. While, from the aspect of simulation domain, the computed boundary data with respect to are calculated through a numerical forward solver (see Appendix A) with the employment of an EIT sensor and sensing strategy. A comparison of the boundary voltages in measurement domain and simulation domain is conducted in Figure 1(c).
2.2. Image Segmentation by Thresholding Method
As mentioned above, the grey-value image (see Figure 1(b)), we want to segment, is a tomographic reconstruction of an actual single-material object (see Figure 1(a)), whose axially sliced representations are compared in Figure 2(a). The grey-level histogram of the tomographic reconstruction is shown in Figure 2(b). An appropriate thresholding value is necessary for approaching the segmented image to the actual material distribution when the tomographic reconstruction is segmented. Such segmentation can be defined aswhere is the input image and is the segmented image. is a thresholding function which only depends on the grey value of at the point. For simplicity, let r and s denote the grey value of and at any point, respectively. The transformation in equation (1) can be expressed by grey-level mapping function (equation (2)), which is depicted in Figure 2(c):
(a)
(b)
(c)
With the processing of grey-level mapping function , the tomographic image is transformed to a binary image where the pixels belonging to objects and background are filled with grey value 1 and 0, respectively. The transformation is defined by
By introducing the conductivities of the object and background, a conductivity vector corresponding to the segmented image can be obtained by filling the object conductivity into object-occupied pixels and filling background conductivity into object-unoccupied pixels, as expressed bywhere and are the object conductivity and background conductivity, respectively. After filling the physical parameters of materials, the computed boundary data corresponding to the segmented image can be calculated by a numerical forward solver (see Appendix A).
However, the selection of appropriate thresholding value is a key issue, which affects the accuracy of object extraction. If the thresholding value is larger than the best one, the extracted object in the segmented image will be smaller than the real one and computed boundary data and measured boundary data will have a larger difference, and vice versa.
2.3. Thresholding Value Selection by Projection Error Minimisation
2.3.1. Formulation of Projection Error Function
As demonstrated in Figure 1(c), there is a little difference in boundary voltages between the measurement domain and simulation domain even for the same setup. A concept of relative voltage change was employed to reduce the difference effect and produced an error function to evaluate the approximation of an EIT inverse solution in the literature [19]. It transforms the relative change of a measurement domain to a simulation domain, which significantly reduces the error effects from measurement noise, electrodes inconsistency, and field dimensions [21]. Referring the same concept, a projection error function is proposed to evaluate the similarity of extracted object in the segmented image and the real object, which can be expressed bywhere is the projection error function with respect to the step segmentation. and are the measured reference voltage and measurement voltage, i.e., boundary data in measurement domain. and are the computed reference voltage and measurement voltage with respect to the conductivity distribution and , i.e., boundary data in the simulation domain. and refer to the number of boundary voltage and the pixel number, respectively.
The projection error function represents the difference between relative changes of measured boundary data and computed boundary data with respect to the segmented image, which should be minimal if the segmented image is close to the actual distribution. From this point, the projection error is used as a quantitative measure of the image segmentation quality.
2.3.2. Size Projection-Based Thresholding
In this paper, we focus on an image segmentation issue that a machine-determined thresholding value is to be approached for extracting the object in EIT reconstruction by minimising the projection error, as illustrated in Figure 3. Given a typical EIT system, the boundary voltages of an undermeasured object are collected and are employed to reconstruct an image representing the object distribution. A changeable thresholding value is used to segment the reconstructed image for eliminating the blurred boundary between the object and background. In order to evaluate the accuracy of such segmentation, the segmented image is converted to conductivity distribution by introducing the conductivities of object and background, and the boundary voltages corresponding to the segmented image are calculated via a numerical forward solver (see Appendix A). Following equation (5), a projection error is then calculated to measure the similarity between the segmented image and the real object distribution. Through changing the thresholding value, that is adjusting the projected size of the object in EIT reconstruction, the extracted object can be better approaching the real one once the minimal projection error is reached. The above process transforms an image segmentation optimisation problem to a projection error minimisation problem which can be solved by an optimisation method.
Considering the meaningful range of thresholding value, should be constrained to [0, 1] relative to the grey-level of pixel in a reconstructed image. Using the above definitions, the central problem of this paper can be summarised. Let be a variable parameter, let be a grey-scale image, and let be a vector of measured boundary voltages. Find an optimal , such that is minimal, which is a nonlinear unconstrained optimisation problem, as expressed by
As illustrated in Figure 2(a), the segmented object size in EIT reconstruction is first close to and then deviates from the real object size when changing the thresholding value from 0 to 1, that is, the projection error first decreases and then increases. Therefore, the projection error function is a unimodal function with respect to the definition domain [0, 1].
Golden-section search (GS) method, an efficient one-dimensional search technique, provides general means of solving the optimisation problem expressed by equation (6) since the range [0, 1] is a unimodal interval of projection error function. The employment process of the golden-section search (GS) algorithm is illustrated in Algorithm 1, where the control factors are set with a minimum convergence error (i.e., measurement error) and maximum number of search steps.
|
By the employment of the GS method, an optimal thresholding value is obtained by approaching the minimal projection error, at which the object size in the segmented image is closest to the real object size. As the thresholding value is machine-determined by comparing the computed boundary voltages corresponding to the segmented object size with the measured boundary voltages of real object size (i.e., the process of minimising projection error), the proposed thresholding value selection method for EIT image segmentation is named size projection algorithm (SPA).
3. Simulation
To investigate the accuracy and the limitation of the size projection algorithm (SPA), typical setups were modelled in COMSOL simulation. The corresponding images were reconstructed by commonly used EIT algorithms and segmented by the size projection algorithm. The segmented images from the SPA method were also compared with the other methods. The typical setups are based on (a) two objects with different size, (b) only a large object at centre, (c) four uniformly distributed small objects, where the conductivities of objects (air) and background (salt water) are 0 mS/cm and 1 mS/cm, respectively. In addition, antinoise capability and time consumption of SPA were also investigated. In the simulation, a typical 16-electrode EIT sensor and adjacent sensing strategy were employed, which could produce 104 independent boundary voltages for solving an EIT inverse problem. All processes of image reconstruction and image segmentation were conducted under a mesh with 1536 triangular pixels.
3.1. Accuracy
In order to compare the accuracy of SPA and other methods, such as the empirical value method [7, 10] and histogram-based method [15–17], relative image error (IE), as expressed by equation (7), is employed to evaluate the accuracy of these image segmentations:where and are the true conductivity distribution and the conductivity distribution of the reconstructed image or the segmented image, respectively.
The simulation results are listed in Table 1. For the two different size objects (see the first column in Table 1), the EIT image was reconstructed by the Tikhonov regularisation algorithm, where the object boundary is blurry and the small object is suppressed in the reconstructed image. Through three image thresholding methods, the clear boundary of objects is extracted. The extracted objects in segmented images are underestimated and overestimated by the empirical value method and the histogram-based method (Otsu’s method), respectively. However, the SPA-based method gives a better estimation of extracted objects with relative image error of 6.19%. For a large object (see the second column in Table 1), the EIT image was reconstructed by the sensitivity back-projection (SBP) algorithm, where only the object shape and location are revealed and the object size is hardly determined. The empirical value method and the histogram-based method (Otsu’s method) extract but underestimate the object in the segmented image, while the SPA method gives a better estimation of the extracted object with relative image error of 5.21%. The jagged boundary in the SPA result is a product of pixel error of image discretization as the setup cannot be represented by a series of pixels. For four small objects with the same size (see the third column in Table 1), the EIT image was reconstructed by the Landweber algorithm, which offers a good imaging performance with revealing approximate object shape, size, and location from human visual perception, but serious reconstruction artefacts exist. The empirical value method extracts approximate objects by removing the reconstruction artefacts, while the histogram-based method (Otsu’s method) mistakenly combines four small objects into a large object by determining a lower thresholding value. The SPA method gives a better estimation of the extracted four object with relative image error of 8.91%.
3.2. Universality
Since the SPA is a method for EIT image postprocessing, its universality is necessary to be investigated. In the simulation results in Table 1, different image reconstruction algorithms were employed to obtain the EIT image of each setup, and then three image thresholding methods were employed and compared on such EIT images. The empirical value method presents a challenge on its accuracy since the fixed empirical thresholding value is not adaptive to different setups and different algorithms. Although the thresholding value of the histogram-based method (Otsu’s method) is adaptive to different EIT images, its accuracy needs to be improved further when it is directly employed in EIT image segmentation. When it comes to the SPA method, the setups are approximately approached by segmenting EIT images with an optimal thresholding value, and even the EIT images vary from different setups to different image reconstruction algorithms.
3.3. Antinoise Capability
The antinoise capability of the size projection algorithm was modelled with the setup of four uniformly distributed small objects, and the results are given in Table 2. In the simulation, different levels of random noise were introduced into the boundary voltages to generate the noise-contained boundary dataset. The Landweber algorithm was employed to reconstruct the EIT images from three noise-contained boundary datasets, and then the size projection algorithm produced a sequence of segmented images from the reconstructed images. It appears to have a good antinoise capability when the level of random noise is less than 2.5%, while a serious distortion occurs if random noise level is more than 5%, which is primarily caused by the distortion of the reconstructed image. Therefore, the antinoise capability of the SPA depends to some extent on the quality of EIT reconstruction image.
3.4. Running Time
The time consumption of these image thresholding methods was tested by an Acer laptop with a 3.2 GHz A10-7300 processor and a 4G running memory. The programs and all related libraries (e.g., mesh generator, sensitivity matrix calculator, forward solver, and golden-section searcher) were coded and run in MATLAB 2016a software. Simulation results are given in Table 3. Based on the input of the EIT image and EIT-measured boundary voltage, the empirical value method directly gives the segmented image with a fixed time due to using a fixed empirical thresholding value, while the histogram-based method (Otsu’s method) and the SPA-based method need several searching steps to approach their “machine-determined optimal” segmented image. The speed of the SPA-based method is much slower than that of common methods, which is the product of solving forward solution in each searching step.
4. Measurement
4.1. Phantom Experiment
As shown in Table 4, a sequence of tests in phantom setups (see the first line in Table 4) were conducted using a P2000 EIT system [22] with 16-electrode sensor ring in 50 mm diameter in the On-Line Instrumentation Laboratory (OLIL) at the University of Leeds. The setups consist of nonconductive Plexiglas rods and 0.466 mS/cm tap water, which are a 15 mm diameter rod in the centre of tap water, three 10 mm diameter rods uniformly distributed in tap water, and two ellipse rods randomly distributed in tap water. The P2000 EIT system utilised adjacent sensing strategy to produce 104 independent boundary voltages, where a pair of 15 mA sinusoidal current terminals at frequency 9.6 kHz were employed. With the measured boundary voltages, the truncated singular value decomposition (TSVD) algorithm was employed to reconstruct EIT images (see the second line in Table 4), and then the EIT images were segmented to extract the rods by the SPA-based method, as listed in the third line in Table 4. The results of phantom experiment demonstrate that the SPA method can extract the objects in EIT images by removing the reconstruction artefacts effectively, and the extracted objects have approximate shape, size, and location in segmented images, which provides a method for quantitative measurement and object extraction using the EIT technique.
4.2. Gas-Water Two-Phase Flow Experiment
A gas-water two-phase flow experiment was conducted on the gas-water flow loop in OLIL at the University of Leeds, where the gas phase is compressed air (0 mS/cm) with 0.5 m/s superficial velocity and the water phase is tap water (0.35 mS/cm) with 0.51 m/s superficial velocity. The V5R [23] system with 312.5 frames per second (fps) sampling speed was utilised to collect the gas-water flow data, and the ITS P2+ software [22] was employed to produce the EIT images with a preloaded MSBP algorithm, then the EIT images were segmented by the SPA machine-determined thresholding value. The experimental schematic is illustrated in Figure 4. The experimental results of gas-water flow are shown in Figure 5, where a gas slug bubble is flowing through the EIT sensor and is being recorded by 90-frame discrete EIT cross-sectional images. At the time of V5R collecting the 1580th image, the slug bubble starts flowing through the EIT sensor. The bullet head of the slug bubble flows through during the time of 1580∼1600 images, and the slug body flows through during the time of 1600∼1640 images. After the tail of the slug bubble (1640∼1660 images) flowing through, the whole slug bubble passes through the EIT sensor. In EIT images (see Figure 5(a)), the boundary of the slug bubble cannot be clearly shown, and the size of the slug bubble is hardly identified. However, the boundary and size of the slug bubble is extracted by the SPA-based image thresholding (see Figure 5(b)), which provides a quantitative method to measure slug bubbles in water.
(a)
(b)
5. Conclusions
Focusing on the thresholding value selection in tomographic image postprocessing, a multistep image segmentation process, size projection algorithm (SPA), is proposed to extract the reconstructed objects in EIT images. This method determines an optimal thresholding value for EIT image segmentation by comparing the computed boundary voltages corresponding to each segmentation image with the measured boundary voltages (i.e., minimising projection error). Results of simulation and experiments demonstrated that the performance of SPA-based image segmentation is better than other methods. It can be used to measure and visualise the large gas bubbles in gas-water two-phase flow, and its application also can be extended to other industrial and medical applications.
However, there are still a few aspects to be addressed in future studies. The SPA-based image thresholding is not very effective in time consumption, which needs to be improved by applying an effective forward solution solver to reduce the time consumption in each searching step or by employing advanced optimisation measures to reduce the searching steps. Since the multistep image segmentation is dependent on the quality of the imported EIT image and boundary data, i.e., the quality of tomographic reconstruction and measurement error, the resultant results could be improved by applying advanced tomographic algorithms or high-performance EIT systems. It is also worthwhile to point out that the proposed method has limited capability to extract too small object as its induced boundary voltages are too weak to be identified by the commercial EIT system. This work also can be extended to other tomographic applications for two immiscible phase measurement and visualisation.
Appendix
A. Numerical Forward Solver
The finite element method (FEM) is an efficient method to solve EIT forward problem, where the solution of measured region is approached by discretizing the measured region into finite elements and calculating the approximate solution in each element. A FEM model presenting a 2D cross section of a measured region can be solved by equation (A.1) [24]:where is the admittance matrix corresponding to the conductivity distribution vector of the segmented image , is the nodal voltage vector for obtaining the computed boundary data , and is the nodal current vector. Since the admittance matrix has been proven as a symmetrical and positively defined matrix [25], equation (A.1) can be solved by the conjugate gradient (CG) method. The CG method utilises the idea of minimising function by up to N-time searching in conjugate directions. The forward solution could be accomplished by the following procedures.(1)Calculate the initial nodal voltage vector and the initial searching direction(2)Carry out the iterations from to until the residual is sufficiently small(3)When the approximate solution of nodal voltage vector is reached, the computed boundary voltages corresponding to the segmented image can be obtained by following a specific sensing strategy (adjacent sensing strategy is employed in our case, i.e., the forward solution of 104 independent voltages).
Data Availability
The experimental data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest or state.
Acknowledgments
The Chinese Scholarship Council (CSC) and the School of Chemical and Process Engineering are acknowledged for supporting Mr. Kun Li’s visit to the University of Leeds. The authors would also like to thank Professor Mi Wang and Dr. Qiang Wang at the University of Leeds for offering helpful ideas in this work.