Medical diagnostic systems 

(Orvosbiológiai képalkotó rendszerek) 

Modelling ultrasound image formation 

(Ultrahang képalkotás modellezése) 
Miklós Gyöngy 



Quick reminder #1: linear time-invariant systems Suppose 
• 
System w(t)›a(t) is linear time-invariant 

• 
.w1(t)+ßw2(t) ›.a1(t)+ßa2(t) 

• 
w(t-.) › a(t-.) 



• 
.(t) ›h(t) (impulse response) 


• 
Fourier transform F{h(t)}= H(f) (frequency response) Then 

• 
a(t)=h(t)*tw(t)=..:[-.,+.] h(.)w(t-.)d. (temporal convolution) 


• A(f)=H(f)W(f) 
Model ultrasonic image formation as linear time-invariant? 




Quick reminder #2: point spread functions 
• 	
Point spread function P is 2-D equivalent of impulse response 

• 	
Assuming P is spatially invariant, F 2D{P*2DT}= F 2D{P}F 2D{T} 

• 	
Typically, P is Gaussian-type blurring function (or low-pass filter) 

• 	
Arises in optics, when structures in original scene are larger than wavelength 

• 	
For sub-wavelength structures illuminated by coherent light or ultrasound, interference of scatterered radiation occurs, resulting in speckle 


"ultrasonic" image original scene "optic" psf (x10) "optic" image "ultrasonic" psf (x10) "ultrasonic" image envelope 

Original photograph in public domain. Credit: Victor B. Scheffer / US Fish and Wildlife Service (www.fws.gov/digitalmedia) 



• 
input waveform temporally filtered by electromechanical response 

• 
beamforming operations temporal filters that act as spatial filters on medium 

• 
above create point spread function (P) of imaging system 

• 
assume small region of interest (ROI) where P is spatially invariant 

• 
beamformed image B convolution of tissue impulse response T with P 

• 
envelope operation not linear (no psf) but there is a characteristic beam size 

• 
scattering within beam (resolution cell) leads to interference effects 






Three approaches to calculating point spread function P 
1. 	
Point element (PE) approximation: array elements as point transceivers 

2. 	
Spatial impulse response (SIR): convolution of impulse response with input waveform yields actual response 

3. 	
Frequency-domain (FD) analysis: consider transmit-receive beam shape at carrier frequency of input waveform and multiply by input waveform 


Calculating beamformed image P * T = B • origin of tissue impulse response T 

• 
models for T 


• 
image speckle – junk or disguised information? 
Brief notes on quantitative ultrasound (QUS) and simulation software 






• 
N point transceivers with uniform frequency response 

• 
Define rab := |rb-ra|; .ab := rab/c 

• 
Pressure wavefronts focussed to rf 

• 	
Total incident pressure at ra sum of point sources at rn p(r,t)=.A/rw(t-(.


annna na -.0a)+.nf) 
where An are apodization weights 

0

• 
Monopolar scatterer at rs with scattering strength ss 

• 	
Scattered pressure recorded at rn p’(r,t)= s/rp(r, t -.)


nssn ssnrn
• 	Beamformed image (using two-way travel time) B(r)= .Arp’(r,2.)
annna nna

• 	Sanity check: r = r{0,a,f,s}; .n = .{0,a,f,s}n; An = 1/N B(r) = 1/N2 ..w(2.-.-.+ .-.)= w(0) 
nn nnn nn





N = 128; D = 19 mm; c = 1500 m/s; f0 = 7.5 MHz; rf = (0,0,75 mm); w(t) = cos(2.f0t) {cos(.f0t/3)rect(f0t/3)} (3 cycles of cosine-modulated 7.5 MHz) 




PE approximation: validity (or otherwise) 
• 
Typical array element: D{x,y} = {0.3 mm, 5 mm}, f0 = 7.5 MHz 

• 
Far-field distance: D{x,y} 2f/4c= 0.1 mm, 31 mm [Olympus 2006] 

• 
In practice, far-field usually brought even closer due to elevational focusing 

• 
By using more elements along elevation direction to model transducer, 

• 
In limit of N›. 

– 
array becomes continuous surface 

– 
sum becomes integration 

– 
expression reflects Huygen’s principle 

– 
analytical solutions sometimes possible 





[Gyöngy 2010] 




Huygen’s principle 

• 	
Pressure from integral of contributions across surface 

• 	
Let acceleration dv(r’)/dt of surface be A(r’)w(t) 

• 	
Then, Rayleigh integral [Jensen 1999] becomes 







p(r,t)=(.0/2.).A(r’)/|r-r’| w(t-|r-r’|/c-.(r’)) dr’ 
• 	
NB: using form of Rayleigh integral without angle-dependent obliquity factor [Goodman 1996, p. 51] 

• 	
Pressure response h1(r,t) to w(t)=.(t)? (SIR) h1(r,t)=(.0/2.).A(r’)/|r-r’| .(t-|r-r’|/c-.(r’)) dr’ 

• 	
Pressure response H1(r,f) to w(t)=exp(j2.ft)? (FD) H1(r,f)=(.0/2.).A(r’)/|r-r’| exp[j2.f(t-|r-r’|/c-.(r’)] dr’ 





Spatial impulse response (SIR) [Jensen 1999] 

• 	
Pressure response h1(r,t) to w(t)=.(t) h1(r,t)=(.0/2.).A(r’)/|r-r’| .(t-|r-r’|/c-.(r’)) dr’ 

• 	
Pressure response to arbitrary w(t) p(r,t)= h1(r,t)*tw(t) 

• 	
h1(r,t) known as SIR: response p(r,t) due to transducer excited by impulse 


• 	
What does spatial impulse response look like? Assume focal point rf and time of arrival tf 


• at point of focal formation: h1(rf ,t) ..(t-tf) 
• outside it, response is blurred by imperfect focusing 

Next: consider SIR of three transducer shapes (all unfocussed and unapodised) 





= .0c[H(t-z/c)-H(t-(a2+ z2)1/2/c)] 

h
1([0,0,z],t) 
• 
For each R, equal contribution to final waveform! .0c 

• 
Integration of equally-weighted . yields rect function 


(a2+ z2)1/2/c

z/ct
• For impulse velocity (rather than acceleration) across 

.0c

surface, two opposing . are generated: one from centre 
(direct wave) and one from edge (edge wave) 



-.0c 
(a2
z2)1/2/c
z/c 

t 

11 



Example 2: Far-field SIR of rectangular transducer [Ellis et al. 2007] 

• 	
Define propagation time values t1,t2,t3,t4 from corners of rectangle to r (in ascending order) 

• 	
Far-field approximations 


1. 
no wavelet from rectangle can reach r until t1 (0 response) 

2. 
no wavelet will reach r after t4 (0 response) 

3. 
maximum and uniform response between t3 and t4 

4. 	
linear increase and decrease during intermediate times 





h1(r,t) 

• 	
Above approximations lead to a trapezoid response 

• 	
Linear arrays consist of rectangular elements (though elevational focusing lens is placed in front of them) 

• 	
Arbitrary shapes decomposed into rectangular elements 


t1 t2 t3 t4 t 





Example 3: Polygon [Jensen 1999] 

• 	
Acoustic reciprocity: pressure recorded by point receiver B due to point source A remains the same if locations of source and receiver are interchanged 

• 	
More generally, if source consists of sum/integration of point sources over surface A, pressure recorded at B is equivalent to sum/integration of pressures over A caused by point source at B 

• 	
Intersection of spherically spreading wave with polygon is set of circular arcs whose total length determines magnitude of pressure received 

• 	
Validity of acoustic reciprocity: arcs form region of polygon that determines p(r,t) 

• 	
Arbitrary shapes decomposed into polygons 







h(t) p(r,t)=w(t)*h(t) h(t)*h(t) psf=w(t)*h(t)*h(t) 
waveformw(t) 1st11t1t1t1
1 1 

50 

15 


4 10 
2 
0 

0.5 
0 5 0 

-1 0 -50 0 
-2 
-10 0 10 30 40 50 60 30 40 50 60 80 100 80 100 

time (µs) time (µs) time (µs) time (µs) time (µs) 





Frequency-domain (FD) method 

• 	
Pressure response H1(r,f) to w(t)=exp(j2.ft) H1(r,f)=(.0/2.).A(r’)/|r-r’| exp[j2.f(t-|r-r’|/c-.(r’)] dr’ 

• 	
Pressure in focal plane is Fourier transform of A(r’)! 

• 	
Proof: 


• 	
Assume transducer is on x-y plane 

• 	
Focusing at t=0, r=(0,0,zf), .(x’,y’)=-(x’2+y’2+z’2)1/2/c 

• 	
x’,y’<<zf; approximation 1/|r-r’|.1/zf 


• 	
x’,y’ << zf; Fresnel approximation (zf2+x’2+y’2)1/2.zf+(x’2+y’2)/2zf 

• 	
H1(x,y,zf) = (.0/2.zf) ..A(x’,y’) exp[j(.t-k.)] dx’dy’ 


• 	
.=((x-x’)2+(y-y’)2+zf)1/2 -(x’2+y’2+z’2)1/2 .(x2-2xx’+y2-2yy’)/2zf 

• 	
|H1(x,y,zf)| = (.0/2.zf) ..A(x’,y’) exp[jk(xx’+yy’) dx’dy’ 









FD method 

• 	
Pressure in focal plane is 2-D Fourier transform of A(r’)! 

• 	
Circular aperture: Airy patterns (cf. optical lenses) 




grating lobe 

• Uniform apodization across rectangular aperture Dx, Dy 
• 	
|H1(x,0,zf)| = (.0/2.zf) sinc(x/f#:x.) sinc(y/f#:x.) 

• 
f#:x = zf/Dx,; f#:x = zf/Dx,(f-number in x and y directions) 


sidelobes 

• 	sinc(±0.443) .1/.2; -3dB beam width (one-way!) is .0.89f#. 
main lobe 

• 	
Increasing D, or decreasing zf or ., all narrow the beam 


• 	
Just as in spectrum analysis, windowing function (typically) 
reduces sidelobes at cost of increasing main lobe width 



• 	Electronic steering of discrete elements 

• 	
Fourier theory: multiplication in space - convolution in frequency 

• 	
If discrete elements modelled as multiplication with dirac deltas 
(point sources as before), this creates grating lobes 







FD: central beam dimensions 

• 
Recap: beam in focal plane 2-D Fourier transform of aperture function A(x,y) 

• 	
One-way beam in focal plane (x,y,zf) for uniform aperture: |H1(x,y,zf)| = (.0/2.zf) sinc(x/f#:x.) sinc(y/f#:y.) 

• 
Beam along z-axis? 	Fresnel transform of A(x,y)! [Mast 2007; Gyöngy 2010] 

• 	
Fresnel integral: Fr(z)= .l:[0,z] exp(j.l2/2) dl 

• 	
In analogy with sinc function, define Fresnc function [Gyöngy 2010]: Fc(z) = Fr(z2)/z2 where z2 = (z/2)1/2 

• 	
Beam in axial direction is approximated by larger aperture dimension: |H1(z,zf)| = Fc(.z/f#2.); f# = min(f#:x, f#:y) 





FD: Gaussian approximation to central beam 

• 
Reciprocity theorem: transmit and receive beams identical 

• 	
One-way (transmit or receive) -3dB beam widths along x, y, z: (.x3, .y3, .z3) = (0.89f#:x., 0.89f#:y., 6.95f#2.) 

• 	
3-D Gaussian function with -3dB beam widths .x3, .y3, .z3: exp{(0.59x/ .x3)2 + (0.59y/ .y3)2 + (0.59z/ .z3)2} 

• 	
Thus, beam can be approximated by Gaussian function exp{-(0.66x/f#:x.)2 -(0.66y/f#:y.)2 -(0.085z/f#2.)2} 

• 	
Two way (transmit-receive) exp{-2(0.66x/f#:x.)2 -2(0.66y/f#:y.)2 -2(0.085z/f#2.)2} 





FD: estimating P for pulsed waveform w(t) 74 
• Typically w(t)/c<< .
z3 75 


• P in axial direction determined by w(t) 

• P in transverse direction determined by 

CW beam of carrier frequency 

• Thus, estimate P as product of CW 

beam by waveform 

axial direction (mm)

76 -5 0 

psf 
74 


75 

76 -5 

05 
envelope of psf 
• Compare result with PE method 74 (simulation parameters identical) 75 


76 
transverse direction (mm) 



Summary of different approaches to calculating P 

Point element (PE) approximation 

• 
Arbitrary accuracy achievable 

• 
Simple calculation 


• 
Large number of points may make computation slow Spatial impulse response (SIR) 

• 
Often simple analytical expression exists for SIR 


• 
SIR specific to imaging array only: no recalculation needed if w(t) changed 

• 
Large sampling rate needed to make convolution accurate Frequency-domain (FD) analysis 

• 
Conceptually simple 

• 
Inaccurate extension to pulsed waveform 






Calculating beamformed image 

• 	
Simple shift-invariant model: beamformed image B=P *2D T 

• 
B-mode image obtained by taking envelope of B 

• 
We have found models to estimate P 


• 
Where does tissue impulse response T come from? 

• 
How do we model it? 


• 
What is image speckle – junk or disguised information? 



B-mode image of bovine liver sample in water 





Origin of tissue impulse response T: scatterers 
A useful classification method for scatterers: 
[Szabo 2004, pp. 244-245; Greenleaf and Sehgal 1992] 
Class 0: << . molecules (absorption, e.g. water, protein) Class 1: <. scatterers (speckle, e.g. cells, e.cellular matrix) Class 2: ~. scatterers (distinguishable scatterer, e.g. bubble) Class 3: >. specular boundary (e.g. organ boundary) 
Consider three scattering models: 
1. 
discrete (point or with point origin) scatterer 

2. 
inhomogeneous continuum 

3. 
discrete-continuum hybrid 


B-mode image of 
bovine liver 
sample in water 




Point scatterer model 

• 	
Use P derived from FD method 

• 	
Projection of scatterers onto 2-D imaging plane 


• 	
Here, using point scatterers with same scattering strength 

• 	
As concentration of scatterers increases, packing constraints make scatterers more regular [Hete and Shung 1993] 

• 	
Instead of making location of scatterers totally random 

• 	
have 10 µm × 10 µm pixels (typical cell size ~10 µm) 


• 	
specify Cs: average number of scatterers per resolution cell 

• 	
calculate probability of pixel containing a scatterer 

• 	
sum contributions from each scatterer 




• 	
How does varying scatterer concentration Cs change scattering properties? 

• 	
qualitative appearance of B-mode image 

• 	
histogram of B-mode values 


• 	
consider power spectrum of average of axial auto-correlations of beamformed image B 







Point scatterer model: Cs = 0.01 




amplitude frequency

1010  histogram (I)  
 autocorrelation PSD 0.02axial auto-correlation power spectrum (B )  
105  
 0.01  
 is therefore simply square of power spectrum of pulse  
100 0  0.2  0.4  0.6  0.8  1  0 0  5  10  15  20  
scattering strength  spatial frequency (1/mm)  




Point scatterer model: Cs = 0.1 




beamformed image B B-mode image I 

distinguished from 
-5 05 


-5 0 5 

each otherx 
transverse direction (mm) transverse direction (mm) 


histogram (I) axial auto-correlation power spectrum (B )
6
10
0.4 
amplitude frequency
1010



as with previous 
histogram, most autocorrelation PSD 
however contains 
4

pixels have a values 
0.2 
peaks from 

near 0, with 1 only 
at the scatterer interference 

between scatterers 
2

locations 
0 
scattering strength spatial frequency (1/mm) 




Point scatterer model: Cs=1 




-5 05 


-5 0 5 

transverse direction (mm) transverse direction (mm) 



beamformed image B B-mode image I 

-5 05 


-5 0 5 

from each other 
transverse direction (mm) transverse direction (mm) 

scattering strength spatial frequency (1/mm) the psf 



Point scatterer model: Cs =10 



•scatterers at distances 

point spread function P 

tissue impulse response T 
cycles of transmitted 
-5 05 


-5 0 5 

pulse 
transverse direction (mm) transverse direction (mm) 

•destructive interference can now 

in histogram 
•Rayleigh-like distribution 

histogram (I) axial auto-correlation power spectrum (B ) 

02468 

30 

20 

10 



a DC component •resolution cell 
with average of 
10 scatterers 
means that most 
values are no 
longer near 0 
has appeared 
0 

scattering strength spatial frequency (1/mm) 






Medical diagnostic systems – Modelling ultrasound image formation 
www.itk.ppke.hu 
Point scatterer model: Cs = 100 
axial direction (mm)axial direction (mm)
axial direction (mm)axial direction (mm) 

point spread function P tissue impulse response T

74 

74 

75 

76 

75 

76 

-5 05 

transverse direction (mm) 


beamformed image B 

-5 0 5 

transverse direction (mm) 

B-mode image I 
74 

75 

76 

74 

75 

76 

-5 05 

transverse direction (mm) 

-5 0 5 

transverse direction (mm) 
•distribution still 

histogram (I) axial auto-correlation power spectrum (B ) 
10000 
1000 
the scatterers are 
beginning to be 

significantly constrained 
5000 
500 

by the 10µm×10µm grid, 
as evidenced by large Rayleigh-like 
•increasing 
number of 
scatterers by 10 
only caused 
amplitude frequency

modest increase 
in mean! 
0  0  DC component  
0  5  10  15  20  0  5  10  15  20  
scattering strength  spatial frequency (1/mm)  

28 



Medical diagnostic systems – Modelling ultrasound image formation 
www.itk.ppke.hu 
Point scatterer model: Cs = 1000 
axial direction (mm)axial direction (mm)
axial direction (mm)axial direction (mm) 
point spread function 
P 



tissue impulse response 
T 

74 

75 

76 

74 

75 

76 

-5 05 

transverse direction (mm) 


beamformed image 
B 


-5 0 5 

transverse direction (mm) 

B-mode (envelope of image) I 
74 

75 

76 

74 

75 

76 

-5 05 

transverse direction (mm) 

-5 0 5 
image dominated 
transverse direction (mm) 
by DC term 


B-mode histogram (I) axial auto-correlation power spectrum (B ) 
10000 
10e4 
shape of non-DC 

component very similar addition of DC 
5000  5e4  to single-scatterer case:  
cf van Cittert-Zernike  
theorem  
0 0  20  40  60  0 0  5  10  15  20 [Anderson and  
scattering strength  spatial frequency (1/mm)  Trahey 2006]  

component creates 
Rician distribution 
amplitude frequency
29 




Medical diagnostic systems – Modelling ultrasound image formation 
Point scatterer model: explanation of observations 

• 
Apart from DC term at high Cs, auto-correlation retains shape of P Cs ~ 1 and lower 

• 
Scatterers do not interfere: histogram of amplitudes reflects histogram of P Cs ~ 10-100 

• 	
Phase of backscatter components independent of each other 

• 	
Total i (in-phase) and q (quadrature) components arise from addition of individual ii and qi components (assumed independent) from each scatterer si 

• 	
Central Limit Theorem: i and q have mean=0 Gaussian distributions 

• 
Backscatter envelope |I| =(i2+q2)1/2 is therefore Rayleigh-distributed Cs ~ 1000 and higher 

• 	
10 µm × 10 µm pixel grid constraint: ii, qi no longer independent 

• 	
scattering has coherent component 


• 	
mean.0 Gaussian i, q: |I| is Rician distributed [Anderson and Trahey 2006] 






Medical diagnostic systems – Modelling ultrasound image formation 
Inhomogeneous continuum model 
• 	
In limit of mean scatterer concentration Cs›., scattering strength=1/Cs 

• 	
spatial distribution of scatterers becomes continuous 


• 	
equivalent to inhomogeneous continuum model with small compressibility changes 

• 	
(Note: density changes are usually neglected due to their angle-dependent scatter, although they can also be readily incorporated assuming backscatter – see earlier lecture, “Fundamental Concepts in Acoustics”, slide title “Sub-wavelength scattering”) 



• 	
Suppose medium completely characterised by spatial autocorrelation R. of its fractional changes in compressibility, .=(.-.0)/. 

• 	
B-mode spatial autocorrelation function RI = R. * r P 

• 	
Random distribution of point or <<. particles (R. . .(r)): RI = P 


[Anderson&Trahey 2006;Hill et al. 2004, p. 217;Mallart&Fink 1991;Wagner et al. 1983] 
• 	RI only dependent on imaging system! (cf “Point scatterer model: Cs = 1000”) 
• 	Other R. tissue models: fluid sphere, spherical shell, Gauss [Insana et al. 1991] 




Medical diagnostic systems – Modelling ultrasound image formation 
Inhomogeneous continuum model: 
its use in quantitative ultrasound (QUS) 

• 	
By choosing a tissue model for autocorrelation R., one can predict the spatial autocorrelation of the image RI 

• 	
If RI is taken along the imaging axis (A-lines), scaled 2/c to get the temporal autocorrelation and Fourier transformed, we obtain the frequency response of the scattering medium [Oelze et al. 2002], which is proportional to the oft-used term form factor (see [Insana et al. 1990] for definition) 

• 	
The chosen tissue model is parameterised by the scatterer diameter ar (or, in the case of a Gaussian auto-correlation, effective diameter) 

• 	
By fitting the frequency response over a gated depth to the theoretical frequency response, the (effective) scatterer diameter may be estimated 

• 	
Example 1: slope of response using Gaussian model finds ar, so identification of several slopes reveals scattering sources in kidney [Insana et al. 1991] 

• 	
Example 2: fit resonance peaks of spherical particles to Faran model to estimate size 


[Condliffe 2009] 



Hybrid scattering model 

[Mo and Cobbold 1992; Lim et al. 1996; Cobbold 2007, pp. 315-319] 
• 	
Consider voxels <./10 where phase is nearly constant throughout 

• 	
Scatterers in each voxel “add up” to create single equivalent scatterer 

• 	
Volume of equivalent scatterer equal to sum of volumes 

• 	
Centre of scatterer equal to centre of phase (for equivalent scatterers, this equals centre of mass [Lim et al. 1996]) 


• 


2-D representation of hybrid model. Left: original scatterers. Right: equivalent scatterers for each pixel. Adapted from [Lim et al. 1996]. 




Medical diagnostic systems – Modelling ultrasound image formation 
Summary of lecture 

• 
Considered methods to estimate point spread function P 

• 
Considered several tissue models 




How realistic are these models? 

• 
Hexagonal close packing in liver [Doyle 2009] gives coherent speckle 

• 
Some tissues agree well with continuum autocorr. models [Insana et al. 1991] 

• 
Blood agrees well with hybrid model [Hill 2004, p. 213] 


Effects we have neglected: 

• 
Reflection (e.g. organ boundaries). Need anatomic info [Dillenseger et al. 2009] 

• 
Attenuation. Simple to incorporate. TGC partially compensates for it. 

• 
Instrumentation noise. More significant at larger depths due to attenuation 

• 
Spatial invariance of P. SIR able to model it; dynamic apodization reduces it. 

• 
Non-linear propagation of sound 

• 
Multiple scattering 


• 
Abersim http://www.ntnu.no/isb/abersim non-linear acoustics! 

• 
Field II http://server.electro.dtu.dk/personal/jaj/field/) 

• 
Ultrasim http://heim.ifi.uio.no/~ultrasim/ 


• 
Delfi http://www.mathworks.com/matlabcentral/fileexchange 

• 
DREAM http://www.signal.uu.se/Toolbox/dream/ 




34 





[Anderson and Trahey 2006] A seminar on k-space applied to medical ultrasound. 
http://dukemil.bme.duke.edu/Ultrasound/k-space/index.htm 
[Cobbold 2007] Foundations of biomedical ultrasound [Condliffe 2009] Particle characterization by acoustic microscopy following needle-free injection [Dillenseger et al. 2009] Fast simulation of ultrasound images from a CT volume [Doyle et al. 2009] Simulation of elastic wave scattering in cells and tissues at the 
microscopic level [Ellis et al. 2007] A spline-based approach for computing spatial impulse responses [Greenleaf and Sehgal 1992] Biologic system evaluation with ultrasound [Gyöngy 2010] Passive cavitation mapping for monitoring ultrasound therapy [Hete and Shung 1993] Scattering of ultrasound from skeletal muscle tissue [Hill et al. 2004] Physical principles of medical ultrasonics [Insana et al. 1990] Describing small-scale structure in random media using pulse-echo 
ultrasound 

... ... 



[Insana et al . 1991] Identifying acoustic scattering sources in normal renal parenchyma from the anisotropy in acoustic properties 
[Jensen 1999] Linear description of ultrasound imaging systems. 
http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ultra12/ref_jaj_1999.pdf 
[Lim et al. 1996] Particle and voxel approaches for simulating ultrasound backscattering from tissue 
[Mallart and Fink 1991] The van Cittert--Zernike theorem in pulse echo measurements 
[Mast 2007] Fresnel approximations for acoustic fields of rectangularly symmetric sources 
[Mo and Cobbold 1992] A unified approach to modeling the backscattered Doppler ultrasound from blood 
[Oelze et al. 2002] Ultrasonic characterization of tissue microstructure 
[Olympus 2006] Ultrasonic transducers technical notes. http://www.olympus-ims.com/data/File/panametrics/UT-technotes.en.pdf 
[Wagner et al. 1983] Statistics of speckle in ultrasound B-Scans