Patent Application
20050152588
Field of the Invention
The present invention relates generally to systems and methods for the computer-aided detection of three-dimensionally extended organ lesions. The present invention also generally relates to automated techniques for the detection of abnormal anatomic regions, for example, as disclosed, in particular, in one or more of U.S. Pat. Nos. 4,907,156; 5,133,020; 5,832,103; and 6,138,045; all of which are incorporated herein by reference. Further, the present invention relates to visualization methods for virtual endoscopic examination of CT colonographic data.
The present invention also generally relates to computerized techniques for the automated analysis of digital images, for example, as disclosed in one or more of U.S. Pat. Nos. 4,839,807; 4,841,555; 4,851,984; 4,875,165; 4,907,156; 4,918,534; 5,072,384; 5,133,020; 5,150,292; 5,224,177; 5,289,374; 5,319,549; 5,343,390; 5,359,513; 5,452,367; 5,463,548; 5,491,627; 5,537,485; 5,598,481; 5,622,171; 5,638,458; 5,657,362; 5,666,434; 5,673,332; 5,668,888; 5,732,697; 5,740,268; 5,790,690; 5,832,103; 5,873,824; 5,881,124; 5,931,780; 5,974,165; 5,982,915; 5,984,870; 5,987,345; 6,011,862; 6,058,322; 6,067,373; 6,075,878; 6,078,680; 6,088,473; 6,112,112; 6,138,045; 6,141,437; 6,185,320; 6,205,348; 6,240,201; 6,282,305; 6,282,307; 6,317,617; 6,335,980; 6,363,163; 6,442,287; 6,466,689; 6,470,092; 6,594,378; 6,738,499; 6,754,380; 6,813,375; as well as U.S. patent application Ser. Nos. 09/759,333; 09/760,854; 09/773,636; 09/816,217; 09/830,562; 10/120,420; 10/270,674; 10/292,625; 10/358,337; 10/360,814; all of which are incorporated herein by reference.
The present invention includes use of various technologies referenced and described in the above-noted U.S. Patents and Applications, as well as described in the references identified in the following LIST OF REFERENCES by the author(s) and year of publication, and cross-referenced throughout the specification by reference to the respective number, in brackets, of the reference:
Colon cancer is the second leading cause of cancer deaths in the United States [41,42,63]. Colonic polyps are seen as potential precursors of colon cancer, and studies show that early detection and removal of the polyps can reduce the risk and the incidence of colon cancer [1,2,43]. Virtual endoscopy of the colon, often called computed tomographic colonography (CTC) or virtual colonoscopy, is a new technique for non-invasive detection of colorectal neoplasms by use of a CT scan of the colon [2,3,4,64]. If CTC can demonstrate satisfactory performance in the detection of polyps, it could be an effective technique for screening of colon cancer in large populations [4,5,65].
To be a clinically practical tool for screening of colon cancer, CTC must be feasible for the interpretation of a large number of images in a time-effective fashion [6,7], and it must facilitate the detection of polyps with high accuracy [4]. However, at present, the interpretation of CTC is a potentially slow process. Despite recent advances in the development of 3-D imaging software and hardware, the current interpretation time of CTC is at least 5-20 minutes per case even for expert abdominal radiologists, and human factors such as fatigue may lengthen the interpretation time, particularly when a primary 2-D interpretation technique based on multiplanar reformatting is used [2]. Moreover, the reported diagnostic performance of CTC varies substantially among readers and among different studies. Although several studies have yielded by-polyp detection sensitivities of 70-90% with corresponding by-patient specificities of approximately 90-100% [8,44], some studies have reported much lower performance results [9].
The accuracy of the detection of polyps is often affected by image display methods because they change the visibility and conspicuity of the polyps. Therefore, a suboptimal implementation of the display methods may increase the perceptual error even among experienced readers [10]. To facilitate a rapid and accurate interpretation of CTC, a number of computer-assisted visualization methods have been reported during the past several years [2,11]. Most of them are concerned with the visualization of the entire colon, so that observers do not miss any lesions. Volume-rendered 3-D endoluminal viewing is a widely accepted method for simulating the appearance of optical colonoscopy [12], and multiple windows may be used for presenting a full 360-degree view of a colonic lumen at a given location [11]. The colon may be unfolded to represent it as a 2-D plane [13,14], which could potentially minimize the blind spots that are often associated with an endoluminal fly-through view [15]. Another method, based on volume rendering, is a semi-transparent view of the complete colon, similar to a double-contrast enema, that can be used for displaying the locations of the findings for a follow-up examination [16].
In these visualization methods, the standard practice has been to color the colon by use of a monochromatic coloring scheme, such as a grey-scale or pink-scale coloring based on the CT values within the data set. A careful, and therefore time-consuming, interpretation of the CTC data is required, because the visual color cues that normally exist in colonoscopy, such as mucosal color changes, are absent. For example, as shown in
Pseudo-coloring (or indexed coloring) is a well-known method for enhancing the visualization of data that do not have inherent color information [17,18]. It can be used, for example, to highlight important values or to identify structures, features, or patterns within the data. The input data are typically colored by substituting each data value by a color indicated by a lookup table. Also in CTC, some previously reported investigations have used such color schemes for the visualization of colonic lesions. For example, Choi et al. [19] used a color spectral imaging scheme to enhance the visualization of axial CTC images. Summers et al. [20] visualized the results of a computerized polyp detection scheme by coloring of the colon based on the curvature of the colonic surface. Wax et al. [21] developed a translucent volume rendering method for the analysis of mucosal opacities, which allowed one to differentiate between contrast-enhanced stool and polyps. Others have developed color schemes for the visualization of the texture differences between polyps and stool [22,23].
Because the diagnosis of a CTC case is made by radiologists based on visual examination, an unambiguous display of the colon is expected to improve the reader performance in CTC.
Moreover, although recent advances in virtual endoscopy allow physicians to examine the colon in CTC interactively with endoluminal perspective views [66], complete navigation through the entire colon with virtual endoscopy can be tedious and time consuming because an inexperienced physician may lose track of the position and orientation of the colon. The difficulty and the time required for navigating through the colon could be reduced substantially by use of the centerline of the colonic lumen. Directing a 3-D camera that creates the endoluminal views of the colon along the centerline allows physicians to navigate easily through the colon. A physician can also veer from the path for closer inspection of the anatomy and later use the central path to re-orient him-or herself in the virtual space. The colon centerline can also aid an advanced visualization of the colon, such as straightening the colon [13,67] and the cylindrical and planar-map projections of the colon [68].
Several approaches have been proposed for computation of the centerline of a colon, some of which are semi-automatic, whereas others are fully automatic. Ge et al. [69, 70] used topological thinning for creating a 3-D skeleton, then removed loops and branches and, finally, approximated the centerline by using B-splines. A user defines a starting seed point. Paik et al. [71] generated an initial path at the object surface by using region-growing. The path was later refined to a center position by use of a thinning algorithm. The user supplied the program with two end points. Samara et al. [72] used a seed point in the rectum as the start for a region-growing that generated a centerline. The same procedure was repeated for another seed point placed in the cecum. The final centerline was the average of the two centerlines. Wan et al. [73] used a distance from boundary field to create a directed, weighted graph. The centerline was calculated by computation of the minimum-cost spanning tree for the graph. One or two pre-defined seed points were needed as input. Bitter et al. [74] computed the gradient field of a distance from boundary field, which was first used for computing the distance from the source and the root in order to find the end points, and then to find the actual centerline. This centerline could later be expanded to a skeleton based on the same basic principle. The method required a single connected component as input, but no seed points were necessary. Deschamps and Cohen [75] used one or two user-defined seed points as a starting (ending) point for a front propagation that was used for finding an initial distance map. A second, inwardly propagated front was then computed, and a minimal path search for this map was performed which provided the final centerline.
Most of these existing methods rely on manually extracted seed points or complete colon segments identification for the computation of a single centerline segment. Manual extraction of the seed points can be a tedious and inaccurate process when the location of the colon is obscured by a large amount of small bowel adhering to the colon, or when the colon is fragmented into multiple disconnected segments due to collapsed regions. For the same reasons, the segmentation of the colon is often incomplete, and it contains extra-colonic segments in addition to the colon. This makes methods that rely on a complete colonic pathway and complete colon segments identification work poorly unless some additional seed selection procedure is applied for finding the different parts of the colon. Fully automated methods can be used even in cases with collapsed regions by iterating through each independently segmented region. However, in this way, centerline segments will also be created for each extra-colonic findings.
The existing algorithms for the computation of a colon centerline are also time-consuming. At best, a computation time of 12-17 seconds per centerline, evaluated on 4 cases with use of an Intel Pentium-based 700 MHz PC, was reported [73]. This computation time excludes the time for the Euclidean distance transform (DT) [76] needed by the algorithm. Because the Euclidean DT calculation itself is a time-consuming process, the actual centerline computation time is considerably longer. Other approaches require even more computation time, such as 5 minutes (Silicon Graphics O2 workstation) [72] and 8 minutes (Silicon Graphics R10000 workstation) [69].
Finally, problems such as the potentially slow interpretation time and the variable diagnostic performance of CTC may also be addressed by use of computer-aided detection (CAD).
CAD could complement CTC by reducing the interpretation time, minimizing the perceptual error, and increasing the diagnostic performance of the interpretation of CTC [2]. Several CAD schemes have been developed for the detection of polyps in CTC during recent years. Most use shape-based features to locate suspicious regions from the colonic wall: Vining et al. [45] used curvature and wall thickness features of the colonic surface, Summers et al. [20,46] used curvature and sphericity features, Paik et al. [47,48] used the Hough transform, Yoshida et al. [33,49] used shape and texture features with a fully automated scheme, and Kiss et al. [50] used surface normals and sphere fitting. In addition, several methods have been developed for the reduction of false-positive (FP) detections. Göktürk et al. [51] developed a random orthogonal shape sections method for the differentiation of polyp candidates by use of support vector machines, Näppi et al.[52] developed several volumetric shape and texture methods for the differentiation of the shape and the texture of polyp candidates, Acar et al. [53] used edge displacement fields to model the changes in consecutive cross-sectional views of the CTC data, Jerebko et al.[54] used a committee of neural networks to classify polyp candidates, and Näppi et al.[35] developed a new method for the extraction and analysis of polyp candidates. The evaluation methods and the results vary among studies; generally, however, the results indicate that 70-100% of visible polyps ≧10 mm could be detected by CAD with an average FP rate of 1-8 detections per data set.
Colorectal masses are generally considered well visualized by a radiologist because of their size, continuation, and secondary findings. Therefore, little effort has been made to develop CAD schemes for mass detection. Nevertheless, the detection of both polyps and masses by CAD in CTC would provide a more comprehensive computer aid than the mere detection of polyps. If masses are not detected by CAD, the radiologist needs to perform a careful and complete review of all CTC cases for the presence of masses, which may increase the reading time. In addition, not all masses are easy to detect, depending on the readers' experience and on how rapidly the readers are reading the cases [55]. Therefore, the application of CAD for mass detection could improve the diagnostic accuracy of CTC by reducing potential reading errors due to reader fatigue, inexperience, or too rapid reading. Without explicit mass detection, the computer could also confuse masses with other types of lesions. Some types of masses have local surface bumps that tend to be detected as suspicious regions by a CAD scheme for the detection of polyps (
Previously, a multi-scale method for the detection of masses in CTC was developed [56]. Cap-like shapes were first detected by use of the shape index and the curvedness features in multiple scales based on a Gaussian pyramid of the CTC volume. The mass regions were extracted by use of a region-growing technique, and false positives were reduced by a quadratic classifier. In an evaluation with 20 patients, including 7 patients with a total of 8 colonoscopy-confirmed masses, all masses were detected by the method with an average FP rate of 0.45 detections per patient. Although the detection performance was satisfactory for that preliminary study, the method was well suited only for the detection of masses with a significant intraluminal component. Furthermore, the method added considerable computational overhead to the CAD scheme because it was implemented separately from the polyp detection scheme.
The present invention provides a new method for enhancing the endoscopic visualization of the colonic lumen. Each unique structure of importance in the colon can be highlighted automatically by a unique color, which makes it easy to perceive the relevant lesions at one glance. The present invention provides a method for coloring of the colon, which is based on a detailed theoretical framework of shape-scale analysis and can be adjusted systematically to delineate specific types of lesions. After the shape-scale specification, the different colonic structures are marked automatically with a unique color to provide an immediate visual differentiation of the lesions. In particular, polyps and diverticula, the visual differentiation of which may require additional time-consuming steps in monochromatically volume-rendered images (
The present invention also provides an algorithm for computation of the colon centerline that is not only robust to collapsed regions and poor segmentation, but is also faster than any previously reported method. The input to the algorithm is a distance map (DM) that is the result of a DT of the segmented colonic lumen. The colonic lumen is segmented by first thresholding of the air in the CTC data set, and then removing the lungs, bones, and skin. The regions representing lungs, bones, and skin are obtained by thresholding, where threshold values are determined adaptively from the histogram of the isotropic volume by identification of the characteristic peaks corresponding to air, fat, and muscle [77]. More precise segmentation algorithms, with respect to complete colon segments identification, exist [78]. However, in order to emphasize on the advantages of the present algorithm, the straightforward threshold-based segmentation approach was chosen.
The centerline algorithm is based on the well-known Kruskal algorithm [79] for generating a minimum spanning tree. In the present invention, the cost function is based on the DM value for the nodes and the Euclidean distance between nodes. Unlike the Kruskal algorithm, the iteration is allowed to stop before a single spanning tree was formed.
The centerline algorithm first extracts local maxima in the DM, divided into two categories—seed points and regular points. Seed points are the subset of local maxima that has DM values larger than a pre-defined threshold level. Regular points are the remaining local maxima. Starting with maxima with the largest distance values, available maxima are iteratively transferred to be single-node graphs in a graph structure. During each iteration, graphs are joined by linking of all pairs of nodes that satisfy a set of connection criteria. After the last iteration, branches are removed from the remaining graphs. The end segments of the remaining graphs are then recovered, after which the actual centerline is represented as the set of graphs containing at least one seed point each.
Further, the present invention provides a CAD method for the detection of intraluminal and non-intraluminal types of colorectal masses. Intraluminal masses are masses with a significant intraluminal component, such as lobulated, polypoid, or circumferential types of masses. Non-intraluminal masses, are masses that are associated with a mucosal wall-thickening type of growth pattern rather than those with a significant intraluminal component. Two mass detection methods were developed, one for each type of mass, because these two types can appear in very different forms in CTC (
According to one aspect of the present invention, there is provided a method, system and computer program product for displaying volumetric image data of a target organ as a visual aid in detecting at least one abnormality in the target organ, comprising (1) obtaining the volumetric image data of the target organ, the volumetric image data including a plurality of voxels having a respective plurality of image data values; (2) determining at least two geometric feature values at each voxel in the plurality of voxels based on the plurality of image data values; (3) determining a colormap by assigning a color to each possible combination of the at least two geometric feature values, wherein each assigned color corresponds to a different region type within the target organ; (4) determining a display color for each voxel in the plurality of voxels based on the determined colormap; and (5) displaying, based on the determined display color for each voxel in the plurality of voxels, a color image representing the volumetric image data of the target organ.
Further, according to one aspect of the present invention, there is provided a system for displaying volumetric image data of a target organ as a visual aid in detecting at least one abnormality in the target organ, comprising: an image acquisition unit configured to obtain the volumetric image data of the target organ, the volumetric image data including a plurality of voxels having a respective plurality of image data values; a processor configured (1) to determine at least two geometric feature values at each voxel in the plurality of voxels based on the plurality of image data values, (2) to determine a colormap by assigning a color to each possible combination of the at least two geometric feature values, wherein each assigned color corresponds to a different region type within the target organ, and (3) to determine a display color for each voxel in the plurality of voxels based on the determined colormap; and a display unit configured to display a color image representing the volumetric image data of the target organ, based on the determined display color for each voxel in the plurality of voxels.
Further, according to another aspect of the present invention, the step of determining the colormap comprises obtaining simulated or clinical volume data representing the target organ; determining, based on the simulated volume data, the at least two geometric feature values for each voxel in the simulated volume data; selecting a region type of the target organ; determining combinations of the at least two geometric feature values that most frequently correspond to the selected region type; and assigning a color to each of the determined combinations of the at least two geometric feature values that correspond to the selected region type.
Further, according to another aspect of the present invention, the step of determining the at least two geometric features comprises: determining, based on the plurality of image data values, a shape index at each voxel, the shape index being indicative of a local shape at each respective voxel; and determining, based on the plurality of image data values, a curvedness index at each voxel, the curvedness index being indicative of a magnitude of the local shape at each respective voxel.
Further, the present invention provides a method of processing a set of volumetric data encompassing a target organ, comprising detecting multiple lobulated regions of the targeted lesion; and reclassifying the detected regions to detect the entire lesion. In addition, the present invention provides that the detecting step comprises computing volumetric features characteristic of lobulated regions at each point in the volume representing the targeted lesion; and determining the possible lobulated regions of the targeted lesion based on the computed volumetric features. Further, the reclassifying step comprises merging the detected lobulated regions by fuzzy merging; and expanding the merged lobulated regions to extract the entire lesion. Moreover, the merging step comprises determining the membership grade for each of the detected lobulated regions; and clustering the lobulated regions that have a high membership grade. The expanding step comprises applying a volumetric region expansion method to the merged lobulated regions, and the step for determining the membership grade comprises formulating a membership function for any two lobulated regions, which has a high value when the two regions belong to a single lesion, based on the CT and gradient values along a ray that connects the two regions; and calculating a membership grade for each of the lobulated regions based on the membership function. Finally, the step for applying a volumetric region expansion comprises using a modified level set method with fast marching algorithm.
A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:
indicates that A and B are located within the same lesion;
b show a perspective rendering of the colonic lumen with computer-generated centerline (33A); and a perspective rendering of the displacement at each location in 3-D between the human-and computer generated centerlines (33B);
The description of the present invention includes a written description of: (1) visualization of the colon by shape-scale signatures; (2) computation of the colon centerline in CT colonography; and (3) detection of colorectal masses in CT colonography based on fuzzy merging and wall-thickening analysis.
I. Visualization of the Colon by Shape-Scale Signatures
Local Shape Index and Curvedness
The local shape index and curvedness measure the shape characteristics of a local surface patch at a point, and they form the basis for defining the shape-scale spectrum which is described in the next section. Both quantities are defined based on the notion of the principal curvature. A traditional approach for computing the principal curvatures of 3-D data is to fit a parametric surface to the data and compute its differential characteristics in a local coordinate system [24]. However, parameterization of surfaces with a complex topology is complicated, and curvature information is often needed before these surfaces are generated. Therefore, the differential characteristics of the iso-surfaces were computed directly from the 3-D data volume, without extraction of any surfaces, based on an implicit representation of the iso-surface [25,26]. That is, the differential characteristics are computed locally at each voxel by use of the differentials of the volume in three directions. The interpretation is that the computed differential characteristics at voxel p correspond to the differential characteristics of the iso-surface passing through p. The definition and the computational method of the local shape index and the curvedness are described briefly in the following.
Let h(p) denote a voxel value (typically, the CT value) at a voxel p=(x,y,z)∈IR3. A 3-D iso-surface P at a level a is given by
P={p=(x,y,z)∈IR3;h(p)=a}.
To obtain an infinitely differentiable 3-D function, the volume is convolved with a Gaussian filter. The implicit-function theorem indicates that, at each voxel p, there exists a small neighborhood U of p in which z can be expressed by a smooth function φ of x and y. By denoting (x, y) with (u, v) in this neighborhood, the iso-surface P can be represented in a parametric 2-D form as
P={(u,v)∈IR2;h(u,v,φ(u,v))=a}.
The principal curvatures are obtained as the eigenvalues of the following matrix [27], called the Weingarten endomorphism:
Here, E, F, and G are the first fundamental forms [27] defined at a voxel p as
E≡Pu·Pu, F≡Pu·Pv, G≡Pv·Pv, (2)
where Pu and Pv are the partial derivatives of P in terms of u and v at p, respectively. In the following, a suffix of u, v, x, y, or z indicates the partial derivative in terms of these variables, and all of the quantities are defined on a single voxel p unless mentioned otherwise. L, M, and N are the second fundamental forms [27] defined as
L≡Puu·Q, M≡Puv·Q, N≡Pvv·Q, (3)
where Q is a normal vector to the iso-surface P defined by Q=Pu{circumflex over ( )}Pv.
To calculate the principal curvatures, the implicit differentiation and the chain rule is used to obtain Pu=∂P/∂u=(1,0,∂φ/∂u)=(1,0,hx/hz). Substituting similar expressions for Pv, Puu, Puv, and Pvv in Eqs. (3)-(7), one can compute the determinant of the matrix W, called the Gaussian curvature, K, and half of the trace of W, called the mean curvature, H, as follows:
Here, Perm(x, y, z) denotes a permutation of (x, y, z), i.e., Perm(x,y,z)≡{(x,y,z),(y,z,x),(z,x,y)}. Because the principal curvatures κ1 and κ2 are the eigenvalues of the Weingarten endomorphism of Eq. (3), the Gaussian curvature and the mean curvature can also be defined by [28]
From Eq. (4), it is easy to derive that the principal curvatures at the voxel p can be expressed as
κi(p)=H(p)±{square root}{square root over (H2(p)−K(p))}, i=1,2.
Let us denote κmin=min{κi,κ2} and κmax=max{κ1,κ2}. Then the local shape index S(p) and the local curvedness C(p) at the voxel p are defined as [29]
The local shape index measures the shape of a local iso-surface patch at a voxel. Every distinct shape corresponds to a unique value of the shape index with a value range of [−1,1] (
The idea of determining the local shape and its “scale” at each voxel is somewhat similar to multi-scale methods, where an input signal (e.g., image data) is typically divided into multiple new signals, each representing the original signal at a unique scale [30,31,32]. The original signal may then be studied and manipulated independently at various scales. However, it is sometimes challenging to establish a useful relationship between the various scale representations of the same signal. In contrast, the local shape index and the curvedness capture the intuitive notion of the local shape of an iso-surface at a voxel, and provide an efficient representation of the colonic surface for our application.
The Shape-Scale Spectrum
The local shape index and curvedness value pair (S(p), C(p)) of each voxel p of a 3-D data set, defined by Eqs. (5) and (6), can be used for constructing a mapping from the voxels of the data set to a Cartesian plane called the shape-scale spectrum (
The differentiation of colonic polyps, diverticula, folds, and the colonic wall was examined. Each of these lesions has some unique characteristics. Polyps are generally bulbous, locally cap-like structures that adhere to the colonic wall or to a fold. Diverticula are small, cup-like structures within the colonic wall. Folds are elongated, ridge-like structures that arise from the colonic wall. The underlying colonic wall can be characterized as a rut-like spinode. The shape index values of these lesions are ideally S≈1 for polyps, S≈−1.0 for diverticula, S≈0.5 for folds, and S≈−0.5 for the colonic wall. The curvature values are expected to be high for polyps and diverticula, relatively low for the colonic wall, and moderately high for folds. The highest values of curvedness do not occur in CTC data because of the lack of sharp biological objects and the smoothing effects of the CTC image formation process.
Generation of Shape-Scale Signatures
There are two methods for generating shape-scale signatures of the colonic lesions of interest. The first method derives the signatures from computer-detected lesions in real clinical data, and the second method derives the signatures from computer-simulated lesions that imitate the shape of the lesions of interest. Real clinical data are expected to produce more realistic shape-scale signatures than do those generated by computer-simulated lesions; however, obtaining relevant clinical data is complicated. On the other hand, the development and evaluation of realistic computer-simulated phantoms can be challenging, but their shape-scale signatures can be obtained accurately because the regions of the simulated lesions are known precisely. If the computer-simulated lesions are sufficiently realistic, they may also be used for testing the effect of a color map on the coloring of colonic structures, and several different examples of a lesion can be simulated with ease by appropriate variation of the simulation parameters.
Use of Clinical Data
The clinical data used in this study were obtained from the CTC examinations performed during 1997-2001 at the University of Chicago Hospitals. Each patient underwent standard pre-colonoscopy cleansing, and the colon was insufflated with room air. The CTC scanning was performed with a helical CT scanner (GE 9800 CTi or LightSpeed QX/i; GE Medical Systems, Milwaukee, Wis.) in the supine and prone positions. The collimation was 2.5-5.0 mm, the pitch was 1.5-1.7, the reconstruction interval was 1.5-2.5 mm, the current was 60-100 mA with 120 kVp, and the matrix size of the axial images was 512×512, with a spatial resolution of 0.5-0.9 mm per pixel. Optical colonoscopy was performed on the same day as the CTC, and radiologists confirmed the presence and location of polyps and diverticula in the CTC data sets by use of colonoscopy reports, pathology reports, and a GE Advantage Windows workstation (GE Medical Systems, Milwaukee, Wis.). Fourteen patients had a total of 21 colonoscopy-confirmed polyps; 20 polyps were 5-12 mm, and one was 25 mm in diameter.
To determine the shape-scale signatures of actual colonic lesions, four CTC data sets which contained pedunculated polyps, sessile polyps, and dozens of diverticula were selected. The polyps were detected by use of our previously developed fully automated CAD scheme for the detection of polyps from CTC data sets, which is described briefly below. Diverticula were detected by use of a new method. Folds were sampled from among false-positive detections as described below.
The CAD scheme for the detection of polyps has been described in detail elsewhere [33,34,35]. Briefly, each input data set is interpolated linearly in the axial direction so that isotropic volumes are obtained, and the region of the colonic wall is extracted by use of a knowledge-guided segmentation method [36]. The polyp-like regions are detected based on the shape index and curvature features with hysteresis thresholding, and extracted by use of conditional morphological dilation [35]. The latter extraction step ensures that the shape-scale signatures of the detected polyp candidates are those of their visually perceived regions rather than those within the range of the predetermined shape index and curvature values for initial detection.
Because the CAD scheme described above was designed to detect only polyps, it was necessary to develop a separate method for the detection of diverticula. Because the cup-like diverticula are opposite in shape to the cap-like polyps, the diverticula can be detected similarly to polyps, but by mirroring the values of the shape index used to search for suspicious regions. That is, whereas polyps are detected by searching for shape index values S≈1, diverticula can be detected by searching for shape index values S≈−1 within the colonic wall. The curvedness values of diverticula are similar to those of polyps, although diverticula are expected to have a narrower range of curvedness values than do polyps. To extract the region of a detected single diverticulum, the detected region is expanded by morphological dilation with a spherical kernel. The size of the dilation kernel depends on the size of the detected region: large regions are dilated with a larger kernel than that used for small regions.
To reduce false-positive detections of diverticula, all air voxels were thresholds within the colon (CT value<−700 HU). Then, a morphological opening of 2.5 mm was applied to the thresholded region. Only the voxels within the detected diverticula that were not part of the air region were retained. Finally, the remaining detections that were smaller than 5 mm or larger than 10 mm were eliminated.
It should be noted that a detailed quantitative analysis of the performance of this detection method was not performed; instead, the detected regions were inspected visually. The visual assessment indicates that the method detected all of the visible diverticula, although it detected some false positives due to cup-like regions appearing at bases between folds which connect smoothly to the underlying colonic wall. These regions remained false positives because of the relative simplicity of the false-positive reduction method described above. Therefore, a more elaborate method would be required to achieve a clinically acceptable performance, which is left for a future study. Nevertheless, highlighting both true- and false-positive diverticula is important for the purposes of the visualization scheme presented in this study, because both types of detections are very different in shape from polyps, and thus they are useful for indicating regions that are unlikely polyps. For generating shape-scale signatures, it sufficed to select manually 10 examples of computer-detected diverticula.
Use of Computer-Simulated Lesions
The computer-simulated lesions were generated in a 384×384×512 volume. The isotropic voxel resolution of the volume was assumed to be 0.7 mm, which is the average image resolution in our clinical CTC database. The approximate range of the shapes and sizes of the simulated polyps, diverticula, and folds were determined based on actual lesions in the 14 cases containing polyps in the CTC database.
The colonic lumen was simulated by a cylinder 1.5-7.0 cm in diameter. The surface of the cylinder represented the surface of the colonic wall and was assigned the default background color. Folds were simulated by placement of several overlapping solid spheres or ellipsoids along the boundary of the simulated colonic wall in the axial plane. Seven thick spherical folds and 7 thin ellipsoid folds with diameters of 3.5-21.0 mm were simulated. Sessile polyps were simulated by placement of the center of a solid sphere or ellipsoid at the boundary of the colonic wall, and pedunculated polyps were simulated by placement of the whole sphere or ellipsoid on the boundary of the colonic wall. In each case, 7 polyps with diameters of 2.8-21.0 mm were simulated, resulting in a total of 28 simulated polyps. A diverticulum was simulated by placing of a zero-valued sphere or ellipsoid at the boundary of the colonic wall to generate a “hole” that represents the diverticulum. Seven spherical and 7 ellipsoid diverticula with diameters of 2.8-7.0 mm were simulated. To simulate the blurring of the clinical data acquired with a CT scanner, we applied a Gaussian smoothing filter with σ=2.0 to the volume. The smoothing also reduces the artificially sharp edges that could appear within the simulated objects.
The Shape-Scale Color Map: Characteristic Signature and Color Interpolation
In practice, for obtaining 2-D color maps, the original shape-scale signatures may not be appropriate for the coloring of the shape-scale spectrum. This is because the shapes and scales of the different types of lesions may share common elements, and thus the shape-scale signatures of the different types of lesions may overlap.
To determine representative, but non-overlapping shape-scale regions for each type of lesion for use in the generation of color maps, the characteristic signature Si of the shape-scale signature for each lesion type i is determined. A characteristic signature is a subregion of the shape-scale signature that contains the most frequently occurring combinations of the local shape index and the curvedness values for the type of lesion that the shape-scale signature represents. Once the characteristic signature has been determined, it is assigned the unique color li that was chosen to represent lesions of type i. The points on the shape-scale spectrum that are not colored in this manner are assigned a default background color. On the other hand, the points in the transition region between a characteristic signature and the background region are colored by use of color interpolation, as described below.
Let (I1,I2,I3) represent the components of the color assigned to the characteristic signature Si. Although these could be the red, green, and blue components of the RGB color space, the use of the RGB space in color interpolation may result in false colors. For a precise result, the color interpolation was performed in the YCrCb (YUV) [37,17] rather than the RGB color space.
Let Di(c,s) denote the shortest Euclidean distance between the boundary of the characteristic signature Si and a point (c,s) outside Si in the shape-scale spectrum. Then, each color component thcalCj(c,s) (j=1,2,3) of the color at (c,s) is calculated as
where R is the number of the different types of lesions to be differentiated by the coloring scheme, Bj is the magnitude of the color component Cj in the color assigned to the background, kn is the sharpness parameter for the color of the lesion of type i, and T is a subsampling factor. The background is considered as an undefined lesion not in i=1, . . . ,R. The sharpness parameter controls the boundary profile of the characteristic signature Si: high values of kn produce sharp boundaries, and low values produce smooth boundaries.
Subsampling of the shape-scale signatures occurs when T>1 (T is set to 1 if no subsampling is desired), and can be used to limit the size of the color lookup table for faster processing. For example, if the shape index and curvedness values are represented by 1,000 entries, the resulting 2-D color lookup table would have 1,000,000 entries. For a setting of T=6, only 27,778 entries would be needed.
The color interpolation scheme provides an effective parameterization of the color map construction, and it can also be used for controlling the overlap between the shape-scale signatures. The important parameters are the sharpness parameters of the characteristic signatures and the set of the unique colors assigned to the characteristic signatures. The effects of these parameters are described below.
To implement the color map scheme, the occurrence frequencies of the local shape index and curvedness values within the shape-scale spectrum were first scaled linearly to [0,255] for each lesion type i. Then, to obtain the characteristic signature, each shape-scale signature was thresholded at the level of 84 (33% of maximum value). To limit the size of the color lookup table, the 2-D color map was subsampled by setting of the subsampling factor to T=6 in Eq. (7).
Shape-Scale Color Map: Coloring of Colon
In the third method, called the hybrid shape-scale coloring method, the shape-scale color map information is combined with the output of our CAD scheme. Although the coloring of the colon by the first two methods enhances the lesions of interest, it is based only on the shape or scale of the lesions. Thus, false positives, such as stool that have a polyp-like shape, may be color-enhanced as well. However, our CAD scheme for polyp detection includes an effective method for the elimination of false-positive detections by an analysis of, for example, the internal texture of the polyp candidates [35]. As a result, the scheme generates only a few, if any, false-positive polyp detections for a given data set. Thus, in the third shape-scale coloring method, we first color the colon by use of either of the first two coloring methods. Next, the regions extracted by the CAD scheme, such as polyps or diverticula, are colored with pre-specified colors. The resulting coloring of the colon is expected to delineate the lesions of interest more precisely than does the application of the first two coloring methods.
In step 3401, a colormap is determined by assigning a color to each possible combination of the at least two geometric feature values. The geometric feature values are, e.g., the shape index and the curvedness index, which can be computed at each voxel in a three-dimensional data set. As described above, each assigned color corresponds to a different region type (e.g., polyps) within the target organ.
In step 3402, the volumetric image data of the target organ is obtained. The three-dimensional volumetric image data includes a plurality of voxels having a respective plurality of image data values.
In step 3403, at least two geometric feature values at each voxel in the plurality of voxels are determined based on the plurality of image data values. For example, the shape index and the curvedness index are calculated at each voxel defining the volumetric image data.
In step 3404, a display color for each voxel in the plurality of voxels is determined based on the determined colormap.
Finally, in step 3405, a color image representing the volumetric image data of the target organ is displayed based on the determined display color for each voxel in the plurality of voxels. As described above, a radiologist can use the color image to identify abnormalities in the target organ. Moreover, a CAD method may be used to identify abnormalities in the target organ, which are then displayed in a predetermined color.
The image segmentation unit 3530 is configured to extract a set of voxels from the total scanned volume, forming a database of voxels representing the segmented target organ, which is stored in the segmented organ database 3502. Abnormalities may be detected using the segmented image data of the target organ by the abnormality detection unit 3540 based on features of the image data determined by the feature calculation unit 3550. Data regarding the abnormalities may be stored in the detected abnormalities database 3503.
The feature calculation unit 3550 also calculates the shape index and the curvedness index at each voxel in the volumetric image data. The colormap determining unit 3520 assigns a color to each combination of the at least two geometric feature values (e.g., shape index and curvedness) calculated by the feature calculation unit 3550, as described above. The color assigning unit 3560 uses the colormap determined by the colormap determining unit 3520 to assign a color to each voxel in the volumetric image data. The display unit 3570 is configured to display a color image of the volumetric image data according to the colors assigned by the color assigning unit 3560.
CAD Workstation for Polyp Detection
To examine the effect of the shape-scale color maps in the coloring of CTC data sets, the color map scheme was implemented in a CAD workstation for the detection of polyps in CTC (
The workstation displays the reconstructed volumetric data sets with two types of display modes: 2-D multiplanar reformatted (MPR) views including three views perpendicular to the axial, coronal, and sagittal planes of the volumetric data, and a 3-D endoscopic view that displays the colonic lumen with a perspective view. The users can examine the data set by scrolling through the images in the MPR views, by zooming into regions of interest, by manipulating the grey-scale display range, and by adjusting the 3-D view at will. The output of the detection of polyps is provided as a list of suspicious locations in a window next to the MPR views. By clicking on the list, the user can view the location of a suspected polyp on the MPR and the 3-D views.
The workstation provides these display modes with a grey-scale mode as well as with a color mode that is based on a shape-scale color map scheme. When the latter is selected, all of the visible region in the 3-D endoscopic view is colored based on the shape-scale color map, whereas in the MPR views, only the region corresponding to the colonic wall is colored with the color map.
Evaluation Method
The effect of the shape-scale coloring scheme was evaluated in two steps. In the first step, the computer-simulated lesions and 14 clinical data sets containing polyps were colored with the application of the shape-scale color maps. Several versions of the color maps were generated by adjustment of the color map parameters (i.e., the sharpness parameters and the color assignment). The effect of the color map parameters and the enhancement of the specified colonic lesions were examined visually with the CAD workstation.
In the second step of the evaluation, an observer study was performed to examine the effect of the hybrid shape-scale coloring method in the detection of colonic polyps in CTC. In this study, 20 CTC cases collected as described above were used. For each patient, either the supine or the prone data set was used. There were 10 abnormal data sets containing 11 polyps of 5-12 mm, and 10 normal data sets with no polyps. The detection performance of the computerized detection scheme was set to a by-polyp sensitivity of 92% with 1.0 false-positive detections per data set on average. The shape-scale color map was constructed by use of a high value of the sharpness parameter (kn=0.5). Folds were assigned a pink color, and the background was assigned a dark red color. The regions detected by the CAD scheme were colored green.
Three radiologists examined the data sets by use of the CAD workstation. Observer 1 was a radiologist experienced in the interpretation of CTC data sets, observer 2 was a radiologist not experienced with CTC, and observer 3 was a gastroenterologist. A sequential test was used, in which an observer first reads a data set without the use of the color map, and rates the confidence level regarding the presence of at least one polyp ≧5 mm on a 0-100% scale. Next, the colon is colored by the hybrid shape-scale color map, and the observer rates the confidence level again. After the observers had examined all data sets in this manner, receiver operating characteristic (ROC) analysis [38] was used to compare the observers' performance in the detection of polyps without and with the coloring scheme. For each observer, a binormal ROC curve was fitted to the confidence rating data of each observer by the use of maximum likelihood estimation. A PROPROC program was used for obtaining a maximum likelihood estimate of “proper” binormal ROC curves from the continuous ordinal-scale rating data. The difference in the observers' performance between the first and the second confidence ratings was characterized by use of Az, the area under the best-fit binormal ROC curve, when the curve is plotted on a unit square. Generally, the larger the Az value (0.0≦Az≦1.0), the higher the performance. The statistical difference in the observer performance between the monochromatic and hybrid shape-scale colorings was evaluated by use of the two-tailed t-test to the pair of the Az values.
Visual Effect of Shape-Scale Color Maps
Observer Study
The results of the observer performance in the detection of polyps are shown by ROC curves in
The results indicate that the detection performance of all observers increased when the colon was colored by the hybrid color map. The two-tailed t-test for paired data showed that the difference of the Az values between the monochromatic and hybrid shape-scale coloring was statistically significant (p=0.027), which confirms the usefulness of this color map. The increment was smallest for the expert radiologist (observer 1). However, the performance of the two other, inexperienced observers increased substantially and reached approximately the performance level of the expert.
The visual evaluation of the application of shape-scale color maps indicates that, as expected, colonic structures are enhanced in the colored colon. The overall perception and visibility of the lesions can be adjusted by the choice of colors. For example, the preciseness of the coloring can be adjusted by the sharpness parameter kn. A low value of kn spreads the color of lesion type i over a wider range of shapes and sizes than does a high value of kn that tends to color only the most typical instances of the target lesion. The shape-scale signatures from clinical data and computer-simulated data were rather similar, and the use of a low sharpness value will further increase their similarity.
The hybrid shape-scale color map construction yielded the most precise delineation of the target lesions. However, errors such as false-positive detections or incorrect segmentations of the detected lesions will show up clearly in the coloring of the colon by this approach, and false-negative detections might be colored as background if the choice of the colors is appropriate. Incorporation of the CAD data also makes this method slower than the two other methods. The observer study was also limited in the sense that only the presence of polyps in the colon was considered. The hybrid color map method appears to be the most effective approach for this purpose because the coloring of the polyps is clear and simple to interpret. On the other hand, the other color map methods may be useful for providing more versatile and rich representations of the information on colonic lumen, thereby potentially enhancing the analysis of the colon. However, an extended training period for radiologists may be required to efficiently interpret such visualizations of the colon. The analysis of the clinical benefits of the more complicated color map representations is a topic for a future study.
The results of the observer performance study must be considered preliminary because of the small number of observers (3) and data sets (20). Nevertheless, the results are similar to those obtained previously in the evaluation of CAD schemes for breast and chest imaging: the performance of a radiologist is increased when a computer aid is used, and the performance of non-experienced readers may be improved to the level of expert readers [39,40]. Therefore, one can expect similar performance results with this method when a larger database and a larger number of observers are employed. The use of the scheme has also other potential benefits that were not addressed directly by the study. If the color map visualization is available as soon as a CTC case has been loaded, the radiologist may be able to find the polyps faster than with the use of the standard screening schemes. Also, the color map visualization may boost radiologist's confidence in the diagnosis, thereby allowing the radiologist to screen CTC cases faster than with standard methods in which the visual representation of the colon can be more ambiguous.
Compared with previous work in which fully automated methods for the detection of polyps in CTC were developed, the present invention is unique because it focuses on the visualization of the colon for radiologists. Various visualization schemes are presented that can be used to color lesions of interest in the colonic lumen based on a systematic method. The required color maps may be prepared by use of either simulated or actual lesions, and the output of the CAD scheme can be incorporated into the visualization scheme by use of the hybrid coloring method. Unlike in previous studies, an observer study was provided to evaluate the effect of the coloring schemes in polyp detection by radiologists.
Flat lesions were not included in the study because the clinical database did not contain such lesions. However, although flat lesions hardly protrude from the colonic wall, they are still expected to have the local cap-like shape of a polyp, and the curvedness values of flat lesions are expected to be lower than those of typical polyps. Because the region that corresponds to flat lesions in the shape-scale spectrum is not expected to have a significant overlap with the other colonic non-polypoid lesions, the visualization of flat lesions may be implemented simply by an adjustment of the curvedness range of the polyp signatures.
A new visualization method for a rapid examination of the colon by use of virtual endoscopy was developed. The colonic structures of interest are detected and colored based on their characteristic signature patterns in the shape-scale spectrum. The method was evaluated by use of computer-simulated lesions, clinical CTC data sets, and an observer study. The results of the visual examination indicate that shape-scale color maps can be used for enhancing the visibility and conspicuity of the target lesions. The results of the observer study indicate that the performance of readers not experienced in the detection of polyps in CTC can be improved substantially by use of the hybrid shape-scale color map scheme. Combined with our CAD workstation for the detection of colonic polyps, the visualization scheme is expected to reduce radiologists' interpretation time and to improve the diagnostic performance of the detection of colonic neoplasms in CT colonography.
The present method of shape-scale signatures can also be implemented more generally by one of ordinary skill in the art for the detection of other abnormalities of other organs. In particular, the present method is applicable to any type of convex (peak-like) or concave (pit-like) abnormalities in N dimensions (N>1). Thus, an embodiment of the present method can be readily applied to 2D/3D aneurysms, embolisms, solid lung cancer (especially those that stick on the chest wall), solid stomach cancer, etc.
II. Computation of Colon Centerline in CT Colonography
According to the present invention, there is provided a method for the computation of a colon centerline that is faster than any of the centerline algorithms previously presented and that relies less than other algorithms on complete colon-segments identification. The method first extracts local maxima in a distance map of a segmented colonic lumen. The maxima are considered to be nodes in a set of graphs, and are iteratively linked together, based on a set of connection criteria, giving a minimum distance spanning tree. The connection criteria are computed from the distance from object boundary, the Euclidean distance between nodes and the voxel values on the pathway between pairs of nodes. After the last iteration, redundant branches are removed and end segments are recovered for each remaining graph. A subset of the initial maxima is used for distinguishing between the colon and non-colonic centerline segments among the set of graphs, giving the final centerline representation.
Detection of Local Maxima in DM
Given a DM of the segmented colonic lumen (air-filled region), a subset that consists of the local 6- and 4-neighborhood maxima within the DM is extracted. The size of a maximum is defined to be its DM voxel value, and the position is defined to be its voxel coordinates. Here, a 6-neighbor maximum is defined as a voxel vx,y,z within the DM at which all six neighborhood voxels (vx±1,y,z, vx,y±1,z, vx,y,z±1) have a value smaller than that of the voxel. It can be interpreted as the center of the locally largest sphere fully enclosed by the object. A 4-neighbor maximum is a voxel vx,y,z.in the DM at which all neighboring voxels in two of the three 6-neighbor directions have a value smaller than that of the voxel (i.e., vx±1,y,z, vx,y±1,z, or vx±1,y,z, vx,y,z±1 or vx,y±1,z, vx,y,z±1). A 4-neighbor maximum can be interpreted as the center of the locally largest cylindrical structure that is fully enclosed by the object.
A point at the center of a segment of the colon usually has a fully enclosed, locally largest, sphere due to the folds within the colon, or a locally largest fully enclosed cylindrical structure due to the cylindrical shape of the colon. Therefore, a good representation of a centerline can be found by extraction of the 6- and 4-neighborhood maxima. Our algorithm uses a subset of these maxima as input. The subset, M, is chosen so that, for every maximum Mi in M, there is no other maximum in M that has a distance to Mi less than the size of Mi scaled by a constant (i.e., the seedscalefactor or regularscalefactor, see Section III).
Seedpoint Extraction
Among the set M of maxima, we select an initial set, MS, of maxima in M consisting of 6-neighbor maxima with sizes larger than a user-defined minimum-size threshold value (i.e., minseedsize, see Section III). MS will serve as a set of seed points for the location of the centerline. The seed points are used as a guide for determining which of the final centerline segments actually belongs to the colon. The lower the minimum size threshold value, the more sensitive to non-colonic objects the algorithm is. The higher the minimum size threshold value, the less likely it is that the full colon centerline is found. A parameter that defines the minimum number of seed points selected (i.e., minseedamount, see Section III) ensures that the algorithm will work even in cases where the colon is extraordinarily narrow, even if this implies that the minimum-size threshold value is violated.
Let G represent a set of graphs. Each maximum in MS is considered to be a node in G, assigned the position of the maximum (i.e., the voxel coordinates in 3-D space) and the size of the maximum (i.e., the DM voxel value of the maximum), and removed from M. Initially, each node in G is equivalent to a graph (see
Graph Connection. Let GA and GB be graphs in G. A link between a node A in GA and a node B in GB is created, joining the two graphs, if the following three conditions are satisfied:
The above criteria 1-3 are applied to all possible nodes in G. The resulting G consists of one or more graphs, each graph consisting of one or more nodes, and representing a minimum distance spanning tree within the colonic lumen (see
Window-Based Iteration
Starting with the remaining maxima in M with the largest size and iterating down to the maxima with the smallest size, maxima are removed from M and instead are considered as new single-node graphs within G (see
Branch Cutting
Each graph in G contains a major centerline and several additional branches, but no loops. Iteratively, for graphs GN within G containing a node with more than two links, all available single-ended nodes in GN (i.e., the end nodes of each branch in GN) are removed. The iteration stops when no graph in G contains a node with more than two links (see
Branch Recovery
The branch-cutting process removes not only unwanted branches, but also the end segments of the actual centerline. These unintentionally removed segments are recovered for each end node in each graph in G. For an existing end node in G, end-node candidates are defined as the nodes in all of the branches that were once connected to the existing end node, but were later removed during Branch Cutting. One of the end-node candidates is selected as the new end node. The branch, or part of the branch, that connects the existing and the new end node is recovered. The selection of the new end node is performed in the following manner: First, we measure the distance from each of the end-node candidates to the existing end node. Second, we measure the distance from each of the end-node candidates to the node with a link to the existing end node. Third, comparing the two distances measured for each end node candidate, nodes that have the longer distance to the existing end node are removed from the list of endnode candidates. This prevents a branch that turns back to the extracted centerline from erroneously being recovered as an end segment. Finally, among the remaining end-node candidates, the new end node is chosen as the end-node candidate with the largest distance to the current end node (see
Save Result
Each graph containing at least one seed point is stored as a resulting centerline segment. All other graphs are removed because a graph without a seed point is considered to be a centerline segment for an extra-colonic structure.
For evaluating the effect of different parameters, we validated our algorithm on a computer-generated colon phantom. The input is an arbitrary curve in 3-D space and the radius, R, of the colon phantom to be generated. The phantom is generated as follows: Traversing the curve given as input, at a uniform path distance of 1 voxel, compute the radius of a cylinder with thickness of 1 voxel, and position it with the main axis in the path direction and centered on the curve. The radius of the cylinder is computed as R with additional random noise with range of 2 voxels. The exception is that, at a uniform distance of 40 voxels of path length, a fold is created. A fold is created by linearly decreasing the radius to 20% of R, after which the radius is increased linearly to its initial value. The path length of the fold-decreasing and -increasing phases is 10 voxels each. At this stage, the phantom contains has the overall shape and randomness, but folds are unnaturally deep and saw-tooth shaped, and the surface is jaggy and synthetic looking. By dilation with a spherical dilation kernel with a radius of 6 voxels, the final 3-D phantom of the colonic lumen was created, with the intended look and fold sizes.
In the following, the terms “reference” and “slave” centerlines are used. The reference centerline is considered the “true” centerline. For quantifying the degree of deviation of a slave centerline from a reference centerline, primarily two types of quantities were calculated.
Displacement for a single point on the reference centerline is defined as the shortest Euclidean distance to the slave centerline with a straight path within the lumen, if existing. Otherwise, displacement is not defined for the point. For a reference centerline, mean displacement is defined as the average of all displacements along the centerline, max displacement is defined as the largest of all displacements along the centerline, and standard deviation is defined as the standard deviation of all displacements along the centerline. The centerline extension is the fraction (computed as percentages) of the number of defined displacement points when compared to the number of both defined and undefined displacement points.
For a set of centerlines, the average extension, the average of the maximum, the average of the mean, and the average standard deviations of the slave centerline when compared to the corresponding reference centerlines, are the averages over all cases of the centerline extensions, max displacements, mean displacements, and standard deviation, respectively (see
The centerlines generated by the algorithm presented, were always regarded as being the slave centerlines in the phantom studies.
For evaluating the accuracy of the centerline generated by our algorithm on clinical cases, a total of 20 CTC cases were selected from CTC examinations performed at our institution between 1997 and 2000. All patients underwent CTC after a precolonoscopy preparation of barium enema bowel cleansing, by use of 48-hour liquid diets and agents such as polyethylene glycol, magnesium citrate, or bisacodyl and phospho soda. Then a soft-tip tube without a training cuff was inserted into the rectum, and the colon was gradually insufflated with room air until the patient experienced mild discomfort (approximately 40 puffs or 1 liter). CT scanning was performed by use of a helical CT scanner (GE9800 CTi and LightSpeed QX/i; GE Medical Systems, Milwaukee, Wis.) with 2.5 mm to 5 mm collimation, and reconstruction intervals of 1.5 mm to 2.5 mm. The matrix size of the resulting axial images was 512×512, with a spatial resolution of 0.5 mm/pixel to 0.7 mm/pixel. A reduced current of 60 mA or 100 mA with 120 kVp was used for minimizing the radiation exposure. CT scanning was performed in both supine and prone positions, and thus 20 cases included a total of 40 CTC scans (i.e., 20 supine and 20 prone CTC scans).
Although the colon was cleansed and distended for all of the patients, the optimality of the cleansing and distension varied among patients, and thus some of the colons had residual stool and fluid as well as non-distended segments of the colon, as often reported by other CTC studies [20]. To evaluate the robustness of the centerline algorithm under a wide spectrum of conditions, we included cases with feces, fluid, a poorly distended colon, and motion artifacts in the 20 cases as described below: (1) 20 CTC scans, classified as “straightforward,” had a well distended colon, with only a small number of narrow sections and thin colonic walls; (2) 10 CTC scans, classified as “moderate,” contained narrow and collapsed regions and thin walls; and (3) 10 CTC scans, classified as “challenging,” had a poorly distended colon with thin walls, and a large number of collapsed regions due to residual stool, fluid, and spasm. Also, motion artifacts were present in some of the cases.
Generally, the segmented cases had a large number of fragments of small bowel, lungs, and stomach remaining in the segmented colon. No case was completely free of extra-colonic components (see
Three radiologists, A, B, and C, manually drew, independent of each other, a centerline for the colon in each of the 40 scans by using in-house 3-D computer software with a custom centerline tool. They were instructed to include all non-collapsed parts of the colon.
Implementation Issues
The following implementation details and settings are used in our implementation of the centerline algorithm. They serve as a guide for efficient implementation, but are not proposed to be the only way for implementing the algorithm.
DM. We use the chamfer 3-4-5 DT [21] for computing the distance map.
Seedscalefactor, Regularscalefactor. Scale factors that scale a spherical neighborhood of a seed point and regular point, respectively, used for removing other maxima within a local neighborhood. If the values for these parameters are set too high, maxima that happen to be close in 3-D space, but are located in different parts of the colon may accidentally be removed. If the parameters are set too low, the computation time increases without any increase in the accuracy of the resulting centerline. In our implementation, the size of the sphere for a seed point (seedscalefactor) is set to 1.5 times the physical sphere radius, as indicated by the DM value at that point. The value was chosen in order for the spherical neighborhood it represents to cover the colonic lumen including parts of the wall. The size of the sphere for a regular point (regularscalefactor) is set to 1.1 times the physical sphere radius. The seed points are considered more accurately positioned in 3-D space than the regular points, hence the larger value for seedscalefactor. Additionally, we set a minimum radius to 10 voxels for any sphere. This is not necessary for the algorithm to work. However, it will give a slight speed improvement.
Minseedsize. Smallest maxima size that will be regarded as seed point. It is typically set to a size within the range of 50-80 for 3-4-5 DT, but this may vary depending on the resolution of the colon in CTC scans. minseedsize should reflect the radius in voxels for the well-distended parts of the colon.
Minseedamount. Minimum number of seed points desired. This parameter may, if necessary, override minseedsize. It is typically set within the range of 30-70. It may be set to zero in order to reduce the number of parameters. The parameter ensures that no less than minseedamount seed points are selected, even if this implies that the parameter minseedsize must be overruled. The parameter is useful only for unusual cases in which a large portion of the colon is collapsed. In such a case, if this parameter is not used (i.e., set to zero), minseedsize would have to be lowered for all cases, which increases the risk of finding the centerline of non colonic structures.
Nondt. Highest voxel value in DM considered to be outside the colonic lumen. We typically use a nondt value of 3 for 3-4-5 DT. If Euclidean DT is to be used, a value of 1 would be appropriate. The partial volume effect can create holes between neighboring colonic walls, and thus a higher value of nondt can avoid the generation of centerlines through the holes. However, setting too high a value for nondt may prevent the extraction of the entire colon centerline. The parameter should be set considering the scale of the DT used to reflect correctly how many layers it peels off the surface in the DM.
Detection of Local Maxima in DM. The removal of close-by maxima can be implemented effectively in a two-step manner: First, store all maxima found in DM in a look-up table (LUT). Also, set up an empty (all zero) volume Vm with the same dimensions as those of the DM. Second, starting with the largest maxima in the LUT, for each maximum Mi check whether the corresponding voxel in Vm is zero. If so, set all voxels in Vm within a desired radius (i.e., the DM value of Mi scaled by seedscalefactor or regularscalefactor) to value 1. If not, remove Mi from the LUT.
Graph Connection. It is important to evaluate the iteration on only the recently added nodes in G (i.e., the nodes that have not previously been evaluated by Graph Connection). In addition, the iteration must not re-visit a node before visiting the next node in the list. Otherwise, the computation time for the Graph Connection process could easily take several minutes. It should be noted that we implemented the path search with a simple Bresenham line algorithm that checks the voxel values for all voxels on the path. If any non dt voxel is found, it returns FALSE, but returns TRUE otherwise.
Results: Evaluation Based on Colon Phantom
Two phantom studies were made. In the first study, a centerline was computed for phantoms with different radius, ranging from 3 to 30 voxels (2-21 mm). As seen in Table I, the centerline computation time is nearly constant (2.3-2.9 seconds), and the extension drops from 100% only in very narrow (3-5-voxel radius) and very well-distended cases (30-voxel radius). For narrow passages, this is because the centerline is split into multiple segments due to sharp turns in the colonic pathway. For a large radius, the decrease occurs because the ends of the original and computer-generated segments are not perfectly aligned. Thus, the similarity computation algorithm will regard these segments as missing. With increased phantom radius, the average displacement-related parameters (i.e., the mean displacement, the max displacement, and the standard deviation) of the computed centerline increases when compared to the path used for creating the phantom. This is because, when the phantom radius becomes large, at sharp turns close-by colonic walls will merge, thus shifting the actual central path of the colon from its original position (see
In the second study, the seedscalefactor was varied within the range of 0.1-2.0. The regularscalefactor was set to 0.7×seedscalefactor. The phantom radius was set to 15 voxels. As seen in Table II, all of the colon centerline was computed regardless of the seedscalefactor. Also, the average displacement-related parameters (i.e., the mean displacement, the max displacement, and the standard deviation) are close to constant. The computation time was close to constant except for the smallest scale factor, when the time consumption increased dramatically. The reason for this is that with such small scales, almost every possible maximum found in the DM will be used by the algorithm (8200 out of 8275). It should be noted that the creation of the graph structure is not linear but, for n nodes, in the worst case proportional to (n/2)2. This is because every node added must be compared to all existing nodes in the graph structure. The complexity of the remaining parts of the algorithm is proportional to the number of voxels in the volume data set or to the number of nodes in the graph structure.
To estimate whether a complexity of (n/2)2 will interfere with the performance of the algorithm under normal circumstances, we compute a reasonable upper estimate of the number of maxima picked up by the algorithm by using the following model: Assume that a colon C is cylinder-shaped, approximately 1500 mm long, and has an average radius of 25 mm [22]. Assume also that we have an isotropic voxel space and 1 voxel=1 mm3. C will consist of Vc=π252×1500≈2.94×106 voxels. Each maximum used by the algorithm will be the single maximum within a spherical neighborhood, and occupy vm=(⅓)×4πr3 voxels where r is the shortest distance to the object boundary. This implies that the seedscalefactor and the regularscalefactor parameters are both set to 1.0. If the number of maxima collected is the inverse of the squared maxima radius, and maxima radii are, on average, 15 mm, an upper limit of the total number of maxima used will be
As seen in Table II, column 2, the centerline can be computed in less than 3 seconds for this number of maxima, which indicates that the computation time is not likely to become large for any colon.
Results: Evaluation Based on Clinical Cases
The computer-generated centerlines for the 40 CTC scans were always regarded as being the slave centerlines when compared to the human-generated centerlines.
Table III shows the results of the comparison between computer-and human-generated centerlines, and between human-generated centerlines. As shown, all quantities for the computer-generated centerline, except for the extension, are comparable with or even better than the pair-wise comparison between human centerlines. The average maximum displacement was relatively large for all cases, both when human centerlines were compared and when human centerlines were compared with the computed centerline. However, most of the large displacements were observed in the end regions of either the cecum or rectum. In these regions, it was difficult to determine the centerlines uniquely because of the large and irregular shapes of the regions.
Table III also shows that the extension of the computer generated centerlines, on average, was 92% of the human-generated centerlines, comparable to 94% when human-drawn centerlines were compared with each other. By computation of the same quantity for the three groups (straightforward, moderate, and challenging) separately, it is shown in Table IV that, for the straightforward cases, the average extension for both computer-generated centerlines and human-drawn centerlines were close to 100% of the entire colon (97% and 98%, respectively). For the moderate cases, the average extension of the computer-generated centerlines was lower than that obtained from the straightforward cases (91% and 97%, respectively). For the challenging cases, the extension of the computed centerlines was, on average, 81% of the colon, compared to 84% extension by the human-drawn centerlines. The computer generated centerlines covered less of the colon relative to the human-generated centerlines because some small, very narrow segments of the colon were found by the radiologists, but were missed by the algorithm. This was expected because these narrow segments were considered to be non-colonic objects by the algorithm. However, even the centerlines drawn by radiologists had poor extension when compared to each other. This is mainly for the same reason as above: the challenging cases often contained a large amount of stool, collapsed regions, and extra-colonic findings which made it difficult for the radiologists to find the true colonic path. Thus, the centerlines drawn by each radiologist covered different parts of the colon.
A two-sample t-test was performed to evaluate the difference in the displacements between human-drawn centerlines and between human-drawn and computer-generated centerlines. The null hypothesis was that there was no difference. To explore this hypothesis, all displacements between human-drawn centerlines used for computing the averages in Table III, were pooled and compared to those of all displacements between human-and computer-generated centerlines. The null hypothesis was not rejected at the =0.05 level. This result indicates that there was no statistically significant difference in the displacements between the human-drawn centerlines and between human-drawn and computer-generated centerlines. However, the average extension of the computer-generated centerline was statistically significantly less than that of the human-drawn centerlines.
For evaluating the speed, we measured the CPU time for computing the centerlines. We did not include the calculation time for the interpolation of the input data or for the segmentation of the colon. Using an Intel Pentium-based 800 MHz computer with 512 MB memory, for the 40 scans, the centerline computation time was in the range of 2.6-7.4 seconds, with an average of 4.8 seconds. The distance transform used in this study was computed, on average, in 5.7 seconds, giving a total centerline computation time of 10.5 seconds, on average, for the 40 CTC scans.
In 4 of the 40 cases, the computer-generated centerline contained a total of 5 centerline segments that did not correspond to any segment of the centerline drawn by a human observer. These segments respectively represented 0.9%, 1.0%, 1.4%, 1.6%, and 6.7% of the total centerline, within a path length range of 18-148 mm.
Visual inspection confirmed that, in no case, the computer generated centerline went through any hole created by the partial volume effect or its equivalent.
The proposed centerline algorithm has several advantages. First, to our knowledge, it is the fastest implementation among the methods for computation of colon centerlines reported in the literature.
Second, our algorithm performs a global search for selection of seed points, without any human intervention. The global search is designed to find seed points in all of the distended segments of the colonic lumen even when these segments are disconnected. Therefore, the algorithm can generate centerlines without being trapped by collapsed regions.
Third, the algorithm is robust to segmentation of the colon of non-optimal quality, with collapsed regions, and stool within the colonic lumen; it computes all available centerlines, regardless of the location and length of individual segments. The set of globally extracted seed points determines which segments belong to the colon.
In addition to these advantages, the algorithm has an important feature. It uses 6- and 4-neighbor maxima of the DM as a guide for positioning the centerline. These maxima are, by definition, well centered in the DM, and thus the resulting centerline is also centered.
The main reason for the speed of the algorithm is data reduction. Initially, a set of maxima representing the centerline pathway is extracted from the DM, reducing the amount of data by a factor of 100,000 compared to the DM. Only the reduced data set is used for creation of the final centerline. The fact that the algorithm does not depend on a specific type of DT, but can be used with an integer-based DT, reduces the total computation time even further [21].
Because the algorithm attempts to find segments of the centerline in different parts of the abdomen independently, application of the algorithm to poorly segmented colons may generate centerline segments for non-colonic structures, such as the small bowel and the lungs. In most of these cases we observed that the entire colon centerline was extracted and could easily be identified as being the longest segment. However, if only parts of the entire colon centerline were generated, anatomy-based post-processing might be necessary for the centerline segments to be connected correctly. This, however, is beyond the scope of the present study.
Several criteria were evaluated for recovering the end segments of a centerline in the Branch Recovery step. Generally, long branches appear to be a good representation of the lost end segments of the centerline. It is, however, not guaranteed that this gives the “correct” end segments. Our method described in Section Branch Recovery was able to recover the most extended branch that does not turn back to the extracted centerline, and thus found to be most appropriate method.
As shown in Table III and IV, for all centerline comparisons, the average maximum displacement was at least 15 mm. Such a large deviation often occurred at the end points. This can be explained in terms of uncertainty. That is, neither the radiologists nor the computer algorithm knew exactly in which direction the end point should be, not even whether the centerline was supposed to be drawn all the way to the colonic wall indicating the end of cecum, rectum, or an in-between segment, or to stop somewhat earlier.
As stated above, the computed centerlines, on average, yielded slightly less extension compared to the human generated centerline because narrow segments between collapsed regions were sometimes completely missed by the algorithm. This was expected, because the algorithm could not find seed points anywhere within such a narrow region, and thus the segments were not regarded as part of the colon. It is possible to lower the minimum-size threshold value for selection of seed points to find these narrow segments although it may potentially increase the number of extra-colonic centerlines.
It should be noted that it do exist centerline algorithms which can produce centerlines that are, on average, more well centered with respect to a DM than the algorithm presented in this paper. The reason for this is, the centerline is based on a set of discrete, centered, nodes within the colon. In cases where the Euclidean distance between connected nodes is large, the in-between connection might be off the central path.
III. Detection of Colorectal Masses in CT Colonography
Because the mass detection methods were designed to be integrated in the inventors' previously developed automated CAD scheme for the detection of polyps, a brief overview of that scheme is provided. The specific details are presented elsewhere [33,34,35,36,49].
Given a 3-D isotropic CTC data set, a thick region encompassing the colonic wall is extracted automatically by use of a knowledge-guided segmentation technique [23]. In this technique, the region outside the body, the osseous structures, and the lungs are first eliminated from further analysis by use of histogram analysis, thresholding, and morphological operations (
The regions with polyp-like characteristics are detected from the extracted region of the colonic wall by use of hysteresis thresholding [33,57]. The hysteresis thresholding is guided by two volumetric features called the shape index (SI) and the curvedness (CV) [29], which can be defined as
where k1(p) and k2 (p) are the principal curvatures of the iso-surface passing through voxel p [29]. The SI classifies the topological volumetric shape at a voxel, and the CV characterizes the scale of that shape. To implement the hysteresis thresholding, one first locates the voxels that have SI and CV values within SI∈[0.9,1.0] and CV∈[6.68,8.52]. The regions represented by such voxels are expanded by addition of 26-connected boundary voxels that have their SI and CV values within SI∈[0.8,1.0] and CV∈[6.40,8.52]. The details of this technique have been presented elsewhere [33].
The regions extracted by hysteresis thresholding, or initial detections, do not necessarily cover completely the visually perceived region of the detected lesions. For analyzing the internal structure of the detections for FP reduction, the final regions of the detections, or the polyp candidates, are extracted by use of a conditional morphological dilation technique [35]. In this technique, a detected region is expanded within the colonic wall by morphological dilation, and the dilation step for which the growth rate of the expanding region is smallest is chosen to represent the final region of the detection [35].
Finally, FP polyp candidates are reduced by use of a quadratic classifier based on quadratic discriminant analysis. The classifier is trained by labeling of each polyp candidate as a true-positive (TP) or FP detection based on the feature values calculated for the candidate region. The quadratic classifier then generates a decision boundary that partitions the feature space by use of a hyperquadratic surface. The result of the partitioning can be represented by a discrimination function, which projects the multi-dimensional feature space into a scalar decision variable. For a polyp candidate, the value of the decision variable determines whether the candidate belongs to the TP or FP class [33]. The polyp candidates in the TP class yield the final set of the polyps detected by the CAD scheme.
Overview of Mass Detection
The first method for mass detection, fuzzy merging, identifies intraluminal masses that protrude from the colonic wall into the colonic lumen and is described below. The second method, wall-thickening analysis, detects non-intraluminal types of masses that appear as an abnormal thickening of the colonic wall. This method is described in more detail below. Once the locations of the mass candidates have been established, we a level set method is used to extract the complete regions of the mass candidates, as described below.
To reduce FP candidates, shape and texture features were calculated from the detected mass regions. A quadratic classifier discriminates the mass candidates into TPs and FPs as in the case of polyp candidates. However, the classifier is trained separately for polyps and the two types of masses, because these lesions have considerably different feature characteristics. Once the final candidates are established, the polyp candidates that overlap the regions of mass candidates are eliminated, because such polyp candidates are likely to be FP detections of a mass. The final output of the integrated CAD scheme for the mass and polyp detection is obtained by combining the outputs of the quadratic classifiers from the analyses of the polyp and mass candidates.
The methods for the detection and extraction of the mass candidates involve assignment of several numerical parameter values. The assignment of these values is discussed briefly below.
Detection of Intraluminal Masses: Fuzzy Merging
As discussed above, the CAD polyp detector may generate several initial detections on the surface regions of large unobstructed intraluminal masses. Therefore, to reduce the computational burden introduced by mass detection, a separate detection step for intraluminal masses is not implemented, but they are identified by use of a new fuzzy merging method that analyzes the fuzzy memberships of the initial polyp detections. The detections that are determined to originate from the same lesion are considered as an initial mass detection and are processed separately from the polyp detections.
Let U be a universal set. A fuzzy set F in U can be defined as a set of ordered pairs as
F={(c,μF(c))|c∈U}, (10)
where μF(.) is the membership function (or characteristic function) of F, and μF(c) is the grade (or degree) of membership of c in F [58]. The membership function μF(.) maps U to the membership space M, where usually M=[0, 1]. In a crisp set, M={0,1}. An α-cut of a fuzzy set F is a crisp set Fα that contains all elements of the universal set U with a membership grade in F greater than or equal to α:
Fα={c∈U|μF(c)≧α}, α∈(0,1]. (11)
In our application, the input data are the set of 26-connected regions of the initial CAD polyp detections. The regions that are located beyond Dmax=50 mm from each other may not be merged, and their membership grade is always zero. Otherwise, for two regions, A and B, within Dmax, the symmetric membership function is defined by
μA(B)=μB(A)=(wCTμA,BCT+wGRμA,BGR)εD(A,B), (12)
where the functions μA,BCT and μA,BGR characterize the membership between A and B based on the CT and gradient values along a 3-D sampling ray R between the centers of A and B. Here, the sampling ray R is defined as a line between the geometric center voxels of A and B. As described below, the value of μA(B) is defined to become high when regions A and B are on a single lesion. The weights wCT and wGR are generally set to (wCT,wGR)=(0.5,0.5), but may be modified in situations where it is obvious that A and B should not be merged. The function εD(A,B) decreases the membership grade as the distance between A and B increases, as discussed below.
To calculate μA,BCT, the minimum CT value along R between the regions of A and B, denoted by minCT(A,B), is first computed (
is scaled linearly to [0,1] in order to characterize the membership of A and B based on the CT value:
The value μA,BCT is high if the CT values along R are relatively uniform, indicating that the detections appear within the same lesion. The value is low if the CT value is attenuated significantly between the regions, suggesting that the detections may not appear within the same lesion.
To calculate μA,BGR, one determines the maximum gradient between the regions A and B, or maxGR(A, B), and divides it by the maximum value of the gradient within A and B, or maxGR[A,B]. The resulting value,
is scaled as follows:
The value μA,BGR is high if the highest gradient between A and B is low, i.e., there is no clear separating boundary between the regions (
The numerical boundary values in Eqs. (13) and (14) were determined experimentally, as discussed in Section 4. For example, in Eq. (13), the threshold value 0.9 implies that if the CT value is reduced less than 10% between A and B, the two regions A and B definitely belong to the same lesion. Similarly, if the CT value is reduced by more than 50% between A and B, the regions definitely do not belong to the same lesion.
If the CT value between A and B is less than −700 HU, it is likely that there is air between the regions and they do not represent the same lesion. However, the value of μA(B) may not always reflect this clearly if equal weights are used. Therefore, in this situation, the effect of the CT value is weighted by modifying the membership weights in Eq. (12) as follows:
The membership grade between A and B is also dependent on their distance. The effect of εD(A,B) is to reduce the membership grade for detections that are far apart as compared with detections that are very close. As the distance D(A,B) between A and B increases, it becomes increasingly likely that they originate from a different lesion, unless the membership grade based on the CT value and gradient is very high (
Once the memberships have been established, the detected regions that have a high membership grade are merged by determining the α-cut of the regions. The threshold for this hard clustering is α=0.5, but the regions with the highest membership are merged first. If, during this merging process, one or both of the regions A and B to be merged already belongs to a previously merged region, e.g., A⊂E, where E is a merged region, it is checked that the distance between all components of E and B (or all components of B, if B is also a previously merged region) are within Dmax. Because of this constraint, the merging process cannot spread over a long “chain” of candidates along the colonic wall. On the other hand, masses with complicated shapes can be detected, because a high membership grade is required only between two detections within a region (
The merging process ends when no mergeable regions can be found. The regions that were constructed by merging of initial detections are tested for their spatial spread along their three perpendicular principal axes for estimating the sizes of the lesions that the merged regions represent. A region is considered a mass candidate if the spread exceeds 25 mm along one of the principal axes, or if the largest of these spreads is 22.5-25 mm and the spreads in the two corresponding perpendicular directions are >50% and >75% from the largest spread. Small polyp-like lesions have a small maximum spread, or their spreads in the directions perpendicular to the maximum spread are small. The detections that are not labeled as mass candidates are considered polyp candidates and are processed as described in below.
Detection of Non-Intraluminal Masses: Wall-Thickening Analysis
As shown in
First, the CT values along the sampling ray L with a length of 20 mm are averaged by a 3-voxel window to reduce noise. Let the start and end locations of L be denoted as Lbegin and Lend (
Once the voxels representing potential regions of abnormal thickening have been detected, a symmetry check is performed by projecting an inverse sampling ray from the detected voxels through the colonic lumen to the opposite side of the colonic wall (
To reduce false positives, additional checks are performed for the extracted, suspicious regions. The total volume of the region should be at least 3000 mm3, and there should be at least 300 voxels associated with the locations where suspicious thickening was originally detected. Furthermore, to eliminate FP detections due to fluid, the directions of the gradient vectors within the suspicious region were analyzed. If more than 50% of the region is associated with unit vectors of gradient with a sagittal component ≧0.9 in magnitude, the region is considered to be fluid and is eliminated from further analysis (
Extraction of the Mass Region by Use of Level Sets
The mass candidate regions, detected by the two methods described in the previous sections, are used as seeds for extracting the final mass regions. The extraction step is based on a level set method [59]. The basic idea of the level set method in extracting regions is to allow the initial surface γ of a smooth, closed object to move along its normal vector field with speed F. In this application, the initial regions of the mass candidates provide the initial conditions for the level set method. To extract the complete region of a mass candidate, solve the following flow evolution equation:
φ+(gI+cI)(1−εκM)|∇φ|−β∇P·∇φ=0, (17)
where φ is the level set function to be solved. The zero level set of φ represents the evolving front, or surface, of the extracted region. Ideally, the surface should pass the minor boundaries in the input data and settle to the major boundaries. The front can be considered to expand with a speed 1−εηM, where ηM is the mean curvature of the surface, and ε controls the effect of the curvature. The additional terms in Eq. (8) represent expansion or contraction forces that control how the level set surface evolves.
The gradient-based expansion force gI is
where Gσ*I denotes the CTC data set convolved with a 3-D Gaussian smoothing filter with a standard deviation σ=0.5. If there is no boundary gradient, the expansion force is gI=1.
In addition to the gradient-based expansion force, we introduce an expansion force cI based on the CT value (in HU):
In the present invention, the use of cI is necessary because the automatically detected initial regions may be located at both sides of the mass boundary. Therefore, without cI, the evolving level set front could spread into the colonic lumen.
The surface tension force, |∇φ|, controls the smoothness of the evolving level set surface, and depends on the mean curvature which may be calculated from the level set function φ as
The force that attracts the evolving surface toward a boundary is implemented as the gradient of a potential field defined by
P(x,y,z)=−(max{0,|∇(Gσ*I(x,y,z))|−20}). (21)
The parameter β in Eq. (15) controls the effect of the force that attracts the evolving surface to a boundary.
The flow evolution of Eq. (15) was implemented similarly to the two-stage approach described by Sethian et al. [59]. First, the flow evolution was solved by use of the expansion forces and a fast marching method. The fast marching method was realized by solving of the Eikonal equation
|∇T|gI=1,Γ={(x,y)|T(x,y)=0}, (22)
where Γ is the initial position of the level set surface. The iteration continues until the solution |∇T| becomes large and, hence, the flow has almost stopped because it reached the desired boundary. To implement the fast marching method, write the Eikonal equation as
where the terms DijkT represent approximate forward and backward operators on T. For example,
The fast marching method can be implemented efficiently by observing that the upwind difference structure of Eq. (23) allows propagation from small to large values of T. Only the values around the thin evolving surface boundary need to be considered, and these can be updated by use of a min-heap data structure. The details of the implementation are discussed in Sethian et al. [59].
Once the final level set surface has been approximated by the fast marching method, a more precise, narrow band level set method could be used for obtaining the final solution [59]. However, in this application, the fast marching method was applied a second time by including the boundary-attracting force in the second term of Eq. (15). Three examples of the final extracted regions of masses are shown in
Selection of the Parameter Values and Features
Most of the numerical parameter values of the methods of mass detection and extraction were determined empirically. For the fuzzy merging method, the numerical threshold values in Eqs. (13)-(15) were determined by examination of the CT and gradient values between the mass and polyp CAD detections in one data set that included five polyps and a circumferential mass. For the wall-thickening analysis, most of the parameter values were estimated visually from one dataset containing a typical wall-thickening type of mass. The requirements that a detected region should have a minimum volume of 3000 mm 3 and be associated with ≧300 voxels indicating wall thickening were established by first calculating these values for both non-intraluminal masses: the corresponding minimum values for all data sets were approximately 3600 mm3 and 380 voxels, respectively.
The values of the parameters ε=0.4 and β=0.66 of the level set method, as well as the threshold values in Eqs. (16), (17), and (19), were adjusted approximately by a visual examination of the resulting extracted mass regions.
A combination of two features was used to reduce FP mass detections by use of a quadratic classifier. The use of more than two features would improve the performance of the mass detection, but would also reduce the generalizability of the results. Because only fourteen masses were used in this study, the use of additional features might cause the classifier to over learn the data. With only two features, it is likely that the results obtained in this study can be reproduced when a larger database is used for the evaluation of our mass detection methods.
The features were chosen by evaluation of the FP reduction for combinations of the mean, variance, and the mean of the 10 maximum values of the CT value, the shape index, and the gradient within the regions of the detected mass candidates. The calculation of these features had been implemented previously in the CAD scheme for the detection of polyps. The features were designed to characterize the shape or internal texture of the detected regions for an effective reduction of FPs.
Performance Evaluation Method
The detection results were characterized primarily with a by-mass analysis, where a mass is considered detected if it is detected in either or both supine and prone positions of the patient. Because each abnormal patient had only one mass, the results of the by-mass and by-patient analysis were identical. In a by-patient analysis, an abnormal patient is considered detected if at least one mass is detected in the supine or prone view of the patient. Furthermore, we calculated the detection sensitivities when only prone data sets or only supine data sets were used in the evaluation.
The performance of the detection of intraluminal masses was evaluated by use of the leave-one-out (round-robin) method, which is expected to provide a good estimate of the unbiased performance of a classifier when the available data are limited. The leave-one-out evaluation was implemented with a by-patient elimination. That is, for a given patient, all candidates detected in the patient were removed from the set of all candidates, the quadratic classifier was trained with the remaining candidates, and the resulting discrimination function was tested on the candidates that were removed. This has the effect of reducing the potential bias caused by the use of the data from the same patient for training and testing. Once this process had been repeated for each patient, the discrimination function was thresholded at various detection sensitivity levels to generate a free-response operating characteristic (FROC) curve [61].
Because of the small number (2) of non-intraluminal masses in our database, it was not feasible to perform a leave-one-out evaluation for the wall-thickening analysis. Therefore, for this method, the quadratic classifier was trained with all detections of the non-intraluminal mass detection method and evaluated the result.
The mass detection performance of the CAD scheme was also investigated for the detection of polyps without the use of the explicit mass detection methods. That is, the inventors investigated how many masses would have been detected if the CAD scheme for polyp detection had been used. The quadratic classifier of the CAD scheme for polyp detection was trained and tested with the 82 cases and a leave-one-out method by treating the polyps (5-25 mm) as TPs and other detections as FPs. Then, the number of masses (>25 mm) that were detected as polyps by this CAD scheme at each detection sensitivity level was determined.
We retrospectively collected 82 cases (patients) from standard CTC examinations performed at the University of Chicago Hospitals (80 cases) and the Beth Israel Deaconess Medical Center (2 mass cases). In preparation for same-day optical colonoscopy, all patients underwent a colon-cleansing regimen with polyethylene glycol. The CTC scanning was performed with helical single- and multi-slice CTC scanners (HiSpeed CTi or LightSpeed QX/i, General Electric Medical Systems, Milwaukee, Wis.) in supine and prone positions with reconstruction intervals of 1.0-5.0 mm (within 1.5-2.5 mm: 91%), slice thicknesses of 1.25-5.0 mm (within 5.0 mm: 80%; within 2.5 mm: 15%), tube currents of 50-160 mA (within 100 mA: 87%) with 120 kVp, and tube rotation times of 1.0 (single scan) and 0.8 (multi-scan) seconds. The CTC data sets represented the entire region of the colon, from the rectum to the lung diaphragm. To yield 3-D isotropic data sets, the original 512×512 CT images were interpolated linearly in the axial direction to a voxel resolution of 0.51-96 mm (average: 0.69 mm).
There were 14 colonoscopy-confirmed masses (12 intraluminal, 2 non-intraluminal) measuring approximately 30-60 mm in diameter (
There were 30 colonoscopy-confirmed polyps measuring 5-25 mm in 19% of the cases. Four cases had both polyps and masses, resulting in a total of 15 polyps and 4 masses. Fifty-six cases (68%) were normal, i.e., a complete optical colonoscopic examination did not detect any colorectal neoplasms in these cases.
The fuzzy merging method detected all but one of the intraluminal masses, or 79% of all masses. After the initial detection step, there were 165 FP initial mass candidates. The two features used to reduce FPs are shown in
The wall-thickening analysis detected both of the non-intraluminal types of masses in the database, as well as one of the intraluminal circumferential masses with apple-core morphology. The latter mass resembled the wall-thickening type of non-intraluminal mass, since it obstructed the colon circumferentially. There were 69 FP initial mass candidates, but most of these could be eliminated by use of the quadratic classifier with two features (
The FROC curve obtained in this manner is expected to represent the combined mass detection performance well, because the two mass detection methods run independently and the number of false positives generated by the wall-thickening analysis is negligible. When only prone data sets of the patients were used in the mass detection, the fuzzy merging method detected 10 of the 12 (83%) intraluminal masses, and the wall-thickening method detected 2 of the 2 (100%) non-intraluminal masses. In combination, the two methods detected 12 of the 14 (86%) masses in prone data sets with 0.10 FPs per patient on average. When only supine data sets were used, the fuzzy merging method detected 8 of the 12 (67%) intraluminal masses, and the wall-thickening analysis detected 100% of the non-intraluminal masses. In combination, the methods detected 10 of the 14 (71%) masses in supine data sets with 0.61 FPs per patient on average.
The mass that was missed by both mass detection methods was located at the rectum and was partially cut off from the CTC data set in both supine and prone positions (
To verify that the mass detection methods work consistently, and to identify potential problems, we also analyzed the sources of the FP detections produced by both mass detection methods. The FP polyp detections generated by the CAD method for the detection of polyps have been analyzed similarly in previous studies [49,62] by the present inventors. The mass detection methods generated only few false positives. The fuzzy merging method generated 13 FP mass detections for the 164 data sets. There was no particular single source of FPs, but 43% of the FPs were associated with a rectal wall with bump-like structures caused by a combination of residual stool, hemorrhoids, and extrinsic compression by the rectal tube. Other sources included a dilated fold with stool or a dilated ileocecal valve (31%) (
Finally, we evaluated the mass detection performance of the CAD scheme developed for the detection of polyps alone, as described below. At a 90% by-polyp detection sensitivity level, only two of the masses (14%) were among the final output of the scheme. At a 97% by-polyp detection sensitivity level, five (36%) of the masses were among the polyp detections. Therefore, the CAD scheme for polyp detection cannot be considered effective for the detection of large masses.
Both mass detection methods missed a mass that was partially cut off from the CTC data in both supine and prone positions. This suggests that a specialized method may need to be developed for checking for abnormalities within the boundary regions of the CT data reconstructed from axial images. On the other hand, both methods were able to detect masses that were only partially visible due to a colonic segment collapsing on the mass. In these cases, the complete region of the mass had been imaged by the CTC scanner, and, thus, an analysis of the internal structure of the mass was performed, leading to the detection of the mass.
In the case of intraluminal masses, the highest detection performance was obtained when both supine and prone data sets were used. This is because some masses that were clearly visible in one position of the patient were invisible in the other position due to a collapsed region or residual fluid, as shown in
The number of masses used in the evaluation was relatively small (14 masses in 14 patients). In particular, most of the masses were of the intraluminal type, and only two of the masses were of the non-intraluminal type which was detected by an analysis of the thickness of the colonic wall. Therefore, even though the evaluation was based on a leave-one-out evaluation with a large number (164) of CTC data sets, it is premature to make precise estimates about the performance of our CAD scheme when it is used with a larger population. Nevertheless, it appears that the detection of masses does not affect the result of the polyp detection, and each type of lesion was identified correctly and processed separately by the corresponding method. In particular, each detected polyp was identified correctly as a polyp, and each detected mass was identified correctly as a mass. The initial results suggest that a high sensitivity and a low FP rate may be obtained with automated detection of masses.
The leave-one-out evaluation was performed by use of by-patient elimination, where all detections of a patient are removed at each evaluation step, one patient at a time. The results are expected to estimate the performance of the CAD scheme in a situation where an unseen case is presented to the scheme. On the other hand, because of the relatively small number of masses and polyps, some of the CTC cases are unique in the sense that similar lesions are not found in the other cases. Therefore, the leave-one-out performance with by-patient elimination might give an overly conservative estimate of the performance of the CAD scheme.
Although the current implementation of the CAD scheme has not been optimized for speed and efficiency, it is relatively fast. The complete process of detecting polyps and masses starting from the original CT images in DICOM format and including disk access takes approximately 10 minutes when a PC workstation (Marquis K120-SR, ASL Inc., Newark, Calif.) with dual processors (Athlon MP 1.2 GHz, AMD, Sunnyvale, Calif.) is used. The detection of intraluminal masses by the fuzzy merging method takes only a few seconds, the detection of non-intraluminal lesions takes about 30 seconds, and the level set extraction takes about 30 seconds per mass candidate.
Two methods for the detection of colorectal masses were developed. The methods can be integrated into a CAD scheme for the detection of polyps for computational efficiency. All masses but one that was partially cut off from the CTC data sets were detected at a 93% sensitivity with an average FP rate of 0.21 detections per patient. No polyp was missed due to mass detection when a polyp detection scheme was combined with the mass detection. Integration of the computerized polyp and mass detection in a single CAD scheme appears to be a feasible approach to provide a comprehensive high-performance computer aid for radiologists in a clinical CTC screening setting.
The present method can also be implemented more generally on other medical images of other organs (e.g., mammographic breast images, or CT scans of the thorax, abdomen, or skeletal system) with respect to some other disease state or state of risk. Nodule or lesion feature values can readily be obtained from other medical images by those of ordinary skill in the art. For example, the detection of nodule or lesion feature values in various medical images is also well known in this art. See, e.g., U.S. Pat. No. 5,881,124 (Giger et al., Automated method and system for the detection of lesions in medical computed tomographic scans), the contents of which are incorporated herein by reference.
For the purposes of this description an image is defined to be a representation of a physical scene, in which the image has been generated by some imaging technology: examples of imaging technology could include television or CCD cameras or X-ray, sonar, or ultrasound imaging devices. The initial medium on which an image is recorded could be an electronic solid-state device, a photographic film, or some other device such as a photostimulable phosphor. That recorded image could then be converted into digital form by a combination of electronic (as in the case of a CCD signal) or mechanical/optical means (as in the case of digitizing a photographic film or digitizing the data from a photostimulable phosphor). The number of dimensions that an image could have could be one (e.g. acoustic signals), two (e.g. X-ray radiological images), or more (e.g. nuclear magnetic resonance images).
All embodiments of the present invention conveniently may be implemented using a conventional general purpose computer or micro-processor programmed according to the teachings of the present invention, as will be apparent to those skilled in the computer art. Appropriate software may readily be prepared by programmers of ordinary skill based on the teachings of the present disclosure, as will be apparent to those skilled in the software art.
As disclosed in cross-referenced U.S. patent application Ser. No. 09/773,636, a computer 900 may implement the methods of the present invention, wherein the computer housing houses a motherboard which contains a CPU, memory (e.g., DRAM, ROM, EPROM, EEPROM, SRAM, SDRAM, and Flash RAM), and other optional special purpose logic devices (e.g., ASICS) or configurable logic devices (e.g., GAL and reprogrammable FPGA). The computer also includes plural input devices, (e.g., keyboard and mouse), and a display card for controlling a monitor. Additionally, the computer may include a floppy disk drive; other removable media devices (e.g. compact disc, tape, and removable magneto-optical media); and a hard disk or other fixed high density media drives, connected using an appropriate device bus (e.g., a SCSI bus, an Enhanced IDE bus, or an Ultra DMA bus). The computer may also include a compact disc reader, a compact disc reader/writer unit, or a compact disc jukebox, which may be connected to the same device bus or to another device bus.
Examples of computer readable media associated with the present invention include compact discs, hard disks, floppy disks, tape, magneto-optical disks, PROMs (e.g., EPROM, EEPROM, Flash EPROM), DRAM, SRAM, SDRAM, etc. Stored on any one or on a combination of these computer readable media, the present invention includes software for controlling both the hardware of the computer and for enabling the computer to interact with a human user. Such software may include, but is not limited to, device drivers, operating systems and user applications, such as development tools. Computer program products of the present invention include any computer readable medium which stores computer program instructions (e.g., computer code devices) which when executed by a computer causes the computer to perform the method of the present invention. The computer code devices of the present invention may be any interpretable or executable code mechanism, including but not limited to, scripts, interpreters, dynamic link libraries, Java classes, and complete executable programs. Moreover, parts of the processing of the present invention may be distributed (e.g., between (1) multiple CPUs or (2) at least one CPU and at least one configurable logic device) for better performance, reliability, and/or cost. For example, an outline or image may be selected on a first computer and sent to a second computer for remote diagnosis.
The invention may also be implemented by the preparation of application specific integrated circuits or by interconnecting an appropriate network of conventional component circuits, as will be readily apparent to those skilled in the art.
Moreover, parts of the processing of the present invention may be distributed for better performance, reliability, and/or cost. For example, an outline or image may be selected on a first computer and sent to a second computer for remote diagnosis.
Moreover, the invention provides computer program products storing program instructions for execution on a computer system, which when executed by the computer system, cause the computer system to perform the methods described herein.
Numerous modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described herein.
A, B, C = human-generated centerlines.
Computer = computer-generated centerline.
A, B, C = human-generated centerlines.
Computer = computer-generated centerline.
The present application is related to U.S. Patent Application Ser. No. 60/514,599, filed Oct. 28, 2003, and U.S. patent application Ser. No. 10/270,674, filed Oct. 16, 2002, the contents of which are incorporated herein by reference.
The present invention was made in part with U.S. Government support under USPHS Grant No. CA095279. The U.S. Government may have certain rights to this invention.
| Number | Date | Country | |
|---|---|---|---|
| 60514599 | Oct 2003 | US |
