Exponential-spline-based Features for Ultrasound Signal ...

Tissue Characterization Ultrasound Signals Modeling Exponential Splines From Continuous to Discrete Conclusions Ac

Exponential-spline-based Featuresfor Ultrasound Signal Characterization

Simona Maggio

Department of Electronics, Computer Science and Systems (DEIS)University of Bologna

Scuola di Dottorato: Scienze e Ingegneria dell’InformazioneCorso di Dottorato: Ingegneria Elettronica, Informatica e delle Telecomunicazioni

19/10/2009

1 / 29

Outline

1 Tissue Characterization for Ultrasound Diagnostic

2 Ultrasound Signals ModelingPredictive DeconvolutionContinuous/Discrete Approach

3 Exponential Splines

4 From Continuous to DiscreteE-spline-based Whitening ModelNon Stationary Discrete Whitening

5 ConclusionsAnalysis ResultsFurther Developments

2 / 29

Outline

2 / 29

Outline

2 / 29

Outline

2 / 29

Outline

2 / 29

Ultrasound Imaging

• Pros: real time, non invasive.

• Limits: low resolution.• Tissue characterization:

• Highligthing features invisible by visualinspection

• Cancer detection and staging• Avoiding unnecessary biopsy

• Specific application: Prostate cancercomputer-aided detection

3 / 29

Ultrasound Imaging

• Pros: real time, non invasive.

• Limits: low resolution.• Tissue characterization:

• Highligthing features invisible by visualinspection

• Cancer detection and staging• Avoiding unnecessary biopsy

• Specific application: Prostate cancercomputer-aided detection

3 / 29

Procedure for Tissue Characterization

Figure: Feature extraction for ultrasound analysis

4 / 29

From Tissue to Echo

• Discrete model: z = HΣs + n = Hx + n• Σ: coherent reflection, macroscopic interactions, mean

value of diffused field.• s: incoherent reflections, interactions smaller than

wavelength, random fluctuations of diffused field.• Scattering as generalized gaussian noise:

p(s|µ, b) = a · e−|s−µ

5 / 29

From Tissue to Echo

p(s|µ, b) = a · e−|s−µ

5 / 29

From Tissue to Echo

p(s|µ, b) = a · e−|s−µ

5 / 29

From Tissue to Echo

p(s|µ, b) = a · e−|s−µ

5 / 29

Deconvolution to recover Tissue Response

• Convolutional model: z[k] = x[k] ∗ h[k] + n[k]

• Deconvolution to restore tissue response: x[k]

• Point Spread Function (PSF) h[n] not known

• Blind adaptive deconvolution approach 1

• Advantages: simplicity, low computational cost,variable PSF.

1[Ng et al., 2007], [Michailovich and Adam, 2005],[Jensen, 1994],[Rasmussen, 1994]

6 / 29

Predictive Deconvolution (PD)

• W(z) predictive filter to remove predictable features of z[n]

• Ideal reconstruction if the transducer is an AR system andx[n] is white Gaussian → whitening

• Recursive Least Squares solution for adaptive whitening• z[n] as non stationary AR process and x[n] generalized

Gaussian• Recovering unpredictable part of RF signal: scattering

7 / 29

Predictive Deconvolution (PD)

• W(z) predictive filter to remove predictable features of z[n]

• Ideal reconstruction if the transducer is an AR system andx[n] is white Gaussian → whitening

• Recursive Least Squares solution for adaptive whitening• z[n] as non stationary AR process and x[n] generalized

Gaussian• Recovering unpredictable part of RF signal: scattering

7 / 29

Improvement due to Predictive Deconvolution

Prostate Images

No Preprocessing RLS Deconvolution

SE 0.69 ± 0.06 0.75 ± 0.09SP 0.94 ± 0.02 0.93 ± 0.01Acc 0.93 ± 0.02 0.93 ± 0.02Az 0.92 ± 0.02 0.95 ± 0.02

8 / 29

Benignant case

SE 0.69 ± 0.06 0.75 ± 0.09SP 0.94 ± 0.02 0.93 ± 0.01Acc 0.93 ± 0.02 0.93 ± 0.02Az 0.92 ± 0.02 0.95 ± 0.02

8 / 29

Malignant case

SE 0.69 ± 0.06 0.75 ± 0.09SP 0.94 ± 0.02 0.93 ± 0.01Acc 0.93 ± 0.02 0.93 ± 0.02Az 0.92 ± 0.02 0.95 ± 0.02

8 / 29

Continuous/Discrete Approach for PD

• Continuous whitening• g(t) continuous RF signal in a stationary case• w(t) continuous uncorrelated unpredictable (scattering)• Link with discrete: discretization in a spline basis• Exponential splines: multiresolution version of L~α

• Identification of continuous model from sampled data• Shift-variant multiresolution description of scattering

signal

9 / 29

signal

9 / 29

signal

9 / 29

Exponential Splines

• E-splines2: s(t) =∑

k∈Zakρ~α(t − tk) + p0(t)

• ρ~α green function of L~α

• s(t) reconstructed in exponential B-spline basis:s(t) =

k∈Zc[k]β~α(t − k) β~α(t) = ∆~α{ρ~α(t)}

• Scale T for generic sampling time: β~αT(t) = ∆~αT{ρ~α(t)}

2[Khalidov and Unser, 2006], [Unser and Blu, 2005a],[Unser and Blu, 2005b]

10 / 29

Exponential Splines

10 / 29

Exponential Splines

10 / 29

E-Spline Multi-Resolution Analysis

• Dyadic scale for multiresolution analysis:T = 2i

• Scaling function: ϕi(t) = β2i(t)/‖β2i‖L2

• At scale i it is not a dilated version of thescaling function at scale 0.

• Mallat’s filter bank algorithms butprecomputation of filters for each scale

• Main property: wavelet coefficients of asignal f (t) are samples of smoothedversion of L~α(f )

11 / 29

E-Spline Multi-Resolution Analysis

• Dyadic scale for multiresolution analysis:T = 2i

• Scaling function: ϕi(t) = β2i(t)/‖β2i‖L2

• At scale i it is not a dilated version of thescaling function at scale 0.

• Mallat’s filter bank algorithms butprecomputation of filters for each scale

• Main property: wavelet coefficients of asignal f (t) are samples of smoothedversion of L~α(f )

11 / 29

E-spline-based Whitening Model

• Continuous: Rgg(t) =(ρα1(t) ∗ ρα1(−t)) ∗ · · · ∗

ραp(t) ∗ ραp(−t))

• Discrete:

R(z) =zp · B(−~α:~α)(z)

∏pi=1 e−αi z(1 − eαi z)(1 − eαi z−1)

∏p−1i=1

−1ζi

(1 − ζiz)(1 − ζiz−1)∏p

i=1 e−αi(1 − eαi z)(1 − eαi z−1)

• Inter-dependence: ζi = f (~α)

• Whitening filter: W(z−1) =

p−1Y

(1 − ζiz−1

(1 − eαi z−1)

(stable)

• ARMA(p, p − 1) discrete model12 / 29

• Discrete:

R(z) =zp · B(−~α:~α)(z)

∏pi=1 e−αi z(1 − eαi z)(1 − eαi z−1)

∏p−1i=1

−1ζi

(1 − ζiz)(1 − ζiz−1)∏p

i=1 e−αi(1 − eαi z)(1 − eαi z−1)

p−1Y

(1 − ζiz−1

(1 − eαi z−1)

(stable)

• Discrete:

R(z) =zp · B(−~α:~α)(z)

∏pi=1 e−αi z(1 − eαi z)(1 − eαi z−1)

∏p−1i=1

−1ζi

(1 − ζiz)(1 − ζiz−1)∏p

i=1 e−αi(1 − eαi z)(1 − eαi z−1)

p−1Y

(1 − ζiz−1

(1 − eαi z−1)

(stable)

• Discrete:

R(z) =zp · B(−~α:~α)(z)

∏pi=1 e−αi z(1 − eαi z)(1 − eαi z−1)

∏p−1i=1

−1ζi

(1 − ζiz)(1 − ζiz−1)∏p

i=1 e−αi(1 − eαi z)(1 − eαi z−1)

p−1Y

(1 − ζiz−1

(1 − eαi z−1)

(stable)

• Discrete:

R(z) =zp · B(−~α:~α)(z)

∏pi=1 e−αi z(1 − eαi z)(1 − eαi z−1)

∏p−1i=1

−1ζi

(1 − ζiz)(1 − ζiz−1)∏p

i=1 e−αi(1 − eαi z)(1 − eαi z−1)

p−1Y

(1 − ζiz−1

(1 − eαi z−1)

(stable)

Exponential Parameters Estimation

• AR(p): Yule-Walker equations.Only poles ai = eαi

• ARMA(p,p-1): No pole-zerointer-dependence.

• Annihilating Polynomial:coeff of AP of Rgg(k) arefunctions of ~α:{∏p

1 − aiz−1)}

Rgg(k) = 0

• MA correction:

• AR to compute poles• pole-zero inter-dependence to estimate zeros•

1MA(z) filtering to cancel zeros

• new AR computation for better pole estimation• possibly iterations

13 / 29

1 − aiz−1)}

Rgg(k) = 0

• MA correction:

13 / 29

1 − aiz−1)}

Rgg(k) = 0

• MA correction:

13 / 29

1 − aiz−1)}

Rgg(k) = 0

• MA correction:

13 / 29

1 − aiz−1)}

Rgg(k) = 0

• MA correction:

13 / 29

Restoring Non Stationarity

• Non stationary RF signal

• Shift-variant whitening to get the non stationary scattering• Two possibilities:

• Sliding window• Recursive parameter estimation

• Recursive Annihilating Polynomial vs LS solution

• Time dependent parameters: ~α[n]

• Shift-variant E-spline wavelet

• Shift-variant multiresolution analysis (sliding window)

• Piece-wise multiresolution description of scattering

14 / 29

Phantoms Targets

• Hypo and hyper echoic target detection

• Linear classifier: learning on ±9 dB, testing on ±6 dB

• Comparison with traditional predictive Deconvolution:improvement 40.3%

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1−Specificity

ROC curves: hypoechoic target detection

ESW − 1 featPD

15 / 29

Phantoms Targets

• Hypo and hyper echoic target detection

• Linear classifier: learning on ±9 dB, testing on ±6 dB

• Comparison with traditional predictive Deconvolution:improvement 2.1%

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1−Specificity

ROC curves: hyperechoic target detection

ESW − 1 featPD

15 / 29

Real Data

• Real data cancer detection: prostate trans-rectal• Dataset: 15 benignant cases, 22 malignant cases• Linear classifier: training 18 cases, testing 19 unknown img• Comparison with traditional predictive Deconvolution:

improvement 16.9%

16 / 29

Conclusions

• Predictive Deconvolution to restore US tissue response

• Continuous approach to predictive deconvolution

• Discretization in Exponential spline basis

• Identifications of E-splines tuned on US signal

• Multiresolution description of tissue response

• Improvement in detection performance: 16.9%

Further Developments• Improve parameter estitmation

• Better MA correction algorithm

• E-spline for texture analysis

17 / 29

Conclusions

• Predictive Deconvolution to restore US tissue response

• Continuous approach to predictive deconvolution

• Discretization in Exponential spline basis

• Identifications of E-splines tuned on US signal

• Multiresolution description of tissue response

• Improvement in detection performance: 16.9%

Further Developments• Improve parameter estitmation

• Better MA correction algorithm

• E-spline for texture analysis

17 / 29

Thank you for your attention!

http://mas.deis.unibo.it/

18 / 29

Bibliography

Jensen, J. (1994).

Estimation of in vivo pulses in medical ultrasound.Ultrasonic Imaging, 16:190–203.

Khalidov, I. and Unser, M. (2006).

From differential equations to the constructio of new wavelet-like bases.IEEE Transactions on Signal processing, 54(4).

Michailovich, O. V. and Adam, D. (2005).

A novel approach to the 2-d blind deconvolution problem in medical ultrasound.IEEE Transactions on Medical Imaging, 24(1):86–104.

Ng, J., Prager, R., Kingsbury, N., Treece, G., and Gee, A. (2007).

Wavelet restoration of medical pulse-echo ultrasound images in an em framework.IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 54(3):550–568.

Rasmussen, K. (1994).

Maximum likelihood estimation of the attenuated uktrasound pulse.IEEE Transactions on Signal Processing, 42:220–222.

Unser, M. and Blu, T. (2005a).

Cardinal exponential splines: Part i - theory and filtering algorithms.IEEE Transactions on Signal Processing, 5373(4).

Unser, M. and Blu, T. (2005b).

Cardinal exponential splines: Part ii - think analog, act digital.IEEE Transactions on Signal Processing, 53(4).

19 / 29

Publications

• 2009 (to be published in next issue)IEEE Transactions on Medical ImagingS. Maggio, A. Palladini, L. De Marchi, M. Alessandrini, N. Speciale, G. MasettiPredictive deconvolution and hybrid feature Selection for Computer-Aided Detection of prostate cancer

• 2009 MarchProceedings International Symposium on Acoustical ImagingM. Scebran, A. Palladini, S. Maggio, L. De Marchi, N. SpecialeAutomatic regions of interests segmentation for computer aided classification of prostate TRUS images

• 2008 NovemberProceedings IEEE IUS2008S. Maggio, L. De Marchi, M. Alessandrini, N. SpecialeComputer aided detection of prostate cancer based on GDA and predictive deconvolution

• 2005 NovemberWSEAS Transactions on SystemsS. Maggio, N. Testoni, L. De Marchi, N. Speciale, G. MasettiUltrasound Images Enhancement by means of Deconvolution Algorithms in the Wavelet Domain

• 2005 SeptemberWSEAS Int. Conf. on Signal Processing, Computational Geometry and Artificial Vision (ISCGAV)S. Maggio, N. Testoni, L. De Marchi, N. Speciale, G. Masetti

Wavelet-based Deconvolution Algorithms Applied to Ultrasound Images

20 / 29

Texture analysis through E-splines

• Global learning of exponential parameters• E-spline wavelet tuned on ultrasound signal• E-spline wavelet transform for texture typing• Comparison with traditional wavelets• Improvement with respect to energy information: 14.9%

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1−Specificity

ROC curves: cancer detection

NakagamiNaka + Variance ESW

21 / 29

Comparison with previous works

Table: Published methods for ultrasound-based prostate tissuecharacterization

Work Ground Truth Technique Results %# ROIs Features SE SP Acc Az

Basset 37 Textural 83 71 - -Huynen - Textural 80 88.20 - -Houston 25 Textural 73 86 80 -Schmitz 3405 Multi 82 88 - -Scheipers 170 484 Multi - - 75 86Feleppa 1019 Spectral - - 80 85Mohamed 96 Textural 83.3 100 93.75 -Llobet 4944 Textural 68 53 61.6 60.1Mohamed 108 Multi 83.3 100 94.4 -Han 2000 Multi 92 95.9 - -Maggio 58602 Multi on RLS 75 93 93 95

22 / 29

Autocorrelation Function Discretization 1/4

In the first order case the ACF of g(t) can be obtained as:

Rgg(τ) = F−1 {Φg}

= F−1

1jω − α

∗ F−1

1−jω − α

= eαtu(t) ∗ e−αtu(−t) = ρα(t) ∗ ρα(−t) (1)

Green function of Lα: ρα(t) = eαtu(t) =

+∞∑

pα[k]βα(t − k)

Rgg(t) =

+∞∑

pα[k]βα(t − k) ∗+∞∑

k′=0

pα[k′]βα(−t − k′)

+∞∑

pα[k]+∞∑

k′=0

pα[k′]βα(t − k) ∗ βα(−t − k′)

=+∞∑

pα[k]+∞∑

k′=0

pα[k′]eαβα(t − k) ∗ β−α(t + k′ + 1) (2)

23 / 29

The expression for ACF simplifies to

Rgg(t) =

+∞∑

pα[k]∑

k′′∈Z

pα[k − k′′ − 1]eαβ(−α:α)(t − k′′)

k′′∈Z

(pα ∗ p′α)[k′′ + 1]eαβ(−α:α)(t − k′′) (3)

where p′α[k] = pα[−k].It turns out that the z-transform of the discrete signal ACF isgiven as:

R(z) = eαzPα(z)P′α(z)B(−α:α)(z)

=eαzB(−α:α)(z)

(1 − eαz−1)(1 − eαz)(4)

24 / 29

The extension to higher-order AR models follows naturally asshown below:

Rgg(t) = F−1

1jω − α1

∗ · · · ∗ F−1

1jω − αp

= (ρα1(t) ∗ ρα1(−t)) ∗ · · · ∗(

k1∈Z

(pα1 ∗ p′α1)[k1 + 1]eα1β(−α1:α1)(t − k1)

∗ · · ·

kp∈Z

(pαp ∗ p′αp)[kp + 1]eαpβ(−αp:αp)(t − kp)

25 / 29

The final expression of the z-transform of Rgg(k) is:

R(z) =

eαi zPαi(z)P′αi

(z) · B(−~α:~α)(z)

eαi z(1 − eαi z)(1 − eαiz−1)

· zp · B(−~α:~α)(z) (6)

26 / 29

Whitening Filter Stability 1/3

The proof that the z-transform of B-spline β(−~α:~α), B(−~α,~α)(z),doesn’t have zeros on unit circle is shown by the followingconsiderations:

β̂(−~α,~α)(ω) = F [β~α ∗ β−~α(t)]

e−αiF [β~α ∗ β~α(t − p)]

e−αi e−jωp∣

∣β̂~α(ω)∣

27 / 29

B(−~α,~α)(z) is linked to the Fourier transform of samples ofβ(−~α,~α)(t), β̂s,(−~α,~α)(ω):

β̂s,(−~α,~α)(ω) =∑

β̂(−~α,~α)(ω + 2πk)

e−αi e−jωp∑

∣β̂~α(ω + 2πk)

e−αi e−jωpA~α(ejω) (8)

28 / 29

where A~α(ejω) is the discrete Fourier transform of the Gramsequence of B-splines, and, for the Riesz basis property, it isalways grater than zero:

0 < r2~α < A~α(ejω) < R2

~α < +∞ (9)

As a consequence β̂s,(−~α,~α)(ω) > 0 for every ω and, sinceB(−~α,~α)(z) = β̂s,(−~α,~α)(ω)|e−jω=z, B(−~α,~α)(z) doesn’t have anyzeros on the unit circle. This proof and the consideration aboutreciprocal roots of B(−~α,~α)(z) guarantee the whitening filterstability.

29 / 29

Exponential-spline-based Features for Ultrasound Signal ...

Documents

Transcript of Exponential-spline-based Features for Ultrasound Signal ...