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

Graphical abstract


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 http://dx.doi.org/10.1016/j.media.2014.12.004 1361-8415/Ó 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 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 . 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 subbands 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.
This paper is organized as follows. State of the art and challenges associated with characterizing speckle tissue texture heterogeneity are summarized in Sections 2 and 3, followed by a detailed explanation of the proposed 3D multifractal Nakagamibased feature descriptor in Section 4. Sections 5 and 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 Fig. 1. 3D multifractal Nakagami feature descriptor algorithm design for ultrasonic tissue characterization introduced in this paper. 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 . In a subsequent work malignant tumors were shown to be more pre-Rayleigh distributed than those from benign counterparts ; 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) AND 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. 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 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. 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. 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 l parameter corresponding to the local concentration of scatterers, and the local backscattered energy represented by the scale parameter x i 0, for x > 0, and C Á ð Þ is the Gamma function. If x has a Nakagami distribution N with parameters l and x, then x 2 has a Gamma distribution C with shape l and scale (energy) parameter x/l.
The Nakagami distribution can model various backscattering conditions in medical ultrasound. By varying l, the envelope statistics range from pre-Rayleigh l < 1 ð Þ , Rayleigh 0 < l < 0:5 ð Þ , and to post-Rayleigh l > 1 ð Þ . The Nakagami parameters are generally estimated by the 2nd and 4th order moments, where given x is the ultrasonic backscattered envelope and E Á ð Þ 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 I i : i ¼ 1; . . . ; z f g subdivided into voxel lattices v i , each having a defined size of m and n, For each v i we assume that for a scaling factor r at a specific spatial scale s, the scaled voxel intensity lattice values v i klr of the RF envelope amplitude A klr such that A klr ¼ v klr ð Þ r2Rs , where the different possible resolution levels Rs : r ¼ 1; . . . ; s; . . . ; j reaching to the maximum level j represent a stochastic pattern, and the envelope amplitude of the scales r of v i klr 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 ( Having generated voxel-based Nakagami parameters, 3D wavelet packet  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.
Daubechies analysis (Mallat, 1999) can be applied at multiple scales. Namely: where v klr 2 L 2 R ð Þ is relative to scaling u klr and wavelet function w klr and W u t 0 ; m; n ð Þdefines an approximation of v klr at scale t 0 , and W i w t; m; n ð Þ coefficients add horizontal, vertical and diagonal details for scales t P t 0 . 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 h 0 k ð Þ and wavelet h 1 k ð Þ filters for image decomposition, and one does not need to express the u klr and w klr in their explicit form.
The decomposition process can be viewed as passing the signal through a pair of lowpass L ð Þ and highpass H ð Þ filters, also known as quadrature mirror filters, having impulse responsesh 0 k ð Þ and h 1 k ð Þ, 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 ash 0 a The decomposition is performed recursively on the output ofh 0 a ð Þ;h 1 a ð Þ andh 0 b ð Þ;h 1 b ð Þ. 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 m; n 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. 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).

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 nonstationary 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 n = 3 for 3-D space, is the voxel pair distances; H is called the Hurst coefficient; and the constant K i 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 N xyd defined for each of the tumor voxels v klr is derived at different range scales r, such that the mean absolute difference of each voxel pair Dv and for each voxel pair distances Dr 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 d ¼ 1; . . . ; j is the resolution limits of matrix N xyd which represents the mean absolute intensity difference to center voxels, and i stands for Nakagami shape l and scale x parametric images. Then each element from each array in N xyd is normalized after taking the logarithm and saved in a mean absolute difference row vector Dv. That is, the first element in all arrays of N xyd will compose vector Dv 1 , and all second elements will compose vector Dv 2 , and so on as shown in (8). This process is illustrated in Fig. 6.
The slope -which corresponds to the Hurst coefficient H -of the least square linear regression line of the log-log plot of Dv versus Dr 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 I:

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 A i > e A m , where A i is the tumor area in slice i ¼ 1; . . . ; m; . . . ; z, and e A m 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 x mle and the observed voxels values in each voxel lattice x v 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. 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 mm 3 to 6.60 mm 3 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.

Feature selection optimization
The ultrasound texture fractal maps representing the FD voxelbased 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 I 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 f i;j for a specific subband j to a certain level of decomposition i represent the average value of the generated M Â N fractal image map I of a volume of interest k as defined in Eq. (12).
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 K 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 8 f k perform least square linear regression as

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.

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 20 Â 54 mm with a resolution of 289 Â 648 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 65 Â 160 mm with a resolution of 225 Â 968 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. 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 < l < 0.5), generalized Rician (l = 0.5), pre-Rayleigh (0.5 < l < 1), Rayleigh (l = 1), post-Rayleigh (l > 1) as illustrated in Fig. 9(c), and from low (0 < x < 3), mid (3 6 x < 7), and high (x > 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 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.) 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.

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 K ðkÞ estimated over multiple scales in each case k in which there are K ðkÞ ¼ 8 Â i 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 (p < 0:05).

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 preclinical 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).

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 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. focus on discussing three main contributions of our work and their importance.
(1) 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.
(2) 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 . 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 regionshighlighted 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. (3) 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 Tables 2-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.