Abstract:
A three dimensional mask of the lungs can be automatically created by thresholding, labeling connected components, selecting the dominant object, and alternately employing dilation and erosion operations. With this mask the lungs can be separated from the other anatomic structures on volumetric CT images. Local density maxima in the lungs are then determined by sequentially decreasing thresholds. As the threshold declines, more and more objects (a 3D object is a group of connected voxels with density values larger than the threshold) become apparent. Geometrically overlapped objects at the subsequent threshold levels are either merged into one object or identified as local density maximum (maxima) and plateau. This process terminates if the threshold reaches a predefined density value. Other information about small lung nodules such as compact shape and size are combined into the algorithm to further remove those detected local density maxima that are not likely to be nodules.

Description:
This application claims the benefit of U.S. Provisional Application No. 60/338,880 filed Oct. 24, 2001, which is incorporated herein by reference. 
    
    
     BACKGROUND OF THE INVENTION 
     The present invention is related to computer-aided analysis of lung nodules in medical images. In particular, the invention is directed to computer-aided review of computed tomography images for initially detecting lung nodule candidates for a subsequent analysis based upon the features of the lung nodule candidates. 
     Lung cancer is the leading cause of cancer deaths among the population in-the United States. Each year there are about 170,000 newly diagnosed cases of lung cancer and over 150,000 deaths. More people die of lung cancer than of colon, breast, and prostate cancers combined. Despite the research and improvements in medical treatments related to surgery, radiation therapy, and chemotherapy, currently the overall survival rate of all lung cancer patients is only about 14 percent. Unfortunately the survival rate has remained essentially the same over the past three decades. The high mortality rate of lung cancer is caused by the fact that more than 80% lung cancer is diagnosed after it has metastasized. Patients with early detection of lung cancer followed by proper treatment with surgery and/or combined with radiation and chemotherapy can improve their five-year survival rate from 13 percent to about 41 percent. Given that earlier-stage intervention leads to substantially higher rates of survival, it is therefore a major public health directive to reduce the mortality of lung cancer through detection and intervention of the cancer at earlier and more curable stages. 
     The development of the computed tomography (CT) technology and post-processing algorithms has provided radiologists with a useful tool for diagnosing lung cancers at earlier stages. However, current CT systems have their inherent shortcomings in that the amount of chest CT images (data) that is generated from a single CT examination, which can range from 30 to over 300 slices depending on image resolution along the scan axial direction, becomes a huge hurdle for the radiologists to interpret. Even though small lung nodules can be captured by helical CT and the images can be viewed in either a traditional film-based mode or cine mode on today&#39;s Picture Archiving and Communication System (PACS) workstations, the potential of overlooking small nodules in the diagnostic process has become major concerns. 
     Accordingly there is a need for a system that automatically identifies small lung nodule candidates from helical CT images to assist radiologists in improving the detection of nodules in the clinical practice. 
     SUMMARY OF THE INVENTION 
     The present invention is a method and apparatus for analyzing volumetric chest computed tomography images for lung nodules. The steps of the method of the invention include initially obtaining the volumetric chest computed tomography images for lung nodules from an image acquisition device. The lungs are then separated from the other anatomic structures on the images to form lung images. The lung images are then processed to detect nodule candidates with a local density maximum algorithm. False-positives among the detected nodule candidates are then reduced by an application of parameters concerning lung nodules. Preferably the parameters concern at least the size and shape of the nodule candidates. 
     The present invention similarly includes an article of manufacture for analyzing volumetric chest computed tomography images for lung nodules. The article includes a machine readable medium containing one or more programs which when executed implement the method of the invention. 
     For a better understanding of the present invention, reference is made to the following description to be taken in conjunction with the accompanying drawings and its scope will be pointed out in the appended claims. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     Preferred embodiments of the invention have been chosen for purposes of illustration and description and are shown in the accompanying drawings, wherein: 
     FIG. 1 is a flowchart illustrating the process of examining a patient using computed tomography (CT) technology in accordance with the present invention; 
     FIG. 2 is an illustration of a slice of a typical chest CT image series obtained with General Electric LightSpeed QX/i; 
     FIG. 3 is an illustration showing the volumetric density histogram of the image series of which one slice is shown in FIG. 2; 
     FIG. 4 is an illustration showing the original slice shown in FIG. 2 thresholded at t=750 gray level; 
     FIG. 5 is an illustration showing an incomplete lung mask for the original slice shown in FIG. 2; 
     FIG. 6 is an illustration showing the complete lung mask for the original slice shown in FIG. 2; 
     FIG. 7 is an illustration showing the extracted lungs from the original slice shown in FIG. 2; 
     FIG. 8 is an illustration showing a one dimensional example of the detection of local density maxima; 
     FIG. 9 is a flowchart of the main loop of the local density maximum algorithm for detecting lung nodules; 
     FIGS. 10A and 10B is a flowchart illustrating details on finding objects at a threshold level; 
     FIG. 11 is a flowchart illustrating details for determining the local maxima at the current threshold value based on the analysis of the relationship between previous and current objects and the criteria of being a local maximum on finding objects at a threshold level; 
     FIG. 12 is an illustration showing the three dimensional relationship of the binary image data; 
     FIG. 13 is an illustration showing an extracted lung image; 
     FIG. 14 is an illustration showing the extracted lung image shown in FIG. 13 thresholded at a level of 145; 
     FIG. 15 is an illustration showing the extracted lung image shown in FIG. 13 thresholded at a level of 90; 
     FIG. 16 is an illustration showing the extracted lung image shown in FIG. 13 thresholded at a level of 70; and 
     FIG. 17 illustrates modification of an object&#39;s bounding box when applying parameters concerning lung nodules. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Referring to FIG. 1, a flowchart illustrating the process of examining a patient using the computed tomography (CT) technology in accordance with the present invention is depicted. The initial step  102  is to obtain the volumetric chest CT images of the patient. Preferably the images are acquired helically using a standard clinical protocol and the image data is transferred to and stored on a standard local computer for further processing. Suitable scanners for obtaining the images include General Electric HighSpeed CT/i and LightSpeed QX/i. The image data can be transferred via an Ethernet and the computer can be a personal computer or a UNIX machine. The next step  104  is to separate the lungs from the other anatomic structures on the volumetric chest CT images. The separation of the lungs  104  is explained in detail below with reference to FIGS. 2 through 7. Similar techniques to extract the lungs from CT chest images can be found elsewhere (Shiying Hu, Eric A. Hoffman, and Joseph M. Reinhardt, “Automatic Lung Segmentation for Accurate Quantitation of Volumetric X-Ray CT images”. IEEE Transactions on Medical Imaging, Vol. 20, No. 6, June 2001; Li Fan, Carol L. Novak, Jiangzhong Qian, Gerhard Kohl, and David P. Naidich, “Automatic Detection of Lung Nodules from Multi-Slice Low Dose CT Images”. Proc SPIE 2001; 4322:1828-1835.). After the lungs are separated, the next step  106  is to detect nodule candidates using a local density maximum algorithm in accordance with the present invention. The final step  108  is to reduce false-positives among the nodule candidates by analyzing the characteristics of the detected nodule candidates. Features that can distinguish spherical shape (nodules) from cylinder shape (vessels and bronchi) such as three-dimensional compactness factor (Binsheng Zhao, Anthony Reeves, David Yankelevitz, and Claudia Henschke, Three-dimensional multi-criterion automatic segmentation of pulmonary nodules of helical CT images. Optical Engineering 1999; 38(8):1340-1347, which is incorporated herein by reference) can be used to delete false positive results. 
     Referring initially to FIG. 2, a slice of a typical chest CT image series obtained with General Electric LightSpeed QX/i is shown while the corresponding volumetric density histogram is shown in FIG.  3 . The histogram illustrates that there are typically four peaks around the mean density values representing, from left to right, the background outside the body but within the scan area, lung parenchyma, fat, and muscles. The histogram also includes a long, flat and low valley between the peaks of the lung parenchyma and the fat. The presence of the valley within the histogram indicates that the separation of the parenchyma from the soft tissues and bones is generally insensitive to a threshold set within the valley. This characteristic is useful in extracting or separating the lungs from the other structures on the volumetric chest CT images in accomplishing step  104  as shown in FIG.  1 . 
     The extraction of the lungs from the other anatomic structures generally involves initially selecting a threshold value lying in the valley close to the peak of the fat to initially separate lung parenchyma and the background from the soft tissues and bones as shown in FIG.  4 . In FIG. 4, a threshold value of 750 gray level (All values of gray level obtained from the GE Lightspeed QX/i machine equate to Hounsfield units HU as follows: gray level=HU+1024, for example 1000 gray level=−24 HU.) was selected to extract the soft tissues and bones. Voxels having density lower than the threshold value will be recognized as lung candidates and assigned the value of one and appear white in FIG. 4, whereas other voxels assigned the value of zero in the resultant image appear black in FIG.  4 . In comparing FIG. 4 with FIGS. 2 and 3, the soft tissues and bones which have voxel densities greater than the threshold value appear black while the background and lung parenchyma which have voxel densities lower than the threshold value appear white. 
     The lung parenchyma is next separated from the background by initialing labeling connected voxel components. The voxel components are labeled by finding geometrically connected voxels that have the value of one in the thresholded binary images and assigning an identical number to each individual group of the connected voxels. The largest component that does not touch any margin of the images is selected next as shown in FIG.  5 . Since the apparent densities of vessels and bronchi in the lungs vary widely, regions in the lungs with higher densities generally get grouped into soft tissues and bones causing an incomplete extraction of the lung mask as shown in FIG.  5 . To obtain a complete lung mask (i.e., to include those high-density segments of vessels and bronchi as well as nodules into the parenchyma group) a morphological dilation followed by an erosion of the result is applied as shown in FIG.  6 . The lung mask is used to extract the lung regions as shown in FIG. 7 from the original chest CT images shown in FIG.  2 . After the extraction, the resulting images are referred to lung images. The lung images are preferably converted from the image type of short (16 bits) to the type of byte (8 bits) by a linear projection for further processing. 
     Referring initially to FIG. 8, the detection of lung nodules using a local density maximum algorithm is illustrated in a one-dimensional example. The illustration in FIG. 8 represents a plot of the voxel density along the length of a cut through one plane of volumetric data which would appear as a single image similar to that which is shown in FIG.  2 . In FIG. 8, six (6) local maxima are represented by the individual peaks in the plot which signify the detection of a nodule candidate for each peak. The local maxima are detected in the volumetric data in a dynamic thresholding manner which is described below with reference to the flowcharts shown in FIGS. 9,  10 A,  10 B, and  11 , and FIGS. 13 through 16. 
     Referring now to FIG. 9, the main loop of the algorithm for the detecting nodule candidate step  106  using a local density maximum algorithm is illustrated. The lung images are segmented using an initial threshold determined at step  202 . At step  202 , the initial threshold level is set to a value equal to the maximum voxel density value of the lung images minus a predefined threshold step value. The step value is illustrated in FIG.  8  and is preferably selected within a range between about 5 and 10 (in gray level). The step value selected for generating the images included in this application was 7 (in gray level). At step  204 , the current threshold value is compared against a predefined stop value to determine whether further processing is required. The stop value can be set as low as 0, but is preferably set at about 15 (in gray level). If the threshold value is greater than the predefined stop value, the processing proceeds to step  206  to find 3D objects at this threshold level. Otherwise, the results of the algorithm will be output for further processing at step  108  as shown in FIG.  1 . When step  206  is processed, the binary image that results from the current threshold value is processed to identify three-dimensional (3D) objects that are represented by geometrically connected voxels having a value of one. The details of step  206  are explained in more detail below with reference to FIGS. 10A and 10B. At step  208 , it is determined whether a previous object had been detected. During the first loop when the threshold value is equal to the initial threshold value set in step  202 , there is never a previous object. If no previous object is detected, the current object or objects are saved to a previous object array (step  212 ) and the current threshold value is reduced by the step value at step  214  and the processing loops back to step  204  for repeating steps  204  through  214  as long as the criteria at step  204  is satisfied as discussed above. If there are previous objects at step  208 , the processing proceeds to step  210  for a determination of local maxima at the current threshold value based on the analysis of the relationship between previous and current objects and the criteria of being a local maximum. The details of step  210  are explained below with reference to FIG.  11 . The processing continues with step  212  as described above. 
     Referring now to FIGS. 13 through 16, the concept of the dynamic thresholding of the volumetric data is illustrated in a two dimensional view of extracted lung images that are subsequently thresholded at different threshold values to form binary images. In particular, FIG. 14 shows the original extracted lung shown in FIG. 13 being thresholded at a voxel density level of  145  (in gray level). At the threshold level of  145 , at least all voxel data having a density level greater than  145  is assigned the value  1  and appear as white spots in FIG.  14 . As the threshold value is reduced in FIGS. 15 and 16, more objects are identified. The actual identification of current and previous objects is a three dimensional analysis which is detailed in the flowchart in FIGS. 10A and 10B and is illustrated in FIG.  12 . 
     Referring now to FIGS. 10A,  10 B, and  12 , in step  206  the binary image that results from the current threshold value is processed to identify three-dimensional (3D) objects that are represented by geometrically connected voxels having a value of one. In step  206 , the volumetric binary data is analyzed for segments which are connected pixels having a value of 1 in a row and assigning an object number to each segment. In FIG. 12, object “p”, which is shown encircled, illustrates a single object that is present on planes k−1, k, and k+1. Where two segments overlap in three dimensions, the object number is unified. 
     Referring now to FIGS. 10A and 10B, the loop in step  206  starts from the first plane (k=0) and the second row (j=1)  304  on the first plane. Next at step  306 , if it is beyond the last row on the first plane, the processing of the next plane (k=k+1) is initiated again at the second row (j=1)  308  on the new plane. If the processing has not completed the last row on the first plane, an analysis  310  is made to check Whether there exists any overlap between segments in row j and row j−1 on the plane k and segments in row j on the plane k+1. The object number is unified when there is an overlap  312 . This process repeats  314  until all of the rows in the current plan are processed. Once the processing is completed for a plane, the next plane is similarly processed starting from the first row on the new plane. In summary, each time in the processing of step  206  the binary volumetric data is iteratively processed to find three dimensional objects. 
     Referring now to FIG. 11, in step  210  a determination of local maxima at the current threshold value is performed based on the analysis of the relationship between previous and current objects and the criteria of being a local maximum. Since the objects are composed of segments, the analysis is actually performed on the basis of segments. For each current object an initial determination  404  is made as to whether it has more than 1 previous object on it. If the object has more than one previous object on it, a determination  410  as to whether local maxima have already been found on this object is made. When local maxima have already been found on this object, the current object is marked as a plateau  416 . Otherwise when there have not been any local maxima found on this object, an analysis  412  is conducted to determine whether any of the previous objects meet the criteria of being a local maximum. If the criteria of being a local maximum is met, the previous object is marked as a local maximum  414 . The criteria for being a local maximum preferably include size and density. The size requirements are preferably defined so that the object must be larger than a range from about 15 to 25 voxels, and most preferably about 20 voxels to be considered a local maximum. The density is preferably defined so that the object must be significantly higher than the current threshold level to be considered a local maximum. Generally, the difference between the two densities must be greater than 12 (in gray level). If a local maximum is determined, the object containing these local maxima will be interpreted as a plateau  416 . 
     If the object at step  404  does not have more than 1 previous object on it, then a determination  406  is made as to whether there is only one previous object on the current object. When there is no previous object on the current object a new object is detected and is registered as an object detected at this threshold value. When there is only 1 previous object on the current object, the object-to-box ratios (defined as the ratio of the size of the object to the size of the rectangular box that contains the object) of current and previous objects are checked  408  to determine whether the ratio of the object-to-box ratio of the previous object to that of the current object is larger than a specified value. The specified value can range between about 25 and 35 and is preferably selected to be 30. The current object will replace the previous one provided the specified value is not exceeded. Otherwise, the previous object will be recognized as a local maximum and the current object will be recognized as a plateau. 
     Referring now to FIG. 1, after step  106  is completed to detect the local density maxima in the lungs, false-positives are reduced in step  108  from the nodule candidates by applying knowledge known about the lung nodules. Only the local maxima falling into the interesting range of nodule size (about 2 mm&lt;diameter&lt;about 20 mm) will be reserved. Although a nodule on a CT slice can appear to be similar to a vessel that is perpendicular to the slice, nodules and vessels have very discriminating shapes in three dimensions; the former have a spherical shape, whereas the latter have a tubular shape. The following ratios are calculated to determine the “three-dimensional compactness factor” of an object. The ratios include at least the following: 
     1) the ratio of the size of the object to the size of the rectangular box that contains the object (always smaller than 1.0); and 
     2) the ratio of the object lengths along x, y and z axes. 
     If the first ratio is larger than about 0.4, the object is considered to be compact. If the second ratio of max(object length along x axis, object length along y axis) to min(object length along x axis, object length along y axis)&lt;about 1.5 and the ratio of object length along z axis to min(object length along x, object length along y axis)&lt;about 2, the object is considered to be compact. The latter ratio is set slightly larger than the former ratio because of the partial voxel effect caused by the non-isotropic CT images in three directions. By calculating the object length along z axis, the non-isotropic character of CT images need to be taken into consideration. 
     In a preferred embodiment, only the local maxima falling into the interesting range of nodule size (about 2 mm&lt;diameter&lt;about 20 mm) is reserved, and the following ratios are calculated to determine the shape compactness of an object. The ratios include at least the following: 
     1) R 1  is defined as the ratio of the volume of the object to the volume of a modified bounding box of the object, i.e., 
     
       
           R   1 =number of object voxels/(max( dx, dy )*max( dx, dy )* dz )  (1) 
       
     
     where, dx, dy and dz are the maximal projection lengths (in pixel) of the object along the axes of x, y and z, respectively, and max(dx, dy)=dx, if dx&gt;=dy; otherwise, max(dx, dy)=dy. The volume of the box containing an object is given by the denominator in Eq. (1) rather than by dx*dy*dz which is the volume of the bounding box of the object. The modification of the bounding box in x-y plane allows R 1  to be able to distinguish between a compact-shaped object from a stick-like object lying in the lungs at a small angle to either x or y coordinate axis. Similar modifications could have been made in x-z and y-z planes. However, they are ignored because with original CT data sets the resolution in axial (z) direction is almost 10-time lower that that in x or y direction and an object&#39;s dz is thus much smaller that dx or dy. FIG. 17 explains the modification of the bounding box of an object in two dimensions. 
     2) R 2  is defined as the ratio of the maximal projection length of the object along the axis of z to the maximal projection length of the object along the axis of x or y whichever is larger. In-isotropic characteristic of CT scan is taken into account by multiplying the length (in pixel) by the corresponding pixel size. 
     
       
           R   2 = dz *psize_z/(max( dx, dy )*psize_x)  (2) 
       
     
     where psize_z and psize_x (=psize_y) are the pixel sizes (in mm) along z and x (y) directions, respectively. 
     3) R 3  is defined as the ratio of the maximal projection length of the object along the axis of x or y whichever is larger to the maximal projection length of the object along the axis of x or y whichever is smaller. 
     
       
           R   3 =max( dx, dy )/min( dx, dy )  (3) 
       
     
     where min(dx, dy)=dx, if dx&lt;dy; otherwise, min(dx, dy)=dy. 
     An object (nodule candidate) will be considered as non-compact or not within the size range of interest, and thus be deleted if the following expression is true. 
       R   1 &lt;about 0.3  ∥R   2 &gt;about 5.0  ∥R   2 &lt;about 0.2  ∥R   3 &gt;about 1.5 ∥(( dx *psize_x&lt;about 2.0)&amp;&amp;( dy *psize_y&lt;about 2.0))∥number of object voxels&gt;about 800  (4) 
     Where “∥” is the logical “OR” and “&amp;&amp;” is the logical “AND”. The threshold values for R 1 , R 2 , and R 3  are determined experimentally. The above-defined parameters and conditions ensure that the remaining objects have relatively compact shapes in three dimensions (R 1 &lt;about 0.3), two dimensions (R 2 &gt;about 5.0, R 2 &lt;about 0.2, and R 3 &gt;about 1.5), and fall into the size range of nodule of interest (about 2 mm&lt;diameter&lt;about 8 mm). Although in-isotropic character of a CT scanner is taken into account while calculating the parameter of R 2 , the threshold level of R 2  is deviated from one (1) (the closer the value of R 2  to one, the more compact shape an object has). This is because the partial volume effect along the axis of z is larger than that along the axis of x or y. 
     In addition to the above-mentioned parameters and conditions, there is at least one more situation that needs to be considered, e.g., if a nodule is attached to surrounding vessel(s). In these cases the object should not be deleted even if Equation (4) is true, provided that the object shows a compact shape on one slice, and its segments projected on the axes x and y on this slice are larger than 5 mm but smaller than 30 mm. All threshold levels in Equation (4) are determined experimentally. 
     The present invention can be implemented using a conventional general purpose digital computer or micro-processor programmed according to the teachings of the present specification, as will be apparent to those skilled in the computer art. Appropriate software coding can readily be prepared by skilled programmers based on the teachings of the present disclosure, as will be apparent to those skilled in the software art. The inventor conducted experimentation with a form of the invention that was written in the C programming language. 
     The present invention includes a computer program product which is a storage medium including instructions which can be used to program a computer to perform processes of the invention. The storage medium can include, but is not limited to, any type of disk including floppy disks, optical discs, CD-ROMs, and magneto-optical disks, ROMs, RAMs, EPROMs, EEPROMs, magnetic or optical cards, or any type of media suitable for storing electronic instructions. 
     Stored on any one of the above described storage media (computer readable media), the present invention includes programming for controlling both the hardware of the computer and for enabling the computer to interact with a human user. Such programming may include, but is not limited to, software for implementation of device drivers, operating systems, and user applications. Such computer readable media further includes programming or software instructions to direct the general purpose computer to perform tasks in accordance with the present invention. 
     The programming of general purpose computer may include a software module for digitizing and storing images obtained from the image acquisition device (computed tomography scanner). Alternatively, it should be understood that the present invention can also be implemented to process digital data derived from images obtained by other means. 
     The invention may also be implemented by the preparation of application specific units such as integrated circuits or by interconnecting an appropriate network of conventional component circuits, as will be readily apparent to those skilled in the art. 
     Thus, while there have been described what are presently believed to be the preferred embodiments of the invention, those skilled in the art will realize that changes and modifications may be made thereto without departing from the spirit of the invention, and is intended to claim all such changes and modifications as fall within the true scope of the invention.