## Figure

3-D segmentation stages: (a) presegmented data, (b) seeds set, (c) segmentation complete, and (d) generated paths.

knowledge of their segmented regions, there is no hint that the vessel should be continued across a complete gap. For the human operator, however, such failures are immediately obvious and can be corrected by manually bridging the gap using the path editing interface.

### 3.7.5 Performance

The performance of the three-dimensional segmentation is data dependent. Both the size of the region to be segmented and the intensity distribution affect the running time, as of course does the seeding pattern given by the user.

Putting those factors aside for the moment, we focus on the real-world performance of the application. The deployment platform is a 1.5-GHz Intel machine with 1 to 2 GB of RAM. The typical CT abdominal study contains 250 to 400 512 x 512 slices. Three-dimensional segmentation of the aorta from above the renal branch to below the iliac bifurcation generally requires less than 5 seconds; the additional path generation and visualization steps add an additional 5 to 10 seconds to the total execution time.

MR studies are faster because of the lower resolution. In CT, other narrow-vessel studies such as coronaries and peripherals require less time because the total volume of the segmented region is usually significantly smaller than that of the large aorta.

3.7.6 From Coarse Three-Dimensional to Fine Two-Dimensional Segmentations

Extraction of vessels in three dimensions using the algorithm described provides the user with a voxel-resolution segmentation of the vessel and a centerline that can be used to accurately position an orthogonal MPR plane. In this plane, the vessel diameter and cross-sectional area can be accurately measured. A goal of our implementation is to provide vessel cross-sectional area measurements automatically once the three-dimensional segmentation is completed. These measurements should be given with sub-voxel precision to support evaluation in cases where the vessel may only have a diameter of 2 or 3 voxels. To accomplish this goal, we use a two-dimensional ray propagation method in conjunction with a mean-shift-based noise modulation technique.

3.8 Cross-Sectional Vessel Boundary Detection by Two-Dimensional Ray Propagation

In this section, we describe a ray propagation technique for detecting vessel boundaries in the planes orthogonal to their centerlines [83].

Let the front be represented by a two-dimensional curve C(s, t) = (x(s, t), y(s, t)), where x and y are the Cartesian coordinates, s is the length parameter, and t is time. The evolution is then governed by [51]:

where Co(s) = (x(s, 0), y(s, 0)) is the initial curve,N is the unit normal vector, and S(x, y) is the speed of a ray at point (x, y).

The approach that we consider in this section is based on explicit front propagation via normal vectors. Specifically, the contour is sampled and the evolution of each sample is followed in time by rewriting the Eikonal equation [69] in vector form, namely:

The speed of rays, S(x, y) depends on image information and shape priors, which is explained in the next sections.

This evolution is Lagrangian because the physical coordinate system moves with the propagating wavefront. However, the applications of ray propagation for curve evolution have been limited. This is because, when the normals to the wavefront collide (formation of shocks), this approach exhibits numerical instabilities due to an accumulated density of sample points, thus requiring special care, such as reparametrization of the wavefront. Also, topo-logical changes are not handled naturally; that is, an external procedure is required.

Previously, ray propagation has been used to implement the full-width at half-maximum (FWHM) technique for quantification of two-dimensional airway geometry [31]. In the FWHM technique, the maximum and minimum intensity values along the rays are computed to determine the "half-maximum" intensity value, which is the half-intensity value between the maximum and minimum. However, stable computation of the maximum and minimum along the rays is quite difficult due to high signal variations. In addition, ray propagation was used by Wink et al. [98] for fast segmentation of vessels and detection of their centerline. In this approach, Wink et al. used the intensity gradients to stop the propagation of rays. However, this approach faces difficulties when the vessel boundaries are not sharp (due to partial volume effects) and also when vessels contain isolated noises (e.g., calcifications in CT images).

The main shortcomings of these approaches stem from the computation of image gradients that are not robust relative to the image noise. In this section, we propose to use mean shift analysis for detecting vessel boundaries efficiently and robustly. Comaniciu and Meer [23, 24] showed that the image discontinuities are robustly revealed by a mean-shift process that evolves in both the intensity and image space.

We first review nonlinear smoothing algorithms; second, describe the mean shift analysis; third, illustrate an approach where the mean shift procedure is applied to select discontinuities in one dimensional signal; and finally, present the mean-shift-based ray propagation approach for segmentation of medical structures (e.g., vessels).

### 3.8.1 PDE-Based Smoothing Methods

Often, smoothing is used to remove noise present in three-dimensional images without affecting the important features, such as vascular structures. Traditionally, images are smoothed via the Gaussian function. However, Gaussian smoothing is not an appropriate method for smoothing medical images containing vascular structures because it shrinks shapes and dislocates boundaries when moving from finer to coarser scales. Instead, nonlinear filters are proposed for removing the noise present in three-dimensional images while keeping the important geometrical structures.

Perona and Malik proposed a nonlinear diffusion filtering method that reduces smoothing at the edges to preserve the contrast information and the location of the object boundaries [70]. Specifically, they proposed that:

fut = div( f ( IV u 12)V u) 1 u(x, y, 0) = I (x, y)

where I (x, y) is the image, f (IVu12) = 1+Ivu^, and \ > 0. This method preserves the strong edges for a long time during the smoothing process by turning off the smoothing at high brightness gradient locations. Thus, this approach smooths regions of low brightness gradient while regions of high gradients are not smoothed. Catte et al. [17] suggest a slight modification in which the gradient in f( Vu 2) is computed from the Gaussian smoothed image to account for the large noise in the image:

A similar scheme was proposed by Whitaker and Pizer [96].

However, these anisotropic diffusion approaches are contrast driven and smooth globally salient but low contrast image features [50]. Thus, geometry-driven diffusion techniques are proposed to deal with low-contrast structures. One particular approach is based on curvature-based smoothing of grey-level intensity images. The curvature deformation of surface $(x, y, t), which is also known as the mean curvature flow,

was studied by Evans and Spruck [34] and by Chen et al. [18]. Kimia and Siddiqi [50] also used this evolution for image smoothing. Similarly, Alvarez et al. [3] studied a class of nonlinear parabolic differential equations defined by:

where G is a smoothing kernel (e.g., the Gaussian). This is similar to Perona Malik's anisotropic smoothing in that curvature smoothing is turned off in the vicinity of a high brightness gradient. Osher and Rudin [68] proposed shock filters for image smoothing and enhancement. Other works on the nonlinear geometric diffusion can be found in [82, 88, 89]. Specifically, Weickert presented an excellent review of the related approaches in [93]. We now present mean-shift analysis, a robust way of measuring image discontinuities and smoothing images at the same time.

3.8.2 Mean-Shift Analysis

Given the set {x; }i=1,.,n of d-dimensional points, the mean-shift vector computed at location x is given by [23, 24]:

where K represents a kernel with a monotonically decreasing profile and h is the bandwidth of the kernel.

It can be shown that (Equation 3.31) represents an estimate of the normalized density gradient computed at location x; that is, the mean shift vector always points toward the direction of the maximum increase in the density. As a result, the successive computation of expression (Equation 3.31), followed by the translation of the kernel K by Mh (x), will define a path that converges to a local maximum of the underlying density. This algorithm is called the mean shift procedure, a simple and efficient statistical technique for mode detection.

The mean-shift procedure can be applied for the data points in the joint spatial-range domain [23, 24], where the space of the two-dimensional lattice represents the spatial domain and the space of intensity values constitutes the range domain. In this approach, a data point defined in the joint spatial range domain is assigned with a point of convergence that represents the local mode of the density in this space, e.g., a three-dimensional space for gray-level images. One can define the displacement vector in the spatial domain as the spatial difference between the convergence point and the original point.

When each pixel in the image is associated with the new range (intensity) information carried by the point of convergence, the algorithm produces discontinuity, preserving smoothing. Conceptually, this process is similar to anisotropic diffusion, nonlinear filtering, or bilateral filtering [9].

However, when the spatial information corresponding to the convergence point is also exploited, one can define a segmentation process based on the displacement vectors of each pixel. The convergence points sufficiently close in this joint domain are gathered together to form uniform regions for image segmentation [23, 24].

In this section, we will exploit the mean-shift-generated displacement vectors to guide active contour models. The robustness of the mean shift is thus combined with a priori information regarding the smoothness of the object contours. This processing is integrated into our computationally efficient framework based on ray propagation. The new algorithm allows real-time segmentation of medical structures.

### 3.8.3 Mean-Shift Filtering Along a Vector

Let us first illustrate the mean-shift procedure applied to a one-dimensional intensity profile obtained from a two-dimensional gray-level image. Specifically, let {x\, f }i=1i...,N and {x*i, I*i}i=1/.../N be the two-dimensional original and filtered N image points in the spatial-range domain. In addition, the output of the mean-shift filter includes a displacement vector {di}i=i/.../N that measures the spatial movement of each spatial point. In our algorithm, each point in this spatial range domain is processed via the mean-shift operator until convergence. Specifically, the algorithm consists of three steps:

1. Initialize k = 1 and (x'k,1'f, di) = (xi, Ii, 0)

until the displacement of spatial points xi are small (i.e., |x ; 1 -x'k | < e) and where ux and uI determine the Gaussian spatial and range kernels size, respectively.

Observe that in the last step of the procedure, the original spatial locations (namely, xi's) are assigned with the smoothed intensity values. Figure 3.10a illustrates an example where 1-dimensional intensity data is obtained from a slice of a CT image. Figures 3.10b and c illustrate the original intensity profile and the smoothed intensity profile, respectively. Observe that the mean shift procedure smooths the intensity data while preserving and sharpening its discontinuities. Similarly, Figure 3.10d depicts the displacement vectors along this 1-D signal. Our boundary detection framework exploits the information contained in these displacement vectors.

## Post a comment