Advanced computer-aided diagnosis of lung nodules

Abstract

Methods are herein provided for decision support in diagnosis of a disease in a subject, and for extracting features from a multi-slice data set. Systems for computer-aided diagnosis are provided. The systems take as input a plurality of medical data and produces as output a diagnosis based upon this data. The inputs may consist of a combination of image data and clinical data. Diagnosis is performed through feature selection and the use of one or more classifier algorithms.

Claims

1. A method of providing decision support in diagnosis of a disease in a subject, the method comprising: receiving a multi-slice computed tomography dataset; extracting, with a hardware processor, a volume of interest from the multi-slice computed tomography dataset; segmenting, with the hardware processor, the multi-slice computed tomography dataset to delineate at least one lesion within the volume of interest from a background of the multi-slice computed tomography dataset; extracting, with the hardware processor, image-based features from the at least one lesion in two and a half or pseudo-three dimensions over a set of two or more slices of the multi-slice computed tomography dataset, wherein the image-based features include a numerical feature; converting, with the hardware processor, clinical information into multiple categories and from a categorical format into a numerical format of the image-based features and; combining the converted clinical information with the image-based features, wherein the combined converted clinical information and the image-based features form a feature pool; and selecting, with the hardware processor, an optimal feature subset from the feature pool, wherein the selecting including selecting both image-based features and converted clinical information based on a training dataset; creating, with the hardware processor, a classifier or a committee of classifiers using the selected optimal feature subset and training dataset; and employing the classifier or committee to provide a prediction of a diagnosis of the at least one lesion.

2. The method according to claim 1, wherein the numerical feature represents one from a group comprising a coefficient, a moment, a gray scale value, a derivative of an image intensity, and an image intensity.

3. The method according to claim 2, wherein the clinical information includes at least two categories, male and female, and each of these categories is converted into the numerical format of the image-based features.

4. The method according to claim 1, wherein selecting an optimal feature subset further comprises using at least one of a correlation filter, a recursive feature elimination, and a random feature selection.

5. The method according to claim 4, wherein selecting an optimal feature subset comprises executing a genetic algorithm-based feature selection for a plurality of iterations with a different randomly selected set of training data and testing data to obtain resulting plurality of feature subsets.

6. The method according to claim 1, wherein the image-based features include the fractal dimension of the shape of the at least one lesion, an intensity and a size of dark internal pixel clusters in the at least one lesion, an intensity and a size of bright internal pixel clusters in the at least one lesion, chain code measurements of the at least one lesion, and at least one neighborhood gray-tone difference matrix of the least one lesion.

7. The method according to claim 1, wherein the classifier or committee of classifiers is a linear discriminant analysis.

8. The method according to claim 1, further comprising the classifier or committee of classifiers determining whether the at least one lesion is malignant or benign.

9. The method according to claim 1., further comprising the classifier or committee of classifiers determining a likelihood of malignancy of the at least one lesion.

10. The method according to claim 8, where the determination by the classifier or committee of classifiers is determined by at least one calculation selected from the group consisting of: a simple mean value, a simple vote, a weighted mean value, and a weighted vote.

11. The method according to claim 1, wherein the clinical information comprises a size of lymph nodes.

12. The method according to claim 1, wherein the clinical information comprises a number of lesions around the at least one lesion.

13. The method according to claim 1, wherein the clinical information comprises a location of the at least one lesion.

14. The method according to claim 1, wherein the clinical information comprises a presence of nodes.

15. The method according to claim 1, wherein the clinical information comprises a gender.

16. The method according to claim 1, wherein the clinical information comprises a past pulmonary disease.

17. The method according to claim 1, wherein the clinical information comprises a present pulmonary disease.

18. The method according to clan 1, wherein the clinical information comprises a prior chest surgery.

19. The method according to claim 1, wherein the clinical information comprises an occupational exposure.

20. The method according to claim 1, wherein the clinical information comprises a recreational exposure.

21. A computer-aided diagnosis system, the system comprising: a processor configured to: receive a multi-slice computed tomography dataset; extract an image-based region of interest from the multi-slice computed tomography dataset; delineate at least one lesion in the region of interest; extract image-based features from the at least one lesion in two and a half or pseudo-three dimensions over a set of two or more slices of the multi-slice computed tomography dataset from the region of interest, wherein the image-based features include a numerical feature; accept and convert clinical information from a categorical format into multiple categories and into a numerical format of the numerical feature and form a feature pool in which the clinical information and the image-based features constitute the feature pool; select an optimal feature subset from the feature pool and containing both image-based features and converted clinical information based on a training dataset; create a classifier or a committee of classifiers using the selected optimal feature subset and training dataset; and employ the classifier or committee to provide a prediction of a diagnosis of the at least one lesion.

22. The computer-aided diagnosis system according to claim 21, wherein the processor further selects, during training of the system, relevant image-based features and clinical features from a feature pool.

23. The computer-aided diagnosis system according to claim 22, wherein the processor utilizes an optimization technique selected from the group consisting of at least one of a correlation filter, a recursive feature elimination, and a random feature selection; and a classifier or committee of classifiers selected from the group consisting of at least one of a support vector machine, a decision tree, linear discriminant analysis, and a neural network.

24. The computer-aided diagnosis system according to claim 21, wherein the image-based region of interest is constructed by performing at least one morphological operation.

25. The computer-aided diagnosis system according to claim 24, wherein the at least one morphological operation is selected from the group of erosion, and dilation.

26. The computer-aided diagnosis system according to claim 21, wherein the image-based region of interest is constructed by selecting the largest contiguous object resulting from the segmentation.

27. The computer-aided diagnosis system according to claim 21, wherein the image-based region of interest further comprises an internal region, an external region, and a boundary region wherein unwanted elements are excluded from the region of interest.

28. The computer-aided diagnosis system according to claim 21, wherein the clinical information further comprises at least one data point selected from the group of age, gender, smoking history, cancer history, family history, occupational exposure, recreational exposure, past pulmonary diseases, present pulmonary diseases, prior chest surgery, number of satellite lesions around the lesion to be diagnosed, size of lymph nodes, presence of suspicious nodes, and location of the lesion in the lung.

29. The computer-aided diagnosis system according to claim 21, wherein the classifier determines if a lesion is malignant or benign.

30. The computer-aided diagnosis system according to claim 21, wherein the classifier determines a likelihood of malignancy of a lesion.

Description

(1) The methods and systems provided herein can either give a binary decision, i.e. whether a nodule is malignant or benign, or can provide a probability or a likelihood (e.g. 80%) that a nodule or lesion is malignant. Clinicians use this information as a second opinion to assist in their diagnostic decision making. The methods and systems provided herein can be used for different diseases (e.g. a colon polyp, liver cancer, or breast cancer) and different modalities, such as MRI or CT.

(2) FIG. 1 is a series of drawings showing a box-counting algorithm for finding fractal dimension of a contour.

(3) FIG. 2 is a graph of a computation of fractal dimension, with log (N [boxes]) on the ordinate and log (1/[box size]) on the abscissa, showing a linear relationship between log (N [boxes]) and log (1/[box size]).

(4) FIG. 3, Panel A is a scan of a nodule.

(5) FIG. 3, Panel B is a corresponding segmented contour drawing and the estimated normal angles.

(6) FIG. 4, Panel A is a graph of the intensity gradient of a nodule, with gradient on the ordinate and position along contour on the abscissa.

(7) FIG. 4, Panel B is a graph of the corresponding power spectrum, with power/frequency on the ordinate and normalized frequency on the abscissa.

(8) FIG. 5 is a block diagram of a CADx system.

(9) FIG. 6 shows a border given by gray pixels (left image) and a template (right image) for an example presented herein.

(10) A CADx system that provides high confidence to clinicians improves clinician workflow by providing fast and accurate diagnosis (fewer false positives and false negatives). A CADx system can be used as a second reader to increase clinicians' confidence in their diagnosis, leading to significant reduction of unnecessary biopsies of lung lesions such as nodules. Furthermore, a CADx system can facilitate lung cancer screening of asymptomatic patients since diagnosis can be reached quickly and accurately. MSCT scanners, exemplified but not limited to the Philips Brilliance series, offer increasing resolution and allow finer structures to be observed while producing increasing amounts of image data to be interpreted by radiologists. However, even the latest CADx systems often fail to incorporate clinical information, to use optimized feature extraction, or to apply machine learning techniques.

(11) Features of interest that have not been developed include analysis for both thick and thin slice CT scans. Other proposed features include only simple 3D features or features that are not optimal to describe the differences between benign and malignant nodules. These features often result in a low rate of accuracy and are not desirable to use in a CADx algorithm.

(12) The methods and systems provided herein are based on state-of-the art machine learning techniques, such as genetic algorithms and support vector machines, and innovative image processing algorithms for pre-processing of images and feature extraction. An aspect of the methods and systems provided herein is the capability of combining image-based and clinical information about a patient and a patient's lesion into the decision making process. The methods and systems provided herein combine features extracted from high quality medical images, for example a CT scan, with non-imaging data from patient health records through the use of machine learning and image processing methodology.

(13) The methods and systems (exemplified by FIG. 5) provided herein for computer-aided diagnosis of lesions includes several processing units. A pre-processing unit 100 processes images, e.g. MSCT scans, to create isotropic volumes of interest each surrounding a lesion to be diagnosed, or to segment or delineate lesions, such as pulmonary nodules. A feature extraction unit 106 extracts two dimensional (2D), two and a half dimensional (2.5D) and three dimensional (3D) features from images to characterize lesions. These features, together with clinical information including patient history, constitute the feature pool. A clinical information processing unit 108 accepts and transforms the clinical information to be used in the feature pool. Feature selection, a step used in the design phase of the system, is based on a genetic algorithm, and is used to select an optimal feature subset from the feature pool. A classifier or committee of classifiers 110 are used in the feature selection process, and are built using the selected feature subset to classify a lesion as malignant or benign, or to determine the malignancy of a lesion, such as a pulmonary nodule.

(14) In many instances, because of the complexity of the 3D shape of the nodules, the 2D slices may reveal several disconnected islands. As many of the features described below rely upon a unique delineation of border pixels, it is often desirable to operate on only a single object. The methods and systems in one embodiment use an algorithm to remove all but the largest connected object, and perform all nodule analysis on that object. In this case, analysis is performed only for 2D and 2.5D features.

(15) The methods and systems provided herein facilitate the use of 2D, 2.5D, and 3D features on a multi-slice dataset. Because of the ease of computation of 2D features, it is often desirable to utilize 2D features even when 3D data is available. This is particularly common in thick slice data (i.e. a slice thickness greater than 3 times the in-plane resolution), where 3D features may not be as robust and sufficient. In such cases, the user of the CADx system may manually identify an individual slice for analysis. The CADx system is then said to be operating in a 2D mode, in which the features extracted from this 2D slice are used in performing the classification. Similarly, only 2D features from the training dataset are used in constructing the optimum classifier when the system is operating in 2D mode. Such systems fail to capture the full range of information present in the multi-slice volume.

(16) To overcome this, the methods and systems provided herein use a 2.5D, or pseudo-3D, mode. The same features are used as in the pure 2D mode, and these features are calculated on each slice of the MSCT image. The slice range may encompass the whole volume, or alternatively can be manually selected by the user. The feature values used in the classifier are then taken as the size-weighted average of these slice-by-slice calculations. For some features, it may be more logical to use the maximum or minimum values across slices. Alternatively, the 2D features are computed on maximum intensity projection (MIP) data taken through the range of slices, or from an arbitrary angle, with the 2D feature extraction run on the projected image.

(17) Pre-Processing of Images

(18) Images, such as MSCT scans, are pre-processed to determine regions of interest (ROIs) or volumes of interest (VOIs) for analysis. These serve as an input into the feature extraction unit.

(19) The methods and systems provided herein include several pre-processing steps. ROIs are constructed that allow feature calculation to be robust to segmentation errors. ROIs are constructed by morphological operations of erosion and dilation on the binary nodule image. This is used to construct an internal, external, and boundary region. An additional ROI consisting of the largest square region that can be embedded in the object is also identified. The features described in these methods and systems can be computed using any one or more of these ROIs.

(20) The chest wall and other irrelevant anatomy are excluded from feature calculation. Many current segmentation algorithms also delineate the borders of the chest wall. Voxels or pixels that belong to the chest wall are explicitly excluded from all of the features to be described. Only those pixels that are labelled as belonging to the lung parenchyma or the nodule of interest are included.

(21) For many features, it is desirable to capture the full range of detail within the nodule, including potential cavities within the nodule. Cavitation, or air pockets within the nodule, is often not identified directly by segmentation. Therefore, the methods and systems provided herein use an algorithm for post-processing nodule segmentation to fill in holes or gaps in the segmentation mask, such that air pockets or cavities within the nodule are explicitly considered during feature calculation. No change is performed on the CT images, but rather, on the segmentation results. These methods and systems relate to the use of the post-processed segmentation masks, as well as the original segmentation masks.

(22) As shown in FIG. 5, the pre-processing unit 100 extracts a VOI surrounding a lung nodule based on its location, which can be provided by a clinician or a CAD system.

(23) Since MSCT scans have higher resolution in the slices than among slices, it is desirable to perform interpolation between the slices to create isotropic voxel representation. Interpolation to isotropic voxels is desirable for segmentation purposes. However, it is also preferable to retain the original data for feature calculation, to avoid the filtering properties of interpolation kernels. The methods and systems provided herein use an interpolation method that keeps the original slices and inserts interpolated slices to make the scan near-isotropic. The methods and systems in one embodiment include a means of limiting the interpolation to a single axis in order to reach near-isotropic voxels. These isotropic volumes can then be used for segmentation, and for computation of 3D features. The original 2D slices can later be extracted from this interpolated 3D dataset.

(24) A segmentation step delineates the pulmonary nodules from the background, generating a binary or trinary image (label VOI), where a nodule, background and lung wall regions are labeled. A segmentation algorithm is run on the interpolated data, resulting in a volume of the same dimension as the interpolated volume. See Philips Press Release, Philips' new CT Lung Nodule Assessment and Comparison Option can enable clinicians to identify and treat lung cancer, 2003; Wiemker et al., Options to improve the performance of the computer aided detection of lung nodules in thin-slice CT, Philips Research Laboratories: Hamburg, p. 1-19, 2003; and Wiemker et al., Computer Aided Tumor Volumetry in CT Data, Invention disclosure, Philips Research Hamburg, 2002. For analysis of 2D features, the slices in the label volume that correspond to the original slices are identified and extracted. For 3D features, the full near-isotropic 3D volume is used.

(25) In one embodiment of the invention, manual selection of segmentation masks from segmentation using varying threshold is used. Thus, various segmentation thresholds and seed placements are tested. Grayscale images overlaid with segmentation contours are presented to a user, who then manually selects the optimum segmentation, a single slice for 2D feature extraction, and a range of slices for 2.5-D feature extraction.

(26) Feature Extraction Unit

(27) Feature extraction is performed to extract 2D, 2.5D, and 3D features images to characterize lesions, such as pulmonary nodules. These features, together with clinical information, constitute a feature pool.

(28) Using a gray-level and labeled VOI, the feature extraction unit 106 calculates different features, such as 2D (using the native slice of the VOI with the largest nodule cross-section), 2.5D (an average of 2D features calculated on all the native slices, weighted by the nodule cross-sectional area), and 3D (based on the near-isotropic VOI) features. The feature extraction step is significant, since the calculated features need to have enough discriminatory power together with the clinical information to distinguish between malignant and benign nodules. Features can be, for example, gray level distributions inside and surrounding the nodules, shape information, texture information inside and outside the nodule, gradient information on the surface of the nodule, or contrast between the inside and outside of the nodule. Each of these features can be calculated in 2D, 2.5D, or 3D.

(29) Clinical Information and its Transformation to Clinical Features

(30) Since clinical information is important in the diagnosis process, the methods and systems provided herein include a unit 108 which converts clinical information into a suitable form so it can be combined with extracted image based features for the feature selection process. For example, clinical information for gender is divided into 2 categories, such as whether the patient is male or whether is the patient female. The clinical information which can be used in the proposed system can include, for example, age, gender, smoking history, cancer history, family history, occupational exposure, recreational exposure, past or present pulmonary diseases (e.g. emphysema), number of satellite nodules around the one to be diagnosed, size of lymph nodes, presence of other suspicious nodules, or location of the nodule in the lung (e.g. upper lobe or lower lobe).

(31) Feature Selection

(32) The feature selection unit 112 finds the most relevant features from the feature pool containing image-based as well clinical features. A GA and SVM-based feature selection process is used. Once the most relevant features are determined based on a training dataset, a classifier or committee of classifiers is built based on the optimal feature subsets and the feature selection unit 112 is no longer required. In the embodiment relating to the committee of classifiers, each committee member can be constructed on a feature subset identified through a separate run of the feature selection algorithm. The diversity of the classifiers in the committee is achieved for example for the GA-based feature selection by giving a different randomly selected set of training and testing data to each GA run. Other feature selection methods, such as statistical difference filters, correlation filters, stepwise linear regression, recursive feature elimination (RFE), and a random feature selection can also be used.

(33) Classifier

(34) Following the supervised learning principle, a classifier is built using the selected optimal feature subset and training data. Possible classifiers are SVMs, decision trees, linear discriminant analysis and neural networks. SVMs are often used since they have shown superior performance with respect to classifiers.

(35) Committees of classifiers can also be used. In this case, several classifiers are built on different feature subsets, which were selected as the best feature subsets via the feature selection process. Each classifier will result in a decision, e.g. whether a nodule is malignant or benign. A majority vote of the classifiers will determine for example, the output of the CADx system for the nodule in question. In likelihood CADx, the output likelihood of the committee may be the fraction of positive votes, or an average (which may be weighted) of the individual likelihood results from each member of the committee.

EXAMPLES

(36) The methods and systems provided herein apply several image processing methods to CADx. Examples of these image processing methods are provided below.

Example 1

Invariant Moments

(37) Moments are a means of mathematically representing the spatial distribution of an object. This includes the shape (binary moments) or density distribution (grayscale moments). Invariant moments are those moments that do not change when the object undergoes some transformation, such as rotation, scaling, and translation. Moment based methods have been used at length in computer vision and optical character recognition. The use of invariant moments, including the first 6 invariant moments, is described by the mathematical formalism of Hu (1962).

(38) For an image I of size NM, the moments m.sub.pq can be given by:

(39) $m_{\mathrm{pq}}=\underset{x=1}{\overset{N}{.Math.}}\underset{y=1}{\overset{M}{.Math.}}{x}^p{y}^{q}I\left(x,y\right)$
where I takes on the grayscale intensity of the image at a pixel (x, y) for grayscale moments, or a value of 0 or 1 for binary moments. The centroid of the object is defined using the following moments.
x=m.sub.10/m.sub.00 and y=m.sub.01/m.sub.00

(40) These allow calculation of centralised moments .sub.pq that are invariant under translation:

(41) ${}_{\mathrm{pq}}=\underset{x=1}{\overset{N}{.Math.}}\underset{y=1}{\overset{M}{.Math.}}{\left(x-\stackrel{\_}{x}\right)}^p{\left(y-\stackrel{\_}{y}\right)}^{q}I\left(x,y\right)$
which can be made scale invariant by computing normalised central moments
.sub.pq=.sub.pq/.sub.00.sup.
where =(p+q)/2+1 for all (p+q)2.
These scale and translation invariant moments can be converted into rotationally invariant moments, using a method described by Hu (IRE Trans Information Theory, IT-8 179-197, 1962). These seven moment invariants are given as the following equations.
H.sub.1=.sub.20+.sub.02
H.sub.2=(.sub.20.sub.02).sup.2+4.sub.11.sup.2
H.sub.3=(.sub.303.sub.12).sup.2+(3.sub.21.sub.03).sup.2
H.sub.4=(.sub.30+.sub.12).sup.2+(.sub.21+.sub.03).sup.2
H.sub.5=(.sub.303.sub.12)(.sub.30+.sub.12)+((.sub.30+.sub.12).sup.23(.sub.21.sub.03).sup.2)+(3.sub.21.sub.03)(.sub.21+.sub.03)(3(.sub.30+.sub.12).sup.23(.sub.21+.sub.03).sup.2)
H.sub.6=(.sub.20.sub.02)((.sub.30+.sub.12).sup.2(.sub.21+.sub.03).sup.2)+4.sub.11(.sub.30+.sub.12)(.sub.21+.sub.03)
H.sub.7=(3.sub.21.sub.03)(.sub.30+.sub.12)((.sub.30+.sub.12).sup.23(.sub.21+.sub.03).sup.2)+(3.sub.12.sub.30)(.sub.21+.sub.03)(3(.sub.12+.sub.30).sup.2(.sub.21+.sub.03).sup.2)

(42) The moments and central moments are computed by iterating over different p and q values and over all pixels in an image, skipping over pixels that are not in the nodule. The image function, I, can be binary or real-valued, and in practice, both are implemented. In this way, each nodule can contribute seven invariant binary moments, and seven invariant grayscale moments. These 14 scalar values can then be used as inputs into a CADx system.

(43) Alternatively, to make the feature extraction robust to segmentation uncertainty, the grayscale moments are calculated by performing the moment calculation on a circular area enclosing the segmented nodule. The pleural wall and irrelevant lung structures are identified via the segmentation and removed. To avoid incorrectly placing 0 values in these removed structures, random noise is inserted into these pixels by sampling from a histogram computed from the retained background.

(44) Further, 3D moments are extracted, and features are derived from the 3D moments in binary mask gray scale data as shown by the following equations.
moment.J11=.sub.200+.sub.020+.sub.002
moment.J21=.sub.200.sub.020+.sub.200.sub.002+.sub.020.sub.002.sub.101.sup.2.sub.110.sub.110.sub.011.sup.2
moment.J31=.sub.200.sub.020.sub.002.sub.002.sub.110.sup.2+2*.sub.110*.sub.101*.sub.011.sub.020*.sub.101.sup.2.sub.200*.sub.011.sup.2

(45) The three derived features above are also computed for moments derived for the grayscale nodule image.

Example 2

Invariant Fourier Descriptors

(46) Fourier descriptors are a means of mathematically describing the shape of an object. Conceptually, the border of an object is described in terms of frequencies. Low frequencies describe smooth edges and give the general shape of the object, whereas high frequencies describe irregularities and sharp changes in the contour. The Fourier descriptors describe the particular mix of frequencies that constitute a particular object shape. In the methods and systems provided herein, Fourier descriptors that are invariant to scale, translation, and rotation are used. For use in the classifier, two means of condensing the Fourier coefficients into a scalar value are provided.

(47) Initially, the pixels on the edge of a nodule are identified by a number of well-known methods. Care is taken to ensure that these N edge pixels are listed in clockwise or counter-clockwise order around the object. Each edge pixel is described with its x and y coordinate, yielding the vectors x={x.sub.1, x.sub.2, . . . , x.sub.N} and y={y.sub.1, y.sub.2, . . . , y.sub.N}. The discrete Fourier transform for each object is computed using any of a number of techniques (denoted generically as FT), yielding
v=FT{x} and w=FT{y}
which can be decomposed into real and imaginary parts
v.sub.n=a.sub.n+ib.sub.n and w.sub.n=c.sub.n+id.sub.n.
The Fourier descriptors for the shape uniquely defined by x and y is then expressed as:

(48) $f_n=\frac{\sqrt{a_n^2+b_n^2}}{\sqrt{a_1^2+b_1^2}}+\frac{\sqrt{c_n^2+d_n^2}}{\sqrt{c_1^2+d_1^2}}.$
This yields a vector f which may vary in length between nodules, thus rendering it difficult to use in comparing nodules in a CADx system. In the methods and systems provided herein, vector f is condensed into two scalar descriptors g.sub.1 and g.sub.2:

(49) $g_1=\left(\underset{n=2}{\overset{N}{.Math.}}\frac{f_n}{n}\right)/\underset{n=2}{\overset{N}{.Math.}}n\mathrm{and}{g}_2=\left(\underset{n=2}{\overset{N}{.Math.}}{\mathrm{nf}}_n\right)/\underset{n=2}{\overset{N}{.Math.}}n.$

(50) The mathematical description given here follows that of Nixon et al., Feature Extraction & Image Processing, Butterworth-Heinemann: Woburn, Mass. pp. 269-278 (2002). This calculation is used as an input into a CADx system.

Example 3

Fractal Dimension of Nodule Shape

(51) The contour of an object can be described using the Minkowski-Bouligand fractal dimension, also known as the box-counting dimension. This describes how the measured length of the contour changes with the resolution at which the contour is viewed. A nodule with a high fractal dimension tends to exhibit irregularities in the surface shape, which may indicate malignancy. The fractal dimension can be calculated for 2D, 2.5D, or 3D features, and in fact, all three may be simultaneously used during classification.

(52) A nodule with an outline, shown in the first panel of FIG. 1 (Scale 1), is considered. This outline consists of all pixels that are nearest to the exterior of the nodule. The total number of pixels in the outline is tallied as N (1). The image is then resampled by a factor of , shown in FIG. 1 as Scale 2. This can be thought of as tiling 22 pixel boxes over the Scale 1 image and seeing how many of these boxes contain and edge pixels. This count of boxes is denoted N (). This is then repeated with 33 pixel boxes to yield N (), and so on. For a fractal object, the value N (1/d) varies with the scale d according to.
N(1/d)=(1/d).sup.FD
where FD is the fractal dimension of the object. By algebraic manipulations, this is changed to:
ln N(1/d)=FD ln(1/d)+ln .
Thus, if the value of N.sub.d is computed for several d, then the FD can be estimated by a linear fit between ln (1/d) on the x-axis and ln N (1/d) on the y-axis. This fit may be a least squares fit or a robust fit, or may be taken by taking the average slope between successive points. An example of such a fit is given in FIG. 2. The slope of this line is used as a feature in a CADx system.

Example 4

Fractal Dimension of Nodule Texture

(53) The texture of an object can be described in terms of fractal dimension. A 2D view of the object is considered as a surface embedded in 3D. The fractal dimension measures the variations in grayscale intensity as the resolution of the image is changed. Higher fractal dimensions suggest greater complexity in the internal structure of the object. As with the border, the fractal dimension of the texture can be calculated for 2D, 2.5D, or 3D features, and in fact, all three may be simultaneously used during classification.

(54) The texture dimension is calculated in a manner similar to that of the border, described above in Example 3, with the following changes: the dd boxes are not tiled, but overlap; that is, all possible dd boxes are considered. Then, instead of counting the number of boxes, N (1/d) represents the differences between the maximum and minimum pixel gray value within each dd box, summed over all boxes of size dd.

Example 5

Edge Gradient-Based Features

(55) Border pixels of a nodule are identified, and the normal direction pointing away from an object is computed. The derivative is then estimated by computing the difference in grayscale intensity between the border pixel and the pixel some finite distance away in the normal direction. This second external value is typically found by interpolating the value of the image two pixels away. The derivative as a function of position along the contour can then be analyzed statistically to yield scalar features using the mean, the standard deviation, and the root-mean-square variation and first moment of the power spectrum of this function.

(56) As with the Fourier descriptors, each edge pixel and its x and y coordinate are considered, yielding the vectors x={x.sub.1, x.sub.2, . . . , x.sub.N} and y={y.sub.1, y.sub.2, . . . , y.sub.N}. For an edge pixel i located at (x.sub.i, y.sub.i), the normal angle is computed. This is given by:

(57) $=\frac{1}{2}\left({\tan }^{-1}\left(\frac{y_{i+k}-y_i}{x_{i+k}-x_i}\right)+{\tan }^{-1}\left(\frac{y_i-y_{i-k}}{x_i-x_{i-k}}\right)\right)+{90}^{},$
taking care to ensure the consistency of the signs. That is, the normal is computed as being perpendicular to the local curve of the surface, where the local curvature is estimated by averaging between two nearby (but not necessarily adjacent) points. An example of these angles as computed on a nodule is given in FIG. 3.

(58) A point d pixels away from the edge is then defined by the coordinates
x=x.sub.i+d sin and y=y.sub.i+d cos ,
with an intensity at that location found by bilinear interpolation on the image. The edge gradient is then simply the difference between the interpolated intensity and the intensity of the original edge pixel. This can be computed for every edge pixel in order, as shown in FIG. 4. A large number of statistical features can be calculated from this, including mean and standard deviation. The power spectrum can estimated through many well known techniques, thus yielding information about the frequency content of the fluctuations around the object. The underlying assumption to this analysis is that high frequency variations in intensity may suggest spiculations, which are believed to be an indicator of malignancy. To assess this, the moments and root mean square variation of the power spectrum can be calculated and used as features.

Example 6

Detection and Characterization of Internal Clusters

(59) The presence of dark or bright clusters in the object may be indicative of calcification or cavitation, respectively. The methods and systems provided herein include means of identifying these regions automatically, and using them as features for CADx.

(60) A threshold is iteratively applied to an image, such that the only pixels that remain are above some intensity t. The value of t is lowered, starting at the maximum value in the object to the minimum value of the object. At each threshold, the remaining pixels are grouped into clusters that are connected on the sides. That is, each cluster consists of a set of pixels where it is possible to move from one pixel to any other pixel by travelling in one of the four cardinal directions, without ever leaving the cluster. If the largest cluster has a size of greater than n pixels, then the value of the threshold t is saved and the algorithm halted. The feature that is extracted is the given by either the threshold value expressed in units of image pixel intensity, or as the number of standard deviations between the mean intensity of the rest of the detected cluster and the critical threshold. Similarly, by thresholding below the iterated intensity threshold, it is possible to detect dark clusters; the same two features can be calculated for these clusters.

Example 7

Chain Code

(61) The pixels on the border of an object image can be identified in a continuous order, such that each pixel has a pair of neighbors. By tracking the relative positions of neighboring pixels, a vector describing the shape of an object is identified.

(62) With reference to FIG. 6, for each pixel on the border, the two neighbors of this pixel are defined using a chain code. Consider a border given by the shaded pixels of the image on the left of FIG. 6. By placing the centre of the template on the right of FIG. 6 on each pixel of the boundary, the chain code description is read out {1, 2, 1, 2, 4, 4, . . . } with a difference between successive values of {1, 1, 1, 2, 0, . . . } or an absolute difference of {1, 1, 1, 2, 0, . . .}.

(63) The distribution of these chain code and chain code difference values can be used to calculate features. The fraction of absolute difference values greater than 1, 2, . . . , 6 can be used as a feature. This fraction can be used to detect the number of sudden changes in direction, thus describing the irregularity of nodules, or concepts such as spiculation and lobulation.

Example 8

Texture: Neighbourhood Gray-Tone Difference Matrix

(64) Several mathematical methods exist for approximating the human perception of texture of an image, including the aforementioned fractal method. An alternative method is based on what is known as the neighborhood gray-tone difference matrix (NGTDM). This method seeks to quantitatively describe the differences between each pixel and its surrounding neighborhood, leading to mathematical descriptions that have been shown in psychometric tests to correlate well with subjective ratings of the abstract qualities including coarseness, contrast, busyness, complexity, and strength.

(65) The description given by Amadsun and King (IEEE Trans Sys Man Cybernetics 19(5): 1264-1274, 1989) is followed. The NGTDM is a matrix formed as follows. The NN image as a whole is quantized to a predetermined number of levels, g. For every pixel in the particular ROI with quantized intensity denoted by f(k,l),

(66) $A\left(k,l\right)=\frac{1}{{\left(2d+1\right)}^2-1}\left[\underset{m=-d}{\overset{d}{.Math.}}\underset{n=-d}{\overset{d}{.Math.}}f\left(k+m,l+n\right)\right]$$\mathrm{where}\left(m,n\right)\left(0,0\right)$
That is, A is the average of a (2d+1)(2d+1) neighborhood around the pixel of interest, excluding the pixel itself. The NGTDM matrix N has one column and as many rows as their are levels of intensity in the image. The row i is then given by

(67) $N\left(i\right)=\underset{k,l}{.Math.}.Math.i-A\left(k,l\right).Math.$
where the summation is taken over all pixels where f(k,l)=i. Care is taken to restrict the calculation to the area within the region of interest (ROI), whether this is a nodule itself or the region outside the nodule.

(68) The probability for each bin, p(i), is defined as the fraction of centre pixels that contributed to the calculation of N(i). Two features are shown here; others are described in detail in Amadsun et al., IEEE Trans Sys Man Cybernetics, 19(5): 1264-1274, 1989. These features are:

(69) $\mathrm{coarseness}={\left[\underset{i}{.Math.}p\left(i\right)N\left(i\right)\right]}^{-1}$

(70) $\mathrm{contrast}=\left[\frac{1}{G\left(G-1\right)}\underset{i=0}{\overset{g}{.Math.}}\underset{j=0}{\overset{g}{.Math.}}p\left(i\right)p\left(j\right){\left(i-j\right)}^2\right]\left[\frac{1}{n^2}\underset{i=0}{\overset{g}{.Math.}}N\left(i\right)\right]$
where n=(N2d) and G is the actual number of the g levels that appear in the image.

(71) The methods and systems provided herein include pre-processing CT images for CADx, several new features that can be used as inputs into a CADx system and advanced machine learning techniques for feature selection and classification. The methods and systems provided herein overcome difficulties in dealing with thick-slice CT volumes, provide robustness to errors in identifying the border of lung nodules, and are optimized to improve diagnostic accuracy of classification systems. The methods and systems provided herein use pre-processing and features to solve this problem by providing 2D, 2.5D, and 3D characterization of lesions, such as benign and malignant lung nodules. Therefore, a machine learning system using uses these features can distinguish more accurately benign nodules from malignant nodules and can achieve higher specificity and sensitivity than the systems without them. Furthermore, the proposed system has the capability of combining image-based and clinical information about the patient into the decision making process.

(72) The CADx methods and systems provided herein can be used with several modalities, for example MRI and CT. The methods and systems provided herein can be used in radiology workstations (exemplified but not limited to Philips Extended Brilliance Workstation, Philips Mx8000, and Philips Brilliance CT scanners) or incorporated into PACS systems (e.g. Stentor iSite). The CADx methods and systems provided herein can be used in diagnosing different diseases, including but not limited to colon polyps, liver cancer, and breast cancer.

(73) It will furthermore be apparent that other and further forms of the invention, and embodiments other than the specific and exemplary embodiments described above, may be devised without departing from the spirit and scope of the appended claims and their equivalents, and therefore it is intended that the scope of this invention encompasses these equivalents and that the description and claims are intended to be exemplary and should not be construed as further limiting. The contents of all references cited herein are incorporated by reference.