Statistical framework for ultrasonic spectral parameter imaging

Copyright

ELSEVIER

l

Ultrasound in Med. & Biol., Vol. 13. No. 9. pp. 117 l-1382. 1997 0 1997 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights rearned 0.101.56?9/97 $17.00 + .OO

PI1 SO301-5629(97)00200-7

Original Contribution

STATISTICAL FRAMEWORK FOR ULTRASONIC SPECTRAL PARAMETER IMAGING FREDERIC *Riverside

L. LIZZI,* MICHAEL ASTOR,* ERNEST J.

FELEPPA,*

MARY

SHAO-~ AND

ANDREW KALISZ* Research Institute, New York, NY 10036; +College of Physicians and Surgeons, Columbia New York, NY 10032 (Received

11 February

1997: in jinal form

30 Jul?;

University,

1997)

Abstract-This study examines the statistics of ultrasonic spectral parameter images that are being used to evaluate tissue microstructure in several organs. The parameters are derived from sliding-window spectrum analysis of radiofrequency echo signals. Calibrated spectra are expressed in dB and analyzed with linear regression procedures to compute spectral slope, intercept and midband fit, which is directly related to integrated backscatter. Local values of each parameter are quantitatively depicted in gray-scale cross-sectional images to determine tissue type, response to therapy and physical scatterer properties. In this report, we treat the statistics of each type of parameter image for statistically homogeneous scatterers. Probability density functions are derived for each parameter, and theoretical results are compared with corresponding histograms clinically measured in homogeneous tissue segments in the liver and prostate. Excellent agreement was found between theoretical density functions and data histograms for homogeneous tissue segments. Departures from theory are observed in heterogeneous tissue segments. The results demonstrate how the statistics of each spectral parameter and integrated backscatter are related to system and analysis parameters. These results are now being used to guide the design of system and analysis parameters, to improve assaysof tissue heterogeneity and to evaluate the precision of estimating features associated with effective scatterer sizes and concentrations. 0 1997 World Federation for Ultrasound in Medicine & Biology. Key Words: Ultrasonic spectrum analysis, Ultrasonic parameter images, Prostate ultrasonography, Liver ultrasonography, Integrated backscatter. INTRODUCTION

tral parameters and derived tissue features. We have collaborated with several medical institutions to evaluate the applicability of these parameter images to diagnosis. treatment planning and treatment monitoring in various organs, including the eye (Feleppa et al. 1986), prostate (Feleppa et al. 1996) and liver (Lizzi et al. 198X). The practical utility of spectral parameter images depends upon their statistical attributes. The accuracy and precision of local spectral estimates affect the reliability of differential diagnosis, and their precision is critical in detecting subtle morphologic alterations induced by either disease progression or successful therapy (Lizzi et al. 1997a). Evaluations of tissue heterogeneity and lesion detectability depend upon the probability density functions (pdf) of spectral estimates. In addition. estimator statistics are key considerations for selecting appropriate ultrasonic systems, implementing optimal analysis procedures and determining reasonable tradeoffs between image resolution and estimator precision.

Several laboratories are investigating alternative signalprocessing techniques to complement the information afforded by conventional medical ultrasonography. Many of these techniques employ spectrum analysis of radiofrequency (RF) echo data to derive such features as attenuation, integrated backscatter and sets of spectral parameters (e.g., spectral slope) (Insana et al. 1990; Thomas et al. 1989; Zagzebski et al. 1993). Spectrum analysis is also employed to derive estimates of tissue properties relating to the effective sizes, concentrations, and relative acoustic impedances of constituent scattering centers in tissue (Lizzi et al. 1987). Our laboratories have developed cross-sectional images that depict spec-

Address correspondence to: Dr. Frederic L. Lizzi, Riverside Research Institute, 330 West 42nd Street, New York, NY 10036 USA. E-mail:

[email protected]

1371

Ultrasound

1377

in Medicine

and Biology

Previous studies have analyzed relevant statistics of tissue spectra with emphasis on attenuation estimation (Kuc and Taylor 1982; Narayana and Ophir 1983). Recently. Chaturvedi and lnsana ( 1996) examined statistical error bounds on ultrasonic scatterer size estimates, and Huisman and Thijssen (1996) derived the standard deviations of several spectral parameters. The latter study, as one of our early reports (Lizzi et al. 1976), assumes conditions appropriate for region-of-interest spectrum analysis, in which spectra are averaged from many adjacent scan lines prior to parameter evaluation. As described below, this averaging process alters the statistics of derived parameters, and results are not directly applicable to high-resolution parameter imaging, in which parameters are estimated prior to optional averaging (i.e., image smoothing). The current study was undertaken to elucidate the statistics of several spectral parameter images that are finding clinical application. This report expands upon our previous analyses (Lizzi et al. 1997b) to include additional spectral parameters (linear-regression slope and intercept) that are being used to assess tissue type and to estimate physical scatterer properties. PDFs for each of these parameters are derived in terms of system characteristics and analysis parameters for the case of statistically homogeneous tissue structures. Theoretical pdf results are compared to relevant histograms derived from clinical parameter images of the prostate and liver; these organs offer distinctly different acoustic and morphological environments for validating theoretical results. The report discusses how these statistical results can be applied to integrated backscatter, a parameter developed by Miller and coworkers (Thomas et al. 1989). It also discusses effects of tissue heterogeneity on spectral parameter statistics.

METHODS

AND SPECTRUM PROCEDURES

ANALYSIS

The clinical data used in these studies were obtained from two ultrasonic instruments interfaced with digital data-acquisition systems developed in our laboratories. Liver data (Lizzi et al. 1988) were acquired using an ATL (Bothell, WA) mechanically scanned system (3MHz center frequency); prostate data (Feleppa et al. 1996) were acquired with a B&K Medical Systems (North Billerica, MA) trans-rectal scanner (5.75-MHz center frequency). In each type of examination, ultrasonic surveys were performed in the conventional manner. Then, RF data (8 bits) were acquired in selected scan planes using a computer-controlled analog-to-digital converter (LeCroy, Chestnut Ridge, NY) as described in the preceding references. Data from individual complete

Volume

23. Number

9, 1997

scans were acquired at real-time scan rates before the application of time-gain compensation (TGC) or logarithmic compression. In liver examinations. a 25-MHz sampling rate was used to acquire 2.5 kbytea of RF data along each of 128 scan lines; histologic evaluations of scanned lesions were performed following subsequent biopsies. In prostate examinations. a 50-MHz sampling rate was used to acquire 2.5 kbyte of RF data along 3 18 scan lines; ultrasonically guided needle biopsies were performed in selected planes immediately following acquisition, so that histological observations were available for known sites in scanned tissue planes. Acquired RF data were processed with a sequence of operations to generate cross-sectional images depicting specific spectral parameters (Lizzi et al. 1983; Feleppa et al. 1986). Along each scan line, a sliding Hamming window (typically 64 points) was used in conjunction with a fast Fourier transform (FFT) algorithm to derive the power spectra of RF-signals at sequential, partially overlapping sites (typical center-tocenter spacing of 8 samples). Each spectrum was divided by a stored calibration spectrum obtained from a waterglass planar interface placed in the transducer’s focal plane. This calibration removes extraneous system-transfer functions and references all spectra to a glass-plate reflector. The magnitudes of these unaveraged, calibrated power spectra were then converted to dB units for further processing over the useful frequency range, as determined from calibration, tissue and noise spectra (measured in an anechoic segment of a water tank). For liver spectra, this range was 2 to 4 MHz (~-MHZ bandwidth); for prostate spectra, it was 3.5 to 8 MHz (4.5MHz bandwidth). Calibrated spectra (in dB) at each tissue site were then analyzed using linear regression techniques to compute the spectral slope (dB/MHz), intercept (dB. extrapolation to zero frequency), and the midband fit (dB, value of the regression line at the center frequency). Although only 2 of these parameters are statistically independent, it is useful to consider all 3 because of their relationships to various tissue properties (Lizzi et al. 1987; Lizzi et al. 1988). Spectral slope is affected by intervening attenuation and is related to the effective sizes (correlation dimensions) of tissue scatterers. Spectral intercept is unaffected by intervening attenuation that is linearly dependent (in dB) on frequency; it is related to the effective sizes, concentrations and relative acoustic impedances of tissue scatterers. The midband fit is affected by attenuation and is related to the above tissue properties; as noted in the discussion, this parameter is closely related to integrated backscatter. In this study, we used the results of local spectral computations to generate cross-sectional images with

Parameter

image statistics

0 F. L. LIZZI ET AL.

1373

gray-scale levels that quantitatively depicted spectral slope, intercept and midband fit, respectively. Histograms of parameter values from selected image regions were compared with theoretical descriptions of the pdfs of each spectral parameter. ANALYTICAL

RESULTS

The previous section described our methods for computing spectral slope, intercept and midband fit. Because most examined tissues exhibit stochastic microstructures, these computed parameters form statistical estimates of their true ensemble values (Lizzi et al. 1983). This section derives the statistics of each estimator when the gated RF signals are received from many closely spaced, independent small scatterers. In this situation. the power spectrum technique decomposes signal components into N independent spectral cells, each spanning a frequency range Af determined by the duration of the analysis gate. We first analyzed the statistics of spectral components in each cell; we then treated the statistics of linear-fit parameters that describe the frequency-dependence of the set of cells in each spectrum. Statistics of single spectral-resolution cell The amplitude of the RF signals received from many closely spaced, independent subresolution scatterers exhibits a Gaussian probability density function (Wagner et al. 1986; Tuthill et al. 1988). Within a single spectral cell II, the Fourier transform of such signal has independent real and imaginary parts, each with a Gaussian pdf. The amplitude IS,,1of each spectral component has a Rayleigh pdf with standard deviation a, (Davenport and Root 1958). The power spectrum I&I is a sum of squares of two Gaussian random variables and, thus, has a Chi-squared distribution with two degrees of freedom (Lindgren 1962, Chapter 8). This distribution, often termed a negative exponential function, is:

(1) where y,, = I&l2 and the subscript designates the nth spectral cell. The mean and standard deviation of this distribution are both numerically equal to p,, = a:. When the power spectrum is expressed in dBs as 2, = 10 log lSn12, the corresponding pdf g(Z,) can be derived by setting g(Z,) dZ, = h(y,) dy,. The result can be expressed as g(Z,) = (1/4.34)exp[z/4.34

- In p,, - exp(z/4.34

- InpJ],

(2)

-30

-20

-10

0

SPECTRAL AMPLITUDE -MEAN

10

(dB)

Fig. 1. Probability density function of individual spectral cell.

where z denotes the values that can be assumed by Z,,. This distribution, which also pertains to logarithmically compressed B-mode images (Kaplan and Ma 1994), is plotted in Fig. 1. Its mean z and standard deviation (SD) (T=, can be evaluated using tabulated functions as described in the Appendix: Z,, = 10 log&)

- 2.5(dB);

~z,, = 5.6(dB) (3)

Thus, the SD of this pdf is constant and its shape is independent of its mean. (Note that the abscissas in all pdf plots represent parameter values minus their mean values.) Equation (2) describes the pdf of the power spectrum (in dB) in a single spectral cell. In the next section. we will use this result to compute pdfs of linear regression fits to N such independent spectral cells. At this point, we note that, because parameter images are of interest, we have derived eqns (2) and (3) for individual (unaveraged) power spectra expressed in dB. Different statistics apply for spectral estimates computed by region-of-interest (ROI) procedures that average spectra from M independent scan lines prior to conversion to dB. In the ROI case, Lizzi et al. (1976) showed that a,,,is equal to 4.34/V% if M is large and power spectra are averaged prior to dB conversion. If spectral magnitudes are averaged before squaring and dB conversion, results of Huisman and Thijssen (1996) and Parker (1988) showed that a-, is equal to 4.54/q/M for large M. Statistics of linear fit parameters As noted previously, the power spectrum of the RF signals being considered can be partitioned into N independent spectral cells, each of which has a pdf equal to

I ill

Ultrasound

in Medicine

and Biology

g(Z,,) ]eqn (2)]. N is equal to the signal bandwidth B divided by AjJ’;the bandwidth of a spectral cell. For a Hamming window of gate length L and a typical propagation velocity (1.5 mm / ps), 4f = l/L (Lizzi et al. 1997b); thus. N = BL where B is specified in MHz and L in mm. (Typically, the value of N ranges from 5 to 10 in clinical examinations.) The center frequency of each cell, ,f,,. is given by the relationship:

Volume

73. Number

9, 1997

as x, and let x = s,, - x,, where X, denotes the expected value of x in spectral cell n; then g(x,,) = g(.r + X;,). Now we restate eqns (7) and (8) in terms of G(M~).the characteristic function (Fourier transform) of g(x). Using eqn (3) and standard scaling relationships, these eqns become:

wr/N)]Nej(fi-fi’)‘“&, P(P) = r[G( (9)

.

f =A + (tz - Nl2)Af

(4)

!I

where ,f,. is the center frequency of the entire power spectrum. Spectral parameters are computed using linear regression techniques that determine the function m(fil j;.) + p that fits 2;; with minimum squared error (Lindgren 1962. Chapter 11). /3 is the estimator of midband fit:

P=;czn !,

(5)

(10)

where j denotes 2/-1

and:

p’ = kNi’ 10 log(&) n=O

- 2.5

(11)

N-I

and ttz is the estimator of spectral slope:

tTz =

m’ = ,zO$[10

(W

$z,, n Af

where (11 - N/2) ‘u = Cln _ N,2)z = 12(n - N/2)l(N3

+ 2~).

log(k.,) - 2.51.

(12)

p’ and m’ are the expected values of parameters @ and m as can be seen by computing the mean values of eqns (5) and (6a) and substituting z from eqn (3). It can also be shown that p’ and m’ represent the actual midband fit and slope values if z does indeed follow a linear trend with frequency. The final probability density function to be computed applies to spectral intercept, which is calculated as

(6b) The probability density functions associated with m and /3 can be derived from the pdf in eqn (2). Because Z,z are independent random variables, the pdf of /3, p(p), and the pdf of m, p(m), are (Papoulis 1965): P(P) = gLWN)*.

. .*g(Z,,-,lN)

p(m)= g(z”gj*.. .*g/zN-,$)

(7)

(8)

where * denotes convolution. Equation (7) is independent of AL and eqn (8) can be resealed for different values of Af by letting ZL = Z,/Af. Thus, these pdfs need not be reevaluated for different values of Aj. We denote the argument of g( ) in eqns (7) and (8)

I = p - rnJ.

(13)

This pdf is (Papoulis 1965): ~(1) = p(P)*(WJp(m)

.

(14)

Theoretical results for spectral parameter pdfs Results described in the preceding section were employed to calculate the pdfs for midband fit, slope and intercept for a range of N. Results for midband fit [eqn (9)] are shown in Fig. 2. For N = 1, the pdf is identical to the plot of Fig. 1. As N increases, the pdf becomes more symmetric about its mean value, approaching a Gaussian shape for large N. Figure 3 compares pdfs for N = 2 and N = 10 with Gaussian functions with the same standard deviations. (These functions correspond to the maximum likelihood, ML, Gaussian fits.)

Parameter image statistics 0 F. L. Lrzzr ETAL.

1375

r

Fig. 4. Probability density functions for spectral slope for various numbers N of spectral cells; Af = 1 MHz. Fig. 2. Probability density functions for spectral midband fit for various numbers N of spectral cells.

Figure 4 shows results for spectral slope [eqn (lo)] for various values of N for Af = 1 MHz. (These results p(m) can be resealed for other values of Afby computing

0.06

P D F A M P L I T U D E

Af p(mlAj)). These pdfs are symmetrical, as expected because statistical deviations from the linear fit are equally likely on either side off,. As indicated in Fig. 5, pdfs for low values of N are more sharply peaked than their ML Gaussian functions, but this dissimilarity decreases with increasing N.

P D F

N=2 0.10

0.05

/ -5 MIDBAND

5

0 FIT AMPLITUDE

;

0.02

;‘.-

I

I

0.04

D E

\ I

A M P L I

10

- MEAN

0

15

-30

(dB)

-15

15

0

SPECTRAL SLOPE - MEAN (dB / MHz)

P D F A M P L I T U D E

I

0.3

I

I

1 .o

I

I

I

P D F

N-10

N-10 0.2

A M P L I T U D E

0.1

0 -8

-4 MIDBAND

t 0 FIT AMPLITUDE

0.5

0 4 - MEAN

8

(da)

Fig. 3. Probability density functions (solid line) and maximum likelihood Gaussian fits (dashed line) for midband fit for 2 (top) and 10 (bottom) spectral cells.

I

-1

0

1

2

SPECTRAL SLOPE - MEAN (dB /MHz)

Fig. 5. Probability density functions (solid line) and maximum likelihood Gaussian fits (dashed line) for 2 (top) and 10 (bottom) spectral cells; Af = 1 MHz.

Ultrasound in Medicine and Biology

I?76

I

P D F A M P L I T U D E

I

I

1

I

I

I

!I,

0.15

MIDSAND FIT

0.10

0.05 SLOPE x CENTER FREQ.

0 -25

0

25

PARAMETER VALUES-MEANS

INTERCEPT-MEAN

50

(dB)

(dB)

Fig. 6. Probability density function for spectral intercept. Top plot shows constituent functions that are convolved to obtain the bottom plot of the intercept pdf; N = 5, Af = 0.5 MHz,f,. = 3 MHz.

Figure 6 shows a pdf result for spectral intercept (N = 5, Af = 0.5 MHz, f, = 3 MHz). For comparison, this figure also (top) plots the constituent pdf functions of midband fit and slope multiplied by center frequency; these curves were convolved to compute p(l) as indicated in eqn (14). These results show that the slope pdf component predominates in determining the intercept pdf because, under most operating conditions, the slope pdf is a much broader function than the midband fit pdf. As a check on these pdf results, standard deviations were numerically computed from the above plots and compared with results (Table 1) derived from classic analysis of linear regression procedures (Lizzi et al. 1997b). In all cases, excellent agreement was found. COMPARISON OF THEORETICAL CLINICAL RESULTS

AND

Theoretical pdf results were compared with histograms of spectral parameter data obtained from 3 clinical cases whose midband fit images are shown in Fig. 7a, b, c. The cases are (a) a healthy liver, (b) a liver with focal metastatic carcinoma, and (c) a prostate with carcinoma proximal to the bowel. (Both carcinomas were defini-

Volume 23. Number 9. 1997

tively diagnosed by subsequent biopsy procedures.) The liver images depict a total range segment of 6 cm; the trans-rectal, transverse prostate scan covers a 3-cm depth. The healthy liver image exhibits a progressive range-dependent decrease due to attenuation. The carcinoma images were compensated for an assumed attenuation coefficient as described below. In each clinical case, parameter images of spectral slope, intercept and midband fit were synthesized as described in Methods. A series of sites were analyzed using digital image analysis software to compute histograms of demarcated region-of-interest (ROI) segments. The ROIs were placed in visually homogeneous regions selected to avoid large reflectors. such as tissue boundaries or blood vessels; the same ROIs were used to analyze the 3-parameter images in each case. For all liver data, N was equal to 4, Af was 1 MHz, andf, was 3 MHz; for prostate data, N was equal to 5, Af was 1 MHz, and ji. was 5.75 MHz. For the healthy liver (Fig. 7a), a set of 5 widely distributed ROIs, comprising a total of 600 pixels, was employed to obtain representative tissue data. Gray-scale histograms from each parameter image were then converted to appropriate units (dB or dB/MHz). The range segment of each ROI was set at 2 mm to minimize attenuation effects within each ROI. The mean parameter value within each ROI was subtracted from its histogram; this operation removed effects of intervening attenuation at different depths and allowed parameter statistics to be evaluated in terms of stochastic perturbations from expected values. The resultant histograms for each region were then averaged for comparison with theoretical pdf results. Figure 8 shows the average midband fit histogram, which is seen to be in close agreement with the appropriate theoretical pdf (N = 4). Figure 8 also shows histograms for spectral slope and intercept measured from the same ROIs and processed in the same manner used with midband tit histograms. Again, these results are in close agreement with relevant theoretical results. Spectral parameter data from the liver metastatic

Table 1. Standard deviations Parameter Midband fit

Spectral slope Spectral intercept

of spectral parameters. Standard

deviation (SD)

Parameterimagestatistics0 F. L. LIZZIETAL.

(4

1371

(b)

Fig. 7. Clinical midband fit images: (a) healthy liver; (b) liver with hypoechoic metastatic carcinoma; (c) prostate with carcinoma (CA) proximal to bowel.

carcinoma (Fig. 7b) were first compensated for attenuation effects before image formation to relieve restrictions on ROI size and location. This compensation was implemented by multiplying midband fit values by 2aLRf where R represents range from the transducer and the attenuation coefficient a was set equal to 0.5 dB/ MHz-cm (Lizzi et al. 1988; Lizzi et al. 1997a). Slope values were compensated by multiplication by 2adz. No compensation was applied to intercept values, which are

theoretically unaffected by attenuation (in dB) that increases linearly with frequency. Figure 9 shows histogram results measured from a rectangular ROI that encompassed virtually all of the metastatic tumor seen in Fig. 7b. All measured histograms are in excellent agreement with corresponding theoretical pdfs. Spectral parameter data from the prostate (Fig. 7c) were also compensated for an attenuation coefficient of

137x

;

Ultrasound

in Medicine

and Biology

Volume

9. 1997

0.15

F

P II F

0.10

0.10

0.05

0.05

0

0 -15

P D F

23, Number

-10

-5 MEA; MIDSAND FIT (de)

5

10

-15

0.10 P D F 0.05

-10

I

I

-15

-10

0.10 L

-5 0 MIDBAND FIT - MEAN (de) I

nI

5

I

IO

I

II

i

I

0.05

0 -15

-10

-5 0 5 SPECTRAL SLOPE - MEAN (dB / MHz)

-20 SPECTRAL INTERCiPT-

IO

15

0 5 -5 0 SPECTRAL SLOPE -MEAN (d6 / MHz)

IO

I3

MEAN (de: SPECTRAL INTERCEPT-MEAN

Fig. 8. Measured histograms (bars) and theoretical pdfs (solid line) for healthy liver of Fig. 7a. Top = midband fit; middle = spectral slope; bottom = spectral intercept.

0.5 dB/MHz-cm prior to parameter image synthesis. Figure 10 shows results from a single ROI, encompassing most of the carcinoma (CA), which appears as a dark region just beyond the rectal wall in Fig. 7c. The midband-fit histogram agrees well with the theoretical pdf and exhibits the asymmetry expected for homogeneous scattering. Poorer agreement with theory is seen for slope and intercept histograms. Figure 11 shows histogram results from a large ROI that included most of the prostate posterior to the carcinoma of Fig. 7c. Inspection of the midband-fit image shows that this is a heterogeneous tissue segment, and the midband-fit histogram in Fig. 11 departs significantly from the theoretical pdf derived for homogeneous scatterers. The measured histogram is more symmetric about its mean, and it conforms very well to the indicated ML Gaussian fit for the data. The slope and, especially, the intercept results also are seen to match the ML Gaussian fits more closely than the theoretical pdfs. Table 2 lists the theoretical and measured SDS for each of the cases shown in Figs. 8-11. In the first 3 cases, ROIs were placed in apparently homogeneous tissue segments, and the theoretical SDS agree well with

(dS)

Fig. 9. Measured histograms (bars) and theoretical pdfs (solid line) for liver metastatic carcinoma of Fig. 7b. Top = midband fit; middle = spectral slope: bottom = spectral intercept.

measured values. The average absolute differences for these cases are 0.3 dB for midband fit, 0.43 dB/MHz for spectral slope, and 2.2 dB for spectral intercept. The fourth case involved visible heterogeneities and the average absolute difference between measured and theoretical SDS is 2.2 dB for midband fit, 0.26 dBNHz for spectral slope and 1.9 dB for spectral intercept. DISCUSSION We first discuss general attributes of our results, including applications to nonlinear spectral shapes, attenuation effects and estimator bias. We then discuss how our results pertain to estimator precision for spectral parameters and integrated backscatter and how they provide a basis for future assays of tissue heterogeneity. The statistical analysis described in this report is generally applicable to spectra when echo signals are received from many statistically homogeneous small scatterers that are randomly positioned within an analysis region. Such scattering fulfills the basic conditions that: 1. Rf signals exhibit Gaussian statistics; and 2. spectral

Parameterimagestatistics0 F. L. Lrzzr

P D F

ET

1379

AL.

0.15

0.10

0.05

0.05 0

0 -15

-10

-5 0 MIDBAND FIT-MEAN

5

10

-15

-10

-5 0 MIDBAND FIT-MEAN

(dB)

1

5

10

(dB)

D F

0.2

0.1

0.l

0

0

-6

I ;

4 I

SPEC;‘, I

0 2 SLOPE -MEAN (dB /MHz) I

I

4

6

-6

1

I

0.04

P 0.04 D F

0.02

0.02

F

0

-4

4

-2 0 2 SPECTRAL SLOPE -MEAN (dB 1MHz) I

I

I

1

I

6

I

0 -20 SPECTRAL INTERZEPT - MEAN (dB?

Fig. 10. Measured histograms (bars) and theoretical pdfs (solid line) for prostate carcinoma (CA) of Fig. 7c. Top = midband fit; middle = spectral slope; bottom = spectral intercept.

cells constitute independent random variables. An important general finding is that theoretical pdfs for linearregression spectral parameters have shapes that are independent of their mean values. As seen in eqn (9) the pdf for midband fit is dependent on p--p’, which is the difference between parameter values p and mean value p’. Similar pdf dependencies were found for slope [eqn (lo)] and intercept. Our analysis did not assume a particular spectral model, and placed no restriction on spectral shape. Often, relatively small bandwidths are analyzed in clinical examinations, and spectra (in dB) appear as linear functions of frequency. However, linear-regression parameters still provide useful and statistically valid features for nonlinear spectral shapes (Lizzi et al. 1987) In these cases, even though the mean value of each spectral parameter (e.g., /3’) will depend on the specific frequency range being examined, the statistical variations of each spectral parameter about its mean will still be described by our theoretical results. Attenuation can affect the statistics of spectral parameters, other than intercept, in several ways. As described in Results, uniform attenuation in intervening

-40

-20 0 20 SPECTRAL INTERCEPT -MEAN (de)

40

Fig. 11. Measured histograms (bars), theoretical pdfs (solid line) and maximum likelihood Gaussian fits (dashed line) for prostate segment posterior to carcinoma of Fig. 7c. Top = midband fit; middle = spectral slope; bottom = spectralintercept.

tissue will progressively alter the mean values of slope and midband fit, but will not affect the statistical dispersion of these parameters about their local mean values. Attenuation within analyzed regions with large depths will affect these statistical dispersions because the mean parameter values will vary within the analyzed segment. If attenuation coefficients are known or measured, these effects can be compensated as described in our analysis of clinical data. Estimator bias requires a definition of “true” parameter values. Our theoretical framework report (Lizzi et al. 1983) described the “true” calibrated power spectrum as an ensemble average that was approximated by averaging calibrated spectra from many adjacent scan lines before conversion to dB. In this case, each spectral cell in the ensemble power spectrum has a mean value of F, [as in eqn (l)], which becomes 10 log pn upon conversion to dB. Equation (3) shows that a single-line spectrum has a corresponding expected value of 10 log pn -2.5 (dB), so that each spectral cell is biased by a frequency-independent factor of -2.5 dB. Accordingly, the midband fit and

i ix0

Ultrasound

in Medicine

and Biology

Volume

23. Number

0. 1997

Table 2. Theoretical and observed standard deviations of spectral parameters. Midband Theoretical Healthy liver (Fig. 7a) Lwer metastatic carcinoma (Fig. 7b) Prostate carcinoma (Fig. 7c) Prostate posterior to carcinoma (Fig. 7~)

fit (dB)

Spectral

Measured

Theoretical

slope

CdBIMHz) Measured

Spectral ‘l‘heoretical

intercept tdBi ..___.Measured

2.7

3.0

4.98

1.x3

15.1

14.3

2.7 2.4

7.8 2.9

4.98 I .72

4.10 1.99

15.3 10.3

13.7 14.0

2.4

4.6

1.72

1.98

IO.3

x.4

intercept also are biased by -2.5 dB, as indicated in eqn (11). Slope estimates are unbiased because the constant 2.5 dB offset does not affect spectral shape; this can be shown by evaluating eqn (12) and noting that the summation of the second bracketed term (containing the 2.5 dB bias) is equal to zero. Results presented in this paper describe how the probability density functions of spectral parameter estimators depend on frequency (f, and B) and analysis (L) parameters. These results were obtained for the conditions used in spectral parameter imaging in which unaveraged parameters are computed and displayed. The results can now be used for designing image-processing procedures, examining tissue heterogeneity and determining the precision of other spectral descriptors (e.g., integrated backscatter) and derived estimates of scatterer properties. Our theoretical results show that the pdf of midband fit depends only on N = BL, and its SD is inversely proportional to a. The pdf of spectral slope depends only on N and Af = l/L, and its SD is inversely proportional to X6% (for BL greater than about 10). Thus, increasing bandwidth is an especially important means for improving estimator precision. The pdf for spectral intercept depends upon B, L and the center frequency f,; if, as usual, B is much smaller than 3.5 f,, the intercept SD is approximately equal to the slope SD multiplied by f;.. As described above, if the spectrum has a nonlinear shape, the mean value of each parameter will depend upon the frequency band being analyzed, but the statistical precision of each estimate will still follow these dependencies on L, B and f,. As N becomes large, all of the pdfs can be approximated as Gaussian functions using the SDS in Table 1. For spectral slope, this Gaussian approximation can be made even at relatively small values of BL, as seen Fig. 5. For midband fit, higher values of N are required, as seen in Fig. 3. The current analysis can be expanded to treat cases in which spatial low-pass filtering is employed to improve the quality of parameter images. If a 2D spatial

filter locally averages parameters from V-independent range cells and M-independent scan lines, the resultant pixel value is the average of KM-independent samples, and the SDS for each parameter are those of Table 1 divided by X6%?. The pdf for each averaged parameter can then be evaluated by raising that parameter’s characteristic function [as in eqn (9)] to the VM power and computing an inverse Fourier transform. Equation (5) indicates that the pdf for averaged midband fit values can be evaluated by simply replacing N with the product VMN and applying the results shown in Fig. 2. Such procedures permit analyses of the trade-offs among parameter image attributes, including image stability (estimator precision), axial resolution (equal to VL/2 for a Hamming window) and lateral resolution (M X the effective beamwidth). The results of the current report are directly applicable to the statistics of integrated backscatter (ZB) (Thomas et al. 1989). Unaveraged values of ZB, suitable for forming images, are computed by dividing unaveraged power spectra lS(f>12by a deterministic frequencydependent correction factor C(f) that accounts for beam cross-section. The resultant function is converted to dB, and is summed over the useable bandwidth B. Accordingly, ZB is the sum over n of -10 log C(f,) + 10 loglS(&)I’. The first term sums to a deterministic constant, and the second term is equal for midband fit (eqn. 5). Thus, the statistics of ZB about its mean value are the same as those we have derived for midband fit (about its mean value) in this report. The theoretical pdfs and SDS agreed well with clinical data from homogeneous segments in liver and prostate (Figs. 8-10; first 3 cases in Table 2). These analyzed segments were selected to test the validity of theoretical results in homogeneous tissue regions, and were not intended to provide general conclusions for liver and prostate scattering. Although agreement with theoretical results was very good, we believe that some of the differences between midband-fit histograms and theoretical pdfs may be associated with local heterogeneity or with the presence of coherent scatterers (whose video

Parameter

image

statistics

signals would produce Rician, rather than Rayleigh, statistics). Variations in attenuation coefficients may also contribute to these differences. We are expanding our analysis to examine these topics, and to determine if other histogram parameters (e.g., skewness) may serve as quantitative indices of such tissue attributes. Results for the posterior prostate segment (Fig. 11 and last case in Table 2) demonstrate the effects of tissue heterogeneity, which is apparent in the image of Fig. 7c. Figure 11 shows that the midband fit pdf is essentially a Gaussian function, rather than the theoretical pdf for homogeneous scatterers. This result may be associated with the presence of different scatterer species with different mean values. In such cases, a knowledge of the theoretical pdf for a homogeneous population may help evaluate the spread of mean values within the scatterer population. The results in this report can also be employed to evaluate the precision of estimates of effective scatterer size and acoustic concentration (product of scatterer concentration and the square of relative acoustic impedance with respect to the surround) (Lizzi et al. 1987). Size is estimated as a function of spectral slope, and its pdf can be evaluated based on eqn ( 12). Acoustic concentration is estimated from slope and intercept, and its pdf can also be evaluated using results of this report. Such analyses would help in separating measurement variance from tissue-property variance, an important factor in tissue classification. We are currently applying our results to the preceding topics to help improve the design of appropriate systems and analysis procedures, to provide guidance in interpreting clinical data from several organs being investigated in our laboratories and to evaluate lesion detectability. AcknoM,lrd~emenu-This work was supported in part by NIH Research Grants ROlCA53561 and ROl-CA38400, awarded by the National Cancer Institute, and ROI-EY01212, awarded by the National Eye Institute. We acknowledge the essential contributions of Drs. William Fair and Donald King and their staffs at the Memorial Sloan-Kettering Cancer Center and the Columbia University College of Physicians and Surgeons. We also gratefully acknowledge the talented assistance of Roslyn Raskin in preparing this manuscript.

0 F. L. LIZZI ET

AL.

1381

Gradshteyn IS. Ryzhik IM. Table of integrals, series and products. New York: Academic Press, 1980576-578. Huisman HJ, Thijssen JM. Precision and accuracy of acoustospectrographic parameters, Ultrasound Med Biol 1996;22:85587 I. Insana MF, Wagner RF, Brown DG, Hall TJ. Describing small-scale structure in random media using pulse-echo ultrasound. J Acoust Sot Am 1990;87: 179-192. Kaplan D, Ma Q. On the statistical characteristics of log-compressed Rayleigh signals. Proc IEEE Ultrason Sympos, 1994:961I974. Kuc R, Taylor K. Variation of acoustic attenuation coefficient slope estimates for in viva liver. Ultrasound Med Biol 1982:8:403-4 12. Lindgren BW. Statistical theory. New York: McGraw-Hill. 1962:3190. Lizai FL, Astor M, Liu T, Deng C, Coleman DJ. Silverman RH. Ultrasonic spectrum analysis for tissue assays and therapy evaluation. Int J Imaa Svst Tech 1997a:8:3-I()” ..,. Lizzi FL, Feleppa EJ. Astor M, Kalisz A. Statistics of ultrasonic parameters for prostate and liver examinations. IEEE Tram UFFC 1997b:44:935-942. Lizzi FL, Greenebaum M, Feleppa EJ. Elbaum M. Coleman DJ. Theoretical framework for spectrum analysis in ultrasonic tissue characterization. J Acoust Sot Am 1983:73: 136661373. Lizzi FL. King DL. Rorke HC, Hui J, Ostromogilsky M. Yaremko MM, Feleppa EJ. Wai P. Comparison of theoretical scattering results and ultrasonic data from clinical liver examinations. Ultrdsound Med Biol 1988:14:377-385. Lizzi FL. Laviola MA, Coleman DJ. Tissue sienature characterization utilizing frequency domain analysis, Proc IEEE Ulnason Sympos 1976:76:714-719. Lizzi FL, Ostromogilsky M. Feleppa EJ, Rorke MC. Yaremko MM. Relationship of ultrasonic spectral parameters to features of tissue microstructure. IEEE Trans UFFC 1987;34:319-329. Narayana P. Ophir J. On the validity of the linear approximation in the parametric measurement of attenuation in tissues. Ultrasound Med Biol I983;9:357-361. Papoulis A. Probability. random variables and stochastic processes. New York: McGraw-Hill, 1965:244. Parker KJ. Attenuation measurement uncertainty caused by speckle statistics. J Acoust Sot Am 1988;80:727-734. Thomas LJ III. Barzilai B, Perez JE. Quantitative real-time imaging of myocardium based on ultrasonic integrated backscatter. IEEE Tram UFFC 1989;36:466-470. Tuthill TA. Sperry RH. Parker KJ. Deviations from Rayleigh statistics in ultrasonic speckle. Ultrason Imaging 1988:lO:S l-89. Wagner RF, Insana MF. Brown DG. Unified approach to the detection and classification of speckle texture in diagnostic ultrasound. Optical Eng 1986;25:732-742. Zagzebski LA. Lu ZF, Xao LX. Quantitative ultrasound imaging: in vi\‘0 results in normal liver. Ultrason Imaging 1993: I5:3?5-35 I

APPENDIX This appendix computes the mean and standard deviation of the pdf for & described in eqn (2). The mean is computed by evaluating the integral

/Z,x(Z,, )dZ,,.

REFERENCES Chaturvedi P, Insana MF. Error bounds on ultrasonic scatterer size estimates. J Acoust Sot Am 1996;100:392-399. Davenport WB. Root WL. An introduction to the theory of random signals and noise. New York: McGraw-Hill, 1958: 113-167. Feleppa EJ, Lizzi FL, Coleman DJ. Yaremko MM. Diagnostic spectrum analysis in ophthalmology: A physical perspective. Ultrasound Med Biol 1986; 12:623-63 1. Feleppa EJ, Fair WR. Kalisz A, Sokil-Melgar JB, Lizzi FL, Rosado A. Shao MC, Liu T, Wang Y. Cookson M, Reuter V. Typing of prostate tissue by ultrasonic spectrum analysis. IEEE Tram UFFC 1996:43:609-619.

Let .Y = exp (Z,,/4.34

- In p,,): then using eqn (2). eqn (Al)

(Al)

becomes:

I 2, = 4.34

I,

In(x)

exp( -x)&

4 4.341npL,,

(A2)

The integral in eqn (A2) is evaluated using Gradshteyn and Ryzhik (1980). eqn (4.352.1) as 4.34+ (I), where clr(M) is the psi function defined using Euler’s constant C = 0.577:

1382

Ultrasound

in Medicine

1/&I4)=[1+1/2+1/3+...+Il(M-1)-C];

and Biology M>l

l)(l)

(A3)

Volume

23, Number

After substituting using Gradshteyn

9, 1997 Y = x /.L~ in eqn (A5), the integral and Ryzhik (1980). eqn. (4.358.2)

= -c. uza2 = 4.34’[(1)(1)

t

- ln( l/p,J)’

Thus, r/~ is evaluated as 4.34 ln p,, - 4.34C, which is equal to 4.34 In CL, - 2.50. The variance of eqn (2) is computed by evaluating

5(2.1)]

~ 3.34’[($(1) mz,, = 4.34.

CI az”2 = z:g(z”)dz, - ‘2:. I -* Again,

letting

x = exp (ZJ4.34

- In CL,), eqn (A4)

(A41

where

5 is Riemann’s

E(~n(xpn))2 I 0

exp(-x)dx

[&2,1

)]I’?

(Ah) iA6a)

(A7)

becomes

- 2:.

--. Int,l/&L.i]~

zeta function.

t(2.1) CL2 = 4.342

can be evaluattxl with the result:

= 7r2/6.

(A5) Thus,

oz. is evaluated

as 4.34

n/W’6

or 5.56

(A7a)