1. technology and computing
2. software
3. databases

ALMA MATER STUDIORUM - UNIVERSITÀ DI BOLOGNA
PhD COURSE: ELECTRONICS, COMPUTER SCIENCE AND TELECOMMUNICATIONS
XXIII CYCLE - SCIENTIFIC DISCIPLINARY SECTOR: ING-INF/01
Analysis and Modeling
Techniques for Ultrasonic
Tissue Characterization
Simona Maggio
S UPERVISOR
Professor Guido Masetti
C OORDINATOR
Professor Paola Mello
DEIS - DEPARTMENT OF ELECTRONICS, COMPUTER SCIENCE AND SYSTEMS
JANUARY 2008 - DECEMBER 2010
A BSTRACT
i
Abstract
This thesis introduces new processing techniques for
computer-aided interpretation of ultrasound images with
the purpose of supporting medical diagnostic. In terms of
practical application, the goal of this work is the improvement of current prostate biopsy protocols by providing physicians with a visual map overlaid over ultrasound images
marking regions potentially affected by disease. As far as
analysis techniques are concerned, the main contributions
of this work to the state-of-the-art is the introduction of deconvolution as a pre-processing step in the standard ultrasonic tissue characterization procedure to improve the diagnostic significance of ultrasonic features.
This thesis also includes some innovations in ultrasound
modeling, in particular the employment of a continuoustime autoregressive moving-average (CARMA) model for ultrasound signals, a new maximum-likelihood CARMA estimator based on exponential splines and the definition of
CARMA parameters as new ultrasonic features able to capture scatterers concentration.
Finally, concerning the clinical usefulness of the developed techniques, the main contribution of this research is
showing, through a study based on medical ground truth,
that a reduction in the number of sampled cores in standard
prostate biopsy is possible, preserving the same diagnostic
power of the current clinical protocol.
Keywords: ultrasound, tissue characterization, deconvolution, prostate biopsy, continuous-time autoregressive moving average, exponential splines.
ii
R ÉSUMÉ
Résumé
Cette thèse introduit de nouvelles techniques de traitement pour l’interprétation guidée d’images ultrasons, dans
le but de soutenir le diagnostic médical. L’objectif pratique
de ce travail est l’amélioration du protocole standard pour la
biopsie de la prostate, en fournissant au médecin une carte
visuelle sur l’écographie qui marque les régions potentiellement malignes. En ce qui concerne les techniques d’analyse,
la contribution principale de cette thèse est l’introduction
de la déconvolution comme étape de pré-traitement dans la
procédure standard de caractérisation des tissus par ultrasons, afin d’améliorer la valeur diagnostique des caractéristiques du signal.
Cette thèse présente, en outre, des innovations dans la
modélisation du signal ultrason, en particulier la proposition d’un modèle autorégressif à moyenne glissante en temps
continu (CARMA), le développement d’un nouvel estimateur
au maximum de vraisemblance de paramètres CARMA fondé
sur les splines exponentielles et la définition des paramètres
CARMA comme des nouvelles caractéristiques du signal ultrason capables de capturer l’information sur la concentration de scatterers.
Enfin, concernant l’utilité clinique des techniques développées, la contribution principale de cette recherche est
dans la démonstration, au travers d’une étude médicale, de
la possibilité de réduire concrètement le nombre d’ échantillons prélevés pendant la biopsie de la prostate standard,
en préservant le même pouvoir diagnostique du protocole
biopsique standard.
S OMMARIO
iii
Sommario
Questa tesi introduce nuove tecniche di elaborazione per
l’interpretazione guidata di immagini ultrasoniche, allo
scopo di supportare la diagnosi medica. Dal punto di vista
applicativo, l’obiettivo di questo lavoro è il miglioramento
del protocollo standard per la biopsia prostatica,fornendo al
radiologo una mappa visiva di regioni potenzialmente malate. Per quan-to riguarda le tecniche di analisi, il contributo principale di questa tesi è l’introduzione della deconvoluzione come passo di pre-processing nella procedura
standard di tissue characterization basata su ultrasuoni, al
fine di migliorare il potere diagnostico di alcune caratteristiche estratte dal segnale.
Questa tesi include, inoltre, delle innovazioni nel modeling del segnale ultrasonico, in particolare la proposta di
un modello tempo-continuo autoregressive moving-average
(CARMA), lo sviluppo di un nuovo stimatore maximum-likelihood di parametri CARMA basato su exponential spline e
la definizione dei parametri CARMA come nuove caratteristiche del segnale ultrasonico capaci di catturare l’informazione sulla concentrazione di scatterer.
Infine, per quanto riguarda l’utilità applicativa delle tecniche sviluppate, il contributo principale di questa ricerca
sta nel mostrare attraverso uno studio su ground-truth medico che è possibile ridurre effettivamente il numero di campioni prelevati durante la biopsia prostatica standard,preservando le stesse capacità diagnostiche del protocollo clinico
standard.
C ONTENTS
Abstract
i
Résumé
ii
Sommario
iii
Introduction
1
1
Background
1.1 Diagnostic ultrasound . . . . . . . . . . . . . . . . .
1.2 Tissue characterization . . . . . . . . . . . . . . . . .
1.3 Application: prostate TRUS . . . . . . . . . . . . . .
5
6
10
19
2
Ultrasound image analysis and interpretation
2.1 Dimensionality reduction . . . . . . . . .
2.2 Feature-based segmentation . . . . . . . .
2.3 Deconvolution as pre-processing step . .
2.4 Feature-based classification . . . . . . . .
2.5 Ground truth database . . . . . . . . . . .
2.6 Experimental results . . . . . . . . . . . .
2.7 Conclusion . . . . . . . . . . . . . . . . . .
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
27
29
32
36
39
40
43
51
A CARMA model for ultrasound signals
3.1 CAR and CARMA models . . . . . . . . . . . . . .
3.2 CARMA processes and exponential splines . . .
3.3 Insights on the exponential spline framework . .
3.4 A new approach to CARMA model identification
3.5 CARMA estimation algorithm . . . . . . . . . . .
3.6 Simulation results . . . . . . . . . . . . . . . . . .
3.7 CARMA for ultrasound tissue characterization .
3.8 Comparative study on phantom images . . . . .
3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
53
55
58
68
70
79
86
92
95
99
v
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
vi
4
C ONTENTS
Improving standard prostate biopsy protocol
4.1 Database collection . . . . . . . . . . . . .
4.2 Processing and learning scheme . . . . .
4.3 Classification results . . . . . . . . . . . .
4.4 The future . . . . . . . . . . . . . . . . . . .
4.5 Conclusion . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
101
102
106
111
111
113
Conclusions
115
A Features for ultrasonic tissue typing
119
B State of the art in tissue characterization
123
C Initial conditions in CARMA estimation
125
D Non symmetric exponential B-splines
127
Notations
131
Publications
137
Bibliography
139
I NTRODUCTION
Research is what I’m doing when I don’t know
what I’m doing.
Wernher von Braun (1912-1977)
U
LTRASOUND imaging is one of the most used technologies for
non invasive medical diagnostic. Resolution limits and speckle
noise affecting ultrasonography prevent the employment of standard image processing techniques to perform tissue characterization supported by computer. On the contrary the interpretation of
ultrasound image requires an approach based on features specific
of the echo signal, able to capture the state of the imaged tissue.
The realization of computer-aided tissue characterization schemes in ultrasound imaging is constrained to the analysis of characteristics of the ultrasound signal and their correlation to the pathological state of tissues. In addition, the study of new parametric
models for the ultrasonic echo allows the definition of new image
attributes which may be likewise correlated to the state of the investigated tissue and useful for characterization. Finally, the design of a medical ground truth and a learning scheme to recognize
pathological tissues is essential for an unbiased approach to tissue
characterization.
In this thesis all the fundamental steps for the realization of
ultrasonic tissue characterization schemes are taken into account
and dealt with through the development of new processing methods. Although the developed techniques are applicable to the investigation of any biological tissue, the final goal of this work is the
design of a computer-aided detection scheme to support physicians in the diagnosis of prostate cancer.
This thesis is subdivided in three stages. First, a preliminary
study of the existing ultrasonic features is performed to select the
most significant attributes for tissue representation. The ground
truth used in this first stage is a collection of ultrasound images
1
2
I NTRODUCTION
of the prostate gland, where the pathological region is outlined by
pathologists. Then, a non linear learning approach is validated on
the collected medical ground truth. Within the first stage the novel
contribution of this thesis is the introduction of shift-variant deconvolution as a pre-processing step in the traditional tissue characterization scheme. In this context deconvolution breaks away
from its traditional goal of improving the visual quality of ultrasound images and is used to enhance the diagnostic power of ultrasonic features. Experimental results show that performing deconvolution before feature extraction yields to significantly increase
the ability of some ultrasonic features to recognize pathological
tissues.
In the second stage, a new model is proposed for the ultrasound echo signal, in order to enable the definition of novel attributes useful for tissue characterization. The ultrasound signal
is modelled as a sampled continuous-time autoregressive moving
average (CARMA) process. Two new algorithms are proposed for
the estimation of continuous-time domain parameters of a CARMA
process from its sampled data. The first is a maximum-likelihood
estimator based on exponential splines suitable for any sampling
condition, even when some aliasing is present. Simulations results
show that this estimator outperforms current estimation methods
based on polynomial splines, which are appropriate only in high
sampling frequency conditions.
Additionally, another estimator is designed to provide good estimates with a lower computational cost when high sampling frequency is employed: this is usually the case in ultrasound images.
Although not maximum-likelihood, this estimator is still based on
exponential splines and can be used to compute the CARMA parameters of the ultrasound signals composing an image. CARMA
parameters are considered as new ultrasonic features and their ability in tissue representation are analysed on a phantom groundtruth and compared to traditional ARMA parameters, often used
to model the discrete ultrasound echo. The results of this study
show that CARMA parameters are able to capture the changes in
scatterers concentration better than traditional ARMA parameters.
In the second stage the main contributions of this thesis are (a)
the proposal of a new model for ultrasound signal, (b) the design
of a CARMA estimator providing correct estimates also in aliasing
conditions and (c) the definition of new features correlated to the
scatterers concentration in ultrasound images.
The third stage of this work focuses on the application of the
I NTRODUCTION
3
studied techniques to the improvement of the current prostate biopsy protocol. A large ground truth database of imaged biopsy cores
and the corresponding histo-pathological analysis are collected to
validate the accuracy of processing and classification techniques.
Improved statistical and textural features are used here in a specific learning scheme, able to use both completely known data and
uncertain biopsy cores, where pathology covers just a percentage
of the sample.
The developed algorithm for real-time computer-aided biopsy
(rtCAB) includes feature extraction and learning techniques and
is implemented exploiting CUDA™ parallel processing in order to
enable real-time support to physicians during biopsy sessions. This
technique is able to reduce the biopsy core number from 8-12 to 7,
preserving the same diagnostic power of a standard protocol and
working in real-time.
The structure of this thesis is formed by 5 chapters. The first
chapter develops the background resuming the basic information
about ultrasound imaging and tissue characterization. Chapter
2 describes the proposed characterization scheme with introduction of deconvolution as pre-processing step. Chapter 3 contains
the description of the new continuous-time domain model for ultrasound signals; the derivation of CARMA parameters estimators
and the application of CARMA parameters to the characterization
of ultrasound scatterers concentration. The fourth chapter reports
the realization of the ground truth and the developed learning algorithm for the improvement of the current prostate biopsy protocol. The conclusive chapter summarizes the main contributions
of this thesis and provides some clues about future research in ultrasonic tissue characterization.
One
B ACKGROUND
The essential is invisible to the eyes.
Antoine de Saint-Exupéry (1900-1944)
T
HE main goal of biomedical image analysis is monitoring inner
physiological systems to early diagnose potential diseases. A
pathological process affecting the investigated system reflects on
measured signals that, because of the disease, present characteristics different from the corresponding normal pattern. The essence of image analysis is to capture these characteristics, which
are often not visible to human eye, and to describe the state of the
imaged physiological system. Image modeling can also serve this
purpose, providing some mathematical descriptors of the measured signals that result to be modified by the state of the imaged
tissues and, thus, useful for their characterization.
The interpretation of biomedical images, traditionally performed by physicians, can be supported by automatic image analysis
tools, making a quantitative study of physiological systems possible. This process is called Computer-aided Diagnosis (CAD) [1]
and it is aimed at improving the reliability of physician judgement
by reducing subjectivity, inter- and intra-observer variability and
potential human errors due to fatigue, boredom or environmental
factors.
CAD tools have always been thought as a support and not a replacement of medical specialists, who integrate their experience
and intuition with objective clinical signs to realize the final diagnosis. Biomedical image analysis and enhancement techniques
constitute the core of CAD systems and can be crucial in healthcare quality.
Although the main target of medical image analysis and improvement techniques is early detection of potential diseases, these
5
6
One. B ACKGROUND
post-processing steps are also important to boost the performance
of cost-effective imaging devices. In fact the software enhancement of medical imagery is necessary both to compensate for quality reduction in video acquisitions of portable devices and to answer to the recent needs of world healthcare market, interested
in producing affordable medical equipments for developing countries.
Among all imaging modalities, ultrasound presents some interesting advantages: ultrasound devices are cheap, they provide
real-time feedbacks, they don’t use ionizing radiations, they are a
well established technology present in all health institutions. On
the other hand, ultrasound video sequences are affected by low
resolution and are characterised by speckle noise, preventing the
employment of traditional image processing algorithms to perform
computer-based analysis and enhancement.
Indeed, while some CAD tools are already clinically used in
more expensive imaging modalities, like CT and MRI, the realization of effective image analysis techniques in ultrasound, matching the clinical standards and possibly real-time, is still a challenge.
1.1 Diagnostic ultrasound
Ultrasound echo imaging [2] is the fastest and simplest medical
imaging technique, based on recording the reflected echo of ultrasound pulses. The main information provided by ultrasonography
is the echo intensity value associated to some spatial coordinates.
Unlike CT and MRI, ultrasound does not provide directly any
information about some tissue parameters, in fact the echo intensity depends in a complex way both on the imaged tissue and on
the acquisition system. For this reason ultrasound are used with
the only purpose of imaging the anatomy, or guiding some clinical
operations, but not to characterize tissues.
Nevertheless the echo intensity is just a small part of the information the reflected pulse can provide. The raw radio-frequency
(RF) signal returned by reflection is neglected by the traditional Bscan imaging modality, which directly computes its envelope (figure 1.1). Nevertheless the study of this signal can lead to the definition of significant features of the imaged tissues, as well as provide
a way to reduce system dependent effects. The RF signal is, thus,
the key to perform ultrasonic tissue characterization.
The reflection of an ultrasound pulse generates a single data
1.1. D IAGNOSTIC ULTRASOUND
7
Figure 1.1: Exemple of an RF signal and its envelope shown as a scan line in a
phantom ultrasound image.
line in ultrasound image along the propagation direction, called
axial direction. The resolution along this direction depends on
both the length of the ultrasonic pulse and its wavelength, still it
can only be increased within strict limits, since ultrasound absorption in tissues is linearly proportional to frequency. The poor quality of ultrasound image is due to the intrinsic resolution limits of
this technology.
The lateral resolution depends on the geometry of the probe
generating the ultrasound beam. The beam divergence angle is
smaller for larger diameters of the probe, making small transducer
suitable for near focusing and large probes suitable for grater depth
focusing. A curved probe surface, an acoustic lens or a delayed
stimulation of the transducer elements gives the possibility to tune
the beam divergence and radial direction.
The geometrical arrangement of several acquisition lines along
different directions produces 2D and 3D ultrasound data. For this
reason most algorithms in ultrasound imaging are developed as
1D processing techniques, and then applied to each image or volume scan line.
The fan scan is one of the most used ultrasound rays configurations and it is advantageous when the organ to be investigated is
protected by structures unpenetrable by ultrasound and requires
a small interface area of the transducer. A fan scan transducer is
used, for example, to image the prostate by trans-rectal ultrasound
(TRUS). The main disadvantage of the fan scan probe is the reduction of lateral resolution along the axial depth. Furthermore the
echo data are recorded in polar (r, θ)-coordinates and need a conversion into Cartesian coordinates and interpolation before dis-
8
One. B ACKGROUND
playing on monitors.
From the point of view of image analysis, the fan scan configuration does not constitute an additional difficulty, because the line
signals are processed directly in the (r, θ)-space, where they appear as a standard rectangular image. The (r, θ)-space is always
the reference, unless the processing algorithms require the knowledge of the real anatomical shapes of the imaged scene or in case
the processing results are forwarded to physicians as a map on the
anatomical image.
The three phenomena contributing to the realization of an ultrasound echo are (a) the pure reflection from a smooth surface, (b)
the diffuse reflection from a rough surface and (c) the refraction,
this latter causing artefacts during image formation which considers possible only line transmission. A fourth factor characterizing
the echo is (d) the diffusion caused by interaction with many rough
surfaces smaller than the ultrasound wavelength. In this case the
pulse is diffused in all directions, determining a random component of the ultrasound signal, called scattering. The presence of
scattering in the recorded echo makes the signal non-stationary.
Nevertheless the RF signal can be considered stationary in a
homogeneous area, where it is characterised by a slowly varying
amplitude and a stationary random component. The RF signal
can, thus, be decomposed into its local mean value, representing
the echogenicity of the area, and a variable component representing the speckle noise, a complex type of disturbance usually modelled as multiplicative.
Although the speckle noise is dependent on the acquisition system properties, it also contains some information about the imaged area, intuitively linked to the texture appearance in ultrasound
images. While image processing, aimed at improving the visual
quality of ultrasonographies, are focused on suppressing speckle
noise, image analysis techniques are often interested in capturing as much information as possible from this random component
which can be useful for tissue characterization.
One more important characteristic of the RF signal is the frequency shift: as signal attenuation increases with frequency, during propagation the high frequency components of the echo signal
are more attenuated than the low frequency ones, globally reducing the value of the central frequency of the received signal. The
frequency shift may also deteriorate the radial resolution, but the
estimation of these changes are useful to study the attenuation of
different types of materials and thus perform tissue characteriza-
1.1. D IAGNOSTIC ULTRASOUND
9
Figure 1.2: Standard signal processing chain in medical ultrasound acquisition
systems. The signal recorded by the probe is the raw RF echo signal, while the final
output displayed on monitor is the envelope of the base-band signal.
tion. The knowledge about the frequency shift can also be used
to compensate the attenuation independently for each ray during
time-gain compensation (TGC) procedures and restore the original echo spectrum.
Since the dynamic range of the RF signal is very large because
of the different amplitude of strong reflection components and
weak scattering, a logarithmic transform is employed to reduce
echo dynamics before digitalization. The final digital raw signal
is obtained by sampling the recorded signal, amplified and processed by TGC and logarithmic transform. The sampling frequency
is chosen in the range of 20-60 MHz and the number of bits used
for quantization is between 16 to 24.
In medical ultrasound imaging the standard signal processing
chain performs a demodulation of the RF signal by Hilbert Transform, providing the in phase-quadrature signal (IQ), and computes
the envelope as the module of the IQ analytic signal. Dealing with
the envelope allows a decimation of the signal according to the reduction of its band and provides a better visualization of the ultrasound image. A scheme of the standard processing chain in ultrasound imaging is shown in figure 1.2.
The base-band signal or B-mode signal is the standard output
of an ultrasound acquisition system, while the original RF signal
is the unprocessed information acquired by the probe. The basic
principle of most tissue characterization techniques in ultrasound
is extracting as much information as possible from both the baseband signal and the RF signal.
10
One. B ACKGROUND
(a)
(b)
Figure 1.3: Rectangular ROI segmentation inside the prostate region is illustrated
in figure 1.3(b), while figure 1.3(b) shows irregularly shaped ROIs.
1.2 Tissue characterization
The final purpose of biomedical image analysis is tissue characterization. In practice the analysis consists in segmenting the image
into several regions of interest (ROIs), computing some characteristics, features or measures related to each ROI and classifying them
into one of few known categories. Any category represents a particular state of the tissue outlined in a ROI; the simplest case is binary classification, aimed at discriminating a normal and pathological tissue, neglecting the presence of benignant abnormalities.
ROIs can be defined as irregular shaped areas in the image, following some anatomical characteristics of the imaged tissues; alternatively they can be chosen as rectangular regions covering the
whole image, with or without overlapping (figure 1.3). While rectangular tessellation can be easily performed by computers, segmentation of ultrasound images into irregular regions requires either a manual procedure or a complex automatic algorithm subject to high degree human supervision.
Ultrasound image segmentation presents some characteristic
artefacts that make the segmentation task complicated, and almost impossible by means of standard image segmentation techniques. In fact, while in other areas of medical imaging (CT, MRI)
the application of standard image processing methods is sufficient
to obtain good segmentations, in ultrasound satisfactory results
can only be achieved by specific methods. A review of segmentation techniques in medical imaging is provided by [3] and [4].
The ROIs consist of a two-dimensional arrangement of recorded
RF signals, which are used to compute some features useful to clas-
1.2. T ISSUE CHARACTERIZATION
11
sify the state of the outlined tissue. Let us consider the extraction of N different features from each ROI: they can then be represented as a point x = x1 , · · · , x N f in an N f -dimensional feature
space. ROIs with similar characteristics will form clusters in the
feature space.
The first step in tissue characterization is the selection of significant features which are able to capture the discriminant characteristics of normal and pathological tissues, ideally providing a
representation of different tissues in the feature space as separate
clusters. In the end, a pattern classification step in image analysis
is aimed at defining the optimal decision boundary in the feature
space to separate normal and pathological tissues.
The construction of the optimal decision boundary for tissue
characterization is possible thanks to supervised pattern classification techniques based on an accurate medical ground truth. The
collection of a large number of ultrasound images from biological
tissues and their corresponding histological examination is essential to create a ground truth, which is used in the design of a classification algorithm and its validation.
Ground truth data are divided into a training set and a testing set. The training set consists of a set of ROI-label pairs, already classified by a chosen gold standard like histological analysis, which are used to perform the learning of a classification algorithm. The testing set instead is used to validate the algorithm.
Test ROIs labels are used only to compute the performance criteria and do not provide any further information to the classifier. In
some cases also a validation set is extracted from the ground truth
in order to tune some parameters of the classification algorithm.
The main problem in ultrasonic tissue characterization is the
absence of a public standard database, also for very common and
highly diagnosed diseases, as prostate carcinoma or breast tumour.
One fundamental contribution of this research is the realization
of a large ground truth database for prostate cancer, collected in
collaboration with the department of Urology of the University of
Bologna.
In ultrasound imaging, system effects and tissue micro structures affect in a complex way the signal and the speckle noise, thus
any conclusion about tissue characterization is valid under the
same conditions and equipment settings. Nevertheless algorithms
designed to make the signal less dependent on acquisition system
adjustments can be employed to provide more efficient tissue characterization.
12
One. B ACKGROUND
Reducing system effects
Ultrasound imaging acquisition systems can affect and deteriorate
tissue response. In order to reduce these undesirable effects and,
hopefully, obtain a higher quality and more significant image, several processing techniques based on image formation linear models were developed.
The most common two-dimensional linear image formation
model is defined through a convolution operator:
x
x[n, m]
=
=
H(w) + µ
P P
k l h[n − k, m − l]w[k, l] + µ[n, m]
(1.1)
where x is the image, w is the tissue response to be recovered, H
is the linear operator accounting for system effects which defines
the point spread function (PSF) h, and µ is some additive term
accounting for both experimental and model noises. The inverse
problem of recovering w is also named as deconvolution [5]. In
medical ultrasound imaging the PSF is usually not known and is
estimated together with the true image w. In this case the problem is referred to as blind deconvolution [6].
The problem with model (1.1) is that the presence of tissues between the transducer and the target changes locally the PSF which
cannot be considered shift-invariant in the whole in vivo image.
The presence of phase aberration and dispersive attenuation in ultrasound imaging lead necessarily to a spatially variant model for
the PSF, which requires a non stationary deconvolution.
The most used technique to deal with shift-variant convolution
kernel is to segment the image in several areas where the signal
can be considered stationary and the local PSF shift-invariant. An
alternative to this approach is the employment of adaptive algorithms based on models able to change their characteristics along
with the signal.
Furthermore, although the image formation model is expressed
in 2D, ultrasound imaging deconvolution algorithms are often developed in 1D and separately applied along the lateral and axial
direction. Anyway it is common practice to perform deconvolution along the axial direction only, as the PSF is wider. In this case
the convolutional model becomes:
X
(1.2)
x[n] = h[n − k]w[k] + µ[n]
k
The separability of the PSF allows one to reduce the problem of
a 2D deconvolution to a sequence of 1D deconvolution problems,
1.2. T ISSUE CHARACTERIZATION
13
applied along each direction independently. Specifically, it is assumed that the 2D convolution kernel can be expressed as product
of two one dimensional convolution kernels along the two directions: the axial PSF includes the blurring effects due to the finite
bandwidth of the transducer and the dispersive attenuation of the
acoustical pressure pulse in tissue, the lateral PSF represents the
convolutional components of lateral blurring due to the complex
beam pattern.
One of the classical methods of blind deconvolution used in
medical ultrasound imaging originates from the theory of system
identification and is based on statistical modeling. An efficient
statistical method for deconvolution assumes the RF signal to be
an autoregressive moving average (ARMA) process and estimates
model parameters to recover the tissue response. Hence, in this
interpretation, the RF-signal x[n] is considered to be the output
sequence of a causal filter h[n], while the tissue response or reflectivity function, w[n], is thereby considered to be an input driving
sequence. The latter is commonly assumed to be a realization of
an independent and identically distributed process of zero mean
and fixed variance.
ARMA models present highly flexible modeling capability, because they integrate properties of AR models, useful for modeling
random signals that possess peaky power spectral densities, and
properties of MA models, more appropriate for signals that have
broad peaks or sharp nulls in their spectra.
For the case of ultrasound imaging the choice of an ARMA model
is well justified. This is because the spectra of the RF lines are normally wide-peaked, and, since the requirement of stationarity imposes limitations on the range of sampling intervals, stationary RF
signals are usually short in duration. ARMA modeling is appropriate in this context, since it provides a well-posed estimation problem with fewer parameters than observations and is characterized
by excellent modeling capability for relatively low orders and for
short data segments.
While deconvolution techniques are generally used to improve
the visual quality of ultrasound image, in this work it is introduced
as a preprocessing step in the standard tissue characterization procedure, in order to enhance the diagnostic significance of ultrasound features. A scheme of the proposed tissue characterization
protocol is shown in figure 1.4.
14
One. B ACKGROUND
Figure 1.4: Proposed characterization procedure with the addition of deconvolution to improve diagnostic performance.
Feature extraction and selection
A large number of features were designed to extract interesting information from the ultrasound RF signal. The B-mode signal envelope itself can be considered as a feature containing all the visual
information available to the physician.
The features derived from the RF signal developed for tissue
characterization usually belong to one of these three classes:
• Spectral: these features gather information about the spectral behaviour of the RF signal. This may include the study
of the frequency shift and the related attenuation; the spectral decomposition of the coherent and random component
of the signal; multi-resolution analysis by wavelets or the estimation of any mathematical model for the spectrum.
• Statistical: the RF image can be modeled as random samples
from some statistical distribution. The parameters of such a
distribution can be estimated locally and provide statistical
features.
• Textural: these characteristics give information about the arrangement of grey levels and are computed from the B-mode
signal. Textural features represent the closest attributes to
what is processed by human brain when analysing an image.
The number of features can increase rapidly as new models for
the RF signal or new attributes are proposed. The rise in dimensionality can only apparently make the discrimination task easier.
In practice not all the features extracted are interesting for tissue
characterization. Furthermore high dimensionality causes more
complex classification models which tend more likely to overfit the
1.2. T ISSUE CHARACTERIZATION
15
training data. For this reason it is essential to select few significant
and non-redundant features.
Feature selection algorithms [7] focus on the realization of the
optimal feature subset for a certain problem, according to a dataset
containing feature values for any sample and the corresponding
category. In the context of image feature selection, the dataset is
constituted by feature values from any ROI and the information
about the pathological state of the tissue imaged by the region. The
selection approaches can be roughly divided into three categories:
filters, wrappers and hybrid methods.
Filters sort features according to a score measure that summarizes the importance of a feature with respect to the state of
the tissue. Typically, the ranking criterion can measure either distance, information, dependency, or consistency between features
extracted from each ROI and the respective class. Filters measures
are independent on any classification algorithm.
On the contrary, wrappers sort attributes from a feature set according to the performance of a classifier trained on that feature
set; they are thus classifier-dependent and require a larger computational cost.
Hybrid methods take advantage both of filter and wrapper methods. The filter measure is used to decide the best subset for a given
cardinality, while the wrapper mining algorithm selects the final
best subset among the best subsets across different cardinalities.
In general a complete search procedure is never adopted in
feature selection because of its huge computational cost; on the
contrary, although they cannot guarantee that the actual best feature subset is selected, forward and backward sequential search
or random search techniques are preferred because of their lower
computational cost.
The above methods are obviously supervised since they take
advantage of the information about tissue state, but also unsupervised methods such as principal component analysis (PCA) can be
useful, above all in order to reduce redundancy in subsets where
features have similar informative value. An unsupervised method
can thus be used as a first step to discard redundant attributes and
a second step of supervised feature selection can subsequently be
applied to further reduce the feature set, eliminating attributes irrelevant for the particular classification problem.
A complete feature selection scheme can be tested by introducing artificially designed irrelevant and redundant attributes in
a benchmark feature set; an efficient feature selection algorithm
16
One. B ACKGROUND
must be able to detect and discard these test attributes.
ROI classification
The goal of image tissue characterization is employing the feature
set to build a classifier able to discriminate between different categories of imaged tissues. Among the several approaches to classification existing in literature [8], non linear discriminant analysis
through kernels is used in this thesis to take advantage of the non
linear mapping of the feature space and ease the classification task
for real data.
In supervised machine learning the data set (xk , c k ) composed
of ROI feature vector and label pairs is used to determine the predictive model representing the relationship c = f (x) between the
input variables x and the output variable c. It is the knowledge
about class labels c k to make the learning problem supervised.
The found predictive model is then used to separate in the feature space and normal and pathological ROIs. In this case the classification problem is binary and the label index c k can only assume
two values, for example 1 for benignant tissue and 2 for pathological tissue.
A linear model f (x) employs a predictive function represented
by an hyperplane in the feature space, whose parameters w and
b are determined using classical criteria such as least squares or
maximum likelihood:
f (x) =< w, x > +b
(1.3)
In this case the predictive model is also called discriminant function. Traditional linear methods give the same importance to all
samples in the dataset and this can easily provide non efficient linear separation, because of the excessive influence due to samples
very far from the separation boundary.
A way to overcome this problem is focusing on those training
examples which lie close to the class they do not belong and, thus,
close to the boundary. These samples are called support vectors
and in this approach they are the only data that explicitly define
the model. A linear model in terms of support vector, sv, is expressed as in the following:
f (x) =
X
i∈sv
αi c i < xi , x > +b
(1.4)
1.2. T ISSUE CHARACTERIZATION
17
where the summation includes only the training examples xi that
are support vectors, and αi are coefficients determined as Lagrange
multipliers in the optimization procedure.
An advantage of this approach is that the classifier concentrates directly on examples sv that are difficult to classify and its
computation scales with the number of support vectors rather than
the dimension of the space which can be very large. In addition
this approach is shown to balance training classification error and
model complexity, thereby avoiding overfitting, a situation in which
the model is too finely tune to the training examples and fails to
perform well on new data.
When real ultrasound data are analysed, linear separation is
generally not possible in the feature space, even if a large number
of features is employed. In this condition it is important to exploit
techniques providing a simple way to obtain a non linear model
from any linear model based on inner products. Even classical
techniques can be turned easily into flexible non linear techniques
via the so called kernel trick.
According to this technique, a non linear function φ mapping
the data into a different higher dimensional feature space can be
applied. After this transformation separability will be enhanced
thanks to the higher dimension of the mapped space, which can
also be infinite. The kernel trick recognises that the classification
in this new space can be obtained without actually performing the
transformation, because the transformation φ applied to samples
in predictive function only appears in the form of an inner product. Therefore non linear approaches can be performed efficiently
by using kernel functions K(., .) which act like dot products in the
re-mapped feature space [9].
It is never actually necessary to compute φ or to define it explicitly. Instead it is sufficient to define the kernel function K(., .),
and it can be shown that any symmetric positive semidefinite function suffices.
For each ROI sample, x, the kernel approach provides its projection onto a vector in F , expressed as an optimized linear combination of the inner products of the mapped sample and a sequence
of mapped support vectors:
X
X
ai < φ(xi ), φ(x) >=
ai K(xi , x).
(1.5)
i∈sv
i∈sv
The combination weights ai are tuned to maximize the inter-class
variance and minimize the intra-class variance of data in the mapped
18
One. B ACKGROUND
space. Intuitively, the effect of the kernel is to measure the similarity between a test vector x and each of the support vectors. In
fact vectors belonging to one of the classes are presumably most
similar to the support vectors belonging to that class, hence these
similarity values provide us with the information to build the predictive function.
In non linear data transformations an appropriate kernel for
representing the similarity of a mapped sample to the support vectors is the radial basis function (RBF) kernel:
Ã
!
−kxi − xk2
K(xi , x) = exp
.
(1.6)
2σ2r b f
The kernel variance σ2r b f is tuned by cross-validation on a validation data set. In particular, a part of the dataset is reserved to perform several tests varying the value of σ2r b f , in order to find the
optimal value for this parameter, which provides the best classification performances.
Cross-validation is one of the techniques for statistical resampling used to evaluate performance and improve robustness of machine learning models and to estimate statistical significance levels. These techniques are as important as the predictive model itself.
After non linear transformation, a linear classifier can be used
in the mapped space to find the discrimination hyperplane defined by the following predictive function:
X
f (x) =
ai K(xi , x) + b.
(1.7)
i∈sv
Linear classifiers are suitable to classify kernel transformed
data, which are characterized by a better separation in the new
feature space. For example, the Fisher Linear Discriminant (FLD)
approach tunes the parameters of the linear model in (1.3), w and
b, by maximizing the Mahalanobis distance and in such a way that
different clusters have the same distance from the hyperplane. The
Mahalanobis distance J between two clusters belonging to different classes is expressed through the intra class mean value µi and
covariance matrix Σi of class i :
J = (µ1 − µ2 )T (Σ1 + Σ2 )−1 (µ1 − µ2 )
(1.8)
The performance indicators computed on the testing set offer a way to assess the diagnostic validity of a CAD scheme. In a
1.3. A PPLICATION : PROSTATE TRUS
19
binary tissue characterization problem, where the two categories
benignant/malignant tissue are available, it is common to define
True Positives (T P ) the pathological samples classified as malignant, False Positives (F P ) the misclassified normal regions, True
Negatives (T N ) the normal samples classified as benignant, False
Negatives (F N ) the misclassified pathological regions.
According to these definitions, besides accuracy (ACC ), the standard performance criteria are sensitivity (SE ) and specificity (SP ),
giving information about the ability in recognizing pathological
and normal tissue respectively. The prediction capability of the
method is also expressed by the probability that a malignant classified sample is actually pathological. This information is provided
by the positive predictive value (P PV ). The performance indicators are resumed in the following:
ACC
SE
SP
P PV
=
=
=
=
T P +T N
T P +T N+F P +F N
TP
T P +F N
TN
T N+F P
TP
T P +F P
(1.9)
The above definitions of performance criteria consider a region as a sample of the dataset, but according to the focused scale
of the problem, also pixel-based, core-based and patient-based
performances can be defined, considering as sample of the dataset
respectively the pixel, the core or the patient.
The variation of the threshold in a linear classification method
allows the complete description of the classifier performances,
which can be visualized in a receiver operating characteristic (ROC)
curve plotting the (SE , 1−SP ) graph. The area under the ROC curve
is also specific of the method diagnostic accuracy which increases
as the area approaches the unit value.
1.3 Application: prostate TRUS
One of the applications where Computer-Aided Detection (CAD)
tools would be extremely valuable is the recognition of prostate
cancer in trans-rectal ultrasound (TRUS) images. Some information about the physiology and pathology of prostate gland will be
presented in the following paragraphs.
The main function of the prostate is to store and secrete a clear,
slightly alkaline (pH 7.29) fluid that constitutes 10 − 30% of the volume of the seminal fluid that, along with spermatozoa, consti-
20
One. B ACKGROUND
Figure 1.5: Schematic representation of the prostate gland.
tutes semen. The rest of the seminal fluid is produced by the two
seminal vesicles.
A healthy adult human prostate is a chestnut shaped gland enveloped in a fibrous capsule. Its base is attached below the urinary
bladder neck and the apex is fixed to the urogenital diaphragm.
It borders on the posterior side with the rectum and on the anterior side with the fibromuscolar stroma connected to the pubis
through the puboprostatic ligaments. On the superior posterior
side it is attached to the seminal vesicles, a pair of simple tubular
glands that secrete a significant proportion of the fluid that ultimately becomes semen. The excretory ducts of seminal vesicles
open into the vas deferens, as they enter the prostate gland, and
they are lined with the epithelium of the transition zone. Within
the prostate, the urethra coming from the bladder is called the prostatic urethra and merges with the two ejaculatory ducts. The prostate is finally sheathed in the muscles of the pelvic floor, which
contract during the ejaculatory process.
The prostate gland is divided in three different zones: the transition zone (TZ), the central zone (CZ), and the peripheral zone
(PZ) (figure 1.5). The transition zone surrounds the urethra and extends from the ejaculatory ducts proximally. The peripheral zone
encompasses the urethra from the base to the apex. The central
zone is composed of tissue immediately surrounding the ejaculatory ducts and it expands inferiorly. The significance of this architecture is based upon the relationship of these three zones to prostatic disease [10]. In the young males the peripheral zone comprises 75% of prostate volume, the transition zone 20% and the
1.3. A PPLICATION : PROSTATE TRUS
21
central zone 5%. These ratios change and after the age of 40 years
the transition zone may enlarge ad occupy most of the gland as
benign hyperplasia is almost inevitable.
The prostate gland is the male organ most often smitten by either benign or malignant lesions and the second leading cause of
cancer death for men [11]. Prostate diseases can be distinguished
in three main categories: prostatitis, benignant hyperplasia and
prostate cancer.
Prostatitis is an inflammation of the prostate gland. It is a benign pathology and in its acute case it is mainly treated with antibiotics.
Benign hyperplasia (BPH) consists in a prostate enlargement.
It is fairly common among ageing men and it occurs mainly in the
transition zone. BPH can be treated with medication, a minimally
invasive procedure or, in extreme cases, with a surgery procedure
that removes the prostate.
Prostate cancer, often referred as prostate carcinoma (or adenocarcinoma), is a malignant pathology and, as all cancers, it is
characterized by an abnormal and uncontrolled cells mutation and
replication. Like BPH, prostate cancer might cause pain, difficulty
in urinating, problems during sexual intercourse, erectile dysfunction and other symptoms. If not detected on early stages and in
presence of more aggressive forms, the disease can advance to
stages characterized by local invasion of the surrounding tissues
(seminal vesicles, bones, rectum), usually resulting in lethality.
The heterogeneous and multifocal nature of prostate cancer lesions poses significant difficulties in its detection. With regards to
heterogeneity, prostate cancer tissue typically reveals a juxtaposition on benign cells, preneoplastic lesions and neoplastic lesions
of varying severity [10].
The current clinical procedure to detect cancer is based on a
combination of different diagnostic tools, because none of these
tools is accurate enough to be used alone. Digital Rectal Examination (DRE), Prostate-Specific Antigen (PSA) evaluation, TRUS
image analysis and biopsy are all part of the medical procedure
for prostate analysis. Each of these tools presents important limitations that make accurate prostate cancer detection still an unsolved problem.
Historically, DRE has been the principal method of prostate
analysis, but it is accurate only at detecting large and superficial
lesions and it is strongly operator dependent.
The first step in public screening for prostate cancer is the mea-
22
One. B ACKGROUND
Figure 1.6: The picture on the left gives a schematic representation of prostate
axial section with orientation inverted respect to its real anatomy. On the right a
TRUS image shows the prostate gland in the same orientation.
surement of the blood PSA level, a glycoprotein produced almost
exclusively in the epithelium of the prostate gland. Unfortunately,
PSA is specific for prostate but not for cancer, since other factors
such as BPH, prostate infection, urethral instrumentation and irritation can cause an increase in the PSA value.
At TRUS, the normal prostate gland has a homogeneous, uniform echo pattern. The appearance of carcinoma on ultrasound
is variable and in its early stages a tumour can appear anechoic,
hypoechoic or isoechoic with respect to the surrounding normal
tissues. Potential hypoechoic regions could also include BPH or
even normal biological structure, thus the specificity of the TRUS
images visual inspection is low.
Standard TRUS images show the axial or the lateral section of
the prostatic gland. The gland anatomical shape is inverted in the
image due to the acquisition setting, as shown in Fig. 1.6, where
the superior hyperechoic boundary marks the interface between
the rectum and the peripheral zone of the prostate. The transition
zone is visible as a darker region below the peripheral zone. The
gland contours are generally well defined and visible in the images.
The histopathological analysis of biopsy samples is the standard for cancer detection confirmation.
The dominant and most reliable method for prostate carcinoma
diagnosis and aggressiveness assessment, in research as well as
in clinical procedure, is the Gleason grading [12]. This method is
based entirely on the histological evaluation of carcinomas cells
arrangement patterns in contrast agents stained prostatic tissue
sections. Specifically, the method is a categorization of glandular
differentiation and pattern of growth of the tumour, at relatively
low magnification (×10 - ×40) in five basic grades patterns [12].
1.3. A PPLICATION : PROSTATE TRUS
23
Due to the heterogeneous nature of prostate cancer lesions,
histopathological inspection of prostate tissues often reveal an ensemble of benign lesions, prenoplastic lesions and neoplastic lesions with different aggressiveness. To take into account this heterogeneity, the five basic grade patterns are used to generate an
histological score, which can range from 2 to 10, by adding the primary grade pattern and the secondary grade pattern. The primary
pattern is the one that is predominant in area, by simple visual inspection. The secondary is the second most common pattern.
The Gleason grading score is therefore and indicator of cancer stage: the higher the grade is, the more advanced the cancer
is. Typically, cancers with Gleason scores lower than 6 are considered well differentiated and associated with a good prognosis.
Those with a Gleason score of 8-10 have the worst prognosis and
the highest risk of recurrence. Gleason score is often combined
with PSA level and clinical stage for risk assessment. The primary
goal of staging is to distinguish patients with organ-confined, locally invasive, or metastatic diseases.
Ultrasounds are used for biopsy guidance, to enable the sampling of all relevant areas of the prostate by means of systematic
sampling protocols [13]. However, the main limitations in this procedure are due to the multifocal nature of cancer and to the sampling process.
In fact, uniquely among image-guided biopsies, prostate biopsy is not lesion-directed but rather based on a systematic sampling of those areas where cancer incidence is higher. As the disease is often multi-focal, different areas throughout the whole gland
are sampled. In particular, since most of the cancers arises in the
peripheral zone, most of the protocols aim to maximize its sampling. The motivation behind this, as discussed before, is the weak
diagnostic significance of imaging inspection techniques due to
the high variability of prostate cancer patterns.
Ultrasound imaging is thus used to guide and monitor the biopsy needles insertions in the different areas of the prostate gland.
Since metal, if compared to soft tissues, has an high acoustic impedance, biopsy needles appears as hyperechoic lines in TRUS scans,
as shown in figure 1.7.
The first landmark sampling technique was the sextant protocol reported in 1989 [14]. As originally described, six biopsies were
obtained in a parasagittal line drawn halfway between the lateral
border and midlines bilaterally, from the base, midgland and apex,
as shown on the left side of figure 1.7.
24
One. B ACKGROUND
Figure 1.7: The picture on the left gives a schematic representation of double
sextant biopsy protocol. On the right a TRUS image shows the insertion of biopsy
needle in the prostate.
Although sextant biopsy protocol was a major advance, with
a 20-25% of positive biopsy rate, with a wider experience it was
found also inaccurate, principally because it under-samples the
peripheral zone [13].
Modifications of the sextant protocol were introduced from the
mid 1990s onwards. For example, in the modified sextant biopsy
protocol a better sampling of the peripheral zone around the lateral margins is obtained moving laterally and angling anterolaterally the biopsy trajectories. This improved in some cases the detection rate from 80 to 89 % [13].
However, even the modified sextant protocol was found to miss
some tumours with time and many alternatives where explored.
Therefore, several extended protocols using more cores (8 or 10)
directed to the peripheral zone were introduced. Although the
available results show that the extended protocols can improve diagnosis accuracy it is still a matter of debate whether extended
protocols are substantially better than the modified sextant protocol.
As regards patients perception of this examination, it was documented that 55% of men report physical discomfort during biopsy. Moreover, this procedure is not completely safe, since it carries the risk of bleeding infection or even urosepsis [15].
Research to improve the efficiency of the standard prostate biopsy protocol is proceeding both from an image processing point of
view and from a purely clinical perspective. In fact, as discussed
before, the positive predictive value p p v 0 of this protocol, that is
the probability of sampling a pathological core given a pathologi-
1.3. A PPLICATION : PROSTATE TRUS
25
cal patient, is quite low, between 20 and 25%. As the probability
to detect the disease in a pathological patient increases with the
number of samples N0 , according to 1−(1−p p v 0 )N0 , the detection
rate of the protocol is about 89%. The inefficiency of this protocol
stems from the fact that most of the sampled tissue is benignant
and, thus, represents an unnecessary biopsy.
In order to improve the efficiency of the prostate biopsy protocol it is essential to increase the P PV , in such a way that the clinical procedure can preserve the same detection rate of the standard
protocol with a reduced number of cores.
Since the processing of TRUS images highlights important characteristics of the investigated tissues, tissue characterization techniques can be employed to increase the P PV , providing the radiologist with a malignancy map to perform lesion-directed biopsy.
In this new scenario the number of unnecessary biopsies would be
reduced, while preserving the same diagnostic power as the standard protocol.
Two
U LTRASOUND
IMAGE ANALYSIS
AND INTERPRETATION
I imagine this is because Nature wants to
ensure that the evils of wisdom shall not
spread further throughout mankind.
Desiderius Erasmus of Rotterdam (1466 - 1536)
C
OMPUTER- AIDED detection schemes are decision making support tools, useful to overcome limitations of problematic clinical procedures. The topic treated in this thesis is the study and the
realization of CAD techniques based on trans-rectal ultrasound images, which may be extremely important to support prostate cancer diagnosis.
First works on TRUS-based CAD schemes for detection of prostate cancer consisted on biopsy ground truth and analysis of rectangular regions around needle insertion points. The main characteristic of these studies is the use of textural features extracted from
TRUS images to discriminate different tissues. Basset et al. [16],
Huynen et al. [17] and Houston et al. [18] realized clinical studies
based on the extraction of first and second order statistics textural
parameters and used simple decision trees to perform classification.
After these works, the main trend shifted to a multi-feature approach, extracting features of different nature from TRUS images
and combining them to obtain higher classification performance.
Schmitz et al. [19] and Scheipers et al. [20] employed textural features and spectral parameters extracted from RF data, while
Feleppa et al. [21] introduced clinical data like PSA value and patient’s age in the feature vector besides spectral parameters.
The classifiers used for these works (self organizing Kohonen
map, neuro-fuzzy systems, Artificial Neural Network) are more
27
28
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
complex with respect to decision trees and define a new direction
in this research field. Furthermore in these works the ground truth
comes not from biopsy but from prostatectomy and histological
analysis of prostate slices.
The work of Mohamed and Salama[22] represents an exception
since the gold standard is based on radiology visual inspection.
The proposed scheme utilizes a pure textural feature vector with
a Support Vector Machine (SVM) classifier. In this case the values
of sensitivity, specificity and area under ROC curve are high, but
they are obtained on a small ground truth.
Similar performance was obtained in latest studies like the extension of Mohamed’s work [23] which includes spectral features
and the study performed by Han et al. [24], where morphological features and multi-resolution textural features are employed
and SVM is used as classifier. In the latter work notable values of
sensibility and specificity are reported, but the proposed method
was only tested on malignant images, thus there is no information
about its behaviour in completely healthy cases.
This chapter presents an effective approach to realize a CAD
scheme for prostate cancer detection, employing a multi-feature
nonlinear classification model and the predictive deconvolution
of the acquired RF signals to reduce system-dependent effects.
The mutual information of feature values and tissue pathological state is used to select, among several analysed ultrasonic characteristics, the features essential for tissue characterization.
A clinical study, performed on ground truth images from biopsy findings, provides a comparison of the classification model applied before and after deconvolution, showing in the latter case a
significant gain in accuracy and area under the ROC curve.
Here the aim is to investigate the ability of the features extracted
from deconvolved US images at discriminating pathological tissues. This issue is analysed in terms of performance comparison
of a nonlinear classification model trained on features extracted
from US images, with and without deconvolution preprocessing.
In addition, while rectangular ROIs are usually employed in
ultrasonic tissue characterization literature, we propose an algorithm for irregular shaped ROIs segmentation and perform analysis to test the usefulness of this segmentation for tissue characterization purposes, without including the deconvolution step.
A large number of features is extracted from US RF signal to
characterize different aspects of biological tissues. The combined
use of features of different nature results in more accurate tissue
2.1. D IMENSIONALITY REDUCTION
29
characterization, but requires a critical operation of feature selection to identify features highly correlated to the pathological state
of the tissue.
In our approach a hybrid feature selection algorithm based on
the mutual information of feature set and ground truth class is
used to prune unimportant features and achieve fast computation.
The ground truth used in this clinical study is based on biopsy
findings and histological analysis of some regions in the images
considered suspicious at visual inspection by expert radiologist.
The proposed CAD scheme was tested on a ground truth containing both benignant and malignant cases.
The analysis of the CAD scheme employing irregular shaped
ROIs shows that non linear classification ensures in this case a sensitivity always greater than 80%, but provides a limited specificity.
Furthermore, the addition of predictive deconvolution in the
proposed CAD scheme increases classification performances, providing a sensitivity of 90%, specificity of 93% and area under the
ROC curve of 95%.
2.1 Dimensionality reduction
Results of recent studies show that combining features extracted
from RF analysis of ultrasound signals and image-based texture
parameters results in more effective classification procedures [25].
Literature about tissue characterization in ultrasound analysis provides a large amount of features of different nature which can be
subdivided according to their contribution in highlighting some
properties of the tissue.
Parameters of statistical distributions, such as Nakagami model
[26], give information about scatterer density, regularity and amplitude.
Spectral features [27] describe fluctuations of physical properties as acoustic impedance, viscosity and elasticity resulting in
backscattering signals. Typical spectral parameters capture the
shifting of RF signal central frequency due to attenuation. Also
the wavelet coefficients of RF signal, their polynomial fitting [28]
and the coherent and diffuse components obtained by signal decomposition [29] belong to the spectral features group and provide
important properties to type tissues.
Features extraction from B-mode images aims mainly at the
detection of textural properties of speckle which represents the
30
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
macroscopic appearance of the scattering generated by tissue microstructures. Different kind of textural parameters are available
in literature. Haralick [30] and Unser features [31] are both based
on the gray levels distribution statistics, while Fractal features [32]
rely on modeling and computation of fractal dimension. In our approach, the characteristic skills of different features are combined
to define a feature set endowed with a high discriminating power
between healthy and cancerous tissues.
A complete feature set of all parameters mentioned before
would have a huge dimensionality of about 140 attributes and is
described in appendix A. For this reason, a first selection step is
performed keeping for each group of features only those correlated
to
the
ground truth class, and discarding the other ones.
The feature set is built from the ground truth dataset that will
be described in section 2.5, computing features on rectangular overlapping ROIs selected from the images. The feature data matrix is
represented as S = [f1 | · · · |fN f ], where fi is the i -th feature vector
with a value for each ROI, such that S is a matrix of size NRO I × N f :
number of ROI × number of features. The ground-truth class vector, containing the index representing the pathological state of the
tissue outlined by a ROI, is expressed by c of length NRO I .
By selecting feature through correlation analysis, the dimensionality is reduced, but synergies between different features are
saved. The defined feature set is constituted by 54 features, as
shown in table 2.1.
Table 2.1: Feature Set
Feature
Origin/Type
#
Wavelet Transform (WT)
Polynomial Fit of WT
Wavelet Decomposition
Central Frequency
Attenuation
B-mode
Nakagami
Statistic
Haralick
Unser
Fractal
RF/Spectral
RF/Spectral
RF/Spectral
RF/Spectral
RF/Spectral
B-mode/Envelope
RF/Statistic
RF/Statistic
B-mode/Textural
B-mode/Textural
B-mode/Textural
2
1
1
1
13
1
4
2
4
9
16
The importance of CAD in diagnostic ultrasound imaging depends on its ability to perform nearly real-time classification in
2.1. D IMENSIONALITY REDUCTION
31
order to give a second opinion to the physician. A large feature
set dimensionality prevents the use of ultrasound images automatic characterization as a diagnostic decision support tool, because feature computation would be too demanding in computational terms. For this reason a further classification-oriented Feature Selection (FS) step is essential.
In particular, in our approach a Mutual Information Hybrid FS
(MIHFS) algorithm is used to rank and prune the whole feature
set. Ranking measures based on distance and dependency [33]
were also tested and, although most of them succeed in discarding irrelevant features, only information based measures are able
to recognize and discard redundant features. In the proposed FS
technique the chosen classifier-independent measure is the minRedundant Max-Relevance (mRMR) criterion proposed by Peng et
al. [34]. The mRMR measure is based on mutual information between the current feature set and ROI class, corrected by the averaged mutual information between features in the feature set. Maximizing this measure allows to define a feature set with maximum
relevance D(S, c) and minimum redundancy R(S), by maximizing
the mRMR measure:
1 X
I (fi , c)
(2.1)
max D(S, c) D(S, c) =
S
|S| fi ∈S
min R(S) R(S) =
S
1
|S|2
X
I (fi , f j )
(2.2)
fi ,f j ∈S
S = argmax Φ(D, R) Φ(D, R) = D − R
(2.3)
where I (X , Y ) is the mutual information of two discrete random
variables X and Y , representing two features or a feature and the
ROI class:
¶
µ
X X
p(x, y)
(2.4)
I (X , Y ) =
p(x, y) log
p 1 (x)p 2 (y)
x∈X y ∈Y
The joint and marginal probability density functions, p(x, y), p 1 (x)
and p 2 (y), can be estimated through the ROI-based histograms of
feature values.
MIHFS is an algorithm in two steps. The filter step consists
in ranking all features according to the mRMR measure, shown in
(2.3), following a sequential forward selection as search technique.
MIHFS wrapper step exploits a Fisher Linear Discriminant (FLD)
[35] as classifier and evaluates the mining performance of the ranked
32
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
feature set at increasing set size. For different cardinalities, the
subset which maximizes the mRMR measure is selected and the
performances of FLD trained on this subset are computed. The
best cardinality and consequently the best subset is chosen as the
one performing the minimum FLD misclassification error. Typically the best cardinality is smaller than the maximum number of
features because the classifier over-fits on training data given by a
large number of features.
Since hybrid FS algorithms require both a time-consuming
search for ranking and several training iterations, FLD was selected
to speed up selection experiments. For the selected subset the MIHFS algorithm produces a ranked list of features which highlights
their predictive skill. A previous analysis performed on ground
truth samples has shown the reliability of MIHFS technique in discarding both irrelevant and redundant features. Obviously, any
kind of preprocessing of ultrasound images (such as deconvolution) could modify the features ranking and, in general, can improve or worsen the discriminating power of a feature.
2.2 Feature-based segmentation
CAD techniques relaying on imaging are typically based on the
segmentation of the given image into ROIs and this task is usually
performed by an expert radiologist. In order to make image classification independent of radiologists skills and give them an additional tool for making diagnosis, ROI segmentation may be possibly automated.
The use of rectangular and possibly overlapping ROIs makes
the automation of segmentation very easy. Nevertheless irregular shaped ROIs may be more interesting, since they can outline
anatomical regions or abnormal structure present in ultrasound
image. Segmentation into irregular shaped ROIs is less easy to automatise, but still some algorithms based on ultrasonic features
are able to automatically highlight significant zones.
The ROI segmentation procedure proposed in this work consists of two phases: feature calculation and image segmentation.
In the feature calculation phase the RF signal is used to compute
spectral, statistical and textural features over a rough selection of
the prostate gland that can be performed also by non trained personnel. These features are then used in the segmentation procedure, based on an unsupervised Bayesian learning technique, to
2.2. F EATURE - BASED SEGMENTATION
33
Figure 2.1: Schematic representation of the proposed procedure for computer
aided classification of TRUS images.
assign each image pixel to a class. The assignment is realized by
a K-means clustering followed by an iterative Bayesian regularization algorithm. Such algorithm imposes spatial smoothness constraints to the solution, minimizing the energy function derived by
modeling the signal features with a multivariate Gaussian and the
pixel segmentation class with a Markov random field (MRF). The
contiguous pixels forming a closed region are then identified as a
ROI.
Ultrasound image segmentation is strongly influenced by signal quality as well as by the tissue imaged. The characteristic artifacts of ultrasound images like attenuation, speckle, shadows and
signal dropouts, missing boundaries due to the orientation dependence of acquisition, make the segmentation task complicated.
While in other areas of medical imaging (CT, MRI) application of
general image processing methods is sufficient, in the case of ultrasound signals more complex and specialized methods are
needed in order to obtain satisfactory results.
In literature, there are a large number of papers describing segmentation procedures applied to ultrasound images. For a complete review of the state of the art of medical and ultrasound image
segmentation methods see [3] and [4].
Most of the techniques presented in literature concern prostate boundary segmentation, while only few works about ROI segmentation have been published. Moreover, most of the methods
depend on human interaction and their performances are user dependent. Therefore, to overcome the limitations of the existing
procedures, new procedures should ideally be: user independent,
independent on training images, focused on ROI segmentation.
According to the review on the published methods for prostate
segmentation and tissue characterization, relying on the paradigm
34
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
Figure 2.2: Schematic representation of ROI segmentation algorithm.
of fig.2.1 the procedure for ultrasound aided prostate cancer diagnosis is based on three steps: automatic segmentation of prostate
into ROIs, feature estimation on the segmented ROIs, characterization of ROIs through a classifier.
The automatic ROIs segmentation procedure, proposed here,
is composed of four sequential steps as summarized in fig.2.2: feature estimation, clustering, Markov random field based regularization and ROIs selection.
Since many boundary estimation algorithms are available in
literature, a rough estimation of prostate boundary performed by
automatic algorithm or done manually by non trained personnel
is considered as given. Feature estimation is then performed pixelwise only inside the prostate gland, by means of a sliding window.
The estimation is done following a multi feature approach, i.e. different classes of features are estimated starting from the radio frequency signal before time gain compensation and envelope detection. B-mode features are estimated after envelope detection
and log compression of the radio frequency signal. Among the
features proposed in literature we employed texture parameters
(Unser features [36]) and B-mode statistical modeling (Nakagami
fitting [26]).
The estimated features are then used to segment the image by
mean of a clustering procedure. Clustering techniques are unsupervised learning technique that easily allows the combination of
multiple description parameters and data mining in a high dimension space, with a feasible computational cost [35]. In this work
we adopted K-means algorithm, a based clustering method which
evaluate the similarity of pixels by means of their distance in feature space. More precisely, pixel portioning is obtained by mini-
2.2. F EATURE - BASED SEGMENTATION
35
mizing the following energy function:
W (C ) =
K
X
1 X
d(Xi , y k )
2 k=1 C (i)=k
(2.5)
where C represents the clustering of all the data points i in K predetermined number of clusters, C i = k stands for the assignment
of the i -th data point to the k-th cluster with center y k , and d(., .) is
the distance. Xi = {X i,1 , · · · , X i,N f } represents the N f -dimensional
feature vector of pixel i ; the whole feature data matrix is X . The K
means algorithm minimize the functional W (C ), with a two steps
iterative procedure: for a given clustering C , W (C ) is minimized
with respect to cluster centers {y 1 , · · · , y k } yielding the mean of the
points currently assigned to each clusters; given a current set of
means {y 1 , · · · , y k } the function W (C ) is minimized assigning each
data point to the nearest cluster center y k , with respect to distance
d(., .). Such procedure is repeated until convergence is reached.
The main drawback of clustering techniques is the lack of spatial information, which in the case of ultrasound images causes
noisy segmentations, with irregular jagged-edges, wrong classified
pixel and small isolated ROIs. Such issue can be mitigated applying a MRF based regularization algorithm on the segmentation obtained with the clustering procedure, as proposed in [37].
Image segmentation is formulated as a maximum a posteriori
(MAP) estimation problem, involving the estimation of the estimated image Σ by maximization of P (Σ|X) ∝ P (Σ)P (X|Σ). The data
to be estimated is the assignment of each pixel i to a class k, indicated as Σi = k and is modelled through the a priori probability P (Σ) as a Gaussian Markov random field on a 2D rectangular
grid. The observed feature j in pixel i , X i,j , belonging to class k
is modelled through the conditioned probability P (X|Σ) as a Gaussian distribution with mean value µki,j and standard deviation σki,j .
Since the MAP estimation by maximization of the true P (Σ|X) is
computationally infeasible, a complexity reduction is mandatory.
One of the most used suboptimal minimization algorithm is the Iterated Conditional Modes (ICM), first proposed in [37] and applied
to different medical image segmentation problems. The ICM minimizes the MAP estimator no more on the entire grid, but pixelby-pixel. In other words, we have to minimize for every pixel with
36
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
respect to the assignment Σi = k the functional:

³
´2 
X i,j − µki,j
X ³ k ´
 X ¡ ¢
U (Σi |X) =
log σi,j +
³
´2  + V Σi
j =1
δi
2 σki,j
Nf
(2.6)
where δi is the neighbourhood system centered on the pixel i and
V is the MRF potential function defined over a second order neighbourhood system. Intuitively the first term of (2.6) imposes that
class intensity is close to observed data, while the second term is a
regularization constrain suggesting the algorithm to associate the
same class to close pixels. This technique requires the minimization of several simple functions, rather than the optimization of a
single, much more complex, functional. This results in an iterative
algorithm which drastically reduces computational requirements,
and provides a local optimum very close to the global one, if initialized with K -means segmentation.
The algorithm was implemented following the adaptive methodology proposed in [37]: conditional mean and conditional variance are estimated over windows whose size is reduced after a certain number of iteration until algorithm convergence is reached,
i.e. either the number of modified pixels is lower than a fixed threshold or the maximum number of iteration has been reached.
Finally, the different clusters are examined and further divided
according their spatial distribution: all the pixel that forms a connected and closed area are identified as a single region of interest.
If ROI of a minimum size are required, all the regions smaller than
a fixed number of pixel are reassigned to the adjacent regions applying a morphological operations on them.
2.3 Deconvolution as pre-processing step
An important step in data preparation for analysis is the compensation for system-dependent effects. In this work, this problem is
faced through the use of deconvolution, which is able to reduce
noise and increase both contrast and quality of US images. Although deconvolution is a well known topic, its usefulness in CAD
schemes was rarely explored.
One of the limitations on image quality is due to the blurring
effect on the back-scattered echoes produced by the transducer’s
point spread function. The acquired RF-signal x[n] can be considered to be the result of the transducer’s impulse response h[n] con-
2.3. D ECONVOLUTION AS PRE - PROCESSING STEP
37
volution with the tissue reflectivity function (or tissue response)
w[n]:
x[n] = w[n] ∗ h[n]
(2.7)
Given the imaging system’s output x[n], the deconvolution algorithm produces an estimate of the tissue response w[n].
Standard techniques for in vivo US signal deconvolution require the knowledge of system PSF and are based on regularized
inverse solution of (2.7), obtained imposing L 1 or L 2 norm constraints on the solution [38]. The system impulsive response h[n]
in vivo can not be properly estimated on phantoms because tissues cause unpredictable phase aberration, non linearity and dispersive attenuation. A minimum phase version of system PSF can
be estimated from the observed signal through homomorphic deconvolution techniques [39]. The main drawbacks of these techinques are the high sensitivity of the obtained solution from the
estimated PSF, which often results in image artifacts [40], loss of
information and high computational cost.
Predictive deconvolution techniques that are used as standards
in astronomy and geophysics are less popular in medical ultrasound. Although blind deconvolution techinques are intrinsically
more limited, their choice for tissue characterization is due to several aspects: they do not require PSF estimation, they are simpler
to implement and have low processing time. Furthermore, blind
techniques based on linear prediction perform a local whitening
of the signal, taking into account time-varying PSF, which is a realistic situation in ultrasound due to the variable focusing and pulse
attenuation.
In this work, techniques like Wiener deconvolution with minimum-phase PSF estimation [41] have been compared to the linear
prediction technique. We verified that Wiener technique provides
over-smoothed solution and features extracted from deconvolved
images are not significant for tissue characterization. On the contrary, predictive deconvolution [42] based on linear prediction is
able to highlight diagnostic significance of extracted features.
Predictive deconvolution uses a predictive filter to discard any
deterministic aspects from x[n] and consider the unpredictable
part, that is the error respect to real values of x[n], as the restored
tissue response. Thus it is necessary to estimate the coefficients of
a filter ǫk to predict future values of x[n]:
X
k
ǫk x[n − k] = x̂[n + α]
(2.8)
38
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
Figure 2.3: Autoregressive convolutional model and predictive deconvolution
iterative algorithm.
with x̂[n + α] being an approximation of x[n] after α samples. The
non predictable part of the input signal is represented by the error:
e[n + α]
=
=
x[n + α] − x̂[n + α]
P
x[n + α] − k ǫk x[n − k]
(2.9)
The Recursive Least Squares (RLS) algorithm [43] provides a recursive way to compute the filter which minimizes the weighted
least squares error function:
E (ǫ, n) =
n
X
λn−i |e[i ]|2
(2.10)
i=1
The information of new data updates the old estimate and the forgetting factor λ decreases the weight of data in the distant past to
follow the statistical variations of the observable data. The imaging system model and corresponding deconvolution model are visualized in figure 2.3. It can be shown that the output of the deconvolution process, converges to the tissue response w[n] if the
following hypotheses are fulfilled: i) Feedback hypothesis: h[n] has
an all pole transfer function, ii) Random hypothesis: w[n] is a white
Gaussian noise process, i.e. the RF signal x[n] is modeled as an autoregressive process [44].
Other deconvolution algorithms such as inverse filters as in
[45], or wavelet-based approaches as in [46] and [47] were considered and tested but the RLS prediction resulted to be the best
compromise in terms of classification performance and computational cost.
2.4. F EATURE - BASED CLASSIFICATION
39
2.4 Feature-based classification
The last stage of the proposed method is a supervised classifier.
In this phase, the selected ROIs are classified by means of a set of
features computed from the regions outlined by these ROIs. The
result of the classification will be displayed over the standard Bmode image: position and extension of all cancerous ROIs are
marked with false colors. Such information can be used by radiologist in guiding biopsy protocols.
Previous studies reported that features of different nature
hardly are linearly inter-correlated and that a nonlinear classification model can be able to extract valuable information from a
mixed feature set and to reach higher level of accuracy [25]. On the
other hand, complex model classifiers fail in preserving physical
significance of features and their true dependence on pathology.
As regards tissue characterization based on irregular ROIs,
among the possible supervised machine learning techniques a Support Vector Machine (SVM) classifier was adopted [8]. As reported
in the literature, nonlinear classifier seems to be preferable for prostate tissue characterization and SVM have been proved to be a
good choice. In order to obtain a nonlinear SVM classifier a Gaussian radial basis function was adopted as kernel. All the simulations were performed using the publicly available software implementation of SVM classifiers (svmlight [48]).
Likewise, in the analysis of deconvolution in tissue characterization, a nonlinear classification model was designed and trained
to discriminate between healthy and cancerous zones. The studied classifier is based on a first nonlinear feature extraction step
and a second linear classification step. Empirical results on the
analysed dataset lead to choose the Generalized Discriminant Analysis (GDA) as feature extraction algorithm. This technique [49]
uses kernel transformation to obtain a new feature space F that
is related to the former by a nonlinear mapping φ, performed on
feature vector x:
φ : RN
f →F
x → φ(x)
(2.11)
The subsequent classifier is a FLD that finds the best hyperplane separating healthy and unhealthy samples of the training
set. The resulting new GDA features, z, are projected in the direction which maximizes the between-class distance and minimizes
the within-class distance of samples of the mapped training set.
This is a linear operation realized by simply combining mapped
40
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
feature vectors:
d(z) = 〈w, z〉 + b
(2.12)
where b is an offset chosen to impose the same distance between
centroids of the two classes of samples and hyperplane, while the
weight coefficients w are found by maximizing the following Fisher
criterion:
wT S B w
(2.13)
D F (w) = T
w SW w
In the previous expression S B is the between-class scatter matrix,
while S W is the within-class scatter matrix, referred to the
remapped samples of the training set.
involves the computation of dot products
­ This maximization
®
φ(xi ), φ(x j ) in F : such calculation can be performed efficiently
by choosing kernel functions K(xi , x j ) which act like dot products
in the re-mapped feature space [9]. The choice of these kernel
functions is crucial. In this study different kind of kernels (namely
polynomial, gaussian and sigmoidal kernels) were tested. The best
results were obtained when a Gaussian kernel is applied.
The effect of GDA feature projection is compared to the results
of a simpler linear (i.e. without nonlinear mapping) discriminant
analysis (LDA) in figure 2.4. In particular, the histograms in figure
2.4(a) show the distributions of healthy (in dark) and unhealthy
(in white) samples after LDA projection obtained combining 5 features, while figure 2.4(b) shows the GDA projection of the same
samples: the improvement in separation is evident.
This classification system was compared to different mining
techniques, both linear and non linear (regularized least squares,
elastic net, SVM), and was selected because the computational time
required for its training is significantly smaller with respect to the
other tested classifiers, while providing the same accuracy.
A scheme of the whole classification procedure of an ultrasound
image is provided in figure 2.5.
2.5 Ground truth database
In this study the analyzed TRUS images of prostate sagittal section
are acquired by a standard US System Esaote Megas™ with transrectal transducer EC123 at sampling frequency of 50 MHz and central frequency of 7.5 MHz. Through this commercial US equipment (combined with FEMMINA [50]) it is possible to have direct
access to RF echo signal before time gain compensation sampled
2.5. G ROUND TRUTH DATABASE
(a)
41
(b)
Figure 2.4: Comparison between linear and non linear projection on real data.
figure 2.4(a) and figure 2.4(b) show the histograms of healthy (in dark) and unhealthy (in white) samples (5-feature vectors) after LDA and GDA projection respectively
Figure 2.5: Classification Scheme: processing steps from RF signal measurements to the derivation of final classification label.
with a resolution of 12 bits. The RF signal is essential to compute
some spectral parameters, useful to characterize prostate tissues.
The clinical study involves 37 patients with a high level of PSA
and undergoing their first biopsy. The age of patients ranges from
53 to 75 years. For each patient 10 US image frames were produced
with probe held in a fixed position. Since the frame rate is about
20 fps, all the frames provide the same information unless for the
measurement noise and micro-movements of radiologist’s hand
driving the probe. Thus in this analysis only the first frame is considered, although the other frames can be employed in the study
of stability of features and possibly to extract features averaging on
several frames. The available dataset for this investigation consists
of 15 benignant cases (normal and hyperplastic tissue) and 22 malignant cases (presence of adenocarcinoma).
The ground truth comes from histopathological analysis from
tissue samples acquired through biopsy. The matching of biopsy
results on TRUS images is performed by a single expert radiologist,
using the exact information of needles location in the images. For
42
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
each case the radiologist has matched biopsy results with precise
regions on TRUS images: both pixels found to be part of cancerous zones and pixels belonging to healthy regions are selected and
labeled in every zone examined through biopsy.
Cancer tissue which has not been sampled during this procedure will be labeled as healthy, thus the accuracy of this process depends on the detection rate of biopsy, gold standard in this
study.
In reality, unless radical prostatectomy is performed, with high
probability the portion of tissue sampled during biopsy is smaller
than the true cancer extension, the available ground truth must
be considered as an incomplete labeling. The learning problem is
thus not fully supervised and specificity of the classification can’t
be evaluated in a reliable way.
However, since the main goal of computer methods for prostate cancer diagnosis is the identification of all the cancer foci for
guiding biopsy, possible overestimations of the size of cancerous
region does not affect the clinical performance of the computer
aided diagnosis procedure. Thus, the main goal is the development of a high sensitivity classification method which is able to
identify all the cancer foci with an average level of specificity.
As regards the first study about segmentation, the whole images database is segmented using the method described in section
2.2, obtaining about 300 ROI of irregular size per image.
On the contrary, in the analysis of the proposed CAD scheme
involving predictive deconvolution, described in section 2.3, the
images of 2500 × 96 pixels are segmented in rectangular windows.
The choice of window size is a tradeoff between statistical significance of samples in a ROI used to compute features and the
resolution of the output feature image. In this study ROIs are rectangular windows of 101 × 7 pixels and of about 9mm2 area.
From each image about 3000 ROIs are generated from the prostate zone and used for feature extraction.
The pairs ROI-label constitute the ground truth for both segmentation algorithm validation and classifier training. The evaluation of features is performed over ROIs lying inside the zones of
different class; ROIs centered on the border of two different zones
are not considered in the analysis.
All the examples of benignant ROI are taken only from the benign cases, because there is no sure knowledge about benign regions from malign glands not examined by histology. Furthermore,
only the ROI marked by radiologist as pathological are used as ex-
2.6. E XPERIMENTAL RESULTS
43
amples of cancerous regions in the learning phase.
2.6 Experimental results
Feature selection
Beside Unser textural features, we adopted also Haralick and Fractal textural features [36], radio frequency signal statistical modeling (generalized Gaussian [26]) and spectral parameters (central
frequency, midband and slope [51], polynomial fitting of wavelet
spectrum [52]). As described in section 2.1, since the complete feature
set
of
all
parameters
mentioned
before
would have a huge dimensionality of about 140 attributes, a first
selection step is performed, keeping for each group of features only
those highly correlated to the ground truth class, and so to the
pathology, and discarding the other ones. In this application a hybrid feature selection algorithm, MIHFS, is used to rank and prune
the initial feature set, as explained in section 2.1.
In the study based on the dataset described in section 2.5, MIHFS algorithm selects a feature subset whose cardinality is equal
to 10.
Feature extraction performed after deconvolution yields to different distributions of feature values. As a consequence, MIHFS
performed on the post-deconvolution feature set provides a new
ranked list and a different optimal subset.
FS output for the case without preprocessing and the case with
predictive deconvolution are given in table 2.2, presenting the features from the most to the less important. Listed attributes in the
selected feature sets are referred to the groups described in table
2.1. A brief description of used features is provided in appendix A.
Table 2.2: Ranked Feature Sets
Rank
No Preprocessing
RLS Deconvolution
1
2
3
4
5
6
7
8
9
10
Variance-Fractal(2) Alfa1
Variance-WDES Diff. Proj. n.8
Mean-Intercept AR(2) burg
Mean-correlation 135
Mean-Nakagami logm
Variance-Fractal(1) Beta 10
Variance-Fractal(2) Beta 10
Mean-Naka w Diff. Proj. n.8
Mean-Nakagami logOmega
Mean-homogeneity 90
Mean-Intercept2 AR(3) lmsd
Variance-Fractal(2) Alfa2
Mean-Nakagami logm
Mean-correlation 45
Variance-Fractal(2) Beta 8
Mean-Haralick Sum of Squares
Mean-contrast 90
Variance-Fractal(1) Beta 10
Mean-Nakagami logOmega
Mean-Naka w Diff. Proj. n.8
44
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
MIHFS results show that RLS deconvolution improves the predictive skills of textural features and worsens the abilities of some
spectral parameters. As a matter of fact: i) deconvolution causes
good speckle reduction aimed at decreasing textures dependence
on the acquisition system and allows an easier edge detection, rising textural pattern significance; ii) some spectral features like the
high frequency diffuse component extracted from RF signals partly
lose their diagnostic value because deconvolution acts as a whitening filter, therefore these features must be computed before PSF
compensation; iii) the diagnostic significance of Nakagami features is assured in both cases, with and without predictive deconvolution.
Irregular ROI-based tissue characterization
Concerning the study about irregular ROI segmentation, as ROI
have different sizes, the classification results must be normalized
in order to take into account the size of each ROI. Therefore, classification performance is evaluated in terms of correctly classified
pixels inside the prostate border. As stated in section 1.2, typical
measures to assess classifier performances are sensitivity, specificity, accuracy defined here pixel-wise: TP are true positive pixels,
TN are true negative pixels, FP are false positive pixels and FN are
true negative pixels. Pixel-based criteria instead of ROI-based criteria are necessary when irregular regions are involved, since ROI
can contain at the same time healthy and unhealthy ground-truth
pixels, thus, unless a threshold is defined, the class of the whole
ROI is ambiguous.
The average classification performances obtained by the proposed procedure on the dataset are resumed in table 2.3. In this
study we employed the ranked feature list of 54 attributes provided
from MIHFS algorithm.
The classification results show that the first 20 features are the
most relevant: while increasing the number of features from 10 to
20 a significant performance improvement is noted, only slightly
improvement are obtained when the complete set of feature is used.
Results confirm that a multi-feature approach is fundamental.
The high value of sensitivity obtained (always > 80%) suggests
that the proposed method is able to identify most of the known
cancer locations. As said before, since only partial information
about the real extension of cancerous areas is available and since
classification was performed on the whole image, the specificity
2.6. E XPERIMENTAL RESULTS
45
(a)
(b)
Figure 2.6: Examples ground truth reference images (left panel) and classification results displayed above B-mode prostate sagittal scans as visual guidance for
biopsy (right panel).
value obtained is not fully reliable and it gives just quantitative information about the method performance. However, analysing the
results obtained on the benign cases we observed that the number
of false positive is extremely low and we can therefore conclude
that the low specificity value is due to the over estimation of cancer size in malignant images. Since the classification results are
meant to be used for guiding biopsy, an over estimation of the tumour size is less problematic than false positive on benign cases,
which would causes unnecessary additional biopsies. Finally, in
figures 2.6(a) and 2.6(b) examples of the classification results displayed above the B-mode scans as visual guidance for biopsy are
given.
46
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
Table 2.3: Classification performance of the proposed tissue characterization
method
Settings
10 features
20 features
54 features
SE
SP
Acc
0.80
0.93
0.95
0.71
0.73
0.74
0.75
0.87
0.89
Deconvolution in tissue characterization
In this study the classification model is trained on 1000 ROIs, sampled randomly from 18 TRUS images (7 benignant and 11 malignant). Training is performed through stratified 10-fold cross validation, so the classification model is generated from 900 samples
and validated on the remaining 100 samples. The testing set consists of the other 19 images, completely unknown to the classification model.
Classification performance is computed through ROI-based
metrics on a total of 58.602 testing ROIs, where 58.286 ROIs are
healthy and 316 are unhealthy. Examples of classified images and
their corresponding ground truth are shown in figure 2.7. These
TRUS images show the axial section of a prostatic gland without
scan conversion. According to the gland anatomical shape shown
in figure 1.6 of section 1.3, the images in figure 2.7 show a prostate
in the same orientation, where the superior hyperechoic boundary marks the interface between the rectum (indicated with “R” in
figure 2.7) and the peripheral zone (PZ) of the prostate.
Classification output is visualized directly on the original TRUS
image, by means of transparent colored patches located over the
ROIs. Ground truth images are visualized on the left side of each
image pair in figure 2.7, where healthy regions are marked in yellow and cancerous regions in red. The classification output is
shown on the right side of each image pair in figures 2.7: unhealthy
classified ROIs are covered in red, while healthy classified regions
are marked in blue. Red and blue colors move toward violet and
cyan respectively when the discriminant function value is near zero,
pointing out zones where classification is more difficult. A color
bar is added to the side of each images set to support an easier
reading of classification results. figure 2.7(a) represents a benignant case, while the others pictures show glandes affected by car-
2.6. E XPERIMENTAL RESULTS
47
(a)
(b)
(c)
(d)
Figure 2.7: Examples of ground truth and classified images. The letter “R” marks
the location of rectum, while the color bar shows the colors used to mark zones
classified as healthy (H) and unhealthy (U). The gradients indicate the discriminating function values and, thus, the colors next to the threshold between the two
classes mark regions where classification is more difficult. figure 2.7(a) represents
a benignant case, while figure 2.7(b), figure 2.7(c) and figure 2.7(d) show malignant
cases.
cinoma. Since some false positives are always present, also in the
best case and always outside biopsy region, the positive predictive
value is low, but this is due to the low prevalence of unhealthy regions respect to healthy regions in a ROI-based metrics. Criteria
independent on disease prevalence are employed to evaluate this
clinical test.
In order to assess classifier performances, mean value and standard deviation of sensitivity, specificity, accuracy and area under
the ROC curve (Az) were estimated over 10 different experiments,
where images in the training and test set are randomly selected. In
table 2.4 performance measures related to the GDA-FLD classifier
with and without RLS deconvolution preprocessing step are collected (mean value ± standard deviation). Each table row refers to
a classifier trained on a different number of features from 5 to 54.
48
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
(a)
(b)
Figure 2.8: ROC of classification with and without RLS predictive deconvolution.
figure 2.8(a) and figure 2.8(b) show the performance for a 5 features and 10 features
classification model respectively.
Figure 2.9: Boxplots of distributions of ROC area for classification without preprocessing (In figure: Az 5 and Az 10 when using 5 and 10 features respectively) and
classification after RLS preprocessing (In figure: Az 5 RLS and Az 10 RLS when using
5 and 10 features respectively).
Classification without deconvolution shows sensitivity ranging
between 0.51 and 0.72, specificity and accuracy ranging between
0.89 and 0.97, and ROC area ranging between 0.87 and 0.92. Classification performances decrease as number of features increase,
when more than 10 parameters are employed. This indicates an
over-fitting on training data when a high number of attributes is
used. The addition of predictive deconvolution improves perfor-
2.6. E XPERIMENTAL RESULTS
49
mances of 3% in terms of ROC area when 5 or 10 features are used,
1% and 4% for the case with 20 and 54 features respectively, where
the effect of overfitting is evident. In the best case classification
performances after RLS preprocessing provide a value of the ROC
area of 98%. These results are better visualized in figure 2.8, where
ROC curves of the two classifiers in the best case are compared
when involving 5 and 10 features.
The improvement is visible in the shifting of Az distributions
mean value toward higher values from the case without preprocessing to the case with RLS preprocessing. Az distributions (5 and
10 features) are shown in figure 2.9 for the cases with and without
RLS preprocessing.
Table 2.4: Classifier Performance
#
No Preprocessing
SE
SP
Acc
Az
0.89 ± 0.02
0.93 ± 0.02
0.95 ± 0.01
0.96 ± 0.01
0.91 ± 0.02
0.92 ± 0.02
0.92 ± 0.02
0.87 ± 0.03
#
0.72 ± 0.08
0.89 ± 0.02
0.69 ± 0.06
0.94 ± 0.02
0.64 ± 0.07
0.95 ± 0.01
0.51 ± 0.06
0.97 ± 0.01
RLS Deconvolution
SE
SP
Acc
Az
5
10
20
54
0.79 ± 0.05
0.75 ± 0.09
0.72 ± 0.07
0.62 ± 0.10
0.89 ± 0.02
0.93 ± 0.02
0.94 ± 0.01
0.96 ± 0.01
0.93 ± 0.01
0.95 ± 0.02
0.93 ± 0.03
0.91 ± 0.02
5
10
20
54
0.89 ± 0.02
0.93 ± 0.01
0.94 ± 0.01
0.96 ± 0.01
The overall result of the introduction of predictive deconvolution is an improvement of classification performances with an addition of complexity that can be negligible with respect to feature
computation complexity.
These processing times are still far from being real-time, but
further work can be done to speed up the whole classification procedure:
• use C or CUDA™ implementation of spectral and textural features
• define larger “natural” regions of interest, based on textural characteristics of the prostate, for example with the algorithm proposed in section 2.2. The additional step for irregular ROI definition does not bring extra computational
complexity, as the cost saved in feature computation due
50
Two. U LTRASOUND IMAGE ANALYSIS AND INTERPRETATION
to the analysis of a smaller number of ROIs is more important than the supplementary cost of segmentation. Thus this
step would provide a reduction of both feature generation
and classification computational time.
In table 2.5 a summary of previously published methods for
prostate tissue characterization (discussed at the beginning of this
chapter) is visualized in order to show research developments in
computer-aided cancer detection and as a comparison with the
study presented here. A more detailed descriptions of the state of
the art methods is presented in appendix B.
Table 2.5: Published methods for ultrasound-based prostate tissue characterization
Work
Basset [16]
Huynen [17]
Houston [18]
Schmitz [19]
Scheipers [20]
Feleppa [21]
Mohamed [22]
Llobet [53]
Mohamed [23]
Han [24]
HistoScanning™ [54]
Ground Truth
# patients
Technique
Features
SE
16
51
25
33
100
200
20
300
20
51
29
Textural
Textural
Textural
Multi
Multi
Spectral
Textural
Textural
Multi
Multi
Statistical
83
80
73
82
83.3
68
83.3
92
95
Results %
SP
Acc
71
88.20
86
88
100
53
100
95.9
-
80
75
80
93.75
61.6
94.4
-
Az
86
85
60.1
-
2.7. C ONCLUSION
51
2.7 Conclusion
Prostate cancer is a common disease among men and early stage
detection is crucial for a successful pathology treatment. The current clinical procedure for diagnosis and grading of the disease, i.e.
biopsy and histopathological images analysis, is still matter of controversy because it is invasive, relatively dangerous and affected
by errors due to the sampling process and the multi-focal nature
of prostate cancer. In this context the availability of a CAD tool
would be very useful to physicians to support their decision and
detect cancer even in its early stages.
TRUS images based CAD systems can be efficient tools to perform more objective prostate cancer diagnosis reducing inter operator variability. The use of ultrasound technology for prostate
cancer detection is motivated by its minimal invasiveness and cost.
Furthermore, TRUS technique is already integrated in the standard
procedure for prostate examination.
This study reported tissue characterization based on a multi
feature approach, where spectral, statistic and textural attributes
were extracted from TRUS images in the prostate zone and used
for ROI classification. A hybrid feature selection algorithm exploits
mutual information between features and pathological state of
ground truth tissue to reduce feature set dimensionality, discarding irrelevant or redundant attributes. In the proposed CAD scheme
the resized feature set is processed by a Gaussian re-mapping, shifting the problem in a sub dimensional space, where linear classification is more effective.
Interesting results were obtained by employing both rectangular regions of interest and irregular zones, selected by automatic
segmentation algorithms discussed in this thesis. This work was
also meant to explore the utility of deconvolution for classification purposes and to provide a comparison between non linear
classification performance with and without deconvolution preprocessing. To this aim, a predictive deconvolution is applied to
ultrasound data before feature computation, producing an average increase in classification performance of about 3% in terms of
area under the ROC curve.
The proposed classification procedure was designed and tested
on a 37 images ground truth and provides a CAD system performing a highly accurate detection inside the prostate zone.
Three
A CARMA MODEL FOR
ULTRASOUND SIGNALS
For my part I have sought liberty more than
power, and power only because it can lead to
freedom.
Memoirs of Hadrian (1951)
I
N ultrasonic tissue characterization several kinds of information
may be inferred from the RF signal by computing some specific
characteristics as textural or spectral features. One more way to
extract information is to find a suitable model for the echo signal
and to estimate its parameters.
Modelling the ultrasound RF signal allows to define new features, the model parameters, which may be dependent on the acquisition setting and/or on the tissue crossed during radiation, thus
providing interesting clues to reduce acquisition system effects or
to type tissues.
The autoregressive moving-average (ARMA) model, which is
classical for time-series modelling, has been employed in a shiftvariant fashion to model ultrasound signal and to perform blind
deconvolution [6], as also proposed in section 2.3, where the signal is supposed to be shift-variant auto-regressive.
The ARMA model is suitable to model a process depending on
a series of unobserved shocks as well as on its past behaviour. Since
in ultrasound images the PSF is spatially variant, ultrasound signals can be modelled as ARMA processes, provided that ARMA parameters are allowed to change over time, i.e. over the axial direction.
An easy way to deal with nonstationary models in ultrasound
is to assume that the spatial variability of the PSF is slowly varying,
hence well approximable by a piecewise constant function. In this
53
54
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
case, the image can be divided into segments, whose size is chosen to be small enough to guarantee that each of the segments is
formed by a stationary convolution with a different PSF. Each segment can thus be modelled as a different standard ARMA model.
An ARMA process x[n], driven by white noise w[n] with variance λ2 , is defined as:
x[n] +
n
X
ci x[n − i ] = w [n] +
i=1
m
X
d j w [n − j ]
(3.1)
j =1
The signal x[n] can thus be seen as the output of the linear timeinvariant casual system H (z), excited by a white noise input:
Pm
j =0
H(z) = Pn
d j z−j
c z
i=0 i
−i
(3.2)
An Autoregressive (AR) process is a special case of an ARMA
process, not depending on its own past samples and, thus, corresponding to an all-pole transfer function. There are several estimation techniques for AR and ARMA parameters in literature, mostly
based on prediction error [55].
In this chapter, a continuous-time ARMA (CARMA) model is
proposed for ultrasound signals. Since all the available information is discrete, estimation techniques for CARMA parameters are
based on sampled data, constituted in this case by the RF signal.
Using a continuous-time model for ultrasound signals means employing some function to perform interpolation of the available
sampled RF signal and move this way to the continuous-domain.
In the next section, continous-time AR (CAR) and CARMA models will be presented, then a new approach for CARMA identification from sampled data, based on exponential splines [56, 57] is
proposed. The presented identification technique is very efficient
also in some aliasing conditions. A comparison study on simulation data is described to show the accuracy of the proposed approach with respect to an algorithm based on the interpolation
of discrete data by polynomial splines. CARMA parameters can
hence be estimated from RF signal and used for its characterization.
A second, simpler algorithm is further proposed to actually extract CARMA parameters from ultrasound signals, assuming high
frequency sampling, which is usually the case in ultrasounds acquisitions. Finally an analysis of these parameters is performed to
assess their skill in capturing information about the concentration
of scatterers in ultrasound signals.
3.1. CAR AND CARMA MODELS
55
3.1 CAR and CARMA models
CARMA processes are widely used in control theory and in signal/image processing and analysis. Typical examples of applications are system identification and adaptive filtering [55, 43],
speech analysis and synthesis [58], stochastic differential equations
and image modeling [59, 60, 61]. In practice, the available data
is discrete and one is usually required to estimate the underlying
continuous domain parameters from sample values. Potential examples are continuous domain structure modeling of physical phenomena; linear time invariant (LTI) system identification as well as
numerical analysis of differential operators.
The sampled version of a CARMA process is a discrete-domain
ARMA process whose zeros and poles are coupled in a non-trivial
way [62, 63, 64, 65]. To overcome this difficulty, currently available methods use various types of approximations. Recent works
on this subject are classified into two broad categories: direct and
indirect [66, 67]. Direct methods consist first of parametrizing a
discrete-domain model by the continuous domain parameters.
The discrete-domain model is then used to minimize a cost function that involves the available data. In this way, the required continuous domain parameters are directly estimated by the minimization process. An example of such a method would be replacing
derivative operators by finite difference operations and replacing
the continuous domain white noise by a discrete-domain white
noise [68, 69, 70, 71, 72].
In [73], a continuous-time AR model was recast into a discretetime linear regression formulation rather than into a discrete-time
ARMA process. The discretization of the continuous-time white
noise was considered then to be part of the regression error. This
error also includes contributions from the weighted finite difference approximation. As the least square solution of the regression
problem might result into a biased estimation, the authors of [73]
suggest two methods for reducing the bias effect. This was done
by imposing constraints on the finite difference weights or, alternatively, by compensating for the bias at the final stage of the estimation process. The advantage of these two methods is that they
require no shifting of the data when approximating the derivative
values, giving rise to a reduced computational complexity over
other least square methods.
In [74] it is suggested to parametrize the autocorrelation sequence of the sampled process by applying the Schur decompo-
56
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
sition numerical method to obtain a state-space representation of
the discrete-domain process. The cost function that has been proposed there minimizes the ℓ2 norm of the difference between a
sampled version of the autocorrelation model and the autocorrelation sample values. Another example of a direct method consists
of power spectrum parametrization [75, 76]. Indirect methods, on
the other hand, rely on standard discrete-domain system identification methods such as prediction error minimization [55]. The
discrete-domain system is then mapped into a continuous domain
one while preserving some required properties, e.g. the bi-linear
transform. Such mappings are common practice for digital filter
designers, but alternative transformations are also available [77].
Current methods are sensitive to aliasing artifacts, meaning that
they require the signals to be sampled at a sufficient rate.
Some mathematical conventions and notations that will be used
throughout this work are provided in the next paragraph.
Notations
The bilateral Laplace transform of a scalar function ϕ(t ) is
Z∞
©
ª
Φ(s) = L ϕ(t ) (s) =
ϕ(t )e −st d t ,
(3.3)
−∞
where s takes complex values that satisfy ϕ(t )e −st ∈ L 1 . The inverse
−1
Laplace transform is ©denoted
ª L . The Fourier transform of this
function is ϕ̂(ω) ©= F ªϕ(t ) (ω) = Φ( j ω). The bilateral z-transform
of the sequence ϕ[n] n∈Z is
∞
X
©
ª
ϕ[n] z −n ,
Φd (z) = Z ϕ[n] (z) =
(3.4)
n=−∞
where z takes complex values for which the infinite sum converges.
The Discrete Time Fourier Transform (DTFT) of the same sequence
is ϕ̂d (ω) = Φd (e j ω ) and for sequences of length N , the Discrete
Fourier Transform (DFT) is given by ϕ̂d [k] = Φd (e 2π j k/N ). The argument ω denotes either radial frequency in units of [rad/timeunit] when used with the Fourier transform, or normalized frequency in units of [rad/sample] when used with the DTFT. continuous domain and discrete-domain convolution operations are
denoted as ∗. The transpose of a matrix Σ is ΣT , its inverse is Σ−1
and its determinant is |Σ|. E denotes expected value.
The dependence of a function ϕ(t ) or a sequence ϕ[n] on a set
of parameters θ is represented by a second input: ϕ(t ; θ), ϕ[n; θ].
3.1. CAR AND CARMA MODELS
57
The same notation is employed when the function or the sequence
are expressed in the Fourier, Laplace or z domain: ϕ̂(ω; θ), Φ(s; θ),
ϕ̂d (ω; θ) and Φd (z; θ).
Model Description
A CARMA model that is driven by white Gaussian noise is specified
in the time domain by the following expression:
A(D)x(t ) = B(D)w(t ),
(3.5)
where x(t ) is the CARMA output signal and w(t ) is white Gaussian noise with zero mean and variance σ2 ; A(D) = D p + a1 D p−1 +
. . . + a p I is the p-th order constant coefficient differential operator that acts on the output of the system (AR part), while B(D) =
D q + b 1 D q−1 + . . . + b q I is the q-th order differential operator that
acts on the Gaussian noise (MA part). In the special case of CAR
processes B(D) = 1. In order to identify this stochastic system, we
parametrize its power spectrum by introducing a vector θ that consists of variance, poles and zeros of the system:
©
ª
θ = σ2 , r 1 , . . . , r q , s1 , . . . , s p
(3.6)
Therefore, the CARMA process x(t ) of order p is characterized by
the power spectrum:
¯2
¯ Qq
¯
j ω − r k ¯¯
k=1
Φ( j ω; θ) = σ ¯ Qp
¯
¯
j ω − sk ¯
k=1
2¯
(3.7)
where s = (s1 , . . . , s p ) are the poles of the system and r = (r 1 , . . . , r q )
its zeroes. Stability of the system (3.7) is ensured if ℜ{sk } < 0. The
one-to-one relation between the poles, the zeroes and the coefficients of the differential operators A(D) and B(D) allows one to
express (3.7) as a function of real
in the
© variables that can be used
ª
optimization process and θo as σ2 , a1 , . . . , a p , b 1 , . . . , b q :
¯
¯
¯ ( j ω)q + b 1 ( j ω)q−1 + · · · + b q ¯2
¯
¯
Φ( j ω; θo ) = σ2 ¯
¯
¯ ( j ω)p + a1 ( j ω)p−1 + · · · + a p ¯
(3.8)
Typically, CARMA identification is performed in the discrete domain where the corresponding model is an ARMA process, derived
by some discretization scheme. The main contribution of this work
58
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
is the use of exponential B-splines (E-B-splines) to establish an exact link between the discrete- and the continuous- time domain
models. Interpolating the discrete autocorrelation sequence with
E-B-splines seems to be a proper choice since the autocorrelation
function of a stochastic model (CAR or CARMA) is a linear combination of exponentials, which admits an exact representation in
an E-B-spline basis.
3.2 CARMA processes and exponential splines
Motivated by the deterministic theory of LTI systems, we exploit
the exponential spline mathematical formulation in this work. Exponential splines provide a formal link between continuous domain convolution operators and their discrete-domain counterparts [56, 57] and they will be shown suitable for describing sampled ARMA processes, too.
Considering an ideal sampling procedure, the autocorrelation
sequence of the sampled process corresponds to sample values
of the autocorrelation function of the original continuous domain
process. It then follows that both autocorrelation measures are of
an exponential type, suggesting an exponential spline framework
for describing the relation between an ARMA process and its sampled version.
While many of the currently available estimation algorithms
are focused on base-band power spectra, it seems possible to derive an estimator that overcomes aliasing. Such an estimator could
be found useful to the area of optical imaging for which the acquisition device has limited resolution, and to the compressed sensing setup which uses a reduced number of measurements. Another potential application is resolution conversion in which low
resolution digital images are displayed on a high resolution device
[78]. From a numerical perspective, Whittle’s likelihood function
plays an important role in deriving frequency-based estimation algorithms. This function is often approximated by means of discrete Fourier transform values and by Riemann sums. Such an approximation is not necessarily optimal and there may exist other
numerical schemes that better approximate it.
This work provides a rigorous derivation of a maximum likelihood estimator of continuous domain ARMA parameters from
sampled data. It utilizes the exponential B-spline framework while
introducing an exact zero-pole coupling for the sampled process.
3.2. CARMA PROCESSES AND EXPONENTIAL SPLINES
59
For that purpose, the relation between the autocorrelation function and the autocorrelation sequence is investigated in both timeand frequency domains. Based on this relation, it is shown that
the Cramér-Rao bound can be made arbitrarily small by considering more sample values where the sampling interval can take an
arbitrary value. The likelihood function of the sampled process is
investigated, too. In particular, it is shown that this function possesses local minima that originate from aliasing. The global minimum, however, corresponds to the maximum likelihood value.
The only assumption that is made throughout this study is that the
number of available samples is relatively large, allowing one to replace the whitening matrix with a digital filter. This approximation
is shown to be valid when considering expected values of the likelihood function. Experimental results indicate that the proposed
exponential-based approach is a preferable choice over currently
available methods that are restricted to relatively high sampling
rates.
A differential LTI system of order p is fully described by a rational transfer function (p > q)
Qq
ĥ(ω) = Qk=1
p
k=1
( j ω − rk )
( j ω − sk )
,
(3.9)
If such a system is driven by continuous domain white Gaussian noise, its output is a CARMA process. The spectral density
function of such a process is ϕ̂(ω) = σ2 |ĥ(ω)|2 and the autocorreσ2 is the varilation function is ϕ(t ) = σ2 {h(τ) ∗ h(−τ)} ©(t ), where
ª
−1
ance of the noise and where h(t ) = F
ĥ(ω) (t ) is the impulse
response of the system. Exponential splines provide a mathematical framework for relating continuous domain LTI systems with
their discrete-domain counterparts. The dependence of ϕ(t ) on
h(t ) suggests that these splines may be equally helpful for relating
continuous domain and discrete-domain ARMA processes, too.
Motivating Example: First Order AR Process
The autocorrelation function of a continuous domain AR(1) process that has a pole at s = s1 and a unit variance innovation is the
symmetric exponential
ϕ(t ; s1 ) = −
1 s 1 |t |
e
.
2s1
(3.10)
60
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
The spectral density function is then
©
ª
Φ( j ω; s1 ) = F ϕ(t ; s1 ) (ω) =
1
.
( j ω − s1 )(− j ω − s1 )
(3.11)
Upon sampling, the autocorrelation sequence of the corresponding discrete-domain process is
ϕ[n; s1 ] = −
1 s 1 |n|
e
,
2s1
(3.12)
where a unit-time sampling interval was assumed for simplicity.
The spectral density function of the sampled process is then
©
ª
e 2s 1 − 1
1
.
Φd (e j ω ; s1 ) = Fd ϕ[n; s1 ] (e j ω ) =
·
s
−
j
ω
1
2s1 (1 − e e
)(1 − e s 1 e j ω )
(3.13)
The important observation is that one is able to link the continuous domain and the discrete-domain autocorrelations via the
Shannon-like interpolation formula (see also figure 3.1)
ϕ(t ; s1 ) =
∞
X
ϕ[n; s1 ]β(t − n; s1 ),
(3.14)
n=−∞
where β(t ; s1 ) is an interpolating basis function whose Fourier expression is
β̂(ω; s1 ) =
2s1
(1 − e s 1 e − j ω )(1 − e s 1 e j ω )
·
.
−1
( j ω − s1 )(− j ω − s1 )
e 2s 1
(3.15)
Observe that the latter expression is also equal to the ratio of (3.11)
and (3.13). The corresponding time-domain expression is
(
sinh[s 1 (1−|t |)]
|t | ≤ 1
sinh s 1
.
(3.16)
β(t ; s1 ) =
0
|t | > 1
The key property that will be exploited in this work is that β(t ; s1 )
is compactly supported, which is not obvious a priori from the
Fourier-domain expression (3.15). In fact, β(t ; s1 ) is an exponential B-spline and the above method generalizes for higher order
systems.
General Case
As stated before, a continuous domain ARMA process is fully characterized by its parameters vector, θ.
3.2. CARMA PROCESSES AND EXPONENTIAL SPLINES
0.5
61
0.5
0.45
0.4
0.4
0.3
0.35
0.3
0.2
0.25
0.1
0.2
0
0.15
0.1
−0.1
0.05
0
−6
−4
−2
0
2
4
−0.2
−4
6
−3
−2
−1
0
t
t
(a)
(b)
1
2
3
4
Figure 3.1: Exponential B-spline decomposition of an ARMA model. Shown is
the autocorrelation function of an ARMA process (black, dashed line). This function can be decomposed by a weighted sum of exponential B-spline functions (various colors, various widths, dashed or solid lines). Shown are only few exponential
B-spline functions. The parameters of the ARMA process are 3.1(a) σ = 1, s = −1
and 3.1(b) σ = 1,r = −1, s = −1 ± i .
The Laplace transform of the corresponding autocorrelation
function is given by
2
Qq
Φ(s; θ) = σ Qk=1
p
k=1
(s − r k )(−s − r k )
(s − sk )(−s − sk )
.
(3.17)
By performing the partial fraction decomposition of Φ(s; θ) (we are
assuming for simplicity that the poles are simple), we find that
Φ(s; θ) = σ2
p
X
αk (θ)
k=1
½
¾
1
1
−
,
s − sk s + sk
(3.18)
and we deduce that the autocorrelation function is a sum of exponentials,
p
X
αk (θ) · e s k |t | .
(3.19)
ϕ(t ; θ) = σ2
k=1
In cases where Φ(s; θ) introduces pole multiplicity, ϕ(t ; θ) would
involve polynomial multiplications as given in table 3.1.
Ideally sampling a continuous domain ARMA process yields
a discrete-domain ARMA process. The autocorrelation sequence
of the discrete-domain process is then given by the ideal samples
of the autocorrelation function. For the partial decomposition of
(3.19) and for a unit-time sampling interval, this sequence is given
62
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
by
ϕ [n; θ] = σ2
p
X
αk (θ) · e s k |n| .
(3.20)
k=1
Autocorrelation functions that consist of multiple poles can be discretized in a similar manner as given in table 3.1.
Continuous Domain
Laplace
1
s−s k
1
(s−s k )2
1
(s−s k )3
1
+ −s−s
+
+
k
1
(−s−s k )2
1
(−s−s k )3
Discrete Domain
Time
index
z-transform
e s k |t |
e s k |n|
− z −1e−e −sk + z −1e−e sk
|t | · e s k |t |
|n| · e s k |n|
1 2 s k |t |
2t ·e
1 2
s k |n|
2n ·e
−sk
·
"
e −sk
z −1 −e −sk
+
e −2sk
(z −1 −e −sk )2
sk
¸ ·
sk
+ z −1e−e sk +
e 2sk
(z −1 −e sk )2
¸
#
e −s k
e −2s k
3
1
e −3s k
+ ¡
−
¢ +¡
¢3 +
2 z −1 − e −s k 2 z −1 − e −s k 2
z −1 − e −s k
"
#
e sk
e 2s k
3
e 3s k
1
+ ¡
+
¢ +¡
¢3
2 z −1 − e s k 2 z −1 − e s k 2
z −1 − e s k
3.2. CARMA PROCESSES AND EXPONENTIAL SPLINES
Table 3.1: conitnuous-to-discrete domain mapping of multiple poles.
63
64
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
ARMA models are closely related to generalized exponential Bsplines. These finite-support functions stem from Green’s functions of rational operators [57], and in our case the Green’s function is the autocorrelation function itself. Unlike [57], however,
this work introduces non-causal symmetric B-splines.
Definition 1 (Symmetric Exponential B-spline). The exponential
B-spline β(t ; θ) with parameters θ is specified by the following inverse Fourier transform
β(t;θ) = F −1
(
q
Y
( j ω − r k )(−j ω − r k ) ·
k=1
)
p
Y
(1 − e j ω+s k )(1 − e − j ω+s k )
(t)
( j ω − s k )(−j ω − s k )
k=1
(3.21)
Observe that this function corresponds to the ratio of the ARMA
power spectrum and of a discrete-domain AR power spectrum with
poles at e s k . The exponential B-splines kernels have the following
properties:
£
¤
• Compactly supported function within the interval −p, p .
• Bounded and symmetric functions.
• Smooth functions: the first 2(p − q) − 1 derivatives are L 2
functions.
• Their integer shifts form a Riesz basis.
• Reproduce exponential functions of the type (3.19).
• The convolution of two exponential B-spline kernels yields
another B-spline kernel of an augmented order. This property allows one to iteratively construct an exponential B-splines
of any order.
Definition 2 (Localization filter). The localization filter with parameters θ is specified by the Laplace transform
∆(s; θ) =
p ¡
Y
k=1
¢¡
¢
1 − e s+s k 1 − e −s+s k .
(3.22)
Proposition 1. The autocorrelation function of an ARMA process
with parameters θ can be written as
X
ϕ(t ; θ) = σ2 p[m; θ] · β(t − m; θ),
(3.23)
3.2. CARMA PROCESSES AND EXPONENTIAL SPLINES
65
where β(t ; θ) is the exponential B-spline with parameters θ and the
sequence p is given in the z-domain by
1
1
.
¢=
¢¡
¡
−1
s
s
∆(ln z; θ)
1−e k z 1−e k z
k=1
P d (z; θ) = Qp
(3.24)
Proof. Take the Fourier transform of (3.23) and substitute (3.21)
and (3.17).
Definition 3. The discrete exponential B-spline kernel with paramPp
eters θ is given by Bd (z; θ) = n=−p β[n; θ]z −n and β̂d (ω) = Bd (e j ω ),
where the sequence β[n; θ] contains the samples of β(t ; θ).
We rely on [65] and assume that ideal sampling preserves the
autocorrelation values of the continuous domain ARMA process.
The power spectrum of the sampled process is then related to the
continuous domain power spectrum through aliasing and it can be
described by means of a rational transfer function in the z-domain
[65]. The following Theorem suggests a computationally efficient
way for extracting the discrete-domain power spectrum from the
continuous domain parameters.
Theorem 1. The discrete-domain ARMA process that is given by the
ideal unit-interval samples of the continuous domain process (3.17)
has the following parameters:
Qp−1
Hd (z; θ) = σd (θ) Qk=1
p
k=1
where
(1 − νk (θ)z −1 )
(1 − ρ k (θ)z −1 )
,
(3.25)
ρ k (θ)
=
e sk
(3.26)
νk (θ)
=
roots of Bd (z; θ) inside the unit circle
(3.27)
and where
Bd (1; θ)
σ2d (θ) = σ2 Q p−1
(1 − νk (θ))2
k=1
(3.28)
is the variance of the discrete-domain innovation process.
Proof. As P d (z; θ) is an all-pole filter, the zeros of
Φd (z; θ) = σ2 P d (z; θ) · Bd (z; θ)
(3.29)
66
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
originate from the zeros of Bd (z; θ) only. Those zeros appear in
reciprocal pairs due to the symmetry property of β(t ; θ). Also, the
DFT of the sampled version of the exponential B-spline satisfies
β̂d (ω; θ) =
∞
X
β̂ (ω + 2πk; θ) = Bd (e j ω ),
(3.30)
k=−∞
and it does not vanish on ω ∈ [0, 2π) due to (3.21). β(t ; θ) is of finite
support so that Bd (z; θ) has no poles. The symmetry property of
∆(s; θ) indicates that the poles of P d (z; θ) appear in reciprocal pairs
and it then follows that Φd (z; θ) can be described by a minimumphase filter. The finite support property of β(t ; θ) and the structure
of ∆(s; θ) also guarantee that this minimum-phase filter has a rational transfer function. We further observe that,
Bd (z; θ)
P d (z; θ)
=
Bd− (z; θ) · Bd+ (z; θ)
(3.31)
=
P d− (z; θ) · P d+ (z; θ)
(3.32)
where ‘−’ denotes a causal filter and where Bd+ (z; θ) = Bd− (z −1 ; θ)
¢
Qp ¡
and P d+ (z; θ) = P d− (z −1 ; θ). In particular, P d− (z;θ) = 1/ k=1 1 − e s k z −1
and Bd− (z; θ) has zeros only; those zeros originate from the roots of
Bd (z; θ) inside the unit circle. It then follows that,
Φd (z; θ) = Hd (z; θ) · Hd (z −1 ; θ),
(3.33)
where Hd (z; θ) = σ · P d− (z; θ) · Bd− (z; θ) is the required filter. Imposing the structure of (3.25) results in the expression for σ2d (θ).
Beside equation (3.28) there is another way of computing the
analytical expression of the variance of the discrete innovation signal based on the classical definition of exponential B-spline [56].
The difference with the previous formulation is a multiplicative
term (in the frequency domain) in the definition of P d (z; θ) (3.24)
and of the E-B-spline £β(t ; θ) (3.21).
¤ The new definition of P d (z; θ)
Qp
differs from (3.24) by z p k=1 e s k , likewise the new definition of
£
¤
Qp
β(t ; θ) differs from (3.21) by (−1)q e − j ωp k=1 e −s k . These differences are clarified in appendix D. Provided these changes, we can
derive the following corollary.
Corollary 1. The variance of the innovation signal of the sampled
CARMA process has the following form,
σ2d (θ) = σ2 β[1; θ]
Y p−1
k=1
(−1/νk (θ))
Yp
k=1
e sk
(3.34)
3.2. CARMA PROCESSES AND EXPONENTIAL SPLINES
67
Proof. From equation (3.29) we derive
Φd (z; θ)
=
σ2 z p
p
Y
k=1
e s k · Qp
k=1
Bd (z; θ)
(1 − e s k z −1 )(1 − e s k z)
(3.35)
Bd (z; θ) is a reciprocal polynomial of order 2p − 2 due to its symmetry and, in order to make its reciprocal roots explicit, it can be
expressed in the following manner
Bd (z;θ)
=
2p
X
β[k,θ]z −k
(3.36)
k=0
=
β[1,θ] · z −1
=
β[1,θ] · z −p
Y p−1
k=1
(1 − νk (θ)z −1 )(1 −
1
z −1 )
νk (θ)
Y p−1 −1 Y p−1
(1 − νk (θ)z −1 )(1 − νk (θ)z)
k=1 ν (θ)
k=1
k
Recalling equation (3.33) and replacing Bd (z; θ) in (3.35) results
in the expression for σ2d (θ).
Corollary 2. Let θ be known. Then the autocorrelation function of
a continuous domain ARMA process is uniquely defined by its samples. Further,
ϕ(t ; θ) =
∞
X
ϕ(n; θ) · η(t − n; θ),
(3.37)
n=−∞
where the interpolation kernel η(t ; θ) is specified by its Fourier transform,
β̂(ω; θ)
.
(3.38)
η̂(ω; θ) =
β̂d (ω; θ)
This is the generalization of (3.14) for arbitrary ARMA(p,q) processes.
Proof. It was shown that Φd (z; θ) = σ2 P d (z; θ) · Bd (z; θ) and that
Φ(s; θ) = σ2 B(s; θ)/∆(s; θ). Dividing the two equations while recalling that P d (z; θ) = ∆(ln1z;θ) yields the required result.
The exponential interpolation function η(t ; θ) provides a mean
for interpolating continuous domain ARMA models. It can also be
interpreted as a spectral weighting function that relates the power
spectrum of the discrete-domain sampled process with the power
spectrum of the continuous domain process it originates from. Unlike the polynomial-based weighting function of [76], this weighting function is parameterized, allowing one to describe band-pass
power spectrum, too.
68
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
3.3 Insights on the exponential spline framework
The motivation for using an exponential spline parametrization
for CARMA models stems from the fact that the ideal sampling
model preserves autocorrelation properties. The autocorrelation
functions of continuous-time ARMA processes are exponential splines and the proposed framework has several advantages in this regard:
1. The spline framework allows one to map continuous domain
parameters to discrete-domain parameters by using simple
algebraic operations. The simple mapping of Theorem 1 can
then be used in estimation algorithms that are based on numerical optimization and that require repeated mapping operations. This mapping involves (a) calculating the autocorrelation function by a partial fraction decomposition, (b) expressing the B-spline of Definition 1 by the FIR localization
filter of Definition 2, and (c) calculating sample values of the
B-spline. As the B-spline is compactly supported in [−p, p],
the z-transform of its sampled version is a polynomial of a
relatively low degree. This in turn allows one to calculate the
zeros of Theorem 1 relatively fast.
2. The proposed framework does not depend on sampling interval values when mapping continuous domain parameters
to their discrete-domain counterparts. It provides a discretedomain process whose autocorrelation sequence is in a perfect match with the point-wise values of the autocorrelation
function for any sampling interval value.
3. The proposed framework provides a closed form expression
for the interpolation kernel of Corollary 2. Considering the
minimum-mean-square-error criterion, this kernel can be
used for interpolating the sampled process in a straightforward manner. Moreover, the interpolation task can be carried out using the finite support B-spline of Definition 1. Finite support functions are extremely beneficial for data interpolation in terms of computational complexity.
4. The proposed spline framework can naturally incorporate
non-ideal sampling procedures such as averaged sampling.
For example the zero-order-hold rectangle function is a
spline, too. In such cases, where the sampling function is
3.3. I NSIGHTS ON THE EXPONENTIAL SPLINE FRAMEWORK
69
by itself a spline, the proposed framework allows one to derive the whitening digital filter in a similar way to Theorem
1. The key point is that there always exists a B-spline that
admits the relation (3.23) with a filter p[m; θ] that depends
also on the sampling model.
An alternative mapping to Theorem 1 was given in [63, 64, 77].
Motivated by a state-space representation, these work introduce
closed-form expressions for extracting the discrete-time power
spectrum from a given continuous-time model. Nevertheless, these
expressions introduce additional computational overload
compared to the proposed approach as they rely on matrix operations. In particular, [77] involves exponentials of a matrix and
integration of matrix elements, and involves polynomial factorization that requires one to represent an inverse of a matrix in a
parametrized form. The proposed spline framework, on the other
hand, does not rely on a matrix-based representation. Rather, it
evaluates the autocorrelation function at a finite number of points
only. Owing to the spline property of the autocorrelation function, the continuous-to-discrete mapping of Theorem 1 involves
the partial fraction decomposition of 3.17 and the evaluation of the
B-spline function at only 2p − 1 points. These points give rise to a
polynomial in the z domain, whose roots are readily computed.
Each point-wise value of the B-spline function corresponds to a
weighted sum of 2p + 1 values of the autocorrelation function.
An additional continuous-to-discrete parameter mapping was
also introduced in [79]. This work derives explicit formulae for the
first terms of the Maclaurin series of the zeros and the variance of
the sampled process as a function of the sampling interval variable
T . This, in turn, allows one to analyze convergence for the limiting case of T → 0. The Maclaurin weights are given by means of
the continuous-time parameters and such results can be used for
mapping continuous domain parameters to their discrete-domain
counterparts. The advantage of the results of [79] is in its extremely
low computational complexity. Unlike the proposed spline framework, however, it does not provide an exact continuous-to-discrete
mapping; and according to Example 2 of [79], such an approximation would better be used for sampling interval values that are
considerably smaller than the characteristic time constants of the
system. As we consider prominent aliasing conditions in our work,
the approach of [79] is less applicable here.
70
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
3.4 A new approach to CARMA model identification
Let a continuous domain ARMA(p,q) process be given by its uniform samples only. Considering a large number of data points
N ≫ 1, the Fisher information matrix is given by [80],
I k,l (θ)
∼
=
N
4π
Z2π µ
0
¶
∂ϕ̂d (ω; θ)
1
·
ϕ̂d (ω; θ)
∂θk
¶
µ
1
∂ϕ̂d (ω; θ)
dω,
·
∂θl
ϕ̂d (ω; θ)
(3.39)
where k, l = 1, . . . , p + q +1 and ϕ̂d (ω; θ) is the DTFT of the sampled
ARMA process autocorrelation sequence. The Cramér-Rao Bound
(CRB) is then given by the inverse of I (θ). Sampling interval dependency can be incorporated in the localization filter (3.22) by
∆(s; θ, T ) =
p ¡
Y
k=1
¢¡
¢
1 − e s+s k T 1 − e −s+s k T ,
(3.40)
where T is the sampling interval. It then follows that ϕ̂d (ω; θ) is
not dependent on N nor does the integrand of (3.39). It then follows that, as the number of available samples becomes larger, the
CRB becomes smaller regardless of T . This implies that aliasing
effects can be compensated for by taking more measurements. It
further suggests that there exists at least an ML (Maximum Likelihood) estimator that overcomes aliasing. Such an estimator has
to take into account the fact that the discrete-domain poles and
the zeros of Φd (z; θ) are coupled, ensuring that the discrete process spectrum corresponds to z-transform of the samples of the
autocorrelation function ϕ(t ; θ). Indeed, one can estimate the digital data by a discrete-domain ARMA model while introducing no
zeros-poles coupling. However, in such a case the estimated parameters will not necessarily correspond to a continuous domain
model. For example, a possible discrete-domain ARMA(1,0) estimation algorithm may yield a pole at z = −0.5; the autocorrelation
sequence is then ϕ[n] = (−0.5)|n| which has no continuous domain
ARMA(1,0) counterpart. Furthermore, non-coupled estimators result in baseband power spectrum, making it difficult to cope with
aliasing effects.
Considering this matter further, there are at least two possible ARMA models that can be associated with the sampled data.
One of them is a discrete-domain ARMA model that introduces no
3.4. A NEW APPROACH TO CARMA MODEL IDENTIFICATION
71
zero-pole coupling. The ML estimator in this case corresponds to
minimizing the variance of the prediction error. The other is a constrained discrete-domain ARMA model which takes the sampling
process into account and which introduces the required zero-pole
coupling. The ML estimator of such a constrained model does not
necessarily minimizes the variance of the prediction error. The
constrained digital whitening filter (3.25) is guaranteed to yield a
discrete-domain innovation process that is independent and identically distributed. This does not hold true for the non-constrained
model as the innovation is not guaranteed to be independent, e.g.
when the estimated power-spectrum does not comply with the required zero-pole coupling, although it may yield lower variance
values. The continuous domain whitening filter (3.25) is always
guaranteed to exist and to whiten the continuous domain process
whereas the continuous domain whitening filter of the non constrained model, if exists [77], is not guaranteed to yield white noise.
The purpose of (3.39) is to provide motivation for deriving an
estimation algorithm that overcomes aliasing. We assume that the
number of samples is very large, allowing one to represent the discrete domain power spectrum by means of a rational transfer function. We then rely on the results of [80] which characterize the
Cramér-Rao bound by means of the power spectrum function. Our
parametrization of the discrete-domain power spectrum allows us,
in turn, to demonstrate that the Cramér-Rao bound converges to
zero for sampled ARMA processes and that it is inversely proportional to N . In this regard, the work [81] derives closed form expressions for the Cramér-Rao bound using the state space representation. Also, it assumes data sets of an arbitrary size. The closed
form expressions, however, involve matrix inversion and eigenvalue decomposition, and it is difficult to implement such computations for large data sets as indicated in section 3 of [81]. It is also
indicated there that, under certain circumstances, the Cramér-Rao
bound becomes inversely proportional to N . We demonstrated
such proportionality for every sampled ARMA process.
Motivated by the CRB, we approximate the log-likelihood function by means of a digital filter. The proposed approximation relies
on
parameters of (3.25), namely
ª
© 2 the discrete-domain
σd , ν1 , . . . , νq , ρ 1 , . . . , ρ p . These parameters can be numerically
calculated from θ in a straightforward manner as was shown in
section 3.2. We, however, do not allow the continuous domain
poles of θ to differ by 2kπ j with k ∈ Z as such poles yield the same
discrete-domain model upon sampling.
72
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
Definition 4. Let θ be known and let x be N uniform ideal samples
of the continuous domain process (3.17) taken on a unit-interval
grid. The probability density function of x is
pd f (x; θ) =
¾
½
1 T −1
exp
−
x
Σ
x
,
N
1
2
(2π) 2 |Σ| 2
1
(3.41)
where Σ [m, n] = ϕ [m − n; θ] is the autocorrelation matrix that corresponds to θ.
Definition 5. The corresponding log-likelihood function, including
a sign inversion, is
l(θ; x) = ln |Σ| + xT Σ−1 x.
(3.42)
Definition 6. Let θ be known. The digital filter g θ is given by its
Z -transform G d (z) = 1/Hd (z) where Hd (z) is given as in (3.25).
Definition 7. Let θ be known. The parameter only function κ(θ) is
given by
∞
X
n · c[n; θ]2 ,
(3.43)
κ(θ) =
n=1
where
c[n; θ] =
o
1n n
ν1 (θ) + . . . + νnp−1 (θ) − ρ n1 (θ) − . . . − ρ np (θ) .
n
(3.44)
The constant κ(θ) can be interpreted by means of an inner product operation. Considering a discrete-domain ARMA power spectrum, we define the Fourier coefficients of its logarithm as:
ln φ̂d (ω) =
∞
X
c[n; θ]e − j ωn
(3.45)
n=−∞
Recalling (3.25):
ln Hd (e − j ω ; θ) =
∞
c[0; θ] X
c[n]e − j ωn
+
2
n=1
(3.46)
it then follows that
κ(θ) = |〈
d
(ln Hd (e j ω )), ln Hd (e j ω ; θ)〉L 2 |
dω
(3.47)
3.4. A NEW APPROACH TO CARMA MODEL IDENTIFICATION
73
Theorem 2. Let θ be known and let x be N uniform ideal samples of
the continuous domain process (3.17) taken on a unit-interval grid.
Then,
³
´
lim E l(θ; x) − l˜(θ; x) = 0,
(3.48)
N→∞
where
°
°2
l˜(θ; x) = N ln σ2d (θ) + κ(θ) + °x ∗ g θ °ℓ2
(3.49)
and ∗ denotes discrete-domain convolution of an N -length output
sequence.
Proof. According to Szegö theorem for infinite Toeplitz matrices
[82, 83, 84],
n
o X
∞
n c[n; θ] c[−n; θ],
lim ln |Σ| − N · c[0; θ] =
N→∞
(3.50)
n=1
where c corresponds to the Fourier coefficients of ln ϕ̂d (ω; θ). In
our case, these coefficients are given by (3.44) and the right-handside of the equation can be replaced by κ(θ). Also, by the Residue
theorem
Z
1 2π
c[0; θ] =
ln ϕ̂d (ω; θ) dω = ln σ2d (θ).
(3.51)
2π 0
It then holds that
n
o
lim ln |Σ| − N ln σ2d (θ) − κ(θ) = 0.
N→∞
(3.52)
The constant κ(θ) also provides a mean for describing the limiting behaviour of the determinant. Writing (3.52) in a different
form [84],
|Σ(θ)| ∝ e κ(θ) · (e c[0;θ] )N
(3.53)
we observe that the right-hand-side of the equation is equal to the
determinant up to a multiplicative constant; the important thing
here is that this constant does not depend on N or on θ.
As for the term x · Σ−1 · x that appears in l(θ; x), it may be approximated by digitally filtering x and by calculating the energy of
the output. A possible filter one can use is g θ . This filter is guaranteed by Theorem 1 to exist, to be stable and to be causal. The
N -length output of such a filter can be described by means of the
lower triangular matrix
L[m, n] = g θ [m − n], m, n = 0, . . . , N − 1,
(3.54)
74
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
°
°2
and it holds that °x ∗ g θ °ℓ2 = xT LL T x. Following [85], the two ma-
trices Σ−1 and LL T are asymptotically equivalent as N becomes
larger. This stems from the fact that they both originate from the
same power spectrum ϕ̂d (ω; θ). Such a power spectrum gives rise
to an infinite Toeplitz matrix T M and the N × N matrix Σ−1 is defined by 1) truncating T M and 2) taking the inverse. The matrix L is
defined by 1) inverting T M , 2) finding its Cholesky decomposition
and 3) truncating the lower triangle matrix.
The asymptotically equivalent matrices are covariance matrices that describe asymptotically equivalent random processes. In
particular, let σ21 = E(xT Σ−1 x) and σ22 = E(xT LL T x) correspond to
the variance values of asymptotically equivalent stochastic processes. Then, limN→∞ (σ21 − σ22 ) = 0, which proves the theorem.
Asymptotic equivalence implies that the norm of the difference
between the two matrices converges to zero with increasing values
of N . Focusing on a Toeplitz matrix T M [m, n] = t M [m −n], one can
associate a discrete-time Fourier transform to the sequence t M [k]
and, for autocorrelation matrices, this transform corresponds to
the power spectrum function. Asymptotic equivalence of two covariance matrices implies that the power spectrum of one matrix
converges uniformly to the power spectrum of the other. In the
context of this proof, the matrix L describes a truncated version
of the filter g θ . As g θ ∈ l 1 , its discrete-time Fourier transform uniformly converges to H (e1 j ω ) . This means that with increasing numd
bers of samples, the expression x T L converges in the l 2 sense to
x∗ g θ . This fact stems from Parseval’s property of the discrete-time
Fourier transform. It then follows that the output is a stochastic process with a power spectrum that converges to a constant
function. This constant is σd (θ). The term xT LL T x is then an
estimation for this value. The expression xT Σ−1 x amounts to decorrelating the process x and to estimating the variance of the uncorrelated process. The value of this variance is σd (θ), too. It then
follows that E{xT LL T x} = E{xT Σ−1 x} = N σd (θ).
The proposed approach is related to Whittle’s approximation
of the log-likelihood function of (3.42). Whittle’s approximation is
given by
N
l˜W (θ; x) =
2π
Z2π
0
1
ln ϕ̂d (ω; θ)dω +
2π
¯2
Z2π ¯¯
x̂d (ω)¯
0
ϕ̂d (ω; θ)
dω.
(3.55)
3.4. A NEW APPROACH TO CARMA MODEL IDENTIFICATION
75
It holds that
N ln σ2d
°
°2
lim °x ∗ g θ °ℓ2
N→∞
=
=
N
2π
Z2π
ln ϕ̂d (ω; θ)dω
¯2
Z ¯
1 2π ¯x̂d (ω)¯
dω,
2π 0 ϕ̂d (ω; θ)
0
(3.56)
(3.57)
where the limit is required due to the use of finite-length convolution in l˜(θ; x). It then follows that limN→∞ l˜(θ; x) = l˜W (θ; x) +
κ(θ). The integrals of Whittle’s likelihood function are often approximated by Riemann sums that involve DFT values. DFT values of a finite length signal, in this case x, indeed coincide with
the samples of its Z -transform on the unit circle. Nevertheless,
this is not the case for infinite sequences that are described by
Φd (z; θ) and by 1/Φd (z; θ). The proposed approach, on the other
hand, does provide an accurate evaluation of these integrals and
the truncation error of the convolution operation is exponentially
decreasing, providing a higher convergence rate than the Riemann
sum method.
The proposed ML approach differs from currently available
methods in several aspects:
• It considers exponential autocorrelation models for both the
continuous- and the discrete-domain processes whereas several previous works considered polynomial models.
• It establishes a link between the continuous- and the discretedomain models by incorporating the sampling process into
the problem formulation. This, in turn, allows for the discrete model to stem naturally from the continuous domain
formulation while no a priori assumptions are made on the
digital data.
• The log-likelihood function suggested here holds true for any
value of sampling interval, rather than describing the limiting case of T → 0.
• The log-likelihood function considers discrete-domain data
for determining continuous domain statistics while no approximation of continuous domain frequency spectrum or
impulse responses is required.
The log-likelihood function (3.49) has several local minima as
demonstrated by figure 3.2 where simulation results are shown for
76
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
a continuous domain AR(2) process having poles s1,2 = −1 ± 5i .
Given a sampled version of such a process, the log-likelihood function (3.49) was repeatedly minimized using different initial conditions; the initial conditions are the real and imaginary part of the
poles one starts with, and they are indicated by the x and y axes
respectively. A detailed description of the five ring-like regions
is given in table 3.2. These local minima originate from aliasing
and there exist several continuous domain processes that result in
similar discrete-domain power spectrum upon sampling. These
very processes generate the local minima. Shown in figure 3.3 is
Φd (e j ω ) for the AR processes of table 3.2 while considering a unit
sampling interval. These processes correspond to the local minima of figure 3.2. The aliased spectrum shown here resembles each
other in terms of their bandpass nature and in terms of the frequency of maximum response. Region 1, however, exhibits a lowpass signal, emphasizing the difference between the proposed ML
estimation approach and other base-band methods. The peak response of these power spectra are distributed along the frequency
axis in distinct bands as shown in figure 3.4; these frequency bands
are π [rad/time-unit] wide. This property suggests that every local minimum can be obtained by minimizing (3.49) while allocating initial conditions that correspond to a peak response at the
required band, e.g. at (k + 1/2 )π where k is the band index [86].
Another way of determining initial conditions will be described
later on. Following Theorem 2, the global minimum of (3.49) corresponds to a ML estimator and it is suggested here to minimize
the likelihood function using several initial conditions.
The work of [87] investigates the maximum likelihood estimator of a discrete-domain ARMA process, too. In particular, it characterizes the local minima values of the likelihood function. The
analysis of [87] relies on the equivalence of the maximum likelihood estimator with the prediction error variance minimization
estimator. According to Theorem 2 of [87], there is a unique global
minimum, provided that enough parameters are used in the fitted model. The ARMA model that originates from sampled data,
however, introduces zero-pole coupling which is not considered
in [87]. Also, the maximum-likelihood estimator for such cases
does not necessarily correspond to minimizing the variance of the
prediction error. For these reasons, the results of [87] can be partially applied here. Similar to [87], we also assume a large number of samples, allowing one to formulate the power spectrum by
means of a ratio of two polynomials. We also demonstrate that
3.4. A NEW APPROACH TO CARMA MODEL IDENTIFICATION
Max
77
10
(5)
9
(4)
7
(3)
6
5
(2)
4
3
(1)
Imaginary part of pole
8
2
1
Min
−10
−8
−6
−4
−2
Real part of pole
Figure 3.2: Local minima of the likelihood function.
the likelihood function of the sampled data has a global minimum.
Nevertheless, due to the zero-pole coupling that is introduced by
the sampling process, this global minimum differs from the global
minimum of [87]. Furthermore, the sampling process introduces
local minima that do not coincide with the global minimum and
this property allows one to select one set of parameters which is
the most probable.
Table 3.2: A detailed description of figure 3.2.
Region
Log-likelihood
Value
Estimated
Poles
Frequencyb
[rad/time-unit]
1
-3656
−2.4 ± 2.1i
2a
-3675
−0.96 ± 5.0i
3
-3557
−0.46 ± 7.8i
7.8
4
-3067
−0.25 ± 11.1i
11.1
5
-2212
−0.16 ± 14.0i
14.0
a Global minimum of the log-likelihood function.
b Frequency of maximum response of Φ ( j ω).
L
0
4.9
78
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−4
Region 1
Region 2
Region 3
Region 4
Region 5
−3
−2
−1
0
1
2
3
4
Radial Frequency [rad]
Figure 3.3: Discrete-domain spectrum of sampled AR signals.
Figure 3.4: Frequencies of maximum response. Shown here by ’×’ are the frequencies of maximum response of each local minimum of figure 3.2. Theses frequencies are located within bands of width π [rad/time-unit]. Additional information is given in table 3.2.
3.5. CARMA ESTIMATION ALGORITHM
79
Figure 3.5: The proposed estimation algorithm. Complementary information is
given in figure 3.6 and in the appendix C.
3.5 CARMA estimation algorithm
The proposed algorithm for estimation of CARMA coefficients exploits the multi-band minimization of the log-likelihood function
parametrized by θ. The vector of parameters θ, which consists of
80
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
the zeros and the poles of (3.17), can also be represented by the
polynomial coefficients of its numerator and of its denominator.
These coefficients are real numbers and the numerical optimization was carried out by using this type of parameterization. For example, the power spectrum of an AR(2) process having two poles
s1 , s2 is given by
Φ(s; θ) =
σ2
,
(s 2 − a1 s + a0 )(s 2 + a1 s + a0 )
(3.58)
where a0 = s1 · s2 and a1 = s1 + s2 . The proposed estimation algorithm aims at finding the coefficients a0 , a1 and the variance of the
innovation process σ2 . The dominant coefficient in this example,
i.e. the one that changes substantially as the values of the poles
change, is a0 . It has the most prominent effect on the geometric distance between the coefficients one starts with (i.e. the initial conditions for extracting every local minimum) and the coefficients of the local minimum one ends up with. In cases where s1 , s2
are complex conjugate, the coefficient a0 corresponds to their module value as demonstrated in figure 3.2. Every location in the image corresponds to a different set of initial conditions. This set was
then used for minimizing the proposed approximation of the likelihood function and the color at that point indicates the value it
converged to upon minimization. Initial conditions were determined by choosing various complex values for the pole s1 and the
radial appearance of the image stems from its module. For the
general case of an ARMA(p,q)
process, the parameters that
©
ª are being optimized are θ = σ2 , a0 , a1 , . . . , a p−1 , b 0 , b 1 , . . . , b q−1 , which
gives rise to the following power spectrum
£ q Pq−1
¤£
¤
Pq−1
s + n=0 b n s n (−s)q + n=0 b n (−s)n
2
Φ(s; θ) = σ £
¤£
¤ . (3.59)
Pp−1
Pp−1
s p + m=0 am s m (−s)p + m=0 am (−s)m
The proposed algorithm finds K local minima of the approximated likelihood function at consecutive frequency bands and
chooses the minimum value among them. Obtaining K different
local minima requires K sets of initial conditions. It is suggested
here to choose bandpass power spectra that have peak responses
at ω0 + πk [rad/time-unit] where k is the band index and where ω0
is determined by fitting the available data with a discrete-domain
AR(p) process and by extracting the frequency that corresponds to
the maximum value of the power spectrum. For arbitrary sampling
interval values, one should use πk
T instead of πk. While there are
3.5. CARMA ESTIMATION ALGORITHM
81
Figure 3.6: A detailed description of the log-likelihood function l˜(θ;x) calculation. This workflow is part of the estimation algorithm of figure 3.5.
82
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
many ways of obtaining such a band pass spectrum, we provide a
constructive way of doing so in appendix C. Initial values for σ can
then be obtained by evaluating (3.19) at t = 0 and equating it to
the variance of the available data. The value of K may be derived
from known physical constraints of the problem at hand. When
there is no such knowledge, this value can be determined during
the execution of the estimation algorithm: following table 3.2, likelihood values of local minima exhibit monotonicity; they decrease
towards the global minimum and then start to monotonically increase as a function of the band index k. Such monotonicity can
then provide a decision rule for setting a value for K . It is noted that
this value does not depend on the size of the data or on the ARMA
order. It is related to the sampling rate value. A single sampling
interval value partitions the frequency axis into non-overlapping
segments which are π/T wide, and, for different sampling interval values, a given continuous domain model can be associated
with different segments. Associating a segment with a power spectrum amounts to allocating the peak response frequency to that
segment. Every local minimum corresponds to a single segment
and the choice of K , which is the number of segment the algorithm
examines, should take this partitioning into account. The algorithm is described in figure 3.5 while complementary information
is given in figures 3.6 and C.1.
The work of [77] suggests three algorithms for computing the
continuous-time counterpart of a known discrete-domain model.
These algorithms can be utilized for indirect estimation
approaches, for which the sampled data is used for estimating the
discrete-domain model first. Nevertheless, [77] points out two safeguards: existence of a continuous domain model is not guaranteed
for all discrete-domain models, and uniqueness of the continuous
domain poles is not guaranteed, either. If one, however, restricts
the imaginary part of the continuous domain poles to be in the interval (− Tπ , Tπ ), as suggested in [77], the poles can then be uniquely
resolved. This uniqueness condition imposes an upper bound on
the sampling interval T .
In our work, the discrete-domain model of Theorem 1 is always guaranteed to have a continuous domain counterpart, as it
originates from the sampled version of the autocorrelation function. Also, our approach to the uniqueness problem is different
than in [77]. Instead of imposing restrictions on sampling interval
values, we are allowing for arbitrary values to be considered while
exploiting the local minima properties of the likelihood function.
3.5. CARMA ESTIMATION ALGORITHM
83
In particular, we find several sets of estimated parameters, each
corresponding to a local minimum, and choose the one that corresponds to the global minimum rather than the one that resides in
the baseband (− Tπ , Tπ ). The work of [88] introduces two estimation
algorithms. The first algorithm relies on the results [77] and it consists of two steps: 1) identifying parameters of a discrete-domain
ARMA process by minimizing the variance of the prediction error,
and 2) associating the discrete-domain parameters with a continuous domain ARMA process by means state space representation. Similar to the comparison we made with [77], when fitting a
discrete-domain ARMA model to sampled data one does not take
the zero-pole coupling of the discrete-model model into account.
Such a fitting may result then in ARMA models that have no continuous domain counterpart. Uniqueness of the continuous domain poles is not guaranteed either. Furthermore, the maximumlikelihood estimator in the sampled case does not correspond to
minimizing the variance of the the prediction error any more. The
second algorithm in [88] consists of two steps: 1) replacing the differentiation operator by a difference operation and approximating the continuous domain innovation by a discrete-domain white
noise with a variance that is proportional to the sampling rate, 2)
optimizing the discrete-domain model to minimize the variance
of the discrete-domain innovation process. The discrete-domain
process that results in using the difference operator of [88] provides an approximated model to the sampled process. In particular, the autocorrelation sequence of this model does not necessarily correspond to the autocorrelation values of the continuous
domain model one started with. Similar to the [88], our approach
to discretizing the continuous domain model is simple too, as it
involves partial fraction decomposition and finite weighted sums
only. Nevertheless, our spline formulation does not depend on
sampling interval values and it provides a discrete-domain process whose autocorrelation sequence is in a perfect match with the
point-wise values of the autocorrelation function for any sampling
interval value.
The work of [89] minimizes a likelihood function for parameters estimation. This function is expressed in terms of the discretetime innovation process, allowing one to write the likelihood function as a product of independent Gaussian variables. The log likelihood function of [89] is then given by a sum. The innovation process corresponds to prediction error values and the work in [89] relies on a linear predictor. This predictor is recursively defined and
84
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
depends on point-wise evaluations of the autocorrelation function. The likelihood formulation we propose in our work is different. We take advantage of the fact that the samples are taken at
fixed intervals and find a digital whitening filter that can be directly
applied to the sampled data rather than using a recursion. Also,
our suggested whitening filter is given in the z-domain, allowing
one to use a Direct Form realization rather than a discrete-time
convolution operation which is computationally more demanding. The recursion of [89] involves weighted sums of the sampled
data, which is even more computationally demanding. Both approaches minimize the likelihood function by means of numerical
optimization. Our formulation, however, allows one to consider
large data sets of more than 100.000 sample values compared to
less than 1.000 that are presented in [89]. Additionally, increasing
the number of data points in our approach has no effect on the
digital filter, whereas the approach of [89] requires one to calculate new additional weights. Another difference is the point-wise
evaluation of the autocorrelation function. Our work relies on partial fraction decomposition of the Laplace transform (3.18). The
time-domain expression for the autocorrelation function is then
analytically derived, as shown in table 3.1. The work of [89] on the
other hand, relies on a state-space representation of the continuous domain process. The point-wise evaluation requires a matrix
exponential operation and several integration operations for every
set of parameters. Such a method is computationally expensive
when using numerical optimization. Furthermore, according to
[90], these matrix operations are approximated by Riemann sums
rather than defined by closed form expressions. The reference [90]
is pointed out in [89]. Our likelihood formulation requires 2p + 1
evaluations of the autocorrelation function and there is no dependency on the size of the data, N . This stems from the fact that we
characterize the whitening filter in the z-domain rather than by its
time-domain values. The formulation of [89], on the other hand,
requires one to evaluate the autocorrelation function at N distinct
points.
The work of [91] relies on the Gaussian property of the sampled
process and the likelihood function in [91] is expressed in terms of
the autocorrelation matrix. The approximation that is suggested in
[91] has the following properties: following the profile likelihood
method, the variance of the innovation process is not included in
the likelihood function formulation of [91]. This amounts to reducing the number of the estimated parameters. The innovation
3.5. CARMA ESTIMATION ALGORITHM
85
variance is estimated instead by the sample variance of the data in
[91]. In our work, the innovation variance appears in the likelihood
formulation and we use the sample variance for initializing the numerical optimization only, as described in figure C.1 of appendix C.
The continuous domain power spectrum model of [91] is assumed
to have distinct roots. Our formulation does not include such an
assumption. Power spectra of multiple poles are likely to arise during the numerical optimization process, and our approach is better designed to cope with such instances of power spectra. The
likelihood function of [91] involves matrix inversion, which can
be calculated by means of Cholesky factorization, as described in
[91]. When the data size is too large, it is suggested in [91] to use
a nearest-neighbor method, which corresponds to de-correlating
the innovation process with respect to neighbouring sample values only. This means that for a fixed sampling interval value, one
has to find the inverse of a relatively small matrix and repeatedly
multiply it by small segments of data. Our work, on the other hand,
uses a digital filter for de-correlating all of the data at once. The
underlying assumption in our work is that the data size is very
large, allowing one to approximate the inverted covariance matrix
by a Toeplitz matrix. Theorem 2 guarantees convergence of the
latter to the former with increasing number of samples. Furthermore, our digital filtering approach is computationally less complex than the nearest neighbour method, as it does not require repeated matrix multiplications. The proposed method takes advantage of the existence of local minima, rather than considering it as
a difficulty. Our algorithm is designed to find several local minima
and to choose the minimum value among them, as demonstrated
in figure 3.2 and table 3.2. These local minima stem from aliasing and they can be found by properly allocating the initial conditions of the numerical optimization. The proposed algorithm
minimizes the likelihood function using several sets of initial conditions and chooses the solution that corresponds to the minimum
value among all the local minima that were found.
Numerical optimization was also used in [76, 89, 91] and our
algorithm repeats the optimization procedures several times. The
number of repetitions, denoted K in figure 3.5, is determined by
the number of the local minima one wishes to find and it is a userdependent parameter. Although our algorithm uses several minimizations of the likelihood function, our spline formulation allows
for fast implementation as it involves the simple parameter mapping of Theorem 1 and a Direct Form implementation of a digital
86
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
filter.
3.6 Simulation results
The proposed approach was implemented in Matlab™ using nonconstrained numerical optimization. It was compared with the
polynomial B-spline estimator of [76]. Sampled signals were generated by filtering a discrete-domain Gaussian white noise with
the digital filter (3.25). Several sets of the parameters θ were considered. For each set, several Monte Carlo simulations were conducted; each simulation corresponds to a different sampling interval. Sampling interval values were chosen so as to capture different aliasing constellations. A single Monte Carlo simulation involved 500 experiments and every experiment was carried out using N = 10, 000 sample values. The variance parameter was set to
σ2 = 1 and was unknown to the estimation algorithm. The value
of κ(θ) (3.43) was calculated using the first 500 terms in the infinite sum and the number of local minima that were examined was
K = 2. The estimation error for a Monte Carlo simulation is the relative MSE (Mean Square Error) between the estimated parameter
and the correct one. For example, the estimation error of a single
parameter, e.g. a0 , is given by
!
à 1 P500
2
500 n=1 ( â 0,n − a 0 )
,
(3.60)
ǫ(a0 ) = −10log
(a0 )2
where â0,n is the estimation of a0 at the n-th experiment.
Experimental results for an AR(2) process are shown in figure
3.7 and in table 3.3. In figure 3.7, Monte Carlo simulation results
for a continuous domain AR(2) process are shown. Sampling interval values were chosen so as to capture different aliasing constellations. The proposed exponential-based approach was compared
with the polynomial model of [76]. The poles of the process are
s1,2 = −0.2 ± 7i and the variance is σ2 = 1. The number of samples
is N = 10, 000. The × mark that appears on the x-axis corresponds
to the sampling interval Tmax for which w max = π/Tmax is the frequency that corresponds to the value of the continuous domain
power spectrum at −10 [dB] respect to the peak value. The CRB is
depicted for comparison purposes while other lines are depicted
for presentation purposes.
Let us consider a generic unit measuring time or spatial extension of a simulated CARMA process: [time-unit]. The maxi-
3.6. S IMULATION RESULTS
87
mum value of the continuous power spectrum of the AR model
of figure 3.7 is located at 5.05 [rad/time-unit] and the 3 [dB] frequencies are 4.09 and 6.10 [rad/time-unit]. A sampling interval
value of T = 0.21 [time-unit] corresponds to a sampling frequency
of 29.92 [rad/time-unit] and it introduces minor aliasing effects.
The value of T = 1.24 [time-unit] corresponds to a sampling frequency of 5.07 [rad/time-unit] and it introduces substantial aliasing effects. The various aliasing configurations are depicted in figure 3.8
The CRB values were calculated according to [92] by using
N = 1000 and approximating its down shift for higher dimensional
signals. Because of the approximated CRB computations and the
limited number of Monte Carlo simulations, the analysed estimators may provide errors lower than the CRB. Figure 3.9 also describes the proposed exponential-based approach from a frequency
domain point of view, emphasizing the flexibility of the exponential B-spline model to adapt to band-pass power spectra. Shown
here is a periodogram of sample values (absolute square value of
the DFT), periodized in the range [0, 4π]. The power spectrum of
the AR(2) process is depicted as a solid red line. While the polynomial based estimator captures a base-band signal (black dashed
line), the proposed approach adequately identifies the proper frequency band (magenta dashed-dotted line). The proposed windowing function corresponds to the Fourier transform of the interpolator η(t ; θ). The poles of the process are s1,2 = −1 ± 10i and the
sampling interval is of a unit value, introducing prominent aliasing
conditions.
An ARMA(2,1) estimation comparison is given in figure 3.10
and in table 3.4. We note that more sample values are required for
the proposed estimation algorithm when the number of parameters increases. Our results indicate that the proposed approach
outperforms the polynomial-based direct method while following
the CRB at various aliasing constellations. It also guarantees that
the estimated parameters correspond to a valid continuous domain model whereas this property does not necessarily hold true
for other discrete-domain-based methods, such as the minimization of the prediction error variance. Unlike the proposed approach,
the prediction error method might also result into a continuous
domain power spectrum that occupies lower frequencies as it does
not take aliasing into account.
88
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
Table 3.3: Estimation comparison for AR(2) processes.
Estimation Error [dB]
Power
Spectrum
Sampling
Interval
[time-unit]
[76]
Proposed
[76]
Proposed
[76]
Proposed
1
0.3370
-19.76
-30.98
-42.03
-55.90
-31.89
-39.60
0.6827
-9.40
-33.25
-1.09
-58.61
-1.46
-39.20
1.1150
-13.42
-33.39
-0.40
-59.63
-0.48
-37.81
0.6059
-34.93
-34.97
-34.94
-45.89
-27.92
-36.99
1.3050
-16.68
-35.79
-4.75
-43.77
-3.73
-35.81
2.2370
6.93
-32.64
-1.86
-35.63
-8.16
-33.92
0.9753
-21.09
-22.38
-17.95
-23.70
-11.73
-17.05
2.3507
-3.65
-4.17
-3.44
-3.68
-1.19
-1.23
3.6011
-1.83
-2.28
-1.31
-1.60
-0.27
-0.34
s 2 +s+9.25
¡
ǫ(σ2 )
ǫ(a1 )
Φ(s; θ)
(s 2 +0.4s+49.04)(s 2 −0.4s+49.04)
¡
¢
s1,2 = −0.2 ± 7i
(
ǫ(a0 )
1
)(
s 2 −s+9.25
s1,2 = −0.5 ± 3i
)
¢
1
(s 2 +6s+5)(s 2 −6s+5)
¡
¢
s1,2 = −1, −5
Table 3.4: Estimation comparison for ARMA(2,1) processes.
Power
Spectrum
Estimation Error [dB]
Sampling
Interval
ǫ(a0 )
ǫ(a1 )
ǫ(b 0 )
ǫ(σ2 )
Φ(s; θ)
[time-unit]
[76] Proposed
[76] Proposed
[76] Proposed
[76] Proposed
(s−3)(s+3)
(s 2 +2s+26)(s 2 −2s+26)
0.3009
-18.52
-45.07
-32.18
-51.37
-12.76
-39.02
-23.37
0.6481
-8.44
-36.85
-18.96
-33.29
-1.65
-17.42
-10.43
-22.50
1.1111
-11.60
-29.69
-0.05
-41.59
-0.26
-12.61
-11.83
-15.16
0.1911
-19.51
-33.71
-32.18
-51.37
-12.76
-39.02
-23.37
-42.17
0.4115
-14.48
-34.78
-18.96
-33.29
-1.65
-17.42
-10.43
-22.50
0.7055
-19.96
-29.00
-0.05
-41.59
-0.26
-12.61
-11.83
-15.16
0.3250
-16.42
-21.08
-18.24
-17.26
-10.73
-20.36
1.30
-24.55
1.3412
-4.63
-5.77
-3.49
-4.44
-4.56
-7.66
-15.31
-7.33
1.8831
-2.97
-3.96
-2.00
-2.51
-3.65
-4.73
-9.11
-5.76
¡
¢
r 1 = −3, s1,2 = −1 ± 5i
(s−4)(s+4)
(s 2 +2s+101)(s 2 −2s+101)
¡
¢
r 1 = −4, s1,2 = −1 ± 10i
(s−3)(s+3)
(s 2 +8s+7)(s 2 −8s+7)
¡
¢
r 1 = −3, s1,2 = −1, −7i
-42.17
3.6. S IMULATION RESULTS
89
50
Proposed
Polynomial
CRB
40
Relative MSE [dB]
30
20
10
0
−10
−20
−30
−40
0.21
X
0.36
0.50
0.65
0.80
0.95
1.10
1.24
0.95
1.10
1.24
0.95
1.10
1.24
Sampling Interval
(a)
0
Proposed
Polynomial
CRB
−10
Relative MSE [dB]
−20
−30
−40
−50
−60
−70
0.21
X
0.36
0.50
0.65
0.80
Sampling Interval
(b)
30
Proposed
Polynomial
CRB
20
Relative MSE [dB]
10
0
−10
−20
−30
−40
−50
0.21
X
0.36
0.50
0.65
0.80
Sampling Interval
(c)
Figure 3.7: Estimation error comparison for an AR(2) process. Figures 3.7(a),
3.7(b), 3.7(c) show the estimation error for parameters a0 , a1 and σ respectively.
90
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
Figure 3.8: Power spectra of sampled processes. Shown here are power spectra of
discrete-domain processes that originate from the same continuous domain process. Every power spectrum corresponds to a different sampling interval value. The
poles of the continuous domain AR process are s 1,2 = −0.2 ± 7i
1
Periodogram
Power−spectrum
Proposed windowing
Polynomial windowing
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
2
4
6
ω [rad/time−unit]
8
10
12
Figure 3.9: A frequency domain description of the proposed estimation approach.
3.6. S IMULATION RESULTS
91
(a)
(b)
(c)
(d)
Figure 3.10: Similar to figure 3.7 for an ARMA(2,1) process. The poles of the
process are s 1,2 = −1±5i , the zero is r 1 = −3 and the variance is σ2 = 1. The number
of samples is N = 100,000. Figures 3.10(a), 3.10(b), 3.10(c) and 3.10(d) show the
estimation error for parameters a0 , b 0 , a1 and σ respectively.
92
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
3.7 CARMA for ultrasound tissue characterization
In ultrasound diagnostic, the acquired image is a rich source of information pertaining to the imaged tissue, even though some of it
is not directly accessible at human visual inspection. To reveal this
hidden information it is necessary to resort to computer-aided tissue characterization. This technique can capture several features
carried by medical sonography and correlate them to the state of
the imaged tissue. In particular, the statistical modeling of discrete
ultrasound signals has been widely employed to derive interesting properties from estimated model parameters that are useful in
the context of tissue characterization. Among them, shift-variant
ARMA processes were found suitable to model ultrasound signals
and images and to make the acquisition more robust with respect
to the estimation of the uncorrupted tissue response [6, 93]. More
recently, 2D ARMA modeling was proposed to improve computeraided detection of breast tumors [94, 95].
Although CARMA stochastic processes have been successfully
exploited for image modeling [78, 61], the continuous domain approach has never been investigated in the processing of ultrasound
images. In this study, we consider the radio-frequency (RF) ultrasound signal as a sampled CARMA process and recover its CARMA
parameters from the traditional ARMA coefficients through an indirect approach.
The ARMA coefficients are easily computed by minimizing the
prediction error [55]. Then, the CARMA parameters can be estimated from their discrete-domain counterpart by mapping the
ARMA model back into the continuous domain [67]. In doing so,
both the exact discretization of CARMA model and the exponential properties of the autocorrelation function in the discrete and
continuous domains must be taken into account. As was shown
in [86], the exact link between ARMA and its corresponding CARMA
model can be established, provided one is willing to interpolate
the discrete signal by exponential B-splines specific to the considered model [56, 57]. In situations where a sufficiently high sampling frequency is employed, it is beneficial to exploit the dependence of the exponential B-spline on the ARMA parameters of the
sampled model. This, along with a numerical optimization scheme,
allows for the accurate estimation of the continuous domain model
parameters.
Here, we analyze the CARMA parameters of ultrasound signals
in terms of their ability to capture information about the concen-
3.7. CARMA FOR ULTRASOUND TISSUE CHARACTERIZATION
93
tration of scatterers in imaged materials. For this purpose, we have
created ultrasound images of phantoms with various concentrations of scatterers. The validation of CARMA and ARMA models in
terms of discrimination performance shows that our continuousbased approach is more efficient than the traditional one. We believe that our unconventional ultrasound modeling is robust with
respect to the effects of the system and has a good potential for
tissue characterization.
We start by assuming that the spatial dependence of the effects
of the acquisition system is weak. This allows us to focus on a segment of the ultrasound image and to consider its characteristics
as locally shift-invariant [41]. As a consequence, the 1D discretetime ultrasound signal x[n] can be modeled as a spatially invariant
sampled CARMA process x(t ), where the continuous-time Gaussian white noise w(t ) with variance σ2 represents the innovation
signal.
For a CARMA model of order (n, m), the spectrum of x(t ) is expressed as
Φ( j ω; θ)
=
=
¯2
¯ Qq
¯
j ω − r k ¯¯
k=1
σ ¯ Qp
¯
¯
j ω − sk ¯
k=1
¯ Pq
¡ ¢k ¯2
¯
b j ω ¯¯
2 ¯ k=0 k
σ ¯ Pp
¡ ¢k ¯ ,
¯
a jω ¯
2¯
k=0
(3.61)
(3.62)
k
where θ = (s, r ) = (s1 , . . . , s p , r 1 , . . . , r q ) is the collection of parameters that contains the poles sk and the zeroes r k of the system,
while ak and b k are the related polynomial coefficients.
In corollary equation (3.25) of section 3.2, it was shown that
sampling the model (3.62) leads to a discrete ARMA model of order
(p, p − 1). For the sake of simplicity, in the sequel we’ll consider a
unit-sampling interval. The discrete spectrum is then given by
Φd (z; θ) = σ2d Hd (z; θ) Hd (z −1 ; θ),
(3.63)
where
Hd (z; θ) =
Qp−1 ¡
k=1
Qp
k=1
¡
1 − νk z −1
1 − ρk z
¢
Pp−1
k=0
¢ = Pp
−1
dk z −k
c
k=0 k
z −k
.
(3.64)
The sampled ultrasound signal, which represents the available data,
is thus modeled as a discrete ARMA (p, p − 1). In this context the
94
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
autocorrelation function of (3.62) is derived from its discrete samples thanks to interpolation by shifts of the exponential B-spline
function β(t ; θ) expressed in equation (3.21) and reported here for
clarity:
β(t;θ) = F −1
(
q
Y
( j ω − r k )(−j ω − r k ) ·
k=1
)
p
Y
(1 − e j ω+s k )(1 − e − j ω+s k )
(t)
( j ω − s k )(−j ω − s k )
k=1
(3.65)
For this reason, the zeroes of the ARMA filter in (3.64) are the roots
of the z-transform of the exponential B-spline samples Bd (z; θ) =
P
−k
that are located on the unit circle:
k∈Z β[k; θ] z
Bd (z; θ) = σ2d
p−1
X
k=0
dk (θ) z −k
p−1
X
dk (θ) z k .
(3.66)
k=0
These expressions present a useful link between the discrete- and
continuous domain ARMA models, and describe the mapping required to extract CARMA parameters from ARMA coefficients.
When the sampling frequency is high enough, the estimation
of the CARMA parameters through indirect methods is efficient [67],
and there should be no need to consider approaches focusing on
several bands as presented in section 3.5. Fortunately, high frequency sampling operations are prevalent in ultrasound diagnostic. Thus, the derivation of a CARMA model by the inverse mapping of ARMA parameters offers good estimation performances.
We assume here to know the ARMA model orders, which can be
estimated as the orders providing a good fitting of the model to the
signal power spectrum. In this selection it is necessary to respect
the constraint that the relative order of the model must be 1 for the
signal to be a sampled CARMA process. The CARMA model order n
is constrained to the ARMA orders, while order m < n can be fixed
to provide an equal number of parameters for CARMA and ARMA
models.
The ARMA parameters (c k , dk ) and (ρ k , νk ) can be computed
by the standard minimization of the prediction error [55], which
provides the estimated values (ĉ k , dˆk ) and (ρ̂ k , ν̂k ). While the link
between the continuous domain poles sk and the discrete-domain
poles ρ k is a trivial exponential mapping (in high frequency sampling conditions), the expression of the discrete-domain zeroes νk
is a complicated function of both continuous domain poles and
zeroes, as can be seen from (3.66).
Since Bd (z; θ) depends solely on the CARMA parameters, it is
possible to estimate the zeros r k by the numerical minimization
3.8. C OMPARATIVE STUDY ON PHANTOM IMAGES
95
Algorithm 3.1 Estimation of CARMA parameters
• Discrete-domain estimation of parameters
Estimate (ĉ , dˆ ) and (ρ̂ , ν̂ ) by minimizing the prediction error from x[n]
k
k
Set σ̂2d = VAR{x[n]}
k
k
• continuous domain estimation of parameters
Set ŝ k = log(ρ̂ k ) ⇒ ŝ
w2
w
Set r̂ = arg minr wd̂ − d(ŝ,r )w
´−1
³
Qp
Q p−1
ŝ k
(3.34)
Set σ̂2 = σ̂2d · β[1;θ] · k=1 −1
ν̂ · k=1 e
k
of the distance between the estimated ARMA coefficients dˆk and
the function dk (θ). The latter can be expressed as dk (s, r ), which
highlights the separation between the collection of poles s and the
collection of zeroes r . The parameters dk (s, r ) can be built from
the coefficients of the polynomial that corresponds to the roots of
Bd (z; θ) in the unit circle. The collection of coefficients is indicated by d(s, r). Likewise, the collection of estimated coefficients
of the ARMA model numerator is d̂.
Setting σ̂2d as the variance of the RF signal samples x[n], we can
estimate σ2 from the link between the variances of the discrete and
continuous innovation signal (3.34).
The indirect method for CARMA model estimation is summarized in Algorithm 3.1. We used this algorithm to compute ARMA
and CARMA parameters from a sequence of phantom images with
various concentrations of scatterers. The two models were compared in terms of classification performance, as described in the
next section.
3.8 Comparative study on phantom images
The ultrasound images were realized using as scatterers ultrafine
polyamide particles (Orgasol© , Arkema, France) of diameter
10 ± 2µm and density 1030kg/m3 . The tissue-mimicking phantoms were prepared by mixing a specific concentration of Orgasol©
particles with distilled water. A magnetic agitator preserved the
solution homogeneous throughout the acquisition. In our experiments, we used five different mixtures having Orgasol© concentrations of 0.5%, 0.75%, 3%, 6%, and 12%. Low and high concentrations mimic random and dense mediums respectively.
The RF signal was acquired by a high-resolution ultrasound
scanner (Vevo 770, Visualsonics, Toronto, Canada) at a central fre-
96
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
quency of 30 MHz and at a sampling frequency of 500 MHz with
8-bit resolution. For each phantom three video acquisitions were
made. A selection of 90 scanlines from the central frame of each
video sequence was used (about 300 samples per concentration).
Along each scanline, the RF signal was used to compute the ARMA
and CARMA parameters, as explained in section 3.7, except that
we found it is more convenient to numerically optimize real polynomial coefficients, as opposed to complex poles and zeros.
In order to ensure a good fit of the available data, the orders of
the CARMA model were chosen as m = 2 and n = 4. Without loss
of generality, we constrain the first polynomial coefficients to take
a unit value, so that a0 = b 0 = c 0 = d0 = 1. The process is thus characterized by the four pole parameters a = (a1 , a2 , a3 , a4 ), two zero
coefficients, and the variance of the innovation signal, which we
represent by the collection b = (b 1 , b 2 , σ2 ). The orders of the ARMA
model are constrained because of the sampling process performed
on the continuous-time model: the number n = 4 of poles is the
same, while there are n − 1 = 3 zeroes. The ARMA model parameters are represented by the collections c = (c 1 , c 2 , c 3 , c 4 ) and d =
(d1 , d2 , d3 ). Finally, the full set of parameters is a seven-attribute
vector for both models.
The Mahalanobis distance along a single feature provides information about the individual discrimination power of the considered feature for distinguishing between two concentrations of
scatterers [8]. We show in table 3.5 the Mahalanobis distance projected onto the direction corresponding to the coefficients on the
rows belonging to the two analyzed models. We chose the most
significant coefficients, separately for each model. The columns
of the table represent the binary classification task of two consecutive concentrations of scatterers, where J [x,y ] is the Mahalanobis
distance between the sample clusters with concentration x and y.
From this first analysis, it is possible to observe that i) the higher
the concentration of the scatterers, the more difficult the discrimination by a single coefficient is; ii) the CARMA coefficients clearly
provide higher distances, which eases the classification task.
To extend these results to the classification of joint parameters, we performed several multiclass discrimination experiments
where all available parameters were involved, for both models. The
learning scheme used was a support-vector machine with Gaussian kernel, whose parameter had been tuned by cross-validation,
separately for each classification model. We report in table 3.6 the
rate of misclassification (averaged on 20 experiments) for the four
3.8. C OMPARATIVE STUDY ON PHANTOM IMAGES
97
Table 3.5: Mahalanobis distance for various concentrations.
Model
CARMA
ARMA
a4
b2
c4
d3
J [.5%,.75%]
14.01
11.30
0.56
2.73
J [.75%,3%]
10.55
39.18
5.22
12.02
J [3%,6%]
2.13
5 · 10−4
0.33
0.19
J [6%,12%]
4.44
0.17
0.37
0.08
multiclass discrimination problems listed in the first column. The
classifier was first trained on CARMA and ARMA pole coefficients
a and c. Then, it was trained on the zero coefficients and variances
b and d of both models. Finally, the learning was performed on all
the available coefficients for each model.
The classification results show that the proposed CARMA coefficients are more sensitive to differences in the concentration of
scatterers, above all when the discrimination is multi-class. continuous domain parameters favour a better recognition score and
represent tissues of different concentrations in more compact and
better discernible clusters, as illustrated in figures 3.11 and 3.12
for the discrimination problem with 5 categories of scatterers concentration. In these figures linear discriminant analysis [8] is used
to reduce data dimensionality and to find the number Ncl asses −
1 dimensional space, where mapped samples present maximum
inter-class distance and minimum intra-class variance. The sevendimensional samples are, thus, projected onto a four dimensional
space and normalized; two planes of the projection space are visualized in figures 3.11 and 3.12. The bivariate distribution of
samples in the projection planes is described by a 2D Gaussian fitting for each category, whose ellipse at a probability density function of 80% is shown in the figures. We observe that the clusters
are better defined in the space of CARMA coefficients than in the
space of the traditional ARMA parameters. Furthermore, as regards higher concentrations [3%, 6%, 12%] samples, although the
discrimination task is more difficult, a good classification of these
categories is provided by CARMA parameters analysis considering
all the projection directions, while ARMA coefficients do not bring
more separation, as shown in figure 3.12.
The difference of performances between CARMA-based and
ARMA-based classifiers of the last two columns in Table 3.6 is also
verified by resampled paired t-tests. These statistic tests, perfor-
98
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
Table 3.6: Rate of misclassification: [mean value ± standard deviation] of 20 experiments (1 CARMA model, 2 ARMA model).
Classes
Concentration
[0.5%, 0.75%]
[0.5%, 0.75%, 3%]
[0.5%, 0.75%, 3%, 6%]
[0.5%, 0.75%, 3%, 6%, 12%]
[0.5%, 0.75%]
[0.5%, 0.75%, 3%]
[0.5%, 0.75%, 3%, 6%]
[0.5%, 0.75%, 3%, 6%, 12%]
[0.5%, 0.75%]
[0.5%, 0.75%, 3%]
[0.5%, 0.75%, 3%, 6%]
[0.5%, 0.75%, 3%, 6%, 12%]
Features for classification
poles
a1
c2
0.007 ± 0.005
0.024 ± 0.006
0.010 ± 0.003
0.024 ± 0.004
0.116 ± 0.006
0.127 ± 0.006
0.126 ± 0.007
0.135 ± 0.006
zeroes
b1
d2
0.013 ± 0.005
0.008 ± 0.003
0.008 ± 0.003
0.032 ± 0.006
0.156 ± 0.007
0.270 ± 0.008
0.273 ± 0.009
0.368 ± 0.010
all parameters
1
(a, b)
(c, d)2
0.007 ± 0.003
0.010 ± 0.003
0.003 ± 0.002
0.012 ± 0.003
0.069 ± 0.006
0.113 ± 0.006
0.085 ± 0.006
0.122 ± 0.006
med on the misclassification errors obtained in 20 experiments,
indicate in all the four classification problems an increment in performance of CARMA-based classifier respect to ARMA-based classifier at a 5% significance level. The improvement of performance
is significant and encourages the application of continuous-time
models on in vivo images for the characterization of ultrasound
tissues.
3.9. C ONCLUSION
(a)
99
(b)
Figure 3.11: Scatter plots of normalized CARMA parameters in figure (a) and
ARMA coefficients in figure (b) for several scatterer concentrations. The plane visualized here is given by the first two directions selected by linear discriminant analysis.
(a)
(b)
Figure 3.12: Scatter plots of normalized CARMA parameters in figure (a) and
ARMA coefficients in figure (b) for concentrations [3%,6%,12%]. The plane visualized here is given by the last two directions selected by linear discriminant analysis.
3.9 Conclusion
In this chapter we introduced a new formulation of CARMA models based on exponential splines. We further proposed a maximum likelihood estimator for continuous domain ARMA parameters from sampled data. It utilizes the exponential B-spline framework while introducing an exact zero-pole coupling for the sampled process. The relation between the autocorrelation function
and the autocorrelation sequence of the sampled process was investigated in both time- and frequency domains. Our approach
100
Three. A CARMA MODEL FOR ULTRASOUND SIGNALS
relies on the fact that the Cramér-Rao bound can be made arbitrarily small by increasing the number of samples while fixing the
sampling interval at an arbitrary value. The likelihood function
of the sampled process was investigated too. In particular, it was
shown that it possesses local minima that originate from aliasing.
The global minimum, however, corresponds to the maximum likelihood value.
The only assumptions that were made throughout this work
are that the number of available samples is relatively large and that
the model orders are known, allowing one to replace the whitening matrix by a digital filter. This approximation was then shown to
be valid when considering expected values of the likelihood function. Experimental results indicate that the proposed exponentialbased approach closely follows the Cramér-Rao bound for various
aliasing configurations, while imposing no restrictions on sampling
rate values.
Additionally we presented a continuous-time ARMA model for
RF ultrasound signals and proposed an indirect approach to compute CARMA parameters based on the properties of exponential
B-splines. This work includes a comparative study of the proposed
CARMA coefficients and of the traditional ARMA parameters that
are widely used for the analysis of ultrasound images. We performed experiments on phantom images specifically designed with
various concentrations of scatterers and we focused on the comparison of CARMA and ARMA features in terms of discrimination
performances.
Experimental results reveal that the new CARMA features are
able to capture the information about the concentration of scatterers more efficiently. Thus, they can be extremely useful for the
characterization of in vivo ultrasound tissues.
Four
I MPROVING
STANDARD
PROSTATE BIOPSY PROTOCOL
It is a very sad thing that nowadays there is so
little useless information.
Oscar Wilde (1854-1900)
T
HE practical application of this thesis is the improvement of
standard prostate biopsy protocols guided by ultrasound, used
in clinical routine as the more reliable tool for prostate cancer diagnosis.
During a biopsy examination the urologist removes from 8 to
12 tissue samples, usually with a needle guided by some imaging
tool. Then a pathologist checks for cancer cells. As the bioptic procedure is generally performed in the doctor’s office with local anesthesia, the most preferred guiding tool is transrectal ultrasound.
The current clinical procedure for prostate cancer biopsy has
several limitations, above all due to the variable appearance of cancer on ultrasound. Because of this variability, TRUS images visual
inspection is characterized by low specificity. For this reason TRUS
guided biopsy is performed by uniformly sampling the gland without any information about the disease probability of the imaged
tissues.
Nevertheless computer-aided analysis of TRUS images can be
exploited to highlight suspicious regions, providing useful information to the physician and guiding the standard prostate biopsy
both as a positioning and detection tool.
As stated in section 1.3, standard image-guided biopsy is not
lesion-directed, but it is based on a systematic sampling of those
regions of the prostate with higher cancer incidence, as the peripheral zone. As the disease may present several foci, this standard
procedure is quite inefficient with a detection rate of 25%.
101
102
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
The uniform sampling of some regions of the prostate aims at
maximizing the efficiency of the procedure, but carries some randomness in the selection of cores which can be compensated for
by lesion-directed methods.
The computer-aided analysis of TRUS images proposed in this
work is aimed at improving the current biopsy protocol by marking the interesting regions and discarding the healthy tissues. As a
consequence, a reduction in the number of unnecessary sampled
cores will be possible by lesion-directed biopsy without losing diagnostic power.
An essential requirement for lesion-directed biopsy to be feasible is that the algorithm performing lesion detection is fast enough
to perform real-time processing. This is necessary to have US biopsy sessions as short as the ones performed today and to preserve
the main advantage of the US respect to others imaging modalities: the real time feedback.
Furthermore in order to maximally exploit the current technology already present in healthcare institutions and prevent extra
costs due to physician training or longer biopsy sessions, radiologists must be provided with a fast and user-friendly tool. Non
real-time detection tools are interesting for research purposes but
not appropriate for clinical employment.
In this chapter a tool for real-time Computer-Aided Biopsy (rtCAB) is presented and the procedure used to train it on a medical ground truth is described. rtCAB enhances TRUS video stream
with a false color overlay image, and suggests the physician where
to sample thus reducing the total number of cores.
The proposed detection tool was implemented using parallel
programming on a specific graphic card, such that the recognition
is processed at the same time as the visualization of the standard
US video sequence, generally acquired at 24-50 fps.
On the contrary, the training of the detection algorithm, which
is not included in the classification task performed during biopsy,
is implemented using Matlab™ .
4.1 Database collection
The main limits preventing the realization of an automatic detection tool in ultrasound imaging are the non standard orientation
of the images and the absence of a public database for algorithms
testing and validation.
4.1. D ATABASE COLLECTION
103
The current TRUS guided biopsy protocol overcomes the problem of free US image orientation, in fact it exploits a standard visualization of the gland which appears in a fixed region of the image
with a fixed position with respect to the biopsy needle. This provides robust clinical references for images interpretation.
In order to train the detection tool for suspicious regions recognition, a large TRUS video sequences database was collected in
collaboration with the department of Urology of S. Orsola Hospital in Bologna. A ground truth database was realized, by recording
radio-frequency US video sequences during biopsy sessions since
November 2009.
The probe used during the sessions is an Esaote™ EC123, employed with US central frequency of 7.5MHz, as in the study described in chapter 2. The dataset employed in this work contains
about 400 video sequences, one for each sampled core, belonging
to 42 patients. This ground truth is continuously updated and increased, alongside the proceeding of the research project, and so
far about 1000 US video sequences belonging to 120 patients have
been recorded.
The histopathological analysis of each core performed by an
expert physician provides information about the pathological state
of the imaged tissue, reporting, in case of disease, the percentage
of tumour in the core and the corresponding Gleason score.
Additionally the PSA value, the DRE result, the prostate volume
and the patient clinical history are likewise recorded to complete
the ground truth database. Table 4.1 shows the list of the different
typologies of tissue collected in the database and the number of
cores belonging to each category.
An example of the form used during biopsy sessions, which will
be filled with clinical information, DRE and TRUS analysis results,
is shown in figure 4.1.
104
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
Table 4.1: Ground Truth Dataset
code
# cores
M
UM
BB
BB*
BM
BM*
U
18
46
143
36
81
39
40
TOT
403
Description
Label
Cores with tumour percentage ≥ 70%
Cores with tumour percentage < 70%
Benignant prostate tissue1
Other benignant tissue1
Benignant prostate tissue2
Other benignant tissue2
Precancerous lesions:
Prostatic intraepithelial neoplasia PIN(26)
Atypical small acinar proliferation ASAP(14)
Total number of cores
42 patients (22 pathological)
Malignant
Unknown
Benignant
Benignant
Benignant
Benignant
Unknown
1 from healthy patient
2 from unhealthy patient
Although several video sequences are available for each examined patient, only a small portion of the ultrasound image is used
in this study. Since the only certain pathological information is
in the cylinder sampled by the needle, it is necessary to select the
very region covered by the needle profile in the frames preceding
the insertion.
The detection of the biopsy needle position is easily performed in the frames where the insertion is visible, due to the high
reflectivity of the needle and its fixed location. This information
is used in the frames preceding the insertion to select the core region, which is then segmented in 20 ROIs of size 101× 11. Only the
ROIs included in the core region are used in this study, since there
is no sure knowledge about the tissues not submitted to pathological analysis.
In particular automatic needle region segmentation is possible
in the scan-converted space by analysing the arithmetical difference between consecutive US frames, after being turned to greylevel images. Applying a suitable threshold to the difference image
we can detect the frame where needle insertion begins and also
highlight the needle zone.
The application of a threshold to the image given by the difference between the insertion frame and the previous frame does not
provide a perfect shape of the needle, because of the noisy appearance of US images. What appears in this image is a cluster of small
unconnected regions grouped in the needle zone.
4.1. D ATABASE COLLECTION
105
Figure 4.1: Exemple of form used during biopsy sessions to collect information
about patient history, clinical data and the results of DRE and TRUS examination
performed before and during biopsy.
106
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
(a)
(b)
(c)
Figure 4.2: Segmentation process for needle region detection. In figure 4.2(a) it
is shown the connection of small regions by Hough transform; the needle edges are
visualised in the scan-converted space in fugure 4.2(b) and in the RF space in figure
4.2(c).
These regions are then connected exploiting the Hough transform [96] and the fact that the needle is a line in the scan-converted
space. Since the needle line is always almost vertical the angle
range for connections search in the Hough transform can be constrained to [−20, 20] degrees.
Once the needle profile is detected, it is regularised by morphological operations of dilation and erosion. Finally the needle
region contour is outlined by edge detection. The found needle
profile is converted in the RF space where it appears curve and finally mapped onto frames preceding the insertion to perform ROI
selection. The segmentation process is illustrated in figure 4.2.
In spite of the histopathological analysis, there is still some uncertainty in the ground truth, due to the unknown distribution of
the tumour in the core. For this reason, only malignant cores with
a tumor covering more than 70% of the sample (category M in table 4.1) are used to train the classification system, as it will be explained in the next section. Furthermore cores diagnosed with
precancerous lesions as PIN and ASAP are ambiguous and cannot be included in the learning process as healthy or unhealthy
regions. These observations highlights the necessity of a complex
learning scheme also able to deal with unlabelled data.
4.2 Processing and learning scheme
In order to attain real-time processing capabilities, rtCAB makes
use of a twofold signal processing datapath, a schematic representation of which is depicted in figure 4.3. Along the first branch, the
incoming RF ultrasound signal undergoes envelope detection and
4.2. P ROCESSING AND LEARNING SCHEME
107
Figure 4.3: Signal processing flow graph of rtCAB: from left to right the incoming
RF signal undergoes envelope detection, statistical estimation, non-linear projection onto the feature space and linear classification. Results are scan-converted
and superimposed to the original ultrasound image.
scan conversion; this path has precedence on the second branch
as during prostate biopsy the physician is mostly concerned about
having a real-time video feedback of the organ he is sampling.
According to this requirement, a fast digital envelope detection algorithm [97] and an high speed digital scan converter [98]
based on bilinear interpolation were chosen. Algorithms were optimized and combined to seamlessly perform detection, logarithmic compression and conversion in one step. As the chosen envelope detection algorithm operates memoryless only on a few input samples, parallel processing techniques were also employed
to further speedup this step. By exploiting the computational resources made available by CUDA™ [99] enabled graphic card, a
processing time around 1 ms was obtained on a GPU with 240 cores
during the conversion of a 2792× 250 RF image to a 300× 300 pixel
US image.
The second branch is devoted to the classification task. First
a list of ROIs is automatically generated by tessellating the RF image with fixed size rectangular regions. Then, an improved version
of the envelope detection algorithm used in the first branch is applied to the incoming data. Comparing the output of this stage to
the output of a gold standard, Hilber Transform based, envelope
detector, an accuracy of −49.2dB nRMSE was measured on standard acquired US data.
Following envelope detection, Nakagami distribution [26] parameter m and scaling factor Ω, introduced in section 2.1, are estimated over each ROI. This distribution can give a good description of RF envelope statistics and is able to well capture different
scattering structures. As discussed in [26], the most discriminant
108
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
1000
2.0
Exact
Power
f(m)
100
1.5
Padé
10
m f(m)
1
0.1
0.05
1.0
0.05
0.10
0.50
m
0.10
1.00
0.50
m
1.00
5.00
5.00
10.00
10.00
Figure 4.4: Comparison of the power series approximation in the ML Nakagami estimator and the corresponding Padé Approximant: the proposed solution
is much more accurate, especially for low values of the shape parameter.
feature is the Nakagami parameter: in fact, different values of m
corresponds to different scattering types. Also the scaling factor is
important: by conveying information about the local scatter energy, its logarithm is strongly related to the log-compressed image
the physicians usually employ to detect prostate cancer.
While estimating Ω is not computational intensive, the ML estimation of m requires solving the nonlinear equation
f (m) = ψ(m)−ln(m) = ∆ where ψ(m) is the digamma function and
∆ is a quantity depending on data. As no closed form solution is
known for this equation, authors in [100] suggest to asymptotically
expand f (m) by means of a power series with strictly negative exponents. However a much better approximation can be obtained
by means of Padé Approximants [101]. In fact, substituting f (m)
with its approximant of order (1,2), a rational function in the form
p +p 1 m
fˆ(m) = q +q0 m+q
2 can be used to obtain a closed form ML esti0
1
2m
mator. Parameters p i and q i were obtained by applying numerical
optimization techniques. As shown in figure 4.4, the proposed solution is much more accurate, especially for low values of m.
The next stage is responsible for the classification of each ROI
as malignant or not malignant basing on both Nakagami distribution parameters. The main blocks of the classification scheme are
4.2. P ROCESSING AND LEARNING SCHEME
109
a non linear feature extraction stage based on Generalized Discriminant Analysis (GDA) with Gaussian kernel[49] and a linear
classification through the Fisher Linear Discriminant (FLD)[35]
with an optimized shift of the decisional threshold.
Classification of each ROI, being each other independent, is
performed exploiting CUDA™ parallel processing acceleration.
The classification of a single 101 × 11 ROI requires a mean time
of 47.6µs while 71.4µs are necessary for a larger 200 × 20 one. Results are then superimposed to the scan converted image to form
a malignancy map which helps the urologist during the sampling
procedure.
The classification model selection is based on a two-step learning scheme. The first step concerns the selection of the training set
and the variance σr b f of the GDA kernel, tuned to obtain the maximum sensitivity. The second step optimizes the linear classifier
threshold in order to provide the classification system with maximum positive predictive value.
The performance statistics considered in the first step are ROIbased, such that the sensitivity to be optimized is the number of
well classified malignant ROIs over the total number of ROIs classified as unhealthy in the validation set. The data set for this step includes only a fraction of the database shown in table 4.1. In particular, it contains the benignant ROIs coming from healthy patients
(143 BB, 2860 ROIs) and the malignant ROIs belonging to the cores
with tumor volume percentage larger than 70% (18 M, 360 ROIs).
The latter choice is necessary to make negligible the uncertainty
on the distribution of the disease in the core.
In the first step of the model selection, 500 random data sets
of 360 unhealthy and 360 healthy ROIs are generated and separated in training set (600 ROIs) and validation set (120 ROIs). For
each data set, a cross-validation over 30 values of σr b f , equally distributed between 0.1 and 6, is performed. The output of the model
selection step is constituted by the training set [X̂t r n , ĉt r n ], including the selected matrix data X̂t r n and the corresponding classes
ĉt r n , and the GDA kernel variance σ̂r b f corresponding to the maximum sensitivity in the classification of the validation set. The process for the first part of the learning scheme is shown in algorithm
4.1.
The choice of sensitivity as ROI-based criterion to select the
model is arbitrary: accuracy, specificity and PPV may be used as
well. Other criteria were also tested and the sensitivity was chosen
as the metrics providing the best core-based classification perfor-
110
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
Algorithm 4.1 Training set and GDA parameter selection
[Xtr n , ctr n ]: training set data matrix and labels
[Xv al , cv al ]: validation set data matrix and labels
for i t er = 1 to 500 do
select random [Xtr n , ctr n ] and [Xv al , cv al ]
for σr b f = 0.1 to 6 do
Train on [Xtr n , ctr n ] using σr b f as GDA kernel variance
Test on [Xv al , cv al ]
Compute ROI-based Accuracy: ACC([Xtr n , ctr n ], σr b f )
end for
end for
Output: ([X̂tr n , ĉtr n ], σ̂r b f ) = arg max ACC([Xtr n , ctr n ], σr b f )
mance as it will be shown next.
In the second step of the learning process the threshold of the
FLD is tuned according to core-based statistics in order to maximize the PPV of the classification model over the whole dataset.
The training set for this step is composed of all the cores (benignant, malignant and unknown) of all the 22 unhealthy patients for
a total of 204 cores.
Using the classification model selected in the previous step, a
cross-validation is performed on the values of the threshold varying between the minimum and the maximum value of the FLD
discriminant function. The number of thresholds considered is
1000, ranging from a model with 100% ROI-based specificity and
a model with 100% ROI-based sensitivity. For any value of the
threshold, the classification model is applied to the 204 cores and
some core-based statistics are computed for each patient. In this
case a core is considered unhealthy if at least one of its ROIs is classified as malignant.
As a consequence a core-based true positives value (TP) is the
number of malignant cores correctly classified; likewise the corebased false positives value (FP) is the number of misclassified benignant cores. As stated in section 1.2, the PPV value is given by
TP
T P +F P and represents the probability of prostate cancer positiveness of a pathological gland according to the detection tool.
The final value of the threshold is selected as the one providing the largest median PPV over all pathological patients. Whether
there are several threshold values corresponding to the maximum
PPV, the highest value is selected in order to maximize the sensitivity SE constrained to the PPV maximization. The second step of
the learning scheme is shown in algorithm 4.2.
4.3. C LASSIFICATION RESULTS
111
Algorithm 4.2 Threshold selection
t h = 0 initial threshold at zero
C classification model defined by ([X̂tr n , ĉtr n ], σ̂r b f , t h)
d : FLD discriminant function
for t h = mi n(d ) to t h = max(d ) do
for p at i ent = 1 to 22 do
Classify by C any core of this patient
Compute core-based PPV(p at i ent , t h)
Compute core-based SE(p at i ent , t h)
end for
end for
PPVm (t h) = med (PPV(p at i ent , t h)) median over patients
Set t˜h = arg max PPVm (t h)
if t˜h not unique then
tˆh = max(t˜h) corresponding to the largest SE
else
tˆh = t˜h → SE constrained to PPV
end if
Output: tˆh
4.3 Classification results
The described learning scheme selects a threshold for the classification model providing a median sensitivity of 100% and a median
PPV of 39%, as shown in figure 4.5.
The standard double sextant biopsy is equivalent to considering all the sampled cores as a priori malignant, thus it is characterized by a PPV over the studied data set of 24%. Therefore
the proposed classification scheme improves the positive predictive value of 65% without losing diagnostic power. In fact, called
p p v 0 the PPV of the TRUS guided biopsy, the probability of correctly diagnosing prostate cancer after extracting N0 cores, can be
computed as 1 − (1 − p p v 0 )N0 . Since the PPV of rtCAB in guiding
biopsies, p p v 1 , is higher than p p v 0 , the number of cores may be
log(1−ppv )
reduced keeping the same diagnostic power: N1 = N0 log(1−ppv 01 ) .
Consequently, according to the results of this study, it is possible
to achieve the same diagnostic value of a standard prostate biopsy
(8-12 cores), with no more than 7 cores in rtCAB guided biopsy.
4.4 The future
The detection tool rtCAB can be further improved by adding more
ultrasonic features to Nakagami parameters and by assuring their
real-time computation. Interesting attributes, which can improve
performance detection by keeping low time consumes, are the
Unser textural features, described in appendix A.
112
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
Positive Predicted Value (PPV = 0.39 @ Thbest )
Sensitivity (Se = 1.00 @ Th best )
1
ASA
BCA
0.9
CCE
CFR
0.8
DCE
DLO
EMA
0.7
FBE
FTA
0.6
GBI
GBO
0.5
GPR
GSI
0.4
GVO
LCO
0.3
MBA
MGU
MTU
0.2
NTO
PNA
0.1
RCO
RRO
−3
−2
−1
−3
0
a)
x 10
−2
−1
−3
0
0
b)
x 10
−3
1
0.8
median PPV
median SE
0.6
0.4
0.2
0
−6
−5
−4
−3
−2
−1
Threshold
0
1
2
3
4
x 10
−3
c)
Figure 4.5: Learning process results. figure 4.5(a) shows the core-based PPV of
the classification system. At right the color bar shows the color coding for the sensitivity and PPV values. figure 4.5(b) shows the core-based sensitivity of the classification system for any patients, listed on the y-axis with their representative code.
The x-axis represents the 1000 threshold values, while its optimal value is marked
by a vertical line. figure 4.5(c) shows the median results over all patients.
When performance results will be clinically satisfactory on a
ground truth dataset containing about 100 pathological patients,
the designed detection tool will be installed on a machine for ultrasound acquisition, giving complete freedom to radiologists of tun-
4.5. C ONCLUSION
113
Figure 4.6: Interface for real-time lesion detection proposed to radiologists to be
used during prostate biopsy.
ing detection characteristics and location through a user-friendly
graphical interface. An example of such interface is shown in figure 4.6.
A double screen or a double window on the machine display
will provide both the malignancy map visualised on the US image and the traditional B-mode signal as imaged in standard machines. Radiologists will thus receive on-line information about
suspected regions and will choose, according to his experience, the
cores of the gland to sample.
Once a prototype machine will be ready to use, a clinical trial
will be performed on new patients, analysing both the standard
double sextant sampled cores and the targeted cores obtained by
rtCAB examination. The results of the clinical trial will show the
actual efficiency of the new US-guided lesion-directed protocol for
prostate biopsy.
4.5 Conclusion
TRUS guided biopsy is a fundamental step in the current clinical
procedure for prostate cancer diagnosis. So far TRUS images have
been used to correctly position the probe during biopsy, while the
information extracted from TRUS images, which can be very use-
114
Four. I MPROVING STANDARD PROSTATE BIOPSY PROTOCOL
ful in detecting suspicious regions, has never been systematically
exploited to guide biopsy. The standard double sextant biopsy protocol equally samples the gland, randomly selecting the tissues to
be analyzed. This procedure does not take into account the risk
of sampling a possibly benignant core, not significant for cancer
detection.
In this work we propose a real-time Computer Aided Biopsy
tool, which employs Nakagami parameters extracted from TRUS
images to detect suspicious regions and guides prostate biopsy.
The proposed detection tool uses a non linear classification process which takes advantage of Nakagami parameters computed
through an original Padé Approximant estimator. The non linear
classification algorithm was trained on a 400 cores ground truth
database, collected thanks to the department of Urology at S. Orsola Hospital in Bologna. In addition, the resulting algorithm was
implemented through CUDA parallel processing and it is characterized by a frame rate of 30 fps, suitable for a clinical application.
The proposed rtCAB is able to reduce the biopsy core number
from 8-12 to 7, preserving the same diagnostic power of a classical
double sextant protocol.
C ONCLUSIONS
T
HE research activity presented in this work was developed during a PhD program at University of Bologna in Italy from 2008
to 2010 and was partly carried out at École Polytechnique Fédérale
de Lausanne in Switzerland.
This research focused on ultrasound image processing for tissue characterization with the purpose to support minimally invasive detection of pathologies. An important clinical target has been
the spine of the different studies proposed in this thesis: the realisation of a lesion-directed ultrasound-guided biopsy to improve
the detection of prostate cancer. For this reason the research activity described in this work was conducted side by side with the
researchers of the department of Urology of S.Orsola hospital in
Bologna.
As the investigated problem is the discrimination by a binary
classification scheme of benignant tissue and a tissue diagnosed
with adenocarcinoma, most of the experiments performed
throughout the research activity were based on medical ground
truth derived from biopsy findings. The classifier takes advantage
of some characteristics of ultrasound signals, examined on medical ground truth, to type imaged tissues. Several aspects of this
problem were investigated in this work: i) the ultrasound image
analysis techniques to improve characterization performance, ii)
the design of an appropriate learning scheme on medical ground
truth to assure a good quality classification, ii) the modeling of ultrasound signal to design new features useful for tissue characterization, iv) the optimization of the developed techniques for realtime computation to allow their clinical employment.
The first part of this thesis concentrated on biomedical image
analysis techniques involved in the tissue characterization procedure. Here it was proposed to include a pre-processing deconvolution step in the traditional procedure before the extraction of ul-
115
116
C ONCLUSIONS
trasonic characteristics from the signal. This additional step was
shown to improve the diagnostic skill of ultrasonic features, allowing a better discrimination between healthy and pathological tissues.
The second part of the work focused on ultrasound signal modeling. In particular it was proposed to use a sampled continuous
domain autoregressive moving average model to fit the ultrasound
data. This modeling provides the possibility to analyze the behaviour of the continuous-time signal which is than sampled in
the acquisition process to provide the actual available data. In order to derive information about the continuous-domain signal, it
was proposed a new mathematical formulation of CARMA process
based on exponential B-splines. These functions, appropriately
designed, allow to restore by interpolation the continuous time
signal from the information about its samples. In addition, a novel
maximum-likelihood CARMA identification algorithm based on exponential B-splines was proposed. This algorithm exploits an exact discretization of CARMA model and does not introduce any approximation or assumption on the sampling interval. It is, thus,
able to provide good estimates of CARMA process parameters also
in aliasing conditions and it can be useful in several image processing applications. Since in ultrasound applications the sampling frequency is quite high, there is no risk of aliasing and a simpler CARMA identification was also proposed for this case. Thanks
to the presented techniques CARMA parameters can than be extracted from the ultrasound signal. A study performed on phantom images showed that these new continuous-time descriptors
are more efficient than traditional discrete-domain ARMA parameters in characterizing the concentration of small particles in ultrasound images, called scatterers and responsible of ultrasound
pulse diffraction.
Finally a clinical study was conducted on a large medical database to assess the feasibility of the improvement of the standard
prostate cancer biopsy protocol by image processing techniques.
The developed techniques were optimized for speed up and implemented on a GPU enabled with CUDA™ architecture in order
to allow the ultrasound image processing to run in parallel with
the real-time visualization of B-mode echo. The result of this part
of the research is the development of a real-time Computer-Aided
Biopsy tool working at a frequency rate compatible with standard
ultrasound machines and providing radiologists with a visual malignancy map, displayed on standard B-mode. This tool allows to
C ONCLUSIONS
117
perform a lesion-detected biopsy with the same diagnostic power
of the standard prostate biopsy protocol, but with a reduction in
the number of cores to sample from 8-12 to 7.
Future research will be targeted to the improvement of the classification scheme on a larger medical database with the addition of
textural features and to the installation of the rtCAB tool on a standard machine in order to start a clinical trial at S.Orsola hospital.
Another aspect to take care of is a faster implementation of the
deconvolution algorithm and its assessment on a large medical database for a future inclusion of this pre-processing step in the rtCAB processing.
In addition, further investigations should be performed in order to validate the usefulness of continuous-domain ARMA parameters on in vivo images to characterize biological tissues.
Finally, a delicate aspect of this research is the presence in the
medical database of regions on which we cannot obtain full knowledge. Precancerous lesions, as well as regions with a very small
tumour percentage, cannot be treated as labelled data. For this
reason, and in order to capture as much information as possible
from the medical ground truth, we are also investigating the properties of semi-supervised classification approaches to improve tissue characterization performances.
A
F EATURES FOR ULTRASONIC
TISSUE TYPING
The multi-feature approach used in this work, based on attributes
of different nature, will be briefly described in the following paragraphs. In general, every single group of features produces several
attributes, but, as described in section 2.1, a correlation analysis
performed on each group provides a selection of attributes from a
certain set. Once an output feature image is obtained, a statistical parameter is computed on feature values in each ROI to have a
local information.
Wavelet Transform: this feature groups the wavelet packet coefficients of the RF signal, decomposed in four bands. According to
the correlation analysis only the first and third bands coefficients
are considered meaningful.
Polynomial Fit of Wavelet Packet Transform: the computation of this feature consists of the extraction of the four bands decomposition WT coefficients of the RF signal. The second step is a
fitting of the computed WT coefficients for each point of RF signal
with a third order polynomial [28]. The result is a projection of the
WT coefficients in a lower dimensionality space. Only the second
order coefficient is considered meaningful by correlation analysis.
Wavelet Decomposition: this feature consists on the decomposition of RF signal in its coherent and diffuse parts [44]. The
RF signal is modelled as the sum of a diffuse component, due to
the interaction of the pulse with resolvable scatterers, and a diffuse component due to the randomly located unresolvable scatterers. When a coherent component is present in time, the scaleaveraged wavelet power is characterized by larger peaks, used to
detect the time location of coherent scatterers through thresholding. The coherent component is then reconstructed according to
the model, through the superposition of Gaussian modulated si-
119
120
A PPENDIX A. F EATURES FOR ULTRASONIC TISSUE TYPING
nusoids. The diffuse component is derived from the difference of
the RF signal and the coherent component. An inverse wavelet
packet transform is applied to both coherent and diffuse signals
and the coefficients in the 8t h and 9t h node of the packet wavelet
tree are used as attributes. The generated feature set is then composed of four features but, according to correlation analysis, just
the mean of the third parameter, WDES Diff. Proj. n.8, is meaningful for classification purposes. The name WDES indicates the used
decomposition algorithm: wavelet-based despeckling [102].
Central Frequency: this feature is an estimate of the mean central frequency of the RF signal. According to the algorithm based
on amplitude spectrum magnitude, the integrated backscatter coefficient is first evaluated as an average of the PSD of the RF signal,
computed by means of FFT. Next, the mean central frequency is
estimated as the
of RF signal PSD, i.e. it is
P
P first order momentum
BW f ∗P SD(f )
BW f ∗P SD(f )
P
=
.
computed as
IBS
BW P SD(f )
Attenuation: this group of features is based on the main idea
that slope and intercept of the linear fitting of the RF signal central
frequency are measures of signal attenuation and thus of backscatter [20]. The first step to compute this feature consists of RF signal PSD evaluation, which can be performed through several algorithms: zero crossing applied on a sliding window, FFT, fitting
with an auto-regressive (AR) model of a certain order. In the last
method the meaningful frequencies corresponding to non trivial
maximums and minimums of the AR spectrum can be computed
analytically since they only depend on AR coefficients. A last way
to compute this attenuation features is based on modified PSD
whose linear trending is removed and a smoothing procedure is
applied to minimize the variance of spectral slope estimation [103].
After PSD estimation, the central frequency over the signal and its
linear fitting are computed. The whole feature set contains 18 attributes but the correlation analysis select 13 of them, privileging
computation based on AR model. In particular, in the ranked feature sets in Table 2.2, only the intercept of central frequency fitting
in the case of 2nd order AR model estimated through Burg method
(Intercept_AR(2) burg) and the intercept of the second frequency
extracted in the 3r d order AR model computed through a LMS procedure (Intercept2_AR(3) lmsd) are selected.
B-mode: this single attribute represents the B-mode value of
the analyzed RF signal, computed through its Hilbert Transform.
Nakagami: this group of features comes from a statistical model
A PPENDIX A. F EATURES SET FOR ULTRASONIC TISSUE TYPING
121
of the RF signal based on the Nakagami distribution. Shape and
scale parameters of Nakagami distributions estimated on the RF
signal are shown to be correlated to the density and intensity of
scatterers and can be computed through estimation of some RF
signal statistical moments [26]. The attributes used in this work
are the logarithm of the shape and scale parameters computed on
the envelope of the RF signal (Nakagami logm and Nakagami logOmega). The same parameters can be extracted from the coherent and diffuse components of the RF signal. Correlation analysis
selects shape and scale parameters when they are extracted from
the RF signal, while only the scale parameter (Naka w Diff. Proj.
n. node_number, Naka w Coher. Proj. n. node_number) is selected
when these attributes are extracted from the RF signal diffuse and
coherent component, thus this feature set contains four different
attributes.
Statistic: both the mean and the variance of the pure RF signal
in each ROI are considered by the correlation analysis as important
attributes for diagnostic purposes. Thus these feature set contains
two different attributes.
Haralick: this textural feature contains statistical attributes generated from co-occurrence matrix of the ROI in the direction of
45◦ [30]. Several different statistical parameters can be computed
from the co-occurrence matrix in each direction, but the correlation analysis selects as meaningful features only four of them: sum
of squares (Haralick Sum of Squares), correlation, entropy and average of the histogram of the sum of grey levels.
Unser: this feature contains statistical attributes generated from
histogram of the sum and difference of grey levels in the ROI [31].
Nine statistical parameters can be computed from sum and difference histogram: mean, variance, contrast, homogeneity, cluster
shade, cluster prominence, energy, correlation and entropy. These
measures can be replicated on all the four angular directions (0◦ ,
45◦ , 90◦ , 135◦ ), giving a total of 36 attributes. Correlation analysis
selects the following attributes: mean for the horizontal direction,
correlation for the direction of 45◦ , mean and contrast for the vertical direction and correlation, cluster prominence, entropy and
homogeneity for the direction of 135◦ , providing a 9 attributes feature set.
Fractal: this group of parameters is based on computation of
progressive binarized versions of the B-mode image [32], [36] obtained through ten different binary thresholding operations. Lattice of different grid size highlights the number of squares contain-
122
A PPENDIX A. F EATURES FOR ULTRASONIC TISSUE TYPING
ing active pixels, expressed through the lattice size and the fractal
dimension. Finally, fractal features are computed through linear
regression of the number of active squares with respect to the lattice size. Both slope (Alpha) and intercept (Beta) are used as features, providing 20 different fractal attributes, named Fractal(1) Alpha/Beta num_attribute. A different kind of fractal feature is based
on a sequence of imaginary extensions of original image toward
upper and lower grey levels, considered as surfaces of a volume
with a certain radius surrounding the original image [104], [36].
The surface area contains information about the difference of two
subsequent extended images and it is expressed through fractal dimension and volume radius. Choosing ten different radius values
and performing a piece-wise linear regression of the surface area,
the derived slopes (Alpha) and intercepts (Beta) are used as textural features. This fractal feature provides 20 more different attributes, named Fractal(2) Alpha/Beta num_attribute. Correlation
analysis selected 7 attributes from the first group of fractal features
and 9 from the second group, providing a set of 16 features.
B
S TATE OF THE ART IN TISSUE
CHARACTERIZATION
The technical characteristics of previously published methods for
prostate tissue characterization, discussed in chapter 2, are summarized in table B.1 in order to show a complete framework about
the state of the art in ultrasound-based computer-aided detection
techniques.
123
Work
Author
Year
Gold
Standard
# images/
# patients
Basset [16]
1993
biopsy
37/16
Huynen [17]
1994
biopsy
98/51
Houston [18]
1995
biopsy
25/25
Schmitz [19]
1999
prostatectomy
200/33
Scheipers [20]
2003
prostatectomy
-/100
Feleppa [21]
2004
Mohamed [22]
2005
Llobet [53]
Mohamed [23]
2007
2008
biopsy and
prostatectomy
physician
visual
inspection
biopsy
-/200
33/20
Malign/
Benign
6/10 patients
46/5 patients
11/14
patients
33/0 patients
40517/129967
ROIs
110/909 ROIs
20/0 patients
4593/289
physician
visual
inspection
33/20
Han [24]
2008
-
51/51
HistoScanning™ [54]
2009
prostatectomy
29/29
Dataset
# ROIs
51/0 patients
29/0 patients
ROI size
2
Train/Test
Segmentation
rectangular
around
needle
rectangular
around
needle
rectangular
around
needle
-
37
512×512
1.45 cm
-/-
-
512×512
3×3 mm2
66/32
25
512×512
121×20
-/-
3400
-
170000
-
1019
-
64×64
180/20
96
-
irregular
80/16
128×16
0.1 cm2
128×16
0.1 cm2
32/1
99/1
manual by
radiologist
-
Techniques
Features
Feature
selection
Textural
no
Textural
no
Textural
no
acoustical
and textural
attenuation,
scattering,
textural
spectral and
clinical data
Textural
Gabor filter
4593
20/0 patients
Image size
108
384 × 288
-
2500 pixels
202/87
patients
irregular
90/18
rectangular
around
needle
Gabor filter
200
350 × 300
25 × 25
5/46
-
3D
0.08 cm2
15/14
thresholded
histogram
equalization
-
SFS on
Mahalanobis
distance
ROC area
based
wrapper
-
mutual
information
Textural
no
PSD,
rotational
invariance
technique
textural,
location,
ellipticity
statistical
Particle
swarm
optimization
no
-
Classifier
SE
Results
SP
Acc
one
parameter
decision tree
multiparameter
decision tree
Decision tree
83
71
-
-
80
88.20
-
-
73
86
80
-
-
82
88
-
-
2 parallel
neuro-fuzzy
systems
artificial
neural
network
SVM
-
-
75
86
-
-
80
85
83.3
100
93.75
-
Az
hidden
Markov
model
SVM
68
53
61.6
60.1
83.3
100
94.4
-
SVM
92
95.9
-
-
95
-
-
-
integration of
3 algorithms
124 A PPENDIX B. S TATE OF THE ART IN TISSUE CHARACTERIZATION
Table B.1: Published methods for ultrasound-based prostate tissue characterization
C
I NITIAL CONDITIONS IN
CARMA ESTIMATION
The task of band-pass power spectrum allocation is repeatedly required by the algorithm of section 3.5. Each such allocation corresponds to a single local minima extraction and it is based on equating the derivative of the power spectrum Φ(s; θ) to zero while substituting s = j ωk , where ωk is the required frequency of maximum
response. The coefficients of ΦL (s) are then obtained by numerical minimizing the derivative value at ωk . The initial value for the
innovation variance is determined by the autocorrelation function
of the model, which is the inverse Laplace transform of the power
spectrum, and by the variance of the available data. A description
of the allocation workflow is given in Figure C.1.
125
126
A PPENDIX C. I NITIAL CONDITIONS IN CARMA ESTIMATION
Figure C.1: Allocating initial conditions. Every local minima of the likelihood
function (3.49) is obtained by a different set of initial conditions. This workflow is
part of the estimation algorithm of Figure 3.5.
D
N ON SYMMETRIC EXPONENTIAL
B- SPLINES
Considering a LTI operator L s,r representing a the whitening model
of a CARMA process of orders (p, q) with poles s and zeroes r, an
exponential spline ([105, 56, 57]) is defined as
s(t ) =
X
a[k]h(t − tk ; s, r)
(D.1)
k∈Z
where h(t ; s, r) is the Green function of L s,r . The function s(t ) can
be reconstructed in exponential B-spline basis, more appropriate
descriptor thanks to their compact support unlike the Green function:
X
c[k]β(t − k; s, r)
(D.2)
s(t ) =
k∈Z
The exponential B-spline is a combination of shifted and scaled
version of the Green function and can be defined by applying a
localization operator ∆s to the Green function itself:
β(t ; s, r) = ∆s {h(t ; s, r)}
(D.3)
The localization operator is defined in the Fourier domain:
∆(ω; s) =
p
Y
(1 − e − j ω+s k )
(D.4)
k=1
As the Green function is, by definition, the impulsive response of
the inverse operator L −1
s,r , its Fourier transform is the CARMA system transfer function (3.9). Thus the Fourier transform of β(t ; s, r)
is
q
p
Y
1 − e − j ω+s k Y
j ω − rk
(D.5)
β̂(ω; s, r) =
j ω − sk k=1
k=1
127
128
A PPENDIX D. N ON SYMMETRIC EXPONENTIAL B- SPLINES
In this context the Green function h(t ; s, r) can be expressed as
an exponential spline:
h(t ; s, r) =
X
p[k; s]β(t − k; s, r)
(D.6)
k∈Z
where the coefficients of the interpolation present the following ztransform:
1
(D.7)
P d− (z; s) = Qp
(1
−
e s k z −1 )
k=1
The autocorrelation function ϕ(t ; s, r) of the CARMA process
can be expressed in function of exponential B-splines:
ϕ(t ;s,r)
σ2 h(t ;s,r) ∗ h(−t ;s,r)
X
X
σ2
p[k;s]β(t − k;s,r) ∗
p[k ′ ;s]β(−t − k ′ ;s,r)
=
=
k ′ ∈Z
k∈Z
2
X
σ
=
σ2
=
p[k;s]
k∈Z
p
Y
e si
i =1
X
′
p[k ;s]β(t − k;s,r) ∗ β(−t − k ′ ;s,r)
k ′ ∈Z
X
p[k;s]
X
p[k ′ ;s]β(t − k;s,r) ∗ β(t + k ′ + p;−s,−r)
k ′ ∈Z
k∈Z
The last passage is possible thanks to the symmetry property
of exponential B-splines [56]. In addition the convolution of two
exponential B-spline is still an exponential B-spline with as coefficients the full set of parameters of the first and the second splines. In this case the resultant exponential B-splines has parameters −s, −r : s, r, where the symbol : indicates the union of coefficients. This property leads to the following development:
ϕ(t ;s,r)
=
σ2
=
2
p
Y
X
X
e si
p[k;s]
p[k − k ′′ − p;s]β(t − k ′′ ;−s,−r : s,r)
i =1
σ
X
k∈Z
′′
k ′′ ∈Z
p̃[k + p;−s : s]β(t − k ′′ ;−s,−r : s,r)
(D.8)
k ′′ ∈Z
The definition of the autocorrelation function of a CARMA process
and its spectrum depends on the symmetric exponential B-splines
β(t − k ′′ ; −s, −r : s, r), which was for this reason introduced in section 3.2. In particular, since considering all parameters (−s, −r :
s, r) is redundant, because they are opposite each other, in section 3.2 we take a single parameter vector θ representing poles
and zeroes of the system. Unlike (3.29), here Φd (z, θ) = σ2 P d (z; θ) ·
Bd (z; θ) uses some slightly different definitions of P d (z; θ)
A PPENDIX D. N ON SYMMETRIC EXPONENTIAL B- SPLINES
129
and β(t ; θ). The definition of P d (z; θ) as z-transform of interpolation coefficients p̃[k] in (D.8) is as follows
P d (z; θ) = Qp
k=1
zp
¡
Qp
k=1
e sk
¢¡
1 − e s k z 1 − e s k z −1
¢
(D.9)
©
ª
As β(t ; θ) = F −1 β̂(ω; s, r) · β̂(ω; −s, −r) , let us consider directly the
definition of β̂(ω; θ) in the frequency domain:
β̂(ω;θ)
=
q
Y
( j ω − r k )( j ω + r k ) ·
k=1
=
(−1)q
p
Y
(1 − e − j ω+sk )(1 − e − j ω−sk )
( j ω − s k )( j ω + s k )
k=1
p
q
Y
Y
e − j ω−sk ·
( j ω − r k )(− j ω − r k )
(D.10)
k=1
·
k=1
p
Y
(1 − e j ω+sk )(1 − e − j ω+sk )
k=1
( j ω − s k )(− j ω − s k )
These definitions allow to define the alternative expression for
the variance of the innovation signal of the sampled CARMA process expressed in corollary 1.
N OTATIONS
Chapter One. Background
Symbol
Meaning
x:
xi :
Nf :
x[n, m]:
w [n, m]:
H(.):
h[n, m]:
µ[n, m]:
xk :
ck :
f (x):
w:
b:
αi :
φ(.):
K(., .):
z:
sv :
ai :
σ2 :
generic ROI feature vector of length N f
i -th feature of a ROI
number of features
2D image data
2D tissue response
acquisition system convolution operator
acquisition system point spread function
additive noise corrupting tissue response
k -th ROI feature vector of length N f
k -th ROI class label
generic predictive model
linear classifier coefficients vector
linear classifier trained threshold
coefficients for support vector linear predictive model
non linear mapping for kernel-based classification
kernel function
remapped ROI feature vector in kernel space
set of reference feature vectors
i -th coefficient for kernel-based classification rule
variance of radial basis function kernel
J:
µi :
Σi :
T P:
F P:
T N:
F N:
ACC :
SE :
SP :
PPV :
N0 :
p p v0 :
Mahalanobis distance operator
intra-class mean value of samples belonging to class i
covariance matrix of samples belonging to class i
number of well classified unhealthy (positive) samples
number of bad classified healthy (negative) samples
number of well classified healthy (negative) samples
number of bad classified unhealthy (positive) samples
number of well classified samples over the total
number of well classified unhealthy samples over all the unhealthy ones
number of well classified healthy samples over all the healthy samples
n. of w. c. unhealthy samples over the samples classified as unhealthy
number of cores in standard prostate biopsy protocol
PPV of standard prostate biopsy protocol
rb f
131
132
N OTATIONS
Chapter Two. Ultrasound image analysis and
interpretation
Symbol
Meaning
S:
fi :
N:
NROI :
N pi xel :
c:
D(., .):
R(.):
I (., .):
C:
K:
d (., .):
Xi :
Xi , j :
X:
W (.):
yk :
Σ:
P(Σ):
P(Σ|X):
P(X|Σ):
Ci = k :
Σi = k :
µki, j :
ROI-based feature data matrix, feature set
i -th feature vector of length NROI
number of features
number of regions of interest
number of pixels
data class vector of length NROI
feature set relevance with respect to class
feature set redundancy
mutual information of feature set with respect to the class
K-means output segmented image, clustering of pixels
number of clusters for segmentation
generic distance between pixels
feature vector of the i -th pixel of length N f
j -th feature value for the i -th pixel
pixel-based feature data matrix of size N pi xel × N f
energy function of K-means algorithm
k -th cluster center
final estimated segmented image
a priori probability of segmented image
a posteriori probability of segmented image, given the observed data
probability of the observed data conditioned to the segmented image
assignment of class k to pixel i by K-means
final assignment of class k to pixel i
mean value of Gaussian distribution modelling j -th feature value
U (.):
V (.):
δi :
x[n]:
w [n]:
h[n]:
ǫk :
x̂[n]:
e[n]:
E (., .):
ǫ:
λ:
F:
φ(.):
x:
z:
w:
b:
D F (.):
SB :
SW :
K(., .):
functional to minimize by segmentation algorithm
Markov Random Field potential function
set of pixels in the neighbourhood of pixel i
RF discrete signal
true tissue response
acquisition system point spread function
predictive deconvolution coefficients
estimated prediction of x[n]
discrete innovation signal, estimated tissue response
predictive deconvolution error function
predictive deconvolution coefficients vector
predictive deconvolution forgetting factor
feature space
non linear mapping
generic feature vector of one ROI
mapped ROI feature vector
linear classifier coefficients vector
linear classifier trained threshold
Fisher criterion
between-class scatter matrix for mapped data
within-class scatter matrix for mapped data
kernel function
σki, j :
standard deviation of Gaussian distribution modelling j -th feature value
Appendix A. Features for ultrasonic tissue typing
f:
frequency
PSD( f ):
power spectral density of RF signal
N OTATIONS
133
Chapter Three. A CARMA model for ultrasound signals
Symbol
Meaning
f (x; y 1 , y 2 ):
continuous-time function f is expressed in x
and depends on parameters y 1 and y 2
discrete-time sequence f is expressed in n
and depends on parameters y 1 and y 2
radial frequency [rad/time-unit] in Fourier transform
or normalized frequency [rad/sample] in DTFT
RF signal, discrete ARMA process
discrete ARMA process samples vector
AR coefficients of x[n]
MA coefficients of x[n]
continuous-time ARMA process, continuous-time RF signal
continuous-time ARMA model innovation signal, CD tissue response
continuous-time derivative operator
LTI operator (AR part)
LTI operator (MA part)
variance of continuous-time innovation signal
order of A(D)
order of B(D)
CARMA process parameters vector
zeros of CARMA model
poles of CARMA model
vector of CARMA model zeroes
vector of CARMA model poles
CARMA process real parameters vector
coefficients of A(D)
coefficients of B(D)
CARMA process autocorrelation function
sampled CARMA process autocorrelation sequence
power spectrum of x(t ); Fourier transform of autocorrelation function
Laplace transform of ϕ(t )
power spectrum of x[n]; DTFT of sampled autocorrelation sequence
z -transform of ϕ[n]
CARMA system impulse response
CARMA system transfer function
coefficients of partial fraction decomposition of Φ(s)
exponential B-spline
samples of exponential B-spline
Laplace transform of localization filter
φ(t ) interpolation coefficients in exponential B-spline basis
z -transform of the sequence p[n]
z -transform of the sequence β[n]
DTFT of the sequence β[n]
Fourier transform of β(t )
ARMA model given by sampling of CARMA model
discrete-time innovation signal variance in sampled CARMA model
poles of ARMA model
zeroes of ARMA model
causal part of Bd (z)
anti-causal part of Bd (z)
causal part of P d (z)
anti-causal part of P d (z)
fundamental interpolation kernel of sequence ϕ[n] to rebuild ϕ(t )
element of Fisher information matrix
sampling interval
number of samples of discrete process
probability density function of x
autocorrelation matrix of x
f [n; y 1 , y 2 ]:
ω:
x[n]:
x:
ck :
dk :
x(t ):
w (t ):
D:
A(D):
B(D):
σ2 :
p:
q:
θ:
rk :
sk :
r:
s:
θ0 :
ak :
bk :
ϕ(t ):
ϕ[n]:
Φ( j ω):
Φ(s):
Φd (e j ω ):
Φd (z):
h(t ):
ĥ(ω):
αk :
β(t ):
β[n]:
∆(s):
p[n]:
P d (z):
Bd (z):
β̂d (ω):
β̂(ω):
Hd (z):
σ2d :
ρk :
νk :
B − (z):
d
B + (z):
d
−
P (z):
d
+
P (z):
d
η(t ):
I k,l (θ):
T:
N:
pd f (x):
Σ:
134
|Σ|:
l (θ; x):
˜ x):
l(θ;
G d (z):
g θ [n]:
κ(θ):
c[n]:
L[m, n]:
TM :
t M [n]:
l˜W (θ; x):
x̂d (ω):
ǫ(a k ):
â k,n :
ωpeak :
ωmax :
Tmax :
a:
b:
c:
d:
dk (s, r):
d(s, r):
J [x,y] :
Ncl asses :
N OTATIONS
determinant of Σ
likelihood function of the model to data x
approximation of l (θ; x)
inverse of digital filter Hd (z)
coefficients of digital filter G d (z)
corrective term of log-likelihood approximation
Fourier coefficients of logarithm of ARMA spectrum
lower triangular matrix describing filter g θ
infinite Toeplitz matrix describing ARMA spectrum
sequence associated to the matrix TM
Whittle approximation of log-likelihood function
DTFT of x
estimation error over parameter a k
estimation of parameter a k in the n -th experiment
frequency at which Φ(ω) is maximum
frequency at which Φ(ω) is at −10 dB respect to Φ(ωpeak )
sampling time corresponding to frequency ωmax
vector of coefficients a k of continuous-domain AR part
vector of coefficients b k of continuous-domain MA part
vector of coefficients c k of discrete-domain AR part
vector of coefficients dk of discrete-domain MA part
coefficients of the casual polynomial of Bd (z)
vector of coefficients dk (s, r)
Mahalanobis distance between samples corresponding to
scatterer concentration x and y
number of scatterers concentrations in the classification problem
Appendix B. Initial conditions in CARMA estimation
k:
band index for multi band likelihood optimization
ω0 :
peak frequency for the base band solution of likelihood maximization
ωk :
estimated peak frequency in the k -th band
A(s):
casual part of denominator of Φ(s)
B(s):
casual part of numerator of Φ(s)
F (s):
first derivative of Φ(s)
σ̂2 :
variance of data x[n]
Appendix C. Non symmetric exponential B-splines
L s,r :
LTI operator of whitening model of a CARMA process,
inverse CARMA system
s(t ):
exponential spline
h(t ):
Green function of L s,r , impulse response of CARMA model (L −1
s,r )
a[k]:
coefficients of a generic spline interpolation by Green function of L s,r
c[k]:
coefficients of a generic spline interpolation by exponential B-splines
∆s :
localization operator of green function in the time/space domain
p[k]:
coefficients for green function interpolation by exponential B-splines
p̃[k]:
coefficients for autocorrelation function interpolation
by exponential B-splines
N OTATIONS
135
Chapter Four. Improving current prostate biopsy protocol
Symbol
Meaning
m:
Ω:
ψ(m):
f (m):
fˆ(m):
pk :
qk :
C:
Xtr n :
Xt st :
ctr n :
ct st :
X̂tr n :
ĉtr n :
σr b f :
σ̂r b f :
d:
t h:
t˜h :
tˆh :
ACC :
T P:
F P:
PPV :
p p v0 :
p p v1 :
N0 :
N1 :
Nakagami shape parameter
Nakagami scale parameter
Digamma function
non linear function to determine m
approximation of f (m)
fˆ(m) numerator coefficients
fˆ(m) denominator coefficients
generic classification model
generic ROI feature matrix of size NROI in training set ×N f
generic ROI feature matrix of size NROI in testing set ×N f
generic training set class label vector
generic testing set class label vector
selected Xtr n
labels corresponding to the selected Xtr n
variable variance of RBF kernel for GDA processing
selected RBF kernel variance
linear classification discriminant function
variable threshold to apply to discriminant function
set of thresholds t h corresponding to maximum PPV
maximum threshold in the set t˜h
ROI-based accuracy
core-based true positives
core-based false positives
core-based positive predictive value
standard protocol PPV
proposed rtCAB PPV
standard protocol number of cores to sample
roposed rtCAB number of cores to sample
136
N OTATIONS
Abbreviations
Acronym
Meaning
ACC:
AR:
ARMA:
ASAP:
BPH:
BW:
CAR:
CAD:
CARMA:
CD:
CRB:
CT:
CZ:
DFT:
DRE:
DTFT:
FFT:
FLD:
FN:
FP:
FS:
GDA:
GPU:
GUI:
IQ:
LDA:
LTI:
MA:
MAP:
MIHFS:
ML:
MLE:
mRMR:
MRF:
MRI:
nRMSE:
PCA:
PIN:
PPV:
PSA:
PSD:
PSF:
PZ:
RBF:
RF:
RLS:
ROI:
rtCAB:
SE:
SFS:
SP:
SVM:
TGC:
TN:
TP:
TRUS:
TZ:
US:
accuracy
auto-regressive
auto-regressive moving-average
atypical small acinar proliferation
benign hyperplasia
band width
continuous-domain auto-regressive
computer-aided diagnosis/detection
continuous-domain auto-regressive moving-average
continuous-domain
Cramér-Rao bound
computerized tomography
central zone
discrete Fourier transform
digital rectal examination
discrete-time Fourier transform
fast Fourier transform
fisher linear discriminant
false negatives
false positives
feature selection
generalized discriminant analysis
graphic processing unit
graphical user interface
in-phase quadrature signal
linear discriminant analysis
linear time invariant
moving-average
maximum a posteriori
mutual information hybrid feature selection
maximum likelihood
maximum likelihood estimator
min-redundant max-relevance
Markov random field
Magnetic risonance imaging
normalized root mean squared error
principal component analysis
prostatic intraepithelial neoplasia
positive predictive value
prostate-specific antigen
power spectral density
point spread function
peripheral zone
radial basis function
radio-frequency
recursive least squares
region of interest
real-time computer-aided biopsy
sensitivity
sequential forward selection
specificity
support vector machine
time-gain compensation
true negatives
true positives
trans-rectal ultrasound
transition zone
ultrasound
P UBLICATIONS
• 2011 (under revision)
IEEE Transactions on Signal Processing
H. Kirshner, S. Maggio, M. Unser.
A Sampling Theory Approach for Continuous ARMA Identification
• 2011 (under revision)
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
M. Alessandrini, S. Maggio, J. Porée, L. De Marchi, N. Speciale, E. Franceschini, O. Bernard, O. Basset.
An expectation maximization framework for improved tissue response characterization
• 2011 May
Proceedings SampTA2011
H. Kirshner, S. Maggio, M. Unser.
Maximum-Likelihood Identification of Sampled Gaussian Processes
• 2011 March
Proceedings ISBI2011
S. Maggio , M. Alessandrini, N. Speciale, O. Bernard, D. Vray, O. Basset, M.
Unser.
Continuous-domain ARMA modeling for ultrasound tissue characterization
• 2011 February
SPIE Medical Imaging
M. Alessandrini, S. Maggio, J. Porée, L. De Marchi, N. Speciale, E. Franceschini, O. Bernard, O. Basset.
An expectation maximization framework for an improved tissue characterization using ultrasounds
• 2010 November
17 Congresso Nazionale SIEUN
N. Testoni, N. Speciale, A. Bertaccini, D. Marchiori, M. Fiorentino, F. Manferrari, R. Schiavina, R. Cividini, F. Galluzzo, S. Maggio, E. Biagi, L. Masotti,
G. Masetti, G. Martorana.
A retrospective study to reduce prostate biopsy cores by a real time interactive
tool
• 2010 October
Proceedings IEEE IUS2010
N. Testoni, S. Maggio, F. Galluzzo, L. De Marchi, N. Speciale.
rtCAB: a tool for reducing unnecessary prostate biopsy cores
137
138
P UBLICATIONS
• 2010 August
Proceedings EUSIPCO 2010
S. Maggio, H. Kirshner, M. Unser
Continuous-time AR model identification: does sampling rate really matter?
• 2010 February
IEEE Transactions on Medical Imaging
S. Maggio, A. Palladini, L. De Marchi, M. Alessandrini, N. Speciale, G. Masetti
Predictive deconvolution and hybrid feature Selection for Computer-Aided
Detection of prostate cancer
• 2009 March
Proceedings International Symposium on Acoustical Imaging
M. Scebran, A. Palladini, S. Maggio, L. De Marchi, N. Speciale
Automatic regions of interests segmentation for computer aided classification
of prostate TRUS images
• 2008 November
Proceedings IEEE IUS2008
S. Maggio, L. De Marchi, M. Alessandrini, N. Speciale
Computer aided detection of prostate cancer based on GDA and predictive
deconvolution
• 2005 November
WSEAS Transactions on Systems
S. Maggio, N. Testoni, L. De Marchi, N. Speciale, G. Masetti
Ultrasound Images Enhancement by means of Deconvolution Algorithms in
the Wavelet Domain
• 2005 September
WSEAS ISCGAV2005
S. Maggio, N. Testoni, L. De Marchi, N. Speciale, G. Masetti
Wavelet-based Deconvolution Algorithms Applied to Ultrasound Images
B IBLIOGRAPHY
[1]
K. Doi, “Computer-aided diagnosis in medical imaging: Historical review, current status
and future potential,” Computerized Medical Imaging and Graphics, vol. 31, pp. 198–211,
2007.
[2]
T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out. Elsevier Inc., 2004.
[3]
J. Noble and D. Boukerroui, “Ultrasound image segmentation: A survey,” IEEE Transaction on Medical Imaging, vol. 25(8), p. 987:1010, 2006.
[4]
D. Pham, C. Xu, and J. Prince, “Current methods in medical image segmentation,” Annual Reviews in Biomedical Engineering, vol. 2(1), pp. 315–337, 2000.
[5]
P. Boccacci and M. Bertero, Introduction to Inverse Problems in Imaging. Taylor & Francis,
1998.
[6]
O. Michailovich and D. Adam, Blind Image Deconvolution: theory and applications, ch. 5:
Deconvolution of Medical Images from Microscopic to whole body images, pp. 169–230.
CRC Press, 2007.
[7]
H. Liu and L. Yu, “Toward integrating feature selection algorithms for classification
and clustering,” IEEE Transactions on Knowledge and Data Engineering, vol. 17, no. 4,
pp. 491–502, 2005.
[8]
A. R. Webb, Statistical Pattern Recognition. Wiley, 2002.
[9]
M. Aizerman, E. Braverman, and L. Rozonoer, “Theoretical foundations of the potential function method in pattern recognition learning,” Automation and Remote Control,
vol. 25, pp. 821–837, 1964.
[10]
C. Abate-Shen and M. Shen, “Molecular genetics of prostate cancer,” Genes & Development, vol. 14(19), p. 2410âĂŞ2434, 2000.
[11]
J. McAninch and E. Tanagho, Smith’s General Urology. McGraw-Hill/Appleton & Lange,
17th ed., 2008.
[12]
P. Humphrey, “Gleason grading and prognostic factors in carcinoma of the prostate,”
Modern Pathology, vol. 17(3), pp. 292–306, 2004.
[13]
J. Raja, N. Ramachandran, G. Munneke, and U. Patel, “Current status of transrectal
ultrasound-guided prostate biopsy in the diagnosis of prostate cancer,” Clinical Radiology, vol. 61, pp. 142–153, 2006.
[14]
J. De La Rosette, R. Giesen, A. Huynen, R. Aarnink, MP Van Iersel, F. Debruyne, and H. Wijkstra, “Automated analysis and interpretation of transrectal ultrasonography images in
patients with prostatitis,” European urology, vol. 27(1), pp. 47–53, 1995.
[15]
M. Essink-Bot, H. de Koning, H. Nijs, W. Kirkels, P. V. der Maas, and F. Schroder, “Shortterm effects of population-based screening for prostate cancer on health-related quality
of life,” Journal of The National Cancer Institute, vol. 90, pp. 925–931, 1998.
139
140
B IBLIOGRAPHY
[16]
O. Basset, Z. Sun, J. Mestas, and G. Gimenez, “Texture analyisis of ultrasonic images of
the prostate by means of co-occurrence matrix,” Ultrasonic Imaging, vol. 15, pp. 218–
237, 1993.
[17]
A. Huynen, R. Giesen, J. de la Rosette, R. Aarnink, F. Debruyne, and H. Wijkstra, “Analysis
of ultrasonographic prostate images for the detection of prostatic carcinoma: the automated urologic diagnostic expert system,” Ultrasound in Medicine and Biology, vol. 20,
no. 1, pp. 1–10, 1994.
[18]
A. Houston, S. Premkumar, D. Pitts, and R. Babaian, “Prostate ultrasound image analysis: Localization of cancer lesions to assist biopsy,” Proceedings of the Eighth IEEE Symposium on Computer-Based Medical Systems, pp. 94–101, 1995.
[19]
G. Schmitz, H. Ermert, and T. Senge, “Tissue-characterization of the prostate using radio frequency ultrasonic signals,” IEEE Transactions on Ultrasonics, Ferroelectrics and
Frequency Control, vol. 46, pp. 126–138, 1999.
[20]
U. Scheipers, H. Ermert, H. Garcia-Schurmann, T. Senge, and S. Philippou, “Ultrasonic
multifeature tissue characterization for prostate diagnostics,” Ultrasound in Medicine
and Biology, vol. 29, no. 8, pp. 1137–1149, 2003.
[21]
E. Feleppa, J. Ketterling, C. Porter, and J. Gillespie, “Ultrasonic tissue-type imaging (tti)
for planning treatment of prostate cancer,” Proceedings of SPIE, vol. 5373, no. 223, 2004.
[22]
S. Mohamed and M. Salama, “Computer-aided diagnosis for prostate cancer using support vector machine,” Proceedings of SPIE, vol. 5744, no. 898, 2005.
[23]
S. Mohamed and M. Salama, “Prostate cancer spectral multifeature analysis using trus
images,” IEEE Transactions on Medical Imaging, vol. 27, no. 4, pp. 548–556, 2008.
[24]
S. Han, H. Lee, and J. Choi, “Computer-aided prostate cancer detection using texture
features and clinical features in ultrasound image,” Journal of Digital Imaging, vol. 21,
no. Suppl 1, pp. S121–33, 2008.
[25]
M. Moradi, P. Mousavi, and P. Abolmaesumi, “Computer-aided diagnosis of prostate cancer with emphasis on ultrasound-based approaches: a review,” Ultrasound in Medicine
and Biology, vol. 33, no. 7, pp. 1010–1028, 2007.
[26]
P. Shankar, “Ultrasonic tissue characterization using a generalized Nakagami model,”
IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 48, no. 6,
pp. 1716–1720, 2001.
[27]
E. Feleppa, W. Fair, T. Liu, A. Kalisz, W. Gnadt, and F. Lizzi, “Two-dimensional and threedimensional tissue-type imaging of the prostate based on ultrasonic spectrum analysis
and neural network classification,” Proceedings of SPIE, vol. 3982, no. 152, 2000.
[28]
L. Masotti, E. Biagi, S. Granchi, L. Breschi, E. Magrini, and F. D. Lorenzo, “Clinical test of
rules (rules: radiofrequency ultrasonic local estimators),” IEEE Ultrasonics Symposium,
vol. 3, no. 23-27 Aug., pp. 2173–2176, 2004.
[29]
G. Georgiou, F. Cohen, C. Piccoli, F. Forsberg, and B. Goldberg, “Tissue characterization using the continuous wavelet transform. part ii: Application on breast rf data,” IEEE
Transactions on Ultrasonics Ferroelectrics and Frequency Control, vol. 48, no. 2, pp. 363–
373, 2001.
[30]
R. Haralick, K. Shanmugam, and I. Dinstein, “Textural features for image classification,”
IEEE Transactions on Systems, Man and Cybernetics, vol. 3, no. 6, pp. 610–621, 1973.
[31]
M. Unser, “Sum and difference histograms for texture classification,” IEEE Transactions
on Pattern Analysis and Machine Intelligence, vol. 8, no. 1, pp. 118–125, 1986.
B IBLIOGRAPHY
141
[32]
W. Lee, Y. Chen, and K. Hsieh, “Ultrasonic liver tissues classification by fractal feature
vector based on m-band wavelet transform,” IEEE Transactions on Medical Imaging,
vol. 22, no. 3, pp. 382–392, 2003.
[33]
I. Witten and E. Frank, Tools Data Mining: Practical Machine Learning Tools and Techniques. Elsevier, 2005.
[34]
H. Peng, F. Long, and C. Ding, “Feature selection based on mutual information: criteria
of max-dependency, max-relevance and min-redundancy,” IEEE Transactions on Pattern
Analysis and Machine Intelligence, vol. 27, no. 8, pp. 1226–1238, 2005.
[35]
S. Theodoridis and K. Koutroumbas, Pattern Recognition. Academic Press, 2008.
[36]
T. Wagner, Handbook of Computer Vision and Applications, Volume 2: Signal Processing
and Pattern Recognition, ch. 12, Texture Analysis, pp. 275–308. Academic Press, 1999.
[37]
T. Pappas, “An adaptive clustering algorithm for image segmentation,” IEEE Transactions
on Acoustics, Speech, and Signal Processing, vol. 40(4), pp. 901–914, 1992.
[38]
J. Ng, R. Prager, N. Kingsbury, G. Treece, and A. Gee, “Wavelet restoration of medical
pulse-echo ultrasound images in an em framework,” IEEE Transactions on Ultrasonics,
Ferroelectrics and Frequency Control, vol. 54, no. 3, pp. 550–568, 2007.
[39]
O. Michailovich and D. Adam, “Phase unwrapping for 2-d blind deconvolution of ultrasound images,” IEEE Transactions on Medical Imaging, vol. 23, no. 1, pp. 1–7, 2004.
[40]
H. Shin, R. Prager, J. Ng, H. Gomersall, N. Kingsbury, G. Treece, and A. Gee, “Sensitivity to point-spread function parameters in medical ultrasound image deconvolution,”
Ultrasonics, vol. 49, no. 3, pp. 344–357, 2009.
[41]
O. V. Michailovich and D. Adam, “A novel approach to the 2-d blind deconvolution problem in medical ultrasound,” IEEE Transactions on Medical Imaging, vol. 24, no. 1, pp. 86–
104, 2005.
[42]
K. Peacock and S. Treitel, “Predictive deconvolution: Theory and practice,” Geophysics,
vol. 34, no. 2, pp. 155–169, 1969.
[43]
S. Haykin, Adaptive Filter Theory. Prentice Hall, 4th ed., 2002.
[44]
G. Georgiou and F. Cohen, “Tissue characterization using the continuous wavelet transform, part i: Decomposition method,” IEEE Transactions on Ultrasonics, Ferroelectrics
and Frequency Control, vol. 48, no. 2, pp. 355–363, 2001.
[45]
L. De Marchi, A. Palladini, N. Testoni, and N. Speciale, “Blurred ultrasonic images as isiaffected signals: Joint tissue response estimation and channel tracking in the proposed
paradigm,” IEEE Ultrasonics Symposium, pp. 1270–1273, 2007.
[46]
R. Neelamani, H. Choi, and R. Baraniuk, “Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems,” IEEE Transactions on Signal Processing, vol. 52,
no. 2, 2004.
[47]
S. Caporale, A. Palladini, L. D. Marchi, N. Speciale, and G. Masetti, “Wavelet-based algorithms for speckle removal from b-mode images,” Proceedings IASTED International
Conference on Biomedical Engineering, vol. 417, 2004.
[48]
[49]
T. Joachims. SVMLight: Support Vector Machine. SVM-Light Support Vector Machine
http: // svmlight. joahims. org/ , University of Dortmund, 1999.
G. Baudat and F. Anouar, “Generalized discriminant analysis using a kernel approach,”
Neural Computation, vol. 12, no. 10, pp. 2385–2404, 2000.
142
B IBLIOGRAPHY
[50]
M. Scabia, E. Biagi, and L. Masotti, “Hardware and software platform for real-time processing and visualization of echographic radiofrequency signals,” IEEE Transactions on
Ultrasonics, Ferroelectrics and Frequency Control, vol. 49, no. 10, pp. 1444–1452, 2002.
[51]
F. Lizzi, E. Feleppa, S. K. Alam, and C. Deng, “Ultrasonic spectrum analysis for tissue
evaluation,” Pattern Recognition Letters, vol. 24(4-5), pp. 637–658, 2003.
[52]
L. Masotti, E. Biagi, S. Granchi, L. Breschi, E. Magrini, and F. D. Lorenzo, “Tissue differentiation based on radiofrequency echographic signal local spectral content,” IEEE
Symposium on Ultrasonics, vol. 1, pp. 1030–1033, 2003.
[53]
R. Llobet, J. Perez-Cortes, A. Toselli, and A. Juan, “Computer-aided detection of prostate
cancer,” International Journal of Medical Informatics, vol. 76, pp. 547–556, 2007.
[54]
J. Braeckman, P. Autier, C. Garbar, M. P. Marichal, C. Soviany, R. Nir, D. Nir, D. Michielsen,
H. Bleiberg, L. Egevad, and M. Emberton, “Computer-aided ultrasonography (histoscanning): a novel technology for locating and characterizing prostate cancer,” BJU International, vol. 101, pp. 293–298, 2007.
[55]
L. Ljung, System Identification - Theory For the User, 2nd ed. PTR Prentice Hall, Upper
Saddle River, N.J., 1999.
[56]
M. Unser and T. Blu, “Cardinal exponential splines: Part i - theory and filtering algorithms,” IEEE Transactions on Signal Processing, vol. 5373, no. 4, 2005.
[57]
M. Unser and T. Blu, “Cardinal exponential splines: Part ii - think analog, act digital,”
IEEE Transactions on Signal Processing, vol. 53, no. 4, 2005.
[58]
I. Lim and B. Lee, “Lossless pole-zero modeling for speech signals,” IEEE Transactions on
Speech and Audio Processing, vol. 1(3), pp. 269–276, 1993.
[59]
H. Kirshner and M. Porat, “On the role of exponential spline in image interpolation,”
IEEE Trans. Signal Processing, vol. 18(10), pp. 2198–2208, 2009.
[60]
J. M. Francos, A. Z. Meiri, and B. Porat, “A unified texture model based on a 2-d wold like
decomposition,” IEEE Trans. Signal Processing, vol. 41, p. 2665âĂŞ2678, 1993.
[61]
S. Nemirovski and M. Porat, “On texture and image interpolation using markov models,”
Signal Processing: Image Communication, vol. 24, pp. 139–157, 2009.
[62]
K. J. Aström, Introduction to Stochastic Control Theory. Academic Press, 1970.
[63]
B. Wahlberg, “Limit results for sampled systems,” International Journal of Control,
vol. 78(3), pp. 1267–1283, 1988.
[64]
B. Wahlberg, L. Ljung, and T. Söderström, “Sampling of continuous time stochastic processes,” Control - Theory and Advanced Technology, vol. 9(1), pp. 99–112, 1993.
[65]
T. Sd̈erstrm̈, Discrete-Time Stochastic Systems, 2nd ed. London: Springer-Verlag, 2002.
[66]
M. M. E. K. Larsson and T. Söderström, “An overview of important practical aspects of
continuous-time arma system identification,” Circuits Systems Signal Processing, vol. 25,
pp. 17–46, 2006.
[67]
H. Garnier and L. Wang, eds., Identification of Continuous-time Models from Sampled
Data. London: Springer-Verlag, 2008.
[68]
R. Johansson, “Identification of continuous-time models,” IEEE Trans. Signal Processing,
vol. 42(4), pp. 887–897, 1994.
[69]
D. Pham, “Estimation of continuous-time autoregressive model from finely sampled
data,” IEEE trans. Signal Processing, vol. 48(9), pp. 2576–2584, 2000.
B IBLIOGRAPHY
143
[70]
H. Garnier, M. Mensler, and A. Richard, “Continuous-time model identification from
sampled data: Implementation issues and performance evaluation,” Int. J. Control,
vol. 76(13), pp. 1337–1357, 2003.
[71]
H. Garnier and P. Young, “Time-domain approaches to continuous-time model identification of dynamical systems from sampled data,” in Proc. Amer. Control Conf., Boston,
MA, 2004.
[72]
G. C. Goodwin, J. Yuz, and H. Garnier, “Robustness issues in continuous-time system
identification from sampled data,” in Proc. 16th IFAC World Congress, Prague, Czech Republic, 2005.
[73]
H. Fan, T. Söderström, M. Mossberg, B. Carlsson, and Y. Zou, “Estimation of continuoustime ar process parameters from discrete-time data,” IEEE Trans. Signal Processing,
vol. 47(5), pp. 1232–1244, 1999.
[74]
M. Mossberg, “Estimation of continuous-time stochastic signals from sample covariances,” IEEE trans. Signal Processing, vol. 56(2), pp. 821–825, 2008.
[75]
R. Pintelon, J. Schoukens, and Y. Rolain, Frequency-domain approach to continuous-time
system identification: Some practical aspects, ch. 8, pp. 214–248. London: SpringerVerlag, 2008.
[76]
J. Gillberg and L. Ljung, “Frequency-domain identification of continuous-time arma
models from sampled data,” Automatica, vol. 45, pp. 1371–1378, 2009.
[77]
T. Söderström, “Computing stochastic continuous-time models from arma models,” International Journal of Control, vol. 53(6), pp. 1311–1326, 1991.
[78]
H. Kirshner, M. Porat, and M. Unser, “A stochastic minimum-norm approach to image
and texture interpolation,” in European Signal Processing Conference, Denmark, 2010.
[79]
E. K. Larsson, “Limiting sampling results for continuous-time arma systems,” International Journal of Control, vol. 78(7), pp. 461–473, 2005.
[80]
B. Friedlander, “On the computation of the cramér-rao bound for arma parameter estimation,” IEEE Trans. Acoustics, Speech, and Signal Processing, vol. ASSP-32(4), pp. 721–
727, 1984.
[81]
E. K. Larrson and E. G. Larrson, “The crb for parameter estimation in irregularly sampled
continuous-time arma system,” IEEE Signal Processing Letters, vol. 11(2), pp. 197–200,
2004.
[82]
G. Szegö, “On certain hermitian forms associated with the Fourier series of a positive
function,” Communications du Séminaire Mathématique de l’Université de Lund, vol. 9,
pp. 228–238, 1952.
[83]
R. E. Hartwig and M. E. Fisher, “Asymptotic behavior of toeplitz matrices and determinants,” Archive for Rational Mechanics and Analysis, vol. 32(3), pp. 190–225, 1969.
[84]
E. Basor, “Asymptotic formulas for toeplitz determinants,” Transactions of the American
Mathematical Society, vol. 239, p. 33âĂŞ65, 1978.
[85]
R. M. Gray, “Toeplitz and circulant matrices: A review,” Foundations and Trends in Communications and Information Theory, vol. 2(3), pp. 155–239, 2006.
[86]
S. Maggio, H. Kirshner, and M. Unser, “Continuous-time ar model identification: Does
sampling rate really matter?,” in European Signal Processing Conference, Denmark, 2010.
[87]
K. J. Aström and T. Söderström, “Uniqueness of the maximum likelihood estimates of the
parameters of an arma model,” IEEE Trans. on Automatic Control, vol. AC-19(6), pp. 769–
773, 1974.
144
B IBLIOGRAPHY
[88]
E. K. Larsson, M. Mossberg, and T. Söderström, “An overview of important practical aspects of continuous-time arma system identification,” Circuits Systems Signal Processing, vol. 25(1), pp. 17–46, 2006.
[89]
H. Tsai and K. S. Chan, “Maximum likelihood estimation of linear continuous time long
memory processes with discrete time data,” Journal of the Royal Statistical Society, Series
B (Statistical Methodology), vol. 67(5), pp. 703–716, 2005.
[90]
“Maximum likelihood estimation of a class of continuous-time long-memory processes,” tech. rep., Department of Statistics and Actuarial Science, University of Iowa,
Iowa City. Tech. Rep. 324, available at http://www.stat.uiowa.edu/techrep/, 2005.
[91]
R. H. Jones and A. V. Vecchia, “Fitting continuous arma models to unequally spaced spatial data,” Journal of the American Statistical Association, vol. 88(423), pp. 947–954, 1993.
[92]
P. Stoica and R. Moses, Introduction to Spectral Analysis. New Jesey: Prentice Hall, 1997.
[93]
S. Maggio, L. D. Marchi, M. Alessandrini, and N. Speciale, “Computer aided detection
of prostate cancer based on gda and predictive deconvolution,” Proceedings IEEE IUS,
2008.
[94]
A. Abdulsadda, N. Bouaynaya, and K. Iqbal, “Stability analysis and breast tumor classification from 2d arma models of ultrasound images,” in Proceedings of the Thirty-First
Annual International Conference of the IEEE Engineering in Medicine and Biology Society
(EMBSâĂŹ09), Minneapolis MN, USA, 2009.
[95]
J. Zielinski, N. Bouaynaya, and D. Schonfeld, “Two-dimensional arma modeling for
breast cancer detection and classification,” in Proceedings of the 2010 IEEE International
Conference on Signal Processing and Communications (SPCOM2010), Bangalore, India,
July 18-21, 2010.
[96]
R. Gonzalez and R. Woods, Digital Image Processing. Prentice Hall, 2007.
[97]
C. Fritsch, A. Ibanez, and M. Parrilla, “A digital envelope detection filter for real-time
operation,” vol. 48, pp. 1287–1293, Dec. 1999.
[98]
J. H. Chang, J. T. Yen, and K. K. Shung, “High-speed digital scan converter for highfrequency ultrasound sector scanners,” Ultrasonics, vol. 48, pp. 444–452, September
2008.
[99]
M. Garland, S. Le Grand, J. Nickolls, J. Anderson, J. Hardwick, S. Morton, E. Phillips,
Y. Zhang, and V. Volkov, “Parallel computing experiences with cuda,” vol. 28, pp. 13 –27,
jul. 2008.
[100] J. Cheng and N. Beaulieu, “Maximum-likelihood based estimation of the nakagami m
parameter,” vol. 5, pp. 101 – 103, mar. 2001.
[101] G. A. Baker and P. Graves-Morris, “Padé approximants second edition,” in Encyclopedia
of Mathematics and its Applications, no. No. 59, p. 746, Cambridge University Press, 1996.
[102] F. Argenti and L. Alparone, “Speckle removal from SAR images in the undecimated
wavelet domain,” IEEE Transactions on Geoscience and Remote Sensing, vol. 40, no. 11,
pp. 2363 – 2374, 2002.
[103] D. Liu and M. Saito, “A new method for estimating the acoustic attenuation coefficient
of tissue from reflected ultrasonic signals,” IEEE Transactions on Medical Imaging, vol. 8,
no. 1, pp. 107–110, 1989.
[104] S. Peleg, J. Naor, R. Hartley, and D. Avnir, “Multiple resolution texture analysis and classification,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 68, no. 4,
pp. 518–523, 1984.
[105] I. Khalidov and M. Unser, “From differential equations to the constructio of new waveletlike bases,” IEEE Transactions on Signal processing, vol. 54, no. 4, 2006.