Skip to main content

Full text of "DTIC ADA193042: Methods for Correction of Refractive Errors."

See other formats


• >p; j . 1 a i r i » i»v a i j a n 11 •{ • j 3 h 1 11 ■ m ■ ! j < a a •; • • m i*i j j { i \-wi i j i -m ji 4|WV|: 


UNIV LAFAVETTE IN SCHOOL OF ELECTRtCAL ENGINEERING 
A C KAN 11 DEC 84 DAMD17-82-C-2819 

UNCLASSIFIED F/6 28/1 NL 





































- r .-\ j ', / --. 

V--.v 

v '/■ .* v* v v “/- v;/ ■/ -'V 































AD-A193 04: 


U 

,« 




% “w 


in: 








* * - 


rv - 1 




I.', 


.v 


& 




t 


AD. 


fSJ 




DTiS.HlE com 


METHODS FOR CORRECTION OF REFRACTIVE ERRORS 


A. C. Kak 


Annual Report 

December 31, 1984 


Supported by 


U. S. Army Medical Research and Development Command 
Fort Detrick, Frederick, Maryland 2170l‘ i *5012 


e* 

4e*»' 


Contract Number DAMD17-82-C-2019 


School of Electrical Engineering 
Purdue University 
West Lafavette. IN 47907 


Approved for public release: 
Distribution Unlimited 



The findings of this report are not to be construed as an 
official Department of the Army position unless so designated 
bv other authorized documents. 


8 8 4 1 12 5 


-'.AA.VA-V 


















SEC^RO' C_ASSiF CAPON Of This PAGE 


REPORT DOCUMENTATION PAGE 


lb RESTRICTIVE MARKINGS 


Form Approved 
OMB No 0704 0188 
Exp Dare lun 30 1986 


'a REPOR T SECURITY CLASSIFICATION 

UNCLASSIFIED 


La SECURITY CLASSIFICATION AUTHORITY 


Lb DECLASSIFICATION DOWNGRADING SCHEDULE 


4 PERFORMING ORGANIZATION REPORT NUMBER(S) 


6a NAME OF PERFORMING ORGANIZATION 6b OFFICE SYM80L 

School of Electrical Engineering (if applicable) 
Purdue University 


6c. ADDRESS (City, State, and ZIP Code) 

vest Lafavette, IX 47907 


3 distribution - availability of report 

Approved tor public release; 
distribution unlimited 


5 MONITORING ORGANIZATION RE p ORT NUM8ER(S) 


7a NAME OF MONITORING ORGANIZATION 


7b ADDRESS (City. State, and ZIP Code) 


8a. NAME OF FUNDING / SPONSORING 

organization U.S. Armv Medical 
Research 4 Development Command 


8c. ADDRESS (City, State, and ZIP Code) 

Fort Detrick 

Frederick, MD 21701-5011 


8b OFFICE SYM30L 9 PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 
(if applicable) 

SU.RD-RMI-S DAMD1 7-82-C-2019 


10 SOURCE OF FUNDING NUMBERS 


PROGRAM 
ELEMENT NO 

62777A 


PROJECT 

task 

WORK UNIT 

NO 3 E L — 

NO 

ACCESSION NO 

62777A878 

BB 

017 


ll TITLE (Include Security Classification) 

(U) Methods for Correction of Refractive Errors 

> 

- 

12 personal AUTHOR(S) 

Kak, A.C. 




13a TYPE OF REPORT 

13b TIME COVERED 

14 DATE OF REPORT (Year, Month, Pay) 

15 PAGE COUNT 

Annua 1 

from 1 / 84 to 12/84 

1984 December 31 

226 

16 SUPPLEMENTARY NOTATION 





17 

COSATI 

CODES 

FIELD 

1 GROUP 

| SUB-GROUP 


18 SUBJECT terms ( Continue on reverse if necessary and identify by block number) 

Biophysics; Bioeffects; Dosimetry, Microwaves 


19 ABSTRACT ( Continue on reverse if necessary and identify by block number) 

The problem of cross sectional (tomographic) imagine of objects with diffracting sources 
is addressed. Specifically the area of investigation is the effect of multiple scattering 
ar.u attenuation phenomena in diffraction imaging. 

The validity of either the Born or the Rytov approximations is the basic assumption 
behind all the inverse scattering techniques in diffraction tomography. To test these 
techniques '.•/hen these assumnt ions a re not satisfied, we have developed a computational 


:>r v 

edure f 

■.' r 

t 

n ca 1 c 

ul 

it: inn of 

tiie 

"true" 

sea- 

ttered fields fre 

'm a mul 

ti-component obhot . 

7T1 :i 


pm 

’e 

d 1 i r ■- , 

til 

e : n r i o r 

ma no 

O O f t \v 

:o av. 

a i 1 a b 1 e d i f f r a e t i 

.on reco 

nstruct ion teciin Lques 

i -> e 

i" i n 

i r. 

r 

a m 

se 

na- nf m 

:u i t i 

vie sea 

tter 

ing effects. The 

1 Simula 

tion results shew 

tile 

- 1 1 • r in 

r i c 


of tile 

S 

vnt'llet i r 

Age 

rture t 

echn 

ique. We have al 

SO Stud 

ied the roK of 

If. 

nun r_ i ■ >:i 

i n 

r 

Uu l\'i' 

on 

st ruct io 

n t o 

<’ r i n l nuo 

s. 

To calculate the 

scat ter 

eci fields from an 

■ ' 0 * L' 

ct in t 

: ; t. • 

p r 

e Seilfe 

o 

f attenu 

at in 

n, new 

com pi 

liter simulation ; 

1 rocrams 

are developed. 

lo > 

S'p B O' 

N A v 

■An 

_ABi Li 

OF 

ABS'RAC* 




Ll A3STRAC' SECu° 'Y 

C.ASS -;CA 

TON 

□ 

yNC.ASS'F 

ED u 

NL 

MlTEO 

Q 

SAVE AS 3 

0 ' 

□ ore ' 

jSERS 

UNCLASSIFIED 



223 ' 


ESPO 

NS 

BlE ND 

V D 





:;b ’EwE d hQ\E ( Include 

Area Code) 

2Lc O c f CE S v V30. 

M, i r v 

; ; ranntu 

s B 

ns 

t i a n 

_ 





301/66323 


SGRD-RMI-S 


00 FORM 1473, 84 mar 


rc r or ^ay oe ^seo exhausted 


$Er _R *Y C.ASS oca'ON Of 










































































TABLE OF CONTENTS 


Page 

LIST OF FIGURES.iv 

ABSTRACT.xxiv 


CHAPTER I. INTRODUCTION. 

List of References. 

CHAPTER 2. FORWARD PROBLEM. 

2.1 Propagation of Angular Spectrum. 

2.2 Scattering from a Single Object. 

2.2.1 Plane Wave Illumination. 

2.2.2 Point Source Illumination. 

2.2.3 Numerical Calculation of Scattered Fields. 

2.3 Scattering from a Multi-component Object. 

2.3.1 Review of the Basic Formulation of Multiple Scattering 

2.3.2 A New Procedure for Computing Multiple Scattering ... 

2.3.3 Numerical Results on the Calculation of Multiple 

Scattering. 

2.3.4 The New Procedure in Attenuating Medium. 

2.4 Test of Numerical Results. 

List of References. 


1 

5 


:io 

.12 

13 

.16 

18 

24 

30 

34 


39 

46 

54 

62 


CHAPTER 3. INVERSE PROBLEM.64 

3.1 Formulation of the Problem.66 

3.1.1 Interaction of Wave with Tissue.66 

3.1.2 An Integral Solution.68 

3.1.3 Born Approximation.70 

3.1.4 Rytov Approximation.70 

3.2 Conventional Diffraction Tomography.71 

3.3 Synthetic Aperture Diffraction Tomography.78 

3.3.1 Transmission Mode.81 

3.3.2 Reflection Mode .87 

3.3.3 Lossy Medium .98 

List of References.103 

































































Pago 


CHAPTER 4. A SIMULATION STUDY OF MULTIPLE 
SCATTERING AND ATTENUATION 

IN DIFFRACTION TOMOGRAPHY'.106 

4.1 Artifacts Caused by Multiple Scattering Phenomena.106 

4.1.1 Artifacts in Conventional Diffraction Tomography.107 

4.1.2 Different Approaches to Reduce the Artifacts.112 

4.1.2.1 First Approach.115 

4.1.2.2 Second Approach.117 

4.1.2.3 Third Approach.118 

4.1.2.4 Simulation Results and Discussion.118 

4.1.3 Artifacts in Synthetic Aperture Diffraction Tomography.131 

4.2 A Simulation Study of Attenuation in Diffraction Tomography.159 

4.2.1 A Simplified View of the Problem.160 

4.2.2 Attenuation in Conventional Diffraction Tomography.164 

4.2.3 Attenuation in Synthetic Aperture Diffraction 

Tomography.17§ 

List of References. I8f 

.. v 

«*- 

CHAPTER 5. SUMMARY.189 

BIBLIOGRAPHY.195 


.APPENDIX 


200 




















LIST OF FIGURES 


Figure Page 

2.1 The scattered field from a cylinder of radius a, illuminated by 

a plane wave, is measured at the point (r,0).14 

2.2 The scattered field from a cylinder of radius a, illuminated by a. 

point source at (r o ,0 o ), is measured at the point (r,0).17 


2.3 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.1 (embedded in a lossless medium) illuminated 

by a plane wave, measured on the receiver line at a distance of 10 X 
Irom the cylinder. The solid line corresponds to the case of a lossless - 
cylinder, while the dashed line represents the field when the 
attenuation constant of the cylinder is 3.5 db/cm.21 

2.4 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.5 (embedded in a lossless medium) illuminated 

by a plane wave, measured on the receiver line at a distance of 10 X 
from the cylinder. The solid line corresponds to the case of a lossless 
cylinder, while the dashed line represents the field when the 
attenuation constant of the cylinder is 3.5 db/cm.22 

2.5 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.1 illuminated by a plane wave, measured on 
the receiver line at a distance of 10 X from the cylinder. 

The solid line corresponds to the case of a lossless cylinder 

and medium, while the dashed line represents the field when the 

attenuation constants of the cylinder and the medium are 

both 2.5 db/cm.23 


I 





















Figure 


v 


Page 


2.6 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.5 illuminated by a plane wave, measured on 
the receiver line at a distance of 10 X from the cylinder. 

The solid line corresponds to the case of a lossless cylinder 

and medium, while the dashed line represents the field when the 

attenuation constants of the cylinder and the medium are 

both 2.5 db/cm.25 

2.7 Magnitude of the scattered field from a cylinder of radius 3 X and 

refractive index 1.1 (embedded in a lossless medium) illuminated 
by a point source 10 X away from the cylinder, measured on the 
receiver line at a distance of 10 X from the cylinder. The solid 
line corresponds to the case of a lossless cylinder, while the dashed 
line represents the field when the attenuation constant of the 
cylinder is 3.5 db/cm.26 

2.8 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.5 (embedded in a lossless medium) illuminated by 
a point source 10 X away from the cylinder, measured on the 
receiver line at a distance of 10 X from the cylinder. The solid 

line corresponds to the case of a lossless cylinder, while the dashed 
line represents the field when the attenuation constant of the 
cylinder is 3.5 db/cm.27 

2.9 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.1 illuminated by a point source 10 X away from 
the cylinder, measured on the receiver line at a distance of 10 X 
from the cylinder. The solid line corresponds to the case of 

a lossless cylinder and medium, while the dashed line represents 

the field when the attenuation constants of the cylinder and 

the medium are both 2.5 db/cm.28 

2.10 Magnitude of the scattered field from a cylinder of radius 3 X and 
refractive index 1.5 illuminated by a point source 10 X away from 
the cylinder, measured on the receiver line at a distance of 10 X 
from the cylinder. The solid line corresponds to the case of 
a lossless cylinder and medium, while the dashed line represents 
the field when the attenuation constants of the cylinder and 
the medium are both 2.5 db/cm. 


29 














2.11 A two component phantom illuminated with incident field Ujff).31 

2.12 This figure depicts the first- and the second-order scattered fields 

for a two component object.33 


2.13 To calculate the second order scattered fields, the scattered 
fields on the line AA’ as generated by cylinder 1 is decomposed 

into plane waves.35 

2.14 A new set of Cartesian coordinates (£,r)) is chosen such that 
the rj axis is along the vector from?, to7 2 and the origin is 

at the midpoint between 7j and 7 2 . K' b the wavenumber vector 

in the new coordinates and (7-7 0 ) b the vector position of 

the point 7 in the new coordinates.36 

2.15 A comparison between the fields on the receiver line when the 
second-order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 6 X and 1.01, respectively, a) Geometry of the simulation. 

b) Magnitudes of U,,, U 12 , U 21 and U 22 components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.42 

2.16 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 6 X and 1.03, respectively, a) Geometry of the simulation. 

b) Magnitudes of U 11 , U 12 , U 2I and U 22 components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.43 














Figure 


Page 


2.17 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 6 X and 1.05, respectively, a) Geometry of the simulation. 

b) Magnitudes of U u , U 12 , U 21 and U 22 components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.44 

2.18 A special case of the summation of different order vector scattered 
fields. Even if the magnitude of second order contribution (U 12 ) 

is large compared with that of (U M ), the magnitude of the 

summation (U n +U I2 ) may still be close to that of (U n ), 

especially when the phase angle 9 is large.45 

2.19 The magnitude of U n in Figure 2.17b is represented by the dashed 
line, while the solid line shows the magnitude of the summation of 

U 11 and 1 Jj 2 -.............45 

2.20 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 9 X and 1.01, respectively, a) Geometry of the simulation. 

b) Magnitudes of U n , U 12 , U 2 , and components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.47 

2.21 A comparison between the helds on the receiver iine when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 9 X and 1.02, respectively, a) Geometry of the simulation. 

b) Magnitudes of U u , U 12 , U 21 and U 22 components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.48 













vv 


Figure 


Page 


2.22 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 9 X and 1.03, respectively, a) Geometry of the simulation. 

b) Magnitudes of U M , U I2 , U 21 and components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.49 

2.23 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 6 X and 1.05, respectively, a) Geometry of the simulation. 

b) Magnitudes of U n , U 12 , U 2 , and components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.50 

2.24 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 6 X and 2.0, respectively, a) Geometry of the simulation. 

b) Magnitudes of Ujj, U 12 , U 2 , and components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.51 

2.25 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 0.4 X and 1.0, respectively, a) Geometry of the simulation. 

b) Magnitudes of Un, U 12 , U 21 and U 22 components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.52 












2.28 A comparison between the fields on the receiver line when the 
second order contributions are included in one case and neglected 
in the other. Diameters and refractive indices of the cylinders 
are 0.4 X and 2.0, respectively, a) Geometry of the simulation, 
b) Magnitudes of U n , U 12 , U 2 , and components of the scattered 
field on the receiver line, plotted one after another. c,d) Magnitude 
and phase, respectively, of the scattered fields are shown by the 
solid line when the second order contributions are included, and 
by the dashed line when these contributions are neglected.53 

2.27 Magnitude (a) and phase (b) of the scattered field from a cylinder 
illuminated by a plane wave, measured on the receiver line at a 
distance of 15 X from the cylinder. Solid line represents the 
magnitude of the field when it is calculated directly, while the 
dashed line corresponds to the case where the field is first calculated 
at a distance of 5 X from the cylinder and then propagated forward 
to the receiver line. Diameter, refractive index and attenuation 
constant of the cylinder are 6 X, 1.1, and 4.0 db/cm, respectively. 

The attenuation constant of the medium is 3.5 db/cm.56 

2.28 Magnitude (a) and phase (b) of the scattered field from a cylinder 
illuminated by a point source 3 X off axis and 10 X away 

from the cylinder, measured on the receiver line at a distance of 
15 X from the cylinder. Solid line represents the magnitude of the 
field when it is calculated directly, while the dashed line 
corresponds to the case where the field is first calculated at a 
distance of 5 X from the cylinder and then propagated forward to 
the receiver line. Diameter, refractive index and attenuation 
constant of the cylinder are 6 X, 1.1, and 4.0 db/cm, respectively. 

The attenuation constant of the medium is 3.5 db/cm.57 

2.29 Magnitude (a) and phase (b) of the scattered field from the 
two-cylinder object shown in Figure 2.15a measured on the receiver 
line at a distance of 20 X from the origin. Solid line represents 

the magnitude of the field when it is calculated directly, while 

the dashed line corresponds to the case where the field is first 

calculated at a distance of 10 X from the origin and then 

propagated forward to the receiver line. The refractive index of 

the cylinders is 1.1.58 












Figure 



2.30 Magnitude (a) and phase (b) of the scattered field from the 
two-cylinder object shown in Figure ?.15a measured on the receiver 
line at a distance of 20 X from the origin. Solid line represents 

the magnitude of the field when it is calculated directly, while 
the dashed line corresponds to the case where the field is first 
calculated at a distance of 10 X from the origin and then 
propagated forward to the receiver line. The refractive index of 
the cylinders is 1.5. 

2.31 Magnitude (a) and phase (b) of the scattered field from the 
two-cylinder object shown in Figure 2.21a measured on the receiver 
line at a distance of 2C X from the origin. Solid line represents 

the magnitude of the field when it is calculated directly, while 
the dashed line corresponds to the case where the field is first 
calculated at a distance of 10 X from the origin and then 
propagated forward to the receiver line. The refractive index of 
the cylinders is 1.1. 

2.32 Magnitude (a) and phase (b) of the scattered field from the 
two-cylinder object shown in Figure 2.21a measured on the receiver 
line at a distance of 20 X from the origin. Solid line represents 

the magnitude of the field when it is calculated directly, while 
the dashed line corresponds to the case where the field is first 
calculated at a distance of 10 X from the origin and then 
propagated forward to the receiver line. The refractive index of 
the cylinders is 1.5.i 

3.1 Illustration of the Fourier Diffraction Projection Theorem. 

3.2 In conventional tomography the field is measured on the receiver 

line which is perpendicular to the direction of the incident 

plane wave.' 

3.3 The Fourier transform of the scattered field on the receiver line 

gives the values of the object's Fourier transform over a 

semi-circle in the Fourier domain.‘ 

3.1 By rotation of the object over 360 degrees the object's Fourier 
transform is recovered over a set of semi-circles inside the 
circle v 2k .i 





























Figure 


3.5 The four points used for the bilinear interpolation in the frequency 


domain can be assumed to be the edges of a parallelogram if the 
rotational angular step is small and the number of points on the 
receiver line is large.79 

3.6 A schematic of the geometry of the Synthetic Aperture Diffraction 

Tomography system.80 

3.7 The direction of the wave vectors K and A relative to the 

space coordinates are depicted.82 

3.8 For a fixed orientation of the object in the transmission mode, the 

object's Fourier transform is obtained in the region shown.85 


3.9 For two rotational positions of the object, 90 degrees apart, the 
object's Fourier transform is recovered inside a circle of radius 


v2k 0 . The coverage in the region C is double.86 

3.10 The size of the object, S, the distance of the object from the 
transmitter line, d, and the distance between the arrays, D, 
are depicted in this figure.88 


3.11 When the transmitter array in figure 3.6 is utilized only in the 

reflection mode, for a fixed orientation of the object, the object's 
Fourier transform is recovered inside the region shown.90 

3.12 The coverage inside the circle of radius 2k 0 is achieved by utilizing 

the system both in the transmission and the reflection modes for 
only one rotational position of the object.91 

3.13 In the reflection mode, the object’s Fourier transform is obtained 

over the semi-circles SC( — k x k y ), one of which is the semi-circle 
centered at point A. By moving this point on the dashed 
semi-circle, the coverage shown in figure 3.11 is achieved.93 

3.1 1 The points A and B are the center of semi-circles SC'(-k x ,-k y ) 
which pass over the point D. Point D is restricted to the region 
shown.96 




































Figure 


Xll 


P.1CTP 


3.15 The curved surface shows the coven go achieved in the object's 

“Fourier-Laplace" domain when the medium is attenuatin g. The 
attenuation of the medium was taken to be 3.0 db/cm.100 

3.16 For a lossy medium, a selective part of the collected data n the 

Synthetic Aperture technique, and many rotational views of the 

object, result in the values of the object's Fourier transform 

over radial lines.102 

4.1 Cross section of a two-component object reconstructed using the 

conventional diffraction reconstruction technique. Diameter of 

the cylinders is 6 X and their refractive index is 1.02. 

a) Only the first-order scattered fields are used for generating 

the data for this reconstruction, b) Doubly scattered fields are 

included in the projection data for this reconstruction.108 

4.2 Cross section of a two-component object reconstructed using the 

conventional diffraction reconstruction technique. Diameter of 

the cylinders is 6 X and their refractive index is 1.03. 

a) Only the first-order scattered fields are used for generating 

the data for this reconstruction, b) Doubly scattered fields are 

included in the projection data for this reconstruction.109 

4.3 Cross section of a two-component object reconstructed using the 

conventional diffraction reconstruction technique. Diameter of 

the cylinders is 6 X and their refractive index is 1.05. 

a) Only the first-order scattered fields are used for generating 

the data for this reconstruction, b) Doubly scattered fields are 

included in the projection data for this reconstruction.110 

4.4 Cross section of a three-component object reconstructed using the 

conventional diffraction reconstruction technique. Diameter of 

the cylinders is 6 X and their refractive index is 1.05. 

a) Only the first-order scattered fields are used for generating 

the data for this reconstruction, b) Doubly scattered fields are 

included in the projection data for this reconstruction. Ill 











Figure 




Page 



4.5 


(a) The condition when the object components block each other is 
illustrated here, (b) With the oblique reception illustrated here 
the “direct blockage” is avoided. Direct blockage occurs when 
the direction of the propagation of the incident field is lined 
up with the line joining the two object components and the 
normal to the receiver array.113 


4.6 In the second approach several plane waves are used to illuminate 
the object. The frequency domain coverage obtained with the three 
plane waves shown here are averaged prior to inversion.114 


4.7 (a) When the incident plane wave is normal to the receiver line, 

the Fourier transform of the projection gives the values of the 
object’s Fourier transform over a semi-circle as shown here. 

(b) When the incident wave is at an oblique angle with the 
receiver line, the coverage is still on a semi-circle, although 
positioned differently.116 

4.8 In the third approach, again a number of plane waves are used to 

illuminate the object, however no independent frequency domain 
coverages are generated from each of these plane waves. For each 
rotational position of the object, the coverage generated by all 
the illuminating waves is extracted first, and then the object 
is rotated till the frequency domain is filled up.119 

4.9 Reconstructions of a two-component object using the first approach. 

The cylinders of refractive index 1.05 are 4 X in diameter, and 
positioned at a distance of 12 X from each other. The projections 
are sampled at 64 points with a sampling interval of 0.45 X. The 
receiver line is 20 X from the center of the coordinates. The angle 
between the direction of illumination and the normal to the receiver 
line is: (a) 0°; (b) 7.5°; (c) 15°; (d) 30°.120 

















Figure 


xiv 


Page 


4.10 Reconstructions of a two-component object using the second 

approach. The cylinders of refractive index 1.05 are 4 X in diameter, 
and positioned at a distance of 12 X from each other. The 
projections are sampled at 64 points over a sampling interval of 
0.45 X. The receiver line is 20 X from the center of the coordinates. 

(a) 3 incident plane waves, with the angle of illumination between 

-15° and 15°. (b) 3 incident plane waves, with the angle of 
illumination between -30° and 30°.123 

4.11 Reconstructions of a two-component object using the third 

approach. The cylinders of refractive index 1.05 are 4 X in diameter, 
and positioned at a distance of 12 X from each other. The projections 
are sampled at 64 points over a sampling interval of 0.45 X. The 
receiver line is 20 X from the center of the coordinates, (a) 3 incident 
plane waves, with the angle of illumination between -7.5° and 7.5° 
and 32 receiver positions, (b) 5 incident plane waves, with the 
angle of illumination between -15° and 15° and 32 receiver positions. 

(c) 5 incident plane waves, with the angle of illumination between 
-30° and 30° and 16 receiver positions.124 

4.12 Reconstructions of a two-component object using the first approach 

consisting of oblique reception. The cylinders of refractive index 
1.08 are 4 X in diameter, and positioned at a distance of 8 X from 
each other. The projections are sampled at 64 points over a 
sampling interval of 0.45 X. The receiver line is 20 X from the 
center of the coordinates. The angle between the direction of 
illumination and the normal to the receiver line is: (a) 0°; 

(b) 10°; (c) 20°; (d) 30°.126 

4.13 Reconstructions of a two-component object using the second 

approach. The cylinders of refractive index 1.08 are 4 X in diameter, 
and positioned at a distance of 8 X from each other. The 
projections are sampled at 64 points over a sampling interval of 
0.45 X. The receiver line is 20 X from the center of the coordinates. 

(a) 3 incident plane waves, with the angle of illumination between 
-10° and 10°. (b) 5 incident plane waves, with the angle of 
illumination between -20° and 20°. (c) 7 incident plane waves, 
with the angle of illumination between -30° and 30° 


128 





















4.14 Reconstructions of a two-component object using the third 

approach. The cylinders of refractive index 1.08 are 4 X in 
diameter, and positioned at a distance of 8 X from each other. 

The projections are sampled at 64 points over a sampling interval 
of 0.45 X. The receiver line is 20 X from the center of the 
coordinates, (a) 5 incident plane waves, with the angle of 
illumination between -20° and 20° and 16 receiver positions. 

(b) 7 incident plane waves, with the angle of illumination 

between -30° and 30° and 8 receiver positions.130 

4.15 Reconstructions of a three-component object using the first 

approach consisting of oblique reception. The cylinders of 
refractive index 1.05 are 4 X in diameter, and positioned at 
a distance of 6 X from the center of coordinates. The projections 
are sampled at 64 points over a sampling interval of 0.45 X. The 
receiver line is 20 X from the center of the coordinates. The 
angle between the direction of illumination and the normal to the 
receiver line is (a) 0°, (b) 10°, (c) 20°, (d) 30°.132 

4.16 Reconstructions of a three-component object using the second 

approach. The cylinders of refractive index 1.05 are 4 X in diameter, 
and positioned at a distance of 6 X from the center of coordinates. 

The projections are sampled at 64 points over a sampling interval 
of 0.45 X. The receiver line is 20 X from the center of the 
coordinates, (a) 3 incident plane waves, with the angle of 
illumination between -10° and 10°. (b) 5 incident plane waves, 
with the angle of illumination between -20° and 20°. (c) 7 incident 
plane waves, with the angle of illumination between -30° and 30°. ...134 

4.17 Reconstructions of a three-component object using the third 

approach. The cylinders of refractive index 1.05 are 4 X in 
diameter, and positioned at a distance of 6 X from the center 
of coordinates. The projections are sampled at 64 points over a 
sampling interval of 0.45 X. The receiver line is 20 X from the 
center of the coordinates, (a) 5 incident plane waves, with the 
angle of illumination between -20° and 20° and 16 receiver 
positions, (b) 7 incident plane waves, with the angle of illumination 
between -30° and 30° and 16 receiver positions.136 









Figure 


XVI 


Page 


4.18 A simple interpolation scheme for the Synthetic Aperture technique. 

Points A and A', B and B', and C and C lying on two different 
semi-circles SC(-k x ,-k y ) posses equal values of respectively. 

(a) small values of k y . (b) values of k y close to k 0 .138 

4.19 Points Ci and C 2 are the intersection points of the semi-circles 

0 4 and 0 2 with the circle centered at the origin going through 
point C. The neighboring points of Cj and C 2 are chosen for the 
kind of interpolation desired.140 

4.20 Magnitude of the reconstructed refractive index of a cylinder with 

diameter of 6 X and refractive index of 1.03, using the Synthetic 
Aperture technique. Transmitter and receiver arrays, positioned 
20 X apart, have 32 elements with inter-element distance of 0.45 X. 

(a) Simple nearest neighbor interpolation, (b) Modified nearest 
neighbor interpolation, (c) Simple bilinear interpolation. 

(d) Modified bilinear interpolation.141 

4.21 Magnitude of the reconstructed refractive index of a cylinder with 

diameter of 6 X and refractive index of 1.03, using the Synthetic 
Aperture technique. Transmitter and receiver arrays, positioned 
20 X apart, have 64 elements with inter-element distance of 0.45 X. 

(a) Simple nearest neighbor interpolation, (b) Modified nearest 
neighbor interpolation, (c) Simple bilinear interpolation. 

(d) Modified bilinear interpolation.143 

4.22 Magnitude of the scattered field from a lossless cylinder, embedded 

in a lossless medium, as a function of transmitter and receiver 

positions. Diameter and refractive index of the cylinder are 6 X 

and 1.03, respectively.146 

4.23 Magnitude of the Fourier transform of the scattered field from a 

lossless cylinder, embedded in a lossless medium, with respect 
to transmitter and receiver position. Diameter and refractive 
index of the cylinder are 6 X and 1.03, respectively, (a) 32 element 
transmitter and receiver arrays, (b) 128 element transmitter and 
receiver arrays.147 










Figure 


Pago 


4.24 Due to the loss of high frequency components in the projections in 

the Synthetic Aperture technique, for one rotational view of the 

object the object’s Fourier transform is obtained only in the 

region shown.148 

4.25 When the high frequency information in the projections is lost, 

larger number of views can increase the resolution and decrease 
the distortion in the reconstructed image. The coverage shown 
corresponds to three rotational views of the object 60° apart.149 

4.26 Magnitude of the reconstructed object’s Fourier transform in a 

lossless medium. The object is a lossless cylinder with diameter 
and refractive index of 6 X and 1.03, respectively, (a) 32 element 
transmitter and receiver arrays, (b) 128 element transmitter and 
receiver arrays.151 

4.27 Real part of the reconstructed object function of a cylinder with 

diameter of 6 X and refractive index of 1.03, using the Synthetic - 

Aperture technique. Transmitter and receiver arrays, with the 
inter-element distance of 0.45 X, are positioned 20 X apart, (a) 32 
element arrays, (b) 64 element arrays, (c) 128 element arrays.152 

4.28 Imaginary part of the reconstructed object function (difference 

in attenuation constant of object and medium) of a cylinder with 
diameter of 6 X and refractive index of 1.03, using the Synthetic 
Aperture technique. Transmitter and receiver arrays, with the 
inter-element distance of 0.45 X, are positioned 20 X apart, (a) 32 
element arrays, (b) 64 element arrays, (c) 128 element arrays.154 

4.29 Magnitude of the reconstructed refractive index of a 

two-component object, using the Synthetic Aperture technique. 

The components are cylinders with diameter of 6 X and refractive 
index of 1.03. Transmitter and receiver arrays positioned 30 X 
apart, have 128 elements with the inter-element distance of 0.45 X. 
a) Only the first-order scattered fields are used for generating the 
data for this reconstruction, b) Doubly scattered fields are 
included in the projection data for this reconstruction.156 

1 

I 

I 








«V«'. 
































Figure 


Page 


xvi ii 


4.30 Magnitude of the reconstructed refractive index of a 

two-component object, using the Synthetic Aperture technique. 

The components are cylinders with diameter of 6 X and refractive 
index of 1.05. Transmitter and receiver arrays positioned 30 X 
apart, have 128 elements with the inter-element distance of 0.45 X. 
a) Only the first-order scattered fields are used for generating the 
data for this reconstruction, b) Doubly scattered fields are 
included in the projection data for this reconstruction.157 

4.31 Magnitude of the reconstructed refractive index of a three-component 

object, using the Synthetic Aperture technique. The components 

are cylinders with diameter of 6 X and refractive index of 1.05. 

Transmitter and receiver arrays positioned 30 X apart, have 

128 elements with the inter-element distance of 0.45 X. a) Only 

the first-order scattered fields are used for generating the data 

for this reconstruction, b) Doubly scattered fields are included in 

the projection data for this reconstruction.158 

4.32 Energy loss of plane wave components resulting from an angular 

spectrum expansion in a lossy medium as a function of directional 
angle and attenuation constant of the medium. 7 is a measure of 
the direction of propagation of these plane waves. 7 = 0 
corresponds to the perpendicular direction to the line of expansion, 
and 7 = 1 corresponds to the parallel direction to that line.162 

4.33 Object’s Fourier domain coverage for the conventional tomography 

in transmission and reflection modes, in a lossy medium. The lower 
resolution in the conventional transmission tomography is due to 
the smaller radius of the circle of coverage as shown here.163 

4.34 When medium is lossy, larger number of views in Synthetic 

Aperture technique can increase the resolution and decrease 
the distortion in the reconstructed image. The coverage shown 
corresponds to three rotational views of object, 60° apart.165 




v. «V«*. 


% A A 









































4.35 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.01, using the conventional technique. Both the object and 
the medium are lossless. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is at a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference in 
attenuation constant of the object and the medium.166 

4.36 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.01, using the conventional technique. The attenuation 
constant is 1.0 db/cm in both the object and the medium. The 
projections are sampled at 64 points over a sampling interval of 
0.45 X. The receiver line is at a distance ol 10 X from the cylinder. 

(a) Refractive index, (b) Difference in attenuation constant of the 
object and the medium.167 

4.37 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.01, using the conventional technique. The attenuation 
constant is 3.0 db/cm in both the object and the medium. The 
projections are sampled at 64 points over a sampling interval of 
0.45 X. The receiver line is at a distance of 10 X from the cylinder. 

(a) Refractive index, (b) Difference in attenuation constant of the 
object and the medium.168 

4.38 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.03, using the conventional technique. Both the object and 
the medium are lossless. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is at a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium.169 

4.39 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.03, using the conventional technique. The attenuation 
constant is 1.0 db/cm in both the object and the medium. The 
projections are sampled at Cl points over a sampling interval of 
0.15 X. The receiver line is at a distance of 10 X from the cylinder. 

(a) Refractive index, (b) Difference in attenuation constant of 

the object and the medium.170 

























Figure 


Page 


4.40 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.03, using the conventional technique. The attenuation 
constant is 3.0 db/cm in both the object and the medium. The 
projections are sampled at 64 points over a sampling interval of 
0.45 X. The receiver line is at a distance of 10 X from the cylinder. 

(a) Refractive index, (b) Difference in attenuation constant of 

the object and the medium.171 

4.41 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1 . 0 , using the conventional technique. The attenuation 
constant of the medium and the object are 0.1 db/cm and 
0.2 db/cm, respectively. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is at a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium .172 

4.42 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1 . 0 , using the conventional technique. The attenuation : 

constant of the medium and the object are 2.0 db/cm and 
2.1 db/cm, respectively. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is at a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium.173 

4.43 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1 . 0 , using the conventional technique. The attenuation 
constant of the medium and the object are 0.1 db/cm and 
1.0 db/cm, respectively. The projections are sampled at f »4 points 
over a sampling interval of 0.45 X. The receiver line is at a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium.174 

4.44 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.0. using the conventional technique. The attenuation 
constant of the medium and the object are 2.0 db/cm and 
3.0 db/cm, respectively. The projections are sampled at 61 points 
over a sampling interval of 0.45 X. The receiver line is a' a distance 
of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium.175 











xx i 

Figure Page 


4.45 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.01, using the conventional technique. The medium is 
lossless while the attenuation constant of the object is 0.5 db/cm. 

The projections are sampled at 64 points over a sampling 

interval of 0.45 X. The receiver line is at a distance of 10 X 

from the cylinder, (a) Refractive index, (b) Difference in 

attenuation constant of the object and the medium.177 

4.46 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.01, using the conventional technique. The medium 

is lossless while the attenuation constant of the object is 

3.0 db/cm. The projections are sampled at 64 points over a 

sampling interval of 0.45 X. The receiver line is at a distance 

of 10 X from the cylinder, (a) Refractive index, (b) Difference 

in attenuation constant of the object and the medium.178 

4.47 Reconstruction of a cylinder with diameter of 6 X and refractive 

index of 1.02, using the conventional technique.'The attenuation 
constants of the medium and the object are 3.0 db/cm and 
3.5 db/cm, respectively. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is at a 
distance of 10 X from the cylinder, (a) Refractive index. 

(b) Difference in attenuation constant of the object and 

the medium.179 

4.48 The object is a lossy cylinder with diameter and refractive index 

of 6 X and 1.03, respectively. The attenuation constants of the 
medium and the object are both 3.0 db/cm. There are 128 elements 
in the transmitter and the receiver arrays, with the inter-element 
distance of 0.45 X. (a) Magnitude of the scattered field, as a 
function of the transmitter and receiver positions, (b) Magnitude 
of the Fourier transform of the scattered field, with respect to 
the transmitter and receiver positions, (c) Magnitude of the 
reconstructed object’s Fourier transform.181 


















4.49 Magnitude of the reconstructed refractive index of a cylinder 

with diameter of 6 X and refractive index of 1.03, using the 
Synthetic Aperture technique. The attenuation constant of the 
object and the medium are both 3.0 db/cm. The transmitter and 
the receiver arrays, positioned 20 X apart, have 128 elements 
with the inter-element distance of 0.45 X. The attenuation 
of the medium has been compensated for.183 

4.50 Magnitude of the reconstructed refractive index of a cylinder 

with diameter of 6 X and refractive index of 1.03, using the 
Synthetic Aperture technique. The attenuation constant of the 
object and the medium are both 3.0 db/cm. The transmitter and 
the receiver arrays, positioned 20 X apart, have 128 elements 
with the inter-element distance of 0.45 X. The attenuation of 
the medium has not been compensated for.184 

4.51 Magnitude of the reconstructed refractive index of a cylinder 

with diameter of 6 X and refractive index of 1.03, using the 
Synthetic Aperture technique. The attenuation constant of the 
object and the medium are both 3.0 db/cm. The transmitter and 
the receiver arrays, positioned 20 X apart, have 128 elements 
with the inter-element distance of 0.15 X. The attenuation of 
the medium has not been compensated for. but four rotational 
views of the object has been used.185 

4.52 Magnitude of the reconstructed attenuation constant of a cylinder 

with diameter of 6 X and refractive index of 1.0, using the 
Synthetic Aperture technique. The attenuation constants of the 
object and medium are 3.0 db/cm and 2.0 db/cm, respectively. 

The transmitter and the receiver arrays, positioned 20 X apart, 
have 128 elements with the inter-element distance of 0.45 X. 

The attenuation of the medium has been compensated for.186 








Figure 


XXI 11 


4.53 Reconstruction of a cylinder with diameter of 6 X and refractive 
index of 1.02, using the Synthetic Aperture technique. The 
attenuation constants of the object and medium are 3.5 db/cm 
and 3.0 db/cm, respectively. The transmitter and the receiver 
arrays, positioned 20 X apart, have 128 elements with the 
inter-element distance of 0.15 X. The attenuation of the 
medium has been compensated for. (a) refractive index of the 
object (b) difference in attenuation constant of the object 
and the medium. 


w . *. '• '* N *. % * - ", , N % *, _**. *• A v ,'> j,'* * -'V 
































CHAPTER 1 
INTRODUCTION 


Tomography refers to the measurement in cross-secto planes uf t:-> 
parameters from either transmission or reflection data. I i most ce_-s, i 
?is : -n« structure i- illuminated from many different direct i. n- and the image 
two-dimensional map of the chosen ti.-.-ue parameter) is rec a .■meted from I 
data •, ollected either in transmission or reflection (called project: ,ns). 

Tfmre is a fundamental diiTerence between tomographic imaging with 
rays on the one hand, and with ultrasound and microwaves <n the other hat 
X-ravs. being non-diffracting, travel in straight lines, and therefore, t 
projection data measure the line integral of some object parameter aln 
straight lines. This makes it possible to apply the Fourier Slice '! h»or< 
[Rosen felds-jj, which say % that the Fourier transform of a pr>T ti : :s e.jual 
a slice of the two-dimensional Fourier transform of the objt •••. tin- 
f' rms the foundation of the filtered-backproj-ctjon algoro m.; ■ ; with r: 

diffracting sources, such as x-rays. 

On the other hand, when either microwaves or ultrasound are used I 
tomographic imaging, the energy often does not propagate along straight 1m 
When the object inhomogeneities arc large compared to a v, iv . length, cmt 
propagation is characterized by refraction and multipath • ■< N*■"r. 

amounts of ray bending induced by refraction can be taken uo account 
combining algebraic reconstruction algorithms [AnderscnMj w;;h iigital r 
tracing and ray linking algorithms (Andersen^'ij. The multipath effects, can-' 
by the arrival at the receiving transducer of more than ou< ray prnpagati 
through different parts of the refracting object, can usually !•’ elmiinat 
through homomorphic processing of the projection data [( rawArds. 

W’mn t J,i- object inhomogeneities become comparable m -mo t<> 
wave eiigth, it is not even appropriate to talk about propagation along l.n 
and rays, and energy transmission must be discussed in terms of wa.- fr a 
and fo lds scfitt.-rcd bv the inhomogeneities. With diffracting illurnmat i u: 
e\e. t -oitition to the problem is made difficult by the complicated nature of * 
’vv itter interaction. These interactions are descri* d by ditfereip; i! 




iTnTn 



















integral equations which, in general, <lu not possess closed-form soimTwo 
types of solutions, called the first-order Born and Rytov, are avail. a! i- :.t this 
time if one can get away with the assumption of weak interactions implying 
small inhomogeneities. By either theoretical reasoning or computer simulations 
on simple cross-sectional distributions, it is usually possible to determine the 
magnitude of the inhomogeneities, which if exceeded would lead to a 
breakdown in the first-order Born or Rytov based solutions. 

Small perturbation (Born and Rytov) solutions for reconstruct'!! n were 
initiated in 1960 by Wolf [WolfSOj. He showed that with an incident 
monochromatic plane wave, the Fourier transform of the scattered data 
measured on a receiving planar aperture and the Fourier transform of the 
object, function are linearly related. Similar results for three-dimensional 
reconstruction were obtained by Dandliker and Weiss [Dandliker'/O], Carter 
used this method for reconstructing refractive index distribution by using 
holograms to record the scattered fields in the Fresnel zone of a 
semitransparent object [CarterTO], In spite of many experimental problems, his 
results were quite promising. Iwata and Nagata addressed some of the 
limitations of this approach and proposed a cylindrical receiver line vis-a-vis 
the straight line (as suggested by Wolf) for small objects [Iwata75], 

Some theoretical procedures for reconstructing weakly inhomogeneous 
objects hate been recently' derived by several researchers. A three-dimensional 
inverse scattering scheme for a spherical recording aperture enclosing the object 
was derived by Ball and Stenger [BallSO]. They considered monochromatic 
spherical wave sources for the object insonification and point receivers for the 
measurement of the scattered field. Norton and Linzer [NortouSl] have 
considered broadband insonification to derive reconstruction methods for 
planar, spherical and cylindrical geometries. 

The method suggested by Wolf has been employed by Mueller et ai. 
[Mueller?*). KavehTOaj to implement a cross-sectional reconstruction scheme for 
imaging the refractive index of objects. They have shown that a coverage of 
the Fourier transform of the object function can be achieved by measuring the 
scattered field on a receiving line for different rotational positions of the object 
'iver -%0 degrees. Mueller et al.’s method gives the values of the Fourier 
tran-form of the object on a set of arcs, while for reconstructing the object 
funeti' n through the discrete Fourier transform operation on a digital 
compiler. the values are required on a two-dimensional rectangular grid. 
J)evnney if), vaneyvjj has derived a method (suggested independently by 
Soumekh e-, -j. [Soumrkh^.'tJ 1, which he has called the Fourier back- 




















propagation technique to solve this problem. This approach i- similar in spirit 
to Fourier Back-Projection for straight ray tomography. A drawback of this 
approach is its extensive demand on computer time. 

Nahamoo and Kak [Nahamoo82] generalized the technique of diffraction 
tomography as presented in [Mueller78j to derive a new scheme which, in 
principle, permits the use of any type of insonification, as opp i to only the 
plane wave type which is experimentally difficult to generate. With this 
approach, the requirement that the object be viewed over 360 degrees is 
replaced by the need for merely two views, which ideally should be 90 degrees 
apart (however, that is not strictly necessary). For each of these two views the 
method requires independent traversal movements for both the transmitter and 
the receiver on two parallel lines in the plane of the cross-section. 

Kenue and Greenleaf [Kenue82] extended the approach of Mueller et al. 
[Mucller79] to the multi-frequency case. They have shown that by measuring 
the projections at more than one frequency one can reduce the number of views 
required and thereby higher speeds in data collection may be achieved. Note 
that the total amount of data collected remains basically the same, because for 
each rotational position of the object, each projection must now be recorded at 
different frequencies. It is possible that in the actual implementation, the 
frequency dependence of tissue attenuation may become a source of difficulty 
for this method. 

Recently Johnson et al. [Johnson83, Tracy83, Johnson84] have -uggested a 
method based on solving systems of nonlinear algebraic equations that are 
derived by applying the method of moments to a Sine basis function expansion 
of the fields and scattering potential. Although it seems that the limitation of 
this approach is quite similar to that of the Born approximation, the domain of 
applicability of this approach still remains to be proved. 

Diffraction tomography has the potential of providing us with techniques 
for the quantitative measurement of tissue parameters that take into account 
the inherently diffractive nature of wave propagation. Diffraction tomography 
methodology is not limited by the assumptions of geometric propagation, since 
concepts such as rays are not invoked, and is capable of providing us with 
distributions for virtually all the tissue parameters. Unfortunately, at this time 
the promise of this approach can only be realized for tissues that are weakly 
scattering; but it is important to bear in mind that the condition of weak 
.-catt ering is not a theoretical limitation, and work is progressing at many 
ro'e'.urch centers toward the development of computational algorithms for 
handling larger tissue inhomogeneities. 










4 


Chapter 2 is devoted to the forward scattering process. Propagation of 
the angular spectrum of a field in a homogeneous medium is presented in 
section 2.1. In section 2.2 scattering from single component objects is 
discussed. The scattering process is considered for the case of plane wave 
incidence as well as point source illumination. The role of attenuation in the 
object as well as in the medium is described and simulation results, based on 
the theoretical formulation presented, are shown. In section 2 3, a new 
procedure for calculation of scattered fields from a multi-component object is 
proposed and computer simulation results presented and discussed. The 
prospect of the new procedure in an attenuation medium is examined at the 
end of that section. Finally in section 2.4, a testing procedure for the 
calculated scattered fields is used to check the results of the previous sections. 

In chapter 3 the inverse scattering problem is considered. The formulation 
of the process is presented in section 3.1, and the two well known Bom and 
Rytov approximations discussed. The conventional Diffraction Tomography 
and related reconstruction techniques are discussed in section 3.2. The 
derivation of the Synthetic Aperture Technique in the transmission mode is 
reviewed in section 3.3, and the corresponding theory for the reflection mode as 
well as the coverage achieved in a lossy medium are derived. 

A simulation study of multiple scattering and attenuation phenomena in 
the inverse problem is presented in chapter 4. In section 4.1, artifacts caused 
by multiple scattering are studied by considering computer simulation results. 
Several approaches to reduce these artifacts are also discussed. A simulation 
study of attenuation in Diffraction Tomography is presented in section 4.2. 

Although the only source of energy we have considered in our study is 
ultrasound, the formulation of the problem is the same for microwave energy. 
The pressure field has to be replaced by the scalar potential in the case of 
microwave. However, the incident electric field has to be polarized either along 
or perpendicular to the axis of the cylindrical objects considered. Any further 
extension of the theory still remains to be explored. 










5 


[Andersen82] 


[Andersen84] 


[BallSO] 


[CarterTO] 


[Crawford82] 


[Dandliker70] 


[Devaney82] 


[Iwata75] 


[Johnson83] 


List of References 


A. H. .Andersen and A. C. Kak, “Digital ray tracing in two- 
dimensional refractive fields,” J. Acoust. Soc. Amer., Vol. 72, 
pp. 1593-1606, 1982. 

A. H. Andersen and A. C. Kak, “Simultaneous algebraic recon¬ 
struction technique (SART): A superior implementation of the 
.ART algorithm,” Ultrasonic Imaging, Vol. 6, pp. 81-94, Jan. 
1984. 

J. J. Ball and F. Stenger, “Explicit inversion of the Helmholtz 
equation for ultrasound insonification and spherical detection,” 
Acoustical Holography, Vol. 9, A. Metherell, Ed., 1980. 


W. H. Carter, “Computational reconstruction of scattering ob¬ 
jects from holograms,” Journal of the Optical Society of America, 
Vol. 60, No. 3, pp. 306-314, March 1970. 

C. R. Crawford and A. C. Kak, “Multipath artifact corrections 
in ultrasonic transmission tomography,” Ultrasonic Imaging, 
Vol. 4, pp. 234-266, 1982. 

R. Dandliker and K. Weiss, “Reconstruction of the three- 
dimensional refractive index from scattered waves,” Optics 
Communications, Vol. 1, No. 7, Feb. 1970. 

A. J. Devaney, “A filtered backpropagation algorithm for 
diffraction tomography,” Ultrasonic Imaging, Vol. 4, pp. 336- 
350, 1982. 


K. Iwata and R. Nagata, “Calculation of refractive index distri¬ 
bution from interferograms using Born and Rytov’s approxima¬ 
tion,” Japan J. Appl. Phys., Vol. 14, Supp. 14-1, pp. 383, 1975. 

S. A. Johnson and M. L. Tracy, “Inverse scattering solutions by 
a Sine basis, multiple source, moment method - Part I: Theory," 
Ultrasonic Imaging, Vol. 5, pp. 361-375, 1983. 























6 


[Johnson84j S. A. Johnson, Y. Zhou, M. K. Tracy, M. J. Berggren, and F. 

Stenger, “Inverse scattering solutions by a Sine basis, multiple 
source, moment method - Part HI: Fast algorithms," Ultrasonic 
Imaging, Vol. 6, pp. 103-116, 1984. 

[Kaveh79a] M. Kaveh, R. K. Mueller, R. Rylander, T. R. Coulter and M. 

Soumekh, “Experimental results in ultrasonic diffraction tomog¬ 
raphy,” Acoustical Imaging, Vol. 9, pp. 433-450, 1979. 

[Kenue82] S. K. Kenue and J. F. Greenleaf, “Limited angie multifrequency 
diffraction tomography,” IEEE Transactions on Sonics and Ul¬ 
trasonics, Vol. SU-29, No. 6, pp. 213-217, July 1982. 

[Mueller78] R. K. Mueller, M. Kaveh and R. D. Iverson, “A new approach 
to acoustical tomography using diffraction techniques,” Acousti¬ 
cal Holography, Voi. 8, A. Metherell, Ed. , New York: Plenum 
Press, 1978. 

[Mueller79] R. K. Mueller, M. Kaveh and G. Wade, “Reconstruction tomog¬ 
raphy and applications to ultrasonics,” Froc. IEEE, Vol. 67, pp. 
567-587, 1979. 

[Nahamoo82] D. Nahamoo and A. C. Kak, “Ultrasonic diffraction imaging,” 
Purdue University, School of Electrical Engineering, Technical 
Report, TR-EE-82-20, 1982. 

[Norton81] S. J. Norton and M. Linzer, “Ultrasonic reflectivity imaging in 
three dimensions: exact inverse scattering solutions for plane, 
cylindrical and spherical apertures,” IEEE Transactions on 
Biomedical Engineering, Vol. BME-28, No. 2, pp. 202-220, 1981. 

[Rosenfeld82] A. Rosenfeld and A. C. Kak, Digital Picture Processing, 2nd ed., 
Vol. 1, Chap. 8, New York: Academic, 1982. 

[Soumokh83] M. Soumekh, M. Kaveh, and R. K. Mueller, “Algorithms and 
experimental results in acoustical tomography Using Rytov's ap¬ 
proximation," Proceeding of International Conference on Acous¬ 
tics, Speech and Signal Processing, Vol. 1, pp. 135.138, 1983. 

[TracyS3] M. K. Tracy and S. A. Johnson, “Inverse scattering solutions by 
a Sine basis, multiple source, moment method - Part II: Numeri¬ 
cal evaluations,” Ultrasonic Imaging, Vol. 5, pp. 376-392, 1983. 











[Wolf69] E. Wolf “Three dimensional structure determination of semi¬ 

transparent objects from holographic data,” Opt. Comm., Vol. 1, 
pp. 153-156, 1969. 




















CHAPTER 2 
FORWARD PROBLEM 


When a wave encounters an obstacle, some part of the wave energy is 
deflected from its original course. It is usual to define the difference between 
the actual wave and the undisturbed wave, which would be present if the 
obstacle was not there, as the scattered wave. When a plane wave, for 
instance, strikes a body in its path, in addition to the undisturbed plane wave 
there is a scattered wave, spreading out from the obstacle in all directions, 
distorting and interfering with the plane wave. If the obstacle is very large 
compared with the wavelength, half of this scattered wave spreads out more or 
less uniformly in all directions from the scatterer, and the other half is 
concentrated behind the obstacle in such a manner as to interfere destructively 
with the unchanged plane wave behind the obstacle, creating a sharp-edged 
shadow there. This is the case of geometric optics; in this case the half of the 
scattered wave spreading out uniformly is called the reflected wave, and the 
half responsible for the shadow is called the interfering wave [Morse68], If the 
obstacle is very small compared with the wavelength, the scattered wave is 
propagated out in all directions, and there exists no sharp-edged shadow. In 
the intermediate cases, where the obstacle is about the same size as the 
wavelength, a variety of curious interference phenomena can occur. 

In our studies, we have treated the forward scattering process, as weii as 
the inverse process, as a two-dimensional problem. In other words, the objects 
under consideration are assumed to be invariant in one dimension. Although 
most of the Diffraction Tomography formulations and reconstruction techniques 
are or can be generalized to the three-dimensional case, we believe that most of 
the current effort should be directed toward resolving the many important 
limitations still facing the two-dimensional problem. 

If the forward problem can be represented as a boundary value problem, it 
is possible to find an exact solution which consists of writing down a coordinate 
system matching the boundaries of the system, and then using the 


.VA’AViv ’ ' ' 



































9 


eigenfunction expansions for solving the partial differential equations. 
Unfortunately, only the circular inhomogeneities can be solved for in this 
manner. In the past, this procedure has been the only tool for the generation 
of simulated data for diffraction tomography [Iwata75, Kaveh79b, Azimi83], 
But note the limitation: the method is suitable only for single component 
objects that are circular in shape. 

However, if the object inhomogepeities are very small such that the well 
known Born or Rytov approximations could be justified, the scattered held can 
be calculated by modeling the scattering process as a linear system [Slanev84]. 
Use of iterative techniques, such as higher order Born or Rytov techniques, can 
increase the domain of applicability of this approach to some extent. 

Since the Angular Spectrum expansion has been extensively utilized in our 
study, either in the calculation of scattered fields or in the development of 
reconstruction techniques, section 2.1 is devoted to the propagation of the 
angular spectrum of a field in either lossy or lossless medium. Section 2.2 is 
devoted to the solution of the scattered field from a single object illuminated 
either by a plane wave or a point source. Related computer simulation results 
are also presented in that section. 

When an object possesses more than one component, due to the 
interaction between the components the totai scattered fields cannot be 
considered as a superposition of the fields due to each of the components 
individually, as generated by the above method. In section 2.3, we present a 
new way of including the interactions between the components by expanding 
the scattered fields from each component as a superposition of plane waves and 
then analyzing the interaction of each one of these plane waves with the other 
object components, and finally synthesizing the total field. The effect of 
component interaction on the scattered fields is discussed in this section by- 
considering the numerous computer simulation results presented. Finally, at 
the end of this section, the prospect of this new technique in an attenuating 
medium is examined. In section 2.1 we have considered a testing procedure for 
checking the results of our computer programs. 

Past efforts on the forward problem as applied to Diffraction Tomography 
have ignored the attenuation present in the medium or the object. Although 
ultrasound undergoes only a slight attenuation., in microwave imaging 
attenuation in the medium and the object can be a limiting factor. : this 
chapter, v- 1 • have also talked about the extension of the algorithms for the 
calculation of scattrrv' fields when both the object and the m bum are f.ssv. 






sVj,-- as,'- 







s A 


10 


2.1 Propagation of Angular Spectrum 

Formulating the propagation of fields in the Fourier domain is a useful 
concept. This is especially true from a computational point of view, since FFT 
is computationally very efficient. In a lossless medium the field and its plane 
wave spectrum are related through a Fourier transform, thus each plane wave 
can be propagated to a new receiver line by the addition of a phase shift. On 
the other hand, in a lossy material the Fourier transform does not represent the 
field in a form that is physically realistic and we must follow a pure 
mathematical representation of the concept. 

If the complex field on a line is Fourier transformed, one can identify each 
Fourier component as a single unique plane wave. The field at any other point 
can then be calculated by adding a phase shift to each plane wave and then 
summing their contribution at that point. In an attenuating media, this 
formulation changes slightly because the phase shift includes a complex 
component to represent the attenuation. 

In a lossy medium, the wavenumber assumes a complex value, the 
imaginary part of which represents the effect of attenuation in the medium. 
We will use the notation k Q for the wavenumber in a lossless medium and k c0 
when the medium is attenuating. k co is defined as: 

k co = K + J f o (2-1) 

If the field energy is attenuated in the medium by <; db/cm (the unit mostly 
used in microwave imaging), ( 0 would be linearly related to <; through the 
following relation, 

_ In 10 

£„ - - c in units of cm . (2.2) 

° orv J ' ' 













11 


u(x 2l y) = / A(k y ,x 2 ) e J k » y dv 

27r -CO 


(2.5) 


Also it must satisfy the following Helmholtz equation 
V 2 u + k co 2 u = 0 


( 2 . 6 ) 


Substituting the expression in Eq. (2.5) for u in Eq. (2.6), we get 


/ 


d 2 A(k y ,x) 


dx‘ 


:=X2 - k y 2 A(k r x 2 ) 4- k c0 2 A(k y ,x 2 ) 


k ’ y dv = 0 . 


(2.7) 


For the above equation to be satisfied, the term inside the integral has to be 
equal to zero. 


d 2 A(k y ,x) 

dx 2 


:=x,-( k co 2 “ k v 2 ) -Mk X 2 ) 


( 2 . 8 ) 


A solution for A(k y ,x 2 ) is given by 

A(k y ,x 2 ) = A(k r x,) e j (2.9) 

This equation shows the exact relationship between the angular spectra of the 
fields measured on two lines separated by the distance (x 2 ~x,). This is exactly 
the same expression for the propagation of the angular spectrum in a lossless 
medium [Gcodman68], except that k c0 is complex in this case and thus each 
Fourier component sees both a phase and an amplitude change. Substituting 
Eq. (2.9) in Eq. (2.5), we have 

u(x„.y) = — 7 A(k y ,x.) e J + k, y | dk ( 2 . 10 ) 

* 2tt _ } 


Now recall that the equation for a unit-amplitude plane wave propagating 
with direction cosines (a.J) is simply 


i7l<>« + d y ) 
B(x,y) = o x 


( 2 . 11 ) 


where X is the wavelength. Thus, when k ro is real, the complex exponential 
function inside the integral in Eq. (2.10) may be regarded as a plane wave 
propagating with direction cosines (a.i). where 


o 




(2.12a) 


and 








































12 




(2.12b) 


When k co takes complex values, a, as given in Eq. (2.12a), also becomes 
complex. Substituting a complex value of a in Eq. (2.11). causes each of the 
plane wave components to be attenuated in the x direction. Thus, while the 
direction of propagation is determined by the direction cosines 
(real part of a , J ), the direction of attenuation is fixed in the x direction. 
Note that these non-uniform plane waves do satisfy the homogeneous 
Helmholtz equation ( Eq. (2.6) ). 

After the above reasoning, the following question may arise. Why not 
consider a different plane wave expansion which is physically more realistic, 
namely a summation of plane waves for which the directions of propagation 
and attenuation are the same? For such an expansion k y in Eq. (2.3) has to be 
replaced by k co sin0, where 0 is the angle between the direction of propagation 
of the plane wave and the x axis. If the expansion exists, the field on the line 
x =x, could be expressed as (change of variable k v by 0 ) 

<X> 

u(x,.y) = — / A(^,x 1 )e jkoSin(< ’ )y e' ( » sin,#ly k co cos(0) d0 . (2.1-3) 

2 7T 


The last exponential term inside the above integral will blow up either for large 
negative or positive values of y depending on the sign of 6. This means that 
the field on the line x=x, can not be represented as a summation of 
attenuating plane waves. 


jc.v. 




2.2 Scattering from a Single Object 

In general, it is not possible to calculate a closed form expression for the 
scattered fields from arbitrarily shaped objects, and the fields must be 
calculated by means of numerical methods such as finite elements. Usually, to 
find a closed form solution to the problem, the scattering process is represented 
as a boundary value problem. The field is expressed in a coordinate system 
matching the boundaries of the pioblern. and then the partial differential 
equation is solved by applying the eigenfunction expansions. The solution to 
the spattering problem is much simpler and easier to calculate numerically 
when the object is considered to be of cylindrical shape rather than of a more 
complicated shape. 

There are two types of illumination used in Diffraction Tomograph) 
techniques: a plane wave or a point source. 


-V 

& 






A /. 


*- «r „• 




i 








13 



N.'y.* 



2.2.1 Plane Wave Illumination 

Assume a plane wave propagating in the direction x is incident on a 
cylinder with the density p e and the compressibility K e centered at the origin, 
embedded in a medium with a density of p and a compressibility of /c (see Fig. 
2.1). The complex envelope of the incident field, defined as the field present 
when the cylinder is removed from the medium, is expressed as 

Uj(x.y) = U 0 e j k °° x (2.14a) 

or in polar coordinates as 

u,(T) = U„ e j k ~ r cos( * ) (2.14b) 

where k f0 is the wavenumber and (r,0) are the coordinates of the point (x,y) 
defined in cylindrical coordinates. The cylindrical boundaries of the object, 
consisting of a single cylinder, suggest the eigenfunctions to be of cylindrical 
type. Therefore, the incident plane wave should be expressed in terms of these 
cylindrical eigenfunctions as [Morse68, \Veeks64] 

'ii(T) = U 0 V £ n j" cos(no) J n (k f0 r) , (2.15) 

n =0 

where £ n is 1 for n=0 and 2 otherwise and J n (.) is the n th order Bessel function. 
When the cylinder with its axis at r=0 is present, the field is no longer given 
by either of Eqs. (2.14a) or (2.14b). Instead there is present, in addition to the 
incident plane wave, a scattered outgoing wave. In general the form of this 
scattered wave is given by the series [Morse68] 



" 5 m = 


v A n cos(no) J n (k c r) 

n =0 


V B n cos(nd>) j .J n (k ro r) + j N n (k P0 r) 

n =0 


r < a 


r > a 


(2.16) 


where k,, is the wavenumber inside the object and N n (.) is the n th order 
Neumann function. The combination I+jN. used for r > a, is to ensure that 
all of the scattered energy travels away from the cylinder, and J, used for 
r < a. is to ensure that the scattered energy travels toward the interior of the 
cylinder. The boundary conditions to be imposed on the problem are the 
continuity of the pro-sure and the axial velocity across the object-medium 
interface for the case of ultrasound energy. 


• ' 

P 





u »# ■ ie 1 1 » ■ w p 1 p 1 p 1 p. 1 pjpj p , p. 1 wr w rr^ w v *r wi 


■ * ■ v."».■ «.**■»." •.■ «■ ■*■ 1 ■•' ^ ''i-,'” -ss 



14 



k’- 


i* 

,* •/ 


£ 


Ki 

ky! 

i 



Figure 2 1 

The scattered field from a cylinder of radius a. illuminated by 
wave, is measured at the point [r,<j>). 


a plane 







r» 


In a visco-elastic rr.edium, the complex envelope of the radial velocitv, V 
and the pressure, P, are related by the following relationship 
1 <9P 

J u V r =-— , (2.1/» 

p dr 

where is the angular frequency of the incident field. We take u s , defined as 
the scattered field in Eq. (2.16), to be the scattered pressure field. By imposing 
the boundary conditions (continuity of V r and P) at the surface of the cylinder, 
one can solve for the coefficients A n and B n , defined in Eq. (2.16). Since we are 
only interested in the scattered fields outside the object, we only bring the part 
of solution for r > a, given by [Morse68j 


u s01 - S 4 j n D n cos(n?) H n l»(kj7 


(2.18) 


where 


Dn = 


Jn(^e a ) J'n(^co a ) a ^ n(^e a ) ^nt^co 3 -) 

aJ' n (k e a) H n “'(k eo a) - J n (k e a) H„* l )'(k co a) 


( 2 . 10 ) 


o- = — , (2.20) 

K p e 

and £ n is 1 for n=0 and 2 otherwise, J n (.) is an n th order Bessel function, H^ 1 '(.) 
is a Hankel function of the n lh order, k co and k e are, respectively, the 
wavenumbers in the medium and the object, a is the radius of the cylinder, 0 is 
the angle between the direction of the incident plane wave and the line joining 
the center of the cylinder and the observation point, and, finally, p and k are 
the density and the compressibility constants of the medium. 

When the object is lossy, k e assumes a complex value and only the 
coefficients D n in Eq. (2.18) will change. However, when the medium is 
attenuating (complex k co ), the Hankel function in Eq. (2.18) is affected as well. 
Because of the attenuation the Hankel function dies out exponentially as a 
function of |t| . This is clear by considering the asymptotic form of this 
function, 

i j ( z - n - ~ 1 ') 

'(z) - -e - (2.21) 

z 

This behavior agrees with <>nr physical intuition of the problem, since th“ 
energy from a source should decay rapidly as a function of distance. 



















' v K V W V" \ »V •YrJT'J V"/ ^ 


’ V 1 V wK."V^T’ 


2.2.2 Point Source Illumination 

In some reconstruction techniques, instead of illuminating the object with 
a plane wave a small sized transducer is used as the transmitter. We will talk 
about these reconstruction techniques in later chapters. A transducer with a 
very small sized aperture can be approximated by a point source. 

Assume a point source located at 7 0 is illuminating a cylinder with the 
density p e and the compressibility K e embedded in a medium with the density p 
and the compressibility k (see Fig. 2.2). The complex envelope of the incident 
field at 7 is expressed by 

u,(7) = U o H o 01(k co |7-7 o | ) . (2.22) 

where k co is the wavenumber. The cylindrical boundaries of the object, 
consisting of a single cylinder, suggest the eigenfunctions to be of cylindrical 
type. Using the expansion of the Hank'd function in cylindrical coordinates 
[Morse53j, we express the incident wave in the following form 


Uj(7) - U 0 V) ‘ n cos n(0-<p o ) 

n = 0 1 


J n( k co r ) H n (1) ( k co r o) r < T 0 


J n (k co r 0 ) H n 0)(k co r) 


r > r 0- 


(2.23) 


When the cylinder with its axis at the origin is present the total field would be 
the incident w'ave plus the scattered wave. In general the form of this 
scattered wave is given by the series in Eq. (2.16). Again, the boundary 
conditions imposed on the problem are the continuity of the pressure and the 
axial velocity across the object-medium interface. 

By imposing the boundary conditions at the surface of the cylinder, we 
solve for the coefficients A n and B n . Since we are only interested in the 
scattered field outside the object, we only bring the result for r > a. 


>i 5 m - 1] J n D n cos n(0-d> o )J H n ,n (k, 0 | r| ) 

n =0 


where 


D„ = 


H n <M r ol 1 MM) I'n( k r„a) - J' n (M) I n (k ro a) 
«J' n (M) H n (,, (ka) - J n (M) H n "i'(ka) 


where the parameters are the same as those defined in section 2.2.1. 
Considering Eqs. (2. IS). (2.19), (2.21) and (2.23), we notice that the coefficients 
I) n are different by a factor of II n (k.„|7„| ) and also the argument of the cosine 
fiinct ion m I ■ | (2 21) includes the posit jonal angle of the point source at 7 (1 . 





















































18 


If the medium is attenuating, the same effects as we discussed in section 
2.2.1 would be observed. As the point source at 7 0 is moved away from the 
cylinder the term H n (k co |7 0 | ) in Eq. (2.25) will decrease and the coefficients D n 
would get smaller and so would the magnitude of the scattered field, as is 
expected. 

2.2.3 Numerical Calculation of Scattered Fields 

We have based our computer simulation programs on the expressions 
derived in Eqs. (2.18) and (2.24). The field is calculated on a straight receiver 
line on the opposite side of the object relative to the direction of illumination. 
All Diffraction Tomographic techniques in transmission mode require such a 
geometry for the receiver positions. Since the amount of data required for 
reconstructing the cross section of an object is relatively large, it is beneficial to 
generate efficient code to speed up the computation. We have used the array 
processor FPS AP120B to carry out the vectorized operations. Most of the 
computational effort in our problem is spent on the calculation of the Bessel 
functions, J n , and the Neumann functions, N n . A Neumann function is the 
imaginary part of a Hankel function of the first kind, H n ^. 

When both the object and the medium are lossless, we calculate the 
functions J 0 , J 1( N 0 and Nj by means of non-vectorized codes available in the 
UNIX mathematical library. Higher order terms can be calculated through the 
following recursive formulas [Bowman58] 

J n + I (x) = J„(x) , (2.26) 

and 

N„ + ,(x) = N„(x)-N„.,(x) . (2.27) 

By forward recursion, with correct values of J 0 (x) and J,(x). one fails to 
produce satisfactory values for J n (x) much beyond n = x. This is due to the 
large factor multiplied by J n (x) in Eq. (2.26), which causes any initial 
uncertainties and rounding errors to propagate in multiples of that factor. 
Knowing that J n (x) goes to zero as n goes to infinity, one can reformulate the 
problem. Starting with J n + m (x) = 0 and the arbitrary J n + m _,(x) = 1, for some 
sufficiently large m, the solution for .J n (x) is sought by backward recurrence 
[Fox68], With a "check" value at some n, say n=0, the '‘trial" solution can be 
scaled to satisfy this condition. It should be noted that recursive calculation of 
the Neumann function does not suffer from this complication. VV e have 
















developed vectorized programs for these recursive calculations (forward 
recursion as well as backward recursion), and in fact, calculation of all the 
required higher order terms takes almost no time in comparison with that of 
zeroth and first order terms which are calculated by the non-vectorized code. 

The derivative of the Bessel functions and the Neumann functions are 
calculated through the following relations [Bowman58j: 

J'„W = {(•>„-.(*) -•>„*.(*)] • (2 28) 

and 

N'.M = 7 ( N n _,(x) - N„ + ,(x) ) . (2 29) 

When the object is lossy (k e is complex), by considering equations (2.18) 
and (2.24), one realizes that only the coefficients D n will be different from the 
case of a lossless object. This means that for all possible positions of the 
receiver, the coefficients D n are fixed and need to be calculated only once for 
each cylinder. Therefore the computation time for the calculation of scattered 
fields from lossy objects would be on the same order as that of the lossless 
objects. According to Eqs. (2.19) and (2.25), calculation of the coefficients D n , 
requires a routine providing Bessel functions of complex arguments, so we made 
use of the “mmbzjn” routine in the well known “IMSL” library. 

When both the object and the medium are lossy, the computational 
difficulties increase severely. This time k co , the wavenumber in the medium 
becomes a complex number. This affects the calculation of coefficients D n , 
because not only Bessel functions with complex arguments are required, as in 
the case of lossy object, but also Neumann functions with complex arguments 
are needed as well. Unfortunately we were not able to find a routine that 
calculates the Bessel functions of the second kind in any of the software 
libraries at Purdue. However, after a long search for computational procedures 
for calculation of Neumann functions with complex arguments, we found 
several useful references [Goldstein59, Marion61, StegunoT]. Finally, we were 
led to the following procedure for obtaining Bessel functions of the second kind 
with complex arguments. 

To compute N 0 (z), the zeroth order Neumann functions with complex 
arguments, we use the following series 







20 


9 7 00 1 

N 0 (z) = f [ C + log(f) J J 0 (z) + 2 V (~1 ) n + * ± J 2n (z) . (2.30) 

l 2 n=l n 


In this equation the constant C is the Euler's constant. N,(z) can be found 
using the relation 

J,(z) N 0 (z) - J 0 (z) N,(z) = ~ . (2.31) 

The functions N n (z), n=2,3,.. may now be obtained from the usual recurrence 
relation. 


N p + i(z) = 2 p z N p (z) - N p _,(z) . (2.32) 

Although the calculation of N 0 is very expensive, higher order terms are simply 
calculated from Eqs. (2.31) and (2.32). Thus, in principle, the calculation of 
scattered fields for attenuating objects and media is not much more expensive 
than that for the case of lossless material, when non-vectorized code is used in 
both cases. 


As mentioned, in the case of a lossless medium, our computer program is 
developed so that the zeroth and the first order Hankel functions are calculated 
using non-vectorized routines and then fed to the FPS-AP120B array processor 
for vectorized calculation of the higher order terms. An attempt to implement 
a similar program for Bessel functions with complex arguments failed due to 
the complex nature of the algorithm. With more work, the calculation of 
Neumann functions could be vectorized but since it can only speed up the 
calculation by a factor of two (the calculation of J n (z) could not be vectorized) 
we do not think it is justifiable. 

The computer program developed was used to generate the scattered fields 
from a cylinder with diameter of 6 wavelengths illuminated by a 5 MHz plane 
wave. In Figs. 2.3 and 2.4, we have shown the magnitude of the scattered field 
from a lossless cylinder with the solid line and that for a lossy cylinder with the 
dashed line. The object's refractive index is of 1.1 and 1.5, respectively. In 
both cases the cylinder is embedded in a lossless material and the receiver line 
is 10 wavelengths from the center of the cylinder. Considering these two 
figures, one observes that while the field magnitude decreases in the presence of 
the attenuation in the object, the overall structure of the scattered field does 
not change appreciably. 

In Fig. 2.5, we have shown the efToct of an attenuating medium on the 
scattered fields from an object. The solid line in Fig. 2.5 represents the 
magnitude of the scattered field from a lossless cylinder with diameter of G 





■ v v » r 1 -.. .Tn ' - ' V 1 . T* T* 


t t * u.*v» 


nj_*-prM_* t*t« 


rr r ^ yr ] \ 


21 



Figure 2.3 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.1 (embedded in a lossless medium) illuminated by a plane 
wave, measured on the receiver line at a distance of 10 X from the 
cylinder. The solid line corresponds to the case of a lossless cylinder, while 
the dashed line represents the field when the attenuation constant of the 
cylinder is 3.5 db/cm. 


• V 

0-J 

\y 

*4 „ - 





position on the receiver line < * X > 


Figure 2.4 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.5 (embedded in a lossless medium) illuminated by a plane 
wave, measured on the receiver line at a distance of 10 X from the 
cylinder. The solid line corresponds to the case of a lossless cylinder, while 
the dashed line represents the field when the attenuation constant of the 
cylinder is 3.5 db/cm. 






V,- v v v' 




>\ •• 



























Figure 2.5 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.1 illuminated by a plane wave, measured on the receiver line 
at a distance of 10 X from the cylinder. The solid line corresponds to the 
case of a lossless cylinder and medium, while the dashed line represents 
the field when the attenuation constants of the cylinder and the medium 
are both 2.5 db/cm. 





























wavelengths and refractive index ef 1.1 embedded in a lossless medium while 
the dashed line corresponds to the case where the attenuation constant of the 
medium and the cylinder are both 2.5 db/cm. In each case an attenuating 
plane wave was incident on the cylinder and the field was measured aiong the 
receiver line, 10 wavelengths away from the cylinder. Fig. 2.6 reflects the same 
information for a larger refractive index of 1.5. These figures show that not 
oulj the attenuation causes a iarge loss of energy but it also causes the field on 
the receiver line to be multiplied by a smooth, bell-shaped window. This 
window is due to the direct relationship between the attenuation of the field 
and the range of propagation. 

Figs 2.7-2.10 show the same information as in Figs 2.3-2.6, while the 
incident plane wave is replaced with a point source 10 wavelengths away from 
the cylinder. 

2.3 Scattering from a Multi-component Object 

One can, in principle, obtain the exact solution to the wave equation for a 
multi-component object provided one is able to solve the boundary value 
problem for the entire object. In practice it is not possible to do so even for 
two- or three-component objects, and one must take recourse to computational 
procedures. 

For calculating the scattered fields, a major source of difficulty with 
multi-component objects is dealing with the interaction between the various 
components. Depending on the interaction between the components, the total 
scattered field may or may not bear any resemblance to the simple sum of the 
scattered fields for each of the components, assuming the others to be absent. 
In this section we have presented a new computational procedure for 
calculating the inter-component interaction. With the computer programs 
developed to date, we now’ are able to generate the scattered fields that are not 
limited by the first order assumption. Although the results shown will only 
include the second order fields for a multi-component object, the computational 
procedure can easily be generalized for higher order scattering effects. 

In what follows, we have first briefly reviewed the basic theory of multiple 
scattering (Section 2.3.1). This is followed by t ho discussion of the 
computational procedure for the calculation of the multiple scattering effects 
(Section 2.3.2). Results on two component objects are shown in Section 2 3 3 

































































v ■ w-* trw 


•w im¬ 


position on the receiver line < t X ) 


Figure 2.8 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.5 (embedded in a lossless medium) illuminated by a point 
source 10 X away from the cylinder, measured on the receiver line at a dis¬ 
tance of 10 X from the cylinder. The solid line corresponds tc the case of 
a lossless cylinder, while the dashed line represents the field when the at¬ 
tenuation constant of the cylinder is 3.5 db/cm. 


y-y v va v-v-v- 


. *• A A .V/. A A‘. 


































position on the receiver line < t X ) 


Figure 2.9 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.1 (embedded in a lossy medium) illuminated by a point source 
10 X away from the cylinder, measured on the receiver line at a distance of 
10 X from the cylinder. The solid line corresponds to the case of a lossless 
cylinder, while the dashed line represents the field when the attenuation 
constants of the cylinder and the medium are both 2.5 db/cm. 


,".\V vV.V v\- v \. v,. . ■> .> . v \. 






















position on the receiver line ( * X > 


Figure 2.10 

Magnitude of the scattered field from a cylinder of radius 3 X and refrac¬ 
tive index 1.5 (embedded in a lossy medium) illuminated by a point source 
10 X away from the cylinder, measured on the receiver line at a distance of 
10 X from the cylinder. The solid line corresponds to the case of a lossless 
cylinder, while the dashed line represents the field when the attenuation 
constants of the cylinder and the medium are both 2.5 db/cm. 




































30 


2.3.1 Review of the Basic Formulation of Multiple Scattering 

We have based our algorithms on basic theory in which the scattered field 
is expressed as an infinite summation of different order terms [Ishimaru78j. 
First order fields are obtained by considering the interaction of the original 
incident field with each component, assuming the others to be absent. First 
order fields caused by one component when incident on the other components 
generate the second order fields, and so on. 

The b asic elements of multiple scattering theory that we will use in our 
computational modeling will now be explained in detail with the help of a two 
component object (The approach is easily generalized to more than two 
components). Fig. 2.11 shows a two component phantom being illuminated by 
a field denoted by U;(7) (In the absence of the phantom the field everywhere will 
be Ujfr}). We will use u(7 a ) to denote the actual field at a point frj as shown in 
the figure. u(T a ) is equal to 

u(TJ = u,fr a j + 0,fr a ) + p,fT a ) . (2.33) 

where dj(7 a ) and <^> 2 (7 a ) are, respectively, the scattered fields at 7 a caused by the 
phantom components 1 and 2. 

Let u inc (T;) be the incident field at the center of the i th component, and let 
°s(7 a 7j) be an operator function such that when it is applied to the fields 
incident ^n the scattering component at it generates the scattered field at 


the observation point 7 a . In terms of o 3 ’s, 0, and <j> 2 are given by 

<*i(T a ) = o s (r v r,)u inc ffi) i = L2 (2.31) 

Substituting (2.31) in (2.33). we get 

uCFJ = u.ffj + o ? ff a |7,)u inc (r,) + o,(r a |r 2 )u inc (To) (2.35) 

The field incident at the site of each scatterer may be expressed as 

u inc( r i) = u i(F|) + o # rr,|7 2 )u inc (7.) . (2.36a) 

u in ,.(T 2 ) = U|fr 2 ) + o (T>jT,)u in .(T|) (2.36b) 

Substituting (2.36a,b) in (2.35). we get 

h(T ) ~ u fT\) <> (T.J 7| )u rf,) <> (7 a j 7 2 )u;(7..) 

+ "K T, K(7, i 7 2 )u i: ,..(r 2 ) ^ ".rr a ; 7 2 ).' (7 2 j T, )u(T*,) (2.371 


If we again substitute (2 36a.b) in (2.37). ue have the following expression 








O u ^ 


Incident 

field 



component I 


Figure 2.11 

A two component phantom illuminated with incident field w (F) 






























32 


u(TJ = u,(Tj + o s fT a |r,)u i (r 1 ) + o s f7 a | FoKfTo) 

+ o s0 r a|'ri)Osf r i|f2)u 1 (T 2 ) + o s (7 a j 7 2 )o s (7 2 j 7,) u ,(7,) 

+ °s(T a J Tj)o s frj j 7 2 )o s (7 2 j T, )ii inc (T,) 

+ o 8 (?;|T 2 )o s (r 2 |T|)o 9 (Fi|T 2 )u inc (F 3 ) (2.38) 

A more compact way to express the above equation is 

U (TJ = u,(f a ) + u fir5t -orderf^a^ ^second-order(T a ) 

+ U higher-order(T a ) (2.39) 

The quantities u first _ order , u second _ order and u higher _ order represent the first order, 
second order, and higher order contributions at the observation point T a . These 
are given by 

Ufirst-orderfa) = °,(T»| 7 1)«i(T,) + O s f7 a | 7 2 )U;(7o) (2.40a) 

^second-orderC^a) = To)o s (7 2 1T, )u i fF 1 ) 

+ Osteal 7,)o g (7, I T 2 )u,(T 2 ) (2.40b) 

U higher-ordertTa) — ^sC^a | r 1 | | r 2 ( r l ) 

+ Os(T a |r 2 )o s (r 2 |T,)o s (T,|r 2 )u inc (r 2 ) (2.40c) 

In Fig. 2.12, we have shown the first and the second order scattering processes 
for a two component object. For the first order scattering, each component 
interacts independently with the incident illumination, being oblivious of the 
existence of the other. To compute the second order scattering terms, each 
component interacts with the fields sent in its direction by the other 
component, and so on. 

To compute the second term in (2.39), which is given by equation (2.40a). 
we require o.(T 3 |Tj)Uj(7j), where 7, is the location of the scatterer and 7 a the 
observation point. When the components are cylinders as in Fig. 2.12. and 
when the illumination field is a plane wave, the operator is given by the 
expression in Fq '2.1^). When the source of illumination is a point source, one 
should consider t expression in Fap (2 21) instead. 






Receivers 

line 



Receiver 

line 



Figure 2.12 

Th is figure depicts the first- and the second-order scattered fields for a two 
component object. 


















34 


W-\ 












2.3.2 A New Procedure for Computing Multiple Scattering 

In this section we will use the notation k 0 for the wavenumber in the 
medium, and assume that there is no attenuation present in the medium. The 
reason will be explained in a later section (2.3.4). We can directly use Eq. (2.18) 
or (2.24) to compute the first order scattering term in (2.40a); it is the 
computation of the second and higher order terms that poses difficulty. The 
computational procedure presented in this section addresses the calculation of 
the second order term as given by (2.40b). Generalization of this procedure for 
computing the higher order scattering effects is straightforward. 

By using (2.18) or (2.24), w'e first calculate the scattered fields as generated 
by the interaction of the incident field with the first component on the line AV 
as shown in Fig. 2.13. For this purpose, we choose a new set of Cartesian 
coordinates denoted by (£ ,rj ) such that the r) axis is along the vector from 7 t to 
7 2 . The origin of the new coordinate system is at the midpoint between 7, and 
7 2 ; the position vector for the origin is therefore given by 


7, + 7 2 



(2.41) 


We will use A(a) to denote the Fourier transform of the scattered field on this 
line. It is given by 


A(o) - / o s (£,f/=0|7,) Uj(7,) e“ JO * d£ (2.42) 

“ X) 




L-v- 










m 

» 

F 

[t%V 

* 

ft 

Ml* - ' 

¥ 

i 


i 


In terms of the Fourier transform A(rv), the field at any point in the half space 
r]>0 can be expressed as 

o s (7j7,) 11,(7,) = f A(o) e jRK ' {f 7n) do (2.43) 

2 K 

where K' — (o, \/k 0 ~- a*) and R is the rotational transformation applied to the 
rj coordinates to position it along the x— y coordinates (see Fig. 2.14). N oto 
that o is the £ component of f\ , and ^/k,,"-o _> the component. 

When the magnitude of o exceeds the wavenumber k . the component of 
hd in the r/ direction becomes imaginary. 'The corresponding wave components, 
called evanescent waves, will be attenuated as they travel in rj direction: the 
larger the value of rv. the greater the attenuation Beyond a distance of about 
10 wavelengths from the source, the evanescent waves are usually unimportant. 


/.v'/.V. 



































































/ u 


cylinder 2 


cylinder 1 


Figure 2.14 

A new set of Cartesian coordinates (£,» 7 ) is chosen such that the r;-axis is 
along the vector from 7 X to To and the origin is at the midpoint between 
T, and To. ^' * s wavenumber vector in the new coordinates and 
(7-7 0 ) is the vector position of the point 7 in the new coordinates. 


* V* " «. * m '** ^ f « J* ,, M A *^V ^ ^ - - 1 " « * " - 1 * • *f » "V" « » „ * . * * “ 



























37 


For digital implementation, we rewrite equations (2.42) and (2.43) in the 
following discrete form 


T-' 


A m = A£ V) o s (nA£, t ]~0 | T,) u,(7,) e N 
_ N 


JTT m n 


(2.44) 


?i) UifF,) - 


*-i 


2ir 

NA£ 


m = — 


\ 

+ 1 

o 


A e 1 


jRK' m . (?-?„] 


f" 


(2.45) 


where A£ is the sampling interval in the £ direction, and N is the total number 
of sampling points. The variable R' m is given by 


, £7mii /. o , 2trm , * . . 

V ~ NA^ ’ V k NA£ "' 46 

The operation in (2.44) can be implemented as a Fast Fourier Transform 
(FFT). When the argument of the square-root becomes negative, the 
corresponding wave becomes evanescent. This evanescent wave propagates in 


the £ direction with a spatial frequency of 


-, while it undergoes 


attenuation by the factor exp ~r/- 


•) 2 ~k 0 2 in the r/ direction. The 


higher the spatial frequency in the £ direction is, the greater the attenuation in 
the i / direction will be. 

The sampling interval is determined by the choice of the extent to which 
on.- wishes to include the evanescent waves in the measurement. The parameter 
7 will be a measure of this inclusion. Its value will range from 0 to 1; when it is 
equal to 0 the sampling interval will be such as to exclude the evanescent 
waves. When it is equal to 1, all evanescent waves will be included. We obtain 
A£ by using the following equality 


, N , 

1 

-•'7 1 

L 2 

[ NAC 

— K . 


e 1 k ' = 7 12.17) 

The left hand side is the attenuation suffered by the evanescent wave, 
characterized by the highest spatial frequency along the £ direction. A£ can be 
calculated from (2.47) to be 




s s 
















For 7 equal to 1 , this expression yields a value of X /2 for the sampling interval, 
a result that is intuitively justifiable. 


It is less possible to be analytic about the choice of X, the number of data 
points on the £ axis. The interval between the samples in the angular 
spectrum is determined by the total number of sampling points N. This 
sampling interval, and therefore the value of X, should depend on how fiat or 
sharply peaked the angular spectrum representation is. Ideally the best choice 
for X is that for which no undersampling occurs in the discrete representation 
of the angular spectrum. 


The second order scattered fields on the receiver line may then be 
determined by applying the formula in (2.18) and ( 2 . 2 - 1 ) to each plane wave 
component. Therefore, 




V \ --iRK'm-ffj-M 
XA£ “' N ' m 


x|expr. in eq.(2.18) or (2.21) with r replaced by Toj (2.49) 

Third order scattering terms can be calculated in the same way; this time the 
second order fields are first calculated on the line AA instead of the receiver 
line; and then decomposed into plane waves. The only harrier to the 
calculation of higher order terms is the computer time involved. 

Although the rate of convergences of the series in Kqs. (2.18) and (2.24) 
are dependent on the change in refractive index, it is much more sensitive to 
the value of real(k ro )a, the higher the value of real(k ro )a is, the slower the 
convergence. In all the computer simulation results shown in this paper, the 
number of terms retained on the right hand side in Kqs. (2.18) and (2.24) was 
100 . this number being more than adequate to render the truncation errors 
negligible in all cases of interest to us. Although we could have used fewer than 
100 terms without sacrificing accuracy for some values of real(k ro )a and other 
parameters, due to the programming requirements of the Floating Point Array 





39 


ol 


o 


Processor (AP120B) that was used to generate the results, we preferred using a 
fixed number of terms. 

If we only compute the second order contributions to the scattered fields, 
the computing effort required is proportional to 0(KW1 2 ) , where K is the 
number of terms retained on the right hand side in Eq. (2.18) or (2.24), N is the 
number of plane wave components in the expansion of the fields on line AA’, 
and M the number of object components. For a two-component object, the 
CPU time was less than 1 minute. For computing higher order scattering 
contributions, the computation time is proportional to 0(KNM q ) , where q is 
the maximum order of scattering included in the simulation. 

2.3.3 Numerical Results on the Calculation of Multiple Scattering 

Although the results shown will only include the second order fields for a 
multi-component object, the computational procedure can easily be generalized 
for higher order scattering effects. 

We will now show the results obtained with the computer implementation 
of the procedure discussed in the preceding section. In this section, all our 
results are on two component objects, and only the scattering contributions up 
to the second order have been included. When the ratio of back-scatter to 
forward-scatter in a two component object i> small, which happens when the 
components are large compared to a wavelength, the scattering contributions of 
orders higher than two become negligible. This condition is satisfied in a 
majority of the simulation results shown; while for others, computational 
expediency was our only reason for not going beyond the second order scatter. 
W e should like to add that even when the object components are large, it may 
still be necessary to include the higher order scattering terms if there are more 
than two components involved. 

The results shown in Figs. 2.15 through 2.17 and 2.20 through 2.26 are for 
different situations with regard to whether the components do or do not block 
each other, their sizes, and the deviations in their refractive index values from 
the background. Our intention is not to inundate the reader with a large 
number of plots. The plots were selected either bacause they illustrate how 
rapidly the higher order scattering contributions become important as the 
controlling parameter (refractive index deviation, diameter, etc.) was changed 
incrementally over one or two steps, or because they represented different 
situations such as the components blocking or not blocking each other, etc. 


V•W'A'VVV'VVV'- V-V'.V, *» 



























There are four plots for each case considered in Figs. 2.15 through 2.17 
and Figs. 2.20 through 2.26. Plots labeled (a) show the scattering geometry, 
such as the position of the cylinders relative to the receiver line (the vertical 
line on the right) and their sizes. The incident plane wave travels from left to 
right parallel to the horizontal axis. Note that the vertical and the horizontal 
scales have been normalized by the wavelength. 

Interpretation of plots labeled (b) requires some care. To explain what 
different segments of these plots stand for, we define l' nm as follows, 

° s fr a [T n ) u,(T n ) 

°/r a | rj o.frjTj u,(7 n ) 

where 7 a is a point on the receiver line. Therefore, U nn is the singly scattered 
wave from object n, and l : nm is the doubly scattered wave generated by the 
incident field first striking the object n, the resulting scattered field then 
hitting object m, and the subsequently scattered field reaching the receiver line. 
Plots labeled (b) show the magnitudes of l T u , l T jo , IN, and l : o 2 over the 
receiver line (the vertical line on the right in plots a), simply plotted one after 
another on the same axis. In plots labeled (c) the solid line shows the 

o o 

magnitude of u 3d = V l' nm (the singly plus doubly scattered field) 

n =1 m = l 

O 

whereas the dashed line represents the magnitude of u 3 = V U nn (the singly 

n = 1 

scattered field). In plots labeled (d), the solid line represents the phase of u. d , 
whereas the dashed line represents the phase of u s . The phase plot is shown 
only over the middle part of the receiver line in order to avoid the oscillations 
near the ends of this line which detract from the effects we wish to highlight. 
It should be noted that the subscript I is always used for the cylinder residing 
on the left or at the bottom (compared to the other one). Also the starting 
point of horizontal axes in plots labeled (b) and (c) corresponds to the top most 
point on the receiver line depicted in (a). 

The cases that we have studied can be divided into two categories. In one 
the cylinders are in line with the direction of illumination, and therefore are 
blocking each other. In the other category, the cylinders are not in the shadow 
of their partners. All the studies in diffraction imaging assume either the Born 
or the Rytov approximations as the first st op. Although the Rytov 
approximation is considered to be superior to the Born approximation to some 
extent, they both impose a severe restriction on the maximum allowable 
















u 


change of refractive index of the object. The Born approximation implies th<it 
the total field inside the object can be approximated by the incident field 
However, the doubly scattered field U 12 (the scattered field of the second 
cylinder when the incident field is assumed to be the scattered field of the other 
cylinder) will be negligible only if the total field incident on the second cylinder 
is close to the original incident field. Therefore we expect that long before the 
Born approximation breaks down for each cylinder, the effect of double 
scattering should become large enough to distort the received field. This is 
borne out by the computer simulations discussed in what follows. 

Figs. 2.15, 2.16 and 2.17 show the results for the refractive index changes 
of 1, 3 and 5 percent, respectively. The cylinder diameters were 6X for all 
these cases. In plots labeled (b), the reader may wonder about the large 
magnitude of U 12 in relation to that of l T ,,. Since the refractive index 
variations in the object are small, the reader might rightfully expect the 
magnitude of the summation U,, +1 )2 t0 be close to that of Ujj, and might 
therefore find a contradiction in the large values of U,o shown. However there 
is no contradiction. Since the summation of U n and U 12 is in the complex 
domain, even if (he magnitude of U 12 is larger than l r n , the magnitude of the 
summation of the two may indeed be very close to that of U n . This is 
illustrated in Fig. 2.18, where we have assumed a large angle, 0, between l',, 
and l’| 2 . This angle can be expected to be large when dealing with object 
components that are large compared to a wavelength, this being caused 
primarily bv the differences in the speed of propagation between the second 
component and the background. To show that indeed this is the case, we have 
chosen Fig. 2.17, which is the worst example amongst Figs. 2.15, 2.16 and 2.17 
in the sense that it is characterized by the largest | U 12 | /| U,,| , and plotted in 
Fig. 2.19 the magnitude of the fields U,, and l’ u +U, 2 on the same scale. As 
was expected, the difference between | U n +U 12 | and | U n | is small. 

In Figs. 2.15, 2.16 and 2.17, a simple comparison between the solid and 
dashed lines (single plus double scattered field versus the single scattered field) 
in plots labeled (c) and (d) reveals that a 3 percent change in refractive index is 
enough for the double scattering effect to change the total scattered field 
completely. Although the magnitude of the field does not change considerably, 

the deviation in phase is as high as -y. Usually the Born approximation is 

assumed to hold when the change in the phase of the incident field going 
through the object is less than 2rr. In figures 2.9, 2.10 and 2 11 the approximate 
phase shift for each cylinder at a point, on the axis joining the center of the 















n : I 01 



•o*:ria« or* «ci!yc* uni *mm»* or* «ct:v«i uni 

(c) (d) 


Figure 2.15 

A comparison between the fields on the receiver line when the second- 
order contributions are included in one case and neglected in the other. 
Diameters and refractive indices of the cylinders are 6 X and 1.01, respec¬ 
tively. a) Geometry of the simulation, b) Magnitudes of f M , C 12 , f.' 21 
and U 22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 























43 


-*l! J l« J «t 



^ogriori a* «ctrva iHC an «ccx/^c» tf»c 

(c) (d) 


Figure 2.16 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 6 X and 1.03, respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of £/,,, U l2 , U 2l and 
0^2 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


























\r -■ a." -j* ~j" _» w 





■ V - » ' > *• V ' V ■ '* w rf ^ V - V — k ' 


" V ” V ' V ’ V w V " i. ’ 


O O 


ii : I 

1 /\ 


t a -« 4 -n i -it • 


*o*iT;af» y* trctiuct .:•< 


V. M / 


• : .4 • 

7:1 i 

hi 


It -Iff -lit ■< X IM ‘ 1 .1* I* « 

•ujptor* a* •cct:^ .:•< 


a -I it ■ > :I i ■ « "* It * 

•c^i’tor* •ccti'Jd »:< 



Figure 2.17 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 6 X and 1.05, respective¬ 
ly. a) Geometrv of the simulation, b) Magnitudes of L’ n , b r 12 , b r 2 i and 
components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


S.'s /. . . . 


: : : i-v'v 

_ .-C- ... . s .«■„ s .W 
























u i 1 (^ r ) +u ,i ( *) 


u«<7) 





Figure 2.18 

A special case of the summation of different order vectors of scattered 
fields. Even if the magnitude of the second-order contribution (f/, 2 ) is 
large compared with that of (£/ n ), the magnitude of the summation 
((/ji + (/ 12 ) may still be close to that of (f/,,), especially when the phase 
angle 6 is large. 


0,00 f - - 

-23 i -19 2 -12 » 


.*0 3 00 


19 2 23 i 


POSITION ON RECEIVER LINE 


Figure 2.19 

The magnitude of T',, in Figure 9.17b is represented by the dashed line, 
while the solid line shows the magnitude of the summation of U n and 
V n . 


^.V.V V VV/Vy.V’,;.-;.- . V 

fc.V_VLV_ y.yJVvV. v. V.V. _L“L« , v vlV-"w2vlvl‘vI ^-1 


--JV*’j* ', 


^ -r "_n .* « *» Vi Vj " * 














46 


cylinders is respectively , —■ and Therefore the total scattered field 

gets '’distorted” long before the phase shift along the length of the objects 
reaches the limit of 2 t. In Figs. 2 . 20 , 2.21 and 2.22 the change in refractive 
index is one, two and three percent, but the objects are taken to be larger in 
size, 9X to be precise. This time a two percent change seems to be the upper 

limit. The phase shift in each cylinder is now - 7 -, — and - 77 -. 

3 3 3 

In Figs. 2.23 and 2.24 we have changed the position of the cylinders to see 
the effect when they are not in each other s shadow. The refractive index of 
the cylinders is taken to be 1.05 and 2.0, respectively. It is quite clear that 
multiple scattering has a very small effect on the result in this case, even for 
very high refractive indices of the object components. This is expected because 
the object components are large and thus forward scattering. However, even 
when the size of the cvlinders is decreased in Figs. 2.23 and 2.26, no significant 
change is observed (refractive indices were respectively, 1.05 and 2.0 in those 
figures). Since the scatterers are not directional in this case, the scattered 
energy spreads out in all directions every time a wavefront hits a component. 
The multiple scattering effects are minimal in this case because the energy jn 
the doubly scattered fields is much lower than that contained in the singly 
scattered fields, due to the two isotropic scattering interactions involved in the 
former, compared with a single such interaction in the latter. 

2.3.4 The New Procedure In an Attenuating Medium 

Although the algorithms for calculating the scattered fields from non¬ 
attenuating cylinders and media were easily extended to lossy objects and 
media, the same is not true for calculating the scattered fields from multiple 
cylinders in a lossy medium. This limitation is due to the special characteristic 
of the angular spectrum expansion technique in a lossy medium, as discussed in 
section 2 . 1 . 

In the case of an attenuating medium, it is not possible to associate the 
Fourier transform of the field with uniTorm attenuating plane waves. Tim 
Fourier components represent non-uniform plane waves which possess different 
directions of propagation and attenuation. The interaction of these non- 
uniform plane waves with cylindrical objects can not be formulated as before. 

Assume the scattered field from a cylinder is calculated on the line x =0 
and is expanded into non-uniform plane waves propagating toward a second 
cylinder. The aim is to formulate the interaction of these non-uniform plane 


V .* Y- v %“ . 

•V- V A'-V'-A'- 




















f «»"r 


"T 


w 

Jv 

L\ \n 






> „% 


**_“ *s‘ 

& 

v 

0 

£;> 

W. 

■1 1 


El 

Ei 


■:•> 
r;-\" 
i-.v," 

!■.;< 

r. -\ 


> - 
v>: 

fo-c- 



vj ".j u ■.» ".• ■• ‘j ■ 


I P'*- f>.« I 


47 



Figure 2.20 

A comparison between the Selds on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 9 X and 1.01, respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of l' u , t', 0 , Uy and 
V 22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 



. 

V. ’ Aji 














\ f H 


t » H 3 4 -It t -It $ -* •** • » • 



»oem9i ■tnciuo» Lire 

(b) 



-a i -if i -it • -* «• « m i a it a if ■ at 

•VITIQfi »t ttctluo lIK 


i »1 I j , 



■ 1T1WI Oh aCTiLC* -ire 


Figure 2.21 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 9 X and 1.02, respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of f r n , L j 2 , b ; o, and 
U 22 components of the scattered field on the receiver line, plotted one 
after another. e,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line wheD these contributions are neglected. 




■ *• a h. * *. 












r^TTTV 




Figure 2.22 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 9 X and 1.03, respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of £/,,, t/ 12 , (/ 21 and 
U 22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


• ' 

|t« r . 


-/ •.* . 

:/,v, ■ / 

JLjlL jl Tj. V,V, l 


JdA 


' , ******** */• 













linen on (ccxiuot Life 




t «• u ■ i* • a ( 


•o*ma* a* «ctu*» Life 


-it i -t i -4 m -i it - m * n « a i « m 


"osiTion w tcctrc* 


Figure 2.23 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 6 \ and 1.05, respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of t'.l- ^ t/ 21 and 
U 22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


->v\* 
















T 






Ti",i 





























(a) 



*o»mon on «ctiuc* u< 


(c) 


•otiTion an tcaiue* Lire 

(b) 



»>« * - 1 ■ — -—— 

-it t - 4 « 4 a -* i» • m ni » » »^ 

potitjon on tctXluCR line 


(d) 


Figure 2.24 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 6 X and 2.0. respective¬ 
ly. a) Geometry of the simulation, b) Magnitudes of (/,,, U { n, t/ 2 . and 
U 22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


V- V'-'y 'i'•-* y" 


is f 










































V w T w V" V» V« V w 1T» S " V * 


Vw; lry -w -j ryi v»v ^ w - w " * - wr ^_r kji - -v - *%• > -v tr* lt* 


^v/ 



52 





f./* 

t. 

r v ••. /. ■'. •\*'r- 

> v ■/vV.-.'/a'-’/"-'/.!''. 

C/- w\Jaj.Y. ilmAji j.. 



Figure 2.25 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 0.4 X and 1.05, respec¬ 
tively. a) Geometry of the simulation, b) Magnitudes of U u , (V 12 , (/„, 
and i ;2 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 










*o«inan ah «ct:uo» l. 


K A i 

i W i it 


« M t « »• • I* * ■ * 


•II* -W J IS »* « * '** 


*o*man »» name* urc 


smoh an waive* l:w 


Figure 2.26 

A comparison between the fields on the receiver line when the second ord¬ 
er contributions are included in one case and neglected in the other. Di¬ 
ameters and refractive indices of the cylinders are 0.4 X and 2.0, respec¬ 
tively. a) Geometry of the simulation, b) Magnitudes of L' u , f/ 12 , i' 2l 
and L '22 components of the scattered field on the receiver line, plotted one 
after another. c,d) Magnitude and phase, respectively, of the scattered 
fields are shown by the solid line when the second order contributions are 
included, and by the dashed line when these contributions are neglected. 


" *** k.* *■ n G/^ •ji. m * 1. TVt * 
























waves with the second cylinder, 
plane waves can be written as 

, i _ r - j ( k.x + k,v 

u,(x.y) - 1' 0 e J 

where k v has a real value. k x 
relationship 


The mathematical expression for one of these 

C2.51) 

would be complex because of the following 


K 




- K r + i k 


If the position of the point (k x ,k v ) is defined by (k,.0) in cylindrical coordinates 
Eq. (2.51) can be rewritten in the following form 


jk,rco«(0~0) j (j kx, r ) oo ?0 


uifr) = l' o e J ■•'■-"• e' ' (2.53) 

where (r.o) are the coordinates of the point (x.y) in cylindrical coordinates. 

To solve the boundary value problem, one has to expand the incident fit 
in terms of cylindrical basis functions. The result is 


(2.54) 



’sc. ft 



ffifr) = <- ; o 

v; £n j n COS n(0-o)l J n (k,r) 

X 

v; ( n C os(nd) I n (k x r) 


n = 0 v ’ 


n =0 


where I n (—j n J n ) is the modified Bessel function. 


Because of the complex nature of the expression in Eq. (2.54), it is quite 
difficult, if not 'mpossible, to impose the boundary conditions on Eqs. (2.16) 
and (2.54) and solve for the coefficients A n and B n . Therefore, we can not 
account for the interaction of the plane waves generated by the angular 
spectrum expansion with the object components. This conclusion suggests that 
the extension of our algorithm for the calculation of scattered fields from 
multi-component objects to the case of a lossy medium is not feasible. 


2.4 Test of Numerical Results 

A simple test of validity or the simulation results is to check whether the 
calculated scattered fields satisfy the wave equation. Due to the Huygen s 
principle when a field is known over a line in a two dimensional space, it is 
known over the whole space. However, to test a numerical technique, one has 
to check the validity of the result over the range of interest and prove the 
numerical stability of the technique. To test our computer programs for 
calculation of scattered fields, we made use of the formulation of propagation 
of the angular spectrum (see section 2.1). We first calculate the scattered fields 








r lf> if .V> irnti i iTm/i »< «> AUmdt 










along the line x^Xj and then propagate it to the line x =Xo, using the 
formulation in Eq. (2.9). This field is then compared to the one directly 
calculated by our program along the line x=x 2 . We choose the two ranges x, 
and x 2 to be at the two extremes of the range of interest, to make sure that the 
field generated by our program in that region satisfies the wave equation. 
Although there is a possibility that the fields calculated do satisfy the wave 
equation but do not correspond to the exact problem the algorithms are 
designed for, one should consider its very low probability . 

The solid line in Fig. 2.27(a) shows the magnitude of the scattered field as 
calculated by our program along a line 15 wavelengths from the center of a 
cylinder when it is illuminated by a plane wave. The dashed line shows the 
field which is first calculated 4 wavelengths from the cylinder and then 
propagated to a distance of 15 wavelengths from the object using the 
formulation in Eq. (2.9). The phase information is reflected in Fig. 2.27(b). 
The attenuation in the medium is assumed to be 3.5 db/cm. Notice that the 
two calculated fields match exactly and the dashed line may not even be 
recognized in this figure. 

In Fig. 2.28, we have shown the same information as in Fig. 2.27, except 
that the incident plane wave has been replaced by a point source. The point 
source was positioned on a line parallel to the receiver line 10 wavelengths 
away from the cylinder and 3 wavelengths off axis. Again the two calculated 
fields are very close. 

For testing the program that calculates the scattered field from a multi- 
component object, we chose the geometry in Figs. 2.15, where the components 
block each other, and the one in 2.23, where no blocking occurs. We chose the 
incident field to be a plane wave; however, a point source illumination could 
have been used as well. In Figs. 2.29 and 2.30 we have shown the same 
information as plotted in Figs. 2.27 and 2.28, while using the geometry of Fig. 
2.15, with the refractive indices of the components being 1.1 and 1.5, 
respectively. Figs. 2.31 and 2.32 are generated for the geometry used in Fig. 
2.23. 

Since the fields calculated by the method of propagation of the angular 
spectrum are very close to the ones calculated directly, we have been persuaded 
that our computer programs are error-free. 




position on ihf recewtr lin* < * X ) 
<b> 


Figure 2.27 

Magnitude (a) and phase (b) of the scattered field from a cylinder il¬ 
luminated by a plane wave, measured on the receiver line at a distance of 
15 X from the cylinder. Solid line represents the magnitude of the field 
when it is calculated directly, while the dashed line corresponds to the case 
where the field is first calculated at a distance of 5 X from the cylinder and 
then propagated forward to the receiver line. Diameter, refractive index 
and attenuation constant of the cylinder are 6 X, 1.1, and 4.0 db/cm, 
respectively. The attenuation constant of the medium is 3.5 db/cm. 































position on th* rectiutr line ( * \ ) 
<a> 



-i* ' -!• * -» Ji -> rr - va j M in im i< i 

position on th* r»criv*r line ( t X > 

<S) 


Figure 2.28 

Magnitude (a) and phase (b) of the scattered field from a cylinder il¬ 
luminated by a point source 3 X off axis and 10 X away from the cylinder, 
measured on the receiver line at a distance of 15 X from the cylinder. 
Solid line represents the magnitude of the field when it is calculated direct¬ 
ly, while the dashed line corresponds to the eas_ where the field is first cal¬ 
culated at a distance of 5 X from the cylinder and then propagated for¬ 
ward to the receiver line. Diameter, refractive index and attenuation con¬ 
stant of the cylinder are 6 X, 1.1, and 4.0 db/'em, respectively. The at¬ 
tenuation constant of the medium is 3.5 db/cm. 








position on ih* receiver line C t X ) 


position on the '•eceiver 1 • <. * a ) 


Figure 2 20 

Magnitude (a) and phase {b) of the scattered field from the two-cylinder 
object shown in Figure 9.15a measured on the receiver line at a distance of 
20 X from the origin. Solid line represents the magnitude of the field when 
it is calculated directly, while the dashed line corresponds to the case 
where the field is first calculated at a distance of 10 X from the origin and 
then propagated forward to the receiver line The refractive index of the 
cvlinders is 1 1 





















position on th* receiver line ( » X ) 

Cj) 


& 4*«0 


l 

•* - ***00 



-n % -J* f -7 31 -3 77 - 223 3 3* k H 10 * 

position on the receiver line < * X ) 

(b> 


Figure 2.30 

Magnitude (a) and phase (b) of the scattered Geld from the two-cylinder 
object shown in Figure 9 15a measured on the receiver line at a distance of 
20 X from the origin Solid line represents the magnitude of the Geld when 
it is calculated directly, while the dashed line corresponds to the case 
where the Geld is first calculated at a distance of 10 X from the origin and 
then propagated forward to the receiver line. The refractive index of the 
cylinders is 1.5. 






































position on th* receiver lm» ( t \ ) 

<b> 


Figure 2.32 

Magnitude (a) and phase (b) of the scattered field from the two-cylinder 
object shown in Figure 9.21a measured on the receiver line at a distance of 
20 X from the origin. Solid line represents the magnitude of the field when 
it is calculated directly, while the dashed line corresponds to the case 
where the field is first calculated at a distance of 10 X from the origin and 
then propagated forward to tbe receiver line. The refractive index of the 
cylinders is 1.5. 





























List of References 


[Azimi83] M. Azimi and A. C. Kak, "Distortion in diffraction tomography 
caused by multiple scattering," IEEE Transactions on Medical 
Imaging. Vo!. MI-2. No. 4. pp. 176-19.'), Dec. 1983. 

[Bowman58] F. Bowman, Introduction to Bessel Functions, Dover Publica¬ 
tions Inc., New York, 1958. 

[Fox68] L. Fox and D. F. Mayers, Computing Methods For Scientists and 

Engtneers, Chap. 3, Clarendon Press, Oxford, 1968. 

[Go!dstein59] M. Goldstein and R. M. Thaler, “Recurrence techniques for the 
calculation of Bessel functions,” Math. Tables and Aids to Comp. 
(Mow Mathematics of Computation), Vol. 13, pp. 102-108, 1959. 

[Goodman68] J. W. Goodman, Introduction to Fourier Optics, McGraw-IIill 
Book Company, 1968. 

[Ishimaru78] A. Ishimaru, Ware Propagation and Scattering in Random 
Media, Vol. 2, Chap. 14, Academic Press, New York, 1978. 

[Iwata75] K. Iwata and R. Nagata, “Calculation of refractive index distri¬ 
bution from interferograms using Born and Rytov’s approxima¬ 
tion,” Japan J. Appl. Phys., Vol. 14, Supp. 14-1, pp. 383, 1975. 

[Kaveh79b] M. Kaveh, R. K. Mueller, and R. D. Inverson, “Ultrasonic to¬ 
mography based on perturbation solutions of the wave equa¬ 
tion,” Computer Graphics and Image Processing, Vol. 9, pp. 
105-116, 1979. 

[Marion6l] G. C. Marion. “Bessel functions of integral order and complex 
arguments,” Communications of the ACM, Vol. 4, No. 4, page 
169, April 1961. 

[Morse53] P. M. Morse and H. Feshbach. Methods of Theoretical Physics, 
Part I, Chap 7, McGraw-Hill Book Co., 1953. 




63 


[Morse68] P. \1. Morse and K. V. Ingard, Theoretical Acoustics, Chap. 8, 
McGraw-Hill Book Co., 1968. 

[Slanev84| Malcolm Slaney, A. C. Kak, and L. E. Larsen, "Limitations of 
imaging with first-order diffraction tomography," IEEE Transac¬ 
tions on Microwave Theory and Techniques, Vol. MTT-32, No. 
8, pp 860-874, August 1984. 

[Stegun57j I. A. Stegun and M. Abramowitz, “Generation of Bessel func¬ 
tions on high speed computers," Math. Tables and Aids to Comp. 
(Now Mathematics of Computation), Vol. 11, pp. 255-257, 1957. 

[Weeks64] VV. L. Weeks, Electromagnetic Theory for Engineering Applica¬ 
tions, ch. 6, John Wiley & Sons, Inc., 1964. 






V V ■-». ji ^ j- 




-; vj *v w-rf »r-r tv ’■Af w > " vit T w T" 




64 


CHAPTER 3 
INVERSE PROBLEM 


In general inverse scattering is a much more complicated problem than the 
forward problem. Inverse scattering techniques are usually accompanied by 
severe restrictions on the different parameters of the system under 
consideration. 

Small perturbation solutions for diffraction reconstruction techniques were 
initiated in 1969 by Wolf [Wolf69]. He proved that when a monochromatic 
plane wave is incident on an object, the Fourier transform of the scattered data 
measured on a receiving planar aperture and the Fourier transform oF the 
object function are linearly related. Dandliker and Weiss obtained similar 
results for three-dimensional reconstruction [Dandliker70], Carter used this 
method for reconstructing refractive index distribution by using holograms to 
record the scattered fields in the Fresnel zone of a semi-transparent object 
[Carter70], In spite of many experimental problems, his results were quite 
promising. Iwata and Nagata addressed some of the limitations of this 
approach and proposed a cylindrical receiver line vis-a-vis the straight line (as 
suggested by Wolf) for small objects [Iwata75], 

Some theoretical inverse procedures for reconstructing weakly 
inhomogeneous objects have been recently derived by several researchers. A 
three-dimensional inverse scattering scheme for a spherical recording aperture 
enclosing the object was derived by Ball et al. [BallSO]. They considered 
monochromatic spherical wave sources for the object insonification and point 
receivers for the measurement of the scattered field. Norton and Linzer 
[Norton81] have considered broadband insonification to derive reconstruction 
methods for planar, spherical and cylindrical geometries. 

Mueller et al. employed the method suggested by Wolf to implement a 
cross-sectional reconstruction scheme for imaging the refractive index of objects 
[Mueller79j. They have shown that, a coverage of the object's Fourier 
transform function can be achieved by measuring the scattered field on a 




- •„* v Vv V V. 











receiving line for different rotational positions of the object over 360 degrees. 
This method provides the values of the object’s Fourier transform on a set of 
arcs, while for reconstructing the object function through the discrete Fourier 
transformation on a digital computer the values are required on a two- 
dimensional rectangular grid. Different interpolation techniques have been used 
for this purpose [Pan83, Kaveh79a]. Devaney [Daveaney82] has derived a 
method, which he has called the Fourier back-propagation algorithm (similar in 
spirit to Fourier Back-Projection for straight ray tomography), to avoid 
interpolation errors. A drawback of this approach is its extensive demand on 
computer time. A similar approach has been followed independently by 
Soumekh et al. [Soumekh83j. 

Nahamoo and Kak (Nahamoo82, Nahamoo84j generalized the technique of 
diffraction tomography as presented in [Mueller79] to derive a new scheme 
which, in principle, permits the use of any type of insonification, as opposed to 
only the plane wave type which is experimentally difficult to generate. With 
this approach, the requirement that the object be viewed over 360 degrees is 
replaced by the need for merely two views, which ideally should be 90 degrees 
apart (however, that is not strictly necessary). For each of these two views the 
method requires independent traversal movements for both the transmitter and 
the receiver on two parallel lines in the plane of the cross section. 

Devaney has recently taken an approach similar to that of Nahamoo and 
Kak in applications dealing with geophysical imaging [Devaney84]. He has 
considered the case where the transmitter and receiver arrays are perpendicular 
to each other and has derived the corresponding point spread functions. 

Kenue and Greenleaf [Kenue82j extended the approach of Mueller et al. 
[Mueller79] to the multi-frequency case. They have shown that by measuring 
the projections at more than one frequency, one can reduce the number of 
views required and thereby higher speeds in data collection may be achieved. 
Note that the total amount of data collected remains basically the same, 
because for each rotational position of the object, each projection must now be 
recorded at different frequencies. It is suspected that in the experimental 
implementation, the frequency dependence of tissue attenuation might become 
a source of difficulty for this method. 

Recently Johnson et al. [Johnson83, Tracy83, Johnson81) have suggested a 
method based on solving systems of nonlinear algebraic equations that are 
derived by applying the method of moments to a bine basis function expansion 
of the fields and scattering potential. Although it seems that the limitations of 
this approach is quite similar to that of the Born approximation, the domain of 


1 v.v, -v-w.v .y'■'.v.Nc/y 





66 


applicability of this approach still remains to be proved. 

In section 3.1, we have reviewed the bascis of the interaction of a wave 
with tissue, the integral representation of the scattered field and the two well 
known weak assumptions: the Born and the Rytov. The conventional 
Diffraction Reconstruction technique is discussed in section 3.2, and the 
extension of the technique to the case of an attenuating medium examined. 
We review' the mathematical derivation of the Synthetic Aperture technique in 
transmission mode [Nahamoo82] in section 3.3. The theory is extended to the 
reflection mode in that section and it is proved that if the reconstruction 
technique is utilized in both the transmission and the reflection modes, for a 
single view of the object one can recover the object’s Fourier transform inside a 
circle of radius 2k 0 (k 0 is the wavenumber). Finally in section 3.3.3 the effect 
of attenuation on this technique is considered. 

3.1 Formulation of the Problem 

The interaction of sound with a medium forces the particles of that 
medium to be displaced from their equilibrium positions. The elastic restoring 
forces combined with the inertia of the particles lead to an oscillatory motion 
of the medium. A mathematical description for the propagation of sound in a 
medium can be derived by formulating the principles of the elastic properties of 
the medium and setting up an equation of motion for the particles in the 
medium. 


3.1.1 Interaction of Wave with Tissue 


We only consider the longitudinal component of the sound wave, assuming 
the traversal components get damped because of the fluid structure of the 
medium and the object under consideration This is usually true for tissue 
because of its high fractional water content. Based on this model, the equation 
of motion is given by [Morse68], 

pm = -vpfr.t) . (3.i) 

at 

where V and P are velocity and pressure fields, respectively and p is the local 
density of the fluid. The equation of state of the fluid near the equilibrium 
state (small amplitude disturbances) is 



- pff.t) Kfr.t) 




where k is the complex compressibility of the viscous fluid. The local variation 






67 


of density is related to the flow of fluid at each point by the equation of 
continuity as given by 

MlA = - p fr t ) V ■ V(7, t) (3 3) 

dt 

By combining the three equations of motion (Eq. (3.1)), state (Eq. (3.2)), and 
continuity (Eq. (3.3)), we come up with the following wave equation. 


pm 


VPff.t) J = /c(f,t) 

dpm. 

-(- dt 

pm 


pm *(F,t) 


pm 


The hist term on the right hand side of this equation is the difference between 
the percentage change in time of p k and p multiplied by the percentage 
change in time of p. This term being small compared to the first term on the 
right hand side of Eq. (3.4) is usually neglected. For a single frequency case 
with P(T,t) = u(T) e^ 1 , we obtain 

V • — Vu(7) + k(7) u)'“ ufr) = 0 (3.5) 

p(n 

The wavenumber is defined as: 

k co “ Po K o 'J » ( 3 - 6 ) 

where p 0 and k 0 are the density and the compressibility of the fluid at 
equilibrium, respectively. When the fluid is lossy the compressibility, k 0 , 
becomes complex and so does k co . 

If we introduce the notations 


- 


ptf) ~ Po 


7 ,-m = 


K’cn - k 0 


(3.7b) 


and make use of definition of k ro in Eq. (3.6), we can rewrite Eq. (3.5) in the 
following form 













(3.8) 


V ' ( ( 1 - 7,(T) ] Vuprl ) + k „ 2 | 1 + 7„(J) ] <i(T) = 0 . 

By massaging the equation above, we get the final form of the wave equation 

(V 2 + IO u(F) = - k„ 2 [n 2 (r)-l| U (r) + -e 0 - V V (T) ■ V uf r) , (3. 

Po 


where the refractive index n(7) is defined as 


n(7) - 


(3.10) 


In most practical situations the second term on the right hand side of (3.9) 
is neglected and the following simpler linear inhomogeneous Helmholtz equation 
is used 


(V 2 + k cc 2 )u(T) = -F(7) ufr) , 
where the source function Ffr) is defined as 
F(f) = k co 2 [n 2 (T) - 1 ] 


(3.11) 


(3.12) 


In all the simulations considered in our study, the density function p(r) is 
assumed to be equal to p 0 , so that 7 (7) becomes zero and the last term on the 
right hand side of Eq. (3.9) vanishes. 


3.1.2 An Integral Solution 

The aim is to convert the differential equation (3.11) to an integral 
equation. An elegant description of this problem can be found in [Morse53], 
while in here we only bring a brief summary of the approach similar to that 
covered in (Nahamoo82|. 

The wave equation as given in Eq. (3.11) does not provide a complete 
definition of the problem unless it is accompanied by the boundary conditions 
which specify the behavior of the field on the surface that encloses the medium. 
The boundary conditions can be expressed in terms of the pressure field and 
the normal velocity at the boundary. Since the velocity is related to the 
derivative of the pressure according to Eq. (3.1), the pressure field and its 
normal derivative at the boundary are usually employed to express the 
boundary conditions. 

In practical cases, the boundary consists of partially rigid or flexible 
materials, while on some part of it a velocity distribution is generated by the 
transmitting transducer. A general representation of such boundary conditions 
is given by 




^B^b) ' °B^b) + z C*b) u Hb) ~ factive(r*B) > (3.13) 

where 7 B is a point on the boundary, 7 B (7 B ) is the unit vector normal to the 
boundary, z(7 B ) represents the reaction impedance between the medium and the 
boundary and factiveOB) is a function which is non-zero for the active part of 
the boundary. For a transmitting transducer, zCie)/jk o /) 0 c 0 represents the 
impedance of the transducer and f ac tiveOB)/j^oPo c o * s ca ^ e( f the driving velocity 
of the transducer. 

The wave equation given in Eq. (3.11) is called inhomogeneous because the 
right hand side of the equation is non-zero and represents a distributed source 
of energy. The wave equation is called homogeneous if the function F(F) is 
identically zero everywhere. Since the function f ac tivetrB) contributes to the 
production of the field almost in the same manner as does F(r), it is 
appropriate to name the boundary conditions similar to the wave equation. 
Homogeneous boundary conditions are associated with f a ctive(*B) being 
indentically zero, while a non-zero f act i ve (rB) represents inhomogeneous 
boundary conditions. 

The unique solution of an inhomogeneous wave equation with the 
associated boundary conditions is the sum of two components. One is the 
solution of the inhomogeneous equation for the homogeneous boundary 
conditions, while the other represents the solution of the homogeneous equation 
for the inhomogeneous boundary conditions. Clearly the sum of these two 
solutions will satisfy both the inhomogeneous wave equation and the 
inhomogeneous boundary conditions. 

The solution for the homogeneous wave equation under inhomogeneous 
boundary conditions is the field generated by the transmitting transducer, 
which we call Uj. The solution for the inhomogeneous wave equation under 
homogeneous boundary conditions is obtained by adding the contribution of all 
the differential elements of the source function, F. The field at point 7 caused 
by a point source at 7 Q is the well-known free space Green’s function g(7[7 0 ) 
which satisfies the following wave equation 

(V 2 + kj) g(7j 7,) = -<5(7-7 0 ) (3.14) 

The field generated by the distributed source function F is the integral of 
F(7 0 ) gfrj7 0 ) over the scattering object Therefore, the integral representation 
of the differential equation in Eq. (3.11). can be written as 










70 


u s (r) = / F(T 0 ) u(f 0 ) g(r]T 0 ) dr 0 , 

(3.15) 

object 


where u 3 is the scattered field defined by 


u s (r) = u(F) - ujfr) 

(3.16) 


Since the total field, u, does appear inside the integral in Eq. (3.15), the 
solution is of transcendental form, which means an inverse reconstruction 
technique based on this solution is not feasible. Some small perturbation 
approximation leads to an expression for the scattered field u s as an explicit 
function of the inhomogeneities. Two of the well known approximations are 
the Born and the Rytov approximations. 

3.1.3 Born Approximation 

In the first-order Bom approximation the total field, u, on the right hand 
side of Eq. (3.15) is approximated by the incident field, u ; . This is just ; dable 
only when the object inhomogeneities are small so that the magnitude of the 
scattered field inside the object is negligible compared to that of the incident 
field. Then the solution is given by the following integral 

u s ffl = / F(7 0 ) Ui(f 0 ) g(r|7 0 ) dT 0 (3.17) 

object 

Second order Born solutions can be obtained by replacing the total field, 
u(7) in Eq. (3.15) by the sum of the incident field and the first-order scattered 
fields calculated using Eq. (3.17). A repetitive application of this procedure 
results in higher order solutions. It is known that the higher order solutions 
converge to the exact solution if the object inhomogeneities are small. Although 
a weak sufficient condition for the convergence of the solution exists, the 
necessary condition is still unknown [Hochstadt73j. 

3.1.4 Rytov Approximation 

An alternative to the Born is the Rytov approximation. In this case the 
incident and the scattered fields are represented by 

u.fr) = e* jn (3.18a) 

and 

U,(r) = e* ,!rl (3.18b) 




Using these representations of the fields, the corresponding integral equation 
can be written as: 

^00 = - 7 - / (l V <*s(n) I 2 + Ffr o )l u,(f 0 ) gff]7 0 ) df 0 . (1.19) 

Under the first-order Rytov approximation, the square of the magnitude of the 
scattered field gradient is assumed to be negligible compared to the object 
function, so that the term in brackets can be approximated by Ffr 0 ). 
Therefore the final equation for the scattered phase subject to the first-order 
Rytov approximation is 

d> s m = I F(7 0 ) U|(7 0 ) g(7|7 0 ) d? 0 . (3.20) 

^iCF) object 

A second order solution is obtained through substitution of the gradient of the 
first-order solution in Eq. (3.19). As in the case of higher order Born solutions, 
by applying the same technique repetitively one can come up with higher order 
Rytov solutions. 

It is important to note that, in spite of the similarities between the Born 
(Eq. (3.17)) and the Rytov (Eq. (3.20)) solutions, the approximations are quite 
different [Keller69, SancerTO]. The Born approximation produces a better 
estimate of the scattered field amplitude for large deviations in the refractive 
index of objects which are small in size. On the other hand the Rytov 
approximation gives a more accurate estimate of the scattered phase for large 
sized objects with small deviations in refractive index. When the object is small 
and the refractive index deviates only slightly from that of the surrounding 
media, it is possible to show that the Born and Rytov approximations produce 
the same results [SlanevSl], 

3.2 Conventional Diffraction Tomography 

The basic concept behind most of the present Diffraction reconstruction 
techniques has been derived by Wolf in 1969 [Wolf69], This concept, called the 
Fourier Diffraction Projection Theorem by Pan and Kak [Pan83], can be stated 

as: 

When an object is illuminated with a plane wave as shown in Fig. 

3 1, the Fourier transform of th«' forward scattered fields 
measured on a line perpendicular to the direction of propagation 
of the wave (line TT' in Fig. 3.1) gives the values of the two- 
dimensional Fourier transform of the object along a circular arc. 

It should be noted that the Fourier Diffraction Projection Theorem is valid 







UNIV LAFAVETTE IN SCHOOL OF ELECTRICAL ENGINEERING 
A C KAN 21 DEC 84 DAHD17-82-C-2819 


UNCLASSIFIED 


F/C 28/1 


NL 





























































































only when the object inhomogeneities are weakly scattering, or in fact when 
either the Born or the Rytov approximations hold. 

In the following, we review the derivation of this concept and extend it to 
the case of an attenuating medium, assuming the Born approximation. 
Consider the geometry shown in Fig. 3.2. We assume an incident plane wave 
propagating along the axis 7, expressed by 

u,(r) = U 0 e ik “ Rr,r , (3.21) 

where f/ is the unit vector along the 7 axis and R is the rotational 
transformation applied to the £-7 coordinates to position it along the x-y 
coordinates. 

The solution for the Green’s function in two-dimensional space is a Hankel 
function of zeroth order [Morse53j 


cr 

D 


F|7 0 ) 


— J_ u (•'/ >- 

~ 4 “o l 




If k co is complex, the Hankel function dies out exponentially as a function of 
|7-7 0 |. This becomes obvious when one considers the asymptotic form of 
this function. 


H o 0)(z) - -ie' 2 ' 2 ’ (3.23) 

z 

This behavior agrees with our physical intuition of the problem, since the 
energy from a point source in a lossy medium should decay rapidly as a 
function of the distance. 

The Green’s function is now expanded in terms of plane waves as 
[Morse53], 

: 00 JRXtf-n) 

g(7|? 0 ) = -j- f - da (T~7 0 ) Rri > 0 , (3.24) 

47r -oo P 


where 


A = (M) 


(3.25) 


P = s/Ko 2 -* 2 . (3.26) 

and n and J axes are in the same directions as f and r; are. The condition 
(7-7 0 ) > 0 in Fq. (3.21) restricts the domain of the plane wave expansion 

expressed by that equation to 7 > 0 . Note that when k ro is complex, 3 would 
also be a complex quantity. 1 'hus each plane wave is attenuated in the 7 





































direction as it propagates through the half space tj > 0. 

Substituting Eqs. (3.21) and (3.24) in Eq. (3.17) and interchanging the 
order of integrations, we have 

u s fr) = “ / f Ff7o) e^ R(X_k ^ ,r<> d? 0 da . (3.27) 

47r -00 P object 

Taking the Fourier transform of the field on the line r)-r] 0} we get 

jU _ \ 

U» = ^ F(r (\-K 0 f,)\ , (3.28) 

where 

CO 

u s (a) = / u s (^ 0 ) d? , (3.29) 

-00 

and 

+ I Fff 0 )e- i? “ (? + ^dr 0 . (3.30) 

object 

the ^ and q[ vectors are in the direction of x and y axes, respectively. 

For a lossless medium (real k co ), F as defined in Eq. (3 30) represents the 
two-dimensional Fourier transform of the object function F (because ^ would 
be zero). Considering the relationship between a and j3 as given by Eq. (3.26), 
it is easy to show that one can recover the Fourier transform of the object 
function on the semicircle shown in Fig. 3.3 by employing Eq. (3.28). As the 
object is rotated over 360 degrees, the Fourier domain is filled up (see Fig. 3.4), 
and then a reconstruction can be obtained by the Fourier transform inversion; 
this being the approach initially suggested by Mueller et al. for the 
conventional Diffraction Tomography (Mueller79). 

When the medium is attenuating, (J is non-zero in Eq. (3.30) and that 
equation defines a Laplace transform relationship between the functions F and 
F. The rotation of the object would provide one with a coverage in the four- 
dimensional Laplace domain over a two-dimensional surface. An inverse 
Laplace transformation is possible only if cj is a fixed vector [Bellman66, 
Krylov69, Krylov77], A close look at Eqs. (3.25), (3.26) and (3.27) reveals the 
fact that p and cf are dependent variables and (j would change whenever p is 
altered Therefore, the conventional reconstruction technique can not be 
extended to the case oi an attenuating medium. 















Direction of ' / 
incident field 


Figure 3.3 

The Fourier transform of the scattered field on the receiver line gives the 
values of the object’s Fourier transform over a semi-circle in the Fourier 
domain. 























kv 


Figure 3.4 

By rotation of the object over 360 degrees the object’s Fourier transform is 
recovered over a set of semi-circles inside the circle v2 k,. 






















The object’s cross section can be recovered by taking the inverse Fourier 
transform of the recovered values in the object Fourier domain. To use an 
FFT based algorithm for the inverse transform, one has to know the values of 
the function over a square grid. This calls for an interpolation process from the 
cicular arcs to the square grid. The resulting interpolation errors lead to 
artifacts in the reconstructed object. Different interpolation skims are 
compared in [Pan83j. If the number of projections and the number of sampled 
points per projection are large enough, bilinear interpolation seems to produce 
relatively small artifacts [Pan83j. This can also be justified by considering Fig. 
3.5 for the case where the cicular arcs and the sampled points are very closely 
spaced. 

It has been shown by Devaney [Devaney83] that the Fourier Diffraction 
Projection Theorem can also be used to derive an algorithm that is 
conceptually similar to the filtered-backprojection algorithm for x-ray 
tomography. This algorithm eliminates the need for interpolation in the 
frequency domain and he called it the filtered-backpropagation algorithm. In 
the backpropagation technique the filtered projections are mapped onto the 
object space via an integral transform that is shown to be the inverse of the 
transform that governs phase propagation within the Rytov approximation. In 
spite of the conceptual similarities between the backprojection and the 
backpropagation algorithms, for its full implementation, for an NxN 
reconstruction with N projections, the filtered-backprojection algorithm needs 
only N Fourier transforms. To contrast, the filtered-backpropagation 
algorithm requires about N 2 Fourier transforms. If the extensive demands on 
computer time are not a factor, one can, of course, use the filtered- 
backpropagation algorithm. In most practical applications, however, one may 
have to resort to procedures that are computationally more efficient. 

3.3 Synthetic Aperture Diffraction Tomography 

Nahamoo and Kak [Nahamoo82] have derived a new synthetic aperture 
imaging technique that permits the use of any type of insonification, as 
opposed to only the plane waves and also requires only two rotational positions 
of the object. For each rotational position, the projection data is collected by 
independent traversal movements for the transmitting and the receiving 
transducers (see Fig. 3.6). Although ideally the angle between the two 
rotational positions should be 90°, theory predicts that valid results should be 
obtainable, albeit with reduced spatial resolution, even when this condition is 
not satisfied. For each rotational position of the object, the data can be 
collected most efficiently by using arrays on both the transmit and the receive 






’our points 
e assumed 

















































81 


sides. The elements of the transmit array could be fired sequentially and for 
each such firing the received field would be measured over all the elements of 
the receiver array. In section 3.3.1, we review the theory behind the Synthetic 
Aperture technique in transmission mode, as proposed by Nahamoo and Kak. 
Later in section 3.3.2, we prove that if the system is utilized in both the 
transmission and the reflection modes, only one rotational position of the object 
will be sufficient for the reconstruction of the object's cross section. We 
examine the Synthetic Aperture technique in the presence of attenuation in the 
last section (3.3.3). 

3.3.1 Transmission Mode 

Consider the geometry in Fig. 3.6. When the transmitter is located at the 
point (x t ,y t ), let the incident field U; and the scattered field u 3 at a point 
T = (x,y) be represented by Uj(7;y t ) and u s (7;y t ) respectively. The incident field 
Uj(T;0) can be represented by the angular spectrum expansion (see section 2.1). 


u i(T;°) = / A t (k y ;0) e jR (7 " x ‘ i) dk, 


(3.31) 


where x is the unit vector along the x axis. The wave vector R is given by (see 
Fig. 3.7 for real k co ) 

fv = (k x ,k y ) (3.32) 


k x = V 1 kco* " V 


(3.33) 


The function A t (k x ;0) represents the complex amplitude of the angular 
spectrum of the incident field at the line x=x t when the transmitter is 
positioned at y= 0 . This function is the Fourier transform of the incident field 
on the transmit line as given by 


A t (k ;0) - / Uj(x t ,y;0) e" jk,y dy 


(3.34) 


The angular spectrum of the field when the transmitter is located at y=y t is 
d('fined as 


A t(k r y t ) - / u i ( x t . y : y t ) e~ jk ' y dy 


Since u i (x l .y-y t ;y l ) is the same as U;(x t .y,0), we have 


(3.35) 










































83 


A t (k y ;y t ) = A t (k y ;0) e jk ^ . (3 36) 

Finally the following expansion for the field is found. Now the field at point 7 
when the transmitter is at (x t ,y t ) can be expressed as 


UitfyJ = “ / A t (k y ) e“ Jlc ' y ‘ e - 


jk y y» 


where we have replaced A t (k y ;0) by A t (k y ) for convenience in presentation. 

Going back to the basic integral equation in Eq. (3.17), the Green’s 
function is represented by an angular spectrum expansion. The expansion for 
x >x Q is given by 


; * JA-( 7 -r.J 

^l 7 " 1 = £ / -T" ' 

- X) 

where the wave vector A is given by (see Fig. 3.7 for real k co ) 
A = (a,/?) 


(3.38) 


(3.39) 


a - \Ae? “ ^ 


(3.40) 


Let u s (y r :y t ) represent the scattered field at the point (x r ,y r ) on the receiver line 
when the transmitter is at (x t ,y t ). We substitute equations (3.37) and (3.38) in 
equation (3.17) to obtain 

i 00 --It 1 JforXr + A'r) 

u s(y r ;yt) = f f (g) II A t(M e J ,y ‘ — dk y d ^ 

~(“tr) object -oo 


f j(A-K)r 0 (Jf 


(3-41) 


Switching the order of integrations and rearranging terms, we get 


■ 00 ■ joxr 

, , _ 1 f r je 1 V 

u ’ (y '-■■>' ^ f ! y . ^ 


AJk ) F(A-R) e lf9yr e~ jk ' y ' d/3 dk 


(3.42) 


where F is the function defined in (3.30). Taking the Fourier transform of both 
sides of this equation with respect to y r and v,, we have 


>T(Tk y ) - 


jk.x, jcxx r 
ie e J 


A t (k ) F(A-K) 


where ud d,k y ) has been defined as 







(3.44) 


u s (^,k y ) = If u s (y r ;-y t ) e' jk ^> dy r dy t 

-00 

Let P s (x r ;x t ) be the measured scattered field when the transmitter and the 
receiver are located at points (x t ,y t ) and (x r ,y r ), respectively. We assume that 
the receiver is a flat transducer with a uniform sensitivity over its aperture. 
The effect of the receiver can then be represented by an aperture function 
I r (y r ). The measured scattered field is then given by 


P s(y r Tt) - / u s (y;y t ) I r (v r ~y) dy 

“OO 

A Fourier transform representation of equation (3.45) is given by 
P s (/?,k y ) - A t {8) u s (/?,k y ) , 

where \(3) is the Fourier transform of the aperture function I r (y r ). 
If w'e substitute equation (3.43) in equation (3.46), we obtain 


(3.45) 


(3.46) 


F(A-K) = 


- 2 jo e jk »*‘ e' iaXr 
A t (k y ) 


P.(Akv) 


(3.47) 


Nahammo and Kak [Nahamoo82] have shown that for a lossless medium, 
for which F is the object’s two-dimensional Fourier transform (see Eq. (3.30)), 
Eq. (3.47) defines a one-to-one relationship between the data u 3 and the object’s 
Fourier transform, F, inside the region shown in Fig. 3.8. The values of the 
object's Fourier transform are recovered over the semicircle drawn by the solid 
line in Fig. 3.8 centered at -lx. The coverage area shown in the figure is 
achieved by moving the center of the solid semicircle over the dashed semicircle 
centered at the origin. Therefore, for two rotational positions of the object, 90 
degrees apart, one can recover the object’s Fourier transform inside a circle of 
radius \/ 2 k 0 , as shown in Fig. 3.9. 

Similar to the conventional approach, the values of the object’s Fourier 
transform are calculated on some circular arcs. Therefore, prior to the inverse 
FFT operation, an interpolation technique should be used to recover the values 
of the function over a square grid. An interpolation-free version of the 
Synthetic Aperture technique is proposed in [Nahamoo82], and a very simple 
implementation is presented in [Nahamoo84]. However, the proposed procedure 
requires no interpolation only if the sampling interval on the transmit, the 
receive line and in the object are all the same. This constraint may be very 
restrictive in a practical implementation of such a system. 




























Figure 38 

For a fixed orientation of the object in the transmission mode, the object’s 
Fourier transform is obtained in the region shown. 




• • * i/o • ' x ’ . <■ . * . •*. « F . ,*. .* ^ «\ , 




.‘•'.‘■‘'a' *.* \* ( 
w.' V l* V 























C=AHB 


>/2ko— circle 



frequency domain 


Figure 3.9 

For two rotational positions of the object, 90 degrees apart, the object's 
Fourier transform is recovered inside a circle of radius v2 k 0 . The cover¬ 
age in the region C is double. 


*** */ v - ■' 






























87 


3.3.2 Reflection Mode 


In this section, we show that in the Synthetic Aperture approach, instead 
of rotation of the object, the backscattered field together with the forward 
scattered field can be used to recover the object’s Fourier transform inside a 
circle of radius 2 k 0 , rather than v/ 2 k 0 achieved in the transmission mode with 
two views. Thus by having two fixed parallel transducer arrays on both sides 
of the object and no mechanical movements, the object cross section can be 
reconstructed. In this approach, in addition to the data collected in the 
transmission mode, for all the positions of the transmitter on both of the 
transducer lines, the scattered field on all positions of the receiver on the same 
line is measured. This new set of data provides a coverage in the Fourier 
domain of the object, inside a circle of radius 2 k 0 , which is the complement of 
the coverage achieved in the transmission mode for that view of the object. 

We propose two ways of implementing a system working in reflection 
mode. One approach is to use two transducers for transmission and reception. 
Obviously the scattered field at those points where the transducers overlap 
physically can not be measured. Perhaps the field at those points can be 
estimated bv the average of the field at the neighboring points without causing 
any major distortion in the reconstruction. A second approach is to replace the 
two separate transducers moving on a line by an array of transducers; however, 
the cross coupling between the elements of the array makes the CW 
transmission too erroneous. To avoid the cross-coupling problem, one can 
transmit a long burst of CW pulse, and range gate the received waveform. The 

2 d 

length of the burst, T, has to be shorter than —, where d is the distance from 

c o 

the object to the array as shown in Fig. 3.10, and c 0 is the speed of 
propagation. Before the first echo reaches the transducers, the transmitter is 
switched to the receiving circuitry and can detect the scattered field. On the 
other hand, since the process should be treated as a broadband process, the 
measured scattered field has to be the summation of the scattered fields from 


all the points of (he object, and not only from a certain range of the object. 
Therefore, T, the length of th burst, 

2 S 

has to be chosen so that it is longer than —^ + t^, where S is the size of the 

object (see Fig. 3.10) and tyy is the time required for the quadrature receiver to 
estimate the magnitude and phase of the signal. Thus we have 







































89 



m 




• - 


G 


✓ 

■ h ,-/lVAA\V, 




— + t w < T < —- (3.48) 

C 0 C 0 

D - S 

In this expression d was replaced by —-—, where D (=x r -x t ) is the distance 

<w 

between the arrays. This imposes the following restriction on the object size 

s < J - J tw , (3.40) 

which means that the object size should be at least three times smaller than 
the distance between the transducer arrays. 

The mathematical steps in the derivation of a relationship between the 
reflected data and the object function is quite similar to that of the 
transmission case. The only difference appears in the expansion of the Green’s 
function. Instead of expanding the Green’s function in the forward direction, 
we expand it in the reverse direction 

g(r)T 0 ) = -j- f - d/? , X<X Q ; (3.50) 

4r 4 a 

a is now defined as 

a = s/kj - f (3.51) 

The processing of the backscattered field measured on the transducer line 
at x=x t results in the coverage of the object’s Fourier transform in the region 
shown in Fig. 3.11. We will prove that this mapping is two-to-one, or in other 
words the value of the object's Fourier transform at any point in the region 
shown in Fig. 3.11 can be calculated from two and only two points of the 
processed data u 3 . When the other transducer array at x=x r is used, the 
coverage area would be the symmetric image of the area shown in Fig. 3.11 
with respect to the k y axis. Therefore, the object’s Fourier transform would be 
recovered inside a circle of radius 2k 0 , when both of the transducer lines are 
used in relflection mode, in addition to the normal transmission mode for that 
view of the object. In Fig. 3.12 we have shown this coverage. Notice that the 
mapping between the data and the object’s Fourier transform is one-to-one in 
the transmission mode. However, repeating the data collection procedure for 
the transmission mode by switching the role of the transmit and the receive 
lines would make it a two-to-one relationship as well. 


L, ,f *. ■*. V", ** l O.V JV \ , r -’"S*4* - ■ . *! - ■ * - V v' • ‘ , 

























Figure 3.11 

When the transmitter array in figure 10.6 is utilized only in the reflection 
mode, for a fixed orientation of the object, the object's Fourier transform 
is recovered inside the region shown. 


•jIoULmJUUX a £mA*bJLmJCmj£mJCjUL +JL*jCa tkjCk iC* £+ ftk 




1 




















>v 









































In the following, we show that by employing Eqs. (3.47) and (3.51) in a 
lossless medium, we can recover the Fourier transform of the object function 
inside the circle of radius 2 k 0 centered at the origin and outside the circles of 
radius k 0 centered at (k o ,0) and (—k o ,0), as shown in figure 3.11. We bring the 
proof for one set of the backscattered data (for the region shown in Fig. 3.10) 
while the derivation for the other region follows from the symmetry of the 


problem. 

Let W = (u,v) be defined by 

W = A - R (3.52) 

Eq. (3.51) can be written in the following form when k co is real. 

a 2 + = k 0 2 (3.53) 

Ignoring evanescent waves, the limits on a are found by considering Eq. (3.51). 

- k 0 < a < 0 (3.54) 

Considering Eq. (3.33), one can write the following equation and inequality 

k x 2 + k y 2 = k 0 2 , (3.55) 

0 < k x < k 0 . (3.56) 

Substituting Eq. (3.52) in (3.53) and (3.54), we get 

(u + k x ) 2 + (v + k y ) 2 = k 0 2 , (3.57) 

and 

- k 0 < u + k x < 0 . (3.58) 


Eqs. (3.57) and (3.58) represent a semicircle of radius k D centred at point 
(—k x ,—k y ). We shall denote this semicircle by SC( -k x , -k y ). On the other 
hand Eqs. (3.55) and (3.56) represent a semicircle of radius k 0 centered at the 
origin. This semicircle that contains all the centers of the SC(-k x ,-k y ) is 
depicted by the dashed arc in figure 3.13. A typical SC(-k x ,-k y ) is shown by 
the solid arc in this figure. The union of SC(-k x ,-k y ) for all possible center 
locations (-k x ,-k y ) will provide a coverage of F(W) in some region in the 
Fourier domain. 

Consider the semicircle SC(-k x ,-k y ) shown by the solid line in Fig. 3.13. 
Let A(k x ,k y ) be the center of this semicircle and D(u,v) a point on this 
semicircle. From equation (3.58) we have u <~k x . Since k x > 0 , u has to be 




£2 








YlWiSsWlYlV 




nSAV.n',' 
























Figure 3.13 

In the reflection mode, the object’s Fourier transform is obtained over the 
semi-circles SC(-k t -k f ), one of which is the semi-circle centered at point 
A. By moving this point on the dashed serai-circle, the coverage shown in 
figure 10.11 is achieved. 













94 


negative, thereby the coverage is restricted to the region u < 0 . In the triangle 
OAD, we have 

OD < AD + OA , (3.59) 

or in terms of u and v 

u 2 + v 2 < (2k 0 ) 2 . (3.60) 


This inequality shows that the coverage will be inside the circle centered at. the 
origin with the radius of 2 k 0 . 

Let the point B be the edge of the semicirlce SC(-k x ,-k y ) in the region 
v < 0 . First we consider the third quadrant ( v <0 ) and make use of the 
symmetry of the problem with respect to the u axis and extend it to the second 
quadrant ( v >0 ) as well. It is clear from Fig. 3.13 that the angle DOE is 
always less than or equal to the angle BOE. The relation between the tangent 
of these two angles can be written as 

tan(DOE) > tan(BOE) . (3-61) 

The angle BOE is given by 


BOE = 


APE 

2 


(3.62) 


We also know that 


tan! AQE 1 = t ~ cos (AOE) 

1 2 1 sin(AOE) 

Using equations (3.61), (3.62) and (3.63), we obtain 
sin(AOE) tan(DOE) > 1 — cos(AOE) . 


(3.63) 

(3.64) 


Replacing AOE and DOE, by their expressions in terms of k x , k y , u and v, we 
get 


K 

K 


_u 

V 



(3.65) 


Making use of Eq. (3.57), the inequality above is rewritten as 

u 2 + (v + k 0 ) 2 > k 0 2 (3.66) 

The inequalities (3.60) and (3.66) define the Fourier domain coverage in the 
region u < 0 , v < 0 to be outside the circle centered at ( 0 ,-k o ) with radius k 0 
and inside the circle at the origin with radius k 0 . 





Since the problem is symmetric with respect to the u axis, one can easily 
generalize the result achieved in the third quadrant to the second quadrant. So 
the inequality (3.66) for the whole half plane u < 0 is written as 

u 2 + (| v| + k 0 ) 2 > k 0 2 u < 0 (3.67) 

Finally the coverage area is defined by the inequalities of equation (3.60) and 
(3.67) and is depicted by the shaded area shown in figure 3.11. 

Up to this point, we have proved that any data point maps into the region 
shown in figure 3.11. To prove that the coverage is indeed that region, we also 
have to prove that any point in that region belongs to a semicircle SC(-k x ,-k y ) 
for some k x and k y . 

In the following, we prove that indeed for any point in the shaded area 
shown in figure 3.11, there exists two and only two semicircles SC(-k x ,-k y ). 
We present the proof only for the third quadrant. One can make use of the 
symmetry of the problem with respect to the u axis to extend it to the whole 
half plane u < 0. Let’s take the point D(u,v) in the coverage area where 
v < 0, as shown in Fig. 3.14. The centers of any circle of radius k 0 passing 
through the point D, has to be located at the intersection of a circle of radius 
k 0 centered at D and the circle of radius k 0 centered at the origin, which we 
call C(0,0). Since the distance between the centers of these two circles is less 
than 2k 0 , there always exists two intersection points. We have to prove that 
these intersection points are always in the region u < 0. and also show that 
the u component of D is smaller than that of the two intersection points, so 
that D lies on the left half of the circles passing through it. 

To prove that the two intersection points A and B lie in the left half 
plane, we show that the distance of the point D to any point P on the right 
half of C(0,0) is larger than k 0 . The coordinates of point D(u,v) and P(up,v p ) 


satisfy the following equations 

u 2 + (v + k y ) 2 > k 0 2 ,u < 0 , v <0 and (3.68) 

up 2 + v P 2 = k 0 2 ,u P > 0 ( 3 . 6 ft) 

The distance DP is given by 

DP 2 = (u-up ) 2 + (v v p ) 2 (3.70) 

Using Eq. (3.69), we can rewrite Eq. (3.70) in the following form 

DP 2 = u 2 + (v +k 0 ) 2 - 2uu P - 2v(k 0 + v P ) (3.71) 


Considering the inequality in (3.68) and the fact that the last two terms in the 




















Figure 3.14 

The points A and B are the center of semi-circles SC{-k t -k f ) which 
pass over the point D. Point D is restricted to the region shown. 






equation above are always positive, we obtain 


DP > k r 


(3.72) 


Therefore, the points A and B can not be on the right half of C(0,0) and are in 
the half-space u<0. 

Let’s call the u components of the intersection points A and B, u A and u B 
respectively. It is left to prove that u is smaller than u A and u B . u A and u B 
can be expressed as 

u A = - k 0 sin(DOE - DOA) and (3.73a) 

u B = - k 0 sin(DOE + DOA) (3.73b) 

Since the point D lies in the third quadrant, the angle DOE is limited to the 
following range 


- > DOE > 0 . 

2 - 

We also have the following restriction on the angle DOA 
DOE > DOA > 0 

Using the inequalities in (3.74) and (3.75), we get 


— > DOE - DOA > 0 

2 ~ 


(3.74) 


( 3 . 75 ) 


(3.76a) 


2 ~ > DOE + DOA > 0 (3.76b) 

Equations (3.73a,b) can be combined to obtain 

u B = u A - 2 k 0 cos(DOE) sin(DOA) (3-77) 

Considering the ranges of DOE and DOA (the inequalities in (3.74) and (3.75)), 
the second term in equation (3.77) is always negative, thus we have u B < u A . 
Thereby one only needs to prove that u is smaller than u B . Equation (3.73b) 
can be rewritten as 

n B = - k„ ( sin(DOE) cos(DOA) - cob(DOE) sin(DOA) ) 


(3.77) 


- k„ Mn(l)OE) cos(DOA) 1 + 


tan(DOA) 


tan(DOE) 


























- u OD 


j + tan(DOA) 
tan(DOE) 


= - k„ ■= 


0 OD 2k 0 


tan(DOA) 

tanjDOE) 


From the inequalities in Eqs. (3.74) and (3.75), we have 
tan(DOE) > tan(DOA) > 0 
This inequality can be written into the following form 


1 > 1 
“ 2 


1 + 


tan(DOA) 


tan(DOE) 


> 1 
~ 2 


(3.78) 


(3.79) 


(3.80) 


The two inequalities in (3.78) and (3.80) show that u < u B . The proof is 
complete. 


3.3.3 Lossy Medium 

The theory presented in section 3.3.1 is general, in the sense that it 
includes the case of an attenuating as well as lossless medium. Consider an 
attenuating medium with the complex wavenumber k c0 . The parameters k x 
and a, as given by Eqs. (3.33) and (3.40), become complex, and as a result the 
first argument of F in Eq. (3.47) (a~k x ) also becomes complex. This means 
that by processing the data according to Eq. (3.47), one can find the values of 
the object “Fourier-Laplace” transform over a certain region. By “Fourier- 
Laplace” transform of a two-dimensional function we mean Laplace transform 
in one direction and Fourier transform in the perpendicular direction (see Eq. 
(3.30)). 

Let’s take the vectors p and t[ in Eq. (3.30) to be represented by 

p = u r + v , (3.8ta) 

cj = U; , (3.81b) 

where u = u r + j U; is now a complex variable (see Eq. (3.53)). Let's define 
W = (u.v) as 

W = A - Fv (3.82) 


Eq. (3.40) can be rewritten as 


O 

a“ 


+ r- 



(3.83) 


Substituting Eq (3.82) in (3.83), we have 


























(u + kj- + (v+ky)- - k C g . (3.84) 

On the other hand, k x and k y are related by Eq. (3.33). 

k, 2 + k, 2 = k„ 2 (3.85) 

Eqs. (3.84) and (3.85) can be solved for u ; . The derivation can be found in 
appendix D, the result is 

7 6 Uj 5 + 05 u i 4 + ^4 u i 3 + 7 3 Uj 2 + 7 2 Ui + 7i = 0 , (3.86) 


where 


7s ~ -u r 


75 - 4k 0 e 0 


7 4 = "2u r (u r 2 - v 2 ) 


7 3 = -4k 0 € 0 (‘2u r 2 ~ v-) , 

7 2 = -u r K 4 + v 4 - 4(k 0 2 -£ 0 2 v 2 + 2u r 2 v 2 ) , 


(3.87a) 

(3.87b) 

(3.87c) 

(3.87d) 

(3.87e) 


7i = 4u r 2 k 0 e 0 (u r 2 + v 2 ) . (3.87f) 

The satisfactory solution of U; should result in positive real values for a and k x 
for the transmission mode. Eq. (3.86) defines the coverage surface in the 
three-dimensional space (u r ,u it v). 

Due to the complexity of Eq. (3.86), one has to take recourse to numerical 
approaches to find the solution)s) of U;. To prove or disprove that the mapping 
to the three-dimensional space (u r ,Uj,v) is one-to-one seems to be very involved 
if not impossible. By applying the well known Newton's method, we solved Eq. 
(3.85) for values of u r and v over a square grid numerically and found the 
mapping to be a one-to-two relationship. However, because of the round off 
errors and the small number of points (over the square grid) examined we can 
not generalize our conclusion from the numerical results. Choosing the solution 
with the smallest magnitude, we have plotted the corresponding coverage 
surface in Fig. 3.15, for a typical attenuation constant of 3.0 db/cm. Those 
values of u ; with their magnitude larger than the maximum scale in Fig. 3.15 
were set to zero. As the attenuation goes to zero, the surface in Fig. 3.15 will 
collapse on the u r v plane and the coverage surface will be the same as that 
shown in Fig. 3.8, as expected. The coverage surface in Fig. 3.15 makes 


























e Q = 3.0 (db/cm) 





o J 


Figure 3.15 

The curved surface shows the coverage achieved in the object’s “Fourier* 
Laplace” domain when the medium is attenuating. The attenuation of the 
medium was taken to be 3.0 db/cm. 
























impossible the inversion process to the object domain. To apply the inverse 
Laplace transform techniques the coverage surface has to be parallel to the u r v 
plane [Krylov69, Krylov77, Beilman66j. 

A possible way of avoiding the “Fourier-Laplace” transform operation is 
to impose the following constraint on the processed data [Nahamoo82] 

|*|=|k x | . (3.88) 

Then the coverage area would be limited to the v axis ( | vj < 2k 0 ), instead of 
the region shown in Fig. 3.8. By rotating the object 180 degrees, the object’s 
Fourier transform can be recovered over radial lines (see Fig. 3.16) and the 
object can be reconstructed by means of Fourier inversion. As a matter of 
fact, the limiting constraint defines the intersection of the coverage surface 
shown in Fig. 3.15 with the plane u r v. This is the only reconstruction 
technique available today for an attenuating medium; however, the amount of 
data required makes it impractical with the present technology. 















































103 




a 


J\ 


List of References 


[Ball80] J. J. Ball and F. Stenger, “Explicit inversion of the Helmholtz 

equation for ultrasound insonifica*ion and spherical detection,” 
Acoustical Holography, Vol. 9, A. Metherell, Ed., 1980. 

[Bellman66] R. Beelman, R. E. Kalaba, and J. A. Lockett, “Numerical 
Inversion of the Laplace Transform: Applications to Biology, 
Economics, Engineering, and Physics,” in Modern Analytic and 
Computational Methods in Science and Mathematics, R. E. Bell¬ 
man and R. E. Kalaba, Eds., American Elsevier Publishing 
Company, Inc., 1966. 

[Carter70] W. H. Carter, “Computational reconstruction of scattering ob¬ 
jects from holograms,” Journal of the Optical Society of America, 
Vol. 60, No. 3, pp. 306-314, March 1970. 

[Dandliker70] R. Dandliker and K. Weiss, “Reconstruction of the three- 
dimensional refractive index from scattered waves,” Optics 
Communications, Vol. 1, No. 7, Feb. 1970. 

[Devaney82] A. J. Devaney, “A filtered-backpropagation algorithm for 
diffraction tomography,” Ultrasonic Imagina, Vol. 4, pp. 336- 
350, 1982. 

[Devaney83] A. J. Devaney, "A computer simulation study of diffraction to¬ 
mography,” IEEE Transactions on Biomedical Engineering ,” 
Vol. BME-30, No. 7, pp. 377-386, July 1983. 

[Devaney84] A. J. Devaney, “Geophysical diffraction tomography,” IEEE 
Transactions on Geoscience and Remote Sensing ,” Vol. GE-22, 
No. 1, pp. 3-13, Jan. 1984. 

[Hochstadt73] H. Hochstadt, Integral Equations, John Wiley & Sons Inc., 1973. 

[Iwata75] K. Iwata and R. Nagata, “Calculation of refractive index distri¬ 
bution from interferograms using Born and Rytov’s approxima¬ 
tion.” Japan J. Appl. Phys., Vol. 14, Supp. 14-1, pp. 383, 1975. 


/ V "w* ■-** “w* \vV V* V V V 'V* "j* jC'J 1 m 

' - <Vv ^■ ■“ - N ■"^ -*• - s - N ' * ^ s •" * * 




























104 


[Johnson83] 

[Johnson84j 

[Kaveh79a] 

[Keller69] 

[Kenue82] 

[Krylov69j 

[Krylov77] 

[\lorse53] 

[Morse68] 

(Muel!er79j 

[Nahamoo82] 


S. A. Johnson and M. L. Tracy, “Inverse scattering solutions by 
a Sine basis, multiple source, moment method - Part I: Theory,” 
Ultrasonic Imaging, Vol. 5, pp. 361-375, 1983. 


S. A. Johnson, Y. Zhou, M. K. Tracy, M. J. Berggren, and F. 
Stenger, “Inverse scattering solutions by a Sine basis, multiple 
source, moment method - Part IE: Fast algorithms,” Ultrasonic 
Imaging, Vol. 6, pp. 103-116, 1984. 


M. Kaveh, R. K. Mueller, R. Rylander, T. R. Coulter and M. 
Soumekh, “Experimental results in ultrasonic diffraction tomog¬ 
raphy,” Acoustical Imaging, Vol. 9, pp. 433-450, 1979. 

J. B. Keller, “Accuracy and validity of the Born and Rytov Ap¬ 
proximations,” J. Opt. Soc. Am., Vol. 59, pp. 1003-1004, Aug. 
1969. 


S. K. Kenue and J. F. Greenleaf, “Limited angle multifrequency 
diffraction tomography,” IEEE Transactions on Sonics and Ul¬ 
trasonics, Vol. SU-29, No. 6, pp. 213-217, July 1982. 


V. I. Krylov and N. S. Skoblya, Handbook of Numerical Inver¬ 
sion of Laplace Transforms, Jerusalem: Israel Program for 
Scientific Translations, 1969. 


V. I. Krylov and N. S. Skoblya, A Handbook of Methods of Ap¬ 
proximate Fourier Transformation and Inversion of the Laplace 
Transformation, Moscow: MIR Publishers, 1977. 


P. M. Morse and H. Feshbach, Methods of Theoretical Physics, 
Part I, Chap 7, McGraw-Hill Book Co., 1953. 


P. M. Morse and K. U. Ingard, Theoretical Acoustics, Chap. 7, 
McGraw-Hill, 1968. 


R. K. Mueller, M. Kaveh and G. Wade, “Reconstructive tomog¬ 
raphy and application to ultrasonics,” Proc. IEEE, Vol. 67, No. 
4, pp. 567-587, 1979. 


D. Nahamoo and A. C. Kak, "Ultrasonic diffraction imaging,” 
Purdue University, School of Electrical Engineering, Technical 
Report, TR-EE-82-20, 1982. 











[Nahamoo84] D. Nahamoo, S. X. Pan, and A. C. Kale, “Synthetic aperture 
diffraction tomography and its interpolation-free computer im¬ 
plementation,” IEEE Transactions on Sonics and Ultrasonics, 
Vol. SU-31, No. 4, pp. 218-229, July 1984. 

[Norton81] S. J. Norton and M. Linzer, “Ultrasonic reflectivity imaging in 
three dimensions: Exact inverse scattering solutions for plane, 
cylindrical and spherical apertures,” IEEE Transactions on 
Biomedical Engineering, Vol. BME-28, No. 2, pp. 202-220, 1981. 

[Pan83j S. X. Pan and A. C. Kak, “A computational study of recon¬ 

struction algorithms for diffraction tomography: Interpolation 
vs. filtered-backpropagation,” IEEE Transactions on Acoustics, 
Speech and Signal Processing, Vol. ASSP-31, No. 5, pp. 1262- 
1275, October 1983. 


[Sancer70j M. I. Sancer and A. D. Varvatsis, “A comparison of the Born 
and Rytov methods,” Proc. of the IEEE, pp. 140-141, Jan. 1970. 

[Slaney84] Malcolm Slaney, A. C. Kak and L. E. Larsen, “Limitations of 
imaging with first-order diffraction tomography,” IEEE Transac¬ 
tions on Microwave Theory and Techniques, Vol. MTT-32, No. 
8, pp 860-874, August 1984. 

[Soumekh83] M. Soumekh, M. Kaveh, and R. K. Mueller, “Algorithms and 
experimental results in acoustical tomography using Rytov’s ap¬ 
proximation,” Proceeding of International Conference on Acous¬ 
tics, Speech and Signal Processing, Vol. 1, pp. 135.138, 1983. 

[Tracy83] M. K. Tracy and S. A. Johnson, “Inverse scattering solutions by 
a Sine basis, multiple source, moment method - Part II: Numeri¬ 
cal evaluations,” Ultrasonic Imaging, Vol. 5, pp. 376-392, 1983. 

[Wolf69] E. Wolf “Three dimensional structure determination of semi¬ 
transparent objects from holographic data,” Opt. Comm., Vol. 1, 
pp. 153-156, 1969. 






























CHAPTER 4 

A SIMULATION STUDY OF 
MULTIPLE SCATTERING AND ATTENUATION 
IN DIFFRACTION TOMOGRAPHY 


All Diffraction Tomography algorithms available today, assume the first 
order Born or Rytov approximations; which implies that they are only good for 
reconstructing weakly scattering objects. In a multi-component object this 
translates into ignoring the interaction between the components. However in 
most cases of practical interest the object would be composed of several 
distinct components. Our objective in section 4.1 is to study the effect of the 
interaction between object components on the reconstructed object, in both the 
conventional and the Synthetic Aperture techniques. In this section, we have 
also looked into different approaches to reduce the multiple scattering effects in 
the conventional technique. 

Although the attenuation in tissue is not severe for ultrasound fields, for 
microwave it is on the order of 3 db/cm at 4 GHz. Since there exists no 
practical reconstruction technique for a lossy medium at the present stage, it is 
of importance to study the performance of the current reconstruction 
techniques in the presence of attenuation. We have presented this study in 
section 4.2. 

The data (projections) required for the computer simulations in this 
chapter were provided by the algorithms discussed in chapter 2. 

4.1 Artifacts Caused by Multiple Scattering Phenomena 

As discussed in chapter 9, the scattering interaction between the 
components of an object can cause the scattered field on the receiver line to 
become quite different from the summation of the first order scattered fields. 
First, order fields are calculated by adding the contributions from each 
component while the others are assumed to be absent. Since the reconstruction 
algorithms are based on the weak scattering assumption, the effect of higher 
order scattering contributions present in the projection should appear as a 




















107 


3 


distortion term in the final reconstruction. We consider cases in which the 
reconstruction of any of the object’s components by its own is satisfactory. 

Both the conventional and Synthetic Aperture techniques will be tested to 
observe the sensitivity of each to multiple scattering phenomena. In section 
4.1.2, we consider several possible approaches to reduce the effect in the 
conventional technique. 

4.1.1 Artifacts in Conventional Diffraction Tomography 

The computing procedure discussed in chapter 9 was used to generate a 
certain number of projections over 360 degrees around the object. Each 
projection is the scattered field measured on a receiver line of finite size when a 
multi-component object is illuminated by a plane wave. It should be noted 
that there has been a huge computational effort required to generate the 
projections. For example, generating 64 projections of an object consisting of 
three cylinders took more than a couple of hours of CPU time on the FPS- 
AP120B array processor. 

In all the simulations presented in this section, 64 projections were 
calculated around the object and a bilinear interpolation based algorithm was- 
then used to reconstruct the object cross section; the result follows. 

In the reconstructions shown in Figs. 4.1-4.4, we have shown the 
magnitude of the deviation of the reconstructed refractive index from that of 
the background, which was assumed to be unity. Plots labeled (a) show the 
reconstruction obtained when the projections were generated by ignoring the 
second and higher order scattered fields. On the other hand, the plots labeled 
(b) show the reconstructions when the second order fields are included in the 
projections. 

In Fig. 4.1 the change in the refractive index of the 6\ cylinders was set to 
2 percent. Although there is some distortion introduced in the direction of the 
line joining the center of the cylinders, it is negligible. However, when the 
refractive index change is increased in Fig. 4.2 to 3 percent, the distortion 
becomes quite noticeable; and in Fig. 4.3 a 5 percent change in refractive index 
is enough to cause the distortion to dominate the reconstruction. When the 
number of cylinders is increased, the distortion is higher, as seen in Fig. 4.4. 
This is expected because in this case there are more projections affected by 
second order scattering. 




■pew 






* 










Figure 4.1 

Cross section of a two-component object reconstructed using the conven¬ 
tional diffraction reconstruction technique. Diameter of the cylinders is 6 
X and their refractive index is 1.02. a) Only the first-order scattered fields 
are used for generating the data for this reconstruction, b) Doubly scat¬ 
tered fields are included in the projection data for this reconstruction. 

























13 3 -16 0 ^ 



13 3^ -16 J 


Figure -4.2 

Cross section of a two-component object reconstructed using the conven¬ 
tional diffraction reconstruction technique. Diameter of the cylinders is 6 
\ and their refractive index is 1 03. a) Only the first-order scattered fields 
are used for generating the data for this reconstruction, b) Doubly scat¬ 
tered fields are included in the projection data for this reconstruction. 
































Figure 4.3 

Cross section of a two-component object reconstructed using the onven- 
tional diffraction reconstruction technique. Diameter of the cylinders is 6 
X and their refractive index is 1.05. a) Only the first-order scattered fields 
are used for generating the data for this reconstruction, b) Doubly scat¬ 
tered fields are included in the projection data for this reconstruction. 























> s s v 


B 








15 5 ~ -16 3 










Figure 41 

Cross section of a three-component object reconstructed using the conven¬ 
tional diffraction reconstruction technique. Diameter of the cylinders is 6 
X and their refractive index is 1.05. a) Only the first-order scattered fields 
are used for generating the data for this reconstruction, b) Doubly scat¬ 
tered fields are included in the projection data for this reconstruction. 


... 

. ** v, * *'• •'*. * . -'5 ►V 

-y^y^y-y^y ■*>. 






112 



4.1.2 Different Approaches to Reduce the Artifacts 

We have experimented with three different approaches for reducing the 
distortions caused by multiple scattering in the conventional Diffraction 
Tomography technique. First a statement about the motivation for this 
experiment: in chapter 9, we showed that even when the object components are 
weakly scattering, the multiple scattering effects become important if the 
components are blocking each other. By blocking we mean the line joining the 
component centers is parallel with the propagation vector (see Fig. 4.5a). In 
tomographic reconstructions, these multiple scattering effects exhibited 
themselves as distortions, especially along the lines joining the components. In 
the first approach for reducing these distortions, we use oblique receptions for 
all projections. In this method, the receive array is made non-parallel with the 
transmit array, as shown in Fig. 4.5b. In this way the blocking effect, which is 
quite severe when the propagation vector and the line joining the object 
components (for a two-component object) are collinear, is eliminated. /Vs will 
be demonstrated in section 4.2.1.1, although this arrangement substantially 
reduces the distortions in spaces between the object components, the price paid 
is a small increase in the general distortion all through the object. This 
approach seems appropriate for two-component objects. 

In a second approach intended for more complex objects, the object was 
illuminated with several non-collinear plane waves for each projection, for each 
rotational position of the receive array (see Fig. 4.6). Each plane wave 
generates a complete frequency domain coverage; these are averaged prior to 
Fourier inversion for object reconstruction. The rationale for this method is 
that the errors caused by blocking vary for the three different illuminating 
waves, one reason being that some of them correspond to the oblique reception 
concept mentioned above. Therefore, averaging in the frequency domain 
should diminish the effect of the errors caused by multiple scattering. To 
substantiate this, we have shown, in section 4.1.2.!, reconstructions obtained 
using this idea. 

We have also considered another approach which is similar in spirit to the 
.Synthetic Aperture technique. For each position of the receiver line the object 
is illuminated by incident plane waves from several direct ion- 1 . (However, unlike 
the second approach, no averaging is performed in the frequency domain.) The 
object is now rotated by a large angle to \ e ld a similar coverage over areas left 
blank before: this procedure is repeated until no holes are left in the frequence 
domain. The rationale here is that with each plane wave at a different 
orientation with respect to the receiver line, the severity of the blocking effect 


i 









113 




Figure 4.5 

(a) The condition when the object components block each other is illus¬ 
trated here. 

(b) With the oblique reception illustrated here the "direct blockage" is 
avoided. Direct blockage occurs when the direction of the propagation of 
the incident field is lined up with the line joining the two object com¬ 
ponents and the normal to the receiver array. 































Figure 4.6 

In the second approach several plane waves are used to illuminate the ob¬ 
ject. The frequency domain coverage obtained with the three plane waves 
shown here are averaged prior to inversion. 


























115 



caused by multiple scattering should be reduced. In this case, the amount of 
data required is close to that for the usual diffraction tomography. It was 
found that this approach does result in better reconstruction of multi- 
component objects. 



4.1.2.1 First Approach 

Examine the reconstruction of a two component object in Fig. 4.3; there 
the blocking effect led to distortion that was relatively large in regions between 
the components. It was, therefore, felt that if we instead used for each 
projection the oblique reception scheme shown in Fig. 4.5b, that would 
eliminate the blocking effect, at least for regions between the components. For 
the projection angle shown in Fig. 4.5b and with this approach, the first order 
fields, Ujj and u 2 o, from each component can propagate unperturbed to the 
receiver line, i.e. without any substantial interaction with the other component 
Although the direct blocking (the components being in each other’s shadow and 
also being lined up with the direction of the incident plane wave) is eliminated 
by this method, this leads to reduced reconstruction distortion in spaces 
between the components. Many more projections experience multiple 
scattering now, albeit at a lower level. While for the direct reception method 
there is only one projection where the normals to the transmit array, the 
receive array, and the line joining the centers of the object components are all 
lined up, for the oblique reception case there will be more directions in which 
the line joining the object centers will be perpendicular to either the receive 
array or the transmit array, ,4s a consequence, the oblique reception method 
reduces the reconstruction distortion between the components, but it raises 
somewhat the overall level of distortion. This is illustrated in the section 4.1.2.4 
with computer simulations. 

With regard to the coverage generated in the Fourier domain with the 
oblique reception approach, it is unaffected. When the incident plane wave is 
normal to the receiver line, the Fourier transform of the projections gives the 
values of the Fourier transform of the object on a semicircle, as shown in Fig. 
4.7a. When the incident wave is at an oblique angle to the receiver line, the 
coverage is still on a semi-circle, but positioned differently. This can be seen in 
Fig. 1.7b. Obviously by rotating the object from 0 to 360 degrees, the arc 0.4 
will cover all the area inside a circle of radius r,. where r ( is 

r, = OA { 1.11 


-'A'T'.-'.l-'A' ‘//■•A’ 


•• v /, /. -\ 




sr-^m 


r V V 





116 


. . . 
cy:v»' 


Direction of incident 
plane wave 




Parallel to 
receiver line 


Direction of incident 
plane wave 




Parallel to 
receiver line 



Figure 4.7 

(a) When the incident plane wave is normal to the receiver line, the 
Fourier transform of the projection gives the values of the object's Fourier 
transform over a semi-circle as shown here. 

(b) When the incident wave is at an oblique angle with the receiver line, 
the coverage is still on a semi-circle, although positioned differently. 




'. ■« % 


-V. 













9 is the oblique angle of incidence. A smaller circle also will be covered by 
rotation of the smaller arc OB; the radius of this circle is r 2 where 

QA°_ A 

r 2 = - k 0 sin(—. (4.2) 

Therefore we get a double coverage inside the smaller circle. 

4.1.2.2 Second Approach 

The preceding method can be made somewhat independent of the number 
of object components by using several illuminating plane waves for each 
angular position of the receive array, as shown in Fig. 4.6, and then averaging 
in the Fourier domain the coverages generated by each of the plane waves. 
The rationale for this approach is that averaging in the frequency domain 
would diminish in importance any distorted values that result from blockages. 
With several illuminating waves, it is unlikely that the received fields for all of 
them would suffer from blocking effects at the same time. 

In this approach the coverage in the Fourier domain would be the union of 
the coverage areas for each incident plane wave. As discussed before in the case 
of a single incident plane wave, there are two circles associated with the 
coverage in the Fourier domain. The points inside the smaller circle get double 
coverage, while those in the area between the two circles are covered only once. 
It is now clear that having several incident plane waves, there will be many of 
these circles, precisely twice the number of existing plane waves. Therefore, the 
number of times a point at the center of these concentric circles is covered 
would be equal to the number of circles (twice the number of incident plane 
waves). As the point is moved away from the center of circles, getting out of 
each circle the multiple coverage is reduced by one. Finally, it should be noted 
that the total coverage area (the area in which all the points are covered at 
least once) lies inside the largest of the circles, with the radius of 

OO° + max(0 ; ) 

r - 2 k„ sin ---- (4.3) 

where max(0j) is the largest oblique angle. 












118 


4.1.2.3 Third Approach 

In this approach, we again use a number of plane waves for illumination 
for each rotational position of the object. However, the coverages are not 
averaged in the frequency domain. Instead, we use the procedure described 
below. 

If we use five plane waves 15° apart, the combined coverage in the 
frequency domain for one rotational position of the object is as shown in Fig. 
4.8. The object is now rotated by a large angle to yield a similar coverage over 
areas left blank by the first rotational position, and so on. With very few- 
rotational positions of the object, the complete frequency coverage can be 
obtained (The method is similar in spirit to the Synthetic Aperture approach). 

We felt that this approach would also yield some protection against the 
multiple scattering distortion due to the fact that for each rotational position 
the severity of the blocking effect would be somewhat reduced owing to the 
plane waves being at different angles to the receive line. 

4.1.2.4 Simulation Results and Discussion 

For the oblique reception method, we first modified our computer 
programs for projection data generation to handle any arbitrary direction of 
the incident field. We used this program to generate the data required for the 
simulations done in this section. Further programming had to be done to code 
the reconstruction algorithms used for the results shown here. 

In Fig. 4.9a-d we have shown the reconstruction of a two-component 
object with each component of refractive index 1.05. The components are t 
wavelengths in diameter and 12 wavelengths apart from each other. The 
receiver line is at a distance of 20 wavelengths from the center of coordinates. 
The sampling interval on this line is 0.45 wavelengths and there are 64 points 
in each projection. Plots a.b.c and d in Fig. 4.9 represent the reconstructions 
of the object when the direction of the incident plane wave is. respectively, 0°. 
7.5°, 15° and 30° from the normal to the receiver line. As we expected, the 
larger the off normal angle, the lower the distortion on the line joining the 
components. However, the general level of background distortion in the 
reconstruction increases with an increase in this angle, as is evidenced by its 
severity when the angle is 30°. The reconstructions in plots 1.10b and 4.10c for 
7.5° and 15° respectively, show some improvements over the normal angle case 
in plot 4.9a. 


, L”'" 


v M " » .■ 


'mjf 



























Incident plane waves 



Receiver line 



Figure -1.8 

In the third approach, again a number of plane waves are used to il¬ 
luminate the object, however no independent frequency domain coverages 
are generated from each of these plane waves. For each rotational posi¬ 
tion of the object, the coverage generated by all the illuminating waves is 
extracted first, and then the object is rotated till the frequency domain is 
filled up. 















Figure 4.9 

Reconstructions of a two-component object using the first approach. The 
cylinders of refractive index 1.05 are 4 X in diameter, and are positioned at 
a distance of 12 X from each other. The projections are sampled at 64 
points with a sampling interval of 0.45 X. The receiver line is 20 X from 
the center of the coordinates. The angle between the direction of illumina¬ 
tion and the receiv'.; line is: (a) 0*; (b) 7.5®; (c) 15®; (d) 30®. 


































































122 


The results obtained with the second approach are illustrated in Fig. 4.10. 
Here the frequency domain coverages for each of the plane waves shown in Fig. 
4.6 are averaged prior to Fourier inversion. Three plane waves are considered 
at. angles of -15°, 0° and 15° in Fig. 4.11a, and -30°, -15°, 0°, 15° and 30° in 
plot 4.11b. Some improvement over the first approach is observed, especially 
in plot 4.10a where the edges are enhanced. 

The reconstructions in Fig. 4.11 were generated using the third approach, 
where we fill a part of the Fourier domain for each receiver position and several 
incident plane waves, then rotate the object and do the same thing until a 
complete coverage is obtained. In plot 4.11a, we have considered three incident 
plane waves at -7.5° and 7.5°, and 32 receiver positions were used. In plot 
4.11b the number of plane waves was increased to 5 with illumination angles 
between -15° and 15°. We see a slight improvement in the reconstructions, 
especially in plot 4.11b. When the range of angles is increased to (-30°, 30°), 
the distortion in the direction of the line joining the cylinders is reduced but we 
see an increase in the general background distortion. 

In the second set of simulations, we changed the distance between the 
components to 8 wavelengths and their refractive index to 1.08. Figs. 4.12a-d 
show reconstructions of the object, when the incident plane wave directions are, 
respectively, 0°, 10°, 20 n and 30° from the normal to the receiver line. We see 
the same changes in the case where the components were farther from each 
other. The interesting case is seen in plot 4 12(1. in which the distortion in the 
region between the components is significantly reduced, but of course some 
distortion has been introduced all over the reconstruction area. Since the 
components are close to each other in this case, the receiver line has to be quite 
a bit oblique in this case so that the components may not block each other 
when they are in line with the direction of propagation of the illuminating field. 
In the reconstruction obtained with the second approach shown in Fig. 4.13 the 
recovered Fourier transforms of the object for different angles of incidence are 
averaged. In plots a, b and c there are, respectively, 3 incident plane waves in 
the range (-10° , 10°), 4 in the range (-20° . 20°) and 7 in the range (-30° , 
30°). In plot 4.13c, the improvement is quite noticeable. We have shown the 
results of the third approach in Fig. 4.14. In plot 4.14a. there were 13 receiver 
positions and 5 incident plane waves in range (-20° , 20°) for each of the 
receiver positions. Comparing these plots with 1.13a, we see no improvement. 
As soon as the range of incident angles is increased to (-30° , 30°) tn plot 4.14b. 
the distortion between the cylinders is reduced. These results lead us to the 
conclusion that for closer spacing of components, the angle of off-normal 
incident directions should become larger. We also see a little improvement in 
























! 


C b i 


Figure 4.10 

Reconstructions of a two-component object using the second approach. 
The cylinders of refractive index 1.05 are 4 X in diameter, and positioned 
at a distance of 12 X from each other. The projections are sampled at 64 
points over a sampling interval of 0.45 X. The receiver line is 20 X from 
the center of the coordinates, (a) 3 incident plane waves, with the angle of 
illumination between -15* and 15*. (b) 3 incident plane waves, with the 
angle of illumination between -30* and 30*. 












Figure 4.11 

Reconstructions of a two-component object using the third approach. The 
cylinders of refractive index 1.05 are 4 X in diameter, and positioned at a 
distance of 12 X from each other. The projections are sampled at 64 
points over a sampling interval of 0.45 X. The receiver line is 20 X from 
the center of the coordinates, (a) 3 incident plane waves, with the angle of 
illumination between - 7 . 5 ° and 7.5® and 32 receiver positions, (b) 5 in¬ 
cident plane waves, with the angle of illumination between -15* and 15® 
and 32 receiver positions, (c) 5 incident plane waves, with the angle of il¬ 
lumination between -30® and 30* and 16 receiver positions. 




. V V V,V V V v ' 





































< b ' 


Figure 4.12 

Reconstructions of a two-component object using the 6rst approach con¬ 
sisting of oblique reception. The cylinders of refractive index 1.08 are 4 X 
in diameter, and positioned at a distance of 8 X from each other. The pro¬ 
jections are sampled at 64 points over a sampling interval of 0.45 X. The 
receiver line is 20 X from the center of the coordinates. The angle between 
the direction of illumination and the receiver line is: (a) 0”; (b) 10”; (c) 
20”, (d) .50”. 






.\AV 

































Figure 4.13 

Reconstructions of a two-component object using the second approach. 
The cylinders of refractive index 1.08 are 4 X in diameter, and positioned 
at a distance of 8 X from each other. The projections are sampled at 64 
points over a sampling interval of 0.45 X. The receiver line is 20 X from 
the center of the coordinates, (a) 3 incident plane waves, with the angle of 
illumination between -10° and lO”. (b) 5 incident plane waves, with the 
angle of illumination between -20' and 20'. (c) 7 incident plane waves, 
with the angle of illumination between -30" and 30 <> . 


















































Figure 4.14 

Reconstructions of a two-component object using the third approach. The 
cylinders of refractive index 1.08 are 4 X in diameter, and positioned at a 
distance of 8 X from each other. The projections are sampled at 64 points 
over a sampling interval of 0.45 X. The receiver line is 20 X from the 
center of the coordinates, (a) 5 incident plane waves, with the angle of il¬ 
lumination between -20® ana 20® and 16 receiver positions, (b) 7 incident 
plane waves, with the angle of illumination between -30® and 30® and 8 
receiver positions. 






















the results from the second approach over that of the third, but the amount of 
data required for the second one is much higher. 

We also tested the different techniques for a three-component object; the 
results are shown in Fig. 4.15-17. There is just one component added to the 
geometry of the first simulation. Plots 16a-d show the reconstruction of the 
object when the incident plane wave direction is, respectively, 0°, 10°, 20° and 
30° from the normal to the receiver line. The improvement in the 
reconstruction is not as noticeable as in the case for a two-component object, 
but still the same effects are seen. The same thing can be said about the 
results of the second approach shown in the next figure. In the reconstruction 
shown in Fig. 4.16. the recovered Fourier transform of the object for different 
angles of incidence is averaged. In plo s a, b and c there are, respectively, 3 
incident plane waves in the range (-10° . 10°), 5 in the range (-20° , 20°) and 7 
in the range (-30° , 30°). We have shown the results of the third approach in 
Fig 4.17. In plot 4.17a, there were 16 receiver positions and 5 incident plane 
waves in the range (-20° , 20°) for ea< h of the receiver positions. When the 
range of angles is increased to (-30° , 30°) in plot 4.17b, the distortion in the 
regions between the components is reduced. 

4.1.3 Artifacts in Synthetic Aperture Diffraction Tomography 

For simulation purposes, the transmitting transducer can be assumed to be 
a point source. Even if the transducer is not an isotropic radiator, the 
deconvolution process in Eq. (3.47) (devision by A t (k y )) would make it function 
as a point source. Therefore by choosing a point source as the transmitting 
transducer, one can avoid the deconvolution computations and the difficulties 
associated therewith. 

Similar to the conventional technique, the Synthetic Aperture technique 
also calls for the interpolation to a square grid in the Fourier domain. An 
interpolation-free version is proposed in [Nahamoo82], and implemented in 
[NahamooSlJ. However, the proposed procedure requires no interpolation only 
if the sampling interval on the transmitter line, receiver line and in the object 
are all the same. Taking recourse to the Fourier domain interpolation 
technique, the values of the object's Fourier transform are obtained only over a 
finite number of semicircles in the Fourier domain, which calls for some kind of 
interpolation scheme. The approach used in [Nahamoo82] was to find the 
neighboring available values of k y , assume the same value of 3 for each of the 
neighboring k y values and find its neighboring available values of 3. In Figs. 
4.18(a,b), the points A and .V . B and TV, and C and C have the same values 






































155 ' -16 0 


Figure 4.15 

Reconstructions of a three-component object using the first approach con¬ 
sisting of oblique reception. The cylinders of refractive index 1.05 are 4 X 
in diameter, and positioned at a distance of 6 X from the center of coordi¬ 
nates. The projections are sampled at 64 points over a sampling interval 
of 0.45 X. The receiver line is 20 X from the center of the coordinates. The 
angle between the direction of illumination and the receiver line is (a) 0*, 
(b) 10*, (c) 20*, (d) 30*. 



























































Figure -*.16 

Reconstructions of a three-component object using the second approach. 
The cylinders of refractive index 1.05 are 4 X in diameter, and positioned 
at a distance of 6 X from the center of coordinates. The projections are 
sampled at 64 points over a sampling interval of 0.45 X. The receiver line 
is 20 X from the center of the coordinates, (a) 3 incident plane waves, 
with the angle of illumination between -lO 9 and lO®. (b) 5 incident plane 
waves, with the angle of illumination between -20" and 20*'. (c) 7 in¬ 
cident plane waves, with the angle of illumination between -30* and SO". 

































Figure 4.16, continued. 


“* k*" «■' . 


* - v /-vat- v' - v 




























Figure 4.17 

Reconstructions of a three-component object using the third approach. 
The cylinders of refractive index 1.05 are 4 X in diameter, and positioned 
at a distance of 6 X from the center of coordinates. The projections are 
sampled at 64 points over a sampling interval of 0.45 X. The receiver line 
is 20 X from the center of the coordinates, (a) 5 incident plane waves, 
with the angle of illumination between -20* and 20® and 16 receiver posi¬ 
tions. (b) 7 incident plane waves, with the angle of illumination between 
-30® and 30® and 16 receiver positions. 


. -V ^*J ✓ J" s w* S J* S .* V .» ,* - ,* 




A. 1 ** 





























of 3, respectively. The values of k Y in Fig. 4.18a and 4.18b were chosen to be 
close to k 0 and zero, respectively. If a bilinear interpolation scheme is used, the 
four neighbors of a point seem to be the edges of a semi-parallelogram, as 
expected. However, a nearest neighbor scheme may result in the completely 
wrong nearest neighbor of the point of interest because of the positions of 
points A and A', B and B' , and C and C 1 relative to each other. In Fig. 4.19 
we have shown a modified interpolation scheme. Point C is the point at which 
the value of the object’s Fourier transform is sought. After finding the 
neighboring values of k y (points O, and 0 2 ), the intersections of the circle 
going through point C and the two semicircles centered at O, and 0 2 are found 
(points C, and C 2 ) and, finally, the neighboring points reached (points A,. B,. 
A 2 and Bo). 

The old and the new interpolation techniques are compared in Figs. 4.20 
and 4.21 for both the nearest neighbor and bilinear schemes. We h ave used the 
computational technique discussed in chapter 9 to generate the required data 
for the object reconstruction. The number of array elements in both the 
receiver and transmitter has been doubled in the simulation result in Fig. 4.21 
compared to that in Fig. 4.20. As expected there is a noticeable improvement 
in the modified nearest neighbor scheme compared to the old nearest neighbor 
technique. But there seems to be no change when the bilinear interpolation 
technique is used. The results of the old and the new techniques become more 
similar when the number of array elements is increased, which is intuitively 
justifiable. 

In the formulation of the Synthetic Aperture technique, the transmit and 
the receive lines were assumed to possess infinite length. However, in the 
simulation environment, the finite size of the transducer arrays may cause some 
distortion in the final reconstruction. Let's first consider the receive line. For 
those transducer points on the edge of the receiver line which are at the far 
field, the change in the field observed by moving to a neighboring point is a 
simple phase change due only to the propagation effect of the wave. If the 
receive line is limited in size to a length several times larger than the object's 
size, the efTect would appear as a lower energy at the spatial frequencies near 
k 0 , which corresponds to the center frequency of the incident field. This 
happens because the high frequency oscillations of the field at the end of the 
receiver line is simply zeroed. Since k 0 is the largest spatial frequency 
considered in the reconstruction technique, thereby one can associate the finite 
size of the receive line with the following limitation on 3 in Fq. (3.47). 






























Figure 4.18 

A simple interpolation scheme for the Synthetic Aperture technique 
Points .4 and A' , B and B , and C and C lying on two different semi 
circles SC(-kj-k ) posses equal values of J, respectively. 

(a) small values of k t . 



























parallel lines 



(b) 


Figure 4.18, continued 



















Figure 4.19 

Points C, and C 2 are the intersection points of the semi-circles O x and 0 2 
with the circle centered at the origin going through point C. The neigh¬ 
boring points of C, and C\ are chosen for the kind of interpolation 
desired. 
















nnb mtarpolat ion . 32 alaaam irr^i 



modified nob interpolation . 32 element arrays 



Figure 4.20 

Magnitude of the reconstructed refractive index of a cylinder with diame¬ 
ter of 6 X and refractive index of 1.03, using the Synthetic Aperture tech¬ 
nique. Transmitter and receiver arrays, positioned 20 X apart, have 32 ele¬ 
ments with inter-element distance of 0 45 X. (a) Simple nearest neighbor 
interpolation, (b) Modified nearest neighbor interpolation, (c) Simple bil¬ 
inear interpolation, (d) Modified bilinear interpolation. 































bilinear interpolation / 32 element arrays 



< c ) 


modified bilinear interpolation . 32 element arrays 



Figure 4.20, continued. 

























rmb Interpolation , 64 alaoant array* 



( a ) 

modified nnb interpolation , 64 element arrays 



< b ) 


Figure 1.21 

Magnitude of the reconstructed refractive index of a cylinder with diame¬ 
ter of 6 X and refractive index of 1.03, using the Synthetic Aperture tech¬ 
nique. Transmitter and receiver arrays, positioned 20 X apart, have 64 ele¬ 
ments with inter-element distance of 0.45 X. (a) Simple nearest neighbor 
interpolation, (b) Modified nearest neighbor interpolation, (c) Simple bil¬ 
inear interpolation, (d) Modified bilinear interpolation. 









bilinear mterpolat ion . 64 element 



modified bilinear interpolation . 64 element arrets 



a l * ^ -tz t ^ 

4 


Figure 4.21, continued. 































| 3\ < ak 0 . (4.4) 

where a is a factor less than unity. 

A similar argument can be used for the transmit array, except this time 
the received field for the transmitter positions at the edge of the transmit array 
are multiplied by a phase factor when one moves to the neighboring 
transmitter position. A limitation similar to that in Eq. (4.4) would be imposed 
on the magnitude of k y 

| k y j < a k Q (1.5) 

The magnitude of the scattered field from a cylinder with refractive index 
of 1.03 and diameter of 6 wavelengths, as a function of the transmitter and the 
receiver positions, is plotted in Fig. 4.22. The effect of the finite size of the 
arrays is seen as the chopping of the signal at the end of array positions. The 
magnitude of the Fourier transform of the field with respect to the transmitter 
and the receiver positions is shown in Fig. 4.23. In plot (a) the number of 
array elements is 32, while in plot (b) the array elements have been increased 
to 128 elements. The clearest advantage in plot (b) relative to plot (a) is the 
recovery of higher frequencies at the end of the diagonal line (higher values of 
k y and 3). We expect a constant value at those points for which k y = 3, 
because that corresponds to the value of the object’s Fourier transform at the 
origin. To correct this effect as observed in Fig. 4.23b, we normalize the result, 
taking the value at k y = 3 = 0 as the reference point. This point was chosen 
because it is the lowest frequency component available and, therefore, is less 
distorted. 

It can be easily shown that the limits in Eqs. (4.4) and (4.5) cause the 
coverage in Fig. 3.8, to shrink to the one shown in Fig. 4.24. Finally the 
coverage achieved for two rotational views of the object 90° apart, would be 
inside a circle of radius 

r = \/2 k 0 ( a - 'Z l~ a 2 ) ; a>— , (4.6) 


r = y/2 k 0 ( a - VI - a 2 ) ; a > , (4.6) 

rather than \/2k 0 . To avoid the loss of resolution, one can collect the data for 
several views of the object. A typical coverage obtained for three rotational 
views of the object 60° apart is shown in Fig. 4.25. The coverage in this case 
would be inside a circle of radius 


r - k r 


a > — 
— o 


• J. * "O 











146 


n = 1 03 , € 0 =c=00 (db'cm) 



Figure 4 22 

agnitude of the scattered field from a lossless cylinder, embedded in a 
■isless medium, as a function of transmitter and receiver positions. Di- 
leter and refractive index of the cylinder are 6 X and 1.03, respectively. 



























v v vi v 


n ■ 103. « Q * c ■ 0 d (db/cm) 32 il tit nt 4rr«^$ 


151 


2 


^ *4 S 

fc 



n » 1 03 t c Q * c * 0 0 Cdb-'cm) # 128 


4 

C 937 



Figure -*23 

Magnitude of the Fourier transform of the scattered field from a lossless 
cylinder, embedded in a lossless medium, with respect to transmitter and 
receiver position. Diameter and refractive index of the cylinder are 6 X 
and 1.03, respectively, (a) 32 element transmitter and receiver arrays, (b) 
128 element transmitter and receiver arrays. 


























































: "Tf ■' 


■ m ■ j i 


f ^ iH7 i 


149 


k-.v. 





• * 



Figure T25 

When the high frequency information in the projections is lost, larger 
number of views can increase the resolution and decrease the distortion in 
the reconstructed image. The coverage shown corresponds to three rota¬ 
tional views of the object 60® apart. 


tc.V.’ 






In Fig. 4.‘26, we have shown the reconstructed Fourier transform of a 
cylinder with refractive index of 1 03 and diameter of (> wavelengths. The 
number of array elements in plots (a) and (b) were taken to be 32 and 12*. 
respectively. The severe distortion in the diagonal direction agree* with the 
concept explained by Fig. 1.24. In Figs. 1.27 and 1.28 we have shown the real 
(refractive index) and the imaginary (difference in attenuation constant) parts 
of the reconstructed object function for different sizes of the array*, 
respectively. Of course, the reconstructed attenuation function should liava* 
been zero. However, due to the large size of the object and large value of the 
refractive index, the large estimation error in the phase causes a transfer of 
energy from the reconstructed refractive index ter the reconstructed attenuation 
function. Since the medium attenuation is zero, the magnitude of the 
reconstructed object function would have been a better estimate of the object's 
refractive index. We have chosen the real and imaginary parts only for a 
better two-dimensional presentation of the distortion. The distortion along the 
diagonal direction in the reconstructed object is again due to the loss of the 
high frequency components of k v and J. One can get rid of the distortion by 
increasing the size of the arrays, as shown in Figs. 4.27c and 4.28c. The 
simulation results show that the loss of coverage, shown in Fig. 4.24, is indeed 
correct and a larger number of elements in the transmitter and receiver would 
solve the problem. 

The simulation study presented in section 4.1.2 persuaded us that the 
performance of the Synthetic Aperture technique in the presence of multi- 
component interactions would be superior to that of conventional diffraction 
tomography. The programs developed, based on the techniques discussed in 
chapter 9, were used to generate the required projections from a multi- 
component object. The computational time for a three cylinder object was on 
the order of 14 hours of CPU time in the array processor FPS-AP120B. The 
reconstructions shown in Figs. 4.29-31, for the Synthetic Aperture technique 
correspond to the ones presented in Figs. 4.2-4 for the conventional technique. 
The number of array elements was taken to be 128, to avoid the effect of the 
finite size of the arrays. The arrays were positioned l’> wavelengths away from 
the center of coordinates and the inter-element distance was set to 0.4) 
wavelengths. The results look very interesting because there is almost no 
distortion along the line joining the center of cylinders as was the case in the 
reconstructions achieved using the conventional technique. The reason is that 
the blocking effect is minimized in this approach. Another interesting point 
about the reconstructions shown in Figs. 4.29-31 is that the components 
blocked are more distorted compared to the ones which are closer to the 




KV* 




n*l03- « o »c*00 (db-'Cm) 32 arrays 




Figure 4.26 

Magnitude of the reconstructed object’s Fourier transform in a lossless 
medium. The object is a lossless cylinder with diameter and refractive in¬ 
dex of 6 X and 1.03, respectively, (a) 32 element transmitter and receiver 
arrays, (b) 128 element transmitter and receiver arrays. 
























S'**.# I 




6H »l*«#m *rr*gj 




Figure 4 27 

Real part of the reconstructed object function of a cylinder with diameter 
of 6 X aDd refractive index of 1.03, using the Synthetic Aperture tech¬ 
nique. Transmitter and receiver arrays, with the inter-element distance of 
0.45 X, are positioned 20 X apart, (a) 32 element arrays, (b) 64 element 
arrays, (c) 128 element arravs. 


































































3 2 el**#ru arrays 


o 

w 




\l ■» 


C a ) 


64 element, arrays 


o 

w 

e 


1« H 


( b > 


Figure 1 28 

Imaginary part of the reconstructed object function (difference in attenua¬ 
tion constant of object and medium) of a cylinder with diameter of 6 X 
and refractive index of 1.03, using the Synthetic Aperture technique. 
Transmitter and receiver arrays, with the inter-element distance of 0.45 X, 
are positioned 20 X apart, (a) 32 element arrays. (b) 64 element arrays, 
(c) 128 element arrays. 







































C & ) 


Figure I 29 

Magnitude of the reconstructed refractive index of a two-component ob¬ 
ject, using the Synthetic Aperture technique. The components are 
cylinders with diameter of 6 X and refractive index of 1.03. Transmitter 
and receiver arrays positioned 30 X apart, have 128 elements with the 
inter-element distance of 0.45 X. a) Only the first-order scattered fu Ids are 
used for generating the data for this reconstruction, b) Doubly scattered 
fields are included in the projection data for this reconstruction. 


























r igure 4 .30 

Magnitude of the reconstructed refractive index of a two-component ob¬ 
ject, using the Synthetic Aperture technique. The components are 
cylinders with diameter of 6 X and refractive index of 1.05. Transmitter 
and receiver arrays positioned 30 X apart, have 128 elements with the 
inter-element distance of 0.-45 X. a) Only the first-order scattered fields are 
used for generating the data for this reconstruction, b) Doubly scattered 
fields are included in the projection data for this reconstruction. 

























Figure 4.31 

Magnitude of the reconstructed refractive index of a three-component ob¬ 
ject, using the Synthetic Aperture technique. The components are 
cylinders with diameter of 6 X and refractive index of 1.05. Transmitter 
and receiver arrays positioned 30 X apart, have 128 elements with the 
inter-element distance of 0.45 X. a) Only the first-order scattered fields are 
used for generating the data for this reconstruction, b) Doubly scattered 
fields are included in the projection data for this reconstruction. 






















159 


transmitter. Note that the transmitter array has been positioned on the 
negative side of the axes x and y in those figures. In the case of two- 
component objects the effect is quite dear. In Fig. 1.51 (threc-oomponet 
object), the cylinder which is farther away from the transmitter array in both 
the projections, is distorted more than the other two, which seem to hold the 
same amount of distortion. 

4.2 A Simulation Study of Attenuation in Diffraction Tomography 

Due to the lack of any practical Diffraction Tomography technique for a 
lossy medium, a study of the performance of the present techniques in the 
presence of attenuation was necessary. Although in ultrasonic imaging the role 
of attenuation is minor, in microwave imaging it can not be ignored, because 
microwaves at 4 GHz undergo almost 3 db of attenuation per centimeter of 
travel in water. 

Let's consider the case where the object and the medium have different 
attenuation constants. Therefore, the wavenumber in the object (k c fr)) and the 
medium (k co ) are both complex. k c (7) is defined as 

k c (T) = k(T) + j c(r) . (4.8) 

The reconstruction algorithms based on the Fourier Diffraction theorem 
estimate the complex object function defined in Eq. (3.12). Since n(7), as 
defined in Eq. (3.10) can also be expressed in terms of the wavenumbers k c (T) 
and k ro as 


nfr) - —— , ( 4 - 9 ) 

kco 

one can rewrite Eq. (3.12) in the following form 

Ffr) = k c 2 - k co 2 . (4.10) 

Since k r0 is known, both the refractive index and the attenuation of the object 
can be estimated simultaneously. 

Substituting the expressions given by Eqs. (2.1) and (4.8) in Eq. (4.10) and 
restricting T to the interior of the object, we get 

bfr) = k 0 2 - k J (7) + c-(T) - c Q 2 + j 2 ( k 0 c 0 - kfr) ((f) J 

— - k n ( k„ ~ k(T) j + j 2 ( k„ c„ - k(T) c(T) ) (4.11) 

The last approximation was applied because the imaginary part of the 


r: -^v?v:v vivvi^ 













wavenumber is much less than the real part in all cases of interest. 
Considering Eq. (4.11), one can see that even if the object inhomogeneities are 
negligible (k = k 0 ), the difference between the attenuation coefficients of the 
object and medium (c 0 - c(7)) may make the imaginary part (and thus the 
magnitude) of the object function large enough so that the Born approximation 
would break down. Therefore, we expect poor reconstructions for large 
differences in the attenuation coefficient, even for very small differences in the 
refractive index of the object and medium. Perhaps this would not happen in 
microwave imaging because the attenuation constants of water-based materials 
are very close. It should be noted that since the estimation parameter in Eq. 
(4.11) is a complex function, any error in either the refractive index (k(7)/k 0 ) or 
the attenuation (c(7)) would affect the estimated value of the other. 

In all the simulation results in later sections, we will refer to the refractive 
index as the ratio of k(7) to k 0 (the ratio of the speed of the wave in the 
medium to that in the object), rather than the complex function nfr) from Eq. 
(4.9). 

4.2.1 A Simplified View of the Problem 

In chapter 10 we showed that due to the lack of an inverse Laplace 
transform matching the domain of the coverage, the current reconstruction 
techniques can not be extended to the case of an attenuating medium. If one is 
not dealing with high attenuation in the medium, we will show that it is 
possible to approximate the Laplace transformation by a Fourier transfomation 
over a sub-domain of the coverage. We will consider both the conventional 
technique and the Synthetic Aperture technique. 

Consider Eq. (3.28) derived for the conventional technique, expressing the 
relation between the Fourier transform of the scattered field on the receiver 
line and the object's Laplace Transform. The Laplace transformation can be 
approximated by a Fourier transformation only if the imaginary part of 

A - k co i/ = ( a , /j - k co ) (4.12) 

is negligible (see Eq. (3.28)). Since a is real, only the imaginary part of J~k ro 
has to be small. Let's define a and f n as fractions of k 0 

o = a k r , , (4.13) 

C = b k Q (4.It) 

For microwave in tissue, the upper limit on b is about 0.07. .?-k co can then be 
approximated by 





















161 


3 ~ k rn - J k 2 - or" - k r 


1 - a- + j 2 b - k c 


(4.15) 


Limiting a to values of less than 0.7, approximately. Eq. (4.15) can still be 
further simplified 


- a - - 1 


3 ~ k co - k o ( V 1 ~ a- - 1 ) + j a' £ 0 (4.16) 

This restriction on a is not required, and the only reason we have applied it is 

to make the mathematics much simpler. The imaginary part of /?-k co in the 
equation above is linearly proportional to a 2 or in fact a 2 . Therefore, for values 
of a smaller than an upper bound, the imaginary part of 3~ k co is negligible and 
the Laplace transformation in Eq. (3.28) can be approximated by the Fourier 
transformation. This point has been mentioned in an unpublished work of 
Mueller [Mueller]. 

In the following we look at an interesting concept, which shows that the 
limitation on ci as suggested in the last paragraph is indeed a physical 
constraint on a in an attenuating medium. As we mentioned in section 2.1, the 

angular spectrum of a field on the line x = x, propagates to line x = x 2 

according to the following relation 


A(k x 2 ) = A(k X|) e J 




(4.17) 


Fig. 4.32 shows a plot of the magnitude of the exponential factor in Eq. (4.17) 
as a function of both the attenuation of the medium and the spatial frequency 
of the plane wave in the y direction. We have plotted the magnitude in db as 
a function of the dimensionless parameter 7 , where k x = 7 k 0 . Thus, for 7=0 
the wave is traveling directly towards the receiver line, while for q~l the wave 
is propagating along the direction of the receiver line. In this plot the 
attenuation factor and 7 have been changed from 0 to 5.0 db/cm and 0 to 1.0, 
respectively. It is clear that the high frequency components (larger 7 ) are more 
attenuated than those components at lower frequencies (those that travel 
directly towards the receiver line). This means that as the field propagates in 
an attenuating medium, it loses its high frequency components. 

Remembering the Fourier domain coverage of diffraction transmission 
tomography, one can now associate this phenomenon with a degradation in 
resolution. This point is illustrated in Fig. 4.33. As in the case of no 
attenuation, the inner circle corresponds to transmission tomography, while the 
outer ring represents the data measured by a reflection tomographic 
experiment. The difference made by attenuation is shown as the shaded area. 


-y -\y 
V.V.V. V.v.V>.V. 


' *’ * •/.. A. ^ . A A **k *"| 




























Figure t32 

Energy loss of plane wave components resulting from an angular spectrum 
expansion in a lossy medium as a function of directional angle and at¬ 
tenuation constant of the medium. 7 is a measure of the direction of pro¬ 
pagation of these plane waves. 7 = 0 corresponds to the perpendicular 
direction to the line of expansion, and 7=1 corresponds to the parallel 
direction to that line. 





























Figure 4.33 

Object’s Fourier domain coverage for the conventional tomography in 
transmission and reflection modes, in a lossy medium. The lower resolu¬ 
tion in the conventional transmission tomography is due to the smaller ra¬ 
dius of the coverage circle as shown here. 



























In this region the attenuation of the medium reduces the amplitude of the 
plane wave components below a minimum tolerable Signal to Noise Ratio 
(SNR) and thus makes this region of the object's Fourier transform 
unmeasurable. 

In the Synthetic Aperture technique, the concept of propagation of the 
angular spectrum in an attenuating medium can be applied again, and this 
results in upper bounds on k y and J. k y is the spatial frequency associated 
with the transmitter, and 3 is the spatial frequency associated with the 
receiver. This is very similar to the effect of the finite size of the transmitter 
and receiver arrays, as shown in Fig. 4.24. This point is illustrated in Fig. 
4.34. The coverage in the reflection mode also shrinks as shown in the figure. 


4.2.2 Attenuation in Conventional Diffraction Tomography 

In the absence of a reconstruction procedure specifically designed for an 
object in a lossy medium a study of the performance of the present algorithms 
in the presence of attenuation is in order. Since the quantitative aspect of the 
reconstruction is of importance to us, we have compensated for the loss of 
energy in the propagation of the field to the receiver array. In the following we 
present reconstructions achieved in this way. based on the projection data 
calculated using the programs discussed in the previous sections. 

In Figs. 4.35-47 we have shown reconstructed values for the refractive 
index and the attenuation for a number of different cylinders and media. In 
each case the first plot (a) represents the difference in the refractive index of 
the cylinder and the surrounding medium and the bottom plot (b) shows the 
reconstructed difference in the attenuation constants. In the first set of Figs. 
4.35-40, the attenuation constant of the object and medium are taken to be 
equal. Considering these figures, we conclude that higher values of attenuation 
result in the loss of higher frequencies in the reconstruction and result, in a loss 
of resolution. The non-zero estimates of the attenuation constant are due to 
the error in the phase of the object function. 

In Figs. 4.41 through 4.44 we have considered different values for the 
attenuation constant of the object am! the medium, while in all cases we have 
kept the refractive index of the object and the medium equal. Again, higher 
attenuation of the medium causes the loss of high frequency components in the 


reconstruction. 


Now, let’s examine a lossy object in a lossless medium. In Figs. 4.45 and 
4.46. we have shown reconstructions of both the refractive index and the 
attenuation constant of cylinders with refractive index of 1.01 and attenuation 




. *v.■'V*. v -VV n/ V 

* - * J 


\ A „*» 

» „*• * J* „ « m r 





















Figure *-3t 

When medium is lossy, larger number of views in Synthetic Aperture tech¬ 
nique can increase the resolution and decrease the distortion in the recon¬ 
structed image The coverage shown corresponds to three rotational views 
of object, 00 ° apart. 


** *“ rf . .r . rfw . fC.V.Y.ir. <.y. y- V v v. vj 




































Figure 4 35 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1 .01, using the conventional technique. Both the object and the medium 
are lossless. The projections are sampled at 64 points over a sampling in¬ 
terval of 0.45 X. The receiver line is at a distance of 10 X from the 
cylinder, (a) Refractive index, (b) Difference in attenuation constant of 
the object and the medium. 























n = 1 0 1 


1 0 i db- cm ) 



15 5 ~ -16 0 



Figure \ .55 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.01, using the conventional technique. The attenuation constant is 1.0 
db/cm in both the object and the medium. The projections are sampled 
at 64 points over a sampling interval of 0.45 X. The receiver line is at a 
distance of 10 X from the cylinder, (a) Refractive index (b) Difference in 
attenuation constant of the object and the medium. 


V •/* . 































■ jii am riiMiii • 


I ija Hll’j 


1 3 ! tl 


UNCLASSIFIED 


UNIV LAFAVETTE IN SCHOOL OF ELECTRICAL ENGINEERING 
A C KAK 21 DEC 84 DAHD17-82-C-2819 

F/G 28/1 














































n 


1 01 


3 0 COD- cm 



'5 5 T -|6 0 


o 



Figure 4.37 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.01, using the conventional technique. The attenuation constant is 3.0 
db/cm in both the object and the medium. The projections are sampled 
at 64 points over a sampling interval of 0.45 X. The receiver line is at a 
distance of 10 X from the cylinder, (a) Refractive index, (b) Difference in 
attenuation constant of the object and the medium. 

































< a 



<■ b ) 


Figure 4 38 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.03, using the conventional technique. Both the object and the medium 
are lossless. The projections are sampled at 64 points over a sampling in¬ 
terval of 0.45 X. The receiver line is at a distance of 10 X from the 
cylinder, (a) Refractive index, (b) Difference in attenuation constant of 
the object and the medium. 




















































I f ife 






15 5 




155 -ISO A 


y “CD 


stu 




- mi -^feg 








15 5^. lt o 


Figure 4.40 

Reconstruction of a cylinder with diameter of 8 X and refractive index of 
1.03, using the conventional technique. The attenuation constant is 3.0 
db/cm in both the object and the medium. The projections are sampled 
at 84 points over a sampling interval of 0.45 X. The receiver line is at a 
distance of 10 X from the cylinder, (a) Refractive index, (b) Difference in 
attenuation constant of the object and the medium. 

























0 1 (db/cm) 


0 2 (db/cm) 


(5»g! 


Uli 


C'li 




§>=£$: \v > , U J 

^ -JjJi 


T/TOfe 


; 'ft? 


fW.i 


nttiii 




Figure -4.41 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.0, using the conventional technique. The attenuation constant of the 
medium and the object are 0.1 db/cm and 0.2 db/cm, respectively. The 
projections are sampled at 64 points over a sampling interval of 0.45 X. 
The receiver line is at a distance of 10 X from the cylinder, (a) Refractive 
index, (b) Difference in attenuation constant of the object and the medi- 


p\ -". / A / 



























n “ 1 u ' = 2 0 (jb>cm) , t - 2 1 


(absent; 



-16 o m 



Figure 4.42 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.0, using the conventional technique. The attenuation constant of the 
medium and the object are 2.0 db/cm and 2.1 db/cm, respectively. The 
projections are sampled at 64 points over a sampling interval of 0.45 X. 
The receiver line is at a distance of 10 X from the cylinder, (a) Refractive 
index, (b) Difference in attenuation constant of the object and the medi- 


























1 0 


*J i <.db-CiH) 


1 0 \.db'Cn> 



< b ) 

Figure 4.43 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.0, using the conventional technique. The attenuation constant of the 
medium and the object are 0.1 db/cm and 1.0 db/cm, respectively. The 
projections are sampled at 64 points over a sampling interval of 0.45 X. 
The receiver line is at a distance of 10 X from the cylinder, (a) Refractive 
index, (b) Difference in attenuation constant of the object and the medi¬ 
um. 





















Figure I II 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.0, using the conventional technique. The attenuation constant of the 
medium and the object are 2.0 db/cm and 3.0 db/cm, respectively. The 
projections are sampled at 64 points over a sampling interval of 0.45 X. 
The receiver line is at a distance of 10 X from the cylinder, (a) Refractive 
index, (b) Difference in attenuation constant of the object and the medi¬ 
um. 






























constants of 0.5 and 3.0 db/cm, respectively. It is only for the ca^c of 3.0 
db/cm attenuation that the reconstruction is really distorted quantitatively, 
while for the value of 0.5 db/cm attenuation we get a reasonably accurate 
reconstruction. Since in most practical cases of interest the difference between 
the attenuation constants of the object and the medium is not larger than 10 
db/cm, perhaps the difference in the attenuation constants is not a limiting 
factor for the reconstruction algorithms in a lossy medium. 

In the final reconstruction (Fig. 4.17) we have tried to consider a case close 
to reality, taking all the parameters of the object and medium to be different. 
The refractive index of the object has been set to 1.02 and the difference in 
attenuation to 0.5 db/cm. The estimation error in the phase of the object 
function causes the estimate of the refractive index to become smaller than the 
actual value, while the difference in attenuation is overestimated. 

The simulation results presented in this section show that the only effect 
of the medium attenuation in the conventional approach is some loss of 
resolution. 

4.2.3 Attenuation in Synthetic Aperture Diffraction Tomography 

As discussed in chapter 10, one can not extend the Synthetic Aperture 
technique to the case of an attenuating medium. Therefore, in the absence of a 
reconstruction procedure specifically designed for a lossy medium we will 
examine the performance of the present algorithm in the presence of 
attenuation. As in the conventional tomography, the quantitative aspect of 
the reconstruction is of importance to us, thus we should compensate for the 
loss of energy in the propagation of the field to the receiver array. However, 
the compensation is quite different from that of the conventional technique. In 
the conventional technique, the attenuation is compensated for by multiplying 
the each projection data by a constant value, while in tne bynthetic Aperture 
technique the multiplication factor is dependent on the transducers positions. 
The scattered field at each position of the transmitter and receiver is amplified 
by the amount of attenuation the incident field has suffered traveling between 
the two transducers. This means that the field at the end of the arrays is 
weighted higher relative to the one at the renter. Of course, in the presence of 
noise, the noise will be amplified as well. In the following, we present 
reconstructions achieved in this way, based on the projection data calculated 
using the programs discussed in chapter 2. 


v 'ey- v rVvvv-.u J 


















n = 1 01 . € Q - 0 0 


= 0 5 C dt>- cm ) 






^ 5 -ib o 4 








-ISO VJ> 


Figure 4.45 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.01, using the conventional technique. The medium is lossless while the 
attenuation constant of the object is 0.5 db/cm. The projections are sam¬ 
pled at 64 points over a sampling interval of 0.45 X. The receiver line is at 
a distance of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium. 












































11 2 T -I* o 


'.5 5 


Figure 4 16 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.01, using the conventional technique. The medium is lossless while the 
attenuation constant of the object is 3.0 db/cm. The projections are sam¬ 
pled at 64 points over a sampling interval of 0 45 X. The receiver line is at 
a distance of 10 X from the cylinder, (a) Refractive index, (b) Difference 
in attenuation constant of the object and the medium. 






mJ?M. K A m. «. <"» a"*- Vk. .V. 


•*’» /*-. O. «*. **_ ft 


















17 <> 


i iu J 



to) 



15 5 




( b > 


Figure 4 17 

Reconstruction of a cylinder with diameter of 6 X and refractive index of 
1.02, using the conventional technique. The attenuation constants of the 
medium and the object are 3.0 db/cm and 3.5 db/cm, respectively. The 
projections are sampled at 64 points over a sampling interval of 0.45 X. 
The receiver line is at a distance of 10 X from the cylinder, (a) Refractive 
index, (b) Difference in attenuation constant of the object and the medi¬ 
um. 


‘ v.v- v.v vv V Vv'v l .- 



























In Fig. 4.48a,b we have shown the magnitude of the scattered field from a 
cylinder with a refractive index of 1.03 and a diameter of 6 wavelengths and its 
Fourier transform. The attenuation constants of the medium and the object 
were set to 3.0 db/cm. As one expects intuitively, the field goes to zero at the 
end of the receiver array and if no attenuation compensation is performed, the 
length of the arrays would not affect the reconstruction. From Fig. 4.48b, the 
attenuation of high frequency components is obvious. The reconstructed 
object’s Fourier transform is shown in Fig. 4.48c. It is quite clear that those 
values in the diagonal directions are attenuated and will cause some distortion 
in those directions. However, if the attenuation is compensated for, this effect 
will be reduced. In Fig. 4.49 we have shown the magnitude of the 
reconstructed a refractive index of a cylinder with refractive index of 1.03 and 
diameter of 6 wavelengths when the attenuation is compensated for. Tim 
attenuation constant of the medium and the object are both taken to be 3.0 
db/cm. In Fig. 4.50 we have shown the reconstruction when no attenuation 
compensation has been performed. Of course the reconstruction is qualitatively 
wrong and some distortion is present along the diagonal directions. When four 
views, 45 degrees apart, are used the distortion is removed as shown in Fig. 
4.51. In Figs. 4.50 and 4.51 the smoothness of the reconstruction is due to the 
loss of high frequency components in the projection due to the medium 
attenuation. 

In Fig. 4.52, the reconstructed difference between attenuation constant of 
the object and medium is presented. The attenuation constants of the medium 
and the object were set to 2.0 and 3.0 db/cm, respectively. The cylinder was 6 
wavelengths in diameter and the refractive index of the object and medium 
were taken to be equal. Other than a small amount of distortion in the 
diagonal directions, the quality of the reconstruction is satisfactory. 

Finally, in Fig. 4.53 we have shown the reconstruction results for the 
object considered in Fig. 4.47. The result obtained under the Synthetic 
Aperture technique (Fig. 4.53) is quite similar to the one obtained using the 
conventional technique (Fig. 4.47). In both of these reconstructions the 
attenuation has been compensated for. 




























r 


181 


n ■ 1 03 , J 0 ■ i . 3 0 (clb'cm) 



n * 1 03 . i ■ < . 3.0 (db'cnO 





< b > 


Figure 4.48 

The object is a lossy cylinder with diameter and refractive index of 8 X 
and 1.03, respectively. The attenuation constants of the medium and the 
object are both 3.0 db/cm. There are 128 elements in the transmitter and 
the receiver arrays, with the inter-element distance of 0.45 X. (a) Magni¬ 
tude of the scattered field, as a function of the transmitter and receiver 
positions, (b) Magnitude of the Fourier transform of the scattered field, 
with respect to the transmitter and receiver positions, (c) Magnitude of 
the reconstructed object’s Fourier transform. 































3.0 <db/cm) 
































e 


€ o ~ 3.0 (db/cm) 


1.60 

/\ 

3- 

I 



Figure 150 

Magnitude of the reconstructed refractive index of a cylinder with diame¬ 
ter of 8 X and refractive index of 1.03, using the Synthetic Aperture tech¬ 
nique. The attenuation constant of the object and the medium are both 
3.0 db/cm. The transmitter and the receiver arrays, positioned 20 X 
apart, nave 128 elements with the inter-element distance of 0.45 X. The 
attenuation of the medium has not been compensated for. 




















3.0 (db/cm) 


Figure 4.oI 

Magnitude of the reconstructed refractive index of a cylinder with diame¬ 
ter of 6 X and refractive index of i.03, using the Synthetic Aperture tech¬ 
nique. The attenuation constant of the object and the medium are both 
3.0 db/cm. The transmitter and the receiver arrays, positioned 20 X 
apart, nave 128 elements with the inter-element distance of 0.45 X. The 
attenuation of the medium has not been compensated for, but four rota¬ 
tional views of the object has been used. 




MB 




























e = 3.0 <db/cm> 


2.0 <db/cm) 


Figure 4.52 

Magnitude of the reconstructed attenuation constant of a cylinder with di¬ 
ameter of 6 X and refractive index of 1.0, using the Synthetic Aperture 
technique. The attenuation constants of the object and medium are 3.0 
db/cm and 2.0 db/cm, respectively. The transmitter and the receiver ar¬ 
rays, positioned 20 X apart, have 128 elements with the inter-element dis¬ 
tance of 0.45 X. The attenuation of the medium has been compensated 
for. 












































List of References 


[Mueller] R. K. Mueller and M. Kaveh, “Diffraction tomography,” Unpub¬ 
lished work. 

[Nahamoo82] D. Nahamoo and A. C. Kak, "Ultrasonic diffraction imaging,” 
Purdue University, School of Electrical Engineering, Technical 
Report, TR-EE-82-20, 1982. 

[Nahamoo84] D. Nahamoo, S. X. Pan, and A. C. Kak, “Synthetic aperture 
diffraction tomography and its interpolation-free computer im¬ 
plementation,” IEEE Transactions on Sonics and Ultrasonics, 
Vol. SU-31, No. 4, pp. 218-229, July 1984. 


I 






























CHAPTER 5 
SUMMARY 


To give the reader a brief overview of the work presented in the second 
part of this thesis, we have devoted this chapter to the summary of the topics 
and results already discussed in chapters 2, 3, and 4. 

The forward problem in diffraction imaging, by which we mean calculation 
cf mattered fields from the inhomogeneities present in an object, was addressed 
in chapter 2. Two distinct types of objects were considered: single-component 
objects and multi-component objects. The object components were assumed to 
be cylinders, because in the formulation of the problem for non-cylindrical 
objects one is faced with severe mathematical difficulties. 

For incident plane waves, we used the well-known mathematical 
expressions for the scattered fields derived in most of the classical books on 
scattering theory. For the case of an incident field generated by a point source, 
we derived similar expressions for the exact scattered fields using the series 
expansion of a Hankel function in cylindrical coordinates. A computer program 
was developed to simulate the process by utilizing the mathematical 
formulations. Some parts of the code was vectorized to speed up the floating 
point calculations in the array processor FPS-AP120B. By using the developed 
computer program, we are now capable of generating the exact scattered fields 
from a cylinder with arbitrary size, refractive index, and attenuation constant 
embedded in a medium with an arbitrary attenuation constant, illuminated by 
either a plane wave or a point source. 

For calculating the scattered fields from multi-component objects, a major 
source of difficulty is dealing with the interaction between the various 
components. Depending on the interaction between the components, the total 
scattered field may or may not bear any resemblance to the simple sum of the 
scattered fields for each of the components, assuming the others to be absent. 
In chapter 2, we presented a new computational procedure for calculating the 
inter-component interaction. We have based our technique on a basic theory in 
which the scattered field is expressed as an infinite summation of different 



190 


order terms. For the first order scattering, each component interacts 
independently with the incident illumination, being oblivious of the existence of 
the others. Second order scattering terms are generated when each component 
interacts with the first order fields sent in its direction by the other 
components, and so on. Second order fields are calculated by expanding the 
first order fields from each component as superposition of plane waves and then 
analyzing the interaction of each of these plane waves with the other object 
components. Higher order fields can be generated similarly. 

This technique was coded in a vectorized fashion to benefit from the speed 
of available array processor. Although our program only calculates first and 
second order scattering terms, extension of the code for calculation of higher 
order terms is trivial. Using this procedure we showed by computer 
simulations that even when each component of a two-component object is 
weakly scattering, the multiple scattering effects become important when the 
components are blocking each other. We further showed that when strongly 
scattering components that are large compared to a component are not 
blocking each other, the scattering effects can be ignored. Both these 
conclusions agree with intuitive reasoning. 

In chapter 3. we investigated the inverse scattering problem. In general 
inverse scattering is a much more complicated problem than the forward 
problem. Inverse scattering techniques are usually accompanied by severe 
restrictions on the different parameters of the system under consideration. We 
have considered those inverse scattering techniques which are based on the 
Fourier Diffraction Projection Theorem. This theorem related the Fourier 
transform of an object with the Fourier transform of the scattered field 
measured on a straight receiver line when the object is illuminated by a plane 
wave. It should be noted that the Fourier Diffraction Projection Theorem is 
valid only when the object inhomogeneities are weakly scattering, or in fact 
when either the two well known Born or Rytov approximations hold. 

The two inverse reconstruction techniques considered are: the conventional 
technique and the Synthetic Aperture technique. In the conventional technique 
the object is illuminated by a plane wave and the scattered field is measured on 
a receiver line which is perpendicular to the direction of the propagation of the 
incident field. Many rotational views of the object provides the necessary data 
to fill up the Fourier space of the object inside a circle of radius v^kr, (k 0 is the 
wavenumber). The object's cross section is then recovered by an inverse 
Fourier transformation 















The Synthetic Aperture technique is a more general technique which 
permits the use of any type of insonification, as opposed to only plane waves 
and also requires only two rotational view's of the object. However, for each 
rotational position, the projection data is collected by independent traversal 
movements of the transmitting and the receiving transducers, over two parallel 
lines on both sides of the object. In this technique, the coverage in the Fourier 
domain of the object is inside a circle of radius v/2k 0 , the same as that of the 
conventional approach. 

We have modified the theory behind the conventional reconstruction 
technique to include the case of lossy objects and media as well. We have 
shown that in a lossy medium one would recover the two dimensional Laplace 
transform of the object over a surface and moreover the inverse Laplace 
transform techniques can not be used to recover the object’s cross section. 
Using the Synthetic Aperture technique in a lossy medium, the object's 
Fourier-Laplace transform is recovered over a surface in the Fourier-Laplace 
domain of the object. By Fourier-Laplace transformation we mean Fourier 
transformation in one direction and Laplace transformation in the other 
direction. Still no inverse transform technique can be used for the inversion 
process. A possible way to avoid the inverse Fourier-Laplace transformation is 
to impose a constraint on the processed data to convert the transform to a 
simple two-dimensional Fourier transform. However, the experimental 
procedure has to be repeated for many different rotational views of the object. 
Of course this requires enormous amount of data to be collected and may be 
quite restrictive in some applications. 

In chapter 3, it has also been proven that if the Synthetic Aperture 
technique is used in the reflection mode as well as the transmission mode, one 
rotational position of the object is enough for recovering the object’s Fourier 
transform. In addition, the coverage would be inside a circle of radius 2k 0 
rather than v/2k 0 for the conventional and the Synthetic Aperture techniques. 
Further study of this technique still remains to be examined by simulation.il 
means. 

Chapter 1 is devoted to a simulation study of multiple scattering and 
attenuation phenomena in diffraction imaging. The programs developed in 
chapter 2 Ur the calculation of scattered fields together with the reconstruction 
schemes presented in chapter 3 were used to generate the computer simulation 
results in this chapter. 




All diffraction tomography algorithms available today, including the 
conventional and the Synthetic Aperture techniques, assume the first order 
Bom or Rytov approximations; which implies that they are only good for 
reconstructing weakly scattering objects. In a multi-component object this 
translates into ignoring the interaction between the components. However in 
most cases of practical interest the object would be composed of several 
distinct components. One of our objectives in chapter 4 has been to study the 
effect of the interaction between object components on the reconstructed 
object, in both the conventional and the Synthetic Aperture techniques. 

We generated two series of reconstructions of a multi-cylinder object using 
the conventional reconstruction technique. In the first series, the projections 
were generated by ignoring the second and higher order scattered fields. In the 
other the second order scattered fields were included in the projections. A 
comparison of these two series of simulations showed that the scattering 
interaction between the components cause some distortion along the line 
joining the components, the distortion increasing with larger values of the 
refractive indices of the components. The distortion also becomes more 
noticeable as the number of components increases. 

We examined three different approaches for reducing the distortions in the 
conventional technique. First we made the receiver array to be oblique to the 
direction of the incident field. In this way the blocking effect, which is quite 
severe when the propagation vector and the line joining the object components 
are collincar, is reduced. Simulation results showed that although this 
arrangement reduces the distortion in the region between the object 
components, the price paid is an increase in the general distortion all through 
the object. This approach seems appropriate for two-component objects. 

In the second approach intended for more complex objects, the object is 
illuminated with several non-collinear plane waves for each projection, for each 
rotational position of the receive array. Fach plane wave generates a complete 
frequency domain coverage, these are averaged prior to Fourier inversion for 
object reconstruction. The rationale for this method is that the error caused 
by blocking vary for the three different illuminating waves. The simulation 
results were similar to that of the first approach. 

We also considered another approach which is similar in spirit to the 
Synthetic Aperture technique For each position of the receiver line the object 
is illuminated by incident plane waves from several directions. However, unlike 
the second approach, no averaging is performed in the frequency domain. The 
object is now rotated by a large angle to yield a similar coverage over areas left 








193 



blank before; this procedure is repeated until no holes are left in the frequ 
domain. The rationale here is that with each plane wave at a difff 
orientation with respect to the receiver line, the severity of the blocking e 
caused by multiple scattering should be reduced. Simulation results show 
this approach does result in better reconstructions of multi-component obje< 

The superiority of the third approach to the others persuaded us that 
performance of the Synthetic Aperture technique in the presence of m 
component interactions should be superior to that of the conventh 
technique. Simulation results proved that our expectation was in fact t 
There was almost no distortion between the components as was the case in 
reconstructions achieved with the conventional technique. The reason is t 
the blocking effect is minimized in this technique. Another interesting p( 
was that the components which were blocked became more distorted compa 
to the ones at closer distances to the transmitter array. Our general conclus 
was that the performance of the Synthetic Aperture technique in the prese 
of multiple scattering is far superior than that of the conventional technique. 

In the second half of the chapter 4, we discussed the effect of 
attenuation on the object reconstructions in diffraction tomography. Althoi 
in ultrasonic imaging the role of attenuation is minor, in microwave imagin' 
can not be ignored, because microwaves at 1 GHz undergo 3 db of attenuat 
per centimeter in water. In the lack of any reconstruction technique specifi 
designed for lossy media, we studied the effect of medium attenuation on 
calculated projections and thereby on the reconstructions themselves. 

We proved that the medium attenuation suppresses the high frequer 
components present in the projections. In the conventional technique this le: 
to a loss of resolution in the final object's reconstruction. This point is back 
by our computer simulation results. However in the Synthetic Aperti 
technique, medium attenuation causes a loss of coverage in the object’s Four 
domain in a manner that the lost region can be recovered only by consider; 
more rotational views of the object. Heuristic compensation of medii 
attenuation in both techniques results is acceptable quantitative obj 
reconstructions. However such compensation techniques may not be advisal 
in the presence of noise in the collected projections. 

Finally simultaneous reconstruction of the refractive index and t 
attenuation constant of an object was discussed in chapter 4. It was sho 1 
that large refractive index of an object could result in an erroner 
reconstruction of the attenuation constant, while typical small differem 
between the attenuation constants of the medium and the object would r 


194 


affect the reconstructed refractive index of the object. As long as the noise level 
is under certain threshold so that the heuristic compensations for the medium 
attenuation can be applied, the attenuation of the medium and object will not 
be a limiting factor. 





















BIBLIOGRAPHY 


[Andersen82j 

A. H. Andersen and A. C. Kak, “Digital ray tracing in two-dimensional re¬ 
fractive fields,” J. Acoust. Soc. Amtr., Vol. 72, pp. 1593-1606, 1982. 

[Andersen84] 

A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction 
technique (SART): A superior implementation of the ART algorithm,” Ul¬ 
trasonic Imaging, Vol. 6, pp. 81-94, Jan. 1984. 

[Azimi83j 

M. Azimi and A. C. Kak, “Distortion in diffraction tomography caused by 
multiple scattering,” IEEE Transactions on Medical Imaging, Vol. MI-2, 
No. 4, pp. 176-195, Dec. 1983. 

[Ball80] 

J. J. Ball and F. Stenger, “Explicit inversion of the Helmholtz equation for 
ultrasound insonification and spherical detection,” Acoustical Holography, 
Vol. 9, A. Metherell, Ed., 1980. 

[Bellman66| 

R. Beelman, R. E. Kalaba, and J. A. Lockett, “Numerical Inversion of 
the Laplace Transform: Applications to Biology, Economics, Engineering, 
and Physics,” in Modern Analytic and Computational Methods in Science 
and Mathematics, R. E. Bellman and R. E. Kalaba, Eds., American El¬ 
sevier Publishing Company, Inc., 1966. 

[Bowman58j 

F. Bowman, Introduction to Bessel Functions, Dover Publications Inc.. 
New York, 1958. 

[Carter70] 

W. II. Carter, ‘Computational reconstruction of scattering objects from 
holograms," Journal of the Optical Society of America, Vol. 60, No. 3. pp. 
306-314, March 1970. 








■ 

























196 


[Crawford82] 

C. R. Crawford and A. C. Kak, “Multipath artifact corrections in ultra¬ 
sonic transmission tomography,” Ultrasonic Imaging, Vol. 4, pp. 234-266, 
1982. 

[Dandliker70j 

R. Danaliker and K. Weiss, “Reconstruction of the three-dimensional re¬ 
fractive index from scattered waves,” Optics Communications, Vol. 1, No. 
7, Feb. 1970. 

[Devaney82] 

A. J. Devaney, “A filtered backpropagation algorithm for diffraction to¬ 
mography,” Ultrasonic Imaging, Vol. 4, pp. 336-350, 1982. 

[Devaney83l 

A. J. Devaney, “A computer simulation study of diffraction tomography,” 
IEEE Transactions on Biomedical Engineering Vol. BME-30, No. 7, pp. 
377-386, July 1983. 


[Devaney84l 

A. J. Devaney, “Geophysical diffraction tomography,” IEEE Transactions 
on Geoscience and Remote Sensing ,” Vol. GE-22, No. 1, pp. 3-13, Jan. 
1984. 

[Fox 68] 

L. Fox and D. F. Mayers, Computing Methods For Scientists and En¬ 
gineers, Chap. 3, Clarendon Press, Oxford, 1968. 


(Goldstein59l 

M. Goldstein and R. M. Thaler, “Recurrence techniques for the calculation 
of Bessel functions,” Math. Tables and Aids to Comp. (Now Mathematics of 
Computation), Vol. 13, pp. 102-108, 1959. 


(Goodman68l 

J. W. Goodman, Introduction to Fourier optics, McGraw-Hill Book Co., 
1968. 


[Hochstadt73] 

H. Hochstadt, Integral Equations, John Wiley & Sons Inc., 1973. 


[Iwata75] 

K. Iwata and R. Nagata, “Calculation of refractive index distribution 
from interforograms using Born and Rytov’s approximation,” Japan J. 
Appl. Phys., Vol. 14, Supp. 14-1, pp. 383, 1975. 






[JohnsonSo] 

S. A. Johnson and M. L. Tracy, “Inverse scattering solutions by a Sine 
basis, multiple source, moment method - Part I: Theory,” Ultrasonic Imaa- 
ing, Vol. 5, pp. 361-375, 1983. 

[Johnson84l 

S. A. Johnson, Y. Zhou, M. K. Tracy, M. J. Berggren, and F. Stenger, “In¬ 
verse scattering solutions by a Sine basis, multiple source, moment method 
- Part III: Fast algorithms,” Ultrasonic Imaging , Vol. 6, pp. 103-116, 1984. 

[Kaveh79a] 

M. Kaveh, R. K. Mueller, R. Rylander, T. R. Coulter and M. Soumekh, 
“Experimental results in ultrasonic diffraction tomography,” Acoustical 
Imag-ng, Vol. 9, pp. 433-450, 1979. 

[Kaveh79b] 

M. Kaveh, R. K. Mueller, and R. D. Inverson, “Ultrasonic tomography 
based on perturbation solutions of the wave equation,” Computer Graphics 
and Image Processing, Vol. 9, pp. 105-116, 1979. 

[Keller69] 

J. B. Keller, “Accuracy and validity of the Born and Rytov Approxima¬ 
tions,” J. Opt. Soc. Am., Vol. 59, pp. 1003-1004, Aug. 1969. 

[Kenue82] 

S. K. Kenue and J. F. Greenleaf, “Limited angle multifrequency 
diffraction tomography,” IEEE Transactions on Sonics and Ultrasonics, 
Vol. SU-29, No. 6, pp. 213-217, July 1982. 

[Krylov 69] 

V. I. Krylov and N. S. Skoblya, Handbook of Numerical Inversion of La¬ 
place Transforms, Jerusalem: Israel Program for Scientific Translations, 
1969. 

[KryIov77j 

V. 1. Krylov and N. S. Skoblya, A Handbook of Methods of Approximate 
Fourier Transformation and Inversion of the Laplace Transformation, Mos¬ 
cow: MIR Publishers, 1977. 

[MarionOl] 

G. C. Marion, “Bessel functions of integral order and complex arguments,” 
Communications of the ACM, Vol. 4, No. 4, page 169, April 1961. 

























198 


[Morse53] 

P. M. Morse and H. Feshbach, Methods of Theoretical Physics, Part I, 
Chap 7, McGraw-Hill Book Co., 1953. 


[Morse68j 

P. M. Morse and K. V. Ingard, Theoretical Acoustics, Chap. 8, McGraw- 
Hill Book Co., 1968. 


[Mueller78] 

R. K. Mueller, M. Kaveh and R. D. Iverson, “A new approach to acousti¬ 
cal tomography using diffraction techniques,” Acoustical Holography, Vol. 
8, A. Metherell, Ed. , New York: Plenum Press, 1978. 


[Mueller79j 

R. K. Mueller, M. Kaveh and G. Wade, “Reconstruction tomography and 
applications to ultrasonics,” Proe. IEEE, Vol. 67, pp. 567-587, 1979. 

[Mueller80] 

R. K. Mueller and M. Kaveh, “Diffraction tomography,” Unpublished 
work. 


[Nahamoo82] 

D. Nahamoo and A. C. Kak, “Ultrasonic diffraction imaging,” Purdue 
University, School of Electrical Engineering, Technical Report, TR-EE- 
82-20, 1982. 

[Nahamoo84] 

D. Nahamoo, S. X. Pan, and A. C. Kak, “Synthetic aperture diffraction 
tomography and its interpolation-free computer implementation,” IEEE 
Transactions on Sonics and Ultrasonics, Vol. SU-31, No. 4, pp. 218-229, 
July 1984. 


[Norton8l] 

S. J. Norton and M. Linzer, “Ultrasonic reflectivity imaging in three di¬ 
mensions: Exact inverse scattering solutions for plane, cylindrical and 
spherical apertures,” IEEE Transactions on Biomedical Engineering, Vol. 
BME-28, No. 2, pp. 202-220, 1981. 


[Pan83] 

S. X. Pan and A. C. Kak, “A computational study of reconstruction algo¬ 
rithms for diffraction tomography: Interpolation vs. filtered- 
backpropagation,” IEEE Transactions on Acoustics, Speech and Signal 
Processing, Vol. ASSP-31, No. 5, pp. 1262-1275, October 1983. 


















m 


[Rosen feld82] 

A. Rosenfeld and A. C. Kak, Digital Image Processing, Academic Press, 
Second Edition, 1982. 

[SancerTO] 

M. I. Sancer and A. I). Varvatsis, “A comparison of the Born and Rytov 
methods,” Proc. of the IEEE, pp. 140-141, Jan. 1970. 


[SlanevSl] 

Malcolm Slaney, A. C. Kak, and L. E. Larsen, “Limitations of imaging 
with first-order diffraction tomography,” IEEE Transactions on Microwave 
Theory and Techniques, Vol. MTT-32, No. 8, pp 860-874, August 1984. 


[Soumekh83] 

M. Soumekh, M. Kaveh, and R. K. Mueller, “Algorithms and Experimen¬ 
tal Results in Acoustical Tomography Using Rytov’s Approximation,” 
Proceeding of International Conference on Acoustics, Speech and Signal 
Processing, Vol. 1, pp. 135.138, 1983. 

[Stegun57] 

I. A. Stegun and M. Abramowitz, “Generation of Bessel functions on high 
speed computers,” Math. Tables and Aids to Comp. (Now Mathematics of 
Computation), Vol. 11, pp. 255-257, 1957. 


[Tracy 83j 

M. K. Tracy and S. A. Johnson, “Inverse scattering solutions by a Sine 
basis, multiple source, moment method - Part II: Numerical evaluations,” 
Ultrasonic Imaging, Vol. 5. pp. 376-392, 1983. 


[Weeks64] 

\V. L. Weeks, Electromagnetic Theory for Engineering Applications, John 
Wiley & Sons, Inc., 1964, ch. 6. 


[Wolf69] 

E. Wolf “Three dimensional structure determination of semi-transparent 
objects from holographic data,” Opt. Comm., Vol. 1, pp. 153-156, 1969. 








200 


APPENDIX 


The aim in this appendix is to solve the following set of equations for u ; in 
terms of u r , v, k 0 , and ( Q 

(li+k,) 2 + (v+k,) 2 =k J , (1) 

k„ 2 + k, 2 = k to 2 , (2) 

<1 - u, + j Uj (3) 

k„ = k „ + j <„ ■ (■») 

Substituting Eq. (2) in (1), we have: 

u 2 + v 2 + 2vk y = ~2uk x . (5) 

Squaring both sides of Eq. (5) and making use of Eq. (2), we get 

u 4 + (v 2 + 2vk y ) 2 + 2u 2 (v 2 + 2vk y ) - 4u 2 (k co 2 -ksyby 2 ) =0 (6) 

Real and imaginary part of Eq. (6) define the following two equations: 

(u r 2 -u 2 ) 2 - 4u r 2 u 2 + (v 2 + 2vk y ) 2 + 2(u r 2 ~u 2 ) (v 2 + 2vk y ) 

- 4(u r 2 -u ; 2 ) (k 0 2 -t 0 2 ~k y 2 ) + 16u r Ujk 0 t 0 = 0 (7) 

u r u ,K 2 " u i 2 ) + u r Ui(v 2 + 2vk y ) - 2(u r 2 -Uj 2 )k 0 e 0 

~2u r u i (jc 0 2 -c 0 2 -k y 2 ) = 0 (8) 

Substitution of the expression for k y from Eq. (8) into Eq. (7), provides us with 

the following fifth order equation for U; in terms of u r , v, k c , and f c . 

76 u i 5 + lb"* + I 4 "? + 73 u i 2 + TjV + 7i = 0 (9) 

where 

76 " . (Wa) 





















72 = -u r (u r 4 + v 4 - 4(k 0 2 -« o 2 v 2 + 2u r V ) 


(lOe 


*71 = 4u r 2 k 0 f 0 (u r 2 + v 2 ) (10f 

t 

\ second order equation for k y is achieved by rearranging the terms in Eq 


k y 2 + vk y + 7(Ur 2 -Ui 2 + V 2 ) - (k 0 2 ~f 0 2 ) - - ■ -“'- koto = 0 • 

i UfU; 


































Y> 



/ s m