Computer-aided detection of colonic polyps with level set-based adaptive convolution in volumetric mucosa to advance CT colonography toward a screening modality.

As a promising second reader of computed tomographic colonography (CTC) screening, the computer-aided detection (CAD) of colonic polyps has earned fast growing research interest. In this paper, we present a CAD scheme to automatically detect colonic polyps in CTC images. First, a thick colon wall representation, ie, a volumetric mucosa (VM) with several voxels wide in general, was segmented from CTC images by a partial-volume image segmentation algorithm. Based on the VM, we employed a level set-based adaptive convolution method for calculating the first- and second-order spatial derivatives more accurately to start the geometric analysis. Furthermore, to emphasize the correspondence among different layers in the VM, we introduced a middle-layer enhanced integration along the image gradient direction inside the VM to improve the operation of extracting the geometric information, like the principal curvatures. Initial polyp candidates (IPCs) were then determined by thresholding the geometric measurements. Based on IPCs, several features were extracted for each IPC, and fed into a support vector machine to reduce false positives (FPs). The final detections were displayed in a commercial system to provide second opinions for radiologists. The CAD scheme was applied to 26 patient CTC studies with 32 confirmed polyps by both optical and virtual colonoscopies. Compared to our previous work, all the polyps can be detected successfully with less FPs. At the 100% by polyp sensitivity, the new method yielded 3.5 FPs/dataset.

Introduction

According to the American Cancer Society statistics, colorectal carcinoma is the third most commonly diagnosed cancer and the second leading cause of death from cancer in the United States (see http://www.cancer.org/docroot/home/index.asp).1 It was estimated that 154,000 new cases would be diagnosed with 52,000 dying from the disease in 2007. Fortunately, most colon cancer probably arises from polyps, which can take 5–15 years for malignant transformation. Thus, early detection and removal of small polyps of less than 10 mm will effectively decrease the incidence of colonic cancer.2–4 As a new minimally-invasive screening technique, computed tomographic colonography (CTC) or CT-based virtual colonoscopy (VC) has shown several advantages over the traditional optical colonoscopy (OC),5 and an annual screening is recommended for people over 50 years old. In order to improve the sensitivity and consistency of CTC in the identification of polyps, computer-aided detection of polyps (CADpolyp) has been suggested as an effective second reader assisting radiologists for finding polyps in the colon.6–8

Typical polyps appear as protrusions out of the colon wall into the lumen. Therefore, in state-of-the-art CAD systems, the colon wall is firstly segmented from the original three-dimensional (3D) CTC image, followed by a shape analysis on the colon wall to generate initial polyp candidates (IPCs). Finally, various learning methods are used to reduce false positives (FPs).9

As the very first step of a CAD system, colon wall segmentation explicitly or implicitly determines the overall detection performance via the successive shape analysis. In this study, we employed a soft segmentation method10 to extract the colon wall as a volumetric mucosa (VM) with a thickness of several voxels, which is considered to be a complete representation of the colon wall without any artificial shrinkage and enlargement. Unfortunately, the relatively-thick VM often appears blurred and noisy, so we employed an anisotropic diffusion filter to smooth the VM. More details about the VM extraction and smoothing are referred to in the section on Volumetric mucosa extraction, below. Upon the extracted colon wall, polypoid-like initial candidates are determined by thresholding the geometric measurements, such as the mean curvature and sphericity ratio,11 the shape index and curvedness,12,13 the convexity or concavity,14 and the surface normal overlap.15 All these geometric measurements are functions of the first- and second-order spatial derivatives of the 3D CTC images. Monga and colleagues16 convoluted the 3D image volume with derivatives of Gaussian kernels (the 3D Deriche filters17) to compute the derivatives. Thirion and colleagues18 just used a central differencing scheme to calculate the derivatives but after filtering the 3D image with a Gaussian kernel. However, the robustness and effectiveness of both methods highly depend on the kernel window size (decided by the parameters of α1 and α2 in Monga’s method, and σ in Thirion’s method) they used. Therefore, Campbell and Summers19 explored the optimal α1 and α2 based on several phantoms, and they found that even with those optimal parameters, spurious computations could be observed on thin flat folds, caused by the interference of different topological structures residing in the same kernel window. To address such a problem, Huang and colleagues20 and Sundaram and colleagues21 built a geometric surface representation of the colon wall instead of the often used volumetric voxel representation, and shape analysis was conducted only on the surface. With the assumption that the surface is accurately extracted from the 3D CTC image, the interference between different topological structures can be avoided. However, Wijk and colleagues22 employed the normalized convolution to calculate the first- and second-order spatial derivatives, where the different structures were differentiated by using the gradient of a distance transform23 based on the air-tissue interface. The mathematically sophisticated method used by Wijk and colleagues22 is time consuming, and the hard segmented air-tissue interface cannot reflect the partial volume effect (PVE). Similarly, Wang and colleagues24 directly eliminated the contribution of voxels of different structures in the convolution by Deriche.17 However, the single-voxel interface layer (referred to as the starting layer (SL) in the following text) used to build up the distance transform was not properly selected either, and due to PVE, their SL would overlook the polypoid structures and incorrectly synthesize nonexisting polypoid structures. The details of this will be addressed in the level set-based adaptive convolution (LSAC) method discussed below. In this study, we employed a LSAC method to calculate the partial derivatives. As such, a more accurate and robust SL was retrieved.

Geometric features aforementioned to characterize colonic polyps have been extensively presented towards CADpolyp. However, most of them were extracted based on a single layer mucosa representation. Necessity still exists for exploring the correspondence among different layers of the mucosa/colon wall in CTC images. Wang and colleagues25 introduced “global” versions of the shape index and curvedness, within spatial domain, by linear smoothing along the principal directions. Nevertheless, their method ignored the correspondence between different layers of the mucosa/colon wall. Therefore, Wang and colleagues24 smoothed the first and second principal curvatures along the normal direction, referred to as weighted integral curvature along normal directions (WICND). In their method, the voxels on the middle-layers of the VM were treated equally as those voxels on the other layers. However, we argue that the middle-layers of the mucosa/colon wall appear smoother, and should be more reliable than the other ones when conducting shape analysis. In this study, we introduced a middle-layer enhanced integration (MEI) strategy to refine the geometric features, where higher weights were put on the voxels at the middle layers of the VM.

The remainder of this paper is organized as follows. In Methodology, we firstly provide a brief overview of the proposed CAD scheme, and then the aspects of the pipeline are demonstrated, especially the LSAC method. We then describe the performance evaluation and give the results. Finally, some conclusions are drawn.

Methodology

Overview

The whole CAD pipeline can be outlined in four stages: VM extraction, IPC detection, FP reduction and detection presentation, as shown in Figure 1.

Figure 1

Pipeline of the presented CADpolyp system.

Abbreviations: CTC, computed tomographic colonography; FP, false positive; IPC, initial polyp candidate; LSAC, level set-based adaptive convolution; SVM, support sector machine; VM, volumetric mucosa.

In the first stage (to be detailed in Volumetric mucosa extraction, below), through the MAP-EM (maximum a posteriori expectation maximization) soft segmentation algorithm,10 the mixture percentage of each tissue type inside each voxel in the CTC images was estimated, from which the VM was extracted and smoothed afterward. In Stage 2 (to be detailed in The level set-based adaptive convolution (LSAC) method, The middle-layer enhanced integration (MEI) strategy, and Polyp candidate extraction, below), we firstly applied the LSAC method inside the VM to conduct the geometric computation. The computed geometric measures were then smoothed by the MEI strategy and used to extract suspicious patches. The suspicious patches were fed to the filter25 to remove noise-induced patches. In what follows, the survived patches grew into the tissue area25,26 to form the IPCs. In the third stage (to be detailed in FP reduction, below), both shape and textural features were extracted for each IPC, and fed into a weighted SVM model27 to reduce FPs. In the last stage, the detection results were displayed in a CTC commercial system and presented to physicians to serve as second opinions (to be detailed in Detection presentation, below).

Volumetric mucosa extraction

For any imaging modality (eg, CT), PVE always exists. Considering the use of positive-contrast tagging agents to opacify the residual fecal for differentiation of the materials from colon wall and polyps,28 PVE became severe and the thickness of the VM varied dramatically. In the case of extracting mucosa for CAD in CTC, a sheer surface ignoring PVE, as segmented by most of the previous colon-wall segmentation algorithms, undoubtedly shrunk or even distorted the true mucosa without any chance of recovery. Figure 2 fully illustrates the limitation of such hard thresholding segmentation without considering PVE. Although by locally adjusting the threshold for a better colon wall representation with less possibilities of shrinkage or enlargement (ie, by a particular threshold value), it was still data-oriented and hardly to be generalized to all the cases. More specifically, the over-fitted colon wall12 may include noncolonic structures with increased FP rates, as shown in Figure 2d. Furthermore, a large amount of unwanted air and tissue areas may also be included, resulting in an increased computational burden.

Figure 2

Segmented colon walls by traditional hard segmentation algorithms: a) part of an original CTC image slice, the red arrow indicates the blurring area near the air-colon interface reflecting the PVE and the yellow arrow flags a 9 mm polyp; b) the colon wall (the cyan curve) in another slice near the same polyp in a) was created by region growing to −300 HU,11 which cuts off part of the polyp in light gray indicated by the yellow arrow; c) the colon wall (the cyan curve) of the same slice in a) was created by region growing to −800 HU,11 where the polyp gets enlarged near the yellow arrow; and d) a thickened colon wall (the area between the two green curves) was created by morphologically dilating the iso-surface at −650 HU (the cyan curve) inward and outward by 5 mm,12 where the yellow arrow indicates that part of the small bowel is incorrectly included.

Abbreviations: CTC, computed tomographic colonography; PVE, partial volume effect; VM, volumetric mucosa.

To address the PVE, a soft image segmentation algorithm10 was developed. The algorithm aims to estimate the tissue mixture percentages inside each voxel within the framework of iterative MAP-EM estimation. By viewing the tagged fecal material as an extra class type, its mixture fractions inside each voxel could be accurately estimated and cleansed by assigning them to those of air class. The output of the MAP-EM algorithm includes: (1) an electronically cleansed image in CT densities (Figure 3b); and (2) the colon lumen (air class) percentage distribution map in the range of [0, 1], , to reflect the PVE (Figure 3c). Finally, the VM with variable thickness was obtained by directly extracting those voxels whose colon percentages are between 0 and 1, which means that they are not fully occupied by either tissue or air (Figure 3d). Compared to the sheer surfaces of Figure 2b and Figure 2c via the traditional hard thresholding segmentation, the extracted VM of Figure 3d fully addressed the PVE without any loss of texture information, and represented the colon wall without any biased shrinkage or enlargement. While compared to the artificially-thickened colon wall of Figure 2d, the extracted VM did not introduce any noncolonic structures, ie, it was more compact (saving computational time).

Figure 3

The process of VM extraction. The red arrow flags an 8 mm polyp. a) Part of an original CTC slice. b) Cleansed image. c) Colon lumen. d) The extracted VM reflecting the PVE.

Abbreviations: CTC, computed tomographic colonography; PVE, partial volume effect; VM, volumetric mucosa.

To smooth the noise and mitigate the blurring inside the VM, we applied a limited anisotropic diffusion (LAD) strategy.29 In this method, the anisotropic diffusion filter30 was applied inside the VM in the colon lumen image , as shown in Figure 3d. The innermost and outermost layers of the VM were taken as the boundaries of the diffusing process, and the Neumann boundary condition31 was applied when solving the diffusion equation at the boundaries. Therefore, the diffusion process was limited just inside the VM. Compared to those fixed-sized kernel window filters, the advantage of such strategy is that the smoothing process is adaptive to the local topology of the VM, so that the interference of different geometric structures would be reduced or avoided in this stage. Based on such smoothed VM, we conducted the following geometric analysis only on the voxels in the VM.

The level set-based adaptive convolution (LSAC) method

It is widely accepted that polyps appear as protrusions out of the colon wall. In previous studies, geometric measures like the mean curvature and sphericity ratio,11 the shape index and curvedness,12,13 the convexity or concavity,14 and the surface normal overlap,15 were explored to characterize such appearance. Interestingly, all these geometric measures are functions of the first- and second-order spatial derivatives of the 3D CTC images. Specifically,16 a convolution procedure was employed, ie, I = (f1 (i) f0 (j) f0 (k))*I for the first-order derivatives, while I = (f2 (i) f0(j) f0(k)) * I and I = (f1 (i) f1(j) f0 (k)) *I for the second-order and mixed partial derivatives, respectively. Here, I represents the image intensity, and i, j and k indicate the three dimensions of the 3D space. The operators f0(x), f1(x) and f2(x) are defined as where c0, c1, and c2 are constants for normalization, and parameters α1 and α2 determine the convolution window. Figure 4a shows a 31 × 31 kernel window in 2D case, which can effectively smooth the CTC image and get meaningful geometric measures. Unfortunately, the colons are densely structured, and different topological structures may be located in one neighborhood area, as shown in Figure 4b. The reliability of the computed spatial derivatives is impaired, and so forth the geometric analysis, which is exactly the reason why spurious computations were observed on thin flat folds.19

Figure 4

a) A 31 × 31 neighborhood of the red pixel when α1 was set to be 0.719 with truncation value to be 0.00005. b) The zoomed image in the neighbor window of the voxel at the red cross position indicating the 8 mm polyp in Figure 3. Voxels A and B belong to another structure and they are interfered with the voxel at the red cross position.

The distance transform maps each image pixel into its smallest distance to a given region of interest (ROI),23 which can help to distinguish the different structures in the local neighborhood. As shown in Figure 5, the object in cyan is encompassed by a thick curve, ie, the boundary of ROI. After distance transform, the number on each pixel indicates the square of the smallest distance from the pixel to the thick curve. The arrows at pixels M, N, and P point to the directions of their distance mapping-based gradients. Obviously, we can claim that if the gradient directions at any two pixels differ significantly, like those at M and N, these two pixels belong to different structures. However, the opposite argument does not always hold, ie, any two pixels having slightly-varied gradient directions are not guaranteed to originate from the same structure, like pixels M and P in Figure 5. Wijk and colleagues22 quantified this by the confidence value of 0 or 1 in a diagonal matrix in their normalized convolution method to address the interference between different structures. However, we noticed that the above criterion of differentiating structures can be adaptively incorporated in the widely used convolution method16 by simply modifying the convolution weights as where l varies among 0, 1, and 2, and f1(x) is formulated in Equation 1. f1(x) denotes the modified version of f1(x). Notation angle (c, x) evaluates the angle between the center pixel c and any neighboring pixel x, in terms of the distance mapping-based gradient directions. 0 is the predefined thresholding angle, like 90 degrees. Therefore, most pixels in the convolution window but belonging to other structures are excluded from center pixel c in the following convolution operation. However, if the convolution window is too big, some farther pixels will still contribute undesirably to the convolution although they belong to different structures from center pixel c. For example, pixel P will contribute to the convolution operation centered at pixel M, if the distance between pixels P and M is shorter than half of the convolution window size. Fortunately, according to Equation 1, a longer distance from the center pixel will lower down the weights exponentially, and the adverse contributions to the convolution will be effectively reduced. Therefore, the robustness of such adaptive convolution method can be guaranteed, and its capability of relieving the interference of different structures has been demonstrated by a phantom study.24

Figure 5

The distance transform used to distinguish different topological structures.32

In this paper, we refer to the thick curve in Figure 5 as the starting layer (SL) of the distance transform. The overall performance of differentiating different structures in terms of the derived distance field is largely determined by the specified SL. However, the previous work22 segmented the air-tissue interface, serving as the SL, simply by using a thresholding method and ignoring the PVE. While the soft image segmentation algorithm used by Liang and Wang10 was later employed by Wang and colleagues24 and led to a volumetric mucosa to represent the colon wall more accurately, the SL was not properly selected yet, which was just the mucosa-tissue interface as shown in Figure 6b by the red dotted curve. The major problem of this mucosa–tissue interface is its failure to detect the 10 mm polyp due to the severe PVE, along with the remaining interference near the polyp. An alternative layer, ie, the mucosa–air layer, might serve as the SL as well. However such layer, as illustrated in Figure 6b by the blue dotted curve, may also deform the colon wall in the sense that the concave structure at the neck of the polyp, as blue arrow indicated, has been artificially distorted towards an opposite direction. In this paper, we are motivated to introduce a level set-based adaptive convolution (LSAC) method, where a better SL, reflecting the colon wall more accurately, is generated to improve the adaptive convolution.

Figure 6

a) One slice of a supine image, where the red arrow indicates a 10 mm tubular adenoma. b) The zoomed lumen image near the 10 mm polyp (pointed with the red arrow). The red and blue dotted curves represent the mucosa-tissue and mucosa-air layers respectively, while the green dotted curve is the new SL generated with the level set method in this paper.

Abbreviation: SL, single layer.

Level set framework is one of the geometric deformable models and is based on the curve evolution theory.33 Compared with other methods, it is able to combine region-based and edge-based information together, make use of global and local information simultaneously and control the geometric property of level set function easily. We intentionally introduce the level set method to retrieve a better SL, from which we build the distance transform to distinguish different topological structures.

Considering the PVE, the SL interface best depicting the shape of colon wall is absolutely not the outermost layer (the mucosa–tissue layer, the red dotted line in Figure 6b) or the innermost layer (the mucosa–air layer, the blue dotted line in Figure 6b) of the VM. Intuitionally, it is located in between the outermost and innermost layers, where the variation of CT intensities across the different layers remains relatively stable. Therefore, we extended the C–V model34 to multi-scale space, where the global and local intensity distribution properties are taken into consideration. Specifically, the gradient of image intensity is used to construct the evolution stopping rule, so that where φ is the Lipschitz function, and I represents the image intensities. The two superscripts in and out indicate the regions where φ > 0 and φ ≤ 0, respectively, while Ī, Ī, and I represent the mean intensity values of the whole image, voxels in the narrow band and the local neighborhood respectively. The notations λ, μ0, μ1, μ2, and μ3 are controlling parameters, and ∇ represents the gradient operator. is the curvature of Lipschitz function, which determines the smoothness of the zero level set surface.

Once Equation 3 converges, the resulting zero level set surface of φ = 0 indicates a layer between the outermost and innermost layers, with slightly-changed CT intensities coming across different layers. As described in Figure 6b by the green dotted curve, the acquired zero level set surface significantly improves the geometric representation of the colon wall.

The middle-layer enhanced integration (MEI) strategy

After applying the above LSAC strategy and the curvature computation method,16 geometric measures, like the gradient, the first and second principal curvatures, can be retrieved for each voxel in the VM. Similar to our previous work,24 we further smooth the geometric measures ci, eg, the principal curvatures, via the following convolution operation where L is the path starting from the concerned voxel along the positive and negative gradient directions inside the VM, as shown in Figure 7. The Gaussian distributed weight W(x) was defined as , x ∈ L with the standard deviation formulated as , where a and b are predefined constants. The term represents the ratio of the local smoothness of to the global one along the path L, whose details are referred to by Wang and colleagues.24 However, following the same arguments as we acquired SL above, we believe the middle layers (eg, layers 2 and 3 in Figure 7) of the VM would be more reliable than the other layers (eg, layers 1 and 4 in Figure 7) for evaluating the geometric features of the colon mucosa. Therefore, we introduced the MEI strategy to highlight the middle layers of the VM through the term

Figure 7

The integral path across the layers of the volumetric mucosa.

Taking Figure 7 as a typical example with totally four layers inside VM and starting from layer 1, the breadth-first dilation strategy marches towards layers 2, 3, and 4. For the presentation simplicity, we simply labeled the dilation step as the layer number when attached to each voxel inside the VM (as shown in Figure 7). Once all the layer numbers are determined, the weight is formulated as: where n is the layer number of voxel on the path, while n̄ is the mean value of the layer numbers of all voxels on the path. Voxels with layer numbers close to n̄ are those in the middle of the path, and therefore they have smaller values of in Equation 5, smaller standard deviations σ in the Gaussian model, and finally higher weights of W(x) in Equation 4.

After applying Equation 4 by means of the MEI strategy, the first and second principal curvatures were refined and denoted by and , respectively, from which the shape index SI (indicates the type of the surface) and curvedness CV (indicates how curved the surface is) were computed35 where . One step further, by reinterpreting the geometric measures c in Equation 4 as SI and CV respectively, the smoothed version, denoted as SI and CV, are acquired to further exploit the correlations between different layers in the VM.

Polyp candidate extraction

As assumed in the previous CAD reports,11–15,25,26 the upper part of the colonic polyps usually appears as a “spherical cap” bulging into the colon lumen. However in our CAD pipeline, the colon lumen image (ie, the percentage distribution map in the range of [0, 1]) views the polyps oppositely. Therefore, the colonic polyps appear as “spherical cups” in the lumen image. Given SI and CV as obtained in The middlelayer enhanced integration (MEI) strategy, the variations on the colon wall with local “spherical cup” shapes can be detected by the use of a simple clustering rule.12 However, the geometric measures, like the shape index and curvedness, are very sensitive to the small variations on the colon wall. Usually, more than one hundred suspicious patches would be generated including a great number of nonpolyp ones. To eliminate the nonpolyp patches as many as we can in this stage, we followed our previous work25 to remove those patches whose sizes were smaller than the clinicallyimportant colonic polyps (eg, larger than or equal to 5 mm in diameter).

Colonic structures, like feces and prominent folds, mimic the spherical shape of true polyps. It’s difficult to identify the colonic structures from polyps by using only geometric information. To deal with such cases, analyzing the internal textural features for each detection would be an effective way.25 For the purpose of conducting textural analysis, a volume of interest (VOI) is extracted for each suspicious patch. Because of the small contrast, it’s hard to identify the border between a polyp and its surrounding normal tissue. Wang and colleagues25 assumed that ellipsoid was a choice to approximate the shape of polyps. However, the shape of the outer part of polyps bulging into the colon lumen varies greatly, like sessile, pedunculated and flat shapes. Therefore, the assumption of elliptical region of interest (eROI)25 might be violated in some cases and consequently include unwanted voxels. Nappi and colleagues26 introduced an iterative conditional morphological dilation method to extract the VOIs. This method proceeded with morphological dilation, which terminated if a minimum growth rate was detected or a predefined dilation distance was reached. Advantages of this method lie in the real reflection of the outer border of polyps. It works fine for pedunculated polyps, as the minimum growth rate can often be detected under this circumstance. Otherwise, the morphological dilation process would reserve the morphological characters of the starting suspicious patch, which frequently happens to the sessile and flat lesions. As shown in Figure 8c, the inner border of the extracted VOI is morphologically similar to the starting region in Figure 8b and still appears as concave, which is generally contradictory to the polyp pathology in general.36 Therefore, in this study, we utilized the inner border generated25 to stop the dilation process when the minimum growth rate cannot be detected. The improved result is shown in Figure 8d.

Figure 8

The VOI extraction for a 6 mm sessile polyp in a sigmoid colon. a) The original polyp (pointed by the red arrow). b) The suspicious patch (the red area). c) The VOI generated with the conditional morphological dilation method.26d) The VOI from our new method.

Abbreviation: VOI, volume of interest.

The grown VOIs are called the IPCs, from which rich textural information can be extracted for the purpose of FP reduction.

FP reduction

For FP reduction, numbers of features were extracted for each IPC to form an associated feature vector. Once we collected all the feature vectors as a sample pool, they were randomly divided into two groups, the training and the testing samples. Given ground truth of the training samples (true positive or false positive, ie, TP or FP), supervised learning machine like LibSVM27 with nonlinear RBF (radial basis function) kernel was trained for maximizing the overall classification accuracy. The classifier was then applied to the testing samples for prediction. In our study, the sample pool was heavily imbalanced, ie, the number of FPs was significantly larger than that of TPs, such that the cost function was overwhelmed by the contribution of FPs, ie, LibSVM favored the decrease in FPs without ensuring the increase in TPs. Zheng and colleagues37 reported that the resulting sensitivity was always less than 0.5, which was very poor for a CAD scheme.

To address the above dilemma, two modifications were made based on the LibSVM package:

Instead of maximizing the overall classification accuracy, we modified the training goal as minimizing the number of FP at a specified detection sensitivity of TP, denoted as ζ. By varying ζ (sensitivity levels), we investigated the classification potentials of different feature combinations through fROC (free-response operator characteristics) analysis.

Those training samples of TP were assigned larger weights to modulate the imbalanced dataset. In the LibSVM package,27 the training process was actually a procedure of 2D grid search optimizing two parameters to meet the specified goal, ie, the cost function value C and the RBF kernel parameter γ. We extended such process to a 3D grid search by adding another search for an optimal weight value for TPs. Zheng and colleagues37 assigned a larger value to the parameter C for TPs and a smaller one for FPs to balance the datasets, sharing a similar goal as ours but by a different approach.

In order to evaluate the performance of LSAC and MEI strategies in CAD applications, we simply employed the four well established geometrical and textual features by Nappi and colleagues,26,38 which were the mean and variance of shape index, the variance of curvedness, and the variance of CT intensity respectively. They were stacked together to form a feature vector for each IPC and fed into the SVM learning package. The output of the SVM would reflect the gain by the use of the LSAC and MEI strategies.

For the purpose of CAD methodology development, three more aspects would be further considered, in addition to the above presented strategies in steps 1–3 of Figure 1. One is the extraction of more features from the VOI and selection of useful ones for the purpose of detection.39 The second one is the enhancement of classifier for improved detection by using the selected features.39 The third one is the presentation of the CAD detection. In what follows, we give a brief narrative regarding the use of CAD to boost radiologist’s detection performance on polyps everywhere throughout the entire colon.

Detection presentation

The SVM-reported final polyp candidates are delivered to a commercially available colon visualization system (ie, the V3DColon Module; Viatronix Inc, Stony Brook, NY). The system is dedicated to realize the concept of virtual colonoscopy (VC) by mimicking the optical colonoscopy (OC) 3D navigation procedure. Once locating a polyp candidate in such virtual environment, a 3D endoscopic view pops out right in the center screen, along with three cross-sectional views displayed on the right side one after another, as depicted in Figure 9. As a strong supportive module, the CADpolyp pipeline was integrated into the visualization system. Radiologists would then pay more attention to the CAD detected locations and spend less time on other areas of the colon mucosa. In other words, the integrated CAD and visualization system plays the role of computer assisted detection of polyps by human experts.

Figure 9

Presentation of a polyp candidate, which is a 35 mm sessile polyp in rectum according to the VC and OC reports.

Abbreviations: OC, optical colonoscopy; VC, virtual colonoscopy.

Performance evaluation

CTC database

To evaluate the performance of the presented CADpolyp pipeline of Figure 1 with the newly developed LSAC and MEI strategies (see the sections on The level set-based adaptive convolution (LSAC) method and The middle-layer enhanced integration (MEI) strategy, above), VOI extraction (see the section on Polyp candidate extraction, above), and modified SVM parameter setting (see the section on FP reduction, above), twenty-six CTC patient studies, aged from 50 to 80 years old, which were acquired from the University of Wisconsin Hospital and Clinics from July 2006 to April 2007, were used. Each study includes two CT scans, one at supine and the other at prone positions, resulting in 52 CTC datasets. A multi-slice CT scanner (Light Speed Ultra; GE Medical Systems, Milwaukee, WI) was used in helical mode to collect data. The scanning protocol includes modulated mAs in the range of 120–216 mA with kVp of 120–140 values.40 A total of 32 clinically significant polyps were confirmed with both optical and virtual colonoscopies: nine polyps sized in the range of 5–8 mm and 23 larger than 8 mm with the largest one being 35 mm. In this evaluation, the supine and prone scans of each patient study were considered as different datasets, and the detection sensitivity was evaluated by per polyp and FP rate by dataset (per image). Therefore, there were 64 polyps in 52 CTC images in the CTC database.

Evaluation of initial polyp candidates (IPC) detection

We set (α1, α2) = (1, 0.8) in Equation 1. After applying LSAC method and MEI strategy, those voxels with SI in [0, 0.06] and CV in [0.05, 0.14] were labeled as seed, and those in [0, 0.25] and [0.02, 0.5] were identified as growable voxels. Via the above empirically-set thresholds, we could make sure that all the polyps can be initially detected successfully in this stage. Seed voxels were firstly clustered together if they were spatially six-connected, grew into suspicious patches by including the six-connected growable voxels. By imposing constraint on the minimal surface area of the found patches, ie, 15 mm2,25 the surviving suspicious patches penetrated into the surrounding tissue area to form IPCs. The locations of the ground truths were determined by physicians, which were represented with the slice numbers in axial, saggital, and coronal directions. These locations were used to specify if a candidate is a TP or FP. If all the 64 polyps were detected successfully and appeared in the resulting group of IPCs, we say that the detection sensitivity in this stage reaches 100% (ie, the by-polyp sensitivity), and the number of false negatives (FNs) is zero.

Evaluation of overall performance of the CADpolyp pipeline

Four features were extracted for each IPC, and stacked together to form a feature vector. The feature vectors of all IPCs were randomly split into two groups of the same size having the same number of TPs, serving as training and testing set respectively. When we conducted the 3D grid search for the optimal parameters C, γ, and ω mentioned in FP reduction, above, the sensitivity level ζ served as the sweeping variable, and varied from 0 to 1. The final performance was measured based on the testing process. In the testing set, there were 32 true polyps. If all the 32 polyps were classified as TPs, we say that the final detection sensitivity was 100% and there was no FN; while if all the 32 polyps were classified as FPs, the final detection sensitivity would be 0%, and the number of FNs would be maximally 32. However, we could not tell how many datasets there were, since the testing candidates were randomly selected from all the IPCs. In order to calculate the number of FPs per dataset, it’s reasonable to assume that there were 26 datasets in the testing set. As a result, the final performance was presented by using fROC curves.

Results

Performance of IPC detection

Under the pipeline shown in Figure 1, three CAD schemes of detecting IPCs were applied to the CTC image database mentioned earlier: traditional method, previous or old adaptive method and the presented LSAC method. In the traditional method (or the first scheme), we just employed the isotropic Gaussian filter to compute the first and second derivatives of the images. In the previous adaptive method (or the second scheme), the adaptive method24 was used. In the third scheme, the presented LSAC method was applied. Based on the first and second derivatives of the three schemes, the successive steps with the same parametric configuration in the pipeline of Figure 1 were applied to extract the IPCs. The IPC extraction performance of these three schemes was shown in Table 1, where the detection sensitivities were evaluated based on datasets. The number of FP per dataset (NFPPD) and the FP reduction percentage (FPRP) based on the traditional scheme were listed in the last two columns.

Table 1

Performance of the CAD schemes

	Sensitivity	NFPPD	FPRP
Traditional method	100%	33.4	–
Previous adaptive method	96.9%	19.0	43.0%
LSAC method	100%	16.4	51.0%

Abbreviations: CAD, computer-aided detection; FPRP, false positive reduction detection; NFPPD, number of false positives per dataset.

Compared to the traditional method, the two adaptive schemes significantly reduced the number of FPs (ie, 43% and 51% respectively). However, two polyps, one of which was 10 mm as shown in Figure 6 and the other was a 6 mm tubular adenoma in a prone image, were missed by the previous adaptive method, because the SL used by Wang and colleagues24 incorrectly depicted the real shape of the colon wall near the polyps and failed to recognize the topological structures there. Such problem was rectified by our new LSAC method and even more FPs were reduced.

Performance of overall CADpolyp pipeline

For the three schemes mentioned in the above section, the generated feature vectors of all the IPCs were fed into the same SVM classifier for the purpose of FP reduction. To evaluate the overall CAD detection performance, three fROC curves in green, blue and red, corresponding to three detection schemes respectively, were generated in Figure 10 as described in the evaluation of overall performance of the CADpolyp pipeline, above, and compared with each other from two different aspects. First of all, the two adaptive methods (red and blue curves) noticeably outperformed the traditional method (green curve). Secondly, the old adaptive method (blue curve) missed two clinically significant polyps in the stage of IPC detection, and consequently it could not achieve 100% by polyp detection sensitivity. The presented LSAC method (red curve) solved this problem, and more FPs were reduced. More specifically, the traditional method yielded 100% sensitivity with 5.9 FPs/dataset, the previous adaptive method reached 96.9% sensitivity with 3.7 FPs/dataset, and the presented LSAC method achieved 100% sensitivity with 3.5 FPs/dataset.

Figure 10

fROC curves generated by the three schemes. The green, blue and red curves represent the overall detection performances of the traditional, the previous adaptive and the new LSAC method respectively.

Abbreviations: fROC, free-response operator characteristics; LSAC, level set-based adaptive convolution.

Discussion and conclusion

In this study, we introduced a LSAC strategy to launch a more meaningful adaptive convolution method, which led to a more accurate geometric analysis on CTC images. Such an improved geometric analysis greatly benefited the CADpolyp scheme of Figure 1. In our experiments, the missed polyps in the IPC detection stage by using the previous method24 were successfully detected with the proposed adaptive convolution strategy, and the number of FPs per dataset was reduced from 19.0 to 16.4, which was a 13.95% reduction. After applying the SVM classifier, our new CAD scheme yielded 100% sensitivity with 3.5 FPs/dataset.

This paper mainly focused on introducing the LSAC strategy to improve the computation of the first- and secondorder derivatives, which was usually the starting point of the geometric analysis. Then, the MEI strategy was applied to exploit the advantage of our volumetric mucosa which was derived from the soft segmentation algorithm.10 The advantages of the presented CAD scheme (utilizing both LSAC and MEI strategies) were demonstrated with four established features26,38 together with the well-known SVM classifier. The CAD scheme could be further improved by considering more useful features and using a more powerful classifier.39

With the above excellent result of 3.5 FPs/dataset at 100% by polyp sensitivity, we noticed that in the final CAD detections, there were still some removable FPs, eg, tubes and floating feces. These survived FPs of CAD could be recognized by experts, eg, radiologists, by using a sophisticated visualization system, eg, the V3D Colon Module system. Therefore, in this paper, we presented an integrated system of CADpolyps and the visualization module which would help the experts further classify the TPs and FPs (see the section on Detection preventation, above). More investigation and validation on the integrated system are under progress.

Inheriting the same assumption regarding the polyps’ spherical-shaped protrusions as the previous CAD schemes did,12,13,25,38 our CAD pipeline also has lower detection rate, as the conventional colonoscopy,12 for irregularly-structured polyps, like those flat ones. A more generalized model, which is applicable to almost all kinds of clinically significant polyps, deserves our great efforts in the future. The gain of the integrated system for detecting flat polyps will be another research topic.

We recognize the limitations of our small-sized database with confirmed and spatially-located polyps. And we are currently devoting ourselves into the access to more datasets. Experienced physicians are invited to locate the (x, y, z) positions of more polyps in our CTC image database. It is expected to conduct a more thorough and statistically-significant evaluation of the presented CAPpolyp pipeline.