Quantification of ultrasonic texture intra-heterogeneity via volumetric stochastic modeling for tissue characterization.

Intensity variations in image texture can provide powerful quantitative information about physical properties of biological tissue. However, tissue patterns can vary according to the utilized imaging system and are intrinsically correlated to the scale of analysis. In the case of ultrasound, the Nakagami distribution is a general model of the ultrasonic backscattering envelope under various scattering conditions and densities where it can be employed for characterizing image texture, but the subtle intra-heterogeneities within a given mass are difficult to capture via this model as it works at a single spatial scale. This paper proposes a locally adaptive 3D multi-resolution Nakagami-based fractal feature descriptor that extends Nakagami-based texture analysis to accommodate subtle speckle spatial frequency tissue intensity variability in volumetric scans. Local textural fractal descriptors - which are invariant to affine intensity changes - are extracted from volumetric patches at different spatial resolutions from voxel lattice-based generated shape and scale Nakagami parameters. Using ultrasound radio-frequency datasets we found that after applying an adaptive fractal decomposition label transfer approach on top of the generated Nakagami voxels, tissue characterization results were superior to the state of art. Experimental results on real 3D ultrasonic pre-clinical and clinical datasets suggest that describing tumor intra-heterogeneity via this descriptor may facilitate improved prediction of therapy response and disease characterization.

Introduction

Analysis of the local characteristic patterns of tissue texture can reveal subtle pathological features deemed important for clinical diagnosis. Spatial variation of textons quantified in terms of image “surface roughness” has been shown to reflect tumor functional heterogeneity, and to lead to a better understanding of disease state (Bae et al., 2013, Chicklore et al., 2013, Davnall et al., 2012, Al-Kadi and Watson, 2008). However, sub-voxel resolution complex and higher order textural features can be difficult to discern by simple observation. These texture signatures may convey significant information about disease progression or regression. However, quantifying these subtle signatures in ultrasound images is challenging.

Our motivation stems from a clinical need to improve the diagnosis and therapy of liver cancer. Approximately 100,000 patients are diagnosed each year with primary liver cancers in the United States and Europe (Cancer Research UK, 2014, American Cancer Society, 2014). When this is compared against worldwide statistics, liver cancer is even more common in developing countries (World Health Organization: International Agency for Research on Cancer, 2012). Although it the sixth most common cancer in the world (Ferlay et al., 2013), incidence varies across the world, and it is the most cancer type in some developing countries (Parkin et al., 2014). Surgery is considered the only curative treatment; however, this is not suitable in the majority of cases due to co-morbidity, extent or location of the cancer, with chemotherapy forming the mainstay of treatment in these patients. Chemotherapy can have significant side effects, and may not be effective in all cases. Development of monitoring techniques during the course of chemotherapy may permit dose adjustment in responders to minimize side effects, while alternative treatments could be offered to non-responders. Current monitoring techniques rely on computed tomography and magnetic resonance imaging, with frequency limited by the potential damage from ionizing radiation and cost consideration respectively. Despite the difficulties of using ultrasound for monitoring disease (e.g. operator dependent, poorly reproducible and non-standardize), it is a technique that is known to be rapid, relatively inexpensive, readily available, with no exposure to ionizing radiation, making it ideal for frequent monitoring of liver tumors during a course of treatment.

Given the advantages of ultrasound, analyzing tissue speckle from a single resolution perspective is limiting, as substantial information that could assist tumor tissue characterization can be hidden at sub-voxel resolution. This is true for the smaller necrotic or functionally low-activity regions that exhibit a hyperechogenic appearance compared to healthy tissue (Ueta et al., 2011). We hypothesize that the difference in echogenicity of the tumor speckle texture can be exploited as an indicator of disease responsiveness to treatment (Czarnota et al., 2013). Nevertheless, the functionally low-activity regions within the tumor texture are relatively small, especially in the early sessions of tumor chemotherapy treatment. They also tend to have low intensity contrast compared to the aggressive or functionally active background of the remaining tumor. Identification of subtle changes in these regions based on visual assessment of the intensity alone can be challenging.

Our approach is motivated by four observations:

Tumors are heterogeneous: most previous work has accounted for functionally active malignant regions rather than the peripheral low activity necrotic regions which may additionally provide key information on disease progression or regression. These subtle variations and deviations within the speckle tissue texture were deemed too chaotic to be characterized in Larrue and Noble (2014), but are important for understanding disease state.

Heterogeneity suggests using a multi-resolution texture analysis: a carefully designed multi-resolution approach which is visually discriminative and geometrically informative could reveal small speckle changes and is better suited to describe the mixture distribution complexity that underpins a heterogeneity model.

Fractal analysis is well-suited to this problem: conventional energy-based wavelet decompositions are susceptible to local intensity distribution variations; the fractal signatures used in our approach, derived from the wavelet representation sub-bands related to physiological properties of texture surface roughness are not; finally.

Analysis should be three-dimensional: performing a 3D texture analysis based on a volumetric Nakagami modeling could facilitate a more reliable estimate of the Nakagami parameters, where the 3D location of each voxel provides a better localization of speckle distribution mixtures.

In this work, a novel multifractal Nakagami-based volumetric feature descriptor that is invariant to local speckle attenuation changes is proposed. A pipeline summarizing the stages of our approach is illustrated in Fig. 1. It is postulated that fractal tissue characteristics locally derived from 3D textural tumor patterns at several scales and from the RF envelope of the ultrasound backscattered volumes can assist in attaining descriptive features that relate to underlying biological structure. These tissue textural fractal characteristics tend to change in cases of therapeutic response, providing an attractive indicator for disease response to treatment during chemotherapy.

Fig. 1

3D multifractal Nakagami feature descriptor algorithm design for ultrasonic tissue characterization introduced in this paper.

This paper is organized as follows. State of the art and challenges associated with characterizing speckle tissue texture heterogeneity are summarized in Sections 2, 3, followed by a detailed explanation of the proposed 3D multifractal Nakagami-based feature descriptor in Section 4. Sections 5, 6 present the experimental results and discuss the potential significance of the work. The paper concludes in Section 7.

Related work

One of the effective statistical techniques used for modeling various backscattering conditions in medical ultrasound is the Nakagami distribution. This probabilistic distribution is known for its analytical simplicity and effective modeling of dense scatterers, accounting for amplitude and spacing, and can be reduced to a Rayleigh distribution under certain assumptions of scatterer density and number (Shankar, 2000). Shankar et al. first proposed the Nakagami distribution for characterizing conditions ranging from pre- to post-Rayleigh existing in ultrasound images, and later for modeling the radio-frequency (RF) envelope of the ultrasound backscattered signal in characterizing B-mode breast masses (Shankar et al., 2001). Others have attempted to tackle the issue of accurate estimation of the Nakagami distribution. For instance (Larrue and Noble, 2011) employed Gamma kernel density estimation to compute a smooth estimate of a distribution within small windows of B-mode ultrasound images, but the mixture of distributions occurring at the boundaries between structures was not accounted for. The impact of morphological parameters and tumor structures on the Nakagami parameters statistics were analyzed in Larrue and Noble (2014). A limitation observed was that there was a need for a robust algorithm to compute the Nakagami parametric images that better delineate the structures and the context in and around the tumor. Characterizing homogeneous tissues via improving the smoothness of the Nakagami parametric images was shown in Tsui et al. (2014). The technique relies on summing and averaging the Nakagami images formed using sliding windows with varying window sizes related to the transducer pulse length. However, a relatively large window size (required for stabilization) may affect the reliability of the estimated Nakagami parameters, and hence degrade the spatial resolution of the resulting Nakagami image.

The Nakagami distribution has been further employed as an image feature in various image analysis contexts. For instance, five contour features and the two Nakagami parameters were used for classification of benign and malignant breast tumors in Tsui et al. (2010). In a subsequent work malignant tumors were shown to be more pre-Rayleigh distributed than those from benign counterparts (Tsui et al., 2010); however, the calculation of the average intensity value in the Nakagami image makes it susceptible to spatial frequency intensity variability. Further, that particular technique was optimized for 2D ultrasound images which may not reliably represent heterogeneous distributions of scatterers (or speckle) encountered within a tumor volume. A random forest based solution to learn tissue-specific ultrasonic backscattering and a signal confidence for predicting heterogeneous composition in atherosclerotic plaques was proposed in Sheet et al. (2014). That technique was developed for intravascular ultrasound and risk assessment of plaque rupture (Zhou et al., 2002). Necrotic core was not considered in that method. Finally, Bouhlel and Sevestre-Ghalila, 2009, Klein et al., 2011 describe a Markov random field model combined with Nakagami distribution estimation to differentiate malignant melanoma from normal tissue. However it was found that the estimated scale model parameter was highly sensitive to image quality, and hence subtle variations could go unnoticed. For an overview of ultrasound tissue characterization we refer the reader to Noble (2010). Previous work on ultrasound texture analysis of tumors has considered both global and local non-uniformity quantification of the tumor texture at only a single analysis scale. Herein we are primarily concerned with tumor intra-heterogeneity (i.e. micro-structures within the tumor speckle texture) which is more challenging.

Challenges in ultrasonic speckle texture characterization

Speckle is a granular-shape stochastic pattern which appears in an image resulting from the scattering of an RF incident signal on an object (Sanches et al., 2011). The spacing and localization of the scatterers in the scanned object structure contribute to the local variation and distribution of the recorded texture pattern. However the characteristic interference patterns, known as speckle, produce an overall reduction in global image contrast (Noble, 2010). As a consequence, the boundaries separating different structures are less well defined, increasing difficulty in delineating regions of interest with a resultant increase in inter-and intra-observer variability for tumor detection.

A means to mitigate against effects such as the beam-tissue physical interaction and other acquisition factors is to characterize the objects via their speckle textural properties (De Grandi et al., 2003, Madabhushi and Metaxas, 2003, Sadeghi-Naini et al., 2013). Textons or texels (texture elements) which are the fundamental components of texture that collectively form the observed speckle pattern texture do not directly correspond to the underlying structure; however, the local intensity textural pattern can reflect the local echogenicity of the underlying scatterers (Anderson and Trahey, 2006), see Fig. 2. This is due to the stochastic nature of the speckle pattern. Viewing the structure locally as a collective texton structure can give information about the underlying scatterer behavior. We hypothesize here this may lead to an improvement in internal structure delineation, and hence tumor characterization.

Fig. 2

Six ultrasound hypoechoic to hyperechoic gray scale target phantoms having 8 mm diameter and 4 cm depth and corresponding simulated B-mode image representing a varying intensity from hypoechoic, −6, −3, +3, +6 dB, and hyperechoic, respectively.

Fig. 3(a) shows a simulated lesion phantom having three different sizes at a fixed depth and with four levels of intensity contrast variability. The small round simulated hyperechoic region (marked with an arrow) resembles the functionally low-activity regions on real B-mode ultrasound images; appearing subtle and challenging to identify. Note that the situation would be even more complicated in real B-mode ultrasound tissue where the functionally active background of the aggressive tumor tissue would not be as uniform as in this example, and the non-aggressive regions do not usually have a constant intensity distribution, see Fig. 3(b).

Fig. 3

(a) Simulated ultrasound B-mode image following the method in Bamber and Dickinson (1980) showing different 4 cm depth of 4, 6 and 8 mm diameter gray scale target phantoms ranging from −6, −3, +3 and +6 dB varying intensity, (b) a real ultrasound B-mode volume of interest of a liver tumor with corresponding fractal slice map in (c) – estimated from the RF envelope of the ultrasound backscattered signal – indicating the subtle low-activity regions.

The tissue characterization solution, as discussed in the paper, is to use a multi-resolution approach that highlights higher order statistical features of the RF envelope. Such features could go unnoticed in B-mode ultrasound and an experienced observer could struggle to identify subtle interval changes in these important texture features.

Fractals and wavelet packet analysis provide effective ways to break down statistical complexity to distinguish between different texture regions, where the invariance to affine speckle intensity changes for the former and the high sensitivity to local features for the latter facilitates effective texture discrimination (Al-Kadi, 2009). Furthermore, according to the uncertainty principle, the wavelet packets can achieve an optimal joint spatial-frequency localization, i.e. simultaneously maintain a good boundary accuracy and frequency response (Rangayyan, 2005), and the estimated fractal dimension can give a quantitative assessment of the surface roughness (Al-Kadi and Watson, 2008, Lopes and Betrouni, 2009). Finally, simultaneous macro and micro scale tumor texture analysis provides a more complete characterization of dense and sparse textural regions within a tumor volume of interest. As demonstrated later, this progressive refinement process optimizes characterization by giving a better fit to the underlying tumor speckle texture.

Methodology

Our goal is to derive a locally-based feature signature based on volumetric generated Nakagami shape and scale parametric voxel lattices, and subsequently to perform an intensity-invariant texture analysis at various spatial resolutions for tissue characterization. This allows us to perform a more complete characterization of tumor texture at the optimal resolution scales compared to single or mono-resolution approaches (Al-Kadi, 2009). The proposed volumetric dense-to-sparse approach can break-down the speckle complexity and provide a robust estimation of model parameters, while having the advantage of simultaneously localizing both large high-contrast and small low-contrast structures at low and high spatial resolution levels.

Nakagami probabilistic distribution

The Nakagami distribution N(x) is an analytically simple distribution that has been proposed as a general model for the ultrasonic backscattered envelope under all scattering conditions and scatterer densities (Shankar, 2000). This distribution has the density function:where x is the envelope of the RF signal, with the shape of the distribution defined by the μ parameter corresponding to the local concentration of scatterers, and the local backscattered energy represented by the scale parameter ω 〉 0, for x > 0, and is the Gamma function. If x has a Nakagami distribution N with parameters μ and ω, then has a Gamma distribution Γ with shape μ and scale (energy) parameter ω/μ.

The Nakagami distribution can model various backscattering conditions in medical ultrasound. By varying μ, the envelope statistics range from pre-Rayleigh , Rayleigh , and to post-Rayleigh . The Nakagami parameters are generally estimated by the 2nd and 4th order moments, where given x is the ultrasonic backscattered envelope and denotes the statistical mean, the two Nakagami parameters can be calculated as:

Volumetric multi-scale Nakagami modeling

A 3D feature signature that operates locally is defined by having each volume V consisting of z acquired slices subdivided into voxel lattices , each having a defined size of m and n, where for , such that . For each we assume that for a scaling factor r at a specific spatial scale s, the scaled voxel intensity lattice values of the RF envelope amplitude such that , where the different possible resolution levels reaching to the maximum level j represent a stochastic pattern, and the envelope amplitude of the scales r of follows a Nakagami distribution. Given the large number of voxel samples to analyze and the known family of probability distributions, the maximum likelihood estimators would tend to have a higher probability of being close to the quantities to be estimated and more often unbiased as compared to moments-based estimation (Cheng and Beaulieu, 2001), therefore the associated shape and scale parameters were estimated via maximum likelihood estimation (MLE) by operating on each voxel lattice region and at different scales. The maximum likelihood estimate for a density function when θ is a vector of parameters for the Nakagami distribution family Θ, estimates the most probable parameters , where is the score function. Having generated voxel-based Nakagami parameters, 3D wavelet packet Daubechies analysis (Mallat, 1999) can be applied at multiple scales. Namely:where is relative to scaling and wavelet function and defines an approximation of at scale , and coefficients add horizontal, vertical and diagonal details for scales . The Daubechies wavelet family can account for self-similarity and signal discontinuities, making it one of the most useful wavelets for characterizing signals exhibiting fractal patterns (Daubechies, 1990). In our case an orthogonal 8-tap Daubechies filter was used to obtain the wavelet packets by expanding the basis having the most significant fractal signature rather than energy.

An octant wavelet transform depends mainly on the scaling and wavelet filters for image decomposition, and one does not need to express the and in their explicit form. The decomposition process can be viewed as passing the signal through a pair of lowpass and highpass filters, also known as quadrature mirror filters, having impulse responses and , while holding the size of the transformed image the same as the original image as we are applying an overcomplete wavelet representation; hence giving a better representation of the texture characteristics at each decomposition. The impulse responses of L and H are defined as and for scaling parameter a, and and for translation parameter b, where . The decomposition is performed recursively on the output of and . Hence, the 3-D wavelet (or octant wavelet packet) can be expressed by the tensor product of the wavelet basis functions along the horizontal, vertical and depth directions. The corresponding filter coefficients can be recursively decomposed by a factor of eight as illustrated in Fig. 4 and expressed in (5), with subscripts indicating the low and high pass filtering characteristics in the and j directions:By decomposing the approximation coefficients of the signal as well, the wavelet transform can be extended in the middle and high frequency channels, providing a more complete partitioning of the spatial-frequency domain, which is known as the octant wavelet packet transform (Coifman and Wickerhauser, 1992). As the textural information about the structural arrangement of surfaces and their relationship to the surrounding neighborhood is spread across the frequency sub-bands, most of the important discriminant features related to the structure terminations and endpoints of surface edges will have a stronger response in higher frequencies (Al-Kadi, 2009). Thereby, this gives an equal opportunity for investigating descriptors of textural features prevailing in the middle and high frequency bands.

Fig. 4

Multiresolution volumetric modeling showing the decomposition up to 3 hierarchical levels by recursive subdivision of volume into octants voxels (left) and corresponding decomposition tree (middle).

From a pattern recognition perspective, the selection of the most suitable wavelet is associated with the understanding of the tissue textural properties and synthesis wavelet. Wavelet analysis using Daubechies wavelet basis functions can achieve a good spatial-frequency localization by having narrow high and wide low frequencies simultaneously. With the increasing number of zero or vanishing moments – which are half the number of filter taps N – this can give a sparse representation for a large class of signal types. Also the Daubechies orthogonal wavelet family consists of purposefully designed filters which account for self-similarity and signal discontinuities, making them one of the most useful wavelets for characterizing signals exhibiting fractal patterns. Besides, they are also considered to be sensitive in recognizing fine characteristic structures, and its application of overlapping windows, unlike other wavelets such as the Haar wavelet, facilitates the capture of all high frequency changes easily (Mallat, 1999). As our work is concerned with the estimation of texture surface roughness from a fractal dimension perspective, the choice of this wavelet is more suitable than other wavelet families (Daubechies, 1992). Therefore an orthogonal 8-tap Daubechies filter (Daubechies, 1990) in a tree structure decomposition is used to obtain the wavelet packets by expanding the basis having the most significant fractal signature, see Fig. 5. This approach gives flexibility to finely tune the signal to the characteristic intrinsic properties of an image (Wang and Yong, 2008).

Fig. 5

Normalized Daubechies’ orthogonal wavelet showing scaling (father) and wavelet (mother) functions with 4 vanishing moments, and corresponding first level Daubechies wavelet decomposition for a liver tumor volume of interest showing from left to right the approximation, horizontal, vertical, and diagonal coefficients, respectively.

Multi-fractal textural model

Fractals can be used in tissue characterization to describe irregular structures that exhibit semi self-similarity at different scales, and can further give an estimation of surface roughness (in our case of the RF envelope surface). There are several fractal models used to estimate the fractal dimension (FD); the FD can be estimated via the fractal Brownian motion (fBm) defined in (6) below, which is a non-stationary model known for its capability for describing random phenomena (Lopes and Betrouni, 2009). Its statistical invariance to dilation, translation and rotation, can mitigate multiplicative speckle scale changes, making it a perfect candidate to be integrated with the Nakagami modeling and multi-resolution decomposition:where is the mean absolute difference of voxel pairs where n = 3 for 3-D space, is the voxel pair distances; H is called the Hurst coefficient; and the constant K 〉 0.

Fractal map estimation

After application of the Daubechies 3D wavelet analysis, the roughness of each voxel lattice surface is determined via estimating its corresponding FD. The estimated voxel-by-voxel array of fractal dimensions for each voxel lattice, which we call a fractal map, provides a basis for characterizing the tissue and for building a bag-of-words of fractal features as a 3D feature descriptor.

A multi-dimensional matrix defined for each of the tumor voxels is derived at different range scales r, such that the mean absolute difference of each voxel pair and for each voxel pair distances are estimated. Thereby the first dimension d represents the voxel after it has been scaled once, and the second dimension represents the voxel at scale 2, and so on until the highest scale j is reached.where M and N are the size of each ultrasound image slice and is the resolution limits of matrix which represents the mean absolute intensity difference to center voxels, and i stands for Nakagami shape μ and scale ω parametric images. Then each element from each array in is normalized after taking the logarithm and saved in a mean absolute difference row vector . That is, the first element in all arrays of will compose vector , and all second elements will compose vector , and so on as shown in (8). This process is illustrated in Fig. 6.

Fig. 6

Overcomplete multi-scale volumetric Nakagami tumor regions; the small voxel lattice centered on the localized voxel , where , and r are the voxel position on the lattice for scale r, is convolved with larger voxels up to j resolution levels for estimation of the mean absolute difference of voxel pairs matrix.

The slope – which corresponds to the Hurst coefficient H – of the least square linear regression line of the log–log plot of versus can be determined by means of sums of squares as in (9).Finally, the slope of the linear regression line defines the textural fractal characteristics, which we call the fractal map :

Volume of interest refinement

It is important to estimate the Nakagami model parameters with good accuracy, but still have a simple model that is easy to interpret. Estimation from small cuboids of interest can provide poor estimation of the Nakagami parameters (Larrue and Noble, 2014). Larger volumes have more data points to fit allowing for averaging of random error, yet this might not be good for tumors with relatively small size. In order to balance the trade off, volume reconstruction was designed to eliminate 2D tumor slices with low information content in order to provide a good characterization of tumor textural patterns.

In practice, as a tumor grows it tends to adopt a non-uniform shape. This will cause sections in the acquired tumor volume to have a relatively small area compared to the whole tumor, where the texture patterns within these small regions cannot be reliably extracted. Therefore removal of these small regions will not only assist in reducing irrelevant features, computational time and memory, but will also direct the efforts of the developed feature descriptor to focus on characterizing the tumor patterns provided in large volume. The selection of volume slices was performed such that , where is the tumor area in slice , and is the slice m with the median area size.

Another important design decision is the selection of the size of the lattice utilized in Nakagami distribution estimation which is ideally performed automatically. To address this, a varying size voxel lattice was introduced as illustrated in Fig. 7 to measure the goodness of fit for the estimated Nakagami model parameters. For error of fit we used the estimated root mean square error (RMSE) between the MLE-estimated Nakagami values and the observed voxels values in each voxel lattice starting at size 2, as there would be no meaning if a lattice had a size of 1. The RMSE was thus defined as:

Fig. 7

Segmented liver tumor volume of interest; and the annotation squares in the enlarged image show the variation of voxel lattice size used in the experiments.

Fig. 8 is an empirical plot of the goodness of fit of the estimated Nakagami parameters versus voxel lattice size for a typical dataset. Sizes varying from 0.03 mm3 to 6.60 mm3 where used in the experiments. The RMSE oscillates as it reaches its minimum, recording a residual error of 0.68 at a lattice size of 7 voxels before the accuracy starts to decrease for larger sizes. Also in the process of generating the Nakagami parametric images and when the voxel lattice happen to be on the border, all voxels laying outside the volume of interest are eliminated from the calculations in order to discard any bias and to maintain a more credible estimation of the parameters.

Fig. 8

Goodness of fit for optimizing the voxel lattice size utilized in the Nakagami distribution fitting.

Feature selection optimization

The ultrasound texture fractal maps representing the FD voxel-based signature for the estimated Nakagami shape and scale parametric images are shown in Fig. 9(e) and (f), respectively. Various wavelet decomposition techniques apply the sub-bands’ energy for decomposition which is susceptible to intensity variations in ultrasound images due to speckle; however, the local density function known as the FD, can instead overcome these local variations in voxel intensities as it gives a representation of texture surface roughness, and hence is employed for the multiresolution analysis. The fractal characteristics are estimated for all sub-bands at each level of the wavelet packet decomposition, where the FD is computed on a voxel-by-voxel basis to produce a fractal map for each sub-band i.e., each voxel in the fractal map has its own localized FD value estimated from its neighborhood as described in the previous section, where the rougher the surface the higher the FD values get, and vice versa. Features at boundaries are computed after assuming that each slice is mirror-like continually extended in both directions. Specifically, the fractal features for a specific sub-band j to a certain level of decomposition i represent the average value of the generated fractal image map of a volume of interest k as defined in Eq. (12).

Fig. 9

Example of a voxel-based tissue characterization for a non-progressive liver tumor case. The tumor 3D volume is reconstructed in (a), and the B-mode middle slice (b) after transforming using the MNF algorithm is shown in (c–f). The Nakagami shape and scale parametric voxels (c) and (d) and the corresponding multi-resolution fractal slice maps (e) and (f) illustrates how the case responds to chemotherapy treatment – the blue color regions in (e) and (f) which correspond with the RECIST criteria. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

This local estimation gives a more reliable estimation compared to a single global value. Finally, the optimized multi-fractal feature vector descriptor consists of all selected sub-band fractal feature signatures f in each volume of interest k, expressed as , where .

In order to save processing time, the dimensionality of the extracted feature vector is reduced by applying a differential threshold which eliminates weak FD signatures. The threshold is defined by the condition such that the FD signature absolute difference of the previous decomposition level is satisfied, then the decomposition should terminate. The new 3D multi-fractal Nakagami feature descriptor is abbreviated subsequently as MNF, and its estimation summarized in the following pseudo code:

Multifractal Nakagami feature descriptor estimation

Experiments

This section describes experiments on pre-clinical and clinical images to illustrate the new MNF algorithm and to compare its characterization performance with previous single scale methods. A tumor was classified as non-progressive if categorized as partial response and progressive if no change or disease demonstrated non-responsiveness. The response evaluation criteria in solid tumors (RECIST) was adopted to categorize the cases into progressive versus non-progressive (Eisenhauer et al., 2009). The baseline cross-sectional imaging was compared against those performed at the end of treatment according to the RECIST criteria to determine response to treatment for each target tumor.

Data

Pre-clinical data: xenograft tumor imaging protocol

RF ultrasound data was acquired using a diagnostic ultrasound system (z.one, Zonare Medical Systems, Mountain View, CA, USA) with a 10 MHz linear transducer and 50 MHz sampling. The output 2D image size was  mm with a resolution of  pixels. A total of 227 cross-sectional images of hind-leg xenograft tumors from 29 mice (20 progressive or stable disease and 9 non-progressive disease) were obtained with 1 mm step-wise movement of the array mounted on a manual positioning device until the whole tumor volume was imaged (Fig. 4). All studies were ethically approved and performed in line with UK Home Office regulations, and in accordance with personal and project licenses.

The 2D images were composed together to create a 3D ultrasound volume. In order to ensure that nearby healthy tissue is not included in tumor tissue characterization, two expert radiologists manually segmented each image prior to applying tissue characterization. The Nakagami distribution was fitted to the distributions in each voxel lattice and parametric volumes were generated for each tumor.

The complete 3D RF ultrasound dataset along with a description of case categorization can be downloaded from the following weblink url: https://ibme-web.eng.ox.ac.uk/livertumour. An example of one of the cases presented as an animated GIF for the fractal slice maps and a video of the corresponding fractal volume map can be found with the dataset.

Clinical data: clinical study imaging protocol

Cross-sectional images of liver tumors undergoing chemotherapy treatment obtained as part of an ethically approved prospective study was used to validate our proposed technique. A total of 394 cross-sectional images (186 from tumors demonstrating partial response categorized as non-progressive, and 208 from tumors demonstrating progressive disease categorized as progressive) were obtained using a diagnostic ultrasound system (z.one, Zonare Medical Systems, Mountain View, CA, USA) with a 4 MHz curvilinear transducer and 11 MHz sampling. Each dataset was acquired prior to commencement of chemotherapy. Response to treatment was determined based on conventional computed tomography follow up imaging as part of the patient standard clinical care according to the RECIST criteria (Eisenhauer et al., 2009).

The transducer beam was initially directed through the target liver tumor in the intercostal imaging plane. Patients were asked to maintain breath hold inspiration, in order to stabilize the tumor target during image acquisition. Using a smooth movement of approximately constant speed, the ultrasound probe was angled whilst maintaining a skin contact position in a cranial to caudal direction to capture sequential 2D cross-sectional images of the target liver tumor. Each output 2D image size was  mm with a resolution of  pixels. Similar to the xenograft tumor dataset, the 2D images were composed together to create a 3D ultrasound volume for each target tumor. The acquisition was repeated in a similar fashion three times at each time point. Manual segmentation of the liver tumor was also performed in a similar fashion prior to image texture analysis.

Fractal maps

Fig. 9 shows the new parametric mapping on an example case. The various scattering conditions related to tissue characteristics within the tumor texture are shown using color mapping, where the various intensity distributions corresponding to the local concentration of scatterers vary from pre-Rician (0 < μ < 0.5), generalized Rician (μ = 0.5), pre-Rayleigh (0.5 < μ < 1), Rayleigh (μ = 1), post-Rayleigh (μ > 1) as illustrated in Fig. 9(c), and from low (0 < ω < 3), mid (3 ⩽ ω < 7), and high (ω > 7) for local backscattered energy as in Fig. 9(d). The MNF fractal slice maps – which correspond to each slice in the reconstructed volume – shown in Fig. 9(e) and (f) correspond to the Nakagami shape and scale voxels of the mid-slice, respectively. Regions with different texture properties (smoother or with lower local fractal dimension values) become more apparent compared to the Nakagami parametric voxels. Another interesting point is that the scale fractal map highlights the low activity regions near to the edge of the tumor, while the shape fractal map complements the characterization by highlighting low activity regions belonging to the inner part of the tumor tissue texture. A holistic view of the overall progression or regression of tumor spread can be effectively revealed via the fractal volume maps where low activity regions which correspond to necrotic tissue are labeled by a dark red color as illustrated in Fig. 10. Also a comparison between a non-progressive and progressive case and for pre and post-chemotherapy treatment is shown in Fig. 11.

Fig. 10

Fractal volume maps: Volumetric rendering of Nakagami parametric scale (first row) and shape (second row) for a progressive and non-progressive liver tumor volume, respectively.

Fig. 11

Pairwise horizontal comparison of texture-based volume rendering of Nakagami scale (a–d) and shape (e–h) multi-fractal volume maps. (a) and (b) Are an example of a non-progressive case in pre and post-chemotherapy treatment, and the (c) and (d) are for a progressive case in pre and post-chemotherapy treatment, respectively; (e) and (f), and (g) and (h) are the corresponding volumetric Nakagami shape cases. Red color labels indicate low local fractal dimension or low-activity regions which correspond to necrotic regions according to RECIST criteria. In first row, it is noticed that the spread of the red voxels has increased in post-treatment as compared to pre-treatment, and vice versa in the second row. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Statistical analysis

Since our primary concern is to demonstrate the texture expressiveness of the new multi-fractal feature descriptor in this subsection we describe experiments conducted to compare the discriminative power of the MNF descriptor with established features for mass classification. To do this we have chosen to use a naïve Bayesian classifier, although SVM or random forests might also have been used.

Specifically, we consider the MNF estimated over multiple scales in each case k in which there are features per voxel, where i is the number of decomposition levels estimated adaptively (see step 7 in Algorithm 1). This feature vector was fed into a simple naïve Bayesian classifier (nBC) determine performance of classifying progressive versus non-progressive cases.

Fig. 12 summarizes classification results for the pre-clinical dataset using the proposed MNF features with results from six other classic filter, model and statistical-based texture analysis methods and for B-mode intensity and Nakagami-based volumes of interest. The compared texture analysis methods are: Gabor filter (GF), fractional Brownian motion (fBm), Gaussian Markov random field (GMRF), gray-level co-occurrence matrix (GLCM), run-length matrix (RLM), and autocovariance function (ACF). Details of the extracted features can be found in Table 1, and a leave-one-out validation approach was employed. MNF-based performance gave the best overall cross-validation accuracy of 98.95%. Also the Wilcoxon Signed-Rank test on paired accuracy in Nakagami and intensity-based of each subject for both two-class classification shows that there is a significant difference ().

Fig. 12

Performance comparison for the MNF method against the filter-based Gabor filter (GF), model-based fractional Brownian motion (fBm) and Gaussian Markov random field (GMRF), and statistical-based gray-level co-occurrence matrix (GLCM), run-length matrix (RLM), and autocovariance function (ACF) texture analysis methods. The blue columns represent the operation of the texture descriptors on volumetric Nakagami parametric volume of interests, while the red columns are results from conventional ultrasonic intensity B-mode volume of interests. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1

Features extracted from the comparative texture analysis methods in Fig. 12.

Method	Texture features
GF	Energy of each magnitude response for five radial frequencies 2/2^6,2/2^5,2/2^4,2/2^3,2/2^2 with 4 orientations 0°, 45°, 90° and 135°
GMRF	Seven features estimated from a third order Markov neighborhood model
fBm	Mean, variance, lacunarity, skewness and kurtosis derived from the generated FD image
GLCM	Contrast, correlation, energy, entropy, homogeneity, dissimilarity, inverse difference momentum, maximum probability statistical features derived in the 0°, 45°, 90° and 135° directions
RLM	Short run emphasis, long run emphasis, gray level non-uniformity, run length non-uniformity and run percentage statistical features derived in the 0°, 45°, 90° and 135° directions
ACF	Peaks of the horizontal and vertical margins values and associated exponential fittings of the ACF

Clinical application

In order to demonstrate the applicability of the MNF algorithm for analysis of data acquired to clinical protocol the new method was applied to a clinical liver tumor dataset. In this case, RF ultrasound data is acquired in a fan-like scanning protocol (cf. the pre-clinical dataset was acquired using a linear transducer) in which a series of 2D images are collected as the transducer is tilted and then reconstructed into a 3D image.

Table 2 summarizes the classification performance for the clinical dataset following the same classifier design as in Section 4.3. The results show a good classification accuracy of 92.90% using a leave-one-out cross-validation approach, and a 92.01% ± 0.50 and 92.60 ± 0.30 for 5-fold and 10-fold cross-validation (results are the mean ± standard deviation of the performance over 60 runs). The texture descriptor was also compared against the RF backscatter signal using the localized voxel-based Nakagami parameters without generating the fractal volume maps (see Table 3), and by only deriving the fractal maps directly from the intensity (B-mode) images as well (see Table 4).

Table 2

Detailed classification performance for the 3D clinical RF ultrasound liver tumor test set using the MNF algorithm.

Classification performance	Cross-validation
Loo	5-fold	10-fold
Recall	0.935	0.92±0.919	0.921±0.931
FP rate	0.065	0.08±0.080	0.069±0.079
Accuracy	0.929	0.92±0.005	0.926±0.003
Precision	0.941	0.93±0.911	0.937±0.914
F-measure	0.928	0.92±0.005	0.925±0.003
J-Index	0.929	0.92±0.005	0.926±0.003
Dice SC	0.963	0.96±0.003	0.961±0.002
ROC Area	0.929	0.92±0.006	0.926±0.003

Table 3

Detailed classification performance for the 3D clinical RF ultrasound liver tumor test set using only the Nakagami parameters.

Classification performance	Cross-validation
Loo	5-fold	10-fold
Recall	0.715	0.73±0.492	0.72±0.491
FP rate	0.514	0.51±0.272	0.51±0.276
Accuracy	0.594	0.60±0.007	0.60±0.006
Precision	0.656	0.56±0.670	0.56±0.666
F-measure	0.595	0.60±0.007	0.60±0.006
J-Index	0.594	0.60±0.007	0.60±0.006
Dice SC	0.745	0.75±0.005	0.75±0.005
ROC Area	0.600	0.61±0.007	0.61±0.006

Table 4

Detailed classification performance using B-mode images of the 3D clinical ultrasound liver tumor test set.

Classification performance	Cross-validation
Loo	5-fold	10-fold
Recall	0.823	0.82±0.870	0.82±0.868
FP rate	0.135	0.13±0.182	0.13±0.180
Accuracy	0.845	0.85±0.005	0.85±0.004
Precision	0.845	0.85±0.842	0.85±0.844
F-measure	0.845	0.85±0.005	0.85±0.004
J-Index	0.845	0.85±0.005	0.85±0.004
Dice SC	0.916	0.92±0.003	0.92±0.002
ROC Area	0.844	0.84±0.005	0.84±0.004

Discussion

An outstanding challenge in medical ultrasound image analysis addressed in the paper is to provide an efficient characterization of subtle speckle textural changes or intra-heterogeneity within tumor tissue. The low signal to noise ratio and imaging artifacts frequently present in clinical ultrasound images make extracting such diagnostically useful information hard. In this section we focus on discussing three main contributions of our work and their importance.

Features of the method which make it novel and how this makes a difference: We have proposed a novel and meaningful fractal feature vector representation of ultrasonic signal characterization across spatial scales. This uses a multi-scale analysis where the voxel lattice is optimized to tumor size. We have also shown that performing the estimation in a volumetric fashion improves classification accuracy.

The proposed MNF approach simplifies analysis of the higher order statistics of the speckle texture via investigating different sub-bands at various decomposition levels using wavelets that tend to exhibit fractal characteristics. The tailored spatial-frequency localization provided via the Daubechies wavelet can facilitate the subsequent measurement of (signal) surface roughness while simultaneously filtering out irrelevant speckle features. The resulting fractal features quantify the speckle texture. Unlike conventional sub-band energy decomposition, sub-resolution level probing via the fBm (which satisfies an affine intensity invariance property) is not susceptible to sudden changes in the speckle intensity spatial frequency distribution.

Significance of the MNF descriptor to describe cancer morphology at the image level: The results in Section 4 are very interesting as they suggest that the MNF representation provides a new way to visualize cancer morphology using an ultrasound-based descriptor, and specifically to partition cancerous masses into necrotic and non-necrotic areas which has potential clinical utility.

In particular, the sparse blue colored clusters – which correspond to the low fractal dimension values shown in Fig. 9(e) and (f) – appear to correspond to necrotic regions in the tumor that are beginning to respond to chemotherapy treatment as defined by the RECIST criteria. This can be explained by first noting that tumorous tissue tends to have a “rougher” appearance than non-tumorous tissue due to the chaotic way that tumors build their network of new blood vessels. These angiogenesis networks tend to be leaky and disorganized, unlike blood vessel vasculature in normal tissue. This heterogeneity introduces a degree of randomness in appearance or “roughness” and gives a higher fractal dimension value compared to necrotic tissue.

Necrotic regions have different echogenicity characteristics (Sadeghi-Naini et al., 2013). In a necrotic region there is no cell growth and hence could be quantified if investigated at the appropriate analysis scale. The fractal volume maps reveal some of the intra-heterogeneity regions that tend to have a different textural characteristics to that of tumor tissue. These observations correspond to the dark red voxels in Fig. 10 and in Fig. 11 in the post-chemotherapy case. The red voxels in the pre-treated tumors in Fig. 11(a) and (c) represent regions with low activity, i.e. non-aggressive, with potential to become aggressive after treatment. We consider these red voxels as suspicious regions within the tumor and not as active as the rest of the malignant tumor tissue. From a pattern recognition perspective, the red voxels are closer in terms of their surface roughness characteristics to necrotic (i.e. low fractal dimension values) rather than aggressive regions (i.e. high fractal dimension values). However, at the pre-treatment stage, the nature of these regions is not yet confirmed since the tumor has not been subjected to chemotherapy treatment. The different shades of red in the pre-treatment figures reflect the varying degree of low activity that exists in the tumor. Subsequently, and after the first session of chemotherapy treatment, we notice that increase in these low activity regions in the non-progressive tumor of Fig. 11(b). Since the tumor has responded to treatment, we can now be confident that the aforementioned low activity regions in fact refers to necrotic regions. Conversely, the red voxels seen on the pre-treatment images for the tumor that progressed post chemotherapy nearly disappears on the post-treatment images of Fig. 11, thus suggesting that they have now become aggressive resulting in progression of the tumor.

We also observe that the scale and shape fractal maps tend to complement each other, in the sense of highlighting different aspects of the analyzed speckle texture pattern. This is evident when examining the Nakagami shape and scale parametric voxels of the non-progressive case shown in Fig. 9. Here we observe that a number of necrotic regions – highlighted in red in Fig. 9(c) and (d) – become more apparent in the corresponding fractal maps of Fig. 9(e) and (f). The scale fractal map highlights low intra-heterogeneity regions on the outer surface and near to the edge of the tumor. The shape fractal map complements this by revealing most of these low intra-heterogeneity regions belonging to the inner part of the tumor.

Expressiveness of MNF as a mass characterization feature: Results in Fig. 12 for the pre-clinical cases show that the MNF performed best based on Nakagami-based parametric images reporting an accuracy of 98.95%, while the fBm worked well on traditional B-mode intensity images with an accuracy of 83.16%. Similarly for the clinical cases, where under all classification performance metrics by comparing Table 2, Table 3, Table 4, we can see that the MNF algorithm outperformed both using the Nakagami parameters alone and when texture analysis was applied to the B-mode images as well. Moreover, we found that using a 2D single slice gave a lower accuracy of 74.74% compared to a volumetric analysis (98.95%). Furthermore, combining the texture-based multiresolution fractal features extracted from the Nakagami maps via the MNF algorithm with the features extracted from the conventional B-mode intensity images resulted in 4.1% degradation in the overall accuracy, as compared if the MNF method was employed alone. Note that in part due the limited data available in the pilot work, we have not looked at the benefits from a mass classification perspective of combining the MNF descriptor with other ultrasound tissue characterization parameters or image texture features which would be a natural topic to explore in the future.

Finally, we comment on three limitations of the current research which can be translated into opportunities for future investigation. Firstly, we are currently relying on the consensus of two radiologists to provide the gold standard and training data. Although we have shown good results with our current strategy, of possible concern is that the training samples may be mislabelled. This would reduce the accuracy of the results. Future work might look at the significance of this and strategies for mitigation. Secondly, fatty livers may result in attenuation of tissue properties and it would be interesting to investigate how this affects MNF classification accuracy. Thirdly, RF characteristics (and hence speckle appearance) tends to differ between ultrasound devices. It would be interesting to investigate whether a training set from one ultrasound scanner can be used for classification of images from a different scanner or results are scanner specific.

Conclusion

A new approach for assessing tumor heterogeneity via 3D multi-fractal multi-scale Nakagami-based feature modeling has been presented which we believe is the first work to consider intra-heterogeneity quantification of a cancerous mass. We estimated volumetric Nakagami shape and scale parameters from which the novel fractal descriptor is estimated. Future work will investigate the use of the method for both staging liver tumors and in longitudinal analysis as an image-based biomarker of tumor growth and therapeutic response.

Input: Volumetric ultrasound images I, volume of interests from each volume I: I_i=x_1,y_1,z_1,…,x_L_n,y_L_n,z_L_n, segmented volume of interests: V_ss=1Z.

Output: multifractal feature descriptor Λf(k)

1: for all Segmented volumes of interest V_1→V_Zdo

2:  {Step 1} //subdivide each volume of interestV_l̲into voxel latticesv^i̲.

3:  for all Voxel lattices v1111→vmnjzdo

4:  Fit with a Nakagami distribution N(x|μ,ω)=2μω^μ1Γμx^2μ-1

5:  {Step 2} //calculate Nakagami shapeμ̲and scaleω̲parameters using maximum likelihood estimation as:

θˆv=argmax_θDθ/v1111,…,vmnjz

where θ is a vector of parameters for the Nakagami distribution family fv1111,…,vmnjz/θ

6:  {Step 3}//construct Nakagami shapeN_μ̲and scaleN_ω̲parametric array.

7:  end for

8: end for

9: for all voxels v_i in N_μ and N_ωdo

10:  {Step 4}//Perform Daubechies wavelet packet transform

W_φt_0,x,y,z=1/mnj∑r=0j-1∑l=0n-1∑k=0m-1v_klrφ_klrt_0,x,y,z

Wψit,x,y,z=1/mnj∑r=0j-1∑l=0n-1∑k=0m-1v_klrψklrit,x,y,z

11:  {Step 5}//multifractal estimation and optimization phase

12:  for all decomposition levels ido

13:  for all voxels vmnri in W_φ and Wψido

14:   for all voxel pair distances Δr in W_φ and Wψido

15:   compute mean absolute difference Δv of each voxel pair q_i,p_i

16:   {Step 6}//construct a multidimensional volume of interest matrixN_d(x,y,d)̲

17:   normalize and take the logarithm Δvˆ=logΔvmnri/Δvmnri,

where m,n and r are the size of a voxel i at a certain scale

18:   normalize voxel pairs distances Δrˆ where Δr=∑i=1nq_i-p_i^2

19:   perform least square linear regression as

S_rr=∑i=1j-1Δrˆi2-∑i=1j-1Δrˆ_i^2/(j-1) ,

S_rv=∑i=1j-1∑k=1j-1Δrˆ_ivˆ_k-∑i=1j-1Δrˆ_i∑k=1j-1Δvˆ_k/(j-1)

20:   estimate the Hurst coefficient H matrix which represents the slope H=S_rv/S_rr

21:   estimate the fractal map I=3-H

22:   end for

23:  end for

24:  extract mean fractal dimension fi,jk←1/MN∑r=1MNI_r

where λ_k=f1,11,…,fi,jk,…,fm,nz

25:  construct feature descriptor from all wavelet sub-bands Λ=⋃_kλ_k

26:  fi+1,jk-fi+1,j+1k

27:  repeat

28:   {Step 7}//determine fractal absolute difference between decomposition level i and subsequent level as:

Dik=fi,jk-fi,j+1k,

Di+1k=fi+1,jk-fi+1,j+1k

29:  untilDi+1k⩽Dik

30:  end for

31: end for

32: return optimized multifractal feature vector: Λf(k)←argmaxλ_k