Semiautomatic Parametric Model-Based 3D Lesion Segmentation for Evaluation of MR-Guided Radiofrequency Ablation Therapy1

Original Investigations

Semiautomatic Parametric Model-Based 3D Lesion Segmentation for Evaluation of MR-Guided Radiofrequency Ablation Therapy1 Roee S. Lazebnik, MD, PhD, Brent D. Weinberg, BS, Michael S. Breen, PhD, Jonathan S. Lewin, MD, David L. Wilson, PhD2

Rationale and Objectives. Interventional magnetic resonance imaging (iMRI) allows real-time guidance and optimization of radiofrequency ablation of pathologic tissue. For many tissues, resulting lesions have a characteristic two-boundary appearance featuring an inner region and an outer hyper-intense margin in both T2 and contrast-enhanced (CE) T1-weighted MR images. We created a geometric model-based semiautomatic method to aid in real-time lesion segmentation, crosssectional/three-dimensional visualization, and intra/posttreatment evaluation. Materials and Methods. Our method relies on a 12-parameter, 3-dimensional, globally deformable model with quadric surfaces that describe both lesion boundaries. We present an energy minimization approach to quickly and semiautomatically fit the model to a gray-scale MR image volume. We applied the method to in vivo lesions (n ⫽ 10) in a rabbit thigh model, using T2 and CE T1-weighted MR images, and compared the results with manually segmented boundaries. Results. For all lesions, the median error was ⱕ1.21 mm for the inner region and ⱕ1.00 mm for the outer hyper-intense region, values that favorably compare to a voxel width of 0.7 mm and distances between the borders manually segmented by the two operators. Conclusion. Our method provides a precise, semiautomatic approximation of lesion shape for ellipsoidal lesions. Further, the method has clinical applications in lesion visualization, volume estimation, and treatment evaluation. Key Words. Parametric Deformable Model Segmentation; Interventional Magnetic Resonance Imaging; Radiofrequency Thermal Ablation; Medical Image Processing. ©

AUR, 2005

Acad Radiol 2005; 12:1491–1501 1 From the Department of Biomedical Engineering (R.S.L., B.D.W., M.S.B, D.L.W.) and the School of Medicine (R.S.L., B.D.W.), Case Western Reserve University, Cleveland, OH; the Department of Radiology, University Hospitals of Cleveland, Cleveland OH (D.L.W.); and the Department of Radiology, The Johns Hopkins Hospital, Baltimore, MD (J.S.L.). Received April 18, 2005; revision received July 20; accepted July 23. Sponsored by NIH grant RO1-CA84433 to DLW. RSL was supported by a Whitaker Foundation graduate fellowship. RSL and BDW are supported by the CWRU Medical Scientist Training Program. We thank Sherif G. Nour for his help with the in vivo rabbit procedure. Address correspondence to: D.L.W. e-mail: [email protected] 2 Present address: Department of Biomedical Engineering, 10900 Euclid Ave., Cleveland, OH 44106

© AUR, 2005 doi:10.1016/j.acra.2005.07.011

Interventional magnetic resonance imaging (iMRI) relies on MR images for guidance during minimally invasive medical procedures (1–3). Among such procedures is thermal ablation of pathologic tissue using a radiofrequency (RF) energy source (4). Imaging of ablated tissue reveals a significant change of the MR signal in the vicinity of the current source. In many tissues, including liver and muscle, the characteristic cross sectional lesion appearance features a hypo-intense inner zone and a hyperintense outer zone (5–7). The lesion also possesses approximately ellipsoidal geometry (8). Studies of ablated gelatin phantoms and ex vivo liver tissue observed spheroidal lesions symmetric about the exposed tip of the current delivery electrode (9,10). The approximately ellipsoidal configuration is also predicted by numerical simula-

1491

LAZEBNIK ET AL

tion of the current, electric field, and resulting temperature field (11). Clinically, the inner zone depicted by a T2-weighted MR image is used to estimate the ablated tissue boundary. Previous studies suggest the hypo-intense zone approximately corresponds to tissue destruction, whereas the hyper-intense zone corresponds to an edematous or hemorrhagic rim surrounding the ablated area. However, it is unknown whether the peripheral cells within the hypointense region and those corresponding to the hyper-intense region are truly nonviable or sufficiently damaged so that apoptosis or osmotic destruction is guaranteed. We are investigating the relationship between histologically determined boundaries of cell destruction and the MR signal (12,13). Image segmentation is often a difficult task because of both variability of object shapes and image quality. One solution is preprocessing of images to remove noise and postprocessing to correct shapes, which may require extensive human observer interaction. To address these limitations, a class of methods known as deformable models was conceived (14). In general, deformable models constrain allowable shapes by accounting for internal and external forces defined within the surface itself and computed from image data, respectively. Internal forces are designed to keep the model smooth during deformation while external forces move and deform it toward shape boundaries or image features. Advantages of deformable model techniques include potential subvoxel accuracy, the ability to mathematically constrain boundaries using a priori information, as well as resistance to image noise and boundary gaps (particularly for global deformation methods). Most generally, approaches may be classified into parametric, defined by a boundary model of with a fixed number of parameters, and nonparametric where the boundary is a mesh with an arbitrary number of control points. Many variations of deformable model algorithms are available with differing level of automation, accommodation of topological complexity, incorporation of a priori constraints, dependence on initial values, and convergence speed. Among the most cited approaches is the parametric snake active contour (15). More recent approaches include implicit geodesic deformable models and Bayesian probability frameworks. Geodesic models are not strictly parameterized, rather defined by a zero level set contour of arbitrary dimensions with well-defined evolution function. Although more accommodating of variations in expected shape, this results in decreased resistance to image

1492

Academic Radiology, Vol 12, No 12, December 2005

noise compared with globally deformable models (16). Bayesian frameworks constrain manipulation of exemplary model templates and allow for parameterized physical transformations such as scaling and bending (17). The most cited application of deformable models to medical imaging is the segmentation of left cardiac ventricle shape. Previous studies utilized a parametric superquadric ellipsoidal approximation of the endocardium and epicardium to model ventricular deformation during contraction (18,19). Shape extraction followed an energy minimization approach, using an internal regularization term and an external attraction potential. Another study used a global deformable ellipsoidal ventricle model with parametric offsets to describe cardiac deformation observed with tagged MR data. The model was applied using the constraint that cardiac wall volume remains constant during contraction, highlighting the potential for inclusion of a priori information in model based methods (20). A study of tagged MR data applied a volumetric deformable model incorporated Lagrangian dynamics and finite element theory (21,22). An alternative approach relied on a kinematic model incorporating myocardial fiber orientation into deformable hexahedral elements (23). These demonstrate the application of biomechanical functions as parameter constraints. In a previous study, we quantitatively validated a double ellipsoid model of lesion shape using digital phantom and in vivo data (8). For all in vivo lesions, the median distance from the model surface to segmented boundaries was ⱕ0.58 mm (interquartile range ⱕ0.89 mm) for both the inner and outer surfaces, values less than the voxel width of 0.7 mm. These values include lesions spanning connective or heterogeneous tissue, and in vicinity of vasculature. We thus showed that the model is sufficiently descriptive of in vivo lesion geometry and accommodates several potential asymmetries. Because of the previously mentioned power of model-based segmentation, we now extend the concept to semiautomatic segmentation. We present a method to segment a lesion in three dimensions (3D) using a globally deformable ellipsoidal model of both boundaries. Global deformations are particularly well suited for this application as low field iMRI is often noisy. To our knowledge, this is the first attempt to apply an automated segmentation method to RF ablation data. The model is semiautomatically fitted to a grayscale volume through minimization of an energy function that simultaneously accounts for lesion shape and grayscale properties. We apply the method to in vivo images of lesions in an animal model and quantify the model fit

SEMIAUTOMATIC PARAMETRIC MODEL-BASED 3D LESION

Academic Radiology, Vol 12, No 12, December 2005

Figure 1. Contiguous contrast enhanced (a) T1- and (b) T2-weighted magnetic resonance image slices and corresponding overlay from the automatic segmentation method. Two characteristic boundaries are visible, one surrounding the inner core and the other surrounding a hyper-intense margin. Note the elliptical cross section of both boundaries.

to actual lesion geometry. There are several potential applications for our method, including pseudo real-time two-dimensional and 3D visualization, geometric dimension measurements, and quantification of imaging signal properties within the treated region as they relate to tissue response (13). Our method can thus provide an interventional radiologist with automated treatment evaluation information.

ALGORITHM Model Geometry An ellipsoid is a quadric surface whose radii occur along the principal axes are determined by the constants a, b, and c as shown by Eq 1. x2 a

⫹

y2 b

⫹

z2 c

⫽1

(1)

tively, radii rx, ry, and rz, to provide x, y, and z values along the ellipsoid surface (24).

冤冥 冤 x

y ⫽ ry · sin u · cos v z

rz · sin v

冥

(2)

An advantage of Eq 2 is the ability to sweep a surface by systematically varying u and v rather than solving Eq 1. Using this form, we can easily translate the center from the coordinate system’s origin along each axis and center by tx, ty, and tz. Further, we can rotate the surface about the coordinate system’s axes (␪x, ␪y, and ␪z), before translation, through multiplication by Givens matrices (Rx, Ry, and Rz). There are thus nine parameters in this formulation of the ellipsoidal surface.

冤冥 x

y ⫽ Rx · R y · Rz ·

Figure 2b depicts an example ellipsoid surface. In general, the ellipsoid can be translated to a new point along the x, y, and z axes or rotated about any of these axes. An alternative parameterization of an ellipsoid surface uses the spherical coordinate system angles u and v, varied between ⫺␲ ⬍ u ⬍ ␲ and ⫺␲/2 ⬍ v ⬍ ␲/2, respec-

rx · cos u · cos v

z

冤

1

Rx ⫽ 0

冤

rx · cos u · cos v

冥 冤冥 tx

ry · sin u · cos v ⫹ ty rz · sin v 0

0

cos ␪x

sin ␪x

0 ⫺sin ␪x cos ␪x

(3)

tz

冥 1493

LAZEBNIK ET AL

Academic Radiology, Vol 12, No 12, December 2005

Figure 2. Typical manually segmented lesion boundaries and corresponding automatic segmentation. Units are in millimeters. Manually segmented boundaries, for every applicable slice, define the inner core and hyper-intense margin (a). Subsequent to model fitting, both the inner and outer ellipsoidal model boundaries are visualized in three dimensions using surface rendering (b).

Ry ⫽

冤 冤

cos ␪y 0

0 sin ␪y 1

0

⫺sin ␪y 0 cos ␪y cos ␪z

sin ␪z 0

Rz ⫽ ⫺sin ␪z cos ␪z 0 0

0

1

冥

冥

volume of signed grayscale edge intensities approximating the 3D gradient of V.

ˆ v ⬵ ⵜV ⫽ G

⭸x

ˆi ⫹

⭸V ⭸y

ˆj ⫹

⭸V ⭸z

kˆ

(5)

(4)

To accommodate both lesion boundaries, we create a model consisting of two concentric ellipsoidal surfaces. This gives us two sets of three radii, for the inner and outer ellipsoid, whereas other parameters are identical for the two surfaces. Thus the final model includes 12 paramout eters: tx, ty, tz, ␪x, ␪y, ␪z, rxin, riny , rinz , rxout, rout y and rz . We will refer to this set of parameters as ␾, its subset defining the inner ellipsoid as ␸in and subset defining the outer ellipsoid as ␸out. Model Energy To fit our model to a gray-scale volume V(x, y, z), we will define model energy such that its minimization will optimize model parameters to fit a grayscale volume. We ˆ v, a vector start by computing a directed edge gradient G

1494

⭸V

ˆ v,(x, y, z) is computed using numerical central differencG ing along each direction, accounting for voxel dimensions (dx, dy, and dz), to approximate the first derivative.

ˆ v(n, m, k) ⫽ G

V(n ⫹ 1, m, k) ⫺ V(n ⫺ 1, m, k) 2 · dx ⫹ ⫹

ˆi

V(n, m ⫹ 1, k) ⫺ V(n, m ⫺ 1, k) 2 · dy V(n, m, k ⫹ 1) ⫺ V(n, m, k ⫺ 1) 2 · dz

ˆj kˆ (6)

ˆ v, we independently compute the compoTo compute G nents Gx, Gy, and Gz, along the x, y, and z directions respectively, using central differencing as indicated above.

SEMIAUTOMATIC PARAMETRIC MODEL-BASED 3D LESION

Academic Radiology, Vol 12, No 12, December 2005

We know a priori that the inner surface is on a positive (dark to light) edge from the lesion center, whereas the outer is along a negative (light to dark) edge. To determine the sign of the edge, we compute the projection ˆ v, on a unit vector û pointing from the of each vector in G model center to the corresponding spatial location.

ˆ v (n, m, k) 䡩 uˆ (n, m, k) Gs(n, m, k) ⫽ G

(7)

The resulting volume Gs is a 3D gray-scale edge map with positive values for rising edges and negative values for falling edges, as one moves radially from the model center. We define the model energy objective function that accounts for model position with respect to gray-scale edges.

E(␸) ⫽ ⫺

1

n

兺

n i⫽1

Sin(ui, vi, ␸in) ⫹

1

p

兺S

p j⫽1

out

(u j, v j, ␸out) ⫹ P(␸) (8)

We choose n and p equally spaced points, along the inner and outer ellipsoidal model surfaces respectively, where we will compute a term in the objective function. We define Sin(um,vm,␾in) as the value of Gs at a point m along the inner ellipsoidal surface defined by {u,v} and model parameters ␾in. There is also a corresponding Sout(um,vm,␾out). Note that the sampling points are not discretized to voxel locations, but rather gray-scale values are obtained at real valued locations using trilinear interpolation of Gs to make parameter space more continuous. Finally, we define a penalty function P(␾) such that if any inner radius exceeds the corresponding outer radius, the function equals the difference times a large constant, else the function is zero valued. P(␸) ⫽ 兵 if (rin ⬎ rout) then P(␸) ⫽ C * (rin ⫺ rout) else P(␸) ⫽ 0 (9) P(␾) could also similarly account for additional a priori constraints such as ensuring that equatorial radii do not differ by more than a specified amount. Equation 8 summarizes the composite energy function E(␾). Note that the first and second terms are negatively

and positively signed, respectively, because minimization of E(␾) results in large positive edge values along the inner boundary and large negative values along the outer boundary. Parameter Estimation We apply an iterative optimization technique to minimize the model energy E(␾) in Eq 8. There are many algorithms suitable for a nonlinear search of parameter space. Typically the choice of algorithm depends on the tradeoff between number of objective function evaluations and resistance to local minima. Because our objective function is computationally inexpensive, such that each evaluation requires only a fraction of a second on a Pentium 4 class computer, we use the Nelder-Mead simplex method (25) for its robustness and resistance to discontinuities in parameter space despite a large number of expected iterations. Initial values for the parameter search are determined from 12 operator-selected points. First, the operator chooses two points that define each of the inner and outer maximum visualized diameters along each axis. Typically, these are obtained along the needle axis perpendicular to the image slices. A set of parameters ␾=0 is computed by assuming the model center is at the centroid of all points, setting all radii to values defined by halving diameters, and setting rotation values to zero. The model is fitted to the operator-selected points using ␾=0 as the initial guess for iterative minimization of the sum of squared distances from the model surfaces to the points. The resulting set of optimized parameters ␾0 is the initial set of parameters for the minimization of Eq 8. To improve robustness, we add some modifications to the procedure. We perform several passes of optimization. ˆ v, with a 3D smoothing kerFor each pass, we convolve G nel of decreasing box size, and use the optimized parameters as the starting point for the next pass. For example, a 5 ⫻ 5 ⫻ 5 followed by a 3 ⫻ 3 ⫻ 3, etc. This makes edges detectable at a greater distance from the current model surface. In each pass, the parameter search is paused when an iteration does not reduce E(␾) by at least 0.1% of its previous best value. We then perturb the current set of parameters by a random displacement less than 1% of their value and repeat until no further reduction in the objective function is observed following three consecutive searches.

1495

LAZEBNIK ET AL

EXPERIMENTAL METHODS We applied the semiautomatic segmentation method to 20 MR volumes, consisting of a contrast-enhanced (CE) T1- and T2-weighted images for each of 10 in vivo lesions. The algorithm was implemented using Matlab 6.0 (Mathworks, Natick, MA) software. Results were compared with multioperator manual segmentation. MR Imaging We used a clinical C-arm 0.2 T MR imaging system (Siemens Magnetom Open, Erlangen, Germany). The rabbits were imaged using a T1-weighted spin-echo sequence (repetition time/echo time/number of excitations [TR/TE/ NEX] 624/26/6, field of view 180 ⫻ 180 mm, matrix 256 ⫻ 256, 13 continuous slices with thickness 3 mm) and a T2-weighted spin-echo sequence (TR/TE/NEX 3362/68/8, field of view 180 ⫻ 180 mm, matrix 256 ⫻ 256, 13 continuous slices with thickness 3 mm). A 0.2 mmol/kg intravenous injection of gadolinium contrast medium (gadopentate dimeglumine: Berlex Laboratories, Wayne, NJ) was administered 5 minutes before the T1weighted image acquisition. All image volumes were approximately centered about the lesion and perpendicular to the electrode path. Rabbit Thigh Ablation Procedure The experiments were performed following a protocol approved by the Institutional Animal Research Committee. In each experiment, a New Zealand White rabbit was anesthetized, its abdomen and right thigh shaved, and its legs secured to a custom-built Plexiglas support to prevent movement of the right thigh. The animal was positioned within the gantry of the MRI system described previously. Two 8 ⫻ 12 cm wire mesh grounding pads (Radionics, Burlington, MA), coated with conductive gel (Aquasonic 100; Parker Laboratories, Orange, NJ), were placed on the rabbit’s abdomen. The thighs were placed within a 12-cm diameter multiturn solenoid receive-only coil. Subsequently, a MR-compatible 17 G copper radiofrequency electrode, with a 10-mm or 20-mm exposed tip (Radionics, Burlington, MA), was inserted percutaneously into the thigh muscle using MR image guidance. Before ablation, T2-weighted images were acquired. After the electrode tip was positioned, a lesion was induced in the muscle tissue by increasing the local temperature, using RF electric current between the electrode tip and ground pads. RF energy was applied for 2 minutes using the

1496

Academic Radiology, Vol 12, No 12, December 2005

10-mm tip or 3 minutes using the 20-mm tip by a 100 W RF generator operating at 500 kHz (model RFG-3C, Radionics). The tip of the RF electrode was maintained at 90 ⫾ 2°C. A thermistor within the electrode tip provided accurate instantaneous temperature information. Before postablation imaging, the RF electrode was removed from the thigh to prevent the appearance of a susceptibility artifact in the lesion. Subsequently, T2-weighted images were acquired followed by CE T1-weighted images. The rabbit was sacrificed approximately 30 minutes after imaging using a barbiturate overdose technique. Manual Segmentation To compare the semiautomatically fitted model to actual lesion shape, we manually segmented MR thermal lesion images. Both the inner and outer boundaries were segmented for each CE T1- and T2-weighted image volume. As a group, investigators established criteria for boundaries under supervision of a radiologist (JSL) specializing in iMRI RF ablation. For example, it was determined that boundaries should be smooth as to avoid occasional long “finger-like” protrusions from the outer boundary, probably corresponding to hemorrhaging or an edematous response. Before actual segmentation, participants segmented a training set consisting of images similar to the experimental data. Training set results were compared among all investigators and further training performed until all group members obtained consistent results. Segmentation was independently performed, over a 4-week period by two of the authors (RSL, BDW), using a freehand region of interest (ROI) tool in Analyze 3.1 (Analyze Direct, Lenexa, KS) image analysis software. Important features of the software include the ability to look at adjacent image slices simultaneously and a ROI editing tool. Persons segmenting data followed a strict protocol. All participants used the same workstation, lighting conditions, and perceptually linear display (ColorVision Optical, Rochester, NY). Window and level settings for each image set were fixed before segmentation such that the qualitative contrast between the inner and outer lesion zones was maximized. The T2-weighted images, acquired before ablation, could be viewed simultaneously on our display and were used to visualize background features such as fat and connective tissue to appropriately eliminate these structures from the segmented lesions. Each participant was blinded from the result of the other participant. In a single, 2-hour session, two to four image volumes were segmented. To minimize bias,

Academic Radiology, Vol 12, No 12, December 2005

SEMIAUTOMATIC PARAMETRIC MODEL-BASED 3D LESION

Table 1 Typical Model Parameter Values During Optimization. Units are in Millimeters and Radians. Phase 1 Indicates Starting Values, Obtained Directly from User Selected Points. Phase 2 is the Refinement of Phase 1 Parameters by Fitting, in a Least Squares Sense, the Model to the User Selected Points. Phase 3 Indicates the Final Values After Fitting the Model to the Gray-Scale Image Volume Data

CE T1- and T2-weighted images for a given rabbit were never segmented during the same day. After segmentation, the ROI boundary coordinates were exported and analyzed using Matlab. Each segmented point was assigned a 3D coordinate, in millimeter dimensions, based on the voxel dimensions. For each volume, the two independently segmented boundaries were “averaged” using a morphologic image processing technique. First, the area between a pair of boundaries for each slice was filled to create a binary image. Second, this image was thinned until only the skeleton remained. Finally, spurs were pruned from the main skeleton using further morphologic processing, leaving only the average boundary. In addition, an automatic function was used to compute the mean distance between the two operator segmentations, where closest points are assumed to correspond. If, for a given slice, the mean distance between the boundaries from the two operators differed by more than two voxel widths (1.4 mm), the slice was discussed and resegmented by both operators. However, among the 260 slices involved in this study, only 4 required repeat segmentation. The averaged set of boundaries was used for error estimation.

{u,v}. Because the distance function was well behaved, a minimum distance was obtained with very little iteration. This method is similar to the one recommended by Hart (26). For both the inner and outer surfaces separately, we obtained a set of minimum signed distances in which manually segmented points within and outside the fitted surface were assigned a negative and positive sign, respectively. We computed the mean and standard deviation of these sets of data as well as the median and quartile values of the absolute distances. We also compared the semiautomatically segmented model with one fitted directly to manually segmented boundaries. We fit the model to segmented points by minimizing the sum of squared distances from segmented boundaries to model surfaces (8), and similarly analyzed distances. This allowed estimation of error from lack of the assumed ellipsoidal geometry. By comparing the semiautomatic and the manual segmentation methods, we could thus better evaluate how much error is attributable to nonellipsoidal lesion geometry rather than failure to find gray-scale edges.

Model Fit Error Estimation To evaluate the quality of the semiautomatically fitted model, we compare it with manually segmented data. For points on the segmented boundaries, we determine the closest points on the corresponding model surface and consider these as corresponding points for error measurements. To start the process, we determined which region of the ellipsoid surface was closest to a manually segmented point by testing several {u,v} pairs (typically 20) spaced equally about the surface. We then used the best pair as the starting point for iterative minimization of the distance using a Gauss-Newton algorithm that varies

For all experiments, the algorithm converged to a solution within a few thousand calls of the objective functions and less than 10 restarts. Typical parameter evolution and corresponding error values are shown by Tables 1 and 2. Virtually identical final model parameters were obtained for the two operators’ initial parameters, resulting in almost identical parameters and error metrics. That is, the mean differences between operators for the median absolute errors were 0.09 ⫾ 0.10 mm (mean ⫾ 1 SD) and 0.12 ⫾ 0.11 mm for the inner surface, 0.07 ⫾ 0.08 mm and 0.08 ⫾ 0.11 mm for the outer surface, where values are given for CE T1- and T2-weighted images, respec-

RESULTS

1497

LAZEBNIK ET AL

Table 2 Typical Model Error Values During Optimization. Units are in Millimeters. Error is Computed by Comparing a Current Model Instance to Manually Segmented Lesion Boundaries. Values Correspond to the Same Experiment as Table 1

Figure 3. Signed distances from the automatically segmented surface to manually segmented in vivo lesion boundaries. For contrast-enhanced (a) T1- and (b) T2-weighted images, the mean and one standard deviation of signed distance errors for the inner (I) and outer (O) lesion boundaries are plotted. Signed three-dimensional distance errors are obtained from manually segmented boundary points to the nearest point on the model surface, where points inside and outside the model surface are negatively and positively valued, respectively.

tively. This suggests the algorithm converged to nearly identical parameters for both operators’ starting points. Thus, in subsequent discussion, we will report only one result per initial parameter set—that associated with the largest error in cases where there was a difference. We also note these differences are much smaller than differences between the two operators manual segmentation of lesions. Mean absolute distances between the borders segmented by the two operators were 0.38 ⫾ 0.13 mm and 0.59 ⫾ 0.20 mm for the CE T1-weighted data, 0.48 ⫾

1498

Academic Radiology, Vol 12, No 12, December 2005

0.14 mm and 0.57 ⫾ 0.16 mm for the T2-weighted data, for the inner and outer surfaces, respectively. These values provide an estimate of the uncertainty in manual lesion segmentation. To evaluate semiautomatic segmentation, we first visually examined an overlay of the segmented border on the applicable MR images. In Fig 1, we show typical overlay images, demonstrating the characteristic elliptical appearance of the inner and outer regions. We also visualize the results in 3D. In Fig 2, we qualitatively examine model fit by comparing manually segmented lesion boundaries with a surface rendering of the corresponding semiautomatically fitted model. We also quantitatively examined the semiautomatic segmentation results both in terms of signed and absolute error. In Fig 3, we examine the signed mean distances from manually segmented in vivo lesion boundaries to the semiautomatically fitted model, where points inside and outside the model ellipsoid are negative and positive, respectively. For all but two lesions (specimens 9 and 10), the mean was ⱕ0.62 mm. For most specimens, the SD was ⬍1.0 mm. We note that the electrode moved slightly during the ablation of specimen 9, resulting in a slightly banana-shaped lesion, and the lesion for specimen 10 bordered a connective tissue boundary on one side that may have acted as a barrier to current flow and heat transfer. Both lesions were less elliptical than normal; consequently, the associated signed error distances were larger and less symmetrically distributed. In Fig 4, we examine the median values of absolute distances. For reference, we also show the same values for a model directly fitted to manually segmented boundaries. For 7 of 10 lesions, the median was less than a voxel width (0.7035 mm) and less than two voxel widths (1.4 mm) for the remaining three lesions. We also note that all median values were significantly smaller than the slice thickness of 3.0 mm. Further, for most lesions, median values are approximately comparable to the previously reported uncertainty in manual segmentation. In general, the error for the semiautomatically fitted ellipsoid error was only slightly greater than that for the model fitted directly to segmented points; the mean differences between corresponding values were 0.30 ⫾ 0.16 mm and 0.17 ⫾ 0.06 mm for the CE T1-weighted data, 0.33 ⫾ 0.26 mm and 0.17 ⫾ 0.14 mm for the T2 weighted data, for the inner and outer surfaces, respectively.

Academic Radiology, Vol 12, No 12, December 2005

SEMIAUTOMATIC PARAMETRIC MODEL-BASED 3D LESION

Figure 4. Absolute distances from manually segmented in vivo lesion boundaries to the fitted model. Values for contrast-enhanced (a) T1- and (b) T2-weighted images are shown along with interquartile ranges (error bars) for both the inner (I) and outer (O) lesion boundaries. For reference, we also show the corresponding median distance to the lesion model directly fitted to the manually segmented boundaries (bars).

DISCUSSION Application of the automated method to in vivo rabbit thigh data was successful. As indicated by Fig 4, all data sets featured a median distance ⱕ1.11 mm from manually segmented boundaries. The in-plane voxel resolution was 0.70 mm and the slice thickness was 3.0 mm; hence, even the largest error in Fig 4 was about one and a half voxel widths. Qualitative inspection of slices for which error was largest suggests that they often contain areas of connective tissue or edematous fluid that infiltrated the lesion. These features distort lesion appearance and made their segmentation more difficult. For these slices, the mean and standard deviation of difference between operator-segmented boundaries were also highest. Hence, for these difficult to segment images, any error from the semiautomatic segmentation is generally within the ability of a human to segment. We also examined discrepancies between manually and semiautomatically segmented boundaries, using the signed mean distance from manu-

ally segmented points to the fitted model. We note that the model surfaces are typically interior to the manually segmented edges by approximately half a voxel width. This is likely because of the tendency of operators to demarcate edges by drawing the boundary just outside of them, whereas the semiautomatic method locates edges with subvoxel resolution. An advantage of operating on the entire volume simultaneously, rather than individual slices, is that the method does not require well-defined lesion boundaries in every slice. In general, slices in which the contrast between the inner core and outer hyperintense regions is largest are those that provide the most useful edge information. In practice, some images are difficult to manually segment because of hemorrhagic and edematous zones that obscure lesion boundaries, particularly in T2-weighted scans. For this study, operators segmented all lesions as completely as possible, such that only the inner zone of one slice in one lesion was not fully demarcated and subsequently not included in the error analysis. By assuming an ellipsoidal shape, we effectively interpolate a lesion boundary through these areas based on nearby edge information. Similarly, because the shape is constrained, the method does not require the entire lesion volume for segmentation (8), allowing for gaps between slices and thus potentially reducing imaging time. Further, unlike the semiautomatic segmentation technique, manual segmentation often requires well chosen window and level setting for the grayscale images. A distinct advantage of model-based segmentation is the ability to parameterize the model using physically meaningful quantities, such as radii, and constrain parameters using a priori information about these quantities. Given that lesion size is predictable and dependent on controlled factors such as electrode tip length, ablation duration, and power output, it is reasonable to constrain parameters to a range that only encompasses realistic lesion sizes. In our experience, this typically decreases the search algorithms convergence time. As with many optimization problems, a reasonable initial set of parameters is required. Our experience suggests that the semiautomatic technique described above, where starting parameters are computed based on 12 operator-chosen points that define lesion maximum diameters, is highly effective and requires less than 30 seconds per lesion. However, fully automated strategies are feasible. Because lesion size is a function of the exposed electrode length, ablation time, and tissue, it is possible to approximately predict lesion dimensions. Given this infor-

1499

LAZEBNIK ET AL

mation, the operator would only need to choose an approximate lesion center. Moreover, if coupled with an electrode tip tracking mechanism, the center is automatically chosen as well (27–30). Consideration of Bayesian region information may further aid initial parameter choice (31). Our method’s segmentation quality depends on two factors. First, error depends on the ability to discern continuous edges. We chose a three-point first derivative approximation to the gray-scale gradient. There are other approaches for edge detection using larger convolution kernels that simultaneously perform noise averaging. However, preliminary experiments with such approaches suggested that, for our in vivo data, the associated edge blurring was sufficiently detrimental such that a first-order approach was superior. Second, segmentation error depends on the ability of an ellipsoidal model to approximate lesion shape. That is, more ellipsoidal lesions allow a lower minimum achievable error, as shown by comparing the error of the model fitted to segmented points to that of the semiautomatic segmentation in Fig 4. Overall, the mean difference between these methods was less than half a voxel width (0.35 mm) for both boundaries and scan types, suggesting the methods are very comparable. In a previous study (8), we determined that when the model is fitted to manually segmented data, the median absolute distance from the model surface to data was ⱕ0.69 mm for the inner surface and ⱕ0.76 mm for the outer surface. These values are approximately a voxel width (0.70 mm) and much less than the slice thickness (3.0 mm). On a per lesion basis, the error due to automated edge fitting is only slightly more than error the due to the model’s approximation of lesion geometry. Thus our semiautomatic method provides a very good geometric description of a lesion generated by a stationary RF electrode in relatively homogeneous tissues, as often encountered in clinical practice. Our anecdotal experience suggests that many clinical treatment scenarios, particularly in liver, are satisfied by such lesions, and are thus appropriate for our method. However, previous studies have commented on the ability to create nonellipsoidal lesions through electrode movement during ablation, or by compounding multiple ablations (9,32). In addition, the presence of significant vasculature, heterogeneous tissue, or anatomic boundaries that impede current conduction, may result in less predictably shaped lesions. Our specific model accommodates significant lesion asymmetry by using three radii for each surface. To accommodate more extreme asymmetries, not described by our model, a

1500

Academic Radiology, Vol 12, No 12, December 2005

hierarchical segmentation approach is possible. Our method can provide a good initial estimate of lesion shape. This will serve as the initial value for segmentation of a more complex surface using other techniques as described in the introduction. Further, by manipulating the radii of a fitted ellipsoidal model, segmentation boundaries can be generated to decrease convergence time of the complex surface. There are many potential clinical applications for our method. Interior and margin volumes are easily computed from the model. Such measurements can be used in animal model studies to optimize RF application, including modifications to cooling (33), gauge and length of exposed tip (34), RF pulsing paradigms (35), application of saline to increase the size of the lesion (36), and multiprobe configurations (32). Another application is comparison of measured volumes to theory (11,37). The ability of the method to estimate good volumes in the absence of missing boundary information is especially helpful. In addition, partial volume effects should be reduced because the model surface is described analytically and cuts through voxel boundaries. By comparing estimated radii to the direction of potential tissue homogeneities such as muscle fiber orientation or the presence of a nearby bone or blood vessel, one can examine such effects. Lesions volume registration for serial studies of thermal lesions might be aided by using the model center as a correspondence point. Although we are using the method for animal model studies, it is equally applicable for clinical studies. There, one might document lesion volumes for correlation to items such as tumor regression/progression, application paradigms, tumor location, tumor type, and even patient survival. In this regard, we are investigating the relationship between the MR signal and histological response to RF ablation (12,13). The results of this study could be combined with the model to provide an interventional radiologist with a 3D visualization of cellular viability probabilities after ablation treatment. REFERENCES 1. Lamb GM, Gedroyc WM. Interventional magnetic resonance imaging. Br J Radiol 1997; 70:S81–S88. 2. Lewin JS. Interventional MR imaging: concepts, systems, and applications in neuroradiology. AJNR Am J Neuroradiol 1999; 20:735–748. 3. Jolesz FA, Nabavi A, Kikinis R. Integration of interventional MRI with computer-assisted surgery. J Magn Reson Imaging 2001; 13:69 –77. 4. Goldberg SN, Gazelle GS. Radiofrequency tissue ablation: physical principles and techniques for increasing coagulation necrosis. Hepatogastroenterology 2001; 48:359 –367. 5. Aschoff AJ, Rafie N, Jesberger JA, et al. Thermal lesion conspicuity following interstitial radiofrequency thermal tumor ablation in humans: a comparison of STIR, turbo spin-echo T2-weighted, and contrast-en-

Academic Radiology, Vol 12, No 12, December 2005

6.

7.

8.

9.

10. 11.

12.

13.

14. 15. 16. 17.

18.

19.

20.

SEMIAUTOMATIC PARAMETRIC MODEL-BASED 3D LESION

hanced T1-weighted MR images at 0.2 T. J Magn Reson Imaging 2000; 12:584 –589. Steiner P, Botnar R, Dubno B, et al. Radio-frequency-induced thermoablation: monitoring with T1-weighted and proton-frequencyshift MR imaging in an interventional 0.5-T environment. Radiology 1998; 206:803– 810. Steiner P, Botnar R, Goldberg SN, et al. Monitoring of radio frequency tissue ablation in an interventional magnetic resonance environment. Preliminary ex vivo and in vivo results. Invest Radiol 1997; 32:671– 678. Lazebnik RS, Weinberg BD, Breen MS, et al. A 3D model of lesion geometry for evaluation of MR-guided thermal ablation therapy. Acad Radiol 2002; . Chung YC, Duerk JL, Lewin JS. Generation and observation of radio frequency thermal lesion ablation for interventional magnetic resonance imaging. Invest Radiol 1997; 32:466 – 474. Goldberg SN, Gazelle GS, Dawson SL, et al. Tissue ablation with radiofrequency using multiprobe arrays. Acad Radiol 1995; 2:670 – 674. Cheng YC, Brown RW, Chung YC, et al. Calculated RF electric field and temperature distributions in RF thermal ablation: comparison with gel experiments and liver imaging. J Magn Reson Imaging 1998; 8:70 –76. Breen, M. S., Lancaster, T. L., Lazebnik, R. S., Nour S.G., Lewin, J. S., and Wilson, D. L., “Three Dimensional Method for Comparing in vivo Interventional MR Images of Thermally Ablated Tissue with Tissue Response,” Journal of Magnetic Resonance Imaging, in press 2003. Lazebnik RS, Breen MS, Fitzmaurice M, et al. Radiofrequency induced thermal lesions: sub-acute MR appearance and histological correlation. J Magn Reson Imaging. In press. McInerney T, Terzopoulos D. Deformable models in medical image analysis: a survey. Med Image Analysis 1996; 1:91–108. Kass M, Witkin A, Terzopoulos D. Snakes: active contour models. Int J Comp Vision 1987; 1:321–331. Niessen WJ. Geodesic deformable models for medical image analysis. Med Imaging IEEE Trans 1998; 17:634 – 641. Joshi S. Multiscale deformable model segmentation and statistical shape analysis using medial descriptions. Med Imaging IEEE Trans 2002; 21:538 –550. Bardinet E, Cohen LD, Ayache N. A parametric deformable model to fit unstructured 3D data. Computer Vision Image Understanding 1998; 71:39 –54. Bardinet E, Cohen E, Ayache N. Analyzing the deformation of the left ventricle of the heart with a parametric deformable model. RR-2797. 2002. O’Donnell T, Boult T, Gupta A. Global models with parametric offsets as applied to cardiac motion recovery. Proceedings CVPR 96. New York: IEEE Press; 1996, 293–299.

21. Park J, Metaxas D, Axel L. Analysis of left ventricular wall motion based on volumetric deformable models and MRI-SPAMM. Med Image Analysis 1996; 1:53–71. 22. Park J, Metaxas D, Young A, et al. Deformable models with parameter functions for cardiac motion analysis from tagged MRI data. IEEE Trans Med Imaging, 1996; 15:290 –298. 23. Sitek A. Deformable model of the heart with fiber structure. Nuclear Sci IEEE Trans 2002; 49:789 –793. 24. Turner DA, Anderson IJ, Mason JC. An algorithm for fitting an ellipsoid to data. RR9803. School of Computing and Mathematics, University of Huddersfield, UK; 1999. 25. Nelder JA, Mead, R. A simplex method for function minimization. Computer J 1965; 7:308 –313. 26. Hart JC. Distance to an ellipsoid. In: Heckbert PS, ed. Graphics gems IV. Boston: AP Professional, 1994; 113–119. 27. Patel KC, Duerk JL, Zhang Q, et al. Methods for providing probe position and temperature information on MR images during interventional procedures. IEEE Trans Med Imaging 1998; 17:794 – 802. 28. Wendt M, Wacker FK. Visualization, tracking, and navigation of instruments for magnetic resonance imaging-guided endovascular procedures. Top Magn Reson Imaging 2000; 11:163–172. 29. Zhang Q, Wendt M, Aschoff AJ, et al. Active MR guidance of interventional devices with target-navigation. Magn Reson Med 2000; 44:56 – 65. 30. Zhang Q, Wendt M, Aschoff AJ, et al. A multielement RF coil for MRI guidance of interventional devices. J Magn Reson Imaging 2001; 14:56 – 62. 31. LaValle SM, Hutchinson SA. A Bayesian segmentation methodology for parametric image models. IEEE Trans Pattern Analysis Machine Intell 1995; 17:211–217. 32. Goldberg SN, Gazelle GS, Dawson SL, et al. Tissue ablation with radiofrequency using multiprobe arrays. Acad Radio 1995; 2:670 – 674. 33. Lorentzen T. A cooled needle electrode for radiofrequency tissue ablation: thermodynamic aspects of improved performance compared with conventional needle design. Acad Radiol 1996; 3:556 –563. 34. Goldberg SN, Gazelle GS, Dawson SL, et al. Tissue ablation with radiofrequency: effect of probe size, gauge, duration, and temperature on lesion volume. Acad Radiol 1995; 2:399 – 404. 35. Goldberg SN, Stein MC, Gazelle GS, et al. Percutaneous radiofrequency tissue ablation: optimization of pulsed-radiofrequency technique to increase coagulation necrosis. J Vasc Interv Radiol 1999; 10:907–916. 36. Boehm T, Malich A, Goldberg SN, et al. Radio-frequency tumor ablation: internally cooled electrode versus saline-enhanced technique in an aggressive rabbit tumor model. Radiology 2002; 222:805– 813. 37. Johnson PC, Saidel GM. Thermal model for fast simulation during MRI guidance of RF tumor ablation. Ann Biomed Eng 2002; .

1501