Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2838620

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. A METHOD TO DETERMINE 3D SURFACES OVER TIME
- III. EXPERIMENTS AND RESULTS
- IV. CONCLUSION
- REFERENCES

Authors

Related links

J Signal Process Syst. Author manuscript; available in PMC 2010 April 1.

Published in final edited form as:

J Signal Process Syst. 2009 April 1; 55(1-3): 139–156.

doi: 10.1007/s11265-008-0219-1PMCID: PMC2838620

NIHMSID: NIHMS171406

Peter C. Tay, Member, IEEE, Bing Li, Member, IEEE, Chris D. Garson, Student, IEEE,, Scott T. Acton, Senior Member, IEEE, and John A. Hossack, Senior Member, IEEE

Peter C. Tay, Dept. of Electrical and Computer Engineering Technology, Western Carolina University, Cullowhee, NC 28723 USA;

Peter C. Tay: ude.ucw.liame@yatp; Bing Li: ude.ainigriv@n3lb; Chris D. Garson: ude.ainigriv@d6gdc; Scott T. Acton: ude.ainigriv@notca; John A. Hossack: ude.ainigriv@kcassoh

A method to perform 4D (3D over time) segmentation of the left ventricle of a mouse heart using a set of B mode cine slices acquired *in vivo* from a series of short axis scans is described. We incorporate previously suggested methods such as temporal propagation, the gradient vector flow active surface, superquadric models, *etc.* into our proposed 4D segmentation of the left ventricle. The contributions of this paper are incorporation of a novel despeckling method and the use of locally fitted superellipsoid models to provide a better initialization for the active surface segmentation algorithm. Average distances of the improved surface segmentation to a manually segmented surface throughout the entire cardiac cycle and cross-sectional contours are provided to demonstrate the improvements produced by the proposed 4D segmentation.

Modern medical imaging systems provide researchers and physicians the ability to view inside living beings and observe organs *in vivo*. This allows observation of useful anatomical and functional information. The ability to observe organs *in vivo* is more desirable when the data used to produce images can be safely, noninvasively, accurately, and cost effectively collected. These and other advantages (portability, real-time or near real-time imaging, *etc.*) make the ultrasound imaging modality an attractive choice for use in cardiology.

An important task for medical researchers and physicians is to be able to distinguish boundaries and segment regions contained within these boundaries. This allows functional information such as the volume of the left ventricle to be calculated over time, which would yield ejection fraction, cardiac output, *etc.* This paper describes the algorithmic steps needed to construct the left ventricle (LV) of a mouse heart from *in vivo* ultrasound scans. This task requires that 3D data of the LV be acquired *in vivo* throughout the systolic cycle. There are numerous ways to create 3D data from 2D slices [1], [2]. A typical method to create 3D data and the one that is used in this paper is to stack 2D slice images along the elevation axis as described by Fenster *et al.* in [2]. This method to create 3D data requires that slices at different elevations be acquired simultaneously or at corresponding times when the cardiac cycles. The simultaneous acquisition of 3D data would require the use of a 2D transducer. Although this technology promises greater functionality, it is still lacking in spatial resolution and is too expensive when compared to a single element or linear array transducer. The scans reported in this paper are the elevational slices that are acquired over different cardiac cycles. The acquisitions are gated to occur at corresponding times during the cardiac cycle.

Previously proposed 3D segmentation methods can be categorized as model fitting methods, model deformation methods, or some combination of both. A compilation of recently proposed 3D cardiac models using a variety of medical imaging modalities have been surveyed in the literature [3], [4]. A comprehensive survey of recently proposed deformable model methods using a variety of medical imaging modalities is provided by McInerney and Terzopoulos in [5]. The use of a globally optimized model that is fitted to the data then locally deforming the model was proposed by Bardinet *et al.* in [6], [7] and also by Montagnat *et al.* in [8] and Chen *et al.* in [9].

The construction of the left and right ventricle chambers using cine magnetic resonance imaging (MRI) has been proposed by Goshtasby and Turner in [10] and Kaus *et al.* in [11]. A level set method that uses *a priori* shape information to track the motion of the LV using MRI and ultrasound imaging is reported by Paragios in [12]. LV segmentation and/or volume estimation methods specific to echocardiography have been reported by Corsi *et al.* [13], Coppini *et al.* [14], Gustavsson *et al.* [15], McCann *et al.* [1], McPherson *et al.* [16], Strickels and Wann [17], and Geiser *et al.* [18].

An automatic or semi-automatic method to distinguish surfaces of interest require edges and borders of objects to be well defined. This can be problematic when working with ultrasound, since ultrasound images can be ambiguously corrupted by speckle and image dropout^{1} commonly occurs. These problems with ultrasound images are discussed by Chu *et al.* in [19]. Speckle can cause ultrasound images to appear mottled, true edges to be indistinguishable from other gradient causing factors, and false gradients to be detected in homogeneous regions. This can result in unusable segmentation when using methods that depends on well defined borders and edges such as with most, if not all, active contour/surface methods [20]–[25]. There are a multitude of proposed methods to reduce speckle [26]–[39]. Application of the gradient vector flow (GVF) active surface (AS) [23] to segment 3D data composed of despeckled 2D slices can provide accurate surfaces except in areas of weak edges and where image dropout occurs. Areas of weak edges and image dropout are evident in Fig. 1, which is a typical image from a single element high frequency transducer attached to the VisualSonics Inc. (Toronto, Canada) Vevo 770 ultrasound system. Parts of the manually drawn endocardial (red) and epicardial (green) contours of Fig. 1 are ambiguously defined. An autonomous 3D segmentation method could erroneously places both surfaces in close proximity to each other or the endocardial surface of the LV extending past the epicardial surface. These errors generally occur in areas where the endocardial border is weak (relative to the epicardial borders) or where the data is ambiguous perhaps due to image dropout. Our 4D (3D over time) segmentation method performs a 3D GVF AS segmentation at each time instant where the initial approximation is achieved through temporal propagation [5]. Temporal propagation uses the final surface from the previous time instant as the initial approximation in the current time instant, which is required by the AS segmentation algorithm. The advantage of temporal propagation is that user interaction is only required for the initial approximation at the first time instant. Thus, this approach is less labor intensive and more objective than requiring user input at every time instances. A problem with temporal propagation is that part(s) of the endocardial and epicardial surfaces where the true borders are ambiguous tend to converge towards each other. The propagation of these errors will cause other parts of both surfaces to converge towards each other. Thus, the possibility of large regions of both surfaces erroneously occupying the same location exists, if left unchecked.

A typical ultrasound cross section of LV with manually segmented endocardial (red) and epicardial (green) contours.

Our solution to the problems of erroneously placed surfaces in areas of weak boundaries, image dropout, and part(s) of the epicardial and endocardial surfaces converging towards each other is to assume the epicardial surface accurately depicts the true epicardial border and locally fit a superellipsoid model [6] to the outer surface. The locally fitted superellipsoid yields a closed form equation, which we use to determine areas where the outer and inner surfaces are in close proximity or where the inner surface extends past the outer surface. The locally fitted superellipsoid is determined by minimizing an objective function. This minimization is determined by a mesh search algorithm, which is capable of finding a global minimum and can avoid being stuck at a local minimum^{2}. The locally fitted superellipsoids are used to adjust the time propagated initial approximation of the inner surface. Once the inner surface is sufficiently contained within the outer surface, the application of the GVF active surface segmentation produces a smooth representation of the endocardial border.

Although more novel methods for segmentation such as the level set methods [40] could be used and incorporation of such methods is an interesting subject for future research, the segmentation method implemented in the proposed algorithm described in this section relies on the GVF AS. A robust implementation of the GVF AS requires edges and borders to be well defined.

Unfortunately, ultrasound images are significantly hampered by a phenomena called speckle. This type of noise make autonomously delineating edges and borders difficult and false edges to be present. As a pre-segmentation step the ultrasound images must be despeckled. A novel despeckling technique that is shown in [39] to provide excellent contrast while smoothing homogeneous regions is used to reduce the effects of ultrasound speckle. The despeckling filter describe in [39] is aptly name the squeeze box filter (SBF) because of the compressing nature of the filtering method. The results of SBF applied to the original six slices of Fig. 2(a) are shown in Fig. 2(b). The filtered images in Fig. 2(b) show that the speckle is smoothed in homogeneous regions and important edge features corresponding with the myocardium and papillary muscles are preserved in each short axis slice.

(a) The original short axis slices acquired *in vivo*.. (b) Speckle reduced slices of the LV of an mouse heart.

The proposed algorithm to reconstruct the mouse LV consists of the following six sequential processes.

- User initialization consisting of defining the endocardial and epicardial borders for all the acquired slices (six in this paper) at time
*t*= 0. - The GVF AS is used to determine approximations of the endocardial and epicardial surfaces.
- The time variable
*t*is incremented. The epicardial surface from the previous time instance*t*− 1 is used as the initial approximation to the GVF AS, which determines the current epicardial surface. - The epicardial surface is modeled by locally fitted superellipsoids.
- The locally fitted superellipsoids are used to define the initial approximation of the endocardial surface, which is the modified propagated endocardial surface from the previous time instant
*t*− 1. - The GVF AS uses the initial approximation to determine the current endocardial border, then the algorithm is repeated from step three.

The flow diagram of our 3D segmentation over time is given in Fig. 3 and each process is explained in the following sections.

A manual slice by slice segmentation is used to initialize the endocardial and epicardial surfaces of the 3D data set at *t* = 0. A dense set of points in 3D space representing the initial surface could be defined by interpolating the points of the slice contours in the elevational direction, as suggested in [41]. The downside of this method is that the number of contour points at different slices must be the same. Since the lengths of the cross-sectional contours at various slices can vary significantly, this method could cause large distances between contour points. Instead, a 3D object is created from a connected component analysis. The connected component in each slice is determined by the manually defined contour of each slice.

The interior of each manually placed closed contour is segmented by a 8-connective component analysis. The 8-connective component analysis takes the binary contour image, which is one if the pixel is on the contour and zero otherwise, and creates a binary image that is one if the pixel is within the interior of the contour and zero otherwise. The closed contour is densely sampled such that adjacent points on the contour are one pixel apart. The 8-connective component analysis creates a binary image by scanning in some fixed direction for example each row from left to right. When the first nonzero pixel is encountered, then the subsequent pixels in that direction in this example each row moving left to right become one. The same is applied when scanning from right to left, top to bottom, bottom to top, top left to lower right, top right to lower left, lower left to top right, and lower right to top left. The binary images of each direction produced by this technique is shown in Fig. 4. A pointwise multiplication, equivalently a logical “and” operation, of the eight binary images in Fig. 4 results in a connected component of the interior defined by the closed contour. Fig. 5 shows the blob created from a pixelwise multiplication (logical “and”) of all eight connect component binary images in Fig. 4.

An application of the 8-connective component analysis to each of the six user defined contours at frames indexed by *t* = 0 produces a 3D object, which approximates the interior volume of the LV. An example of the 3D connected component of all six slices is shown in Fig. 6(a). Each connected component slice is downsampled and elevationally interpolated in exactly same manner as the interpolated despeckled 3D data to produce the 3D object used as the initial approximation at time *t* = 0 (shown in Fig. 6(b)). An isosurface of this object is shown in Fig. 7(a).

The user specified initial approximation shown by the isosurface in Fig. 7(a) suffers from striations or banding effect [2], [5]. The application of AS was used by Taubin in [42] as a method to smooth a surface while being constrained by three dimensional borders. In the same spirit the GVF AS is applied to the vertices of the striated approximation to attain a surface that varies smoothly. The effects of the GVF AS applied to the vertices of the striated surface of Fig. 7(a) are shown in Fig. 7(b).

Up to this point in the algorithm, the endocardial and epicardial surfaces are defined in the same manner. Fig. 8 show the side, the top, and the bird’s eye view of the endocardial (red) and epicardial (green) surfaces super imposed on the same image. These 3D surface renderings give an approximate reconstruction of the myocardial borders of the LV. It should be noted that at this early stage both surfaces are well separated and distinct. Segmentation results of this nature are desired, since it characterizes the true anatomy of the myocardial walls. Repetition of the user involved segmentation process to this point for the remaining 4D sequences is labor intensive and strongly discouraged.

The 3D segmentation of the epicardial border at each time increment *t* is independent of the 3D endocardial segmentation. The segmentation of the epicardial surface is needed so that corrections to the time propagated initial approximation of the endocardial surface can be made. The epicardial border in our 3D data is assumed to be more reliable and easier to detect than the endocardial border. This assumption is due in part to the shape of the epicardial border. The epicardial border can be observed to be more regular in shape and more smoothly varying (the concepts of regular and smooth are both defined in the Gestalt [43] sense) than the endocardial border, which includes the protruding papillary muscles. From experimental observation it was discovered that due to its regularity and smoothness the epicardial border is more reliably segmented by the GVF AS. Thus, based on observation the epicardial surfaces produced by the GVF AS are more accurately defined.

The resulting epicardial surface from the previous time increment *t* − 1 is propagated as the initial approximation. The GVF AS is applied to determine the epicardial surface at the current time increment *t* ≠ 0. At *t* = 0 the initial approximation of the epicardial surface is determined from user inputs. The final epicardial surface at *t* = 0 is the result of the GVF AS applied to the user dependent initial approximation of this surface.

Finding the best fitting superellipsoid to a given set of point is an minimization problem that is determined with a search algorithm. An exhaustive search would guarantee an optimal solution, but is prohibited by high computational cost. Greedy methods like a gradient descent algorithm may converge only to a locally optimal solution with no chance of attaining a globally optimal solution. Methods such as simulated annealing, genetic algorithms, or pattern searches compromising computation time with the possibility of finding a globally optimal solution. The need for computational efficiency and accuracy motivates our use of a pattern search method to determine a set of local superellipsoids that is fitted to the epicardial surface.^{3} This set of local superellipsoids will yield a close form equation in which to evaluate discrete samples of the approximated endocardial surface.

Let ${D}_{\mathit{\text{EP}}}^{t}$ be a finite set of vertices that is a discrete sampling of points on the estimated epicardial surface at time *t*. A set of samples of the epicardial surface at elevation *z* is denoted as

$${D}_{\mathit{\text{EP}},z}^{t}=\{(x,y,z)\phantom{\rule{thinmathspace}{0ex}}\in \phantom{\rule{thinmathspace}{0ex}}{D}_{\mathit{\text{EP}}}^{t}\phantom{\rule{thinmathspace}{0ex}}\text{for}\phantom{\rule{thinmathspace}{0ex}}\text{some}\phantom{\rule{thinmathspace}{0ex}}\text{fixed}\phantom{\rule{thinmathspace}{0ex}}z\}.$$

(1)

Clearly, the union of all such sets defined in equation (1) provides a discrete sampling of the approximated epicardial surface, that is ${D}_{\mathit{\text{EP}}}^{t}={\displaystyle \underset{\forall z\in \mathbb{Z}}{\cup}{D}_{\mathit{\text{EP}},z}^{t}}$. A pattern search method is used to find the parameters of the locally optimal superellipsoid, which is denoted as the vector ${\mathbf{v}}_{z}^{t}=({\widehat{x}}_{c},{\widehat{y}}_{c},z,{\widehat{r}}_{x},{\widehat{r}}_{y},1,{\widehat{\epsilon}}_{1},{\widehat{\epsilon}}_{2})$, so that the energy function used in [6], [7] defined in equation (2)

$${E}_{\mathit{\text{SE}}}({\mathbf{v}}_{z}^{t})={\displaystyle \sum _{\forall (x,y,z)\in {D}_{\mathit{\text{EP}},z}^{t}}{({S}_{{\mathbf{v}}_{z}^{t}}(x,y,z)-1)}^{2}}$$

(2)

where

$${S}_{{\mathbf{v}}_{z}^{t}}\left(x,y,z\right)={\left\{{\left[{\left(\frac{x-{x}_{c}}{{r}_{x}}\right)}^{\frac{2}{{\epsilon}_{2}}}+{\left(\frac{y-{y}_{c}}{{r}_{y}}\right)}^{\frac{2}{{\epsilon}_{2}}}\right]}^{\frac{{\epsilon}_{2}}{{\epsilon}_{1}}}+{\left(\frac{z-{z}_{c}}{{r}_{z}}\right)}^{\frac{2}{{\epsilon}_{2}}}\right\}}^{\frac{{\epsilon}_{1}}{2}}.$$

(3)

is minimized.

The method to determine a close form approximation for the epicardial surface represented by a set of discrete 3D points when the time index *t* is equal to zero is slightly different than for other time instances. When *t* = 0 and given a significant region of the epicardial surface is between *z _{min}* and

- If
*z*=*z*, then a user defined v_{min}_{0}where the superellipsoid*S*_{v0}is required as an initial model approximation of ${D}_{\mathit{\text{EP}},{z}_{\mathit{\text{min}}}}^{0}$ . Otherwise, ${\mathbf{v}}_{0}={\mathbf{v}}_{z-1}^{0}$. - A set of twelve mesh points at distance λ > 0 about
**v**_{0}are defined as$$\begin{array}{cc}{M}_{k}=\hfill & \{{\mathbf{v}}_{k}\pm (\mathrm{\lambda},0,0,0,0,0,0,0),{\mathbf{v}}_{k}\pm (0,\mathrm{\lambda},0,0,0,0,0,0),{\mathbf{v}}_{k}\pm (0,0,0,\mathrm{\lambda},0,0,0,0),\hfill \\ \hfill & {\mathbf{v}}_{k}\pm (0,0,0,0,\mathrm{\lambda},0,0,0),{\mathbf{v}}_{k}\pm (0,0,0,0,0,0,\mathrm{\lambda},0),{\mathbf{v}}_{k}\pm (0,0,0,0,0,0,0,\mathrm{\lambda})\}.\hfill \end{array}$$(4) - If there is a
*M*such that_{k}*S*_{}<*S*_{v}for all**v***M*_{0}{**v**} and all $({x}_{i},{y}_{i},z)\phantom{\rule{thinmathspace}{0ex}}\in \phantom{\rule{thinmathspace}{0ex}}{D}_{\mathit{\text{EP}},z}^{0}$, then_{k}**v**_{k+1}= and λ is doubled that is λ ← 2λ. The process is iterated from step two and*k*is incremented by one. Otherwise**v**_{k+1}=**v**and λ is halved that is $\mathrm{\lambda}\leftarrow \frac{\lambda}{2}$._{k} - If λ >
*T*for some predefined threshold*T*, then the search is iterated from step two and*k*is incremented by one. Otherwise, the search is terminated and ${\mathbf{v}}_{z}^{0}={\mathbf{v}}_{k}$.

When *t* is not equal to zero, a model composed of optimal locally fitted superellipsoids ${S}_{{\mathbf{v}}_{z}^{t}}(\xb7)$ is determined similarly by an iterative pattern search method that scan within a varying mesh in the following manner:

- The search is initialized by setting ${\mathbf{v}}_{0}={\mathbf{v}}_{z}^{t-1}$.
- For a vector
**v**_{k}^{8}, the set of twelve mesh points at distance λ > 0 about**v**is denoted as_{k}*M*and defined in equation (4)._{k} - If there is a
*M*such that_{k}*S*_{}<*S*_{v}for all**v***M*_{0}{**v**_{0}} and all $({x}_{i},{y}_{i},z)\phantom{\rule{thinmathspace}{0ex}}\in \phantom{\rule{thinmathspace}{0ex}}{D}_{\mathit{\text{EP}},z}^{t}$, then**v**_{k+1}= and λ is doubled that is λ ← 2λ. The process is iterated from step two and*k*is incremented by one. Otherwise**v**_{k+1}=**v**and λ is halved that is $\mathrm{\lambda}\leftarrow \frac{\mathrm{\lambda}}{2}$._{k} - If λ >
*t*for some predefined threshold*t*, then the search is iterated from step one and*k*is incremented by one. Otherwise, the search is terminated and ${\mathbf{v}}_{z}^{t}={\mathbf{v}}_{k}$.

The flow diagram of the superellipsoidal model fitting search is illustrated in Fig. 9.

The result of this search is a superellipsoid ${S}_{{\mathbf{v}}_{z}^{t}}(\xb7)$ that is well fitted to the set of local 3D points ${D}_{\mathit{\text{EP}},z}^{t}$. After application of this search for all *z* , a global closed form model *S ^{t}*(

$${S}^{t}(x,y,z)={S}_{{\mathbf{v}}_{z}}^{t}(x,y,z)$$

(5)

for all $(x,y,z)\phantom{\rule{thinmathspace}{0ex}}\in \phantom{\rule{thinmathspace}{0ex}}{D}_{\mathit{\text{EP}}}^{t}$. The optimal locally fitted superellipsoids at time *t* = 0 and at time *t* = 25 for elevations *z* = 25, 35, 45 are shown super imposed on the epicardial surface in Figs. 10(a)–10(c) and Figs. 10(d)–10(f), *resp*.

The previous set 3D discrete points of endocardial surface denoted as ${D}_{\mathit{\text{EN}}}^{t-1}$ is temporally propagated and is used as the initial approximation of the GVF AS algorithm, which is used to determine the final endocardial surface of the current time sample. This set of endocardial points is denoted as ${D}_{\mathit{\text{EN}}}^{t}$. A downfall of temporal propagation is that portions of the endocardial and epicardial surfaces converge towards each other in regions where the data is ambiguous (possibly due to image dropout, shadowing, or other artifacts). This is evident when the epicardial and endocardial surfaces are superimposed on each other as in Fig. 11(a). The red surface shown in Fig. 11(a) reflects the approximated endocardial border at time *t* = 4. The green surface shown in Fig. 11(a) is the approximate epicardial border at *t* = 5. Both surfaces were determined by the GVF AS algorithm where the initial approximation was determined by temporally propagating the final GVF AS from the previous time sample *t* = 4. The estimates of the two surfaces are reasonable in that the endocardial surface (red) is mostly contained well within the epicardial surface (green). It is shown in Fig. 11(a) a region where the endocardial surface (red) resides erroneously close to the epicardial surface (green). Except at the top of the LV, the endocardial and epicardial borders should be clearly separated. The top slice corresponds to the cross-section of the LV closest to the base of the heart. The endocardial and epicardial connected components contained in this slice is taken to be the top surface of the 3D endocardial and epicardial surfaces, *resp*. This results in the red connected component seen in the surface renderings (red) shown in Fig. 11. Likewise, the connected components in the lowest elevation closest to the apex of the heart are take to be the bottom surface of the LV. The inaccuracies caused by the overlap of endocardial and epicardial surfaces at the highest and lowest elevational slices are not alleviated by the proposed segmentation algorithm. Rather, our segmentation algorithm proposes to contain the endocardial surfaces within the epicardial surfaces between the lowest and highest elevations during the entire systolic cycle.

Endocardial and epicardial surfaces at *t* = 5 (a) without and (b) with endocardial surface correction.

In Fig. 11(a) the surface approximation of the endocardial, epicardial, or both surfaces in a region between the base and apex is clearly incorrect, in that the endocardial surface of the LV is in close proximity to the epicardial surface. This is problematic not only in the visual rendering of the two surfaces but estimates of the LV based on such a segmentation will produce erroneous results, which may lead a cardiologist to incorrectly diagnose a thin area of the myocardium or other abnormalities. If left uncorrected, then the regions where the two surfaces coincide or where the endocardial surface is erroneously close to or possibly extends past the epicardial surface will be temporally propagated and may increase over time.

This problem is alleviated by using the epicardial model to correct the temporally propagated initial approximation of the endocardial surfaces. Evaluation of each point of the endocardial initial approximation with the close form model of the epicardial surface *S ^{t}*(·) defined in equation (5) allows us to determine which points are erroneously close to or extend beyond the epicardial surface. Additionally, the epicardial model

Let 0 < ε ≤ 1 be pre-defined and ${D}_{\mathit{\text{EN}}}^{t-1}$ be a set of 3D discrete points taken from the final GVF AS segmentation of the endocardial border of the previous time *t* − 1, and *S ^{t}*(·) be the locally fitted superellipsoid model of the epicardial surface defined by the set of points on the estimated epicardial surface ${D}_{\mathit{\text{EP}}}^{t}$. A new set of endocardial surface points to be used as discrete samples of the initial approximation required by the GVF AS algorithm is denoted as ${\widehat{D}}_{\mathit{\text{EN}}}^{t}$
. The set ${\widehat{D}}_{\mathit{\text{EN}}}^{t}$ is determined by evaluating all the points in ${D}_{\mathit{\text{EN}}}^{t}$ with the model

Fig. 11(b) shows the corrected initial approximation (prior to applying the GVF AS) at *t* = 5, of the temporally propagated surface ${\widehat{D}}_{\mathit{\text{EN}}}^{4}$. A comparison of the initial approximation of endocardial (red), which was temporally propagated, and the current epicardial (green) surfaces in Fig. 11(a) exhibits a region where the two surfaces are erroneously close to each other. Fig. 11(b) show the corrected initial approximation of the endocardial surface with the epicardial surface. It is clearly evident that the two surfaces shown in Fig. 11(b) are well separated.

The corrected temporally propagated initial approximation of the endocardial surface shown in Fig. 11(b) cannot be expected to accurately follow the true endocardial border. Additionally, moving points toward the centroid of the locally fitted superellipsoid cause the surface to be unrealistically rough. Since AS have been demonstrated to accurately find borders [22] and is an excellent method to smooth surfaces [42], the result of a few iterations of the GVF AS algorithm to the points of the corrected initial approximation contained in set ${\widehat{D}}_{\mathit{\text{EN}}}^{t}$ is used to determine the final endocardial surface at time *t*. Fig. 11(a) shows the resulting endocardial surface (red) after applying the GVF AS to the non-corrected temporally propagated endocardial surface shown as the red surface in Fig. 11(a). It is evident on the middle left side that the result produced by the GVF AS of endocardial (red) and the epicardial (green) surfaces are erroneously close in proximity. This error is avoided when we use the corrected initial approximation shown as the red surface in Fig. 11(b). The result of applying the GVF AS using the corrected initial approximation is shown as the red surface in Fig. 11(c) along with the epicardial surface (green). A comparison of the surface renderings in Fig. 11(c) with the surface renderings in Fig. 11(a) shows the benefits of correcting the initial approximation of the temporally propagated endocardial surface with the optimal locally fitted superellipsoid model, in which region(s) where the two surfaces coincide or are in close proximity are separated while other regions are not affected by this correction. Additionally, application of the GVF AS to the corrected initial approximation provides a smooth approximation of the endocardial border of the LV.

This section describes the experiments performed on three mouse heart scans and provides the results of our experimentation. The mouse heart scans are acquired from a VisualSonics Vevo 770 ultrasound scanner operating at 25 MHz using an single element mechanical sweeping RMV-707B scan head. This system is electrocardiographically gated to capture ultrasound data between two successive R-waves. This echocardiographic system captures a sequence of images beginning at approximately the same time during the cardiac cycle, near end diastole^{4}. Our acquisition technique consists of attaching a transducer to a calibrated mechanical arm, starting from the base of the heart, and acquiring a sequence of short axis slices at a fixed elevation for one entire cardiac cycle. The transducer is moved one millimeter towards the apex along the elevational axis. Another sequence of images over one cardiac cycle is acquired that is parallel to the previous elevational slice plane. This process is repeated for six elevations. This scanning method produces a sequence that consists of approximately 100 frames of six elevation short axis movies. This acquisition system simulates a frame rate of approximately 1000 frames (images) per second. The description of the simulated frame rate of this proprietary system is given by Chérin *et al.* in [45]. Since the mouse heart rates in our experiments are approximately 500 cycles per minute, a high frame rate is needed to promote smooth motion. It should be noted that the simulated frame rate is fully performed in the VisualSonics Vevo 770 scanner and beyond the control of the authors. The frames from the all the sequences are despeckled using the SBF method describe in [39]. The 3D data is composed of elevationally interpolating the six short axis despeckled slices to 80 slices. This elevational interpolation factor was chosen to approximately match the axial and lateral resolution.

The frames indexed by *t* = 0 from each of the six elevation are shown Fig. 2(a). The corresponding despeckled frames are shown in Fig. 2(b). Fig. 12(a) shows the slices elevationally stacked. Fig. 12(b) show an example of the 3D data composed of despeckled and elevationally interpolated slices at time *t* = 0. The 4D data of despeckled and elevationally interpolated slices for all 100 time instances is the input of our proposed 4D segmentation algorithm.

(a) 3D matrix composed of short axis slices. (b) The 3D matrix created from interpolating the SBF slices.

In this comparison of the model corrected segmentation method with the non-corrected method, all the GVF AS parameters are fixed and identical for both the non-corrected and the corrected segmentation algorithms. The model correction parameter ε is set at 0.75. The same user supplied initialization at *t* = 0 is used for both segmentation methods.

The average distances of each segmentation to the manually segmented surface of each mouse left ventricle are shown in Fig. 13 in units of voxels. The average distance of each point on model corrected segmented surface to the manually determined surface at each time instance is plotted with the black line in each graph of Fig. 13. The red line in each graph of Fig. 13 is the average distance of each point of the non-corrected segmented surface. In all three cases the model corrected surfaces are closer to the manually segmented surface except for a few instances with mouse 1.

The average distances of the model corrected (red) and the non-corrected (blue) surfaces to the manually segmented surface over one cardiac cycle using (a) mouse 1, (b) mouse 2, and (c) mouse 3 scans.

Cross-sectional slice contours at three elevational frames corresponding to two data captured before and after end systole, *t* = 30 and 70, *resp*. are shown to provide a visually meaningful validation to the robustness of incorporating the model correction method. Slices at different three elevations are taken from the final segmentation results. These cross-sectional slices are illustrated in Fig. 14 using the surface attained with mouse 1.

Final endocardial surface of mouse 1 using model correction to the initial approximation at frame *t* = 30 and elevations 60, 30, and 1.

Fig. 15 and Fig. 16 show the cross-sectional contours at elevations 60 (closest to the base), 30, and 1 (closest to the apex) of the manual (gold), the model corrected (red), and non-corrected (blue) segmentation of the endocardial surfaces at frames *t* = 30 and *t* = 70, *resp*. It can be seen that in all but two cases the model corrected endocardial contours more closely resemble the manually segmented contours. The cross-sectional contours of the LV segmentation using scans of mouse 1, elevation 1, frames indexed by *t* = 30, are shown in Fig. 15(g). In this image, it can be seen that the contour from the model corrected segmentation (red) produces a results that better resembles the manually segmented contour (gold). Additionally, the model corrected contour is smoother than the cross sectional contour produced by the segmentation without correction. The contours of the LV segmentation using the scans of mouse 1, elevation 30, frames indexed by *t* = 70, are shown in Fig. 16(d), the model corrected and non-corrected contours are equally different from the manually segmented contour. Again, the contour from the model corrected surface is much smoother than the contour from the non-corrected surface. Thus, in all but two cases the model corrected segmentation produces cross sectional contours that more accurately resemble manual segmentation. In the two cases where the two segmentation methods were equally unlike the manually segmented contour, the model correction produced a much smoother contour than the contour produced by omitting this correction.

The manually segmented (gold), the non-corrected (blue), and the model corrected (red) contours at frame *t* = 30.

Segmentation of the LV to visualize anatomic information and to calculate functional information is an important and difficult task. This paper provide a 4D segmentation method using ultrasound short axis slices, which are despeckled and downsampled. The despeckled and downsampled slices are elevationally interpolated to form 3D data at each time increment or instance. The time varying 3D data is used to segment the LV. The proposed method requires manual contour segmentation at the beginning of the 3D sequence to initialize the six short axis slices. This initialization determines the initial endocardial and epicardial surfaces. A starting superellipsoid for a slice at *t* = 0 is also required by the user, which initializes the search for locally fitted superellipsoids that provide a model of the time varying epicardial surfaces. These locally fitted superellipsoids allow us to correct initial approximation errors due to temporally propagating the final surface from the previous time to be used as the initial approximation of the endocardial surface in the current time. It was observed that temporal propagation causes parts of the endocardial and epicardial surfaces to converge to the same regions. An optimally locally fitted superellipsoidal model of the epicardial surface allowed corrections to the time propagated initial approximation of the endocardial surface to be made. The application of the GVF AS to the corrected initial approximation of the endocardial surface results in a smooth surface that is well separated from the epicardial surface.

The average distance to the manually segmented surface is decreased when applying model correction to the temporally propagated initial surface of each time instance. This supports an increase in the robustness of our proposed model corrected 4D segmentation method. Cross-section of the segmented surface at two different time instances at various elevations for three different mouse scans show that model corrected segmentation provide contours that better reflect manually segmented contours than contours produced by non-corrected surfaces.

This work was supported in part by NIH NIBIB grant EB001826, US Army CDMRP grant (W81XWH-04-1-0240), and NIH NCRR RR022582 (for purchase of VisualSonics Vevo 770 scanner).

**Peter C. Tay** is an Assistant Professor at Western Carolina University in Cullowhee, North Carolina, USA. He received a B.S., M.A. both in Mathematics and a Ph.D. in Electrical and Computer Engineering from the University of Oklahoma.

**Bing Li** received the B.S. degree in Electrical and Computer Engineering from Peking (Beijing) University, Beijing, China in 2003. He received the M.S. and Ph.D. degrees in Electrical Engineering from the University of Virginia, Charlottesville, VA in 2005 and 2007, respectively. He is working on mask inspection systems in the RAPID division for KLA-Tencor Corporation. His research interests include image analysis, segmentation, active models, registration, biometrics, and target tracking.

**Chris D. Garson** is a Ph.D. student in the Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA. He recieved a B.S. and M.S. in Biomedical Engineering from Duke University.

**Scott T. Acton** is a Professor in the Dept. of Electrical and Computer Engineering and the Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA. He recieved a B.S. in Electrical Engineering from Virginia Tech University, a M.S. and a Ph.D. in Electrical and Computer Engineering from the University of Texas.

**John A. Hossack** is an Associate Professor with the Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA. He recieved a B. Eng. and a Ph.D. from University of Strathclyde, UK.

^{1}Area in the image where features or structure exist but were not captured in the image.

^{2}Although the speed and robustness of the algorithm could be improved by establishing a optimal search method, the optimality of the search method is beyond the scope of this paper.

^{3}Determination of which search algorithm provides the best tradeoff between computational speed and accuracy is beyond the scope of this paper.

^{4}The point during the cardiac cycle when the myocardium is most relaxed and the volume of the LV is at its maximum. See [44] for a concise definition.

Peter C. Tay, Dept. of Electrical and Computer Engineering Technology, Western Carolina University, Cullowhee, NC 28723 USA.

Bing Li, Dept. of Electrical and Computer Engineering, University of Virginia, Charlottesville, VA 22904 USA.

Chris D. Garson, Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA 22908 USA.

Scott T. Acton, Dept. of Electrical and Computer Engineering and also the Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA 22904 USA.

John A. Hossack, Dept. of Biomedical Engineering, University of Virginia, Charlottesville, VA 22908 USA.

1. McCann HA, Sharp JC, Kinter TM, McEwan CN, Barillot C, Greenleaf JF. Multidimensional ultrasonic imaging for cardiology. Ultrasound in Med. & Biol. 1988 Sept.vol. 76(no. 9):1063–1073.

2. Fenster A, Downey DB, Cardinal HN. Three-dimensional ultrasound imaging. Physics in Medicine and Biology. 2001;vol. 46:R67–R99. [PubMed]

3. Frangi AF, Niessen WJ, Viergever MA. Three-dimensional modeling for functional analysis of cardiac images: a review. IEEE Trans. Med. Imag. 2001 Jan.vol. 20(no. 1):2–25. [PubMed]

4. Jasjit SS. Computer vision, pattern recognition and image processing in left ventricle segmentation: the last 50 years. Pattern Analysis and Applications. 2000 Sept.vol. 3(no. 3):209–242.

5. McInerney T, Terzopoulos D. Deformable models in medical image analysis: a survey. Medical Image Analysis. 1996;vol. 1(no. 2):91–108. [PubMed]

6. Bardinet E, Cohen LD, Ayache N. A parametric deformable model to fit unstructured 3D data. Computer Vision and Image Understanding. 1998;vol. 71(no. 1):39–54.

7. Bardinet E, Cohen LD, Ayache N. Tracking and motion analysis of the left ventricle with deformable superquadrics. Medical Image Analysis. 1996;vol. 1(no. 2):129–149. [PubMed]

8. Montagnat J, Delingette H. Globally constrained deformable models for 3D object reconstruction. Signal Processing. 1998;vol. 71(no. 2):173–186.

9. Chen CW, Huang TS, Arrott M. Modeling, analysis, and visualization of left ventricle shape and motion by hierarchical decomposition. IEEE Trans. Pattern Anal. Machine Intell. 1994 April;vol. 16(no. 4):342–356.

10. Goshtasby A, Turner DA. Segmentation of cardiac MR images for extraction of right and left ventricular chambers. IEEE Trans. Med. Imag. 1995 Jan.vol. 14(no. 1):56–64. [PubMed]

11. Kaus MR, von Berg J, Niessen W, Pekar V. Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer; 2003. Automated segmentation of the left ventricle in cardiac MRI; pp. 432–439.

12. Paragios N. A level set approach for shape-driven segmentation and tracking of the left ventricle. IEEE Trans. Med. Imag. 2003 June;vol. 22(no. 6):773–776. [PubMed]

13. Corsi C, Saracino G, Sarti A, Lamberti C. Left ventricular volume estimation for real-time three-dimensional echocardiography. IEEE Trans. Med. Imag. 2002 Sept.vol. 14(no. 21):1202–1208. [PubMed]

14. Coppini G, Poli R, Valli G. Recovery of the 3-D shape of the left ventricle from echocardiographic images. IEEE Trans. Med. Imag. 1995 June;vol. 14(no. 2):301–317. [PubMed]

15. Gustavsson T, Pascher R, Caidahl K. Model based dynamic 3D reconstruction and display of the left ventricle from 2D cross-sectional echocargiograms. Comput. Med. Imag. and Graphics. 1993;vol. 17(no. 4/5):273–278. [PubMed]

16. McPherson DD, Skorton DJ, Kodiyalam S, Petree L, Noel MP, Kieso R, Kerber RE, Collins SM, Chandran KB. Finite element analysis of myocardial diastolic function using three-dimensional echocardiographic reconstruction: application of a new method for study of acute ischemia in dogs. Circulation Research. 1987 May;vol. 60(no. 5):674–682. [PubMed]

17. Strickels KR, Wann LS. An analysis of three-dimensional reconstructive echocardiography. Ultrasound in Med. & Biol. 1984;vol. 10(no. 5):575–580. [PubMed]

18. Geiser EA, Lupkiewicz SM, Christie LG, Ariet M, Conetta DA, Conti CR. A framework for three-dimensional time-varying reconstruction of the human left ventricle: sources of error and estimation of their magnitude. Computers and Biomedical Research. 1980;vol. 13:225–241. [PubMed]

19. Chu CH, Delp EJ, Buda AJ. Detecting left ventricular endocardial and epicardial boundaries by digital two-dimensional echocardiography. IEEE Trans. Med. Imag. 1998 July;vol. 7(no. 2):81–90. [PubMed]

20. Kass M, Witkin A, Terzopoulos D. Snakes-active contour models. International Journal of computer Vision. 1988 Jan.vol. 1(no. 4):321–331.

21. Cohen LD, Cohen I. Finite element methods for active contour models and balloons for 2-D and 3-D images. IEEE Trans. Pattern Anal. Machine Intell. 1993 Nov.vol. 15(no. 11):1131–1147.

22. Xu C, Prince JL. Snakes, shapes, and gradient vector flow. IEEE Trans. Image Processing. 1998 March;vol. 7(no. 3):359–369. [PubMed]

23. Xu C, Prince JL. Gradient vector flow deformable models. In: Bankman IN, editor. Handbook of Medical Imaging. Orlando, FL, U.S.A.: Academic Press; 2000.

24. Li B, Acton ST. Active contour external force using vector field convolution for image segmentation. IEEE Trans. Image Processing. 2007 Nov.vol. 16(no. 8):2096–2106. [PubMed]

25. Li B, Acton ST. Automatic active model initializations via Poisson inverse gradient. IEEE Trans. Image Processing. 2008 in press. [PubMed]

26. Nagao M, Matsuyama T. Edge preserving smoothing. Computer Graphics and Image Processing. 1979 April;vol. 9(no. 4):394–407.

27. Lee JS. Digital image enhancement and noise filtering by use of local statistics. IEEE Trans. Pattern Anal. Machine Intell. 1980 March;vol. PAMI-2(no. 2):165–168. [PubMed]

28. Frost VS, Stiles JA, Shanmugan KS, Holtzman JC. A model for radar images and its application to adaptive digital filtering of multiplicative noise. IEEE Trans. Pattern Anal. Machine Intell. 1982 March;vol. PAMI-4(no. 2):157–166. [PubMed]

29. Kuan DT, Sawchuk AA, Strand TC, Chavel P. Adaptive noise smoothing filter for images with signal-dependent noise. IEEE Trans. Pattern Anal. Machine Intell. 1985 March;vol. PAMI-7(no. 2):165–177. [PubMed]

30. Kuan DT, Sawchuk AA, Strand TC, Chavel P. Adaptive restoration of images with speckle. IEEE Trans. Acoust., Speech, Signal Processing. 1987 March;vol. ASSP-35(no. 3):373–383.

31. Loupas T, McDicken WN, Allan PL. An adaptive weighted median filter for speckle suppression in medical ultrasonic images. IEEE Trans. Circuits Syst. 1989 January;vol. 36(no. 1):129–135.

32. Karaman M, Kutay MA, Bozdagi G. An adaptive speckle suppression filter for medical ultrasonic imaging. IEEE Trans. Med. Imag. 1995 June;vol. 14(no. 2):283–292. [PubMed]

33. Yu Y, Acton ST. Speckle reducing anisotropic diffusion. IEEE Trans. Image Processing. 2002 November;vol. 11(no. 11):1260–1270. [PubMed]

34. Hao SGX, Gao X. A novel multiscale nonlinear thresholding method for ultrasonic speckle suppressing. IEEE Trans. Med. Imag. 1999 Sept.vol. 18(no. 9):797–794. [PubMed]

35. Achim A, Bezerianos A, Tsakalides P. Novel Bayesian multiscale method for speckle removal in medical ultrasound images. IEEE Trans. Med. Imag. 2001 Aug.vol. 20(no. 8):772–783. [PubMed]

36. Lopes A, Touzi R, Nezry E. Adaptive speckle filters and scene heterogeneity. IEEE Trans. Geosci. Remote Sensing. 1990 Nov.vol. 28(no. 6):992–1000.

37. Michailovich OV, Tannenbaum A. Despeckling of medical ultrasound images. IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 2006 Jan.vol. 53(no. 1):64–78. [PMC free article] [PubMed]

38. Dutt V, Greenleaf J. Adaptive speckle reduction filter for log-compressed B-scan images. IEEE Trans. Med. Imag. 1996 Dec.vol. 15(no. 6):802–813. [PubMed]

39. Tay PC, Acton ST, Hossack JA. Proc. IEEE Int’l. Symp. on Biom. Imag. Arlington, Virginia, U.S.A.: 2006. Apr 6–9, A stochastic approach to ultrasound despeckling; pp. 907–910.

40. Sethian J. Level Set Methods and Fast Marching Methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. New York, NY: Cambridge University Press; 1999.

41. Klingensmith JD, Shekhar R, Vince DG. Evaluation of three-dimensional segmentation algorithms for the identification of luminal and medial-adventitial borders in intravascular ultrasound images. IEEE Trans. Image Processing. 1996 June;vol. 5(no. 6):996–1011. [PubMed]

42. Taubin G. A signal processing approach to fair surface design. SIGGRAPH. 1995:351–358.

43. Spelke ES. Principles of object perception. Cognitive Science. 1990;vol. 14:29–56.

44. Agnew LRC, Aviado DM, Brody JI, Burrows W, Butler RF, Combs CM, Gambill CM, Glasser O, Hine MK, Shelley WB, Daly LW, editors. Dorland’s Illustrated Medical Dictionary. 24th ed. Philadelphia: W. B. Saunders company; 1965.

45. Chérin E, Williams R, Needles A, Liu G, White C, Brown AS, Zhou Y, Foster FS. Ultrahigh frame rate retrospective ultrasound microimaging and blood flow visualization in mice *in vivo*. Ultrasound in Medicine & Biology. 2006 May;vol. 32(no. 5):683–691. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |