# 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