Multimodality quantitative ultrasound envelope statistics imaging based support vector machines for characterizing tissue scatterer distribution patterns: Methods and application in detecting microwave-induced thermal lesions

Highlights • A multi-QUS envelope statistics imaging SVM method for scatterer distribution.• Ultrasound log10(α), m, and hNSE as the input features of the SVM classifier.• Phantom simulations and MWA experiments of porcine liver ex vivo were conducted.• The proposed method outperformed single QUS envelope statistics imaging method.

Ultrasound envelope statistics imaging is an important group of QUS techniques, including ultrasound Nakagami imaging, homodyned-K (HK) imaging, and information entropy imaging, as they can characterize scatterer distribution patterns [3,[6][7][8][9][10][11][22][23][24][25][26], such as local scatterer concentrations and arrangements, which are useful in disease diagnosis, therapy evaluation or monitoring, and prognosis prediction.Ultrasound Nakagami and HK imaging methods are model-based envelope statistics imaging techniques, where the Nakagami or HK models are used for describing the probability density functions (PDFs) of the envelopes of ultrasound radiofrequency (RF) signals.The Nakagami [27,28] and HK [29][30][31] distributions are generalized statistical models for ultrasound backscatter envelope statistics [31].The Nakagami distribution has two parameters, m and Ω, which can describe the shape and scale of envelope PDFs, respectively.The shape parameter m is connected with scatterer concentrations and arrangements.The HK distribution has three parameters of physical meanings, α, ε, and σ, where α is the scatterer clustering parameter related to the effective number of scatterers per resolution cell, ε 2 denotes the coherent signal power, and 2ασ 2 represents the diffuse signal power [31].A derived parameter of the HK distribution, k, can describe the ratio of the coherent to diffuse signal, where k = ε/(σ ̅̅̅ α √ ).It should be noted that the Nakagami distribution can be considered as an approximation of the HK distribution [31].Ultrasound information entropy imaging methods are non-model-based envelope statistics imaging techniques, where entropy provides an insight into RF signal uncertainty and allows for the analysis of backscattered statistics without considering the distribution of data or the statistical properties of ultrasound signals [1].A histogram is usually used to estimate the probability value of each histogram bin, and these probability values are used to calculate the Shannon entropy.It has been demonstrated that ultrasound Shannon entropy imaging can reflect the changes of scatterer concentrations [32].
Thermal tumor ablation, such as microwave ablation (MWA) and radiofrequency ablation (RFA), is a minimally invasive therapy technique [33], where an applicator is inserted into a tumor to heat the tumor, causing coagulative necrosis of the tumor cells to kill the tumor in situ.High-intensity focused ultrasound (HIFU) does not need an applicator and enables heating of a target tissue volume from sources placed outside the body which produce a focused ultrasound field.The heating induced gas bubbles and tissue necrosis in a coagulation zone or ablation zone causes tissue microstructural changes which can be reflected in scatterer distribution patterns.These bubbles can play the role of "natural" contrast agents, contributing to additional acoustic scatterers in the ablation zone.The interaction of the bubbles with incident ultrasound waves results in additional signal components in the RF signals [34][35][36].Ultrasound envelope statistics imaging techniques have shown potential in detecting or monitoring heating-induced thermal lesions [11,[22][23][24], as they are capable to characterize local scatterer distribution alterations in the ablation zone, such as the formation and growth of heating-induced bubbles.Zhang et al. [37,38] explored the feasibility of using Nakagami imaging to evaluate the ablation zone induced by HIFU exposures in tissue-mimicking phantoms.Huang et al. [39] conducted ex vivo porcine liver HIFU experiments to compare Nakagami imaging to conventional B-mode ultrasound imaging in identifying thermal lesions.Wang et al. [40] demonstrated that Nakagami imaging could visualize thermal lesions of porcine liver ex vivo induced by RFA.Zhou et al. [41] proposed a real-time Nakagami imaging method based on frequency and temporal compounding for monitoring RFA and measuring the areas of ablation zones in porcine livers and muscle tissues ex vivo.Zhang et al. [42,43] demonstrated that Nakagami imaging had potential to monitor and characterize thermal lesions induced by MWA through ex vivo and in vivo porcine livers experiments.Our previous study [22] presented an algorithmic scheme of Nakagami imaging based on Gaussian pyramid decomposition, with MWA coagulation zone detection based on polynomial approximation (PAX).Our previous studies demonstrated the feasibility of HK imaging in monitoring thermal lesions in porcine livers ex vivo induced by MWA [11,24].Monfared et al. [44] demonstrated that Shannon entropy imaging could detect thermal lesions in porcine muscle tissue ex vivo induced by HIFU, and exhibited higher lesion-to-normal contrast comparing to the Nakagami imaging.Li et al. [45] proposed a horizontally normalized Shannon entropy (hNSE) imaging method for monitoring MWA-induced thermal lesions of porcine livers ex vivo and in vivo.
In this study, we focused on QUS envelope statistics imaging as it is one of well-established QUS techniques associated with tissue microstructures.Although ultrasound Nakagami imaging, HK imaging, and information entropy imaging methods have been used for characterizing tissue scatterer distribution patterns, an approach that integrates the strength of multimodality envelope statistics imaging techniques have not been investigated, particularly in terms of MWA.We hypothesized that using a machine learning model to combine the features of ultrasound Nakagami imaging, HK imaging, and information entropy imaging may enhance the capability to characterize tissue scatterer distribution patterns.In this study, we proposed a support vector machine (SVM) approach to combining multimodality QUS envelope statistics imaging techniques for characterizing tissue distribution patterns.We conducted simulation experiments and ex vivo porcine liver MWA experiments to evaluate the performance of the proposed method.The major contributions of this work are as follows.
• A multimodality QUS envelope statistics imaging based SVM method for tissue scatterer distribution characterization was proposed.• Ultrasound HK-log 10 (α ) , Nakagami-m, and hNSE parameters were extracted form parametric maps to construct the feature vectors, which were inputs of the SVM classifier.
• Phantom simulations and MWA experiments of porcine liver ex vivo were conducted to evaluate the performance of the proposed method.• The proposed method was compared to single QUS envelope statistics imaging method with adaptive threshold segmentation in identifying the boundary of coagulation zones.

Phantom simulations
The heterogeneous phantom simulations were conducted using the ultrasound simulation toolbox Field II [46,47].A 3.5-period Gaussian excitation pulse (central frequency = 7.5 MHz) was used as the incident wave and a 192-element linear-array transducer with 32 active elements was used to generate backscattered echo signals.Acoustic attenuation and noise were neglected.The width, elevation height, and kerf of the element was 0.3 mm, 4 mm, and 12.5 μm, respectively.The speed of sound was 1540 m/s and the sampling frequency was set at 30 MHz.The central frequency and sampling frequency were set to conform to the physical transducer in the ex vivo experiments.The three-dimensional (3D) resolution cell corresponding to the smallest resolvable detail was considered as the full width at half maximum of the point spread function in the axial, lateral, and elevation directions, which was estimated using the method in [48].A relatively rough approximation was used in this paper to regard the volume of resolution cell as a cube with the dimensions equal to each resolution [49].Thus, we obtained a 3D resolution cell of 0.3713 mm 3 (i.e., 0.7175 mm × 0.7194 mm × 0.7194 mm), and generated various phantoms according to different number densities of scatterers per resolution cell.The volume of the S. Li et al. heterogeneous phantoms was x × y × z (lateral × elevational × axial) = 40 mm × 1.5 mm × 40 mm, and there was a cylindrical dense scattering inclusion with a diameter of 10 mm in the center.A total of five types of heterogeneous phantoms were generated.When the scatterer concentration of the background media was 2 scatterers/resolution cell, the scatterer concentrations of the inclusion media were 16, 24, and 32 scatterers/resolution cell, respectively.When the scatterer concentration of the background media was 10 scatterers/resolution cell, the scatterer concentrations of the inclusion media were 32 and 64 scatterers/resolution cell, respectively.Thirty groups of experiments were conducted for each type of heterogeneous phantoms.Each group of experiments yielded a frame of backscattered signals, corresponding to a single phantom.The backscattered signals with a size of 1558 (sampling points) × 256 (scan lines) were obtained from each phantom.The inclusion with higher scatterer concentrations and the background with lower scatterer concentrations were used to simulate different scatterer concentrations.In this study, a SVM classifier was used to classify these two scattering media.

MWA experiments of porcine liver ex vivo
We conducted 60 cases of MWA experiments of porcine liver ex vivo (n = 60) at 80 W for 1 min (n = 20), 70 W for 2 min (n = 20), and 60 W for 3 min (n = 20).Porcine livers slaughtered and freshly excised within 12 h were purchased from a local market, and they were immersed in 0.9% saline before the experiments.Fig. 1 shows the experimental setup for MWA and ultrasound backscattered signal acquisition.The thermal lesions were induced in porcine livers ex vivo using a clinical MWA system (KY-2000, Nanjing Kangyou Co., Ltd., Jiangsu, China).Before cutting the porcine liver sample into an appropriate shape and placing it in the polycarbonate box (5 cm × 5 cm × 5 cm), the B-mode image was used to avoid large vessels or cavities in the area to be ablated [45].Fig. 2 shows the polycarbonate box for holding porcine liver samples ex vivo.Fig. 2(a = 0.7 mm] just happened to be stuck by the opening left by the upper cover and caused the transducer surface to sink to a certain distance from the inner box (in order to nearly contact the porcine liver samples).The insertion depth of the ablation needle was measured and set in advance, and the transducer scanned vertically from top to bottom on the plane each time where the central heating point of the needle tip was located.During the ablation, ultrasound backscattered data [1558 (sampling points) × 256 (scan lines)] were acquired by a C++ program [41] developed in-house on the ultrasound scanner (Model 3000, Terason, Burlington, MA, USA) at a sampling frequency of 30 MHz and saved as.bin format files at 2 frames/s for off-line processing.For each case of experiments, only the frame corresponding to the end time of MWA duration was used in this study.After the MWA, the transducer and upper cover were removed, and the ablation needle was pulled out.The inner box was then taken out from the outer box using the handle.There were slits with a width of 0.5 mm left on both sides of the inner box (on the scanning plane of the transducer) [Fig.2(c)].The porcine liver samples were sliced by a blade from top to bottom in the inner box (corresponding to the imaging plane through the center of the lesion) in preparation for gross pathology photograph.

Ultrasound envelope statistics parameter estimation and imaging
The two-dimensional envelope data were obtained from backscattered signals using the Hilbert transform.A sliding window with a specific window side length (WSL) was used to traverse the entire envelope data to obtain the parametric maps.For the local envelope data in each window, a specific estimator was used to calculate the parameter value, which was assigned to a new pixel located in the center of the window.A 50% overlap was set for the sliding window and cubic spline interpolation was adopted to make the parametric map have the same size as the envelope.Finally, the parametric images were obtained by using pseudo-color mapping.Specifically, each element in the resized parametric map was assigned a specific color through the MATLAB colormap() function, so that the color of a pixel in the yielded parametric image corresponded to a specific value of the envelope statistics parameter.A colormap was a matrix of values that define the colors for images and consisted of three columns.Each row in the matrix defines one color using an RGB (red, green, blue) triplet.An RGB triplet is a three-element row vector whose elements specify the intensities of the red, green, and blue components of the color.We used the MATLAB colormap(jet) function to gradually change the pixel color from blue to cyan, green, yellow, and finally to red, indicating the variation of parameter values from 0 to 1.In addition, the B-mode image can be yielded by logarithmic-compression and grayscale mapping of the envelope data.
For the selection of the estimators for envelope statistics parameters, we considered fundamental ones that have been frequently used in the field of QUS, because the focus of this study was to evaluate the feasibility of multimodality envelope statistics imaging based SVMs in enhancing scatterer distribution pattern characterization.Specifically, the XU estimator [50] was used for HK parameter estimation because of its low root mean square error, the moment-based estimator [28] was utilized for Nakagami parameter estimation because of its high computational efficiency, and the probability distribution histogram (PDH) based estimator [45] was adopted for hNSE parameter estimation because of its sensitivity to disorder of data.

The HK parameter estimation
Under the HK distribution, the PDF of envelope data is modeled as [31,50] where A⩾0 denotes the envelope amplitude; α > 0 denotes the scatterer clustering parameter and reflects the effective number of scatterers per resolution cell; ε 2 (ε⩾0) and 2ασ 2 (σ > 0) denote coherent and diffuse signal powers, respectively; μ = ε 2 +2ασ 2 is the mean intensity; and J 0 (.) denotes the Bessel function of the first kind, zeroth order.The HK-α parameter was demonstrated to be feasible as a parametric imaging and was better than the HK-k parameter (i.e., the ratio of the coherent to diffuse signal) for MWA thermal lesion monitoring [24].From the experiments, we found that the value of the HK-α parameter had a large dynamic range and parametric imaging of HK-log 10 (α ) had a better visualization capability than that of HK-α, so HK-log 10 (α ) was used as a QUS feature input to the SVM model in this study.The estimation method of HK distribution parameters based on the mean intensity and X-and U-statistics, namely the XU estimator, was used in this study [50]: (2) X = IlogI/I − logI; (3) where I = A 2 is the intensity of envelope signals.Given a sample set of envelope signals {A 1 , …, A N }, the notation g(A) represents the mean value of the function g(.) over the sample set.For instance, I is equal to In order to solve the nonlinear system of equations for the variables ε 2 , σ 2 , and α, it would be convenient to adopt an algorithmic parameter γ, γ = ε 2 / ( 2σ 2 ) .Consequently, the analytical expressions X HK (α, γ) and U HK (α, γ) for the two log-moments as a function of the model parameters were obtained, and the estimation was formulated as a minimization problem [50]: where α max denotes the upper bound of α (α max = 40.5 in this work), and the subscript 'HK' denotes analytic (i.e., theoretical) values.The parameter α which optimally satisfies the equations were estimated using a binary search algorithm [50].In this study, four times the transducer PL [11] was adopted as the WSL to obtained the HK-log 10 (α ) parametric map.

The Nakagami-m parameter estimation
The Nakagami distribution has an envelope PDF as [28]: where Γ(.) denotes the gamma function; m is the shape parameter, and Ω is the scale parameter.It has been shown that the variation of m parameter between 0 and 1 could describe the PDF of envelope amplitude ranging from pre-Rayleigh to Rayleigh distributions [28].This characteristic allows m parametric imaging to detect the variation of local scatterer concentrations in a tissue [51].It was demonstrated that a sliding window of 3 PLs is sufficient to obtain stable parameter estimation [51].In this study, the m parameter was estimated by the moment-based estimator [28]: where E[.] denotes the expectation operator.

The hNSE parameter estimation
The Shannon entropy of the backscattered envelope signal r can be estimated in a discrete form [52]: where w(r) denotes the probability distribution of envelope data, which can be calculated using the PDH.The bin number, bins, was set to 60 in accordance with the phantom simulations and ex vivo experiments.The PDH consists of bins that divide the bin limits [r min , r max ] into equal intervals, so r and w(r) are the center and height of bins, respectively.The height of a bin depends on the number of data falling into the bin divided by the total amount of envelope data.For the hNSE parameter estimation, ultrasound backscattered envelope data were normalized into [0, 1] through The bin limits of the hNSE estimator are always [0, 1], rather than varying with envelope data.Therefore, the entropy value was estimated by [45] A sliding window with the WSL of 1 PL [32] was adopted to construct the hNSE parametric images.

Feature extraction and classifier training
The QUS-based multi-parameter SVM classifier algorithm scheme is shown in Fig. 3(a).In the preprocessing step, the HK-log 10 (α ) , Nakagami-m, and hNSE parametric maps were obtained through the estimation methods described above.Then, each parametric map was normalized to [0, 1] by For the simulated phantom data, feature combinations were randomly selected from 20 of the 30 experiment groups for each of the three types of inclusion scatterer concentrations to construct the training set, so the total number of feature combinations was 4000.Note that each group corresponded to one frame of data.For the data collected from MWA experiments of porcine livers ex vivo, feature combinations were selected randomly from 45 of 60 cases (i.e., 15 cases were selected for each combination of treatment power and duration) to form the training set, so the total number of feature combinations was 9000.We used the LIBSVM toolbox [53] to train a SVM classifier.The SVM classifier was trained using Gaussian radial basis functions and a sequential minimal optimization (SMO) algorithm [54] to classify the two classes of data.The soft margin parameter C and kernel parameter γ should be set for model training, which were optimized using a grid search method [55].A 10-fold cross validation was performed to find the optimal hyperplane for classification.Briefly, the training set was randomly divided into 10 equally sized subsamples.The classifier was trained on 9 of the sample sets and tested on the remaining sample set.The process was then repeated 10 times.The accuracy of cross validation was computed and the optimal C and γ parameters were obtained when the accuracy was the highest.The optimal C and γ parameters were used to obtain the final classification model.

Classifier performance evaluation
The remaining 10 groups of five separate phantoms and 15 cases of MWA data were used to provide testing samples to evaluate the performance of the SVM classification model.Fig. 4 shows the process of testing sample construction and classification.Unlike the training set, the testing samples were constructed pixel by pixel using the estimated three QUS parametric maps, that is, the HK-log 10 (α ) , Nakagami-m, and hNSE parameters corresponding to the same pixel point were composed into a feature combination according to the original size of ultrasound backscattered data matrix.Then, the classifier predicted the classes of the feature combinations pixel by pixel into a 0/1 label, and a binarized image was obtained.The performance of the classifier models was evaluated by classification accuracy ACC, sensitivity SEN, specificity SPE, and the area under receiver operating characteristic (ROC) curve (AUC): where TP (true positive) and FP (false positive) denote prediction classes that were judged as thermal lesions by the classifier whereas in fact they were thermal lesions and normal tissues, respectively.FN (false negative) and TN (true negative) denote prediction classes that were judged as normal tissues by the classifier whereas in fact they were thermal lesions and normal tissues, respectively.The ROC curve analysis was performed using the MATLAB perfcurve() function, in which the AUC is equal to the probability that a classifier will rank a randomly chosen positive instance higher than a randomly chosen negative one [56].

Thermal lesion detection
After obtaining classification results using the QUS-based SVM classifier in Fig. 4(d) for testing samples, binarized segmentation images were yielded.In general, the thermal lesion region should be an independent connected region, because the temperature field caused by power radiation was a continuous field during MWA [57].Therefore, independent spots would be filtered out and holes of predicted ablated regions would be filled.To this end, the MATLAB imfill() function was used to fill holes in the classification results [i.e., the binarized image in Fig. 4(d)], and the MATLAB regionprops() function was used to compute the areas of each connected component in the image.Then, the maximum connected region in the final segmentation image was regarded as the thermal lesion region [Fig.4(e)], and the smaller connected components were all deleted.In addition, binarized segmentation images were also obtained from HK-log 10 (α ) , Nakagami-m, and hNSE parametric imaging by the Otsu algorithm [58] to compare with the proposed multi-parameter SVM classifier in terms of the performance of thermal lesion detection.Since Otsu [58] is a non-parametric and unsupervised method of automatic threshold selection for image segmentation, we took advantage of its simple and convenient calculation to achieve thermal lesion detection rapidly.In each binarized image, the target (ablated) region was denoted by 1, and the background (non-ablated) region was denoted by 0.

Thermal lesion detection performance evaluation
We used the ImageJ software (National Institutes of Health, USA) to outline the thermal lesion region from the photo of the tissue section (gross pathology) and binarized it into a mask as the classification ground truth.The images after binarized segmentation, spots filtering, and hole filling were resized to the same pixel size of the mask so that the area detection accuracy of thermal lesions could be calculated as [22] ACC area = where Area truth is the ground truth area of thermal lesions obtained from the mask, and Area detected is the predicted area calculated from the binarized image, which was obtained from multi-parameter SVM classifier or single parametric imaging.Besides, the Dice score DS [59] and Hausdorff distance HD [60] were used to evaluate the similarity between the finally obtained binarized images (i.e., the predicted region of thermal lesions) and the ground truth (the actual region of thermal lesions), which were calculated as HD where a and b are pixels of the predicted binarized image I pred and the ground truth image I truth , respectively, and dist(a, b) is the Euclidean distance.The larger the DS (with the maximum of 1) and the smaller the HD (with the minimum of 0), the higher the similarity between the prediction and the ground truth.The signal processing, parameter estimation, SVM training and testing, and the calculation of evaluation indicators in the five types of phantom simulations and MWA experiments were all implemented offline using MATLAB R2020b (MathWorks, Natick, MA, USA).

Statistical analysis
The statistical analysis of evaluation indicators was implemented using IBM SPSS Statistics 20 (IBM Corp., NY, USA).A Shapiro-Wilk normality test was performed on area accuracies, Dice scores, and Hausdorff distances to determine whether they followed a normal distribution [61].One-way analysis of variance (ANOVA) test [62] and Kruskal-Wallis test [63] were used to analyze the difference in means for normally distributed groups and difference in medians for non-normally distributed groups, respectively.Thus, we checked whether there were significant differences between the proposed method and single QUS envelope statistics imaging methods.

Results
Fig. 5 shows B-mode images of five kinds of heterogeneous phantoms and boxplots of HK-log 10 (α ) , Nakagami-m, and hNSE parametric features extracted from the circular inclusion and background regions of parametric images, respectively, which were normalized according to Eq. (11).Each features was expressed as the median and interquartile range (IQR).The median HK-log 10 (α ) was 0.74 (IQR: 0.66-0.81),0.75 (IQR: 0.67-0.83),0.78 (IQR: 0.71-0.86),0.83 (IQR: 0.78-0.90),and 0.86 (IQR: 0.78-0.93)for the circular inclusions in the five kinds of heterogeneous phantoms, respectively, and the median HK-log 10 (α , Nakagami-m, and hNSE features extracted from circular target regions were substantially higher than those extracted from background regions.However, there was a certain degree of overlap between the IQRs of HK-log 10 (α ) and Nakagami-m features extracted from circular target regions of 32 scatterers/resolution cell and those extracted from background regions of 10 scatterers/resolution cell.
Fig. 6 shows a case of heterogeneous phantom and its B-mode, HKlog 10 (α ) , Nakagami-m, and hNSE parametric images with corresponding Otsu segmentation and binarized images after spots filtering and hole filling, comparing with a SVM classifier with the features input of [log 10 (α ) , m, hNSE] to identify the circular inclusions.The binarized images segmented from log 10 (α ) , m, and hNSE parametric imaging exhibited much interference in the background, while the binarized image predicted by SVM classifier seemed to yield least interference.The final binarized images obtained from both the hNSE parametric imaging and SVM classifier exhibited similar contour to the ground truth.
Table 1 shows the performance of the SVM classifier on the testing set of simulated heterogeneous phantoms.For different types of heterogeneous phantoms, the prediction accuracy exceeded 90%, and high sensitivity, specificity, and AUC were achieved except the case of 32 scatterers/resolution cell in the inclusion and 10 scatterers/resolution cell in the background region.Table 2 shows mean area accuracies, Dice scores, and Hausdorff distances of circular inclusion detection from five kinds of heterogeneous phantoms using single HK-log 10 (α ) , Nakagamim, and hNSE parametric imaging with Otsu segmentation and the QUSbased multi-parameter SVM classifier.Expect for the case of 32 scatterers/resolution cell in the inclusion and 10 scatterers/resolution cell in  the background, the SVM classifier exhibited highest mean area accuracies (> 93%) in identifying the circular inclusions, and was significantly higher than Nakagami-m and hNSE parametric imaging; moreover, the SVM classifier exhibited smaller mean Hausdorff distances (< 13) and significantly higher mean Dice scores (> 0.93) than the HK-log 10 (α ) , Nakagami-m, and hNSE parametric imaging.For the case of 32 scatterers/resolution cell in the inclusion and 10 scatterers/ resolution cell in the background, the SVM classifier exhibited significantly higher mean area accuracy and Dice score, and smaller mean Hausdorff distance than the HK-log 10 (α ) and Nakagami-m parametric imaging.These indicated that the SVM classifier has a higher accuracy in both area detection and contour detection of target inclusions than single parametric imaging with Otsu segmentation.Fig. 7 shows the boxplots of HK-log 10 (α ) , Nakagami-m, and hNSE parametric features extracted from thermal lesions and normal tissues of parametric images, respectively.The median HK-log 10 (α ) was 0.65 (IQR: 0.56-0.75)for thermal lesion regions, and 0.52 (IQR: 0.43-0.63)for normal tissues.The median Nakagami-m was 0.34 (IQR: 0.20-0.50)for thermal lesion regions, and 0.09 (IQR: 0.04-0.20)for normal tissues.The median hNSE was 0.54 (IQR: 0.41-0.69)for thermal lesion regions, and 0.30 (IQR: 0.17-0.44)for normal tissues.On one hand, the IQRs of Nakagami-m and hNSE features extracted from the thermal lesion regions were substantially higher than those extracted from normal tissues.On the other hand, although the IQR of the HK-log 10 (α ) feature extracted from thermal lesion regions was basically higher than that extracted from normal tissues, there was a certain degree of overlap in the value range.) and m parametric imaging were significantly larger and smaller than the ground truth, respectively, whereas the ablated region in the binarized images segmented from hNSE imaging and predicted by SVM classifier was more similar to the ground truth.Meanwhile, the binarized image predicted by the SVM classifier exhibited less interference in the background than that obtained from hNSE imaging.
Table 3 shows the performance of the SVM classifier on the testing set of porcine liver ex vivo induced by MWA at 60-80 W, 1-3 min (n = 15).The mean classification accuracy, sensitivity, specificity, and AUC of the testing samples were 0.85, 0.84, 0.85, and 0.84, respectively.Table 4 shows the mean area accuracies, Dice scores, and Hausdorff

Table 2
The mean area accuracies, Dice scores, and Hausdorff distances of circular inclusion detection from five kinds of heterogeneous phantoms (10 groups of testing sets for each type of heterogeneous phantoms) using single homodyned-K log 10 (α ) , Nakagami-m, and horizontally normalized Shannon entropy (hNSE) parametric imaging with Otsu segmentation and the quantitative ultrasound based multi-parameter support vector machine (SVM) classifier.S/RC: scatterers per resolution cell.ADA: Area detection accuracy.DS: Dice score.HD: Hausdorff distances. 1Significant difference between log 10 (α ) and m. 2 Significant difference between log 10 (α ) and hNSE. 3Significant difference between log 10 (α ) and SVM. 4 Significant difference between m and hNSE. 5Significant difference between m and SVM. 6Significant difference between hNSE and SVM.distances of thermal lesions detection from porcine liver ex vivo induced by MWA (n = 15) using single HK-log 10 (α ) , Nakagami-m, and hNSE parametric imaging with Otsu segmentation and the QUS-based multiparameter SVM classifier.The SVM classifier achieved the highest mean area accuracy and Dice score, which were significantly higher than those obtained from HK-log 10 (α ) , Nakagami-m, and hNSE parametric imaging.Meanwhile, the mean Hausdorff distance obtained from the SVM classifier was significantly smaller than that obtained from HK-log 10 (α ) , Nakagami-m, and hNSE parametric imaging.

Discussion
In this study, a QUS-based SVM classifier was proposed for tissue characterization, specifically for characterization of tissue scatterer distribution patterns.The proposed method was compared with HKlog 10 (α ) , Nakagami-m, and hNSE parametric imaging through five types of phantom simulations and MWA of porcine livers ex vivo.Furthermore, the feasibility in detecting thermal lesions by the proposed method was evaluated.
Currently, there are several ultrasound imaging techniques for detecting thermal lesions or ablation zones, such as ultrasound elastography [64], contrast-enhanced ultrasound (CEUS) [65], QUS [23], and deep learning methods [66].Among them, heat-induced gas bubbles generated during the ablation can cause artifacts in ultrasound elastography of ablated regions [67], so ultrasound elastography is suggested for the evaluation of thermal lesions after ablation.Pohlman et al. [68] attempted to characterize post-ablation zones using electrode displacement elastography for 13 patients with hepatocellular carcinoma or liver metastasis and compared the binarized images of coagulation zones with clinical-standard-of-care segmentation, yielding an average Dice score of 0.88.CEUS provides highly sensitive visualization of blood flow in tissues because of the utility of intravascular microbubble contrast agents [69], which makes it the primary clinical method for evaluating the treatment effect after ablation.Liu et al. [70] used CEUS to quantitatively evaluate the thermal lesions of the livers of 10 young living swine induced by RFA; compared with histopathology,

Table 3
Performance of the support vector machine classifier on the testing set of porcine liver microwave ablation ex vivo (15 groups of testing sets).  1 Significant difference between log 10 (α ) and m. 2 Significant difference between log 10 (α ) and hNSE. 3Significant difference between log 10 (α ) and SVM. 4 Significant difference between m and hNSE. 5Significant difference between m and SVM. 6Significant difference between hNSE and SVM.
CEUS had an accuracy of 81.4%, a sensitivity of 83.3%, and a specificity of 76.9%.Methods of deep learning overcomes the difficulty of finding and computing the features related to lesions, but it needs a large amount of training data, which puts strict requirements on data acquisition and labeling.Zhang et al. [66] proposed an ultrasound imaging method based on a convolutional neural network (CNN) architecture for the detection and monitoring of thermal lesions induced by MWA in porcine livers, where RF data were acquired as the input layer of trained network, and segmentation images based on CNN were reconstructed to display the predicted thermal lesions, yielding an AUC of 0.89 and a Dice score of 0.87.In addition, single QUS envelope statistics imaging methods have been proposed to detect MWA-induced thermal lesions.In our previous work [11,22], the mean coagulation area detection accuracies of 20 cases of MWA-induced porcine liver ex vivo based on Nakagami-m and HK α and k parametric imaging combined with PAX were 89.26%, 89.47%, and 85.27%, respectively.The H-scan-based HK contrast-weighted summation α and k parametric imaging combined with Otsu segmentation, spots filtering, and hole filling achieved mean MWA coagulation area detection accuracies over 87% and 66%, respectively, with mean Dice scores being over 0.73 and 0.57, respectively [11].Li et al. [45] proposed hNSE imaging for MWA-induced thermal lesion detection, with an accuracy of 79.03%, a sensitivity of 79.19%, a specificity of 78.72%, and an AUC of 0.87.
The QUS method quantitatively extracts acoustic characteristic parameters from ultrasound backscattered echo signals.HK, Nakagami, and entropy parametric imaging is an important group of QUS imaging techniques, which uses the statistics in the backscattered echo signals to reflect changes of tissue microstructures.HK and Nakagami imaging methods are model-based techniques, which requires that the envelope of the backscattered signals obey a specific distribution model, whereas entropy imaging is a non-model-based method, which has better adaptability to the data.Note that prior to estimating parametric models, one could classify pixels into a few labels (corresponding to distinct echo envelope distribution) in order to meet the assumption of a single distribution on each sample set [71].Typically, a normal liver parenchyma can be considered as a collection of considerable scatterers [72].Most of the liver volume is occupied by hepatocytes, which contribute to diffuse scattering [73].The hepatocyte is very small; it is approximately cubical, with a side length of 20-30 μm and a volume of 3.4 × 10 − 6 mm 3 [74].Since the strength of ultrasound backscattering depends on the size of scatterers, the hepatocytes may be considered as weak scatterers exhibiting low echogenicities.During MWA, the rapid increase of temperature around the tip of ablation needle led to coagulative necrosis of the tissue in this region, where protein denaturation, cytoplasm disintegration, organelle dysfunction, nucleus and cell membranes destruction, and the presence of vapor bubbles all occurred [37,38].These might act as additional acoustic scatterers causing the increases of scatterer concentrations in the ablation zone.The HK α and Nakagami-m parameters that determine the statistical distribution of ultrasonic backscattered envelope were used to characterize the scatterer concentrations, since both α and m vary with the scatterer concentration.Moreover, ablation led to a higher uncertainty of backscattered signals in ablated tissues than in normal tissues.Entropy accurately captured and characterized the uncertainty and complexity of backscattered signals in ablated tissues which carried information of microstructure changes through calculation of PDHs.
Although the feasibility of deep learning in monitoring MWA has been demonstrated [66], the acquisition of large amounts of training data and labels becomes a major obstacle to its clinical application.In contrast, conventional machine learning algorithms, such as SVMs, can use less training data and labels to achieve the construction of classification models, and may be more explainable.The use of QUS parameters as input features of machine learning models has been studied in the lesion diagnosis and tumor classification [71,75,76].To the best of our knowledge, this paper is the first to use QUS parameters as the input features of the SVM classifier for MWA thermal lesion detection.In addition, our study found that the proposed method exhibited better performance in identifying coagulation zones than the single parametric imaging (HK-log 10 (α ) , Nakagami-m, or hNSE), which indicated that the combination of QUS-based multi-parameter features effectively enhanced the discrimination of ablated tissues.The SVM classifier can classify the image pixel by pixel after ablation, outperforming the method of Otsu image segmentation for single parametric images, which can effectively avoid the impact of the segmentation algorithm on the detection accuracy of the coagulation zones.
In the phantom experiments, we simulated five types of heterogeneous phantoms with higher concentrations in the circular inclusions and lower concentrations in the background media.When the concentration of the weak scattering background media was 2 scatterers/resolution cell and the inclusion concentrations were 16, 24, and 32 scatterers/resolution cell, the effectiveness of the proposed method in characterizing different scatterer distribution patterns was validated.We used scattering media of 10 scatterers/resolution cell (i.e., the echo amplitude follows the Rayleigh distribution) to represent normal liver tissues because the hepatocytes were explained to act as scattering sources at the Rayleigh scattering level [77].Since we did not know how many scatterers are added to the ablation zone due to MWA-induced bubbles, circular inclusions of 32 and 64 scatterers/resolution cell were used to represent the coagulation zone.It was found in Table 1 that the greater the difference of scatterer concentrations between the inclusion and background regions, the better the classification performance of the SVM classifier.
The optimal input feature vector, i.e., [log 10 (α ) , m, hNSE], was determined by experiments.For instance, we compared the accuracy, sensitivity, specificity, and AUC obtained using different feature combinations as inputs of the SVM classifier in the MWA experiments of porcine liver ex vivo (n = 15), and the results are shown in Table 5.Overall, the three-parameter feature combination [log 10 (α ) , m, hNSE] and four-parameter feature combination [log 10 (α ) , k, m, hNSE] both exhibited the best performance of classification.The results showed that there was no need to add the HK k parameter into the input feature of the SVM classifier.Moreover, it has been demonstrated that k parametric imaging was not as good as α parametric imaging in characterizing thermal lesions [11,24].This might be due to that the k parameter represents the ratio of the coherent to diffuse signal, whereas the characterization of ablation zones was more related to the scatterer concentration and arrangement, which could be better described by α,m, or hNSE.The feature vector [log 10 (α ) , m, hNSE] used in this study yielded better performance than the combinations of two features or one single feature.Besides, when the ratio of the scatterer concentration of inclusion-to-background was smaller, the performance of the proposed method would be lower (Tables 1 and 2).Since the input features of SVM, namely, m, log 10 (α ) , and hNSE, needed to distinguish different scattering media by characterizing scatterer concentrations and signal uncertainty, when the difference between the concentrations of the

Table 5
The accuracy, sensitivity, specificity, and AUC obtained by using different feature combinations as inputs of the SVM classifier in MWA experiments of porcine liver ex vivo (15 groups of testing sets).AUC: area under the operating characteristic curve; SVM: support vector machine; MWA: microwave ablation.background and inclusion media is small, the m, log 10 (α ) , and hNSE values extracted from these two regions may be nearly the same, which may not act as effective input features for SVM models.Therefore, the proposed method is suggested to characterize tissue scatterer distribution patterns for different scatterer concentrations.
This study has limitations.First, the simulated phantoms did not consider acoustic attenuation and noise, as the primary goal of this study was to validate the feasibility of multimodality QUS envelope statistics imaging based machine learning models (specifically, SVMs) in enhancing the capability of characterizing scatterer distribution patterns.Note that the echo envelope distribution depended on the settings of focused transducer [78] and the total attenuation of the intervening tissues, even when considering non-model-based methods, which also affected the parameter estimation.Besides, in order to keep the settings of the phantom simulations consistent with or close to those of the MWA experiments, the element width of the transducer we set in simulations was higher than half of the central wavelength, which might cause clutter and degrade the parameter estimation.The 100 pixels randomly selected from each simulated phantom acquisition are not guaranteed to be uncorrelated since some may be close to one another.Second, the types of simulated inclusions and the types of MWA experiments of porcine liver ex vivo are limited.In addition, only MWA experiments of porcine liver ex vivo was carried out to investigate the feasibility of the proposed method.In vivo animal experiments and clinical studies may be conducted in future work.Third, outliers in Nakagami parameter estimates in particular may be due to using the moments-based estimator rather than lower-variance methods like the maximum likelihood estimator.The better Nakagami parameter estimator can be used to replace the moments-based estimator in future work.In addition, only HKlog 10 (α ) , Nakagami-m, and hNSE were combined as the input features for the SVM classifier.Other parameter estimators may be used and more QUS parameters can be studied as the input of the machine learning model in future work.
The typical central frequency of the curved-array ultrasound transducer commonly used for clinical in vivo liver imaging is around 3.5 MHz.The central frequency of the linear-array ultrasound transducer used in this study (phantom simulations and ex vivo experiments) was 7.5 MHz.A higher central frequency corresponds to a smaller size of the resolution cell, which can make the ultrasound system more sensitive to changes in scatterer concentrations, and may provide a more optimistic performance of the SVM.However, a lower central frequency of the 3.5-MHz curved-array transducer corresponds to a larger penetration depth, which is inconvenient for MWA experiments of porcine liver ex vivo because the thickness of porcine liver samples ex vivo is limited.We used a linear-array transducer with a central frequency of 7.5 MHz to conduct ex vivo experiments, whose scanning depth was 4 cm.In order to be consistent with the configuration of the ex vivo experiments, the central frequency of the ultrasound transducer in the phantom simulations was also set to 7.5 MHz.Nevertheless, in accordance with the experimental results of this study, we hypothesize that for the 3.5-MHz curved-array transducer, the SVM model also can enhance the capability of scatterer distribution pattern characterization over the single ultrasound envelope statistics imaging method.This hypothesis needs to be validated in future work.

Conclusion
In this paper, a QUS-based multi-parameter SVM classifier was proposed for characterizing scatterer distribution patterns and for detecting MWA-induced thermal lesions.Through phantom simulations and MWA experiments of porcine livers ex vivo, the proposed method was demonstrated to outperform the HK-log 10 (α ) , Nakagami-m, hNSE parametric imaging in identifying the higher-scatterer-concentration zones and MWA-induced coagulation zones.
) shows that the polycarbonate box consists of three components: A. Outer box, B. Inner box (just enough to fit into the outer box without leaving any gaps), and C. Upper cover (with an opening of 18 mm × 50 mm, and just enough to cover the outer box without leaving any gaps).Fig. 2(b) shows the assembled object of three components.A water-cooled MWA needle (KY-2450B) was inserted into the sample horizontally through the circular hole with a diameter of 2 mm left on the side of the outer and inner boxes, and a linear-array transducer (12L5A, Terason) with a central frequency of 7.5 MHz [pulse length (PL)

Fig. 2 .
Fig. 2. The polycarbonate box for holding porcine liver samples ex vivo.(a) It consists of three components: A. Outer box, B. Inner box, and C. Upper cover.(b) The assembled object of three components.The ultrasonic transducer just happened to be stuck by the opening left by the upper cover.(c) There were slits with a width of 0.5 mm left on both sides of the inner box where a blade could be used to slice from top to bottom.

Fig. 3 .
Fig. 3. (a) The block diagram for the quantitative ultrasound-based multi-parameter support vector machine (SVM) classifier for training and predicting (testing).Ultrasound backscattered signals are preprocessed to obtain homodyned-K log 10 (α ) , Nakagami-m, and horizontally normalized Shannon entropy (hNSE) maps.Thus, each input backscattered signal matrix can provide feature combinations [log 10 (α ) , m, hNSE] ∈ R 3 as the inputs of the SVM classifier.With the training sets, SVM constructs the optimal hyperplane to separate two classes: thermal lesions and normal tissues.(b) Illustration of features extracted from each parametric map, that is, randomly selecting 100 pixel points from the target (thermal lesions, indicated by red crosses) and background (normal tissues, indicated by green points) regions, respectively.Each parametric map provides a feature vector containing 200 points, where the points from the background region are labeled as 0, and the points from the target region are labeled as 1.

)Fig. 4 .
Fig. 4. Illustration of the workflow of quantitative ultrasound-based multi-parameter classification and imaging.

Fig. 6 .
Fig. 6.A case of heterogeneous phantom and its (a) B-mode, (b) homodyned-K log 10 (α ) , (c) Nakagami-m, and (d) horizontally normalized Shannon entropy (hNSE) parametric images with corresponding Otsu segmentation and binarized images after spots filtering and hole filling; (e) a support vector machine (SVM) classifier with [log 10 (α ) , m, hNSE] as the input features was used to identify the circle inclusion.

Fig. 8
shows a case of MWA-induced thermal lesions in porcine livers ex vivo at 80 W, 1 min.It involves the B-mode [Fig.8(a)], HK-log 10 (α ) [Fig.8(b)], Nakagami-m [Fig.8(c)], and hNSE [Fig.8(d)] parametric imaging with corresponding Otsu segmentation and binarized images after spots filtering and hole filling.Fig. 8(e) shows a SVM classifier with [log 10 (α ) ,m, hNSE] as the input features to identify the thermal lesions, in which a 30 mm × 30 mm region of interest (ROI) selected from the tissue section image of the ablated porcine liver was binarized into a mask as the classification ground truth.The ablated regions in the binarized images segmented from HK-log 10 (α

Fig. 8 .
Fig. 8.A case of microwave ablation of porcine liver ex vivo at 80 W, 1 min: (a) B-mode, (b) homodyned-K (HK) log 10 (α ) , (c) Nakagami-m, and (d) horizontally normalized Shannon entropy (hNSE) parametric imaging with corresponding Otsu segmentation and binarized images after spots filtering and hole filling; (e) a support vector machine (SVM) classifier with [log 10 (α ) ,m, hNSE] as the input features was used to identify the thermal lesions, in which a 30 mm × 30 mm region of interest (ROI) selected from the tissue section image was binarized into a mask as classification ground truth.The square enclosed by dotted red line in (a) is the corresponding ROI.

Table 1
Performance of the support vector machine classifier on the testing set of five kinds of heterogeneous phantoms (10 groups of testing sets for each type of heterogeneous phantoms).S/RC: scatterers per resolution cell.

Table 4
The mean area accuracies, Dice scores, and Hausdorff distances of thermal lesion detection from porcine livers ex vivo induced by microwave ablation (15 groups of testing sets) using single homodyned-K log 10 (α ) , Nakagami-m, and horizontally normalized Shannon entropy (hNSE) parametric imaging with Otsu segmentation and the quantitative ultrasound based multi-parameter support vector machine (SVM) classifier.S/RC: scatterers per resolution cell.ADA: Area detection accuracy.DS: Dice score.HD: Hausdorff distances.