IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

2393

Renal Segmentation From 3D Ultrasound via Fuzzy Appearance Models and Patient-Specific Alpha Shapes Juan J. Cerrolaza,∗ Nabile Safdar, Elijah Biggs, James Jago, Craig A. Peters, and Marius George Linguraru Abstract— Ultrasound (US) imaging is the primary imaging modality for pediatric hydronephrosis, which manifests as the dilation of the renal collecting system (CS). In this paper, we present a new framework for the segmentation of renal structures, kidney and CS, from 3DUS scans. First, the kidney is segmented using an active shape model-based approach, tailored to deal with the challenges raised by US images. A weighted statistical shape model allows to compensate the image variation with the propagation direction of the US wavefront. The model is completed with a new fuzzy appearance model and a multi-scale omnidirectional Gabor-based appearance descriptor. Next, the CS is segmented using an active contour formulation, which combines contour- and intensity-based terms. The new positive alpha detector presented here allows to control the propagation process by means of a patient-specific stopping function created from the bands of adipose tissue within the kidney. The performance of the new segmentation approach was evaluated on a dataset of 39 cases, showing an average Dice’s coefficient of 0.86 ± 0.05 for the kidney, and 0.74 ± 0.10 for the CS segmentation, respectively. These promising results demonstrate the potential utility of this framework for the US-based assessment of the severity of pediatric hydronephrosis. Index Terms— Active shape model, alpha shape, collecting system, Gabor filters, hydronephrosis, image segmentation, kidney, ultrasonography.

I. I NTRODUCTION

U

LTRASOUND (US) imaging is the preferred medical imaging modality in pediatric care, due to its non-ionizing and non-invasive properties, its real-time nature,

Manuscript received April 11, 2016; revised May 18, 2016; accepted May 19, 2016. Date of publication May 24, 2016; date of current version October 27, 2016. This project was supported in part by a philanthropic gift from the Government of Abu Dhabi to Children’s National Health System. Asterisk indicates corresponding author. ∗ J. J. Cerrolaza is with the Sheikh Zayed Institute for Pediatric Surgical Innovation, Children’s National Health System, Washington DC 20010 USA (e-mail: [email protected]). N. Safdar is with the Department of Radiology and Imaging Sciences, Emory University School of Medicine, Atlanta, GA 30322 USA (e-mail: [email protected]). E. Biggs is with the Sheikh Zayed Institute for Pediatric Surgical Innovation, Children’s National Health System, Washington DC 20010 USA (e-mail: [email protected]). J. Jago is with the Philips Healthcare, Bothell, WA 98021 USA (e-mail: [email protected]). C. A. Peters is with the Department of Pediatric Urology, Children’s Medical Center of Dallas, Dallas, TX 75235 USA (e-mail: [email protected]). M. G. Linguraru is with the Sheikh Zayed Institute for Pediatric Surgical Innovation, Children’s National Health System, Washington DC 20010 USA, and also with the Department of Radiology and Pediatrics School of Medicine and Health Science, George Washington University, Washington DC 20037 USA (e-mail: [email protected] ). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMI.2016.2572641

Fig. 1. Renal segmentation in US. (a) US longitudinal view of a hydronephrotic kidney. (b) Diagram describing the proposed framework.

safety, and low cost. In particular, US is the mainstay of imaging for the evaluation of the kidney and the urinary tract in children. One of the most common finding in abdominal US studies is hydronephrosis (HN), affecting 2–2.5% of children [1]. HN is defined as the dilation of the renal collecting system (CS, which comprises the renal pelvis and calices) and the distortion of the renal parenchyma as a result of obstruction of the outflow of urine (see Fig. 1(a)), and can manifest in a wide spectrum of severity. An early diagnosis is very important in order to identify those cases with severe obstruction of the urinary tract that require immediate surgical intervention, and prevent permanent kidney damage. The Society for Fetal Urology (SFU) proposed a 5-level scale for classifying the severity of HN [2] based on the visual assessment of the degree of dilation of the CS. However, the potential of this grading system as diagnostic tool is limited by the subjective interpretation of radiologists [3]. Shapiro et al. [4] proposed one of the first quantitative measure of HN severity, the HN index (HI), defined as the ratio of the CS area to the total area of the kidney. More recent works [5] have shown the potential clinical value of automated US imaging processing for the quantitative, reproducible, and standardized assessment of pediatric HN, using manually segmented contours to correlate renal morphology to renal function. In this context, an accurate parameterization and segmentation of the kidney anatomy becomes essential to guarantee optimal performance and robustness. Despite considerable improvements in the image quality in recent years, US imaging remains, arguably, the most difficult medical imaging modality upon which to perform segmentation. Low signalto-noise ratio, signal attenuation and dropout, and missing boundaries due to variability with the orientation of the probe during image acquisition, make the analysis of organs and objects of interest particularly challenging in US. There has been increasing interest in developing new segmentation methods for sonographic images [6]–[11] in

0278-0062 © 2016 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

2394

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

application areas like echocardiography or transrectal US. However, there is a dearth of literature on the segmentation of renal structures, especially in 3DUS. Xie et al. [7] presented a 2D segmentation method based on texture and shape priors, using an expectation-maximization Gaussian-mixture approach and a shape model. Mendoza et al. [8] presented one of the first frameworks to segment both the kidney and the CS from 2DUS images also using appearance and shape models. However, the segmentation, and thus the computation of HI can be significantly affected by the selection of a single 2D sagittal slice. The 3D extension of the HI, 3DHI, recently presented in our earlier work [9] using manual segmentations of the kidney, showed the potential of this new imaging biomarker for the objective and accurate assessment of severity of HN. However, in spite of the superiority of 3DHI over the 2DHI, its computation demands the development of new tools that alleviate the tedious and subjective manual segmentation process. Martín-Fernández and Alberola-López [10] proposed one of the few kidney segmentation methods in 3DUS. However, a slice-by-slice approach is used to estimate the contour of the kidney. Recently, Ardon et al. [11] presented an interactive tool for the detection and segmentation of the kidney in 3DUS. In spite of the promising results reported, further validation with pediatric hydronephrotic cases is required. To the best of our knowledge, the 3D segmentation of CS has not yet been addressed. Despite the fact that CS can be partially identified in US in hydronephrotic kidneys thanks to its hypoechoic nature, its segmentation is particularly challenging. The lack of a definite shape, its blurred boundaries and the presence of heterogeneous structures of different shapes and intensities inside the kidney (i.e., renal pyramids visible in young children) can easily confuse even the experienced eye (see Fig. 1(a)). This paper presents a complete framework for the semiautomatic segmentation and quantification of renal structures (kidney and CS) in 3DUS. The framework was divided in two parts to deal effectively with the particular challenges that arise when working with US renal images. First, Section II introduces a new kidney segmentation algorithm via a new Gabor-based fuzzy appearance model (FAM). The segmentation process incorporates shape priors and appearance model tailored to deal with the high intensity variability and inhomogeneity of US images. Once the kidney is segmented, Section III details a new active contour-based formulation that mimics the evolution of the HN within the kidney. In particular, we combine a new positive delta detector [21] to identify the renal fat in the kidney. This important anatomical constraint allows us to define a new patient-specific stopping function using alpha shapes [26] for the accurate segmentation of the renal CS. Finally, complete segmentation framework is tested in Section IV on a set of 39 cases, both healthy and pathological. The main elements of this framework are depicted in Fig. 1(b). II. K IDNEY S EGMENTATION VIA F UZZY A PPEARANCE M ODELS Based on the well-known active shape models (ASM) introduced by Cootes et al. [12], we proposed a new tailored

Fig. 2. Effect of wavefront orientation on the contrast in US images. (a) Slice of a 3DUS image of the kidney. The estimated propagation directions of the wavefront ( di ) at landmark l i is depicted as green vectors; the normals to the kidney surface (ni ) are shown in cyan; the center of the probe (C) is shown in red. The white dots indicate areas where the propagation direction is almost tangent to the contour of the kidney (i.e., normal to ni ), causing fading effects of the edges. (b) Weights assigned to the landmarks according to the orientation of the kidney surface relative to the US probe (γ = 0.65).

approach to the imaging physics of US imaging in [13]. Inspired by the promising results obtained in this paper, we incorporate a weighted statistical shape model to address the contrast dependency with the propagation direction of the US wavefront. Additionally, we introduce a new weighted fuzzy appearance model (FAM) to deal efficiently with the intensity variability of medical images in general, and US scans in particular. A. Wavefront-Corrected Statistical Shape Model Each kidney, x, is defined as a set of K ∈ N+ 3D landmarks distributed across the surface. Using principal component analysis over the aligned training set, {x a }, it is possible to define a subspace of allowed-shapes by means of the linear equation y = x + P · b, where x is the mean shape, and P is the (3K × t) matrix formed by the t ∈ N+ principal eigenvectors required to explain the 98% of the total variance in the training set. Generally, all the landmarks are considered equally relevant when calculating the corresponding shape parameters b of a new instance y. In particular, b = P T ( y−x) can be considered as the principal component projection of y that minimizes the squared error function E R R =  y − x a 2 . Whereas these expressions provide satisfactory results in most medical applications, in US imaging it is known that those edges tangent to the propagation direction of the wavefront can be affected by fading effects [32] (see Fig. 2(a)). These effects hinder the correct localization of the corresponding landmarks, and hence the segmentation process. Suppose W is a (3K × 3K ) diagonal matrix defining the weight (or reliability) of each landmark. Thus, the error function can be redefined as E R RW = ( y − x a )T W( y − x a ), and the corresponding shape parameters, bW , that minimizes it can be → − obtained as bW [ P T W P]−1 P T W( y−x). Suppose now that d l is an unitary direction vector that represents the propagation direction of the wavefront at the i -th landmark, l i . Defining → − the center of the US probe as C, d l can be approximated as → − −−→ −−→ −−→ d l ≈ C, l l /C, l l , where C, l l  represents the Euclidean −−→ norm of vector C, l l . Therefore, the weight wi associated with → − l i in the shape model can be defined as a function of d l , and → − nl , the unitary vector normal to the kidney surface at l i ,   2 → → γ − a cos |0,π/2 ( d l · − wi = 1 − nl ) , (1) π

CERROLAZA et al.: RENAL SEGMENTATION FROM 3D ULTRASOUND VIA FUZZY APPEARANCE MODELS AND PATIENT-SPECIFIC ALPHA SHAPES

2395

where a cos |0,π/2 (·) represents the inverse cosine function in the range [0, π/2], and γ ∈ R is a configuration parameter of the power law function. From (1) it can be observed how → − → nl are the weight is 0 for those landmarks where d l and − orthonormal, and 1 if parallel, as Fig. 2(b) illustrates. B. Fuzzy Appearance Model (FAM) While the shape prior model described in Section II.A ensures the legitimacy of the shapes obtained during the segmentation process, an adequate texture model is another cornerstone of most US image segmentation methods [6]. Traditionally, ASM appearance models are based on the normalized first derivative of the gray profiles normal to the boundary of the object and centered at each landmark [12]. Despite the popularity of this simple model, more sophisticated alternatives have emerged over time trying to overcome some of the limitations of the original model [14], [15]. However, all these approaches use a single statistical model to characterize the appearance and texture around each landmark, assuming that there exists a structural consistency among the dataset. While this is just an approximation, differences in the surrounding tissue between patients, depth-dependent attenuation in US images, or inaccuracies in the correspondence between landmarks can lead to noisy or even uninformative appearance models. Fig. 3(b) illustrates the intensity profiles of the six different subjects shown in Fig. 3(a) at the same anatomical location between the right kidney and the liver. Typically, the liver is slightly brighter than the kidney (case 3). However, the existence of a highly echogenic bands of adipose tissue between both organs (cases 5 and 6), or the location of the probe (case 4) can generate different patterns of appearance. The resulting average profile (μ) is an uninformative intensity pattern as Fig. 3(c) illustrates. As alternative to the classic single-model approaches, we propose the creation of multiple appearance models for each landmark in order to capture the inherent differences between datasets, whether they are due to anatomical variability or imaging parameters. We use fuzzy clustering theory to identify the different appearance patterns of the anatomical region around the landmark l i . Suppose {a is }s=1...S represents the training set to model the appearance around l i , where S ∈ N+ is the number of training samples. In the most general scenario, ai,s ∈ Rn represents the S th n-dimensional training sample (e.g., a n-dimensional vector containing the intensity profile normal to the contour and centered at l i ). Thus, given {ais }s=1...S , we use a tailored version of the fuzzy clustering algorithm proposed by Gustafson and Kessel [16] (see Appendix) to identify Ti different appearance patterns for each landmark, l i . Each one of these patterns is defined by {μi j , M i j } j =1 , . . . , Ti , where μi j ∈ Rn and M i j ∈ Rn×n represent the mean (center) and the norm-inducing matrix of the j -th cluster respectively. Unlike other fuzzy clustering approaches, {μi j , M i j } defines the inner-product norm for each cluster, (ai,s − μi j )T M i j (a i,s − μi j ), which allows to generate clusters of different geometrical shapes. Note that in the particular case where Ti = 1, M i and μ i represent the sample covariance matrix and the sample mean, respectively, as commonly used in ASM [12].

Fig. 3. Appearance inconsistency between datasets. (a) Slices of 3DUS scans of six different cases, showing a longitudinal section of the right kidney. The red point identifies the same anatomical location for all the images; the green lines indicate where the intensity samples were extracted. The liver (yellow ‘L’) and the right kidney (cyan ‘K’) are also identified in the pictures for clarity. (b) Intensity profiles extracted from each one of the six images depicted in (a). (c) Average intensity patterns obtained when averaging all the samples simultaneously (μ:Avg. Profile), and when using the weighted fuzzy clustering approach described in Section II-B (μ1 : Avg.Cluster 1 and μ2 : Avg.Cluster 2).

As discussed in Section II-A, the contrast, and thus the appearance, around landmark l i can be significantly affected by the propagation direction of the wavefront at that location. Thus, for each sample, ai,s , it is possible to define a weight Wis (see (1)) that controls the impact of that sample when building the appearance models. The Appendix extends the original fuzzy clustering framework proposed by Gustafson and Kessel [16] for the more general scenario of a weighted set of samples {ai,s , Wi,s }. The optimal number of clusters, Ti , is automatically defined by means of a weighted version of the validation index proposed by Tang et al. [17] (see (A18) in Appendix). In order to increase the robustness of the appearance model against possible inaccuracies in the correspondence between landmarks (i.e., small differences in the anatomical location defined by l i through the training set), we expand the training set {ais } by including the appearance information from adjacent landmarks. In particular, the appearance training set of landmark l i , {ais }, consists of texture samples extracted from l i and its neighboring vertices within a ring of length C.1 When segmenting a new image using an ASM-based algorithm, the location of each landmark is updated to that position that maximizes the probability of coming from the learned distribution, i.e., minimizing the Mahalanobis distance to the mean profile learned from the training set (see [12] for details). Similarly, having now Ti different models, {μi j , M i j } j =1,...,Ti , the optimal location for each landmark will be defined by 1 A ring center at l of length C ∈ N+ , R C , is formed by all the vertices l whose face-based distances to l is ≤ C. Thus, Rl1 is defined by all the vertices

from those faces of which l is member (i.e., face-based distances = 1). Similarly, Rl2 is defined by all the vertices whose face-based distance to Rl1 is equal to 1; and so on.

2396

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

Fig. 4. Gabor-based appearance model. (a) Omnidirectional Gabor filter estimation. Vector nl (blue) represents a random direction in the 3D space. G f,θi ,ϕi is approximated by the four nearest discrete directions considered when computing the Gabor filter (red and green). (b) Example of the original US image of the kidney. The normal direction of the kidney contour at the position of landmark l i is defined by the vector nl . (c) Imaginary component of G f,θi ,ϕi at frequency f = 0.0625. (d) Imaginary component of G f,θi ,ϕi at a higher frequency, f = 0.1768.

that model that minimizes the Mahalanobis distance, using the corresponding fuzzy mean and covariance matrix. C. Omnidirectional Gabor-Based Appearance Descriptor Traditional intensity-based appearance models are particularly inefficient when dealing with US images, due also to the aforementioned inherent challenges, such as speckle and low contrast between areas of interest, among others. Here we propose a Gabor filter-based appearance model, as alternative to the classic intensity-based approaches. Gabor filters have been used in US image processing for edge detection, texture representation and discrimination, mainly in 2D [7], [18]. Notably, Zhan and Shen [19] used a Gabor filter bank to extract and characterize texture features in 3DUS image of the prostate. However, only 2D Gabor filters located at two orthogonal planes were used. In this paper we propose an omnidirectional Gabor-based appearance model for 3DUS images. A Gabor filter can be expressed mathematically as: 1 2 2 2 2 e−(x +y +z /2σ ) g f,θ,ϕ (x, y, z) = 3/2 3 (2π) σ · e− j 2π(ux+v y+wz), (2) where u = f sin θ cos ϕ; v = f sin θ sin ϕ; and w = f cos θ ; f is the central frequency of the sinusoidal plane wave, and ϕ and θ are the orientation parameters that together with the Gaussian scale parameter, σ , determine the Gabor filter in 3D. The number of filters, and thus the computational cost, increases significantly with the number of orientations, especially in 3D. Typically, only a discrete number of orientations are considered (i.e., θm = mπ/M|m=0,...,M−1 and ϕn = nπ/N|n=0,...,N−1 ), which limits the capacity of the filter to extract texture features in any direction (θ, ϕ) [18]. Suppose G f,θm ,ϕn represents the filtered 3DUS volume, I , using one of the M · N Gabor filters that sample the entire 3D space. Suppose now the angles θi and ϕi represent respectively the → zenith and azimuth of the vector − nl , the unitary vector normal to the kidney at landmark l i . Thus, G f,θi ,ϕi can be estimated as G f,θi ,ϕi = (1 − ηi )(1 − βi )G f,θmi ,ϕni + (1 − ηi )βi G f,θmi ,ϕni+1 + (3) ηi (1 − βi )i G f,θmi+1 ,ϕni + ηi βi G f,θmi+1 ,ϕni+1 , where m i = θi /(π/M) , n i = ϕi /(π/N) , ηi = (θi /(π/M))−m i , and βi = (ϕi /(π/N))−n i ·G f,θi ,ϕi are omnidirectional approximations of the discrete filter bank G f,θm ,ϕn

computed for each landmark, l i , during the ASM-based iterative segmentation [12] (see Fig. 4(a)). More specifically, we use the imaginary component to create a new texture model for each landmark. Figs. 4(b)–4(d) show how the omnidirectional Gabor filter is able to identify the contour of the kidney in the vicinity of each landmark. Additionally, the multiscale nature of Gabor filter banks allows to characterize textures with different dominant sizes (i.e., using different central frequencies, f ), and thus, to improve the robustness of the segmentation method to local minima. Starting with the lowest frequency (Fig. 4(c))), the coarse Gabor features are used in the initial stages of the segmentation. As the algorithm evolves, the resulting shape becomes closer to the target, using higher values of f (Fig. 4(d))). This framework enables us to hierarchically focus on different image features at different stages of the algorithm. Since the texture information is different at each resolution, a new fuzzy appearance model (Section II-D) is specifically created for each frequency, f . Thus, the Gabor-based fuzzy appearance model for the i -th landmark at frequency f is defined by {μ f i j , M f i j } j =1,...,T f i (see Section II-B). D. Kidney Segmentation Algorithm Algorithm 1 summarizes the key elements of our kidney segmentation method. Typically, most segmentation methods for US imaging require user interaction to initialize the algorithm, or impose hard positional constraints on the location of the target object [11], [18]. Our algorithm requires minimal user intervention by selecting point clicks to roughly define the principal axis of the kidney. These principal axes can be defined as the main axis of the ellipsoid circumscribing the kidney. The effect of the precision of these initial points over the segmentation accuracy is discussed in Section IV. III. C OLLECTING S YSTEM S EGMENTATION Once the contour of the kidney has been delineated, the segmentation of the CS inside the kidney is addressed by means of the active contour formulation presented in Section III-A. In particular, the evolution equation introduced here incorporates a new patient-specific stopping function (S F(·)) using the renal fat as anatomical constraint. The renal fat surrounding the CS is automatically detected thanks to the new positive delta detector described in Section III-B. This fat of the renal sinus is used in Section III-C to define the alpha shape-based

CERROLAZA et al.: RENAL SEGMENTATION FROM 3D ULTRASOUND VIA FUZZY APPEARANCE MODELS AND PATIENT-SPECIFIC ALPHA SHAPES

Algorithm 1 Kidney Segmentation Algorithm 1. Initialization. 2. Detect the center of the probe, C. 3. for f = f min to f max // Multiscale loop. 4. while (NOT CONVERGENCE) or (MAX ITERATIONS) do 5. for i = 1 to K // Landmarks updating process. → 6. Calculate normal vector, − n . − → −−→ −−l→ 7. Calculate d l = C, l l /C, l l , and wi using (1). 8. Calculate G f,θi ,ϕi using (3). 9. Extract texture samples. 10. Update location using {μ f i j , M f i j } j =1,...,T f i // Using the Mahalanobis distance; see [12] for details. 11. end 12. Define y by concatenating the updated landmarks and mapping it to the normalized shape space. 13. Calculate and constraint bw . 14. Define the new kidney shape. 15. end 16. end patient-specific stopping function that controls the evolution of the segmentation process. A first description of the CS segmentation method was recently presented in [29] and validated on a limited dataset of 13 cases. A revised and detailed description of the method is presented below. A. Active Contour Formulation Based on the original active contour formulation proposed by Caselles and Kimmel [20], we propose an energy functional that combines contour and intensity-based terms, and incorporates a new patient-specific positional map as additional stopping criteria. Suppose I : ← − R+ represents a 3D gray level image in the image domain ⊂ R3 , and U : (t, ) ← −R is an implicit representation of the CS at time t, i.e., the CS coincides with the set of points U(t, .) = 0. Here we define the evolution equation of U as ∂U = S F(I)(β · Cont (Iκ, c, U) + (1 − β) · I nt (I , U)) (4) ∂t where S F(I ) represents the new aforementioned stopping function (see Section III-B), and β ∈ [0, 1] is a constant that balances the contour- and the intensity-based terms, Cont and I nt respectively. In particular, Cont (I, κ, c, U) = g(I)|∇U|(κ + c) + ∇g(I) · ∇U,

(5)

where κ = di v(∇U/|∇U|) is the curvature term computed on the level set of U, c ∈ R+ is a constant velocity term, and g(I) is an inverse edge indicator function of the image I. Typically, g(I) is a gradient-based edge detector, e.g., g(I) = 1/(1 + |∇ I|), where  I is a smooth version of I. However, these very simple edge detectors generally perform poorly in US. Alternatively, we use a local phase-based step function detector, the feature asymmetry (FA) detector [22]. The mathematical formulation of FA, and the resulting edgebased stopping function, are detailed in g(I) Section III-B. Moreover, a formulation based exclusively on the expansive

2397

forces described in (5) would turn inefficient when segmenting objects with weak or missing boundaries. Here, we combine the original gradient-based model with I nt (I , U), a minimal intensity variance term [23], defined as   I nt (I , U) = λout (I − μout )2 + λin (I − μin )2 ) |∇U|; (6) μout , and μin are the mean intensities in the exterior and the interior of the CS, respectively, and λout and λin are two control parameters generally defined as λout = λin = 1. Intuitively, this new intensity-based term looks for the best separating contour in I, and the optimal expected values μout and μin . Given the hypoechoic nature of the CS in US images (i.e., μin ≈ 0), the second term of (6) prevents the evolution of the contour into brighter areas, whereas the first term acts as expansive force toward dark areas (i.e., toward the CS). B. Local Phase-Based 3D Positive Delta Detection With the evolution of HN, part of the renal fat originally located in the renal pelvis is displaced into the kidney, surrounding the dilated CS. These bands of fat constitute a key visual clue for the radiologist to differentiate the CS from other hypoechoic structures like the renal pyramids (see Fig. 1(a)). The aim of the new positive delta detector (PDD) presented here is to identify these echogenic regions and to incorporate this anatomical information into the new patient-specific stopping function S F(I ), defiend in Section III-C. To define PDD, we use local phase analysis of the monogenic signal [24], a n-dimensional generalization of the Hilbert transform-based analytic representation of 1D signals, whose potential in echocardiography images was recently reported [22]. Using the Riesz transform, the monogenic signal of a 3DUS image, I, is defined as the 4D vector, I M = (I B P , I R ); I B P is the resulting image of band-pass filtering I, and I R = (I Rx , I Ry , I Rz ) =(I B P ∗h x , I B P ∗h y , I B P ∗h z ) represents the three Riesz filtered components. The spatial representation of the Riesz filters is defined as hk = −k/(2π(x 2 + y 2 + z 2 )3/2 ), where k represents one of the three coordinates of the 3DUS image, i.e., k = x, y or z. Using an isotropic log-Gabor filter with central frequency ω0 , g L G , ω0 the monogenic signal can be represented in polar form as even ω0 = g L G,ω0 ∗ I,  (g L G,ω0 ∗ I Rk )2 . oddω0 = k=x,y,z

(7) (8)

= The local phase can thus be defined as ω0 atan(even ω0 /oddω0 ). The local phase contains structural information of the image, such as the location and orientation of image features, transitions and discontinuities. Typically, these properties are used to detect feature asymmetries in images, e.g., step edges, identifying those points whose absolute value of the local phase is 0 (positive edges), or π (negative edges) [21], [30]:  |oddω | − |even ω | − Tω

, (9) FA = ω oddω 2 + even ω 2 + ε where the · operator zeros the negative values, ε is a small positive constant to prevent division by zero, and Tε is a scale

2398

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

Fig. 5. Patient specific positional map. (a) Longitudinal view of a hydronephrotic kidney extracted from the 3DUS scan. (b) Example of the positive delta detection for the image shown in (a). (c) 3D alpha shape (α = 40; the wireframe is shown in green) and resulting positional-mapbased stopping function from the fat points inside the kidney. (d) 3D alpha shape (light yellow) surrounding the CS (blue). The detected fat points are shown in dark yellow and the kidney in light red.

specific noise threshold defined as Tω = exp(mean[log((oddω 2 + even ω 2 )1/2 )]).

(10)

As alternative to the traditional intensity-based approaches, we define the new edge detector in (5) as g(I) = 1 − F A, whose satisfactory performance as edge detector in US images was emphasized by previous works [22]. Here, we also exploit local phase properties to detect symmetrical image features, positive deltas identified with points whose local phase is close to +π/2 (i.e., points where |even ω |  |oddω | and sign(even ω · oddω ) > 0. These points indicate relatively thin echogenic regions (i.e., ridges) inside the kidney, corresponding to fat tissue. Since positive deltas are scale-dependent features, the multi-scales PDD is defined as  |even ω |−|oddω |−Tω ·sign(even ω ·oddω)

. PDD = ω oddω 2 + even ω 2 + ε (11) It can be observed that PDD ∈ [0, 1] takes values close to 1 near bright bands (i.e., positive delta features), and close to zero otherwise (see Fig. 5(b)). C. Patient Specific Positional Maps Using Alpha Shapes Using the fat areas identified inside the kidney, we create an anatomically-justified stopping criteria for the active contour formulation able to prevent the leakage of the contour outside the region delimited by the fat bands. Mathematically, this region can be defined as the interior of the continuous surface circumscribing the fat points located via PDD (i.e, those points where PDD ≈ 1). In our case, we define those points as PDDT h (I) = PDD ≥ 0.8. However, a densely-sampled point cloud is required by most of the existing point setbased surface reconstruction techniques [33]. That is not the case in most of hydronephrotic kidneys, where the scattered distribution of the inner thin fat results in an unstructured set of dispersed points (Fig. 5(b)). Therefore, we propose the use of 3D alpha shapes as alternative [26]. The concept of alpha shapes is a generalization of the convex hull that formalized the intuitive notion of shape for any random spatial point set data, including non-convex and even non-connected sets of points in 3D. Given a set of

points, an alpha shape, Sα , is a family of 3D polyhedrons (i.e., geometrical volume with flat polygonal faces) defined by the configuration parameter α ∈ R+ . An edge of Sα is defined between two members of the finite set of points if there exists a generalized sphere of radius 1/α containing the entire point set and which has the property that the two points lie on its boundary. In particular, S0 represents the convex hull defined by the points. Given Sα defined by the set of fat points located inside the kidney, the new stopping function S F(I) can be defined as 1 S F(I) = , (12) (1 + D(Sα (PDDT h (I))) )τ where D(Sα (PDDT h (I))) is the signed distance to the alpha shape Sα (PDDT h (I)), taking negative or positive values inside and outside the region, respectively; τ ∈ [1, +∞) is a control variable. The value of S F(I ) will be 1 inside the alpha shape, and close to zero as we move away from it, thus gradually penalizing the leaking of the contour outside of Sα(PDD T h ) (see Fig. 5(c)). Note that the computation of PDD (11), and thus S F (12), can be performed offline for the entire image domain . D. Initialization of CS Segmentation The initialization of the active contour formulation described in Section III-A is fully automated by selecting the darkest region within the Sα closest to the ureteropelvic junction (UPJ) of the kidney. UPJ is an anatomical region that can be automatically identified thanks to the landmark correspondence between cases required to create the statistical shape model of the kidney. Thus, the UPJ can be identified with the position of a predefined landmark in the kidney at the junction between the ureter and the renal pelvis of the kidney. The automatic initialization process provides a valid seed in most of the cases (see Section IV). A simple interactive process (1 mouse click) can be used to refine or correct invalid seeds if necessary (∼ 6% of cases in our experiments). IV. R ESULTS AND D ISCUSSION A. Materials In our experiments, we used a proprietary database of 39 pediatric 3DUS images (B-mode images; 36 patients; 32 right and 7 left kidneys; 16 female and 23 male; mean age: 51.5 months) acquired under an IRB approved protocol. The database includes 21 healthy and 18 hydronephrotic cases of variable severity. The SFU grades of the pathological kidneys ranges from 1(mild dilation) to 4 (severe dilation), with an average value of 2.2. The images were acquired from a Philips iU22 system with an X6–1 xMATRIX array transducer. The average volume size is 484 × 404 × 256 voxels, and the resolution ranges from 0.15 mm to 0.82 mm. All the images were initially transformed to a world space coordinate system (in mm). For each image, the reference shapes (kidney and CS) were delineated manually by two experienced members of our team, under the supervision of an expert radiologist. The accuracy of our segmentation framework was evaluated in terms of the symmetric point-to-surface distance (SP2S), Dice’s coefficient (DC), and relative-volumedifference (RVD) [31], using leave-one-out cross-validation.

CERROLAZA et al.: RENAL SEGMENTATION FROM 3D ULTRASOUND VIA FUZZY APPEARANCE MODELS AND PATIENT-SPECIFIC ALPHA SHAPES

2399

Significance was assessed using the Wilcoxon rank sum test and a p-value of 0.05. B. Kidney Segmentation The dense point correspondence between all shapes in the training set was established via mesh-to-mesh non-rigid registration [27], using 600 landmarks to describe each shape, and obtaining an average post-registration Dice’s coefficient of 0.99 (i.e., Dice’s coefficient between the original ground truth and the moving mesh). The statistical shape model of the kidney was built to explain the 98% of the shape variability observed in the training set, restricting the deformation space to twice the standard deviation in each deformation vector. The configuration parameter γ was defined experimentally as γ = 0.65, although no significant difference was observed for values between 0.5 and 0.8. The frequency values for the Gabor filter bank (see (2)) were f = {0.176, 0.125}, approximating σ ∼ = 1/ f as proposed in [28]. As described in Algorithm 1, the frequency of the Gabor filters was progressively adapted as the segmentation algorithm evolves. Thus, starting with the lowest frequency, the value of f was modified when convergence was reached, or after 10 iterations. The number of orientations considered in the filter bank was 26: θm = mπ/4|m=0,...,3 , ϕn = nπ/4|n=0,...,3 . The appearance model of each landmark was built from profiles normal to the surface of 5 mm length, while the search space during the landmark updating process was set to 2 mm on each side. The maximum number of iterations was set to 10, though convergence is typically achieved after a few iterations (depending on the quality of the image). The performance of the new algorithm based on FAM was compared to four different alternatives: ASM-F, a version of ASM including the new fuzzy appearance model presenting in Section II-B (i.e., without using the wavefront-corrected shape model), ASM-W, a weighted version of ASM (i.e., no fuzzy appearance model is used), the classical ASM, and the elliptical approximation of the kidney (EA), a rudimentary approximation of the kidney region usually employed in clinical practice. Table I shows the segmentation error for the four algorithms compared here. As described in Section II-D, the initialization of the algorithms requires minimal user intervention to define the semi-axis of the kidney. To study the effect of the initialization in the segmentation process, Table I also shows the accuracy of our approach with added Gaussian noise of different mean value and standard deviation to the location of the initial points. When no noise was added in the initialization process, it can be observed how the new FAM provided, not only significantly better average results, but also lower deviations than ASM-F, ASM-W, ASM and EA, for the three metrics studied: DC (0.86 ± 0.05 vs. 0.82 ± 0.07, 0.80 ± 0.06, 0.81 ± 0.14 and 0.70 ± 0.16, see Fig. 6), SP2S (1.78 ± 0.53 vs. 2.15 ± 0.60, 2.35±0.43, 2.29±0.77 and 3.07±0.89), and RVD (0.18±0.09 vs. 0.21 ± 0.10, 0.23 ± 0.10, 0.22 ± 0.17 and 0.42 ± 0.22), respectively. The better performance of FAM was also statistically verified for the two types of cases under study, healthy and pathological kidneys. Unlike ASM-W, ASM and EA, where a decrease in accuracy was observed for hydronephrotic

Fig. 6. Boxplots of Dice similarity coefficient for the segmentation methods.

kidneys, no significant difference was observed when using a fuzzy appearance model in FAM and ASM-F (DC: 0.86±.0.04 vs. 0.85 ± 0.07 and 0.83 ± .0.05 vs. 0.82 ± 0.10 for healthy and hydronephrotic cases, respectively). In addition to the inherent inter-subject appearance variability, the dilation of the CS in hydronephrotic cases introduces an additional distortion in the appearance model. This effect is particularly relevant in those cases with a severe dilation, i.e., SFU grade of 3 and 4, where parts of the dilated CS are close to the surface of the kidney. Compared to the unique appearance model of ASM-W and ASM, the new fuzzy appearance model deals more effectively with the appearance diversity of the pathological population. FAM also provided a better performance than the four alternatives when a moderate Gaussian noise (mean: 2.0 mm; std: 1.0 mm) was added to the initial points. With the exception of the healthy population, the performance of FAM is statistically better for most of the metrics considered. Interestingly, no significant difference between the original initialization and the noisy one was observed for the hydronephrotic cases when using FAM or ASM-F. However, the performance of FAM decreased when a higher Gaussian distortion (mean: 4.0 mm; std: 1.5 mm) was used in the initialization process. In this case, the difference between FAM and the alternative approaches was not significant, likely because the initialization of the model was further from the kidney surface than the search space in FAM. The performance of FAM was significantly better than EA in all cases R implementation of the FAM kidney segmentation Matlab algorithm was tested on a 64-bit 2.8 GHz processor, with an average execution time of 120 s. C. CS Segmentation The initialization of the CS segmentation process was automatic, as described in Section III-D. A valid initialization seed in the form of a 3×3×3 region inside the CS was detected in all cases but one particularly noisy case with low image quality, when a manual correction of the seed was required. For the monogenic signal analysis, the same configuration parameters described in [24] and [25] were used during the experiments. The rest of parameters were defined empirically, using an iterated grid search approach. In particular, c = 1, λout = λin = 1, β = 0.3, α = 60, τ = 9, and ε = 0.01. Table II shows the segmentation accuracy of the proposed alpha shapes-based active contour formulation (αC) for the three accuracy metrics described above (DC, SP2S, and RVD)

2400

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

TABLE I K IDNEY S EGMENTATION R ESULTS N O N OISE A DDED

TABLE II CS S EGMENTATION A CCURACY

trast between the CS and the renal parenchyma, when the seeds distorted by noise fell close to a local minimum. Although the definition of multiple seeds could alleviate this problem, the single-seed approach described in Section III-D provided satisfactory results in most cases. The CS segmentation algorithm R implementation, with an average was tested on a Matlab execution time of 225 s. The overall segmentation time for the kidney and CS was 346 s, a significant improvement over the several hours it takes for the manual delineation by an operator. D. Hydronephrosis Index

on the subjects with hydronephrosis. In particular, the performance of αC was tested for three different initializations: the automatic initializations defined in Section III-D, and adding Gaussian distortion of different intensities to the initial seed. Additionally, the performance of αC was compared with our 3D implementation of the graph cut-based method (GC) described by Mendoza et al. [8], originally conceived for 2DUS images. It can be observed how the new αC provided significantly better accuracy than GC for all the metrics (DC: 0.74 ± 0.10 vs. 0.56 ± 0.21, SP2S: 1.45 ± 0.66 vs. 3.44 ± 3.00; RVD: 0.28 ± 0.16 vs. 0.58 ± 0.27, respectively). Interestingly, the average DC of the new αC method (0.74 ±0.10) is similar to the inter-user accuracy reported by [8] (0.74 ± 0.18) as assessed by means of a two-sample t-test (p-value > 0.2). When Gaussian noise was added to the initial location of the seed, the segmentation accuracy decreased (see Table II). In particular, the DC decreased to 0.69 ± 0.12 and 0.64 ± 0.20 when adding Gaussian noise (mean: 2 mm-std: 1 mm, and mean: 4 mm-std: 1.5 mm, respectively). Under-segmentation was observed in some particularly noisy cases with low con-

Currently, the SFU grading system is one of the most widely used systems for grading HN sonographically. However, the subjectivity of this grading method results in significant interobserver variation [3]. The segmentation of the renal structures allows to define more objective and quantitative measures of HN severity, such as the HI [4]. Originally conceived as a 2D area-based biomarker, our 3D segmentation framework allows the computation of a new volumetric HI (3DHI, defined as the ratio of the volumes of the CS and the kidney), obtaining an average estimation error of 2.8 ±1.5 in absolute percentage points compared to the manual standard. Compared to 2DHI, the use of the new 3DHI has potential to provide a more objective quantitative estimation of the HN severity [9], suppressing the subjective selection of a single 2D sagittal slice required to compute 2DHI, and allowing objective longitudinal monitoring of HN patient by US imaging. The correlation coefficient between the HI and the SFU grading system provided by an expert was 0.92 (3DHI) and 0.82 (2DHI). Recently that some studies [5] suggested the use of a more comprehensive set of 2DUS-based morphological features for the assessment of HN, including area, parenchymal thickness, and curvature of the kidney and collecting system. In this context, the semi-automatic segmentation framework presented in

CERROLAZA et al.: RENAL SEGMENTATION FROM 3D ULTRASOUND VIA FUZZY APPEARANCE MODELS AND PATIENT-SPECIFIC ALPHA SHAPES

this manuscript represents an important contribution to the 3D extension of the morphological analysis of pathological kidneys, where a manual delineation is impracticable in the clinical practice.

This paper presents a new comprehensive segmentation tool of renal structures in 3DUS, of particular clinical relevance in the diagnosis and assessment of pediatric hydronephrosis. The framework was tailored to deal with the specific challenges raised by the segmentation of US images. First, the kidney was segmented using a new ASM-based approach. The weights defined for each landmark based on the position relative to the probe allowed to integrate the image dependency on the propagation direction of the US wavefront into the segmentation framework. The model was completed with a new fuzzy appearance model for a better and more efficient characterization of the different intensity patterns. The CS was delineated using a tailored approach of the active contour formulation. In particular, the algorithm incorporates patient-specific anatomical constraints that control the contour evolution process and mimics the evolution of hydronephrosis inside the kidneys. For this purpose, we created an alpha shape-based personalized positional map from the bands of fat surrounding the dilated CS. These echogenic bands of adipose tissue were automatically identified thanks to the new positive delta detector for US images presented here. The promising results, in terms of segmentation accuracy and volumetric HI estimation, demonstrate the ability of the proposed framework to objectively assess kidneys with HN and support better informed clinical decision making. The system is able to provide objective and accurate information of the renal parenchymal status, estimating the HI value as reliably as the clinicians, but without the high time and human costs of manual segmentation. Some of the new image processing techniques proposed here were developed in the particular context of the segmentation of renal structures. However, the generality of the ideas proposed here might be applicable to other US segmentation problems, or for the analysis of images with poorly defined edges and landmarks.

j =1

i=1

is a smoothing parameter where m i j = m j (x i ), and α ∈ that controls the fuzziness of the cluster. Note that even though J is also a function of w, the weights are constants for a given set of samples. Constraint (A1) can be integrated by setting m i j = si2j , with si j ∈ R. Constraint (A2) can be adjoined to J by means of Lagrange multipliers as follows:  N T J (s, θ, w, λ) = s 2α di j wi i=1 j =1 i j N T + λi ( si2j − 1). j =1

i=1

By setting ∂ J /∂s = 0, we have the necessary conditions: di j wi + λ∗i ) = 0, ∀i, j si∗j (αsi∗2(α−2) j T si∗j = 1, ∀i,

The mathematical framework presented here is based on the original work on fuzzy clustering presented by Gustafson and Kessel [17]. Given a set of N ∈ N+ n-dimensional samples, {x i }i=1,...,N , a fuzzy partition can be defined as the T-tuple of membership functions, m(·) = {m 1 , . . . , m T } on the feature space  ⊂ Rn such that (A1) (A2)

Suppose 0 ≤ wi ≤ 1 represents the weight associated to the sample x i , and that di j = d j (x i ) = d(x i , θ j ) > 0 represents the distance from the sample x i to the j -th class

(A5) (A6)

j =1

where (∗ ) represents the optimal value. Assuming si∗j = 0, ∀i, j 1/α−1  −λ∗i ∗ mi j = , (A7) αdi j wi and from (A2) it follows that 1 m ∗i j = (α di j wi )−1/α−1

T  1 1/α−1 αdiz wi

z=1

=



T z=1

1

 di j 1/α−1 diz

.

(A8)

From (A3), the corresponding extremum for any θ is J ∗ (θ ) = minm J (m, θ, w)  N T wi = i=1

1/1−α d j =1 i j

1−α .

(A9)

To find the optimal parameter set θ ∗ = (θ1∗ , . . . , θT∗ ) we have ∂di j ∂ J (m, θ, w, λ)  N = m αi j wi = 0. i=1 ∂s ∂θ j

(A10)

Defining di j = (x i − v j )T M j (x i − v j )

A PPENDIX W EIGHTED F UZZY C LUSTERING

j =1

parameterized by θ j . The weighted fuzzy clustering (i.e., the optimal membership function) is defined by the cost function  N T J (m, θ, w) = m αi j di j wi , (A3) R+

V. C ONCLUSIONS

0 ≤ m j (x i ) ≤ 1, 1 ≤ i ≤ N, 1 ≤ j ≤ T T m j (x i ) = 1, 1 ≤ i ≤ N.

2401

(A11)

for M j symmetric and positive-definite, J is linear in M j , giving a singular problem. The cost function can be made as small as desired by simply making M j less positive definite. To get a feasible solution, the determinant |M j | can be constrained as |M j | = ρ j , with ρ j > 0. The resulting cost function J (m, θ, w, λ, β) is  N T m αi j wi di j (v j , M j )+ Jˆ(m, θ, w, λ, β) = i=1 j T T N λi ( m i j − 1) + β j (|M j | − ρ j ). i=1

j =1

j =1

The first necessary condition to minimize Jˆ N ∂ Jˆ = −2 m α wi M j (x i − v j ) = 0, i=1 i j ∂v j

(A13)

2402

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 35, NO. 11, NOVEMBER 2016

leads to the definition of the fuzzy weighted mean μ j :

N m ∗α ( i=1 i j wi x i ) ∗ μ f j = v j = N . (A14) ( i=1 m ∗α i j wi ) Using (A11) in the second necessary condition2 N ∂ Jˆ =0= m α wi (x i −v j )(x i −v j )T +β j |M ∗j |M ∗−1 j , i=1 i j ∂Mj (A15) and the previous constraint |M j | = ρ j , we have M ∗−1 = (ρ j | P f j |)−1/n P f j j where P f j is the fuzzy weighted covariance matrix

N α T i=1 m i j wi (x i − μ f j )(x i − μ f j ) . P fj =

N α i=1 m i j wi

(A16)

(A17)

Now, given a set of data {x}i=1,...,N , and the initial guess (0) (0) (0) θ j = {μ f j , P f j }, the algorithm proceeds as follows: for k = 1, 2, . . . (until convergence) (i) compute di(k) j using (A11). (k)

(ii) compute m i j using (A8). (k+1)

(iii) compute new θ j

using (A14), (A16), and (A16)

end Typically, the sample class mean and the identity matrix are P (0) used to initialize μ(0) f j and f j , respectively. For a given par tition μ f j , P f j , {m i j } , it is possible to defined a weighted version of the validation index defined by Tang et al. [18] as

T N 1 T T 2 2 j =1 i=1 m i j wi di j + T (T−1) j =1 k=1 v j −v k  VT =

min j =k v j − v k

2

+

k = j 1 T

.

(A18) The first component of the numerator indicates the compactness of the fuzzy partition, and the second one is a penalizing term that compensates the decreasing tendency as T ← − N. The denominator measures the interclass difference. Thus, the optimal number of clusters, Topt can be found by minimizing (A18) (e.g., using a simple iterative process). R EFERENCES [1] I. Singh et al., “Pathophysiology of urinary tract obstruction,” in Campbell-Walsh Textbook of Urology, A. Wein, L. Kavoussi, A. Novick, A. Partin, and C. Peters Eds., 10th ed. Philadelphia, PA: Elsevier, 2012. [2] S. K. Fernbach et al., “Ultrasound grading of hydronephrosis: Introduction to the system used by the society for fetal urology,” Pediatric Radiol., vol. 23, no. 6, pp. 478–480, 1993. [3] M. A. Keavs et al., “Reliability assessment of society for fetal urology ultrasound grading system for hydronephrosis,” J. Urol., vol. 180, no. 4, pp. 1680–1683, 2008. [4] S. R. Shapiro et al., “Hydronephrosis index: A new method to track patients with hydronephrosis quantitatively,” Urology, vol. 72, no. 3, pp. 536–539, 2008. [5] J. J. Cerrolaza et al., “Quantitative ultrasound imaging of obstruction severity in pediatric hydronephrosis,” J. Urol., vol. 195, no. 4, pp. 1093–1099, 2016.

2 Using the following identities: ∂(x T A x)/∂ A = x x T , ∂| A|/∂ A = | A| A−1 .

[6] J. A. Noble and D. Boukerroui, “Ultrasound image segmentation: A survey,” IEEE Trans. Med. Imag., vol. 25, no. 8, pp. 987–1010, Aug. 2006. [7] J. Xie et al., “Segmentation of kidney from ultrasound images based on texture and shape priors,” IEEE Trans. Med. Imag., vol. 24, no. 1, pp. 45–57, Jan. 2005. [8] C. S. Mendoza et al., “Automatic analysis of pediatric renal ultrasound using shape, anatomical and image acquisition priors,” in Proc. MICCAI, 2013, vol. 8151. pp. 259–266. [9] J. J. Cerrolaza et al., “The value of 3D hydronephrosis index in the assessment of pediatric hydronephrosis,” in RSNA, 2014. [10] M. Martín-Fernández and C. Alberola-López, “An approach for contour detection of human kidney from ultrasound images using Markov random fields and active contour,” MedIA, vol. 9, no. 1, pp. 1–23, 2005. [11] R. Ardon et al., “Fast kidney detection and segmentation with learned kernel convolution and model deformation in 3D ultrasound images,” in Proc. IEEE Int. Symp. Biomed. Imag., 2015, pp. 268–271. [12] T. F. Cootes et al., “Active shape models-their training and application,” Comp. Vis. Image Understand., vol. 61, no. 1, pp. 38–59, 1995. [13] J. J. Cerrolaza et al., “Segmentation of kidney 3D-ultrasound images using Gabor-based appearance models,” in IEEE Int. Symp. Biomed. Imag., 2014, pp. 633–639. [14] B. van Ginneken et al., “Active shape model segmentation with optimal features,” IEEE Trans. Med. Imag., vol. 21, no. 8, pp. 924–933, Aug. 2002. [15] F. M. Sukno et al., “Active shape models with invariant optimal features: Application to facial analysis,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 29, no. 7, pp. 1105–1117, Jul. 2007. [16] D. E. Gustafson and W. C. Kessel, “Fuzzy clustering with a fuzzy covariance matrix,” in IEEE Conf. Decision Control, 1978, pp. 761–766. [17] Y. G. Tang et al., “Improved validation index for fuzzy clustering,” in Am. Control Conf., vol. 2. 2005, pp. 1120–1125. [18] D. Shen et al., “Segmentation of prostate boundaries from ultrasound images using statistical shape model,” IEEE Trans. Med. Imag., vol. 22, no. 4, pp. 539–551, Apr. 2003. [19] Y. Zhan and D. Shen, “Deformable segmentation of 3-D ultrasound prostate images using statistical texture matching method,” IEEE Trans. Med. Imag., vol. 25, no. 3, pp. 256–272, Mar. 2006. [20] V. Caselles and R. Kimmel, “Geodesic active contours,” Int. J. Comp. Vis., vol. 22, no. 1, pp. 61–79, 1997. [21] P. Kovesi, “Image features from phase congruency,” J. Comput. Vis. Res., vol. 1, no. 3, pp. 1–26, 1999. [22] K. Rajpoot, V. Grau, and J. A. Noble, “Local-phase based 3D boundary detection using monogenic signal and its application to real-time 3-D echocardiography images,” in Proc. Int. Symp. Biomend. Imag., 2009, pp. 783–786. [23] T. F. Chan and L. A. Vese, “Active contour without edges,” IEEE Trans. Image Process., vol. 10, no. 2, pp. 266–277, Feb. 2001. [24] M. Felsberg and G. Sommer, “The monogenic signal,” IEEE Trans. Signal Process., vol. 49, no. 12, pp. 2–8, Dec. 1997. [25] P. Kovesi, “Symmetry and asymmetry from local phase,” in 10th Austral. Joint Conf. Art. Int., 1997, vol. 190, pp. 2–4. [26] H. Edelsbrunner et al., “On the shape of a set of points in the plane,” IEEE Trans. Inf. Theory, vol. 29, no. 4, pp. 551–559, Jul. 1983. [27] J. A. Schnabel et al., “A generic framework for non-rigid registration based on non-uniform multi-level free-form deformations,” in Proc. MICCAI, 2001, pp. 573–581. [28] R. L. Devalois et al., “Spatial frequency selectivity of cells in macaque visual cortex,” Vis. Res., vol. 22, pp. 545–559, 1982. [29] J. J. Cerrolaza et al., “Positive delta detection for alpha shape segmentation of 3D ultrasound images of pathologic kidneys,” in Proc. MICCAI, 2015, pp. 711–718. [30] M. Mulet-Parada and A. Noble, “2D + T acoustic boundary detection in echocardiography,” MedIA, vol. 4, no. 1, pp. 21–30, 2000. [31] K. O. Babalola et al., “An evaluation of four automatic methods of segmenting the subcortical structures in the brain,” Neuroimaging, vol. 47, pp. 1435–1447. [32] T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out. New York: Academic, 2004. [33] M. Berger et al., “State of the art in surface reconstruction from point clouds,” in Eurographics, 2014, pp. 161–185.