- Email: [email protected]
A two-step method for preprocessing volume data
Computer Vision and Image Understanding 95 (2004) 150–164 www.elsevier.com/locate/cviu
A two-step method for preprocessing volume data Cheng Bing,a,b,* Wang Ying,c Zheng Nanning,a Bian Zhengzhong,b and Zhang Yongpinga a
Artificial Intelligence and Robotics Institute, Xi’ an Jiaotong University, Xi’an, 710049, China b Department of Biomedical Engeering, Xi’ an Jiaotong University, Xi’an, 710049, China c Medical experiment center of No.1 hospital, Xi’ an Jiaotong University, Xi’an, 710049, China Received 10 August 2002 Available online 20 April 2004
Abstract A method to preprocess volume data for rendering and analyzing is proposed here. We segment the data by constructing a threshold super-surface, which is sensitive to noise. So first we reduce the noise using level sets. The data is decomposed into upper level sets and lower level sets, and the small connected components of each level set are removed. In the second step, the segmentation problem is modeled as a partial difference equation (PDE). The super-surface is obtained by solving this PDE. Finally, the performance of our method is demonstrated with an experiment. Ó 2004 Elsevier Inc. All rights reserved. IDT: Segmentation; Noise reduction; Level sets; Volume data; Visualization
1. Introduction Since modern imaging modalities deliver huge amounts of data, which cannot be assessed easily, visualization techniques are utilized to emphasize the structures of interest. Visualization of volume data in 3D space is a rapidly growing application of computer graphics. This is especially true in medical volume data, where medical imaging modalities such as computed tomography (CT), magnetic resonance *
Corresponding author. E-mail address: [email protected] (C. Bing).
1077-3142/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.cviu.2004.03.003
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
151
imaging (MRI), and positron emission tomography (PET), etc., have been developing very fast recently [1,2]. The volume data we obtained cannot be visualized directly, because there are a lot of noise during acquiring the data and we do not know which tissues the voxels belong to. So the volume data need to be preprocessed before rendering. The segmentation allows extraction of the interested object before 3D representation is constructed. Due to a variety of factors affecting the data acquisition, volume data segmentation is a challenging task. In particular, the data is often corrupted by noise and sampling artifacts [3], which can cause considerable difficulties when applying classical segmentation techniques such as edge detection and thresholding [4–6]. As a result, these techniques either fail completely or require some kind of postprocessing step to remove invalid object boundaries in the segmentation results. Some issues that are important in a general purpose segmentation system can be ignored in medical volume data segmentation, such as illumination and pose determination. Volume data captured by devices such as biomedical scanners have inherent inaccuracies. The degree of this inaccuracy depends on a number of factors including limitations in spatial, temporal, and parametric resolutions and other physical limitations of the device. The objects to be segmented from medical imagery are actual anatomical structures. The anatomical structures in the human body are often nonrigid and complex in shape. They exhibit considerable variation from person to person and we have not the explicit shape models that capture the deformations in anatomy. The data would also have some motion effects and noise because of movement of the subject during the scan and the effects of the imaging facility. In this paper, we present a method, which is composed of two steps, to preprocess the volume data. The noise in the data is reduced and the interested object, which is segmented in a sense into a binary image, is obtained by the method for the future 3D operations of visualization, manipulation, and analysis, etc. The evaluation of segmentation methods is also a difficult problem because of the lack of a unified framework for assessment of the performance. Methods published expressly for addressing segmentation evaluation are rare and are very restricted in scope [7]. In the experiment section, we give a careful evaluation of the method in a systematic way according to the framework proposed by Udupa and co-workers [8,9].
2. Noise reduction The raw data typically are a spatial sequence of image matrices of 2562 pixels acquired through CT, MRI, PET, etc. To save storage space the gray-level data of the original volume is compressed to a dynamic range of 256 gray values. To achieve cubic volume elements, a linear interpolation of the intensity values between the original slices is performed. This results in a 3D array (voxel model). The volume data can be realistically modeled as a function uðxÞ defined on a cube: X ¼ ð0; W Þ ð0; H Þ ð0; DÞ ! R . Where X is the support set of the function. R is the available set of the value of voxels corresponding to a physical property of the tissues or organs of the body, e.g., the X-ray attenuation coefficient in CT.
152
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
A real function f defined on Rn is upper semi-continuous, if it satisfies: for any x 2 Rn and any real k with k > f ðxÞ, there exists a neighbor field Bx of x, and for any y 2 Bx ; f ðyÞ < k holds. The function uðxÞ can be looked on as an upper semi-continuous function because it is digitalized. We cannot tell the difference between an arbitrary function and a semi-continuous function using our eyes. A function can be decomposed into level sets if it is upper semi-continuous. Upper level set with level k is defined by [10] vk u ¼ X k ¼ fx 2 R3 juðxÞ P kg; k 2 R :
ð1Þ
Lower level set is defined in exactly the same way by changing ‘‘ P ’’ into ‘‘<.’’ vk u ¼ Xk ¼ fx 2 R3 juðxÞ < kg; k 2 R :
ð2Þ
The level sets give a complete account of the function. Indeed, the image u can be recovered from its upper level sets vk u by the formula ð3Þ uðxÞ ¼ supfkjx 2 X k g: The reconstruction formula corresponding to lower level sets is as follows: uðxÞ ¼ inffkjx 2 Xk g:
ð4Þ
The level sets of a function have three properties: (i) Xk XTl , if k > l, (ii) Xk ¼ Xl , l
(iii) If ðXk Þk2R is a family of subsets of Rn satisfying (i) and (ii), then the function with values in R, uðxÞ ¼ inffl; x 2 Xl g satisfies vk u ¼ Xk . The value of the volume data represents the physical property, so we have no concern for the value itself. What we care about is the contrast between the values and the ordering of the values that can help us to distinguish the different tissues and organs. A striking property of the level sets is their global invariance to contrast changes. Two functions v and u have globally the same level sets if for every k there is l such that vl v ¼ vk u, and conversely. If u is applied to a contrast change understood as a continuous increasing function g, then it is easily checked that v ¼ gðuÞ and u have globally the same level sets. We study operators that remove peaks in the volume data. Such peaks are often created by impulse noise. Let X R3 be a closed set. We say that X is connected if it is not the union of two disjoint nonempty open sets. We call the connected component of a point x in X , ccðx; X Þ, the maximal connected subset of X containing x. It is a closed S set. A level set can be seen as the union of all its connected components:. vk u ¼ i ccðxi ; vk uÞ. Every connected component is given a measure that is the volume the connected component occupies. The measure is defined by measðccðxi ; vk uÞÞ ¼ volumeðccðxi ; vk uÞÞ:
ð5Þ
Fix a scale parameter a > 0. We define a set denoising operator on F , the set of all compact sets of R3 in the following way: Let X 2 F be a compact set. Then X is the union of all its connected components, X ¼ [i ccðxi ; X Þ and this decomposition is
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
153
uniquely defined. Then all connected components of measure strictly less than a are removed from X . We therefore, define a set operator: [ Ta ðX Þ ¼ ccðxi ; X Þ: ð6Þ measðccðxi ;X ÞÞ P a
First, the function uðxÞ is decomposed into level sets. Then the level sets are processed using the set operator Ta to get rid of the small connected components. Finally, the volume data is reconstructed using the processed level sets and the noise is removed. The process of denoising is shown in Fig. 1. The number of level sets that need to be processed is 2n, if the gray-level of the volume data is n. A large amount of memory and a lot of computing time are needed if we decompose the data into upper and lower level sets and get rid of the small connected components in every level set. Therefore, a fast algorithm is designed using a heap data structure to reduce memory consumption and to improve computing efficiency. We scan the whole volume data from left to right, top to bottom, and front to back to find the local extremum. The voxel of local extremum is put in a buffer, and the value of the voxel is recorded as Current_Value. The voxels in its 6-neighbor are put in the heap data structure, then the voxel having the largest value at the top of the heap are put in the buffer and Current_Value is updated by its value. The process is continued in the new voxelsÕ 6-neighbors until the number of voxels in the buffer exceeds a threshold that corresponds to the scale parameter a. Finally, the value of voxels in the buffer is updated using Current_Value. The following simple pseudo-code clarifies how the removal of the small connected components in the upper level sets can be applied: [1] Begin; [2] set threshold a and clear Buffer and Heap; // Begin to scan the whole volume data [3] compare uði; j; kÞ with the voxels in its 6-neighbor to find Local_MaxVoxel; // uði; j; kÞ is Local_MaxVoxel if its value is not // less than the value of voxels in its 6-neighbor. [4] put Local_MaxVoxelði; j; kÞ into Buffer; [5] Current_Value ¼ Local_MaxVoxelði; j; kÞ; [6] put the voxels in 6-neighbor of Local_MaxVoxel into Heap; [7] get MaxVoxel at the top of Heap; // The value of MaxVoxel is not less than the
Fig. 1. The process of denoising.
154
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
// value of the other voxels in Heap. compare Current_Value with MaxVoxel; if Current_Value < MaxVoxel goto [14]; Current_Value ¼ MaxVoxelði; j; kÞ; get Number of the voxels in the Buffer; if (Number < aÞ { put MaxVoxel into Buffer; put the voxels in its 6-neighbor into Heap; } [13] goto [7]; [14] update the value of voxels in the Buffer with Current_Value; // The value of the voxels in small connected // component is replaced by the value of their // neighbors. [15] clear Buffer and Heap; [16] if the whole volume data is not scanned goto [3]; [17] end; A heap data structure is used in the above process. The elements in the heap data structure are ordered from large to small automatically when they are added into the structure. So the voxel at the top of the heap is larger than the others in the heap. The complexity of the algorithm of adding a voxel and getting the MaxVoxel is Oðlog nÞ, where n is number of elements in the heap. The process of removing the small connected components in the lower level sets is similar to above process.
[8] [9] [10] [11] [12]
3. Segmentation The noise is reduced and the edges are preserved after removing the small connected components in the original data. In the second step, we construct a threshold super-surface that varies by spatial position. The super-surface is obtained using edge information, i.e., the points corresponding to zero-crossings serve as the interpolation points of the threshold super-surface. And the interpolation problem is solved by minimizing an energy function. The points corresponding to zero-crossings, where Laplacian of u is zero, are chosen as interpolating points. The condition of interpolation is to have the supersurface pass these points and to minimize the following energy function: ZZZ LðuÞ ¼ ð7Þ jDuj2 dx dy dz; X o2 u ox2
2
2
where Du ¼ þ ooyu2 þ ooz2u is the Laplacian of u. The edges of objects are usually broken if directly counting the zero-crossings of the data. The data is convoluted with a symmetric smooth kernel u and the energy function of (7) is replaced by:
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
ZZZ
LðuÞ ¼
2
jDðu uÞj dx dy dz ¼
X
ZZZ
155
2
jðDuÞ uj dx dy dz:
ð8Þ
X
The smooth kernel u is chosen as follows c exp x2 þy 21þz2 1 x2 þ y 2 þ z2 < 1; uðx; y; zÞ ¼ 0 otherwise;
ð9Þ
where c is a constant and is defined by c ¼ RRR
1 : uðx; y; zÞ dx dy dz X
ð10Þ
The kernel function u is infinitely smooth with compact support, and its Laplacian is as follows: Wðx; y; zÞ ¼ Duðx; y; zÞ:
ð11Þ
The W is symmetric from u being symmetric, i.e. Wðx; y; zÞ ¼ Wðx; y; zÞ:
ð12Þ
A constraint is added to make the super-surface smooth enough and the energy function is translated into the following form: ZZZ ZZZ 2 SðuÞ ¼ a jðDuÞ uj dx dy dz þ b f ðkrukÞ dx dy dz; ð13Þ X
X
where ru represents the gradient vector of function u, k:k represents the norm of a vector, f is an increasing positive function, a and b are positive parameters which can balance the two terms on the right of (13). Such a process of optimization preserves the edge information. We now seek to minimize SðuÞ by considering the first Gateaux variation of SðuÞ at u in the direction h, as follows: dSðu : hÞ ¼ lim k!0
Sðu þ khÞ SðuÞ : k
ð14Þ
If u is a solution of the above optimization problem, dSðu : hÞ ¼ 0 holds for every small function h. The Gateaux variation of SðuÞ is calculated as follows: dSðu : hÞ Sðu þ khÞ SðuÞ k Z Z Z ZZZ 2 ¼ lim a ðW ðu þ khÞÞ dx dy dz þ b f ðkrðu þ khÞkÞ dx dy dz ¼ lim k!0
k!0
X
X
ZZZ ZZZ 2 a ðW uÞ dx dy dz þ b f ðkrukÞ dx dy dz k X
¼a
ZZZ X
ðW uÞðW hÞ dx dy dz þ b
X
ZZZ
f 0 ðkrukÞ X
ru rh dx dy dz kruk
156
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
¼a
ZZZ
ZZZ
ðW uÞ Wðx s; y t; z pÞhðs; t; pÞ ds dt dp dx dy dz X Z ZXZ ru rh dx dy dz þb f 0 ðkrukÞ kruk X ZZZ ZZZ ¼a hðs; t; pÞ ðW uÞWðx s; y t; z pÞ dx dy dz dsdt dp X Z ZXZ ru rh dx dy dz þb f 0 ðkrukÞ kruk X ZZZ ZZZ ¼a hðs; t; pÞ ðW uÞWðs x; t y; p zÞ dx dy dz ds dt dp X Z ZX Z ru rh dx dy dz þb f 0 ðkrukÞ kruk ZZZ X ¼a ðW ðW uÞÞhðs; t; pÞ ds dt dp Z ZX Z ru rh dx dy dz þb f 0 ðkrukÞ kruk ZZZ X ¼a ðW ðW uÞÞhðx; y; zÞ dx dy dz X ZZZ ru rh dx dy dz: þb f 0 ðkrukÞ kruk X
By assuming symmetric boundary conditions, i.e. *
ruðx; y; zÞ n ¼ 0;
for ðx; y; zÞ 2 oX;
ð15Þ *
where oX is the boundary of image domain X and n is the outward normal to oX, and from Green theorem, we have ZZZ ðW ðW uÞÞh dx dy dz dSðu : hÞ ¼ 2a Z ZXZ ru 0 div f ðkrukÞ h dx dy dz ¼ 2hrSðuÞ; hi; ð16Þ 2b kruk X where h ; i is the inner product defined by hf ; gi ¼
ZZZ
f ðx; y; zÞgðx; y; zÞ dx dy dz;
ð17Þ
X
and ru rSðuÞ ¼ aW ðW uÞ b div f 0 ðkrukÞ ; kruk
ð18Þ
is the gradient of SðuÞ at u. The functional SðuÞ is minimized by moving in the negative direction of the gradient through the following parabolic partial differential equation (PDE):
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
ou ru 0 ¼ rSðuÞ ¼ aW ðW uÞ þ b div f ðkrukÞ : ot kruk
157
ð19Þ
The energy function (13) has unique solution if function f is a convex function. In 2 this paper, we let f ðkrukÞ ¼ kruk . A numerical solution of the partial differential equation can be obtained by the following iteration: ukþ1 ¼ uk a W ðW uk Þ þ b Duk ;
ð20Þ
Eq. (20) is the mathematical formula of continuous models. But uðxÞ is digitalized and the voxel model is a discrete one. So the discrete form of Eq. (20) needs to be derived for digital model. The digital form of uðxÞ is the same as in Section 2. The smooth kernel u is digitalized as follows: 1 c exp i2 þj2 þk i2 þ j2 þ k 2 < 1; 2 1 uði; j; kÞ ¼ ð21Þ 0 otherwise; where jij < Nx ; jjj < Ny ; jkj < Nz ; Nx ; Ny ; Nz are radii in three directions, c is a constant and is defined by 1 c¼ P P P
uði; j; kÞ
:
ð22Þ
jij
The discrete three dimension Laplacian operator is defined as DLði; j; kÞ ¼ Lði þ 2; j; kÞ 2Lði þ 1; j; kÞ þ Lði; j þ 2; kÞ 2Lði; j þ 1; kÞ þ Lði; j; k þ 2Þ 2Lði; j; k þ 1Þ þ 3Lði; j; kÞ:
ð23Þ
After the threshold super-surface is acquired, we compare the value of every voxel with that of the corresponding point in the threshold super-surface u and segment the data by the following procedure: 1 if u0 ði; j; kÞ P uði; j; kÞ; Oði; j; kÞ ¼ ð24Þ 0 otherwise:
4. Experiment results and evaluation In this section, we present several examples to illustrate the behavior of our method. The first example, illustrated in Fig. 2, comes from a volume data set of a normal head. A slice of the head is shown in Fig. 2A. We reduced the noise in the data set using the level set method first, and the result is shown in Fig. 2B. And then the white matter region of the data set is segmented by Eqs. (20)–(24), the intermediate result at the fourth step of the iteration is shown in Fig. 2C. Fig. 2D is the final segmentation result, which iterated 15 steps. The other three examples are shown in Fig. 3. They are (A) white matter region in another head, (B) pathological region in a head with pathology, and (C) bones, tendons, and muscles in a knee. Fig. 4 is an example of noised and more complex images. We compare our denoising method
158
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
Fig. 2. Segment the white matter region of a normal head. (A) A slice of the head. (B) The result after denoising. (C) The intermediate result at the fourth step of the iteration. (D) The final segmented result.
Fig. 3. More examples, the first row represents the raw data of different slices, the second row indicates the corresponding segmented results: (A) another head, (B) a head with pathology, and (C) a knee.
with median filtering and Gaussian filtering. The first row indicates two slices of T2weighted images (TE ¼ 100 ms; TR ¼ 2266 ms), which show the injury in the region of the myotendinous junction of the right biceps femoris muscle. Note the feathery appearance in the biceps due to infiltration with blood and edema. The second row represents the denoised results using our method. The third row and the fourth row are the denoised results produced by median and Gaussian filter, respectively. The final segmented results are illustrated in the bottom row, which show the effectiveness of our algorithm. We will also present a quantitative comparative evaluation of this algorithm with the kMRFOE algorithm [11] using the first example. The MRI data is processed by some operations such as filtering to reduce image noise or enhance the edges, interpolating, and resampling to produce an isotropic volume
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
159
Fig. 4. An example of noised and more complex images. The first row indicates two slices of MRI images. The second row represents the denoised results using our method. The third row and the fourth row are the denoised results produced by median and Gaussian filter, respectively. The fifth row is the segmented results.
160
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
Fig. 4. (continued)
from a stack of cross-sectional images, and selection of a volume of interest from a large volumetric image, to obtain a new data set (with extended name *.vox). The data set is also used in the platform, Win-Volvis [12], which was developed in our lab. The manipulation and visualization are also done on the platform. The segmentation algorithm was implemented in C++ language and executed on a Pentium IV PC. The segmentation is done in three dimensions, although the results are shown for one section. To assess the effectiveness and utility of our segmentation method, we adopt an evaluation method proposed by Udupa. Some notations defined in UdupaÕs method are presented as follows for our description [9]. Application domain: evaluation must be performed for each application domain separately because an algorithm, that signals high performance for a given application domain, may tell nothing at all for a different application domain. Our application domain is defined as the segmentation of white matter region of brain using the data set mentioned above.
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
161
Scene: scene is a volume data of three dimensions, denoted by S ¼ ðC; f Þ, where C is a rectangular array of voxels, and f ðcÞ denotes the scene intensity of any voxel c in C. Fuzzy object: the output of any segmentation algorithm represents the region occupied by an interested object O in a scene S, which is a set of voxels in the scene O 2 S. The fuzzy object defined by O in S a scene S0 ¼ ðC; f0 Þ, where, for any c 2 C; f0 ðcÞ ¼ gðf ðcÞÞ. Here g is a function that assigns a degree of objectness to every voxel c in O depending on the scene intensity f ðcÞ. Some operations on fuzzy objects are defined as follows. Let S1 ¼ ðC; f1 Þ; S2 ¼ ðC; f2 Þ, and S3 ¼ ðC; f3 Þ be any fuzzy objects defined by the same object P O in a scene S. Then, the cardinality jS1 j of the fuzzy object S1 is defined as jS1 j ¼ c2C f1 ðcÞ. Fuzzy set union S1 ¼ S2 [ S3 is defined by, for any c 2 C; f1 ðcÞ ¼ maxðf2 ðcÞ; f3 ðcÞÞ. Fuzzy set intersection S1 ¼ S2 \ S3 is defined by, for any c 2 C; f1 ðcÞ ¼ minðf2 ðcÞ; f3 ðcÞÞ. The efficacy of any segmentation method in an application domain is to be measured in terms of three factors: precision, accuracy, and efficiency. Precision refers to the reproducibility of the segmentation results after taking into account all subjective actions that enter into the segmentation process. Accuracy relates to how well the segmentation results compare with the true delineation of the objects. Efficiency pertains to the practical viability of the method, which is determined by the amount of time required for performing computations and providing any user help needed in segmentation. Precision: in our application, precision is affected mainly by the way of selecting the interested region, the balance parameter a; b in Eq. (20) and the steps of iteration. To assess the precision of our method, four operators repeated the segmentation experiment five times. There are two situations should be considered here. The first situation, intra operator, is that the same operator segments the same object in the same scene twice by using our method. The second situation, inter operator, is that two operators segment the same object in the same scene once by using our method. The segmentation results are shown in Table 1, which lists the volumes of segmented white matter region. The mean intra operator coefficients of variation is 0.414%. The inter operator coefficients of variation is 0.571%. Let S1 and S2 be segmentations of the same object in two repeated trials. A measure of precision, differ in percentage overlap, is given by DðO1 ; O2 Þ ¼ jS1 \ S2 j=jS1 [ S2 j. Where O1 and O2 are segmented objects produced by a method in a trial. The measure represents the total amount of the tissue that is common to both O1 and O2 as a fraction of the total amount of tissue in the union of O1 and O2 . The precision of any two segmentation methods for each situation can
Table 1 Segmented volumes of white matter region
Operator Operator Operator Operator
1 2 3 4
Trial 1
Trial 2
Trial 3
Trial 4
Trial 5
8668 8451 8496 8539
8652 8527 8537 8602
8595 8519 8516 8581
8604 8426 8473 8593
8691 8472 8526 8487
162
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
be compared by comparing their DðO1 ; O2 Þ values. In the first example, the intra-operator and inter-operator percentage overlaps of segmentations of our method are 99.1 and 98.4%, respectively. Their counterparts in kMRFOE algorithm are 98.2 and 97.5%, respectively. Accuracy: of the three indexes used to describe the efficacy of a segmentation method, accuracy is the most difficult to assess. This is due mainly to the difficulty in establishing the true delineation of the object of interest in the scene. Consequently, an appropriate surrogate of truth is used in place of true delineation. The segmentation resulting from manual tracing of the object boundary by a knowledgeable operator is commonly used as this surrogate. We have taken this approach in our study. Let S0 be the scene representing the fuzzy object defined by an object O in S obtained by using method M, and let St be the corresponding scene of true delineation, under the same application domain. The following measures are defined to characterize the accuracy of method M for delineation: FNVF ¼
jSt S0 j ; jSt j
ð25Þ
FPVF ¼
jS0 St j ; jSt j
ð26Þ
TPVF ¼
jSt \ S0 j ; jSt j
ð27Þ
where fuzzy set difference S 0 ¼ S1 S2 is defined by, for any c 2 C, f1 ðcÞ f2 ðcÞ if f2 ðcÞ P 0; f 0 ðcÞ ¼ 0 otherwise:
ð28Þ
These measures are all expressed as a fraction of the volume of delineation. FNVF, which means false negative volume fraction, indicates the fraction of tissue defined in St that was missed by method M in delineation. FPVF, which means false positive volume fraction, denotes the amount of tissue falsely identified by method M as a fraction of the total amount of tissue in St . TPVF, which means true positive volume fraction, describes the fraction of the total amount of tissue in St with which the fuzzy object S0 overlaps. Note that the three measures are independent, that is none of them can be derived from the other two. The comparison of accuracy between our algorithm and the kMRFOE algorithm is shown in Table 2. Efficiency: efficiency refers to the practical viability of a segmentation method. Two factors need to be considered to fully characterize efficiency: computational time and the human operator time required to complete segmentation of each study in a routine setting in the application domain. Few of interaction is needed in our method, so the human operator time is ignored. To evaluate efficiency, we measured the computational time in several trials. The mean computational time study for the software package running with Windows XP and a Pentium IV PC (with 1.2 GHz CPU and 256 MB memory) is about 5.3 s. The mean computational time of the kMRFOE algorithm in the same condition is 2.1 s.
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
163
Table 2 The comparison of accuracy between our algorithm and kMRFOE algorithm FNVF (%)
Operator Operator Operator Operator
1 2 3 4
FPVF (%)
TPVF (%)
Our algorithm
kMRFOE algorithm
Our algorithm
kMRFOE algorithm
Our algorithm
kMRFOE algorithm
1.01 0.77 0.38 0.42
0.98 0.82 0.47 0.69
0.42 0.58 0.21 0.18
0.47 0.65 0.36 0.31
96.71 98.15 99.32 99.39
96.70 97.28 97.92 98.15
From above detailed assessment and comparison, we can get that our method is a little better than the kMRFOE algorithm in precision and accuracy, but our method need much computational time than the kMRFOE algorithm.
5. Concluding remarks In this paper, we have introduced a volume data preprocessing method for volume rendering, visualization, and analysis. Our method has two steps: denoising and segmentation. A partial difference equation (PDE) is solved by a recursive procedure to find a threshold super-surface using edge information. The edge information is not reliable due to the large amounts of noise in typical situations. So we use a level set based method to reduce the noise first, which can preserve the structure of the tissues and organs after denoising. Finally, a quantitative assessment of the method is given. The experimental results show that the method can segment the medical volume data efficaciously. Our next goal is to reduce the computational time and develop a preprocessing workshop, wherein, for any given application the denoising and segmentation protocol can be selected quickly and adaptively.
Acknowledgments The authors thank Dr. Quan Wei and Dr. Yuan Zejian for helpful discussions that greatly improved the paper from its original form. The authors thank Nick Pease for checking the paper for mistakes. This work was supported in part by National Nature Science Foundation of China (Grant No. 60003011) and The National High Technology Research and Development Program of China (Grant No. 2001AA114202).
References [1] R. Acharya, R. Wasserman, J. Sevens, C. Hinojosa, Biomedical imaging modalities: a tutorial, Comput. Med. Imag. Graph. 19 (1) (1995) 3–25. [2] M. Blank, W.A. Kalender, Medical volume exploration: gaining insights virtually, Eur. J. Radiol. 33 (2000) 161–169.
164
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
[3] H. Tang, T. Zhuang, Ed.X. Wu, Realizations of Fast 2-D/3-D image filtering and enhancement, IEEE Trans. Med. Imag. 20 (2) (2001) 132–140. [4] N. Papamarkos, C. Strouthopoulos, I. Andreadis, Multithresholding of color and gray-level images through a neural network technique, Image Vis. Comput. 18 (2000) 213–222. [5] J.R. Parker, Gray level thresholding in badly illuminated images, IEEE Trans. Pattern Anal. Mach. Intell. 13 (8) (1991) 813–819. [6] W. Tao, H. Burkhardt, An Effective Image Thresholding Method Using a Fuzzy Compactness Measure, in: The 12th International Conference on Pattern Recognition, Israel, October 1994. [7] F. Mao, J. Gill, A. Fenster, Technique for evaluation of semi-automatic segmentation methods, SPIE Proc. 3661 (1999) 1027–1036. [8] P.K. Saha, J.K. Udupa, Scale-based fuzzy connected image segmentation: theory, algorithm and validation, Comput. Vis. Image Understand. 77 (2) (2000) 145–174. [9] J.K. Udupa, V.R. LeBlance, et al., A methodology for evaluating image segmentation algorithms, SPIE Proc. 4684 (2002) 266–277. [10] F. Guichard, J.M. Morel, Partial differential equations and image iterative filtering, An Tutorial Int. Conf. Image Processing, Washington, DC, 1995. [11] P.K. Saha, J.K. Udupa, Relative fuzzy connectedness among multiple objects: theory, algorithms, and applications in image segmentation, Comput. Vis. Image Understand. 82 (1) (2001) 42–56. [12] L. Xinxiao, J. Li, Z. Nanning, Win-VolVis: an expandable visualization system in Win95/NT, J. Comput. Appl. Software of China (1) (2002).
A two-step method for preprocessing volume data Cheng Bing,a,b,* Wang Ying,c Zheng Nanning,a Bian Zhengzhong,b and Zhang Yongpinga a
Artificial Intelligence and Robotics Institute, Xi’ an Jiaotong University, Xi’an, 710049, China b Department of Biomedical Engeering, Xi’ an Jiaotong University, Xi’an, 710049, China c Medical experiment center of No.1 hospital, Xi’ an Jiaotong University, Xi’an, 710049, China Received 10 August 2002 Available online 20 April 2004
Abstract A method to preprocess volume data for rendering and analyzing is proposed here. We segment the data by constructing a threshold super-surface, which is sensitive to noise. So first we reduce the noise using level sets. The data is decomposed into upper level sets and lower level sets, and the small connected components of each level set are removed. In the second step, the segmentation problem is modeled as a partial difference equation (PDE). The super-surface is obtained by solving this PDE. Finally, the performance of our method is demonstrated with an experiment. Ó 2004 Elsevier Inc. All rights reserved. IDT: Segmentation; Noise reduction; Level sets; Volume data; Visualization
1. Introduction Since modern imaging modalities deliver huge amounts of data, which cannot be assessed easily, visualization techniques are utilized to emphasize the structures of interest. Visualization of volume data in 3D space is a rapidly growing application of computer graphics. This is especially true in medical volume data, where medical imaging modalities such as computed tomography (CT), magnetic resonance *
Corresponding author. E-mail address: [email protected] (C. Bing).
1077-3142/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.cviu.2004.03.003
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
151
imaging (MRI), and positron emission tomography (PET), etc., have been developing very fast recently [1,2]. The volume data we obtained cannot be visualized directly, because there are a lot of noise during acquiring the data and we do not know which tissues the voxels belong to. So the volume data need to be preprocessed before rendering. The segmentation allows extraction of the interested object before 3D representation is constructed. Due to a variety of factors affecting the data acquisition, volume data segmentation is a challenging task. In particular, the data is often corrupted by noise and sampling artifacts [3], which can cause considerable difficulties when applying classical segmentation techniques such as edge detection and thresholding [4–6]. As a result, these techniques either fail completely or require some kind of postprocessing step to remove invalid object boundaries in the segmentation results. Some issues that are important in a general purpose segmentation system can be ignored in medical volume data segmentation, such as illumination and pose determination. Volume data captured by devices such as biomedical scanners have inherent inaccuracies. The degree of this inaccuracy depends on a number of factors including limitations in spatial, temporal, and parametric resolutions and other physical limitations of the device. The objects to be segmented from medical imagery are actual anatomical structures. The anatomical structures in the human body are often nonrigid and complex in shape. They exhibit considerable variation from person to person and we have not the explicit shape models that capture the deformations in anatomy. The data would also have some motion effects and noise because of movement of the subject during the scan and the effects of the imaging facility. In this paper, we present a method, which is composed of two steps, to preprocess the volume data. The noise in the data is reduced and the interested object, which is segmented in a sense into a binary image, is obtained by the method for the future 3D operations of visualization, manipulation, and analysis, etc. The evaluation of segmentation methods is also a difficult problem because of the lack of a unified framework for assessment of the performance. Methods published expressly for addressing segmentation evaluation are rare and are very restricted in scope [7]. In the experiment section, we give a careful evaluation of the method in a systematic way according to the framework proposed by Udupa and co-workers [8,9].
2. Noise reduction The raw data typically are a spatial sequence of image matrices of 2562 pixels acquired through CT, MRI, PET, etc. To save storage space the gray-level data of the original volume is compressed to a dynamic range of 256 gray values. To achieve cubic volume elements, a linear interpolation of the intensity values between the original slices is performed. This results in a 3D array (voxel model). The volume data can be realistically modeled as a function uðxÞ defined on a cube: X ¼ ð0; W Þ ð0; H Þ ð0; DÞ ! R . Where X is the support set of the function. R is the available set of the value of voxels corresponding to a physical property of the tissues or organs of the body, e.g., the X-ray attenuation coefficient in CT.
152
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
A real function f defined on Rn is upper semi-continuous, if it satisfies: for any x 2 Rn and any real k with k > f ðxÞ, there exists a neighbor field Bx of x, and for any y 2 Bx ; f ðyÞ < k holds. The function uðxÞ can be looked on as an upper semi-continuous function because it is digitalized. We cannot tell the difference between an arbitrary function and a semi-continuous function using our eyes. A function can be decomposed into level sets if it is upper semi-continuous. Upper level set with level k is defined by [10] vk u ¼ X k ¼ fx 2 R3 juðxÞ P kg; k 2 R :
ð1Þ
Lower level set is defined in exactly the same way by changing ‘‘ P ’’ into ‘‘<.’’ vk u ¼ Xk ¼ fx 2 R3 juðxÞ < kg; k 2 R :
ð2Þ
The level sets give a complete account of the function. Indeed, the image u can be recovered from its upper level sets vk u by the formula ð3Þ uðxÞ ¼ supfkjx 2 X k g: The reconstruction formula corresponding to lower level sets is as follows: uðxÞ ¼ inffkjx 2 Xk g:
ð4Þ
The level sets of a function have three properties: (i) Xk XTl , if k > l, (ii) Xk ¼ Xl , l
(iii) If ðXk Þk2R is a family of subsets of Rn satisfying (i) and (ii), then the function with values in R, uðxÞ ¼ inffl; x 2 Xl g satisfies vk u ¼ Xk . The value of the volume data represents the physical property, so we have no concern for the value itself. What we care about is the contrast between the values and the ordering of the values that can help us to distinguish the different tissues and organs. A striking property of the level sets is their global invariance to contrast changes. Two functions v and u have globally the same level sets if for every k there is l such that vl v ¼ vk u, and conversely. If u is applied to a contrast change understood as a continuous increasing function g, then it is easily checked that v ¼ gðuÞ and u have globally the same level sets. We study operators that remove peaks in the volume data. Such peaks are often created by impulse noise. Let X R3 be a closed set. We say that X is connected if it is not the union of two disjoint nonempty open sets. We call the connected component of a point x in X , ccðx; X Þ, the maximal connected subset of X containing x. It is a closed S set. A level set can be seen as the union of all its connected components:. vk u ¼ i ccðxi ; vk uÞ. Every connected component is given a measure that is the volume the connected component occupies. The measure is defined by measðccðxi ; vk uÞÞ ¼ volumeðccðxi ; vk uÞÞ:
ð5Þ
Fix a scale parameter a > 0. We define a set denoising operator on F , the set of all compact sets of R3 in the following way: Let X 2 F be a compact set. Then X is the union of all its connected components, X ¼ [i ccðxi ; X Þ and this decomposition is
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
153
uniquely defined. Then all connected components of measure strictly less than a are removed from X . We therefore, define a set operator: [ Ta ðX Þ ¼ ccðxi ; X Þ: ð6Þ measðccðxi ;X ÞÞ P a
First, the function uðxÞ is decomposed into level sets. Then the level sets are processed using the set operator Ta to get rid of the small connected components. Finally, the volume data is reconstructed using the processed level sets and the noise is removed. The process of denoising is shown in Fig. 1. The number of level sets that need to be processed is 2n, if the gray-level of the volume data is n. A large amount of memory and a lot of computing time are needed if we decompose the data into upper and lower level sets and get rid of the small connected components in every level set. Therefore, a fast algorithm is designed using a heap data structure to reduce memory consumption and to improve computing efficiency. We scan the whole volume data from left to right, top to bottom, and front to back to find the local extremum. The voxel of local extremum is put in a buffer, and the value of the voxel is recorded as Current_Value. The voxels in its 6-neighbor are put in the heap data structure, then the voxel having the largest value at the top of the heap are put in the buffer and Current_Value is updated by its value. The process is continued in the new voxelsÕ 6-neighbors until the number of voxels in the buffer exceeds a threshold that corresponds to the scale parameter a. Finally, the value of voxels in the buffer is updated using Current_Value. The following simple pseudo-code clarifies how the removal of the small connected components in the upper level sets can be applied: [1] Begin; [2] set threshold a and clear Buffer and Heap; // Begin to scan the whole volume data [3] compare uði; j; kÞ with the voxels in its 6-neighbor to find Local_MaxVoxel; // uði; j; kÞ is Local_MaxVoxel if its value is not // less than the value of voxels in its 6-neighbor. [4] put Local_MaxVoxelði; j; kÞ into Buffer; [5] Current_Value ¼ Local_MaxVoxelði; j; kÞ; [6] put the voxels in 6-neighbor of Local_MaxVoxel into Heap; [7] get MaxVoxel at the top of Heap; // The value of MaxVoxel is not less than the
Fig. 1. The process of denoising.
154
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
// value of the other voxels in Heap. compare Current_Value with MaxVoxel; if Current_Value < MaxVoxel goto [14]; Current_Value ¼ MaxVoxelði; j; kÞ; get Number of the voxels in the Buffer; if (Number < aÞ { put MaxVoxel into Buffer; put the voxels in its 6-neighbor into Heap; } [13] goto [7]; [14] update the value of voxels in the Buffer with Current_Value; // The value of the voxels in small connected // component is replaced by the value of their // neighbors. [15] clear Buffer and Heap; [16] if the whole volume data is not scanned goto [3]; [17] end; A heap data structure is used in the above process. The elements in the heap data structure are ordered from large to small automatically when they are added into the structure. So the voxel at the top of the heap is larger than the others in the heap. The complexity of the algorithm of adding a voxel and getting the MaxVoxel is Oðlog nÞ, where n is number of elements in the heap. The process of removing the small connected components in the lower level sets is similar to above process.
[8] [9] [10] [11] [12]
3. Segmentation The noise is reduced and the edges are preserved after removing the small connected components in the original data. In the second step, we construct a threshold super-surface that varies by spatial position. The super-surface is obtained using edge information, i.e., the points corresponding to zero-crossings serve as the interpolation points of the threshold super-surface. And the interpolation problem is solved by minimizing an energy function. The points corresponding to zero-crossings, where Laplacian of u is zero, are chosen as interpolating points. The condition of interpolation is to have the supersurface pass these points and to minimize the following energy function: ZZZ LðuÞ ¼ ð7Þ jDuj2 dx dy dz; X o2 u ox2
2
2
where Du ¼ þ ooyu2 þ ooz2u is the Laplacian of u. The edges of objects are usually broken if directly counting the zero-crossings of the data. The data is convoluted with a symmetric smooth kernel u and the energy function of (7) is replaced by:
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
ZZZ
LðuÞ ¼
2
jDðu uÞj dx dy dz ¼
X
ZZZ
155
2
jðDuÞ uj dx dy dz:
ð8Þ
X
The smooth kernel u is chosen as follows c exp x2 þy 21þz2 1 x2 þ y 2 þ z2 < 1; uðx; y; zÞ ¼ 0 otherwise;
ð9Þ
where c is a constant and is defined by c ¼ RRR
1 : uðx; y; zÞ dx dy dz X
ð10Þ
The kernel function u is infinitely smooth with compact support, and its Laplacian is as follows: Wðx; y; zÞ ¼ Duðx; y; zÞ:
ð11Þ
The W is symmetric from u being symmetric, i.e. Wðx; y; zÞ ¼ Wðx; y; zÞ:
ð12Þ
A constraint is added to make the super-surface smooth enough and the energy function is translated into the following form: ZZZ ZZZ 2 SðuÞ ¼ a jðDuÞ uj dx dy dz þ b f ðkrukÞ dx dy dz; ð13Þ X
X
where ru represents the gradient vector of function u, k:k represents the norm of a vector, f is an increasing positive function, a and b are positive parameters which can balance the two terms on the right of (13). Such a process of optimization preserves the edge information. We now seek to minimize SðuÞ by considering the first Gateaux variation of SðuÞ at u in the direction h, as follows: dSðu : hÞ ¼ lim k!0
Sðu þ khÞ SðuÞ : k
ð14Þ
If u is a solution of the above optimization problem, dSðu : hÞ ¼ 0 holds for every small function h. The Gateaux variation of SðuÞ is calculated as follows: dSðu : hÞ Sðu þ khÞ SðuÞ k Z Z Z ZZZ 2 ¼ lim a ðW ðu þ khÞÞ dx dy dz þ b f ðkrðu þ khÞkÞ dx dy dz ¼ lim k!0
k!0
X
X
ZZZ ZZZ 2 a ðW uÞ dx dy dz þ b f ðkrukÞ dx dy dz k X
¼a
ZZZ X
ðW uÞðW hÞ dx dy dz þ b
X
ZZZ
f 0 ðkrukÞ X
ru rh dx dy dz kruk
156
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
¼a
ZZZ
ZZZ
ðW uÞ Wðx s; y t; z pÞhðs; t; pÞ ds dt dp dx dy dz X Z ZXZ ru rh dx dy dz þb f 0 ðkrukÞ kruk X ZZZ ZZZ ¼a hðs; t; pÞ ðW uÞWðx s; y t; z pÞ dx dy dz dsdt dp X Z ZXZ ru rh dx dy dz þb f 0 ðkrukÞ kruk X ZZZ ZZZ ¼a hðs; t; pÞ ðW uÞWðs x; t y; p zÞ dx dy dz ds dt dp X Z ZX Z ru rh dx dy dz þb f 0 ðkrukÞ kruk ZZZ X ¼a ðW ðW uÞÞhðs; t; pÞ ds dt dp Z ZX Z ru rh dx dy dz þb f 0 ðkrukÞ kruk ZZZ X ¼a ðW ðW uÞÞhðx; y; zÞ dx dy dz X ZZZ ru rh dx dy dz: þb f 0 ðkrukÞ kruk X
By assuming symmetric boundary conditions, i.e. *
ruðx; y; zÞ n ¼ 0;
for ðx; y; zÞ 2 oX;
ð15Þ *
where oX is the boundary of image domain X and n is the outward normal to oX, and from Green theorem, we have ZZZ ðW ðW uÞÞh dx dy dz dSðu : hÞ ¼ 2a Z ZXZ ru 0 div f ðkrukÞ h dx dy dz ¼ 2hrSðuÞ; hi; ð16Þ 2b kruk X where h ; i is the inner product defined by hf ; gi ¼
ZZZ
f ðx; y; zÞgðx; y; zÞ dx dy dz;
ð17Þ
X
and ru rSðuÞ ¼ aW ðW uÞ b div f 0 ðkrukÞ ; kruk
ð18Þ
is the gradient of SðuÞ at u. The functional SðuÞ is minimized by moving in the negative direction of the gradient through the following parabolic partial differential equation (PDE):
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
ou ru 0 ¼ rSðuÞ ¼ aW ðW uÞ þ b div f ðkrukÞ : ot kruk
157
ð19Þ
The energy function (13) has unique solution if function f is a convex function. In 2 this paper, we let f ðkrukÞ ¼ kruk . A numerical solution of the partial differential equation can be obtained by the following iteration: ukþ1 ¼ uk a W ðW uk Þ þ b Duk ;
ð20Þ
Eq. (20) is the mathematical formula of continuous models. But uðxÞ is digitalized and the voxel model is a discrete one. So the discrete form of Eq. (20) needs to be derived for digital model. The digital form of uðxÞ is the same as in Section 2. The smooth kernel u is digitalized as follows: 1 c exp i2 þj2 þk i2 þ j2 þ k 2 < 1; 2 1 uði; j; kÞ ¼ ð21Þ 0 otherwise; where jij < Nx ; jjj < Ny ; jkj < Nz ; Nx ; Ny ; Nz are radii in three directions, c is a constant and is defined by 1 c¼ P P P
uði; j; kÞ
:
ð22Þ
jij
The discrete three dimension Laplacian operator is defined as DLði; j; kÞ ¼ Lði þ 2; j; kÞ 2Lði þ 1; j; kÞ þ Lði; j þ 2; kÞ 2Lði; j þ 1; kÞ þ Lði; j; k þ 2Þ 2Lði; j; k þ 1Þ þ 3Lði; j; kÞ:
ð23Þ
After the threshold super-surface is acquired, we compare the value of every voxel with that of the corresponding point in the threshold super-surface u and segment the data by the following procedure: 1 if u0 ði; j; kÞ P uði; j; kÞ; Oði; j; kÞ ¼ ð24Þ 0 otherwise:
4. Experiment results and evaluation In this section, we present several examples to illustrate the behavior of our method. The first example, illustrated in Fig. 2, comes from a volume data set of a normal head. A slice of the head is shown in Fig. 2A. We reduced the noise in the data set using the level set method first, and the result is shown in Fig. 2B. And then the white matter region of the data set is segmented by Eqs. (20)–(24), the intermediate result at the fourth step of the iteration is shown in Fig. 2C. Fig. 2D is the final segmentation result, which iterated 15 steps. The other three examples are shown in Fig. 3. They are (A) white matter region in another head, (B) pathological region in a head with pathology, and (C) bones, tendons, and muscles in a knee. Fig. 4 is an example of noised and more complex images. We compare our denoising method
158
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
Fig. 2. Segment the white matter region of a normal head. (A) A slice of the head. (B) The result after denoising. (C) The intermediate result at the fourth step of the iteration. (D) The final segmented result.
Fig. 3. More examples, the first row represents the raw data of different slices, the second row indicates the corresponding segmented results: (A) another head, (B) a head with pathology, and (C) a knee.
with median filtering and Gaussian filtering. The first row indicates two slices of T2weighted images (TE ¼ 100 ms; TR ¼ 2266 ms), which show the injury in the region of the myotendinous junction of the right biceps femoris muscle. Note the feathery appearance in the biceps due to infiltration with blood and edema. The second row represents the denoised results using our method. The third row and the fourth row are the denoised results produced by median and Gaussian filter, respectively. The final segmented results are illustrated in the bottom row, which show the effectiveness of our algorithm. We will also present a quantitative comparative evaluation of this algorithm with the kMRFOE algorithm [11] using the first example. The MRI data is processed by some operations such as filtering to reduce image noise or enhance the edges, interpolating, and resampling to produce an isotropic volume
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
159
Fig. 4. An example of noised and more complex images. The first row indicates two slices of MRI images. The second row represents the denoised results using our method. The third row and the fourth row are the denoised results produced by median and Gaussian filter, respectively. The fifth row is the segmented results.
160
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
Fig. 4. (continued)
from a stack of cross-sectional images, and selection of a volume of interest from a large volumetric image, to obtain a new data set (with extended name *.vox). The data set is also used in the platform, Win-Volvis [12], which was developed in our lab. The manipulation and visualization are also done on the platform. The segmentation algorithm was implemented in C++ language and executed on a Pentium IV PC. The segmentation is done in three dimensions, although the results are shown for one section. To assess the effectiveness and utility of our segmentation method, we adopt an evaluation method proposed by Udupa. Some notations defined in UdupaÕs method are presented as follows for our description [9]. Application domain: evaluation must be performed for each application domain separately because an algorithm, that signals high performance for a given application domain, may tell nothing at all for a different application domain. Our application domain is defined as the segmentation of white matter region of brain using the data set mentioned above.
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
161
Scene: scene is a volume data of three dimensions, denoted by S ¼ ðC; f Þ, where C is a rectangular array of voxels, and f ðcÞ denotes the scene intensity of any voxel c in C. Fuzzy object: the output of any segmentation algorithm represents the region occupied by an interested object O in a scene S, which is a set of voxels in the scene O 2 S. The fuzzy object defined by O in S a scene S0 ¼ ðC; f0 Þ, where, for any c 2 C; f0 ðcÞ ¼ gðf ðcÞÞ. Here g is a function that assigns a degree of objectness to every voxel c in O depending on the scene intensity f ðcÞ. Some operations on fuzzy objects are defined as follows. Let S1 ¼ ðC; f1 Þ; S2 ¼ ðC; f2 Þ, and S3 ¼ ðC; f3 Þ be any fuzzy objects defined by the same object P O in a scene S. Then, the cardinality jS1 j of the fuzzy object S1 is defined as jS1 j ¼ c2C f1 ðcÞ. Fuzzy set union S1 ¼ S2 [ S3 is defined by, for any c 2 C; f1 ðcÞ ¼ maxðf2 ðcÞ; f3 ðcÞÞ. Fuzzy set intersection S1 ¼ S2 \ S3 is defined by, for any c 2 C; f1 ðcÞ ¼ minðf2 ðcÞ; f3 ðcÞÞ. The efficacy of any segmentation method in an application domain is to be measured in terms of three factors: precision, accuracy, and efficiency. Precision refers to the reproducibility of the segmentation results after taking into account all subjective actions that enter into the segmentation process. Accuracy relates to how well the segmentation results compare with the true delineation of the objects. Efficiency pertains to the practical viability of the method, which is determined by the amount of time required for performing computations and providing any user help needed in segmentation. Precision: in our application, precision is affected mainly by the way of selecting the interested region, the balance parameter a; b in Eq. (20) and the steps of iteration. To assess the precision of our method, four operators repeated the segmentation experiment five times. There are two situations should be considered here. The first situation, intra operator, is that the same operator segments the same object in the same scene twice by using our method. The second situation, inter operator, is that two operators segment the same object in the same scene once by using our method. The segmentation results are shown in Table 1, which lists the volumes of segmented white matter region. The mean intra operator coefficients of variation is 0.414%. The inter operator coefficients of variation is 0.571%. Let S1 and S2 be segmentations of the same object in two repeated trials. A measure of precision, differ in percentage overlap, is given by DðO1 ; O2 Þ ¼ jS1 \ S2 j=jS1 [ S2 j. Where O1 and O2 are segmented objects produced by a method in a trial. The measure represents the total amount of the tissue that is common to both O1 and O2 as a fraction of the total amount of tissue in the union of O1 and O2 . The precision of any two segmentation methods for each situation can
Table 1 Segmented volumes of white matter region
Operator Operator Operator Operator
1 2 3 4
Trial 1
Trial 2
Trial 3
Trial 4
Trial 5
8668 8451 8496 8539
8652 8527 8537 8602
8595 8519 8516 8581
8604 8426 8473 8593
8691 8472 8526 8487
162
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
be compared by comparing their DðO1 ; O2 Þ values. In the first example, the intra-operator and inter-operator percentage overlaps of segmentations of our method are 99.1 and 98.4%, respectively. Their counterparts in kMRFOE algorithm are 98.2 and 97.5%, respectively. Accuracy: of the three indexes used to describe the efficacy of a segmentation method, accuracy is the most difficult to assess. This is due mainly to the difficulty in establishing the true delineation of the object of interest in the scene. Consequently, an appropriate surrogate of truth is used in place of true delineation. The segmentation resulting from manual tracing of the object boundary by a knowledgeable operator is commonly used as this surrogate. We have taken this approach in our study. Let S0 be the scene representing the fuzzy object defined by an object O in S obtained by using method M, and let St be the corresponding scene of true delineation, under the same application domain. The following measures are defined to characterize the accuracy of method M for delineation: FNVF ¼
jSt S0 j ; jSt j
ð25Þ
FPVF ¼
jS0 St j ; jSt j
ð26Þ
TPVF ¼
jSt \ S0 j ; jSt j
ð27Þ
where fuzzy set difference S 0 ¼ S1 S2 is defined by, for any c 2 C, f1 ðcÞ f2 ðcÞ if f2 ðcÞ P 0; f 0 ðcÞ ¼ 0 otherwise:
ð28Þ
These measures are all expressed as a fraction of the volume of delineation. FNVF, which means false negative volume fraction, indicates the fraction of tissue defined in St that was missed by method M in delineation. FPVF, which means false positive volume fraction, denotes the amount of tissue falsely identified by method M as a fraction of the total amount of tissue in St . TPVF, which means true positive volume fraction, describes the fraction of the total amount of tissue in St with which the fuzzy object S0 overlaps. Note that the three measures are independent, that is none of them can be derived from the other two. The comparison of accuracy between our algorithm and the kMRFOE algorithm is shown in Table 2. Efficiency: efficiency refers to the practical viability of a segmentation method. Two factors need to be considered to fully characterize efficiency: computational time and the human operator time required to complete segmentation of each study in a routine setting in the application domain. Few of interaction is needed in our method, so the human operator time is ignored. To evaluate efficiency, we measured the computational time in several trials. The mean computational time study for the software package running with Windows XP and a Pentium IV PC (with 1.2 GHz CPU and 256 MB memory) is about 5.3 s. The mean computational time of the kMRFOE algorithm in the same condition is 2.1 s.
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
163
Table 2 The comparison of accuracy between our algorithm and kMRFOE algorithm FNVF (%)
Operator Operator Operator Operator
1 2 3 4
FPVF (%)
TPVF (%)
Our algorithm
kMRFOE algorithm
Our algorithm
kMRFOE algorithm
Our algorithm
kMRFOE algorithm
1.01 0.77 0.38 0.42
0.98 0.82 0.47 0.69
0.42 0.58 0.21 0.18
0.47 0.65 0.36 0.31
96.71 98.15 99.32 99.39
96.70 97.28 97.92 98.15
From above detailed assessment and comparison, we can get that our method is a little better than the kMRFOE algorithm in precision and accuracy, but our method need much computational time than the kMRFOE algorithm.
5. Concluding remarks In this paper, we have introduced a volume data preprocessing method for volume rendering, visualization, and analysis. Our method has two steps: denoising and segmentation. A partial difference equation (PDE) is solved by a recursive procedure to find a threshold super-surface using edge information. The edge information is not reliable due to the large amounts of noise in typical situations. So we use a level set based method to reduce the noise first, which can preserve the structure of the tissues and organs after denoising. Finally, a quantitative assessment of the method is given. The experimental results show that the method can segment the medical volume data efficaciously. Our next goal is to reduce the computational time and develop a preprocessing workshop, wherein, for any given application the denoising and segmentation protocol can be selected quickly and adaptively.
Acknowledgments The authors thank Dr. Quan Wei and Dr. Yuan Zejian for helpful discussions that greatly improved the paper from its original form. The authors thank Nick Pease for checking the paper for mistakes. This work was supported in part by National Nature Science Foundation of China (Grant No. 60003011) and The National High Technology Research and Development Program of China (Grant No. 2001AA114202).
References [1] R. Acharya, R. Wasserman, J. Sevens, C. Hinojosa, Biomedical imaging modalities: a tutorial, Comput. Med. Imag. Graph. 19 (1) (1995) 3–25. [2] M. Blank, W.A. Kalender, Medical volume exploration: gaining insights virtually, Eur. J. Radiol. 33 (2000) 161–169.
164
C. Bing et al. / Computer Vision and Image Understanding 95 (2004) 150–164
[3] H. Tang, T. Zhuang, Ed.X. Wu, Realizations of Fast 2-D/3-D image filtering and enhancement, IEEE Trans. Med. Imag. 20 (2) (2001) 132–140. [4] N. Papamarkos, C. Strouthopoulos, I. Andreadis, Multithresholding of color and gray-level images through a neural network technique, Image Vis. Comput. 18 (2000) 213–222. [5] J.R. Parker, Gray level thresholding in badly illuminated images, IEEE Trans. Pattern Anal. Mach. Intell. 13 (8) (1991) 813–819. [6] W. Tao, H. Burkhardt, An Effective Image Thresholding Method Using a Fuzzy Compactness Measure, in: The 12th International Conference on Pattern Recognition, Israel, October 1994. [7] F. Mao, J. Gill, A. Fenster, Technique for evaluation of semi-automatic segmentation methods, SPIE Proc. 3661 (1999) 1027–1036. [8] P.K. Saha, J.K. Udupa, Scale-based fuzzy connected image segmentation: theory, algorithm and validation, Comput. Vis. Image Understand. 77 (2) (2000) 145–174. [9] J.K. Udupa, V.R. LeBlance, et al., A methodology for evaluating image segmentation algorithms, SPIE Proc. 4684 (2002) 266–277. [10] F. Guichard, J.M. Morel, Partial differential equations and image iterative filtering, An Tutorial Int. Conf. Image Processing, Washington, DC, 1995. [11] P.K. Saha, J.K. Udupa, Relative fuzzy connectedness among multiple objects: theory, algorithms, and applications in image segmentation, Comput. Vis. Image Understand. 82 (1) (2001) 42–56. [12] L. Xinxiao, J. Li, Z. Nanning, Win-VolVis: an expandable visualization system in Win95/NT, J. Comput. Appl. Software of China (1) (2002).