This application claims the benefit of U.S. Provisional Application No. 60/637,805 filed Dec. 21, 2004, which is incorporated herein by reference.
BACKGROUND OF THE INVENTION The present invention relates generally to airway segmentation, and more particularly to using 3-dimensional (3-D) morphology operators for airway detection and segmentation.
A Computed Tomography (CT) scan uses x-ray equipment to obtain image data from different angles around the human body and then uses computer processing of the information to produce cross-sectional images of the body tissues and organs. The image can then be analyzed by methods using morphological operators to highlight specific areas so that radiologists (or other medical personnel) can more easily diagnose problems associated with the patient such as cancers, cardiovascular disease, infectious disease, trauma and musculoskeletal disorders. Current high-resolution CT offers high resolution images of the chest and airways.
Virtual bronchoscopy (VB) systems are systems that evaluate an image of a patient's airways for the purposes of stenosis detection, stent design, guidance, or quantitative analysis. One basic precursor that can be vital to these systems is the segmentation of the airways. Segmentation of an airway is the isolation of the airway from the rest of an image. Segmentation offers the ability to later generate paths for guidance and forms a basis for quantitative measurements.
As shown inFIG. 1, thehuman airway tree100 appears on a CT cross-section as a set of connected, dark, branching tubular structures that tend to decrease in diameter as the branching progresses.
Airway tree segmentation is traditionally a challenging problem. While airway voxels have a density near −1000 Hounsfield units (HU), noise and partial volume effects make it virtually impossible to use a simple threshold to identify all airway voxels within an image. Whenever mixtures of different tissue types comprise a voxel, intermediate gray-level values are the result. Voxels straddling air and airway walls typically have values above −1000 HU. These partial voxels can be resolved into tissue components through statistical methods. Moreover, due to the size of the voxel, thin or stenosed airways can appear broken or discontinuous. Finally, image reconstruction artifacts may cause the airways to appear discontinuous. Such discontinuities may cause problems during the segmentation, potentially resulting in both under- and oversegmentation errors.
Previously proposed airway segmentation methods have employed four strategies: (a) knowledge-based technologies, (b) region growing, (c) central-axis analysis, and (d) mathematical morphology. Hybrid algorithms combining two or more of these strategies also exist. Knowledge-based technologies describe structural relationships between airways and neighboring pulmonary vessels. Initially, 3-D seeded region growing is used to identify large airways. The knowledge-based rules are applied to the image on a section-by-section basis.
Region growing methods use voxel connectivity and a threshold to identify regions. In particular, region growing includes merging an initial set of points iteratively according to aggregation criteria. An arbitrary seed point is chosen and compared with neighboring points. A region is grown from the seed point by adding neighboring points that match aggregation criteria. This entire process is continued for each added point until no further additions are possible resulting in a single connected region.
Although 3-D region growing is typically extremely fast, it suffers from partial volume effects and noise due to the global threshold used during segmentation. The “optimal” thresholds differ for large versus small airways because of these factors. The resultant segmentation tends to lack finer details of the airways and contains rough edges. Further, these methods conventionally lose details, depict incomplete structures, or suffer from parenchymal leakage (“explosion”) to varying degrees.
Segmentation algorithms based on central axis analysis depend on central axis estimates for computing the segmentation. The disadvantage of this method is the critical dependence on the central axis analysis results, which may be imperfect or fail.
The field of mathematical morphology involves image-processing operations that focus on shape and gray-scale properties. Airway segmentation methods that draw on mathematical morphology tend to have two or more processing phases. First, candidate airways are detected by means of various morphologic operations. Next, 3-D relationships and shape characteristics help determine the true airways from false candidates. Typically, two-dimensional (2-D) operators of varying sizes are applied to each section of the image to identify candidate airways (i.e., candidates). Next, false candidates are eliminated through 3-D reconstruction. In one form of reconstruction, a region growing is performed on thresholded candidates with the seed point placed in the trachea. One potential pitfall with this approach is oversegmentation of the airways via false candidates. Further, most morphology-based methods require heavy computation times.
BRIEF SUMMARY OF THE INVENTION The present invention is a method and system for detecting a structure (e.g., an airway) within an image. The method and system generate a segmented image from the image. The method and system then iteratively apply three dimensional morphological operators of varying sizes to the image to obtain a resulting image. The method and system generate a candidate mask, apply the candidate mask to the resulting image to generate a masked resulting image, and update the segmented image using the masked resulting image to detect the structure.
In one embodiment, region growing is performed to obtain a preliminary segmented image from the original image. The image is received (e.g., from a CT scan) and thresholded to generate the candidate mask. The structure may be, for example, an airway along with airway-like objects. The iterative application of three-dimensional morphological operators further includes performing a gray scale closing with the image and an operator to obtain a detection image. A gray scale erosion is then performed with the image, the operator, and the detection image to obtain the resulting image. The size of the operator is changed to detect objects of corresponding scale.
The (original) image is then subtracted from the resulting image to obtain a difference image and the difference image is thresholded. In one embodiment, applying the candidate mask includes performing an AND operation between the candidate mask and the result of the thresholding of the difference image. In one embodiment, updating the segmented image further includes performing an OR operation between the segmented image and the result of the AND operation. True candidates are then isolated from the segmented image. In one embodiment, the isolating of the true candidates includes taking the largest component in the segmented image. In one embodiment, this candidate isolation can occur in between applications of operators, while in another embodiment, this step can occur as a single final step.
In another embodiment, a region of interest comprising at least a portion of a segmented structure (e.g. airway) is determined from a segmented image. A three-dimensional structure segmentation is performed on an image. Next, a coronal projection of the three-dimensional structure segmentation is generated. Boundary lines defining the region of interest are generated by determining profiles of the coronal projection. The profiles may include an upper profile, a lower profile, an upper cumulative profile, and/or a lower cumulative profile. The boundary lines further include upper boundary lines, lower boundary lines, and a middle boundary line. The middle boundary line divides the segmented structure into a right structure and a left structure.
In one embodiment, the upper cumulative profile is determined by determining the highest level of the upper profile. The lower cumulative profile is determined by determining the lowest level of the lower profile. In one embodiment, a shaping window is created. The lower profile may then be multiplied by the shaping window. In one embodiment, a coronal diameter of the coronal projection is determined. The shaping window may then have a non-zero value for the span of the coronal diameter. The lower profile may be determined by traversing columns of the coronal projection and determining the lowest pixel in each of the columns. The upper profile is determined by traversing the columns of the coronal projection and determining the highest pixel in each of the columns.
These and other advantages of the invention will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings.
BRIEF DESCRIPTION OF THE DRAWINGSFIG. 1 shows a structural diagram of a human airway tree;
FIG. 2 is a high level block diagram of a computer in accordance with an embodiment of the invention;
FIG. 3 shows exemplary images corresponding to the transition from an original image to a final image of candidates using a prior art approach to segmentation;
FIG. 4 shows a flowchart of the steps performed by the computer to perform the 3-D operator algorithm in accordance with an embodiment of the invention;
FIG. 5A is a rendering of a segmented image produced by region growing in accordance with an embodiment of the invention;
FIG. 5B is a rendering of an image of a segmented airway in accordance with an embodiment of the invention;
FIG. 6A is a block diagram showing the initialization steps before the iterative algorithm is performed in accordance with an embodiment of the invention;
FIG. 6B is a block diagram of an iterative algorithm to produce a final image of a segmented airway in accordance with an embodiment of the invention;
FIG. 7 shows a flow chart illustrating the steps performed by the computer to crop the volume of a segmented image in accordance with an embodiment of the invention.
FIG. 8 is a coronal projection image of a segmented airway tree with automatically determined regions in accordance with an embodiment of the invention;
FIG. 9 is a graphical representation of upper and lower profiles of the segmented airway tree in accordance with an embodiment of the invention;
FIG. 10 is a graphical representation of a shaping window in accordance with an embodiment of the invention;
FIG. 11 is a graphical representation of an upper cumulative profile in accordance with an embodiment of the invention;
DETAILED DESCRIPTION The following description describes the present invention in terms of the processing steps required to implement an embodiment of the invention. These steps may be performed by an appropriately programmed computer, the configuration of which is well known in the art. An appropriate computer may be implemented, for example, using well known computer processors, memory units, storage devices, computer software, and other components. A high level block diagram of such a computer is shown inFIG. 2.Computer202 contains aprocessor204 which controls the overall operation ofcomputer202 by executing computer program instructions which define such operation. The computer program instructions may be stored in a storage device212 (e.g., magnetic disk) and loaded intomemory210 when execution of the computer program instructions is desired.Computer202 also includes one ormore interfaces206 for communicating with other devices (e.g., locally or via a network).Computer202 also includes input/output208 which represents devices which allow for user interaction with the computer202 (e.g., display, keyboard, mouse, speakers, buttons, etc.). One skilled in the art will recognize that an implementation of an actual computer will contain other components as well, and thatFIG. 2 is a high level representation of some of the components of such a computer for illustrative purposes. In addition, one skilled in the art will recognize that the processing steps described herein may also be implemented using dedicated hardware, the circuitry of which is configured specifically for implementing such processing steps. Alternatively, the processing steps may be implemented using various combinations of hardware and software. Also, the processing steps may take place in a computer or may be part of a larger machine (e.g., a medical imaging machine).
A brief introduction to a prior art approach to the application of morphological operators is described below. In one embodiment, the prior art approach is described in A. P. Kiraly, E. A. Hoffman, G. McLennan, W. E. Higgins, and J. M. Reinhardt, “3D Human Airway Segmentation Methods for Clinical Virtual Bronchoscopy,” Academic Radiology, Vol. 9, No. 10, pp. 1153-1168 (2002), which is incorporated herein by reference. The prior art approach proceeds in two major steps. First, candidate airways (i.e., candidates) are selected via 2-D morphological operators that are of different sizes for detecting airways of different sizes and appended to an initial region growing segmentation result. Next, a reconstruction process is used to eliminate false candidates and produce the final segmentation.
Candidate airway locations are determined on a 2-D basis by means of gray-scale reconstruction followed by thresholding. The basis operator for the reconstruction is given by a binary four-connected neighborhood, B4, the smallest element used in the process. Larger elements are computed with repeated dilations of B4as follows:
B4b=bB4=B4⊕B4⊕ . . . ⊕B4,
where bB4is the bth-order homothetic of B4(i.e., bB4is a geometric expansion or contraction of B4, depending on b).
An operator B4b is applied to each individual (2D) section as described below. Given the original image I, each section z=1, . . . , Z of I is windowed (−1000, 0) to form a new 2-D image S, where
Windowing is used to eliminate the effect of variations in more dense structures, such as bone. A marker image for gray-scale reconstruction, J1, is then obtained from the gray-scale closing of S with structuring element B4b,
J1b=S·B4b=(S⊕B4b)ΘB4b.
Hence, J1bis computed from the gray-scale dilation by means of B4bfollowed by an erosion by the same operator. Next, J2bis computed from J1bwith the following:
Jk+1b=max(JkbΘB4,S),
where the right side of the equation computes the voxel-by-voxel gray-level maximum.
The above equation is repeated until no further changes occur; that is,
J∞b(x,y)=Jk+1b(x,y)=Jkb(x,y),∀(x,y)εJkb,
where J∞brepresents the final gray-scale reconstructed image with structuring element B4b. In this image, local minima smaller than B4bin S are filled in with a gray-level value proportional to the difference between the maximum and minimum gray levels computed within a B4b-sized neighborhood of the minima.
A gray-scale difference image is then computed between j∞band S and thresholded:
The difference image is bright where local minima exist in the image S. In this case, since the data windowed was between 0 and −1000, the threshold value is 200.
Finally, the union of all the candidates obtained from all of the operators is taken to create the final image C:
The image C is then reconstructed to eliminate false candidates and later fused with the region growing segmentation for the final result.FIG. 3 shows exemplary images corresponding to the transition from the original image304 to the J image308 after five iterations (using one operator) to the J infinity image312 to an image Cbobtained from a specific operator size b316 using the prior art approach to segmentation. The white regions (circles, spots, speckle) are potential airway candidates. The collection of these 2-D results will then be considered in unison to identify the true airway.
Unlike the prior art approach, which uses 2-D operators, the present invention uses 3-D morphological operators. The use of 3-D morphological operators allows for more image information to be used in context in order to detect airway candidates. Additionally, the threshold is varied based on operator size and the reconstruction stage is merged with the first stage and is iteratively applied as the results from each operator are computed.
The 3-D operator algorithm is based on simultaneously considering multiple adjacent sections (i.e., slices) of the original image. Thus, if a dark spot is present on one slice of the original image (e.g., a dark spot associated with a packet of fat) but is not present on other slices of the image, the 3-D operator algorithm determines that this dark spot is not an airway. Unlike a 2-D operator algorithm, which may determine that the dark spot is an airway because of the examination of that single slice of the image, the 3-D operator algorithm is applied to multiple slices of the image so that a dark spot present on a single slice will not be considered as an airway.
FIG. 4 shows a flowchart of the steps performed for the 3-D operator algorithm. First, an original image I of the chest (i.e., airways) is obtained from a CT scan. To create an initial estimate segmentation of the major airways image (and thereby create a segmented image), region growing is performed at a seed point (e.g., at the trachea) instep404. As shown inFIG. 5A, theregion growing image504 typically only grows to a particular area and does not cover the smaller airways (e.g., airway510) present in asegmented image508 ofFIG. 5B.FIG. 5B shows a rendering of a segmented airway and therefore does not show the rest of the image (that may have other dark areas, such as correlating to emphysema).Image504 depicts a surface rendering of the results obtained from region growing image which is assigned to the “current segmentation” (also referred to below as image H) that provides an initial segmentation that will be completed during later steps.
An airway candidate mask (also referred to as image G) is then obtained instep408. The airway candidate mask is obtained by performing a threshold on the original image I:
G(x,y,z)=1 ifI(x, y, z)≦AirThreshold, 0 otherwise
In one embodiment, a threshold of −675 HU is used. Alternatively, a threshold of −775 HU is used. This thresholding step may be performed at any time during the algorithm before image G is required—it is only performed once and used by all iterations.
Candidates for airways are then determined by a specific-sized operator. Note that these candidates may include airways already defined in image H, but the goal is to find additional airways. This determination is made by using a 3-D, b-sized operator on the original image I instep412 and as shown below:
J1b=I·B4b=(I⊕B4b)ΘB4b
Jk+1b=max(JkbΘB4,I)
Note in this case that the operator of size b is not necessarily obtained from repeated dilations and can be uniquely obtained. For example, b could correspond to the radius of a sphere that is digitized to create the operator. It can also be elliptical in the case of anisotropic data.
The difference image between the image Jk+1b(also referred to below as the resulting image) and the original image I is then calculated instep416. The difference image is then thresholded. The threshold varies for different b, and typically decreases as b decreases.
Cb(x,y,z)=1 ifJ∞b(x,y,z)−I≧Tcb, 0 otherwise
The thresholded difference image is then masked (i.e., ANDed) with the airway candidate mask instep420 to generate a masked resulting image.
Cb(x,y,z)=Cb(x,y,z)·G(x,y,z)
Image H is then updated with the image generated in the previous step (the masked resulting image) instep424 by ORing.
H(x,y,z)=H(x,y,z)|Cb(x,y,z)
The above steps are iterated for different sized structuring elements. For example, thecomputer202 may empirically determine that the largest operator that had useful additions to the region growing was a sphere with a discrete radius of 3 (0.5 mm scale). The above steps are iteratively performed using operators of different sizes (e.g., ofradius 3, 2, and 1 in that order). Thus, the iterations progress for operator size b⊂{M, M−1, . . . 1}, where M is the determined maximum operator size that had useful additions to the region growing image H.
The next step is to isolate the true candidates from the image instep426. In one embodiment, the largest connected component is extracted instep426. Image H is then updated with the true candidates.
The size of the operator is then compared to 1 instep428. If the size b is greater than or equal to 1, then another iteration is performed after the size is decreased (e.g., by one) instep427. If the size b is less than 1, then the last iteration has been performed. Hence, unlike previous methods, the candidate isolation occurs iteratively along with each operator.
FIGS. 6A and 6B show more detailed flow charts illustrating an embodiment of the algorithm used to generate the airway segmentation.FIG. 6A shows the initialization steps before the iterative algorithm is performed. As described above, an original image I604 is obtained via a CT scan of the chest. In this example, showing a synthetic cross-section of an airway, theoriginal image604 includes agray background608, a grayouter ring612, and a blackinner circle616. In one embodiment, a 2-D slice of an airway is shown as a sample of a 2D section of the 3D image. Region growing is then performed at a seed point in step620. Asegmented image H622 is set to the region growing result.
An air threshold is also performed on the original image I604 instep624 to capture air-like regions of theoriginal image I604. The resultingimage G625 has ablack background626 and a whiteinner circle627 for this example. In general, it is white for all voxels that were less than the threshold in the original image. The purpose of this image is to define potential candidates based on gray level values.
FIG. 6B shows a block diagram of an iterative algorithm to produce a final image of a segmented airway. A gray scale closing is performed on the original image instep628 using the operator described above,Bb630. This closing628 producesdetection image J1632 for the first operator. A gray scale erosion with theoperator Bb630 is then performed inblocks634 followed by a MAX operation with the original image I on thedetection image J1632 as well as theoriginal image604. This loop is repeated until the input image into theMAX operation636 is unchanged. The result is a resultingimage J∞b638 which is referred to asJ∞638 for the remainder of the iteration. The resultingimage J638 has a gray background, a gray outer ring and a gray inner circle in this example. In general, some darker portions of the image become brighter.
Thecomputer202 subtracts theoriginal image604 from the resultingimage638 to produce a subtracted ordifference image639. The subtractedimage639 has a black background and a gray inner circle in this example. The subtractedimage639 is then thresholded to generate athresholded image640. The computer masks (i.e., performs an AND operation on) thethresholded image640 with the airwaycandidate mask G625, as shown withstep642, to generate a masked resulting image. The masked resulting image is ORed with the currentsegmented image H622 to produce a newsegmented image644. This image is iteratively updated as the operator size b is adjusted. Once the iterations are complete, the true candidates are isolated from thesegmented image H622. In one embodiment, the largest component646 in the image is taken to create asegmented image650. Thesegmented image650 has a black background and a white inner circle in this example, thereby segmenting the airway (i.e., the white inner circle). Thisimage650 is used to update image H. If necessary, the process is iterated more than once. Further, this operation can occur during each iteration of operator sizes.
Because all operations are performed in 3-D, the initial segmentation is iteratively grown by decreasing the operator size and updating the segmentation for each operator size. This algorithm enables the avoidance of more false positives than if thecomputer202 used only 2D operators as done in the prior art algorithm. Further, morphological operations with much larger structuring elements are avoided since the initial segmentation created from region growing already contains the larger airways.
The main parameters involved for candidate selection are:
- 1. M, the maximal size of the operators (e.g., set to 3 but will increase as the resolution of the input image increases).
- 2. Tcb, the threshold for the detection difference image at each iteration. This is different for every iteration and typically decreases along with b.
- 3. The connectivity constraint for extracting the largest connected component in the last step. This may be different for different iterations.
In one embodiment, the 3-D operator algorithm requires a significant amount of the computer's memory to execute. To solve this problem, thecomputer202 may crop the entire thoracic volume to a region of interest (ROI) which, in this case, is the lung region that contains identifiable airways. If further partitioning were necessary, separate bounding boxes for the left and right airway tree may be identified.
FIG. 7 shows a flow chart illustrating the steps for determining, from a segmented airway image, a region of interest having at least some of the segmented airway (i.e., volume cropping). Thecomputer202 performs a 3-D airway segmentation on an image as described above instep704. Thecomputer202 obtains a coronal projection of the 3-D airway segmentation instep708. Thecomputer202 then determines boundary lines defining the region of interest by determining profiles of the coronal projection followed by sagittal projection instep712.
The profiles include an upper profile, a lower profile, an upper cumulative profile, and a lower cumulative profile. The boundary lines can include upper boundary lines, lower boundary lines, and a middle boundary line that divides the segmented airway into a right airway and a left airway. The determination of an upper cumulative profile includes determining the highest level of the upper profile. The determination of the lower cumulative profile includes determining the lowest level of the lower profile. The region of interest is defined by the boundary lines and may include any portion of the airway.
FIG. 8 shows a diagram of a segmented airway tree obtained using region growing. The diagram is a coronal maximum intensity projection (MIP)800 ofsegmented airway tree802.Vertical line804 separates the left from the right airway, while the upperhorizontal lines806,808 provide the upper bounds of the left and right airway tree separately.
In one embodiment, thecomputer202 determines the upper and lower profiles of the coronalairway segmentation MIP800.FIG. 9 shows a graph of the upper profile904 and thelower profile906. Thelower profile906 is obtained by traversing the columns of theMIP800 and noting the height at which the lowest airway pixel is located in each pixel column. The upper profile904 is obtained by traversing the pixel columns of the MIP and locating the highest airway pixel of each column.
The highest point of the lower profile is typically where the trachea splits. With a more extensive segmentation, however, a branch at the leftmost or rightmost end may have the highest elevation. Therefore, the lower profile is multiplied with a shaping window such that the highest point in the shaped profile coincides with the point at which the trachea splits into the left and right airway. InFIG. 8, the left boundary and right boundary of the airway tree are shown asvertical lines810,812. The distance between the left boundary and the right boundary is the coronal diameter of the airway tree. In one embodiment, the shaping window is set to a non-zero value for the span of the coronal diameter.FIG. 10 shows an exemplary shaping window1000 having alinear ramp1004,1008 for the outer 20% of the coronal diameter on either side.
To obtain thehorizontal boundary lines806,808, thecomputer202 analyzes the upper profile904 ofFIG. 9. Thecomputer202 creates an upper cumulative profile by traversing theupper profile804 from the outside towards the left/right splitting boundary802. At each point, the value of the upper cumulative profile is the highest level that thecomputer202 sees thus far when traversing the upper profile904: i.e., it is the maximum level seen considering the section of the upper profile904 that is bounded by the relevant outer boundary and the current point. The traversal of the upper profile904 results in an upper cumulative profile, shown as theupper curve1100 inFIG. 11. The upper cumulative profile is flipped and the top of the trachea is set at level 0. Thus, the obtained boundary levels will be relative to the top of the trachea. The desired boundary is the line that touches the upper most distal airway of the airway tree of either side. Referring tocurve1100, it is a rather long run in the upper cumulative profile that is proximal to the trachea, or equivalently, toboundary802. Similar to the motivation that led to the shaping window inFIG. 10, the region of theupper curve1100 that is to be used to obtainboundary lines806 and808 can be specified in terms of percentages of the airway tree width and relative toboundary802. Thelower curve1104 shows the first difference ofcurve1100—specifically, the spikes in1104 denote a change in level incurve1100. This information is typically simpler to analyze in light of prior knowledge—to define the region of1100 that dictates theboundary lines806 and808. Once the boundary values that can be found using a coronal MIP have been obtained, the next step is to obtain the anterior and posterior boundaries for the left and the right airway. Thecomputer202 uses the boundaries (especially the left and right splitting boundary) to obtain separate sagittal Maximum Intensity Projection (MIPS) for the left and right airway trees. Once these sagittal MIPS are available, the required boundaries are obtained.
The boundary requirements can be subdivided further into quadrants or even further subdivisions for stricter memory requirements of thecomputer202. These further subdivisions can be used to set custom parameters for region growing, morphology, or other operations to provide better results.
Thus, the advantage of volume cropping is the ability to focus separately on different regions of the airway segmentation for speed and memory efficiency.
Although the above algorithms are described in terms of an airway, they may equally be applied to any structure in an image.
The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.