Spatio-Temporal Segmentation in 3D Echocardiographic Sequences using Fractional Brownian Motion

12/21/2019 ∙ by Omar S. Al-Kadi, et al. ∙ Yale University The University of Jordan 26

An important aspect for an improved cardiac functional analysis is the accurate segmentation of the left ventricle (LV). A novel approach for fully-automated segmentation of the LV endocardium and epicardium contours is presented. This is mainly based on the natural physical characteristics of the LV shape structure. Both sides of the LV boundaries exhibit natural elliptical curvatures by having details on various scales, i.e. exhibiting fractal-like characteristics. The fractional Brownian motion (fBm), which is a non-stationary stochastic process, integrates well with the stochastic nature of ultrasound echoes. It has the advantage of representing a wide range of non-stationary signals and can quantify statistical local self-similarity throughout the time-sequence ultrasound images. The locally characterized boundaries of the fBm segmented LV were further iteratively refined using global information by means of second-order moments. The method is benchmarked using synthetic 3D+time echocardiographic sequences for normal and different ischemic cardiomyopathy, and results compared with state-of-the-art LV segmentation. Furthermore, the framework was validated against real data from canine cases with expert-defined segmentations and demonstrated improved accuracy. The fBm-based segmentation algorithm is fully automatic and has the potential to be used clinically together with 3D echocardiography for improved cardiovascular disease diagnosis.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 5

page 6

page 7

page 9

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

The automatic segmentation of the left ventricle (LV) of the heart is still considered an open challenge in the field of medical image segmentation. The interest in understanding this largest and main pumping chamber of the heart refers to the key role it plays in blood circulation. Clinicians mainly assess the extent of heart muscle damage by measuring the LV ejection fraction [lai16]. The accurate and reliable segmentation of the LV shape, and being able to early characterize ischemic myocardial damage, is an important prerequisite for further quantitative analysis of cardiac function.

In clinical practice, the segmentation task involves the delineation of LV endocardium and epicardium contours. This process, however, is tedious and time consuming and prone to intra- and inter-observer variability [nbl06]. An automatic and robust approach for segmenting cardiac ultrasound time-sequences would highly facilitate the routine clinical work [frg01]. Several key challenges in the automated segmentation of the LV in cardiac ultrasound datasets exist. Namely, speckle intensity heterogeneity: as in LV cavity (blood pool) due to blood flow or the dynamic motion of the heart; obscureness: close proximity of the papillary muscles tend to show speckle intensities similar to that of the myocardium, and thus affecting endocardium segmentation; spatial complexity of anatomy: the separating border between right and left ventricle, and the low contrast between the myocardium and lung air makes segmentation of the epicardium especially difficult. Other factors as partial volume effects due to limited resolution and inter-variability in shape and speckle intensity of the heart chambers across patients due to pathology may pose additional challenges for LV segmentation. To address these challenges, it is advantageous to have an efficient segmentation algorithm that is objective and reproducible to accelerate and facilitate the process of diagnosis. Regarding ultrasound B-mode segmentation efforts, the endocardium and epicardium boundaries are delineated using a variety of strategies. In particular, statistical models which encode high-level knowledge, as the parametric or geometric deformable models, can handle topological changes robustly [vrg16, bkf16, stg15, cro13]

. However, the boundary finding process requires a model to be initialized sufficiently close to the object to converge and is sometimes prone to local minima. Also, many machine learning approaches were proposed for voxel-based classification

[dmg14, brd16], yet the feature engineering process is not straightforward and computationally expensive. Taking advantage of labelled data for detecting data-driven features has drawn increased attention in prior probabilistic maps [hns14, pdr17]

and deep learning techniques

[cro12, oky18]

. The unsupervised learning by deep networks can reduce the need for feature engineering – one of the most time-consuming parts of machine learning practice, and have recently shown very promising results for improving image classification and segmentation. However, limited training data is a common obstacle in the latter techniques, and regularization of the training data with a large amount of human-annotated data or anatomical models from large datasets is not always available. A review on ultrasound image segmentation methods with techniques mainly focusing on B-mode images can be found in

[nbl06]

. An intuitive approach would be to incorporate the spatio-temporal domain for improving structure and inter-dependencies of the output. Dealing with the LV segmentation problem from a spatio-temporal perspective can give further information on the shape boundaries. Due to the nature of the speckle pattern, it is hard to draw conclusions about the boarder of the LV from still frames. Thus, cardiologists usually examine videos of the deformation of the LV wall during the echocardiographic examination. It is logical to assume the speckle pattern structure is better localized when the spatio-temporal coherence is considered. In this regard, Huang et al. exploits the spatio-temporal coherence of individual data for cardiac contour estimation

[hng14]. Others embarked on introducing temporal consistency in the tracked speckle pattern using optical flow [sdr17]

or gradient vector flow approach

[zgk05]. Nevertheless, the spatio-temporal structures of the speckle patterns are usually stochastic in nature. Neglecting the inherent heterogeneity might not characterize the structure in space and time efficiently. The fractional Brownian motion (fBm), which is a non-stationary stochastic process, integrates well with the stochastic nature of ultrasound echoes [alk15]. It has the advantage of representing a wide range of non-stationary signals and can quantify statistical self-similarity in time-sequence ultrasound images.

To address the aforementioned challenges, we present a novel physically motivated stochastic model for improved epicardium and endocardium boundary segmentation. To our knowledge, this is the first time a spatio-temporal fBm process is used for LV segmentation in 3D ultrasound sequences. Motion roughness, reflected in the speckle pattern, is characterized by fBm for local boundary delineation – which is theoretically invariant to intensity transformations [alk17]. Using second-order moments for LV shape global information complements the local characterization of the fBm process.

The direct delineation of the endocardium boundaries is challenging, even for the trained eye, where poor edge information result in some regions becoming obscured, as in the case of the upper left region of Fig. 1(a). The 3D shape modeling using surface parameterization in Fig. 1(b) naturally reflects the physics of the LV wall motion and clearly highlights the blood pool and its associated myocardium boundary. The similarity between speckle patterns is evaluated by analyzing their statistical behavior, in particular, by measuring voxel-by-voxel the distance between their statistical distributions. Context around each voxel of interest is found and optimized in a pair-wise pattern detection approach. The fBm herein serves as a similarity criterion to identify such patterns. The self-similarity between the endocardium voxel patches – quantified as changes in specular reflections – appear more prominent in the generated fBm parametric image. Similarly applies to the region inside the left ventricular cavity where the speckle patterns are considered more homogeneous. By incorporating the spatio-temporal coherence, an effective segmentation can be achieved. The derived fractal dimensions would reflect the spatio-temporal coherence of the LV boundaries that exhibit natural elliptical curvatures. In this context, the fBm process would adapt to both myocardium and ventricular cavity structures that have sufficient spatio-temporal coherence of varying scale speckle patterns.

In this paper, we introduce a fully-automatic method for robustly segmenting the endocardium and epicardium surface from 3D cardiac ultrasound sequences. In Section II, a 3D stochastic fractal model approach using surface parameterization is adopted for capturing the complex anatomical structure of the LV (see subsection II-B). Both local and global information of the LV shape boundaries are characterized with improved precision (see subsection II-D). Experimental results are reported in Section III. Interpretation and analysis of the results – with suggestions to future directions, and summary of major findings are discussed in Sections IV and V, respectively.

Fig. 1: B-mode cardiac ultrasound image of a canine subject showing (a) the left ventricle, (b) corresponding fractional Brownian motion surface image, and (c) 3D reconstruction and projections in different planes. Individual voxels with green-like color represent radial strain during end-systole.

Ii Methodology

This sections describes a method for improving LV segmentation, namely, endocardium and epicardium boundaries in 3D cardiac ultrasound time series data.

Ii-a Automatic Initialization

Formally, let be the acquisition space-time of the reference time-series of images , , . The fBm segmentation maps , representing the fractal dimension computed for each voxel in each 2D slice (derived in subsection (II-B)), is a combination of a spatial transformation and a temporal transformation . We make the reasonable assumption that the temporal transformation is only time dependent : . The initialization method only requires to store a database with sample patches describing the possible variations found in LV ventricular cavity and myocardium. The left ventricular cavity is divided into a stack of ellipsoidal-like discs (conducted on average images per sequence). A two-step postprocessing stage is applied to local information: the first step is aimed at filling voxel gaps in segmented ventricular cavity, while the second step is aimed at removing falsely detected isolated myocardium voxels, cf. subsection III-C. Techniques such as thresholding, mathematical morphology and correlation are combined for this purpose. Then the method entails the identification of base slice centroid acting as a reference for the fBm segmented LV on 3D echocardiography over time. The centroid is tracked throughout the entire fBm segmented left ventricular cavity, where the local shape and orientation varies as the elliptical model parameters are computed towards the apex. This is iterated twice for the endocardium and epicardium boundaries. The volume of the left ventricular cavity represents the sum of the volumes of each of these disks. The volume of each disc is calculated as the cross-sectional area of the disc multiplied by its height.

Fig. 2: Multiscale parametric mapping based on fractional Brownian motion () for voxel-based segmentation. (a) Object having a spatial support in source image , and (b) constructed fractal parametric voxel map with scaling factor at different scales (), s.t. each parametric value represents a localized Hurst index ().

Ii-B Fractional Brownian Motion

Brownian motion is a Markovian

process, whose conditional transition density function is time-homogeneous. That is, the probability

of being in state at time , given all states up to time , depends only on the previous state, , at time . Therefore, the prediction of what states will occur in the future depends only on the current state. Also, the standard

is a Gaussian process, since it has a normal distribution with specific first two moments. Based on the Lagrangian representation, the generalization of a standard

is a Fractional Brownian motion , which is a continuous Gaussian self-similar process in with stationary increments [MandelbrotNess1968]. can be modeled via stochastic integral equation, given by

(1)

where and are the gamma function and standard Brownian motion, respectively, and is called the parameter or index which describes the scaling behavior of the process [mly12] and the roughness of the resultant motion or trajectory [mdl04]. For our case, characterizes the deformation of the left ventricle wall motion, with lower values leading to a heterogeneous motion and vice versa. From (1), the standard Brownian motion (also often called Wiener process) is recovered when . But in contrast to Brownian motion, fBm has dependent increments when . By allowing to differ from , a fBm process is achieved, where for increments are positive correlated, and for increments are negatively correlated. The subtracted term allows the kernel to vanish quickly when , ensuring the convergence of the integral in the range .

There are several ways to estimate the fractal dimension of a stochastic process modeled by a discrete-time representation of fractional Brownian motion [alk17]. All of them are based on the formula

(2)

where

is proportional to the standard deviation

for samples apart. The fBm and power-law variogram fits were used to estimate as a measure of self-similarity. Particularly, the linearly-related and corresponding fractal dimension can be calculated from the slope of the average absolute difference plotted as a function of the increments with sampling interval (or step-size) on a - plot. Since is a real number with values varying between a random walk process and a smooth process (), the sampling interval can be divided by any arbitrary positive value and the result rescaled in the ratio , cf. property 3 in Appendix V. Then the new semivariogram will be identical to the initial one. In this sense is a self-similar or fractal process. Fig. 2 illustrates the process of constructing a fractal parametric voxel map for an object , representing here the blood pool and associated myocardium boundary, from a source image . The resolution for each image slice is investigated for self-similarity patterns by probing for higher resolutions, i.e. searching for voxel pair structures that exhibit self-similarity at different scales in the range , where represents the maximum probed scale. Then the variogram estimates the mean absolute difference of each voxel pair to the scaling factor , such that the resulting localized Hurst indices , .

Fig. 3: Illustrative example of different steps of multiscale parametric mapping based on fractional Brownian motion () for voxel-based segmentation.

The fractal dimension (FD) of an -dimensional fBm is related to the Hurst index [alk17] by

(3)

where is the Euclidean dimension, or the number of independent variables, and quantifies the self-affinity of the process. Therefore, the closer is to one, the lower the FD, and the smoother the process becomes and vice versa. Hence, quantifying how smooth or regular the cardiac wall motion via fBm might assist in detecting abnormal left ventricular relaxation and increased stiffness related to myocardial ischemia and infarction [thn06].

Finally the regional (or parametric) fractal dimension volume () is generated for each voxel in the time-sequence ultrasound images. In (4), each voxel has its own localized Hurst index () after being subtracted from Euclidean dimension as in (3), and hence a fractal surface would represent the volume. The practical implementation is represented in Algorithm 1.

(4)

The volume elements, namely , and the coordinate index for the third dimension – representing the slice position in the processed volume , are defined for different volumes within the ultrasound time-sequence . The scaling factor is the resolution limits of , which represents the mean absolute intensity difference to center voxels based on local range dependence, cf. Fig. 2. The framework as conducted through voxel-by-voxel is illustrated in Fig. 3.

Input: , I: 3D cardiac ultrasound sequences, in space , and time
1 while  do /* :ultrasound time-sequence */
2       foreach LV volume  do
3             Step 1: 3D shape modeling using surface parameterization
4             forall voxels in cuboid lattice  do
5                   Step 2: Initialize self-similarity stochastic process
6                   forall voxel pair distances in  do
7                         Compute mean absolute difference of each voxel pair ;
8                         Normalize and take the logarithm ; /* , and are size of voxel at certain scale */
9                         Normalize voxel pairs distances ; /* s.t. */
10                     

    Perform least square linear regression as:

,
11                         ;
12                         Step 3: Hurst index H matrix
13                         Estimate localized Hurst indices (): ; /* local linear regression */
14                         Generate fractal map: ; /* for time-sequence and scaling factor */
15                   end forall
16                  
17             end forall
18            Estimate surface heterogeneity from ;
19             Construct feature vector for Volume : ;
20       end foreach
21      
22 end while
23return ;
Output: Fractal parametric voxel map and features vector
Algorithm 1 Fractal parametric voxel map construction
Fig. 4: Endocardium segmentation: (a) B-mode 2D ultrasound image, (b) corresponding fractal parametric voxel map (c) fBm segmented binary image, (d) moment-based endocardium region based on enclosing rectangle and (e) ellipsoidal approximation, respectively.

Ii-C Classification

Patches of mm – representing myocardium vs blood pool – were selected randomly from a single ultrasound volume transformed to , and used for training data. Such that, the training set, , is composed of a set of patches for which the feature vector , and the classification result ( or : myocardium or blood pool) are known

. The Hurst index is a useful statistical method for inferring the properties of a time series without making assumptions about stationarity. It is most useful when used in conjunction with other techniques. Thus, features representing the FD mean and variance values over

, lacunarity () – which defines the sparsity of the fractal pattern in terms of the ratio of the variance over the mean of as in (5) [alk17]:

(5)

and higher order statistics, namely, skewness (asymmetry of the probability distribution) and kurtosis (fourth standardized moment) are defined over localized histograms derived from

and normalized to form a 5–D vector in feature space . Then in the classification procedure, one of the classes (myocardium) or

(blood pool) would be assigned to each candidate voxel when its representation is known. Since our primary concern is to demonstrate the robustness of the new multiscale parametric mapping based on fBm, we describe experiments conducted to compare the discriminative power – which relies naturally on the spatio-temporal dependencies of the local descriptors – for classification. A Bayesian classifier, being a simple probabilistic and commonly used machine learning benchmarking method, which performs well even with possible presence of dependent attributes

[dmg97] was selected. This particularly suits fractal patterns in ultrasound images, where not all locally estimated features are conditionally independent given the class.

The sample training patches were collected from manually labeled LV (performed by 2 experts to avoid errors due to ambiguous interpretation of structures) cavity–blood pool and myocardium voxels in the synthetic and real training images. Specifically, around and voxel samples referring to 62 and 48 voxel patches for the synthetic and animal dataset, respectively, were fairly divided into myocardium and blood pool voxel patches and used for training. A total of 62 and 48 sample voxel patches were used for myocardium training, for synthetic and animal dataset. Similarly applies for the blood pool training. The training voxels represent 0.41% and 0.29% of the total number of synthetic and animal dataset test image voxels, respectively. To reduce the risk of introducing errors, and thus noise in the classification stage, training samples were carefully selected to cover all possible regions related to blood pool, endocardium, epicardium, and myocardium. It is worth noting that the LV segmentation was performed voxel-by-voxel in 3D by employing a multiresolution fBm-based kernel operating as a sliding-window (i.e. cuboid ) in a 3D volume sequence (); revisit Algorithm 1. This would give smoother segmentation in homogeneous tissue regions while preserving fine details in heterogeneous regions, viz. high changes in tissue specular reflections.

Ii-D Ellipsoidal model assignment to fBm maps

For an improved refinement of the LV contour boundaries, each of the fBm segmented regions are fitted by means of image moments to an elliptical model. Moments of images can provide efficient local descriptors and have been used extensively in image analysis applications [kls07, sit14]. Their main advantage is their ability to provide invariant measures of shape, which can better characterize the heterogeneity of the LV wall surface.

Image moments can be defined as weighted averages of voxel intensities. For the time-sequence 3D ultrasound images, at time and depth , the raw -moment for an fBm segmented object (LV for our case) is given by:

(6)

where and are the size of a 2D-image slice of a 3D object and are the labeled binary output of the calculated FD values.

The first-order moments and , when normalized by give the coordinates of the binary object – and of the endocardium or epicardium. Accordingly, second-order moments describe the “distribution of mass” with respect to the coordinate axes and define the orientation of . In order to extract the parameters of the equivalent ellipse from the second-order moments , , and , the central moments can be defined as

(7)

such that

(8)

where and are the coordinates of the centroid of having a size of . These moments are invariants to translation. Then the covariance matrix of the binary object would be:

and the eigen vectors of this covariance matrix correspond to the major and minor axes of the equivalent ellipse.

From (7) the summation extends over all elements in , such that represents the area of the pattern, i.e., the number of white pixels in the middle of Fig .4(c). The coordinates of the centroid can be calculated combining with the image moments of the first degree and . The binary image of the equivalent rectangle111The farthest pixel from the centroid of is specified as the upper-left corner of the rectangle and perpendicular lines are projected in each dimension. Then the smallest bounding box enclosing in is found iteratively. in Fig .4(d) has the same zeroth, first and second moments. Using the moments of the second degree , and , the final formulae giving the major axis orientation and the respective major and minor axis lengths and , are calculated as follows:

(9)
(10)
(11)

Having these parameters we can infer the equivalent ellipse, where was first identified on the fBm segmented LV base image, and used afterwards for alignment – throughout the rest of image sequence – towards the apex. Fig .4 illustrates the above concepts, where and . In this paper, we employ enclosing elements with dimensions .

Iii Experimental Results

Iii-a Experimental Setup

The framework robustness was validated on five different synthetic 3D ultrasound sequences simulating normal and ischemic LV conditions, and furthermore on real 3D cardiac ultrasound sequences acquired from two canine subjects under resting-state and stress-state conditions.

Iii-A1 Synthetic Data

Validation was done on 3D ultrasound time-sequences from the KU Leuven synthetic dataset [asd15], representing one normal and four ischemic cases. Therein, different cases related to LV healthy and pathological conditions were generated. In the mechanical simulations, normal contractility and stiffness values – according to the standard American Heart Association [crq02] – were assigned for the LV normal geometry. The mechanical parameters were tuned in order to match ejection fraction measured on the corresponding template acquisition [asd15]. From the healthy geometry, four ischemic simulations were generated by altering contractility and stiffness in diseased segments. The four ischemic simulations corresponded to: a distal and proximal occlusion of the Left Anterior Descending artery (LADdist and LADprox respectively); occlusions of Right Coronary Artery (RCA) and Left Circumflex (LCX).

Fig. 5: Measuring the discrepancy by means of residual sum of squares – bright regions indicate higher error –between a sample cardiac ultrasound sequence image with a voxel size of mm (shown in Fig. 3(a)) and fBm estimation model with different localized range dependence ( mm, mm, mm, mm, mm, mm, respectively).

Each image in the dataset had voxels of size mm. On average there were 35 images per sequence. For each sequence, ground truth motion trajectories were provided at 2250 mesh points. The endocardium and epicardium meshes were then converted into segmentation masks and used for benchmarking.

Iii-A2 Animal Data

3D ultrasound sequence images were acquired from two acute canine studies (open chested) following a severe occlusion of the left anterior descending coronary artery. A Philips iE33 ultrasound imaging system (Philips Health Care, Andover, MA) with a X7-2 probe at 4.4 MHz suspended in a water bath over the heart was used for image acquisition. Acquisition time points included baseline and one hour and 6 weeks after surgical occlusion of the left anterior descending coronary artery. Images typically had voxels of mm with an average of 23 temporal frames. All experiments were conducted in compliance with the Institutional Animal Care and Use Committee policies.

Iii-B Evaluating Model Goodness-of-Fit

The residual sum of squares (RSS) was used to measure the discrepancy in the estimated indices of the parametric volume maps (). A low value indicates the model has a smaller random error component. Accordingly, the surface of was initially estimated using different fBm localized range dependence by adjusting the range of the scaling factor in (4). Higher RSS was encountered in longer range dependence – i.e. higher resolutions for ; revisit Fig. 3– as the introduced error may get accumulated due to echo artifacts, and hence the estimation becomes unstable. In order to graphically represent the method estimation discrepancy, as this can give a better understanding of the relationship between the fBm model and the surface of the speckle pattern, Fig. 5 shows where errors most likely occur when locally computing , and hence the estimation of . Selecting the best fBm dependence range can avoid unnecessary computational time and give better accuracy.

Fig. 6: Left ventricle endocardium and epicardium initial segmentation result for a sample cardiac ultrasound sequence image by varying the localized fBm dependence range. Best delineation of the endocardium versus epicardium boundary is shown in (c).

Iii-C Patch size vs localized fBm

The local estimation of the fBm process by indicating how far the resolution of can be deeply probed is an important factor for improving segmentation quality. Unlike mathematical fractals where they tend to scale infinitely, the localized fBm dependence range – calculating how deeply the resolution limits of can be probed – is essential for an improved segmentation of LV boundaries. Fig. 6 empirically shows how the values vary at different fBm local dependence ranges, with an indication of the best possible range. To insure consistency, the dependence range has been varied in Fig. 6 in a similar ratio to how the patch size was changed. The best localized range that better delineates the blood pool from the epicardium with the least segmentation errors is Fig. 6(c). An example of a LV base segmentation after undergoing fBm stochastic modelling and ellipse fitting for boundary refinement is illustrated in Fig. 7.

Fig. 7: Segmentation of the left ventricle of the heart from (a) 2D time sequence ultrasound image demonstrating in (b) and (c) epicardium ground-truth vs automatic fBm segmentation, (e) and (f) endocardium ground-truth vs automatic fBm segmentation, respectively. A joint comparison between ground-truth vs automatic left ventricle fBm segmentation is shown in (d), where white color means perfect match, while magenta and green refer to dissimilarity of the ground truth vs automatic fBm segmentations, respectively. (For qualitative interpretation of segmentation quality, figures (b)-(f) best viewed in color)

Iii-D Quantitative Evaluation

The proposed fBm segmentation method was benchmarked using state-of-the-art method in [hng14], which is based on employing spatio-temporal coherence under a dynamic appearance model, abbreviated herein as C-DAM. The evaluation of the segmentation quality was performed using three different segmentation measures, which are: Dice coefficient (DC), Hausdorff distance (HD), and mean absolute distance (MAD), and computed for both end-diastolic and -systolic and averaged over all cases, see Table I and II. Fig. 8 shows a slice-by-slice LV epicardium segmentation quality evaluation corresponding to the DICE metric in Table II. Also the results of the resting/stress state of the canine dataset showed an improvement of , mm, mm and , mm, mm for the DC, HD, and MAD segmentation quality metrics of the epicardium and endocardium borders for the fBm-based as compared to the C-DAM method, respectively.

For clinical relevance, as numerical measures tend to be more narrowly focused on a particular aspect of the data and often try to compress information into a single number, the LV volume segmentation quality is represented as well graphically in Fig. 9 and Fig. 10 for the synthetic and animal dataset, respectively. Results show that the method in [hng14] tends to over-segment the epicardium boarder by extending into the surrounding area of myocardium (i.e. LV appearing dilated as compared to ground truth), and the epicardium surface is more irregular as compared to the proposed fBm-based segmentation method, cf. Fig. 9(f) and (k). Also, the mid-ventricular cavity shows improved delineation in the complex canine dataset using the fBm-based segmentation method, cf. Fig. 10(g)-(i). Furthermore, the absolute difference is measured between the different voxel patches in search for self-similarity properties, i.e. invariance under a suitable change of scale. It describes the different statistical distribution variations found in the LV ventricular cavity and myocardium. To this end, the estimated fBm local range dependence serves in determining how far the image resolution can be deeply probed. Terminating the search at the optimized scaling factor contributes in improving segmentation of the LV boundaries, and further saves unnecessary computational time. The implemented fBm algorithm in this paper had a computational complexity of , with running time nearly similar for both methods, i.e. around 1 minute per frame. However, the proposed fBm segmentation method is fully automatic and does not rely on training data.

Left Ventricle
Condition
Method

Segmentation Evaluation Metric

DICE (%) HD (mm) MAD (mm)
Normal
C-DAM
Our fBm
RCA
C-DAM
Our fBm
LCX
C-DAM
Our fBm
LADdist
C-DAM
Our fBm
LADprox
C-DAM
Our fBm
TABLE I: Endocardium segmentation for normal and ischemic cardiomyopathy
Left Ventricle
Condition
Method Segmentation Evaluation Metric
DICE (%) HD (mm) MAD (mm)
Normal
C-DAM
Our fBm
RCA
C-DAM
Our fBm
LCX
C-DAM
Our fBm
LADdist
C-DAM
Our fBm
LADprox
C-DAM
Our fBm
TABLE II: Epicardium segmentation for normal and ischemic cardiomyopathy
Fig. 8: Comparison of slice-by-slice dice coefficient segmentation quality for fBm and C-DAM [hng14] in (a) normal, (b) RCA, (c) LCX, (d) LADdist, and (e) LADprox left ventricle conditions.

Iv Discussion

The fBm process is a useful stochastic method for inferring the properties of a time series without making assumptions about stationarity. It relates to how strong the autocorrelations of the spatio-temporal coherence, and the rate at which these decrease as the lag between pairs of scatterer patterns increases. Although it might seem that searching for self-similarity properties by probing for higher resolutions, i.e. selecting higher values for in (4), could improve the segmentation quality, results show that reliable estimate of , and hence the , is valid only at a certain cutoff scale where there can be no more details. That is, the echo patterns being approximate rather than deterministic, as the characteristics of the pattern tends to scale in a statistical fashion. The self-similarity property in this sense means invariance in distribution under a suitable change of scale . Therefor the local linear regression relation for estimating (revisit Algorithm 1) becomes non-linear at higher values of – due to exceeding the actual resolution of the ultrasound image itself [alk08], and eventually resulting in error accumulation. On the other hand, complementing the self-similarity property of the fBm process with the lacunarity measure can give a better characterization of the LV shape heterogeneity. Lacunarity analysis is a technique introduced to deal with fractal objects of the same dimension with different textural appearances [alk09]. The lacunarity parameter describes the local complexity of the speckle pattern, on different scales, based on the spatial distribution of gaps of a specific size. Namely, it is a measures of the sparsity of the fBm process, providing additional information on the how irregularity fills the space.

From a clinical perspective, one of the key properties of the fBm process is that it can exhibit persistence () or anti-persistence (). Persistence is the property that the LV wall motion tends to be smooth, e.g. near to normal. Anti-persistence is the property that the relative stochastic process is very noisy, and hence LV wall motion trajectories tend to be heterogeneous. The latter case is an example of cardiac ischemia, such that displacements over one temporal or spatial interval are partially cancelled out by displacements over another time interval. From Table I, the fBm segmentation method showed improved performance for nearly all tested LV endocardium and epicardium deformation conditions except for the ischemic–LADdist condition. The occlusion of the distal and side branches could have a negative side-effect on LV myocardial function, introducing heterogeneity in the LV wall motion velocity pattern [azv96]. Therefore the spatio-temporal structure of speckle patterns could be obscured or difficult to discern; affecting the fBm stochastic modelling for LADdist segmentation. Besides the possible presence of extensive epicardial fat may also complicate the LV segmentation. The endocardium of ischemic–LADprox showed nearly equivalent performance in both methods. For practical relevance, results are also demonstrated visually in Fig. 9 and Fig. 10. Improved delineation of the epicardium and endocardium boundaries can be seen, especially for the LV base slice.

Fig. 9: 3D segmentation volumes of the left ventricle for 3D echocardiographic sequences illustrating [upper-row] ground truth for normal, ischemic-RCA, ischemic-LCX, ischemic-LADdist, and ischemic-LADprox; [mid-row] C-DAM method in [hng14]; [bottom-row] proposed fully-automated fBm approach. (Video examples can be referred to in supplementary materials)

In the ellipsoidal modelling phase, where the fBm segmented image is combined with the corresponding shape information, the relative difference in magnitude of the eigenvalues

and are an indication of the eccentricity of the LV. This gives an elliptical ventricle with a cone-shaped apex rather than having a round cylindrical-shape. It is known that the contraction of the LV cavity is less symmetrical because of the systolic increase in wall thickness [kds95]. A way to view the LV is to consider it as a composite of adjoining structures; therefore a circular or elliptical shape at one point may not represent its entire structure. This can be attributed to the heterogeneity between the major axis diameter and cross-sectional area for the different regions of the LV during contraction, which is associated with the descent of the base and rotation of the apex. As the ventricular shape changes from an elliptical to a more spherical form, the assumption of a uniform structure of LV that is localized spatially would be misleading. The 3D spatio-temporal imaging would be advantageous for assessing the LV shape and size change at specific regions of interests. The purpose of the fBm-based method is to account for changes in shape from circular to elliptical as shown in Fig. 9, where the difference in the major axis diameter from the apex to the base of the LV is refined by image moments. Nevertheless, other geometric primitives could be used, as well as deformable surfaces for improving LV boundary refinement of the fBm-based segmentation method.

Methods relying on evolving surfaces, e.g. [hng14], can allow for flexible topology changes and does not assume a priori knowledge of object’s shape; however, they cannot effectively segment surfaces that break apart or intersect – which is a common condition in ultrasound imaging. Ultrasound image segmentation is challenging due to the inherent speckle and presence of artifacts such as shadows, attenuation and signal dropout. This often leads to missing edges, making it difficult for such methods to deal with structure discontinuity. Results in this work support that irregular structure characterization by locally investigating self-similarity patterns appearing at different scales using the fBm process can better quantify the heterogeneity in LV wall motion. The improved localization in both the temporal and spatial domains gives the fBm approach the advantage of accurate segmentation for fine details, and possible compensation of edge disconnection. Moreover, the non-stationarity of the fBm process integrates well with the stochastic nature of ultrasound echoes.

Finally, due to the nature of how the ultrasound volume is segmented – i.e. slice-based approach, the current implementation is not free of several limitations: 1) A disc-like structure is usually assumed in magnetic resonance (MR) data due to the imaging limitations for cine MR images. However, this approach might hamper the segmentation of a closed myocardium – since true 3D ultrasound segmentation algorithms as in [brd16] consider the LV as a closed structure, although not always assumed as such for the basal side, but typically for the apex; 2) Previous work has reported that different 3D echocardiography software packages show variability when used in predicting response to cardiac resynchronization therapy [aly12], therefore comparing performance with similar commercial software packages would be interesting to investigate; 3) Performing the fBm segmentation based on the radio-frequency envelope detected echos may better reflect the intrinsic properties of the myocardium tissue [alk16, alk16b], and assist in overcoming issues related to ultrasound B-mode settings.

Fig. 10: [Left] 3D segmentation of left ventricle for 2 canine subjects illustrating (a)-(c) manual, C-DAM, fBm in resting state condition, and (d)-(f) manual, C-DAM, fBm in stress state condition; [Right] comparison of endocardium segmentation accuracy of LV mid-ventricular cavity in resting state for (g)-(i) manual, fBm-based, C-DAM, respectively. (For qualitative interpretation of segmentation quality, figures (g)-(i) best viewed in color)

V Conclusion

This work presents a novel method for improving LV segmentation by addressing the problem of speckle pattern heterogeneity, where a) the fBm classification-based segmentation method relies naturally on the spatio-temporal dependencies of the local features; b) the 3D sequence of Hurst indices, which was used to derive fractal dimension volume maps, are invariant to intensity transformations; c) global information about the LV shape using second-order moments complements the local characterization of the fBm process. Both local and global boundary information about the LV shape boundaries was captured with improved precision.

-a Covariance of fractional Brownian motion

Given a Gaussian process that is characterized by associated parameter , then it follows that the covariance function is given by

(12)

for , where denotes the expectation operator with respect to probability space, and is the Euclidean norm of .

-B Discrete fractional Brownian motion

A discrete-time representation of fBm can be obtained by sampling the continuous-time fBm. When approximating (1) by sums, the first integral should be truncated, say at -. The approximation is for given by

(13)

where with respect to are mutually independent vectors.

-C Fractal dimension estimation

There are several ways to estimate the fractal dimension of a stochastic process modeled by a fBm [alk17]. All of them are based on the formula

(14)

for the variogram of fBm, which follows from (-A). The fBm and power-law variogram fits were used to estimate as a measure of self-similarity. For the discrete-time fBm process defined in (-B), the following properties hold [dch97]:

Property 1: The mean of the fBm increments, which are samples apart, is zero, i.e.,

(15)

The covariance of the fBm increments, samples apart, is given by

(16)

When and , the variance of the unit increments is , which is the variance of independent identically distributed samples. From (16), the following two properties of fBm can be obtained.

Property 2: For unit increments , the covariance of the increments becomes

Property 3: The variance of the increments when is,

The self-similarity characteristic of fBm can be recognized in property 3 with a scale transformation

where . Such a result indicates that fBm is statistically indistinguishable under a scale transformation and implies scale-invariance. Combining property 1 and 3 to obtain

(17)

which is generally expressed in the form

(18)

where is proportional to the standard deviation . Taking the logarithm of (18) yields the linear equation

(19)

Acknowledgment

The author would like to thank Dr. Albert Sinusas’s lab for provision of the canine dataset, Dr. James Duncan for very helpful discussion and making available the C-DAM method for comparison, and Allen Lu and Dr. Nripesh Parajuli for assisting with the delineation of the ground truth. This work was supported by Fulbright Scholarship Grant No. G-1-00005.

References