Medical diagnostic systems 

(Orvosbiolσgiai kιpalkotσ rendszerek) 

Ultrasound tomography 

(Ultrahang tomogrαfia) 
Miklσs Gyφngy 



Motivation 
• 	
Interpretation of B-mode images subjective and operator-dependent 

• 	
There is a need for quantitative ultrasound images 

• 	e.g. “we have a speed of sound of 1600 m/s here. This is outside the normal range for this type of tissue” rather than “hmm.., this area has a bit less speckle than usual, let me turn up the gain to see better.. aha, that’s alright” 

• 	
Elastography (presented in previous lecture) is a solution, with values of shear moduli provided by many of its implementations 

• 	
However, what about using tomography techniques in say, X-ray, that use information gained from several angles to generate (in theory, at least) quantitative images? 

• 	
After all, B-mode angular compounding already used to reduce speckle 

• 	
Information from several angles also increases resolution 



(tcos.-ssin.,tsin.+scos.)


Conventional tomography 
[Kak and Slaney 1987] 
t 
.

• 	Intensity plots or images p from several angles are projections of object function f 
t

• 
Each pixel of p arises from line integral over fp(t,.) 

• 
In case of 2-D to 1-D projections 

• 
line described by L=(tcos.-ssin.,tsin.+scos.) 

• 
p(t,.)= .Lf(x,y)ds 




• 
Projection is called the Radon transform 

• 
Task is to reconstruct f from p 


t 



. Radon transform p(t,.) of f(x,y) 



Overview of lecture 
• 
Tomographic projections in X-ray CT 

• 
Reconstruction of object function f 

• 
Ultrasound tomography: challenges, opportunities 

• 
Modifications of reconstruction algorithm for ultrasound 




X-ray computed tomography (CT) [Kak and Slaney 1987, Ch 4.1] 
• 
X-rays transmitted at one end travel in straight lines through attenuating tissue 

• 
Remaining photons received at other end 


• 
For every unit distance, proportion ΅(s) of photons is lost 

• 
Suppose Nt photons transmitted, Nr photons received 


dN= -N΅(s)ds

NN (1–΅ds) 
photons photons 


ln(Nr/Nt) = .L –΅(s)ds 
ds 

• 
Object function is ΅(r), projection (Radon transform) is ln(Nr/Nt) 

• 
Projection easily estimated from transmitted and received X-ray intensities 

• 
How to recover object function? 







Reconstruction of object function 

[Kak and Slaney 1987, Ch 3.1] 
• 	
“Brute force” approach: sum backprojections fbp(x,y)= ..:[-.,.] p(t,.)d.; t= (x2+y2)1/2 

• 	
However, this causes blurring of object function 

• 	
If projections are appropriately filtered before being backprojected, reconstruction more accurate 

• 	
What is mathematical justification for filtering? 

• 	
What is optimal filter to use? (Hint: blurring effectively highlights lower frequency components) 



1,4,10,180 backprojections 



1,4,10,180 filtered backprojections 


Fourier slice theorem (aka projection slice theorem) [Meglicki 1999] 
• 
When taking the 1-D Fourier transform (FT) G of a function g, G(0) 

corresponds to the integral of g(x) 
G(0) = .x:[-.,.] g(x) e-2.j0x dx 



• 	
Similarly, when considering 2-D Fourier transform (FT) of object function: x:[-.,.] f(x,y) e-2.j(xX+yY)dx dy


F(X,Y)= .y:[-.,.].F(X,0)= .x:[-.,.] e-2.jxX (.y:[-.,.] f(x,y)dy)dx = .x:[-.,.] p(t=x,.=0) e-2.jxX dx= P(T,.=0) 
FT of object function at 0° (passing through origin) is FT of projection at 0° 
• 
More generally, the projection-slice theorem states that projection FT at . 

and slice of object function FT at . (running through origin) are equal P(T,.)= F(T,.) 

• 
How to use Fourier slice theorem to recover object function f ? 




Using Fourier slice theorem to recover f 

An idea: frequency-domain reconstruction 


Take FT of each projection 	“Place” projections on (X,Y) Fourier grid. Interpolate to find missing values. According to Fourier slice theorem, we are reconstructing FT of object function 


The idea is implementable, but there is an even more efficient method! 



Using Fourier slice theorem to recover f So far: 
f(x,y)= .X:[-.,.].Y:[-.,.] F(X,Y) e2.j(xX+yY)dX dY 
F(T,.)= P(T,.)= .t:[-.,.] p(.,t) e-2.jtT dt 
• 	How to convert F(X,Y) into F(T,.)? Use Jacobian determinant X = T cos.; Y = T sin.; 
J=|.X/.. .Y/.T-.X/.T .Y/..| = |-Tsin2. -Tcos2.|=|T| dX dY = |T| dT d. 
f(x,y)= ..:[-.,.]{.T:[-.,.] P(T,.)|T| e2.jtT dT}d. = ..:[-.,.] p’(t,.)d. where p’ is the ramp-filtered version of the projection 


Using Fourier slice theorem to recover f 
The simple relationship on previous slide shows an even simpler method to recover the object function f, termed the filtered backpropagation method: 
1. 	
Filter projections with ramp filter p’(t,.)= .T:[-.,.](. t:[-.,.] p(t,.) e-2.jtT dt) |T| e2.jtT dT 

2. 	
Backproject f(x,y)= . .:[-.,.] p’(t,.)d.; t= (x2+y2)1/2 


In essence, ramp filtering counteracts the blurring observed on the conventionally backfiltered projections 
Notes: 

• 
full-spectrum information may be missing, limiting the resolution 

• 
integration over 360° unnecessary: all information contained over 180° 






Ultrasound tomography – parameters of interest 

[Kak and Slaney 1999, Chapter 4.3] 
Attenuation 

• 
Arises in ultrasound as in X-rays, recovered similarly Attenuation slope with frequency 

• 
Higher contrast than attenuation image 


• 
Attenuation slope causes shift of central frequency of a Gaussian pulse 

• 
Frequency shift accumulative: proportional to integral of attenuation slope 

• 
Hence, attenuation slope may be recovered by generating attenuation images at 

several frequencies or by observing the frequency-shift Refractive index/speed of sound 


• 
Lower frequency than X-ray: propagation time/phase may be measured 

• 
Propagation time is line integral of refractive index (inverse of sound speed) Scatter 

• 
Inhomogeneities scatter ultrasound in many directions 




US tomography – challenges [Kak and Slaney 1999, Ch. 4.3] 
Wave distortion modifies the line integral; compared to X-rays: 
• 
greater changes in impedance: significant reflections 

• 
greater changes in propagation speed: significant refraction 

• 
longer wavelength: significant diffraction 



Difficulty of encircling the body with tomographic sensors 
• 
limited angles of access due to wave distortion (e.g. reflection from lungs) 

• 
significant attenuation from bone limits SNR for intervening tissue 

• 
coupling medium necessary between transducer and tissue 

• 
the future: wearable sensors? 




Solutions to challenges in US tomography 
1. Algebraic reconstruction [Kak and Slaney 1999, Ch. 7] 
• requires small and spatially slow variations in ., c (weakly refractive medium) 
2. Diffraction tomography [Kak and Slaney 1999, Ch. 6] 
• 
assume that scattering arises from incident field only (Born approximation) 

• 
less restrictive condition is Rytov approximation (beyond scope of lecture) 


3. Pulse-echo tomography [Cobbold 2007, pp. 556-557; Norton 1980] 
• also makes Born approximation but uses pulse rather than narrowband excitation 
• 	
Following discussion will focus on estimating scattering/reflectivity function 

• 	
However, using knowledge of the scattering/reflectivity function, the path of the incident wave may be estimated, which allows the expression of the other parameters as integral projections (see 2 slides back), and paves the way for their estimation 


See [Lavarello 2009] for more advanced methods e.g. with multiple scattering 

Algebraic reconstruction [Kak and Slaney 1999, Ch. 7] 
• 	
Object function is discretised to a grid of pixels f(x,y)= f(xi,yj) where (i,j) = round(x/.x,y/.y) 

• 	
Rays are also discretised according to the sampling available in the projection and assumed to be of finite width 

• 	
Each pixel’s contribution to line integral becomes weighted sum according to area of intersection of finite width ray with pixel 

• 	
Algebraic techniques used to solve system of linear equations 

• 	
One such technique is method of projections: iteratively project guess onto subsequent linear equation until convergence 

• 	
One technique to accommodate refraction (bent-ray tomography) is combining algebraic reconstruction with ray tracing in an iterative manner: generate system of linear equations (SLE) ignoring refraction; get image by solving SLE; trace rays through image to generate more accurate SLE; etc. 



Diffraction tomography[Kak and Slaney 1999, Ch. 6] 
• 	
Consider a parameter such as attenuation: in conventional straight-ray tomography, attenuation integrates over the ray path to yield a reduced received amplitude 

• 	
In the case of refracting boundaries, the line integral is modified to account for the bending of the rays 

• 	
In the case of wavelength or subwavelength scatterers, however, the concept of a ray (line parallel to planar wave, in direction of propagation) is invalidated 

• 	
In this case, it is easier to consider the distribution of scattering sources, and attenuation and refractive index are assumed constant (reducing to factors instead of integration) 

• 	
Is there a nice theorem that relates projections to the object function, similar to the Fourier slice theorem? 



straight-ray tomography 


bent-ray tomography 


diffraction tomography 





The Fourier diffraction theorem [Kak and Slaney 1999, Ch. 6.3] 
• 	Instead of a line (as in the Fourier slice theorem), the projection maps onto a semicircular slice of the object function FT 
forward scattered wave 


Adapted from [Kak and Slaney 1999, Ch 6.3, p. 219]. 



Fourier diffraction theorem: some observations 

[Kak and Slaney 1999, Ch. 6.3] 
• 	
The radius of the semicircle is the wavenumber k0 of the incident wave 

• 	
In the limiting case of k0›., the semicircle becomes a line crossing the origin, and we arrive at the diffractionless Fourier slice theorem 

• 	
Additional information is obtained by recording the backscatter, which projects onto the other half of the circle 

• 	
Again for k0›., the other semicircle becomes the same line through the origin 


(i.e. in the diffractionless case, forward and backprojection are identical) 




Fourier diffraction theorem: 
reconstruction of object function f(x,y) varying angle . 

[Kak and Slaney 1999, Ch. 6.4] 
• 	To cover the Fourier plane one can either 
. 
• 	
vary the insonation frequency (and thereby 
wavenumber k0) AND/OR 


• 	
vary the insonation angle . 



• 	As in diffractionless case, interpolation of projection values onto square grid is one k0 possibility 
• 	It turns out that this is simpler than the analogue of filtered backprojection, namely filtered backpropagation 
Adapted from [Kak and Slaney 1999, Ch 6.3, pp. 228-229]. 
varying wavenumber k0 




Pulse-echo tomography 
[Cobbold 2007, pp. 556-557; Norton 1980] 
• 	
Consider array of transceivers placed over a surface such as a cylinder 

• 	
Assume constant speed of sound c 

• 	
The two-way travel time of a pulse scattered at a distance d is 2d/c 

• 	
Therefore, if a cylindrical wavefront (pulse) is transmitted at t=0, the signal received at t=2d/c is a line integral of reflectivity over a cylinder of radius d 

• 	
Solution: decompose pulse-echo signal p(.,t) into circular harmonics and carry out a transformation that yields the circular harmonics of the object function 




Tomography in 3-D [Kak and Slaney 1999, pp. 99-107, 219] 
• 	
So far, we have assumed 2-D object functions (or 3-D object functions invariant in one of the dimensions) 

• 
What about the third – let us call it z – axis? Two approaches: 

• 	
Focus beam tightly in the z axis, i.e. over a “slice” of tissue 

• 	
Use one of the 2-D reconstruction algorithms discussed so far 

• 	
In the case of a 1-D array, the array will need to be mechanically scanned to cover the entire region of interest 

• 	
Averaging of tissue properties will occur over slice 


OR 


• 	
Extend 2-D reconstruction algorithm to 3-D 

• 	
Simple extension in straight-ray tomography. Receiver becomes plane and results on one value of z are unaffected by the other (decoupled) due to parallel rays 

• 	
Fourier diffraction theorem, pulse-echo tomography: circle becomes sphere 




The breasts: a viable subject of US tomography 
Breasts are an excellent organ to image with ultrasound tomography: 
• 	
easy to surround with gel-coupled sensors 

• 	
soft tissue with little refraction 

• 
no bones or large air pockets to prevent penetration of ultrasound Growing commercial interest in ultrasound tomography of the breast: 

• 	
Techniscan (look it up on Google and Youtube) solves inverse scattering problem using graphical processing units (GPUs) 

• 	
Delphinus Medical technologies is a spin-off of research [Li et al. 2009] (here discussing bent-ray tomography) at Karmanos Cancer Institute, Wayne State University 




[Cobbold 2007] Foundations of biomedical ultrasound 
[Kak and Slaney 1987] Principles of computerized tomography. 

http://www.slaney.org/pct/pct-toc.html 

[Lavarello 2009] New developments on quantitative imaging using ultrasonic waves. 
http://www.brl.uiuc.edu/Downloads/lavarello/Lavarello_PhD.pdf 
[Meglicki 1999] Lecture notes on Advanced Scientific Computing. 
http://beige.ucs.indiana.edu/B673/node19.html 

[Norton 1980] Reconstruction of a two-dimensional reflecting medium over a circular domain: Exact solution