A quantitative evaluation of pleural effusion on computed tomography scans using B-spline and local clustering level set

Abstract

Estimation of the pleural effusion’s volume is an important clinical issue. The existing methods cannot assess it accurately when there is large volume of liquid in the pleural cavity and/or the patient has some other disease (e.g. pneumonia). In order to help solve this issue, the objective of this study is to develop and test a novel algorithm using B-spline and local clustering level set method jointly, namely BLL. The BLL algorithm was applied to a dataset involving 27 pleural effusions detected on chest CT examination of 18 adult patients with the presence of free pleural effusion. Study results showed that average volumes of pleural effusion computed using the BLL algorithm and assessed manually by the physicians were 586 ml±339 ml and 604±352 ml, respectively. For the same patient, the volume of the pleural effusion, segmented semi-automatically, was 101.8% ±4.6% of that was segmented manually. Dice similarity was found to be 0.917±0.031. The study demonstrated feasibility of applying the new BLL algorithm to accurately measure the volume of pleural effusion.

Keywords

CT pleural effusion volume B-spline local clustering level set

1 Introduction

Pleural fluid is physiologically present in the human’s pleural cavity [1, 2]. In normal conditions, the unilateral volume is 8.4±4.3 ml [3]. As the fluid amount increasing, the lung compliance reduces due to the increased lung capillary pressure. It may lead to patient dyspnea [4]. However, pleurocentesis provides a relief from dyspnea. The estimation of the pleural effusion’s volume may directly affect the therapeutic evolution of the disease and supports the decision whether to perform pleurocentesis or not [5]. Meanwhile, the medical treatment will be adjusted based on the change of pleural effusion volume. Therefore, accurate calculation of pleural effusion volume plays an important role in patient management.

There are several approaches to detect and evaluate the volume of pleural effusion. A typical method is posteroanterior chest radiograph. Pleural effusion is visible on the posteroanterior projection as a meniscus at about 200 ml. The hemidiaphragm will obliterate at approximately 500 ml. Effusion which is lower than 200 ml can be detected only when a lateral chest radiograph is taken. Pleural effusion becomes visible as a meniscus in the posterior costophrenic sulcus at approximately 50 ml [5].

Sonography is a good choice for detecting pleural effusions in supine patients, especially for the patient in intensive care unit. The volume of pleural effusion can be quantified by the equation y = 47.6×x-837, where y is the predicted pleural effusion volume in milliliters and x is the maximal thickness of the pleural effusion lamella in millimetres [6].

Computed tomography (CT) is the most suitable method for detection and evaluation of pleural effusion [7 –10]. On CT scans, the effusion dimensions can be obtained easily and the volume of pleural effusion can be evaluated using a variety of methods. Mergo. PJ, Helmberger. T et al. assessed the volume of pleural effusion by the equation V = d² × l, where d is the greatest depth of the effusion on a single CT scan, l is the greatest length of the effusion [8]. Hazlinger. M, Ctvrtlk. F, et al. devised a method with single measurement, the greatest depth of effusion, to evaluate the volume of effusion. The volume can be quantified by their recommended tabular [10].

However, all of the methods may be not accurate enough when the fluid amount is over 300 ml. Large deviation may be occurred under this condition. For example, in M. Hazlinger’s research [10], though the volume of pleural effusion is approximately 900 ml in a patient, the predicted volume is only about 600 ml using the recommended equation. This evaluation may lead the inappropriate treatment. The major limitations of these methods are: firstly, there is an implicit supposition on regression equations. The shape of pleural effusion has geometric similarity. For example, the shape of pleural effusion was often described as meniscus while the volume of pleural effusion can be evaluated by several measurements. However, this supposition may be inappropriate in some situations. The shape of pleural effusion may deform irregularly and cannot be described as meniscus. Thus, the dispersion of the dataset will be increased. The unavoidable errors may be occurred. This deviation would be significant while the measurement of the equations is limited. Secondly, the deviation will be induced when the physician(s) get the equation measurement(s). The shape of pleural effusion may be irregular in some conditions. It is difficult to obtain the measurements precisely. As the result, the deviation will occur. It will be enlarged while there are more measurements required.

Above problems can be avoided using the image segmentation technology [11 –15]. There are several prior efforts to segment and measure the volume of pleural effusion in CT scans [16 –21]. In most of these segmentation methods, the region in inner chest wall and diaphragm was countered in the first step and some of the organs and tissue contained in this region (e.g. lung, heart, etc.) was excluded based on the gray-level of the CT image. After that, the remaining section was the area of the pleural effusion. These methods can be named as “excluding segmentation methods”. However, these methods may be not accurate enough under the following conditions:

First, there is large volume (e.g. over 900 ml) of pleural effusion in the pleural cavity. The shape of the lung will deformed obviously. In some of the CT scans which near the diaphragm, all the air in the lung will be compressed out by the pleural effusion. The gray-level of the lung will be changed significantly. As a result, the boundary of lung and diaphragm will be delineated improperly. The volume of the pleural effusion will be calculated imprecisely, especially when the boundary of lung is an important input parameter for the method.

Second, the patient may have some other disease (e.g. pneumonia) at the same time. The gray-level of the pneumonia is much higher than lung and approximate to the effusion. The region of the pneumonia may be separated as a part of the effusion. Then the errors will be occurred in the calculations.

This study aims to develop and test a novel semi-automatic segmentation algorithm based on B-spline and local clustering level set (LCLS) method, namely BLL, to evaluate the volume of pleural effusion. Comparing with the “excluding segmentation methods”, BLL will focus on the area of pleural effusion and delineate the boundary directly. The above limitations are removed in the present study. The procedure of BLL is as follows:

First, the B-spline method generates the inner chest wall boundary by the interior points of the thoracic cage image. Then images in the boundary were extracted. Based on these images, the LCLS method segments the pleural effusion region into four phases [22]. The region of pleural effusion is assembled by the novel criterions that include analyzing the intensity, location and the shape of the image within each phase. Finally, the pleural effusion volume is calculated using the CT scans parameters.

2 Materials and methods

2.1 Patient population and device

A dataset involving 27 free pleural effusions detected by physicians from 18 adult patients on chest CT images during the period January 2015 to December 2015 was used in this study. There were 11 men and 7 women (mean age: 51 years, age range: 18–75 years). 9 patients were with unilateral pleural effusions (4 right-sided, 5 left-sided) and another 9 patients were with bilateral. Philips Brilliance CT Big Bore Systems (Philips Medical System, Cleveland, OH) performed the CT examinations with 64 detector rows. The tube voltage and tube current were fixed as 120 keV and 200mAs, respectively. The thickness of scans was 5 mm. Meanwhile, the resolution was 512×512 pixel, almost 1mm×1 mm for each pixel, roughly.

2.2 Extraction of the image within the inner chest wall

Thoracic cage consists of ribs, sternum, and thoracic vertebrae. Suppose that the inner chest wall boundary can be delineated by the interior points of the thoracic cage. Then the images within the boundary can be obtained. There are three major steps for this purpose:

Extract the interior points of thoracic cage from the CT scans.

Delineate the inner chest wall boundary based on the thoracic cage images.

Extract the images within the inner chest wall boundary.

2.2.1 Extraction of interior points of thoracic cage

In a patient’s CT scan, the effective area of patient can be countered by a cubic, with the help of three parameters, height (H), width (W) and length (L). We can get the scans and represent the image within the cubic as I = {I_s|s = 1, 2,..., N} (N is the last scan number). Figure 1 shows the patient CT scans.

Fig.1

The patient CT scans.

Because the gray-level of the bones is much higher than the others, almost 1400 HU in average, the images of bones $I_{B}^{*}$ can be easily extracted from the original CT scans. However, there are some unexpected bones (e.g. humeri, shoulder blade, etc.) images I_uB = {I_uB,s|s = 1, 2,..., N} presented in the images and should be deleted. On the other hand, during the chest CT examination, more often practiced poster of patient, with his/her hands settled behind his/her head, is used. Therefore, humeri and shoulder blade are only in the top slices and backside of the scans. Then, the regional growth (RG) method is employed to remove these bones’ image.

The RG method is used twice from two directions and with different mask, I_mask1 and I_mask2, respectively. In the top to bottom processing, I_mask1 is with two triangle areas in the sides of the image, as shown in Fig. 2(a). The humeri and shoulder blade image in the first scan $I_{uB, 1}^{+}$ is found by this mask. Setting $I_{uB, 1}^{+}$ as origin seed image, repeat the two steps:

Use the ith seed image and the RG method to find the intersection image $I_{uB, i + 1}^{+}$ in the i + 1th slice.

Set $I_{uB, i + 1}^{+}$ as new seed image.

Fig.2

The green areas in (a) and (b) are I_mask1 and I_mask1 for the RG method, respectively. (a) h is the distance between the top of the clavicle and the top line of the scan.

The humeri and shoulder blade images $I_{uB}^{+}$ are selected. However, some of the ribs image are also in $I_{uB}^{+}$ . Then the RG method uses again to save the ribs image. This time, a rectangle area is set in the top of I_mask2, as shown in Fig. 2 (b). The humeri and shoulder blade image $I_{uB}^{-}$ is extracted in following processing. The union of last chest CT scan I_N and I_mask2 set as origin seed image. Repeating the following steps:

Use the ith seed image and the RG method to find the intersection image $I_{B, i - 1}^{-}$ and set $I_{B, i - 1}^{-}$ as new seed image.

$I_{uB, i - 1}^{-} = I_{i - 1} - I_{B, i - 1}^{-}$ .

It is worth nothing that the shoulder blades are contacted with the ribs in the middle of shoulder blades only in few of the CT scan. As a result, no more than three slices in $I_{uB}^{+}$ and $I_{uB}^{-}$ are the same. Comparing with $I_{uB}^{+}$ and $I_{uB}^{-}$ , the contacted slices I_BS = {I_BS,s|s = k₁, ..., k₂} (k₁ is the first same scan number, k₂ is the last same scan number) are obtained. Then thoracic cage image I_B is obtained following equations (1) and (2), as shown in Fig. 3. $I_{B, i}^{*} = I_{B, i}^{*} - I_{B, i}^{+} k = 1, . . ., k_{1} - 1$ (1) $I_{B, i} = I_{B, i}^{*} - I_{B, i}^{-} k = k_{2}, . . ., N$ (2)

Fig.3

The thoracic image.

If the thorax is divided into two parts, namely the left one and the right one, the centroids can be set as P_CL = {P_CL,s|s = 1, 2, ..., N} for the left, and P_CR = {P_CR,s|s = 1, 2, ..., N} for the right. It is assumed that the left centroid in the nth CT scan is P_CL,n, and there are M bones images in one CT scan. Setting P_CL,n as the origin point, the coordinate system can be established, as shown in Fig. 4. The boundary point of the mth bone is written as z_n,m(d, θ), where d is the distance between the boundary point, θ is angle of the point. Suppose that, the boundary points in the qth angle section [θ_q, θ_q + dθ] are $Z_{n, m}^{q}$ . Then the boundary point $Q_{n, m}^{q}$ , which is closest to the centroid P_CL,n, is selected. Repeat it for all the angle sections, the inner chest wall boundary points in the nth CT scan Q_B,n are obtained, as shown in Fig. 5.

Fig.4

Schematic representation of extracting interior points of thoracic cage. (a) shows the established coordinate system. (b) shows the detail of selecting process. $Q_{n, m}^{q}$ is closest to the centroid P_CL,n in angle section dθ.

Fig.5

The extracted bone boundary points.

2.2.2 Delineation of the inner chest wall boundary

After Q_B is obtained, the thorax boundary C_Th can be approximated by the B-spline surface [23, 24]. It is widely used in medical field [25 –27]. The B-spline surface S(u,v) can be expressed as: $S (u, v) = \sum_{i = 0}^{n} \sum_{j = 0}^{m} N_{i, p} (u) N_{j, q} (v) P_{i, j}$ (3)

Where P_i,j is control point. m and n are the number of the control points in u and v directions, respectively. N_i,p and N_j,q are basic functions of order p and q, which with the knot vectors U and V. N_i,p can be expressed as: $\begin{matrix} N_{i, 0} = {\begin{matrix} 1, u_{i} \leq u < u_{i + 1} \\ 0, others \end{matrix} \\ N_{i, p} = \frac{u - u_{i}}{u_{i + p} - u_{i}} N_{i, p - 1} (u) + \frac{u_{i + p + 1} - u}{u_{i + p + 1} - u_{i + 1}} N_{i + 1, p - 1} (u) \end{matrix}$ (4)

The knot vector U and V can be expressed as: $U = {\underset{p + 1}{0, . . ., 0,} u_{p + 1}, . . ., u_{r - p - 1}, \underset{p + 1}{1, . . ., 1}}$ (5) $V = {\underset{q + 1}{0, . . ., 0,} u_{q + 1}, . . ., u_{s - q - 1}, \underset{q + 1}{1, . . ., 1}}$ (6) where the number of the knots in U and V are r = n + p + 1 and s = m + q + 1, respectively.

The surface approximation process can be implemented in two steps:

The inner chest wall boundary in each CT scan is approximated by the control point $P^{*} = {P_{s, i}^{*} | s = 1, 2, . . ., N, i = 1, 2, . . ., n + 1}$ .

P^* is set as data point, P, the control point of S(u,v), is obtained. Substituting P into (1), S(u,v) is obtained.

In the first step, the weighted least squares curve approximate Q_B,i in ith scan. The data points are separated into two parts: constraint data $Q_{B, i}^{c}$ and unconstraint data $Q_{B, i}^{u}$ . The constraint data will be exactly approximated. The unconstraint data may not. The weight of each unconstraint point is ω_k (define as ω_k = 1). Then the approximating problem is solved as a constraint minimum processing. The Lagrange multiplier method is employed to solve it. Lagrange multipliers, A = [λ_k], are used in the method, The unconstraint data is described as $N P_{i} = S$ (7) where N is the basic function of the unconstraint data point, P_i is control point, S is unconstraint data $Q_{B, i}^{u}$ .

The constraint data is described as $M P_{i} = T$ (8) where M is the basic function of the constraint data point, T is constraint data $Q_{B, i}^{c}$ .

The error vector of the unconstraint data is S– N P_i. The error under the condition of M P_i = T will be minimized. According to the Lagrange multiplier method, the following function should be minimized. $(S^{T} - P_{i}^{T} N^{T}) W (S - N P_{i}) + A^{T} (M P_{i} - T)$ (9) where W is a diagonal matrix with the element {ω_k}.

If take the derivative of (9), and define it equal to zero, P_i is shown as $P_{i} = (N^{T} WN)^{- 1} N^{T} W S - ((N^{T} WN)^{- 1}) M^{T} A$ (10)

Premultiplying (10) with M , A is solved. $A = (M (N^{T} WN)^{- 1} M^{T})^{- 1} (M (N^{T} WN)^{- 1} N^{T} W S - T)$ (11)

Substituting A to (10), P_i is obtained.

In this study, we used 13 control points to approximate the inner chest wall boundary in each scan, as shown in Fig. 6.

Fig.6

The distribution of the control points in a CT scan. The last control is the same as the first one.

The second step of approximation the surface is similar with the first one. The rth points in all scans P_r = {P_i,r|i = 1, 2,..., N} are set as the new data point Q_r. Q_r is approximated by m + 1 control points P_j. All of the data points are defined as unconstraint data Q^u. The origin weight of these points is $ω_{k}^{(1)} = 1$ . The process of approximation is as follows:

According to (9, 10), the control point $P_{j}^{(t)}$ is obtained under the condition of $ω_{k}^{(t)}$ . Then the data point Q^(t) is solved.

Calculate the distance $d_{k}^{(t)}$ between Q^(t) and Q, $d_{k}^{(t)} = | | Q_{k}^{(t)} - Q_{k} | |$ . Find the maximum distance $d_{\max}^{(t)} = max (d_{k}^{(t)})$ .

Determine whether $d_{\max}^{(t)}$ meet convergence condition. Output $P_{j}^{(t)}$ if it is satisfied. Otherwise, generate $ω_{k}^{(t + 1)}$ and repeat step (1-2).

The convergence condition is set as $d_{\max}^{(t)} \leq d$ , where d is the dimension of 4 pixels. $ω_{k}^{(t + 1)}$ is generated as follow. $ω_{k}^{(t + 1)} = {\begin{matrix} 1, d / d_{k}^{(t)} \geq 1 \\ d / d_{k}^{(t)}, others \end{matrix}$ (12)

After all, the control points P of surface S (u, v) are obtained. Then the inner chest wall boundary in each scan C_Th is calculated. Image I_Th, which contained within the boundary, can be extracted. Figure 7 shows the thorax boundaries in different scans.

Fig.7

The inner wall chest boundaries in different CT scans.

2.3 Segmentation of the pleural effusion

As described above, the image within inner chest wall boundary I_Th is obtained (as shown in Fig. 8). However, many tissues and organs, heart, lung, liver, artery vessels, etc., are still in I_Th. Furthermore, the distribution, component and volume of pleural effusion are quite different from one patient to the other. Thus, the LCLS method is employed to segment the pleural effusion image. Two assumptions about the CT scans are introduced as follows:

The CT scans has suitable signal to noise ratio. It is not necessary to reduce the image noise.

The gray-scale of pleural effusion is approximate for same patient.

Fig.8

The images with the inner chest wall boundary in different CT scans.

In this paper, we consider the CT image is intensity inhomogeneous. I_Th is divided into multiphase. Each phase has intensity constant and contains several image regions. All the regions are assessed by their intensity, location, and the geometric information. Hence, the effusion region is assembled by the part of the selected region. The procedures of the algorithm are as follows:

Segment I_Th by the image intensity, preliminarily. Then the region of interest (ROI) is obtained.

Divide ROI into four phases using the LCLS method [22].

Select and assemble the regions in all phases. Then, extract the region of effusion.

2.3.1 Extraction of the ROI

The ROI is obtained in two steps:

I_Th is segmented by their intensity preliminarily. The region will be removed if their intensity is quite different from the effusion.

the intensity of image is enhanced for the following segmentation.

In the ith scan I_Th,i, Physician(s) will delineate one or several sub-regions of effusion in only one CT scans, namely Ω_m = {Ω_m,j|j = 1, 2,..., n}, where n is the last number of the sub-region (shown as Fig. 9a). The delineation can be easily finished in few second. Let Ω_M = ⋃ _iΩ_m, i, the mean gray-scale Y_m and standard deviation Y_std will be calculated among the points in Ω_M. Because the component of effusion is almost same, the mean gray-scale Y_m and standard deviation Y_std will change a little when the different Ω_M is obtained. Then extracting the points which gray-scale is in the range of [Y_m - 100, Y_m + 100], the origin ROI $I_{ROI}^{*} = {I_{ROI, i}^{*} | i = 1, 2, . . . N}$ is obtained. To separate the tissue region effectively, the ROI image was piecewise enhanced as the follow.

Fig.9

ROI of pleural effusion. (a) shows image within the inner chest wall. The ellipses are the sub-regions delineated by the physician(s). (b) shows the image which enhanced by (14).

$Y_{n} = {\begin{matrix} Y, Y \leq Y_{m} \\ Y_{m} + C \times (Y - Y_{m}), Y > Y_{m} \end{matrix}$ (13) where Y is origin gray-scale, Y_n is new gray-scale after enhanced, C is a constant, C = 2.5. Then, the ROIs in all scans I_ROI = {I_ROI,i|i = 1, 2,... N} are obtained, shown as Fig. 9b. The new mean gray-scale $Y_{m}^{*}$ and standard deviation $Y_{std}^{*}$ will be calculated among the same points in ROI.

2.3.2 Segmentation of the ROI

We hope that each ROI can be segmented into four kinds of regions. The first one’s gray-scale is greater than the pleural effusion’s. The second one’s gray-scale is approximated. The third one’s gray-scale is less than the effusion. And the last one is the background, in which gray-scale should be close to zero. As the result, the tissues and organs around the effusion (e.g. pneumonia) are distinguished by these regions. The LCLS method is employed to achieve it.

2.3.2.1 LCLS formulation

As aforementioned, the CT image is intensity inhomogeneous. The observed image I is defined on a continuous domain Ω, and modeled as $I = bJ + n$ (14) where J is the true image measures an intrinsic physical property of the objects being imaged, b is the component that accounts for the intensity inhomogeneity and assumed to be slowly varying, n is assumed to be zero-mean Gaussian noise image noise.

We believed that, the true image J can be separated into disjoint regions Ω = Ω₁, …, Ω_N by N distinct constant values c = (c₁, …, c_N), respectively. The segmentation process is achieved by minimizing the energy function $F = E + ν L + R_{p}$ (15) where E is the data term, which forces the constant to be close to the image. νL is the length regularization term, which serves to smooth the contour. R_p is signed distance regularization term, which forces the constant to be smooth in the separated regions.

In (15), the data term is of the most important. The local intensity cluster property is introduced in this term. The standard K-means clustering is allowed to classify these local intensities, which is defined as $E_{y} = \sum_{i = 1}^{N} \int_{Ω_{i}} K (y - x) | I (x) - b (y) c_{i} |^{2} dx$ (16) where K (y–x) is a Gaussian function with standard deviation σ. b (y) is varying bias field. The smaller the value of E_y, the better classification. K (u) is defined as $K (u) = {\begin{matrix} \frac{1}{a} e^{- | u |^{2} / 2 σ^{2}}, | u | \leq ρ \\ 0, others \end{matrix}$ (17) where a is a normalization constant such that ∫K (u) = 1. ρ is the radius of classification of the intensities neighborhood. Note that, if the image is more localized intensity inhomogeneity, the bias field b varies faster, thus the smaller ρ should be used as the radius.

Minimizing E_y for all y in Ω, the energy function of the data term in (16) is defined as E = ∫E_ydy. In order to find the solution of the energy E, E = ∫E_ydy convert to a level set formulation. In level set methods, the disjoint regions Ω are partitioned by the level set functions Φ = (φ₁, ..., φ_k). The energy E is written as a multiphase level set formulation $E (Φ, c, b) = \int \sum_{i = 1}^{N} e_{i} (\times) M_{i} (Φ (x)) d x$ (18) where M_i(Φ) is the membership functions, e_i is defined by $e_{i} (x) = \int K (y - x) | I (x) - b (y) c_{i} |^{2} d y$ (19)

The regularization terms in (15) is rewritten as the level set function, $L (Φ) = \sum_{j = 1}^{k} L (φ_{j})$ and $R_{p} (Φ) = \sum_{j = 1}^{k} R_{p} (φ_{j})$ , where L(φ_j) is defined by $L (φ) = \int | \nabla H (φ) | d x$ (20) which computes the arc length of the zero level contour of φ and serves to smooth the contour by penalizing its arc length [28, 29]. The energy term R_p(φ_j) is defined by $R_{p} (φ) = \int p | \nabla φ | d x$ (21) where function p is defined by p (s) = (1/ - 2) (s - 1) ². The energy R_p(φ) will be minimized when | ∇ φ|=1.

2.3.2.2 Energy minimization

Then the energy function F in multiphase level set formulation is defined by $F (Φ, c, b) = E (Φ, c, b) + ν L (Φ) + R_{p} (Φ)$ (22)

The minimization of the energy F (Φ, c, b) is achieved by an iterative process: in each iteration, the energy F (Φ, c, b) is minimized with respect to each of the variables Φ, c and b, given the other uploaded in previous iteration. The solution to the energy minimization as follows.

1) Energy minimization with respect to Φ: For fixed c and b, the minimization of F (Φ, c, b) with respect to Φ is performed solving the following gradient flow equations: $\frac{\partial Φ}{\partial t} = - \frac{\partial F}{\partial Φ}$ (23) where ∂F/∂Φ is the Gâteaux derivative [30] of energy F. By calculus of variations, ∂F/∂Φ is expressed as the following gradient flow equations: $\begin{matrix} \frac{\partial φ_{1}}{\partial t} = - \sum_{i = 1}^{N} \frac{\partial M_{i} (φ)}{\partial φ_{1}} e_{i} + ν δ (φ_{1}) div (\frac{\nabla φ_{1}}{| \nabla φ_{1} |}) + μ div (d_{p} (| \nabla φ_{1} |) \nabla φ_{1}) \\ ⋮ \\ \frac{\partial φ_{k}}{\partial t} = - \sum_{i = 1}^{N} \frac{\partial M_{i} (φ)}{\partial φ_{k}} e_{i} + ν δ (φ_{k}) div (\frac{\nabla φ_{k}}{| \nabla φ_{k} |}) + μ div (d_{p} (| \nabla φ_{k} |) \nabla φ_{k}) \end{matrix}$ (24) where ∇ is the gradient operator, div(·) is the divergence operator, and the function d_p is defined as d_p = p′(s)/s. The variable Φ can be uploaded by solving (24)

2) Energy minimization with respect to c: For fixed Φ and b, the constant vector c is minimized by the energy E (Φ, c, b), denoted by $\hat{c} = ({\hat{c}}_{i}, . . ., {\hat{c}}_{N})$ , is given by ${\hat{c}}_{i} = \frac{\int (b * K) {Iu}_{i} d y}{\int (b^{2} * K) {Iu}_{i} d y}, i = 1, . . ., N$ (25) where u_i is the membership function, u_i(y) = M_i (Φ (y)).

3) Energy minimization with respect to b: For fixed Φ and c, the bias field b is minimized by the energy E (Φ, c, b), denoted by $\hat{b}$ , is given by $\hat{b} = \frac{({IJ}^{(1)}) * K}{J^{(2)} * K}$ (26) where $J^{(1)} = \sum_{i = 1}^{N} c_{i} u_{i}, J^{(2)} = \sum_{i = 1}^{N} c_{i}^{2} u_{i}$ .

2.3.2.3 Numerical implementation

In this study, the ROI is separated into four disjoint regions Ω₁, …, Ω₄ with the intensity constant c = c₁, … c₄. As described in Sec. 2.3.2, the intensity constant c is constrained as follow: $c_{1} \geq Y_{m}^{*} + 4 \times Y_{std}^{*}, c_{2} = Y_{m}^{*}, c_{3} \leq Y_{m}^{*} - 3 \times Y_{std}^{*}$ . And the initial value are $c_{1} = Y_{m}^{*} + 4 \times Y_{std}^{*}, c_{2} = Y_{m}^{*}, c_{3} = Y_{m}^{*} - 3 \times Y_{std}^{*}, c_{4} = 0$ . As a result, one should use at least two level set functions φ₁, φ₂ to define the membership functions M_i of the regions Ω_1,..., Ω₄, such that $M_{i} (φ_{1} (y), . . ., φ_{k} (y)) = {\begin{matrix} 1, y \in Ω_{i} \\ 0, else \end{matrix}$ (27)

For the four-regions case, M_i is defined as M₁ (φ₁, φ₂) = H (φ₁) H (φ₂), M₂ (φ₁, φ₂) = H (φ₁) (1 - H (φ₂)), M₃ (φ₁, φ₂) = (1 - H (φ₁)) H (φ₂), and M₄ (φ₁, φ₂) = (1 - H (φ₁)) (1 - H (φ₂)).

The level set evolution in (24) is implemented using the finite difference scheme as for the distance regularized level set evolution formulation provide in the previous study [31]. The Heaviside function H is replaced by a smoothed Heaviside function H_ɛ, which is defined by $H_{ɛ} (x) = \frac{1}{2} [1 + \frac{1}{π} arctan (\frac{x}{ɛ})]$ (28) where ɛ = 1. The dirac delta function δ is the derivative of H_ɛ, which is defined by $δ_{ɛ} (x) = H_{ɛ}^{'} (x) = \frac{1}{π} \frac{ɛ}{ɛ^{2} + x^{2}}$ (29)

The parameters in the level set function are chose as follow: the parameters μ and the time step Δt can be fixed as μ = 0.5 and Δt = 0.1. The parameter ν is set to 0.0002×255². The parameter σ is set to 4. And the kernel K is constructed as a ω×ω mask, with ω being the smallest odd number such ω≥4*σ+1. When σ = 4, the mask size is 17×17. Fig. 10 shows the regions segmented by the level set method.

Fig.10

The regions segmented by the level set method.

2.3.3 Selecting and assembling criterion

In Sec.2.3.2, the ROI is divided into four regions. Each region contains several sub-regions, as shown in Fig. 11a. Then we select the sub-regions for assembling the pleural effusion region by the following experiential criterion:

Fig.11

The result of the selecting process in each step. (a) is the image region segmented by LCLS method. (b) is the chosen image sub-region which adjoined the lower boundary of inner chest wall. (c) shows the assembled region and the cross in (c) are the separated points. (d) shows the final region, which will be obtained.

choose one sub-region Ω_S_,0, which intensity constant c is close to the mean gray-scale of the pleural effusion Y_m, from the sub-regions adjoined the lower boundary of the inner chest wall, as shown in Fig. 11b.

find all sub-regions adjacent to Ω_S_,0, denoted by Ω_S_,1,...,Ω_S_,n, where n is last number of the adjacent sub-region. The mean gray-scale of Ω_S_,i, i = 1,...,n is calculated, denoted by Y_s,1, … Y_s,n. The Ω_S_,i is selected to assemble the pleural effusion area if its gray-scale is in the range of [Y_m–3·Y_std, Y_m–2.5·Y_std]. Then the origin area of pleural effusion $Ω_{E}^{*}$ and the boundary of the area $C_{E}^{*}$ are obtained, as shown in Fig. 11c.

$Ω_{E}^{*}$ is separated based on the distance between the boundary points. i.e. if the distance between two points P_i and P_j is less than 5 pixel, the image will be divided into two parts by the line between two points. The bottom area of $Ω_{E}^{*}$ will be selected.

Two separated lines are used to divide the effusion area. We define the line with the top point of the spinal column as base line. Two separated lines are defined upon base line 15 and 30 pixels. Both of the lines are half wide of the image, as shown in Fig. 11c. From line 1, if the image can be divided into two parts by it, the upper part will be eliminated and the lower part will be selected as the new $Ω_{E}^{*}$ . Otherwise, line 2 will repeat this step, and $Ω_{E}^{*}$ will be uploaded.

If there were holes in the image domain, the mean gray-scale of the image within the holes are evaluated. And the hole will be filled while its intensity is in the range of [Y_m–3·Y_std, Y_m–2.5·Y_std].

After all, the area of pleural effusion Ω_E can be obtained, as shown in Fig. 11d.

Counting the number of pixel in Ω_E, the volume of pleural effusion is quantified by the parameters of the scans. Figure 12 shows the segmentation of the pleural effusion.

Fig.12

The segmentation of pleural effusion.

3 Results

An experienced radiologist manually contoured the pleural effusions on each CT scan using the Philips EBW 4.5.3.40140 image workstation. Based on the manual boundaries, volume of pleural effusion was calculated in milliliters and the resulted range of the effusion volume was 64 ml to 1475 ml, whereas the effusion volume (mean±standard) was 586 ml±339 ml. Meanwhile, the regression equation used here was designed by Martin Hazlinger, Filip Ctvrtlik, et al. [10], assessed the effusion volume, the effusion volume (mean±standard) was about 646 ml±455 ml. (Table 1 showed the conversing tabular recommended by Martin Hazlinger’s method. One may quantify the volume of pleural effusion basing on the thickness of the effusion, shown as in Fig. 13). Using BLL, the effusion volume (mean±standard) was 604 ml±352 ml. Comparing the results from the assessment, the volume of the pleural effusion segmented by our method was 101.8% ±4.6% of that segmented manually. The maximum error was 11.2% and the minimum error was –6.68%. The Dice coefficient of the segmentation between BLL and the manual one was 0.917±0.031. The volume of pleural effusion segmented by the regression equation was 118.9% ±94.4% of that segmented manually. The maximum error was 462.6%, and the minimum error was –35.3%. The details of the results are shown in Table 2.

Table 1
Conversion of pleural effusion’s thickness to volume*

Thickness d (cm) Pleural effusion volume (ml)

1 67

1.5 142

2 216

2.5 288

3 360

3.5 431

4 501

4.5 572

5 643

5.5 714

6 786

6.5 859

7 933

7.5 1009

8 1086

8.5 1166

9 1248

9.5 1333

10 1421

Thickness d (cm)	Pleural effusion volume (ml)
1	67
1.5	142
2	216
2.5	288
3	360
3.5	431
4	501
4.5	572
5	643
5.5	714
6	786
6.5	859
7	933
7.5	1009
8	1086
8.5	1166
9	1248
9.5	1333
10	1421

*: Martin Hazlinger, Filip Ctvrtlik, et al. [10].

Fig.13

Assessment procedures designed by Hazlinger et al. d is the greatest depth of the effusion on a single CT scan.

Table 2

The results of the pleural effusion volume

Effusion number	Radiologist (ml)	Regression equation (ml)	Our algorithm (ml)	Dice coefficient
1	480	459	506	0.925
2	1280	–^*	1352	0.917
3	760	564	740	0.958
4	453	374	506	0.855
5	825	933	848	0.932
6	269	281	284	0.906
7	193	127	209	0.890
8	1475	–^*	1379	0.882
9	492	572	506	0.939
10	64	360	68	0.921
11	456	501	442	0.935
12	591	459	586	0.951
13	337	259	324	0.943
14	1047	1117	1078	0.947
15	370	346	389	0.926
16	792	757	871	0.844
17	133	157	142	0.910
18	695	757	755	0.867
19	552	600	571	0.941
20	565	402	550	0.947
21	228	316	214	0.903
22	324	430	302	0.893
23	320	302	309	0.926
24	867	844	864	0.963
25	1077	1150	1155	0.890
26	791	700	752	0.917
27	637	845	618	0.934

^*: The greatest depth of the effusion is out of the evaluation range of recommended tabular. The effusion volume cannot be assessed.

The computer with Intel(R) Core(TM) i7-2600 was used for the calculations. The frequency of the CPU is 3.40 GHz and the memory is 8 G. It would spend no more than 8 second for calculating one CT scan. For one patient with 45 scans, the whole processing will be finished in almost 6 minutes.

4 Discusson

The analytical results showed that the algorithm proposed in this study is a feasible way to assess the volume of pleural effusion. The algorithm can assess the volume of the pleural effusion effectively and precisely. In most of the tests, the volume error between BLL and the manual one is less than 10%.

In the published papers, the regression equations may imprecise because of insufficient sample space and the dispersion of the samples. Meanwhile, the pleural effusion images vary from patient to patient, the shape of effusion deformed irregularly (shown in Fig. 14) which is difficult to extract precise measurements for the equations. Therefore, additional error may occur. The existing (semi-)automatic methods also could not segment the image accurately, especially when the fluid amount is large and/or the patient having some other disease (e.g. pneumonia) at the same time. Thus, it is necessary to develop the new quantification algorithm for the pleural effusion calculation.

Fig.14

Images of the pleural effusion. The images of pleural effusion vary from patient to patient. It is then difficult to extract the measurement precisely. However, the algorithm in this paper can delineate their boundaries to evaluate the volume.

The images of pleural effusion vary from patient to patient. It is then difficult to extract the measurement precisely. However, the algorithm in this paper can delineate their boundaries to evaluate the volume.

In this study, we developed a novel semi-automatic algorithm BLL to assess the volume of pleural effusion. The algorithm is designed in two major steps: In the first step, the image within the inner chest wall was obtained by the B-spline method. In the second step, the region of pleural effusion was segmented by the LCLS method and obtained by the novel selecting and assembling criterion. Compared to existing methods, the volume of pleural effusion was calculated precisely by BLL. Moreover, the variation of the effusion volume was estimated within a period of time. So that physician(s) can assess patient’s new condition and adjust their medical treatment.

The limitations of our study consist of quantifying pleural effusion by several factors. First, in this study, authors believed that the CT scans were with good signal to noise ratio, because the physicians selected and detected the scans at the very beginning. As a result, scans with serious image problems which could not be detected had been dropped. However, in the everyday medical practice, CT scans may have different signal to noise levels and/or some “little” image problem. Authors did not consider these issues in this period of work. Thus, it should be completed in future. Second, this algorithm was semi-automatic and it still needed a manual input. It may lead some manual deviations. Therefore, it should be solved in following studies. In summary, the BLL algorithm developed and tested in this study can be applied in different areas to evaluate the volume of pneumothorax after several changes.

5 Conclusion

This study demonstrated the feasibility and advantages of applying a new semi-automatic algorithm namely, BLL, to segment the volume of pleural effusion on CT scans. This algorithm can also quantify the volume of pleural effusion precisely. It is widely applicable in patients with free pleural effusion.

Footnotes

Acknowledgments

This research is supported a Grant No.201402013 from National Health and Family Planning Commission of the People’s Republic of China.

References

Miserocchi

, Physiology and pathophysiology of pleural fluid turnover, Eur Respir J 10 (1997), 219–225.

Lai-Fook

S.J.

, Pleural mechanics and fluid exchange, Physiol Rev 84 (2004), 385–410.

Noppen

, De Waelle

, Li

et al., Volume and cellular content of normal pleural fluid in humans examined by pleural lavage, Am J Respir Crit Care Med 162 (2000), 1023–1026.

Estenne

, Yernault

J.-C.

et al., Mechanism of relief of dyspnea after thoracocentesis in patients with large pleural effusions, Am J Med 74 (1983), 813–819.

Craig Blackmore

, Black

W.C.

, Dallas

R.V.

and Crow

H.C.

, Pleural fluid volume estimation: A chest radiograph prediction rule, Acad Radiol 3 (1996), 103–109.

Eibenberger

K.L.

, Dock

W.I.

, Ammann

M.E.

et al., Quantification of pleural effusions: Sonography versus radiography, Radiology 191 (1994), 681–684.

Kalokairinos

M.I.

, Vassiliadis

S.G.

, Vossou

C.G.

et al., Estimation of the volume of the pleural effusion by computed tomography, Eur Sco Radial ECR 2010 (2010), C-0902.

Mergo

P.J.

, Helmberger

, Didovic

et al., New Formula for Quantification of Pleural Effusions From Computed Tomography, J of Thorac Imaging 14 (1999), 122–125.

Guo

Z.J.

, Lin

, Liu

H.T.

et al., The preliminary exploration of 64-slice volume computed tomography in the accurate measurement of pleural effusion, Acta Radiol 54 (2013), 765–769.

10.

Hazlinger

, Ctvrtlik

, Langova

and Herman

, Quantification of pleural effusion on CT by simple measurement, Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub 158 (2014), 107–111.

11.

Zhen

, Tavares

João Manuel R.S.

, Natal Jorge

Renato

and Mascarenhas

, A review of algorithms for medical image segmentation and their applications to the female pelvic cavity, Computer Methods in Biomechanics and Biomedical Engineering 13(2) (2009), 235–236.

12.

Zhen

M.A.

, Tavares

João Manuel R.S.

and Natal Jorge

R.M.

, A review on the current segmentation algorithms for medical images, 1st International Conference on Imaging Theory and Applications (IMAGAPP) (2009), 135–140.

13.

Jodas

Danilo Samuel

, Pereira

Aledir Silveira

and Tavares

João Manuel R.S.

, A review of computational methods applied for identification and quantification of atherosclerotic plaques in images, Expert Systems with Applications 46 (2016), 1–14.

14.

Zhen

M.A.

, Natal Jorge

Renato

, Mascarenhas

and Tavares

João Manuel R.S.

, Novel Approach to Segment the Inner and OuterBoundaries of the Bladder Wall in T2-weighted Magnetic ResonanceImages, Annals of Biomedical Engineering 39(8) (2011), 2287–2297.

15.

Zhen

M.A.

, Natal Jorge

Renato M.

, Mascarenhas

, Tavares

João Manuel R.S.

, Segmentation of Female Pelvic Organs in AxialMagnetic Resonance Images using Coupled Geometric DeformableModels”, Computers in Biology and Medicine 43(4) (2013), 248–258.

16.

Anderson

S.A.

, Danelson

K.A.

et al., A semiautomatic method of quantifying pleural effusions using computed tomograph, Biomed Sci Instru 49 (2013), 13–19.

17.

Bliton

, Yao

, Bi

, et al., Improved 3D automatic segmentation and measurement of pleural effusions, Proceedings (2011), 1954–1957.

18.

Kido and Tsunomori, Automated extractions of pleural effusion in three dimensional thoracic CT images, Proc. of SPIE, Vol. 7260, 72600D-1.

19.

Donohue et al., Constrained region-based segmentation of pleural effusion in thin-slice CT, Proc IEEE 13th Int Mach Vis Image Process Conf (2009).

20.

Qian et al., New method for quantification of pleural effusions from CT imaging, Proc ISECS Int Colloq Comput Commun Control Manage (2008).

21.

Yao et al., Automatic segmentation and measurement of pleural effusions on CT, IEEE Trans Biomed Eng 60 (2013), 1834–1840.

22.

, Huang

, Ding

, et al., A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI, IEEE Trans Image Process A 20(7) (2011), 2007–2016.

23.

Gordon

W.J.

and Riesenfeld

R.F.

, B-spline curves and surfaces, Eds. New York, Academic Press (1974), 95–126.

24.

Riesenfeld

R.F.

, Applications of B-spline approximation to geometric problems of computer aided design, Ph.D. dissertation, Syracuse Univ, 1973.

25.

Z.J.

, Zhang

Y.L.

, Ma

Z.C.

, et al., A Low-Pass Differentiation Filter Based on the 2nd-Order B-Spline Wavelet for Calculating Augmentation Index, Med Eng Phys 36 (2014), 786–792.

26.

, Heng

and Wong

, Cardiac motion recovery using an incompressible B-solid model, Med Eng Phys 35 (2013), 958–968.

27.

Sun

, Chang

, Dormer

K.J.

, Dyer

R.K.

and Gan

R.Z.

, An advanced computer-aided geometric modeling and fabrication method for human middle ear, Med Eng Phys 24 (2002), 595–606.

28.

Chan

and Vese

, Active contours without edges, IEEE Trans Image Process 10(2) (2001), 266–277.

29.

, Kao

, et al., Minimization of region-scalable fitting energy for image segmentation, IEEE Trans Image Process 17(10) (2008), 1940–1949.

30.

Aubert

and Kornprobst

, Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations, New York: Springer-Verlag, (2002).

31.

, Xu

, et al., Distance regularized level set evolution and its application to image segmentation, IEEE Trans Image Process 19(12) (2010), 3243–3254.

A quantitative evaluation of pleural effusion on computed tomography scans using B-spline and local clustering level set

Abstract

Keywords

1 Introduction

2 Materials and methods

2.1 Patient population and device

2.2 Extraction of the image within the inner chest wall

2.2.1 Extraction of interior points of thoracic cage

2.3.2.1 LCLS formulation

Table 1 Conversion of pleural effusion’s thickness to volume* Thickness d (cm) Pleural effusion volume (ml) 1 67 1.5 142 2 216 2.5 288 3 360 3.5 431 4 501 4.5 572 5 643 5.5 714 6 786 6.5 859 7 933 7.5 1009 8 1086 8.5 1166 9 1248 9.5 1333 10 1421

Footnotes

Acknowledgments

References

Table 1
Conversion of pleural effusion’s thickness to volume*

Thickness d (cm) Pleural effusion volume (ml)

1 67

1.5 142

2 216

2.5 288

3 360

3.5 431

4 501

4.5 572

5 643

5.5 714

6 786

6.5 859

7 933

7.5 1009

8 1086

8.5 1166

9 1248

9.5 1333

10 1421