Acoustic imaging device

Abstract

This invention provides an acoustic imaging device which utilizes scattered acoustic wave information to produce an image. The system consists of an array of acoustic transducers which at least partially encircle a medium to be imaged. The medium is sequentially insonified by each transducer with an acoustic signal comprised of at least one discrete frequency. The acoustic signal scatters from the medium and is detected by the remaining transducers. A data set representing the phase and amplitude of each discrete frequency at each transducer is acquired and used by a computer to calculate an image of the medium. Images acquired and calculated at each identifiable discrete frequency are combined in various ways in order to improve the quality of the final image of the medium.

Description

BACKGROUND OF THE INVENTION

In conventional fan beam transmission tomography, such as utilized in x-ray computed tomography (x-ray CT), the attenuation of a fan beam of x-rays is measured as this beam probes medium along many different trajectories. The information contained in these attenuation projections is then used to reconstruct a tomographic image of the medium. The success of x-ray CT (manifested in the resolution and clarity of the images) is fundamentally linked to the very short wavelength of the incident x-ray beam (.apprxeq.1 .ANG.). The mathematical reconstruction techniques assume a geometrical ray approximation and are quite simplified.

In ultrasonic tomography, the acoustic wavelengths (.apprxeq.1 mm) are on the order of the scale size of the inhomogeneities of the medium and the diffraction effects of the acoustic wave are not negligible. In this case, the attenuation of a sound wave is substantially affected by scattering effects such as diffraction, refraction, and reflection, as well as absorption. The simplified mathematical reconstruction algorithms used in conventional x-ray CT are not applicable in this case.

Ultrasonic diffraction tomography (UDT) techniques attempt to mathematically reconstruct a tomographic image of a medium from acoustic wave data with full consideration to the scattering effects associated with acoustic wavelengths involved. This is done by considering the full wave equation, a much more difficult problem to develop and implement than the geometrical wave approximation reconstruction methods used in x-ray CT. UDT techniques attempt to determine the internal structure of an object which is semi-transparent to acoustic waves from a set of scattered wave data detected outside at the boundary of the object. A number of these mathematical techniques have been theoretically considered, including simplified algorithms which use the Born approximation (E. Wolf, "Three-dimensional Structure Determination of Semi-Transparent Objects from Holographic Data", Opt. Comm. 153-156 (1969)) and Rytov approximations (A. J. Devaney, "Inverse Scattering Within the Rytov Approximation", Opt. Lett. 6, 374 (1981)). The Born and Rytov approximations assume that the medium which is imaged is a weak scatterer of acoustic waves and that the acoustic wave does not experience large phase shifts within the medium. Computer intensive full-wave reconstruction algorithms (S. Johnson and M. Tracy, "Inverse Scattering Solutions By a Sinc Basis, Multiple Source, Moment Method", Ultrasonic Imaging 5, 361-375 (1983) and references therein) have also been proposed for imaging mediums under conditions which are not in the range of validity of the Born or Rytov approximation.

Prior art patents included Devaney (U.S. Pat. Nos. 4,598,366 and 4,594,662) and Johnson (U.S. Patent No. 4,662,222). The patents of Johnson disclose iterative algorithms which form an initial estimate of the sound speed at all points within the object being imaged, calculate a sound speed map (a type of image), and then update the map. This process is repeated until a residual error parameter is small enough. The patents of Devaney form an image directly from the data using Born and Rytov inversions with a technique called filtered back propagation.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide a system and method to insonify a medium with ultrasonic waves, record information relating to ultrasonic waves scattered from the medium and use that information to quickly and efficiently construct images of the medium.

This invention provides an acoustic imaging device which utilizes scattered acoustic wave information to produce an image. The system consists of an array of acoustic transducers which at least partially encircle a medium to be imaged. A preferred embodiment of the device consists of a ring of acoustic transducers which encircle a medium to be imaged. The medium is sequentially insonified by each transducer with an acoustic signal comprised of at least one identifiable discrete frequency. The acoustic signal scatters from the medium and is detected by the remaining transducers. A data set representing the phase and amplitude of each identifiable discrete frequency at each transducer is acquired and is used by a computer to calculate an image of the medium. Images acquired and calculated at each indentifiable discrete frequency are combined in various ways in order to improve the quality of the final image of the medium.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram for use in illustrating the mathematical reconstruction algorithm.

FIG. 2 is a block diagram of the electrical circuitry used in an actual prototype device of the invention.

FIG. 3 is a diagram of the computer simulated object used to demonstrate some of the features of the present invention.

FIG. 4 are the real (right) and imaginary (left) components of a single frequency image of the simulated object.

FIG. 5 are sound speed c(r,.theta.) (left) and attenuation .mu.(r,.theta.) (right) maps of the simulated object.

FIG. 6 are the real (right) and imaginary (left) components of the simulated object. The images are the summation of single frequency images as described in FIG. 5 which have been phase aberration corrected.

FIG. 7 is a composite image of a simulated object. This composite image combines the sound speed map described in FIG. 6 with real component of the image described in FIG. 7.

FIG. 8 displays nine images of a female breast produced with the preferred algorithms using the prototype device described in FIGS. 1-3.

FIG. 9 displays nine images of a porcine kidney produced with the preferred algorithms using the prototype device described in FIGS. 1--3.

FIG. 10 displays images of nine slices of a normal female breast produced with the prototype device described in FIGS. 1--3.

FIG. 11 displays images of nine slices of a cystic female breast produced with the prototype device described in FIGS. 1-3.

FIG. 12 displays an image of a porcine kidney produced with the prototype device described in FIGS. 1-3.

FIG. 13 is an illustration of the invention utilized for breast imaging showing how slices at various locations are obtained.

FIG. 14 is an illustration of the invention utilized for abdomen imaging.

FIG. 15 is an illustration of an application of the invention utilized for biopsy needle guidance.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

This invention, provides a system for insonifying a medium with acoustic waves at discrete frequencies from a plurality of locations, and for recording data representing acoustic waves scattered from the medium. Algorithms are provided for calculating images of the medium from the recorded data.

A. Description of Instrumentation

A schematic diagram of the cross-section of a portion of a preferred embodiment and an object being imaged is shown in FIG. 1. In this embodiment, there are 1024 acoustic transducers 10 or 10.1, . . . , 10.1024 evenly arranged on the locus of a circle of radius r.sub.o =102 mm. The transducers 10, each which can act as either a transmitter or a receiver of acoustic waves, are located in a thin slice approximating a plane which we call the (r,.theta.) (or (x,y)) plane. The z-direction is perpendicular to this plane. The transducers circumscribe object 9 which is assumed to be semi-transparent to the acoustic waves. In the preferred embodiment, a coupling fluid 11 resides in the remaining volume (shown as an area) inside the circle to effectively couple acoustic waves between the transducers 10 and object 9. The object 9 and the coupling fluid 11 in the remaining volume is herein referred to as medium 12. Water is a preferred fluid for imaging human tissue because water and human tissue have approximately equal sound speed and density. When imaging other objects, a coupling fluid 11 should be chosen which preferably has a density .rho..sub.o and sound speed c.sub.o which are approximately equal to the average density and sound speed of the object 9 being imaged (salt water is sometimes better than pure water for use with human tissue).

In a preferred embodiment of the present invention, each acoustic transducer 10 has a 1 MHz center frequency, a 60% bandwidth, and each performs at different times as a receiver or a transmitter of acoustic waves. At a 1 MHz frequency, the acoustic wavelength .lambda..sub.o =1.5 mm in water and approximately this wavelength in human tissue. The spacing between transducers is 0.65 mm which is somewhat less than .lambda..sub.o /2 in water in order to provide adequate sampling of the scattered wavefront (i.e. to avoid aliasing)., This facilitates the construction of unique images from the acquired data. In this preferred embodiment, each transducer has a dipole radiation pattern in the (r,.theta.) plane. The vertical height of the,, transducers 10 is 12 mm and each transducer produces a beam which is approximately collimated in the z-direction for distance in the z direction of approximately 12 mm and propagates from the transducer in the (r,.theta.) plane. We thus insonify an effective thickness across the medium roughly equal to the vertical height of the transducers so that any image obtained could be expected to represent an average of scatter information over the height of the transducers (i.e. a slice about 12 mm thick).

A schematic diagram of our existing laboratory system is displayed in FIG. 2. The transducer ring 10 consists of 1024 acoustic transducers 10.1, . . . , 10.1024 arranged in a circle having a radius of 102 mm. The transducer ring 10 is divided into eight octants designated as octants 1 through 8, each octant containing 128 transducers. Connected to the ring are eight multiplexer sets 70 each multiplexer set consisting of a two-channel 64-to-1 receiver multiplexers 72, two amplifiers 74 each amplifier connected to one channel of the two-channel multiplexer, a two-channel 64-pin transmit multiplexer 76, two drivers 78 each driving one channel of the two channel transmit multiplexers, and an address decoder unit 79. The data is collected by eight two-channel analog-to-digital converter units 22 each having 16 megabytes of buffer memory. The system timing is accomplished by a programmable master clock 28 which drives programmable timing and control unit 42 and a phase-locked arbitrary wave form generator 26. The entire system is controlled by a supervisor computer unit 30 the supervisor having associated with it a high resolution monitor 32, a keyboard 34, a 3 gigabyte hard drive 38, and an ethernet link 31. Timing commands are communicated to the transducer multiplexer sets 70 through a command interpreter 50. The transmitter signal is generated by the arbitrary waveform generator 26 through a transmit buffer 52 and to the transducer multiplexers 70.

The individual components of the described preferred system are described in detail here:

1) Transducers

The transducers 10 used in our system are constructed from composites of piezoelectric and non-piezoelectric materials. Each acoustic transducer has a 1 MHz center frequency, a 60% bandwidth, and each performs at different times as a receiver or a transmitter of acoustic waves. At a 1 MHz frequency, the acoustic wavelength .lambda..sub.o =1.5 mm in water and approximately this wavelength in human tissue. The spacing between centedines of adjacent transducers is somewhat less than .lambda..sub.o /2 in water in order to provide adequate sampling of the scattered wavefront (i.e. to avoid aliasing). In addition, because the width of each transducer is less than .lambda..sub.o /2, each transducer has a dipole attenna pattern in the (r,.theta.) plane. This means that the relative amplitude of the transmitted acoustic energy follows a cos .psi. function of the angle .psi. measured relative to the vector normal to the surface of the transducer face. There are two acoustic impedance matching layers 17 between the transducers 10 and the coupling fluid 11. These layers help to couple acoustic energy from the transducers 10 to the coupling fluid 11. In addition, these layers 17 help to minimize the amount of acoustic energy reflected from the faces of the transducers 10 back into the medium 12. This reflected acoustic energy manifests itself as image artifacts in the final images of the medium 12. Several variations in transducer fabrication could be utilized instead of composites of piezoelectric and non-piezoelectric materials. These include transducers constructed from piezoelectric zirconate titanate (PZT5A, for example), piezoelectric film (PVDF), or co-polymers.

2) Transmitter and Receiver Multiplexers

Transmitter and receiver multiplexer and units are preferably fabricated from standard off-the-self multiplex equipment commonly used in communication equipment. We use Analog Device Inc. ADG506 devices and Com-Lines CLC522 devices.

3) Address Decoder

Our address decoder 79 is a model ED910 programmable logic device supplied by Altera Corporation. This unit is programmed to operate multiplexers 72 and 76.

4) Command Interpreter

Command interpreter 50 receives commands from programmable timing and control module 42 and transmits those commands to the address decoders 79 in multiplexer units 70 in order to activate both the transmitting and receiving transducers. For the command interpreter 50, we use type EP1810 programmable logic devices supplied by Altera Corp.

5) Programmable Master Clock

The programmable master clock 28 is a computer programmable generator capable of precisely generating a known clock frequency for the system. This component 28 generates the master time base that phase locks every element of the system and is set when the system is initialized for the data collection algorithm in use. The programmable master clock 28 is a VXI module of model number 1391 supplied by WaveTek Corporation.

8) Programmable Timing and Control Module

The programmable timing and control module 42 is built around an Advance Micro Devices Corporation 1DT 2910 micro sequencer. The module consists of a programmable bit slice micro-sequencer and micro instruction RAM that is clocked by the programmable master clock 28. This device is designed to generate a unique bit pattern on each successive clock cycle. It is used to initialize the multiplexers, control the command interpreter, produce sample clocks to the analog-to-digital converters 22, trigger the arbitrary wave form generator and enable data collection. Since the module is locked to the time base generator, all of these functions are phase locked operations. The module is designed to interface with the VME bus chassis of supervisor 30 as described below.

6) Arbitrary Waveform Generator

Our arbitrary waveform generator 26 is a phase locked computer programmable generator capable of producing the various arbitrary analog output waveforms necessary for various data acquisition methods. It is used to synthesize the waveform which is transmitted through the transducers 10 into the medium 12. Since it is phase locked to the programmable master clock 28 and triggered by the programmable timing and control module 42, the transmitted wave form is phase locked to the data collection sequence. Our arbitrary waveform generator 26 is a VXI module of model number 1395 supplied by Wavetek Corporation.

7) Analog-to-Digital Converters

The design of the analog-to-digital converters 22 is centered around the model SPT 7922 12 bit, 30 MHz analog-to-digital converters supplied by Signal Processing Technologies. Each channel contains 8 megabytes of local buffer memory. The converters 22 are programmed to collect the analog data received via receive multiplexers 72 and temporarily store this data in the digital buffer memory. The preferred system has 16 channels of analog-to-digital conversion which has a capacity to store 128 megabytes of data.

9) Computer Supervisor

The overall control of the imaging device is provided by computer supervisor 30. Supervisor 30 is used to initialize the system, compile and download timing and control algorithms, perform test functions, collect raw data and perform the data reduction algorithm producing the final image. Interface to the machine is provided by a keyboard 34 and monitor 32. A 3 gigabyte hard drive 38 is provided for data storage. An ethernet link 31 is provided from the machine to download or upload data and algorithms for subsequent off line processing and reduction. The supervisor 30 is a series of board level products from Radisys Corporation based on the Intel 80486 computer family and a compiler supplied by Zortech Corporation.

10) Important Elements

Images obtained with prototype devices built by the inventors and their co-workers are reproduced in FIGS. 8 through 11. These images and the specific algorithms we used to obtain the images are discussed in detail in subsequent sections of this patent. Applicants believe that their images are the best, by far, tomographic images ever obtained with an acoustic devices.

The principle advantages of the present system that has led to these improvements over the prior art are the combination of the following hardware elements:

1. A large number of transducers spaced generally less than .lambda./2 apart on a circle surrounding the object where .lambda., is the acoustic wavelength in the medium 12. 2. A choice of acoustic frequencies between 100 kHz and 1.5 MHz which allows the acoustic wave to traverse human tissue such as the female breast with a greatly reduced amount of distortion to the phase front of the incident acoustic wave (compared to acoustic waves used in the prior art with frequencies in the range of 3 to 10 MHz). The acoustic frequencies in the preferred embodiment are low enough to enable the reconstruction algorithms to more accurately calculate a final image from the acquired data. However, the acoustic frequencies are high enough to provide adequate spatial resolution in the final image. Also, the use of acoustic frequencies in the range of 100 kHz to 1.5 MHz results in a minimal amount of attenuation of the incident acoustic wave thus enabling the production of a single large format image of a large portion of tissue. 3. The use of discrete frequency, steady-state insonification of the medium 12 which allows an accurate measurement of the phase and amplitude of the acoustic wave for each discrete frequency at each receiver transducer. 4. The capability of recording quickly and storing a very large amount of acoustic data in order to develop the matrices of data from which the images can be calculated. 5. An electronic system for phase locking all aspects of the acoustic data generation and collection.

B. Description of Data Acquisition Process

The data acquisition process involves the systematic insonification of the medium 12 by having each transducer 10.j (j 32 1, . . . , N) sequentially perform as transmitter of continuous wave, single frequency sound (frequency .omega..sub.1) with subsequent reception of the scattered waves by the remaining transducers 10.k (k=1, . . . , N; k.noteq.j). The acquired information is in the form of an N.times.N complex matrix with (N-1).sup.2 coefficients which we call m.sub.jk (.omega..sub..alpha.) (j,k=1, . . . , N) for data acquired at insonification frequency .omega..sub..alpha.. The (N-1).sup.2 =1023.sup.2 coefficients m.sub.jk (.omega..sub..alpha.) (j,k=1, . . . , N; j.noteq.k) may be collected as described in the following example.

The object 9 to be imaged such as a female breast is placed within transducer ring 10. A coupling fluid 11 is placed between the object 9 and transducer ring 10. With an operator signal from keyboard 34, individual components of the system are initialized by computer supervisor 30. The operator initiates a command to the programmable timing and control module 42 to start the acquisition sequence. Timing and control module 42 sends a sequence of commands to command interpreter 50 causing a reset of the transmitter and receiver multiplexers 60 to their initial state.

Arbitrary wave form generator 26 is then triggered which causes the first transmitter 10.1, through transmit multiplexers 76, to broadcast a single frequency (.omega..sub.1 =687 kHz) sinusoidal acoustic signal into the medium 12. The transmitted acoustic wave scatters from the inhomogeneities of medium 12 and impinges upon the other N-1 transducers 10.2, . . . , 10.1024 which transform these acoustic signals into single frequency electrical signals. After a time T=270 microseconds corresponding to twice, the transit time of the acoustic wave across the diameter 2 r.sub.o of the device, the acoustic signals at each transducer have reached a steady-state condition. At this time, the electrical signals from each transducer 10.2 to 10.1024 are routed to the eight two-channel analog-to-digital converters 22 through receive multiplexers 72. The exact data acquisition process for each transducer 10 proceeds as follows. The supervisor 30 instructs the programmable timing and control module 42 to take four quadrature measurements from each transducer 10.2 to 10.1024. These four measurements labeled A,B,C, and D are of the transducer voltage output at four short equally spaced intervals during one acoustic period which is phase locked to the signal from the arbitrary waveform generator 26. For example, if the input frequency is 1 MHz (with a period of 1.times.10.sup.-6 second) samples are taken at 0.25.times.10.sup.-6 second intervals. The real part of complex amplitude is determined by subtracting measurement B from D and the imaginary part is determined by subtracting measurement A from C. The calculated complex amplitude at each transducer contains the phase and amplitude information of the acoustic wave at each transducer 10.2 to 10.1024 as compared to the input signal from 10.1. The data is stored temporally in the buffer memory units of the analog-to-digital converter units 22.

The system has eight two-channel analog-to-digital converters so signals are recorded from the transducers 10 in parallel. First signals from sixteen transducers 10.2, 10.66, 10.129, . . . , 10.961 are recorded simultaneously. After acquiring data from the sixteen transducers 10.2, 10.66, 10.129, . . . , 10.961, the receive multiplexers 72 rout signals from the sixteen transducers 10.3, 10.67, 10.130, . . . , 10.962 to the analog-to-digital converters 22 which acquire data from these transducers. The acquisition process continues until each analog-to-digital converter 22 has acquired data from sixty-four transducers 10. This process takes 270 microseconds for the signals to reach steady-state and 512 microseconds for each analog-to-digital converter 22 to acquire samples from sixty-four transducers 10. Digital values representing the phase and amplitude of each receiving transducer 10.2, . . . , 10.1024 relative to the transmitting transducer 10.1 are recorded in the buffer memory of the analog-to-digital converters 22 as N-1 (i.e. 1023) complex coefficients m.sub.jk (.omega..sub.1) (j=1; k=2, . . . , 1024). The data collection procedure continues by routing the .omega..sub.1 =687 kHz electrical signal to transducer 10.2 which now acts as a transmitter of the acoustic energy. After another time T=270 microseconds required for the system to reach a steady-state condition, electrical signals are sequentially measured from the transducers 10.1,10.3, . . . , 10.1024 as described above in order to collect 1023 complex coefficients m.sub.jk (j=2; k=1,3, . . . , 1024). The measurement procedure continues until each transducer 10.1, . . . , 1.0.1024 has acted as a transmitter of the acoustic wave. At this point, a complete set of data representing complex coefficients m.sub.jk (.omega..sub.1) (j,k=1, . . . , N; j.noteq.k) has been recorded for the frequency .omega..sub.1 =687 kHz. A full set of data at a single acoustic frequency is obtained in approximately 1 second. In the preferred embodiment, the entire preceding measurement procedure is repeated ten times in order to acquire a separate m.sub.jk (.omega..sub..alpha.) at each of ten discrete frequencies .omega..sub..alpha. (.alpha.=1,2, . . . , 10) spaced at 62.5 kHz intervals from 687 kHz to 1.250 MHz.

In order to initially calibrate the device, we acquire a complete set of data from the medium 12 with no object 9 (only the coupling fluid 11). This data will only be used to calibrate the transducers as described in Section D.1. The object 9 is removed and the entire measurement procedure described in the previous paragraph is repeated to acquire m.sub.jk.sup.incident (.omega..sub..alpha.) (j,k=1, . . . N.sub.j,j .noteq.k) for the medium 12 with no object 9 (only the coupling fluid 11) at each of ten discrete frequencies .omega..sub..alpha. (.alpha.=1,2, . . . , 10) spaced at 62.5 kHz intervals from 687 kHz to 1.250 MHz.

C. Calculation of an Image from the Acquired Data

The data m.sub.jk (.omega..sub..alpha.) (j,k=1, . . . , N;.alpha.=1,2, . . . , 10) acquired for medium 12 is uniquely determined by the sound speed c(r,.theta.) and attenuation coefficient g(r,.theta.) of the medium. Our objective is to use this set of data m.sub.jk (.omega..sub..alpha.) to calculate mathematically with a digital computer an image of a 12 mm thick slice of medium 12. This image is based on the calculation of an approximation of the complex scattering potential, S.sup.(.alpha.) (r,.theta.), at all locations (r,.theta.) throughout a 12 mm thick slice of the medium (see FIG. 1).

In the preferred embodiment, we calculate an image of the medium 12 from the acquired data m.sub.jk (.omega..sub..alpha.) at each frequency .omega..sub..alpha.. This image, defined as S.sup.(.alpha.) (r,.theta.), is a two-dimensional map of the medium 12. The third dimension is averaged over the vertical extent of the transducers 10.

Our algorithm for S.sup.(.alpha.) (r,.theta.) is derived from the well known acoustic wave equation known as the Helmholtz equation:

.gradient..sup.2 f(r,.theta.)+S.sup.(.alpha.) (r,.theta.)f(r,.theta.)=0(1) where ##EQU1## represents a single frequency acoustic pressure wave .rho.(r,.theta.) at frequency .omega..sub..alpha.. .rho.(r,.theta.) is the density distribution in the medium. S(r,.theta.)=k.sup.2 (r,.theta.)-k.sub.o.sup.2 -.rho.(r,.theta.).sup.1/2.gradient..sup.2 .rho.(r,.theta.).sup.-1/2 is the complex scattering potential of the medium. .omega..sub..alpha. =angular frequency of acoustic wave. k.sub.o =.omega..sub..alpha. /c.sub.o c.sub.o =sound speed in coupling fluid 11 k(r,.theta.)=[.omega..sub..alpha..sup.2 /c.sup.2 (r,.theta.)-i .omega..mu.(r,.theta.)].sup.1/2 is the complex wavevector which accounts for the spatially varying sound speed c (r,.theta.) and includes an acoustic energy loss term .mu. (r,.theta.) which accounts for viscous heating of the medium 12 and diffraction of acoustic energy from the field of view of the transducers 10.

A solution of equation (1) for f (r.sub.o,.theta..sub.k) at transducer element 10.k (at r.sub.o,.theta..sub.k) due to single frequency insonification (.omega..sub..alpha.) of the medium 12 by transducer 10.j (at r.sub.o,.theta..sub.j) can be expressed in terms of a two dimensional integral equation involving:

(i) the complex scattering potential S.sup.(.alpha.) (r,.theta.), (ii) the position of the transducers (r.sub.o,.theta..sub.j) and (r.sub.o,.theta..sub.k), and (iii) tabulated Hankel (H) and Bessel (J) functions. That solution is ##EQU2##

The term f.sup.(.alpha.) (r.sub.o,.theta..sub.k) represents the pressure wave (frequency .omega..sub..alpha.) at the position of transducer 10.k expressed in terms of amplitude and phase, relative to a pressure wave originating at the position of transducer 10.j. This is precisely equivalent to our acquired data m.sub.jk (.omega..sub..alpha.) which we described in Section B. Thus:

f.sup.(.alpha.) (r.sub.o,.theta..sub.k).tbd.m.sub.jk (.omega..sub..alpha.).(3) Therefore, we can write equation (2) as: ##EQU3## where: S.sup.(.alpha.) (r,.theta.) is a complex image of medium 12 determined from data obtained at frequency .omega..sub..alpha..

It is a two dimensional matrix representing the acoustic scattering potential at (1024).sup.2 (r,.theta.) points within a 12 mm thick slice through medium 12 in the (r,.theta.) plane. Each of the (1024).sup.2 complex values of S.sup.(.alpha.) (r,.theta.) includes a real and an imaginary component.

H'.sub.m (k.sub.o r.sub.o) is the first derivative of a Hankel function representing the antenna pattern of transmitter transducer 10.j and is the correct mathematical representation for the dipole transducers used in the preferred embodiment. Hankel functions are discussed and tabulated in, for example, "The Pocketbook of Mathematical Functions," by Abramowitz and Stegun, Chapter 9, Verlag Harri Deutach 1984. The attenna pattern of any arbitrary transducer can be modeled as a weighted sum of monopole, dipole, quadrapole, etc, attenna patterns. These antenna patterns are discussed, for example, in "Classical Electrodynamics", J. D. Jackson, Chapter 16, I. Wiley and Sons, New York, 1962.

H'.sub.n (k.sub.o r.sub.o) is the first derivative of a Hankel function representing the antenna pattern of receiver transducer 10.k. For the dipole transducer used in the preferred embodiment the functions are the same as for the transmitter transducer. J.sub.m (k.sub.o r) and J.sub.n (k.sub.o r) are Bessel functions determined for all values r within the (r,.theta.) plane based on tabulated Bessel function values.

represents the Fourier transform of the indicated functions given by

g(.kappa.,.phi.).tbd. {f(r,.theta.)}.tbd..intg.f(r,.theta.)e.sup.-i.kappa.rcos(.theta.-.phi.) rdrd.theta. (5)

.sup.-1 represents the inverse Fourier transform of the indicated functions given by

f(r,.theta.).tbd. .sup.-1 {g(.kappa.,.phi.)}.tbd..intg.g(.kappa.,.phi.)e.sup.i.kappa.rcos(.theta.-.phi.) .kappa.d.kappa.d.theta. (6)

Equation (4) can be transformed by appropriate mathematical manipulations into the following form which gives S.sup.(.alpha.) (r,.theta.) in terms of the measured acoustic data m.sub.jk (.omega..alpha.). ##EQU4##

D. Creating the Images

The full process of producing the S.sup.(.alpha.) (r,.theta.) image in this embodiment is described here and with reference to equation (7).

Step (1). The first step is to calibrate the transducers. We do this by first acquiring a set of data m.sub.jk.sup.incident (.omega..sub..alpha.) (j,k=1, . . . , N;j.noteq.k) as described in Section B at each of 10 discrete frequencies .omega..sub..alpha.. (.alpha.=1,2, . . . , 10) with object 9 removed so that medium 12 is only the coupling fluid 11. (Frequencies are spaced at 62.5 kHz intervals between 687 kHz to 1,250 kHz.) Calibration factors are determined by asking the computer to create weighting factors to produce a uniform grey scale image when the water bath only data is processed as given by steps 2 through 9 of this Section. Addition calibration factors are determined to remove the frequency dependence of the transducers and to globally adjust each complex grey scale image of the water bath so as to provide the same image at each discrete frequency .omega..sub..alpha.. These calibration factors are stored as a separate 1024.times.1024 calibration matrix for each of the 10 discrete frequencies .omega..sub..alpha. (.alpha.=1,2, . . . , 10). These calibration matrices are stored in the memory of computer 72.

Step (2). Data representing m.sub.jk (.omega..sub..alpha.) of the medium 12 are acquired as described in Section B for each of the 10 frequencies equally spaced between 687 kHz to 1,250 kHz. For each frequency .omega..sub..alpha., m.sub.jk (.omega..sub..alpha.) is an 1024.times.1024 complex matrix of phase and amplitude data.

Step (3). For each discrete frequency .omega..sub..alpha., computer 72 multiplies the m.sub.jk (.omega..sub..alpha.) data by the calibration matrix determined in Step (1) of this Section to provide calibrated m.sub.jk (.omega..sub..alpha.) data.

Step (4). Using the calibrated m.sub.jk (.omega..sub..alpha.) data, computer 72 calculates the quantity ##EQU5## for each frequency .omega..sub..alpha.. This manipulation constitutes a fast Fourier transform of the calibrated complex matrix m.sub.jk (.omega..sub..alpha.). The result of this transformation is that the calibrated m.sub.jk (.omega..sub..alpha.) data for each frequency .omega..sub..alpha. is now in the form of a 1024.times.1024 matrix in an azimuthal mode space dependent on azimuthal mode values designated by the azimuthal mode numbers p and q.

Step (5). We now calculate the antenna patterns H'.sub.p (k.sub.o r.sub.o) and H'.sub.q (k.sub.o r.sub.o) for each of the transducers 10. This involves tabulation of the first derivatives of Hankel functions for azimuthal mode numbers p,q=-512 to 512. These values are stored in the memory of computer 72.

Step (6). In this step we simply divide the quantity determined in Step (4) by the product of the antenna patterns calculated in Step (5) for each transmitter-receiver combination at frequency .omega..sub..alpha. and each mode combination (p,q) to calculate the following quantity: ##EQU6##

Step (7). We now calculate a weighting function J.sub.p (k.sub.o r)J.sub.q (k.sub.o r)e.sup.-i .theta.(p+q) at each point r and .theta. which is a 1024.times.1024 dimensional matrix. J.sub.p and J.sub.q are determined from a Bessel function table stored in the memory of computer 72. Step (8). Next we multiply the quantity calculated in Step (6) by the weighting function calculated in Step (7) and sum over azimuthal modes p and q to obtain the quantity: ##EQU7## This step is done at each point (r,.theta.) in the medium 12. The result of this step is S.sub.conv.sup.(.alpha.) (r,.theta.). This is an image of S.sup.(.alpha.) (r,.theta.) but it is convolved with [J.sub.o (k.sub.o r)].sup.2 which presents a smoothed image of S.sup.(.alpha.) (r,.theta.). For a more accurate image we need to deconvolve the quantity [J.sub.o (k.sub.o r)].sup.2 from S.sub.conv.sup.(.alpha.) (r,.theta.).

Step (9). Convolution and deconvolution methods are standard mathematical calculations and are discussed, for example, in "Discrete-Time Signal Processing", pg. 58, by Oppenheim and Schafer, Prentice-Hall (1989). Convolutions and deconvolutions of two spatial functions can be calculated in the spatial domain or in the Fourier wavevector domain. The preferred embodiment calculates the convolutions and deconvolutions in the Fourier wavevector domain. However, persons skilled in the art will recognize that these convolutions and deconvolutions can easily be done in the spatial domain. The Fourier transform of [J.sub.o (k.sub.o r)].sup.2 is ##EQU8## Thus we deconvolve [J.sub.o (k.sub.o r)].sup.2 from S.sub.conv.sup.(.alpha.) (r,.theta.) to calculate an image S.sup.(.alpha.) (r,.theta.) as follows:

S.sup.(.alpha.) (r,.theta.)= .sup.-1 {.pi..kappa.(4k.sub.o -.kappa.).sup.1/2 {S.sub.conv.sup.(.alpha.) (r,.theta.)}}.(11)

Using equation (11) we can produce an image S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. of medium 12: To show the images on a digital computer monitor, computer 72 is programmed to convert from polar coordinates (r,.theta.) to cartesian coordinates (x,y) to produce S.sup.(.alpha.) (x,y). Each of the S.sup.(.alpha.) (x,y) is a matrix of values. Each value of the matrix contains two numbers, a first number representing a real part of the value and a second number representing an imaginary part of the value.

The magnitude of these numbers can be represented as a grey scale on a computer monitor at a grid pixel (x,y). (The larger the number the whiter the pixel or vise versa.) When this is done we get two separate images for each frequency, one from the real part of the value of S and one from the imaginary part of S.

E. Combining the Single Frequency Images

As stated above the values of S.sup.(.alpha.) (r,.theta.) at discrete frequencies coot produce good images of medium 12. Over a certain frequency range such as 687 kHz to 1.25 MHz, for example, in which the sound speed and attenuation coefficient of the medium 12 are approximately frequency independent, the images of S.sup.(.alpha.) (r,.theta.) acquired at different frequencies .omega..sub..alpha. ideally would present the same images of the medium 12, independent of frequency .omega..sub..alpha.. However, the images can be complicated by frequency dependent artifacts resulting from reflections of acoustic energy from the faces of the receiver transducers 10.k and from multiple scattering events within the medium 12. We have found that combining images of S.sup.(.alpha.) (r,.theta.) reconstructed at different frequencies coax reduces the image artifacts and improves the image quality. The following discussion centers on seven preferred ways of combining single frequency images S.sup.(.alpha.) (r,.theta.) reconstructed at different frequencies .omega..sub..alpha..

For the purposes of combining images S.sup.(.alpha.) (r,.theta.) formed at different frequencies .omega..sub..alpha., we have found it useful to define a term S.sup.(.alpha.) (r,.theta.) which we calculate as: S.sup.(.alpha.) (r,.theta.)= .sup.-1 {.kappa..sup.-1 {S.sup.(.alpha.) (r,.theta.)}}. (12)

In equation (12), S.sup.(.alpha.) (r,.theta.) is calculated by first calculating an FFT of S.sup.(.alpha.) (r,.theta.), then dividing by .kappa. and calculating the inverse FFT of the result.

1) Sum of the complex images S.sup.(.alpha.) (r,.theta.).

From the acquired data m.sub.jk (.omega..sub..alpha.), we reconstruct separate images of S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . , 10) as given in Steps (1)-(9) of Section D. We then sum these complex images to produce: ##EQU9## We convert to cartesian coordinates (x,y) to produce an image S.sub.total (r,.theta.) which has less artifact and more clarity than each single frequency image S.sup.(.alpha.) (r,.theta.).

2) Sum of magnitudes of S.sup.(.alpha.) (r,.theta.) Images

From the acquired data m.sub.jk (.omega..sub..alpha.), we reconstruct separate images of S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . , 10) as given in Steps (1)-(9) of Section D. We then calculate the magnitudes of these complex images as .vertline.S.sup.(.alpha.) (r,.theta.).vertline.=[S.sup.(.alpha.) (r,.theta.)S.sup.(.alpha.)* (r,.theta.)].sup.1/2 where S.sup.(.alpha.) *(r,.theta.) is the complex conjugate of S.sup.(.alpha.) (r,.theta.). We then sum the magnitudes of these images to produce: ##EQU10## Again we convert to cartesian coordinates to produce an image .vertline.S.sub.total (r,.theta.).vertline. which has less artifacts and more clarity than each image .vertline.S.sup.(.alpha.) (r,.theta.).vertline..

3) Sound Speed and Attenuation Maps Calculated from S.sup.(.alpha.) (r,.theta.) Images

Sound speed maps and attentuation maps of medium 12 can be calculated from S.sup.(.alpha.) (r,.theta.) as follows:

Step (1). From the acquired data m.sub.jk (.omega..sub..alpha.), we reconstruct separate images of S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . , 10) as given in Steps (1)-(9) of Section D. We then calculate the maps S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. as given in equation (12).

Step (2). Computer 72 then calculates the following function ##EQU11## For each point (r,.theta.) in the medium 12, P(r,.theta.,.tau.) represents a synthesized acoustic pulse P(.tau.) focused at the point (r,.theta.).

Step (3). For each point (r,.theta.) in the medium 12, P(.tau.) has a maximum value P.sub.max at a value of .tau. which we call .tau..sub.max. The computer 72 is programmed to find, for each point (r,.theta.), P.sub.max (r,.theta.) and to determine the value of .tau..sub.max (r,.theta.).

Step (4). The computer 72 then creates a sound speed map c(r,.theta.) of medium 12 by the following calculation

c(r, .theta.)= .sup.-1 {.kappa. {2r.sub.o /.tau..sub.max (r,.theta.)}}c.sub.o. (16) In equation (16), r.sub.o is the radius of the transducer ring and c.sub.o is the sound speed of the coupling fluid 11.

Step (5). The computer then creates an attenuation map .mu.(r,.theta.) of medium 12 by the following calculation ##EQU12## In equation (17), .mu..sub.o is the attenuation coefficient of the coupling fluid 11. These values of c(r,.theta.) and .mu.(r,.theta.) are 1024.times.1024 matrixes in r and .theta.. They are converted to cartesian coordinates (x,y) and displayed on a computer screen with the values of c and .mu. depicted as gray scale values.

4) Refined Sound Speed and Attenuation Maps Calculated from S.sup.(.alpha.) (r,.theta.)

The sound speed c(r,.theta.) and attenuation .mu.(r,.theta.) maps calculated in equations (15)-(17) represent approximate images of the medium 12. This section discusses refined calculations of c(r,.theta.) and .mu.(r,.theta.) from the single frequency images S.sup.(.alpha.) (r,.theta.) which represent more accurate images of the medium 12.

Step (1). From the acquired data m.sub.jk (.omega..sub..alpha.), we reconstruct separate images of S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . , 10) as given in Steps (1)-(9) of Section D. We then calculate the maps S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. as given in equation (12).

Step (2). We then calculate the maps S.sub..DELTA.p,.DELTA.q.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . , 10) as such: ##EQU13## The maps S.sub..DELTA.p,.DELTA.q.sup.(.alpha.) (r,.theta.) represent coupling of acoustic energy from azimuthal modes p and q to modes p+.DELTA.p and q+.DELTA.q.

Step (3). Computer 72 then calculates the following function: ##EQU14## We have had excellent results using sums of .DELTA.p=-6 to +6, .DELTA.q=0; .DELTA.p=0, .DELTA.q=-6 to +6; and .DELTA.p=-6 to +6, .DELTA.q=-6 to +6. Other combinations of .DELTA.p and .DELTA.q could be used.

For each point (r,.theta.) in the medium 12, P.sub.corr (r,.theta.,.tau.) represents a synthesized acoustic pulse P.sub.corr (.tau.) focused at the point (r,.theta.). P.sub.corr (r,.theta.,.tau.) is a refined version of the function P(r,.theta.,.tau.) given in equation (15).

Step (4). For each point (r,.theta.) in the medium 12, P.sub.corr (.tau.) has a maximum value P.sub.max at a value of .tau. which we call .tau..sub.max. The computer 72 is programmed to find, for each point (r,.theta.), P.sub.max (r,.theta.) and to determine the value of .tau..sub.max (r,.theta.) as was done in Step (3) of Section E.3 except we use P.sub.corr (r,.theta.,.tau.) instead of P(r,.theta.,.tau.)).

Step (5). Using the refined values of .tau..sub.max (r,.theta.) and P.sub.max (r,.theta.), we now calculate the sound speed c(r,.theta.)) and attenuation .mu.(r,.theta.) maps as we did in equations (16) and (17).

5) Sum of Phase Aberration Corrected S.sup.(.alpha.) (r,.theta.) Images

The complex images S.sup.(.alpha.) (r,.theta.) reconstructed as given in Steps (1)-(9) of Section F represent approximate images of the medium 12. The complex phase of these images can be distorted at each point (r,.theta.) due to limitations of the reconstruction algorithm as given in Steps (1)-(9) of Section D. This Section discusses a correction method for the phase aberrations in the images S.sup.(.alpha.) (r,.theta.) and the combination of these images to reduce image artifacts.

Step (1). From the acquired data m.sub.jk (.omega..sub..alpha.), we reconstruct separate images of S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. (.alpha.=1,2, . . . ; 10) as given in steps (1)-(9) of Section D. We then calculate the maps S.sup.(.alpha.) (r,.theta.) at each frequency .omega..sub..alpha. as given equation (12).

Step (2). We calculate the map .tau..sub.max (r,.theta.) as described by the procedures given in Section E.3 or section E.4.

Step (3). Computer 72 calculates the map S.sub.total (r,.theta.) as such: ##EQU15##

Step (4). The final image S.sub.total (r,.theta.) is determined by the following calculation:

S.sub.total (r,.theta.)= .sup.-1 {.kappa. {S.sub.total (r,.theta.)}}(21) This image displays only high spatial frequency components (i.e., fine details) of medium 12. These fine details include sharp edges or discontinuities in the sound speed in the medium 12 and also small objects with a scale size comparable to the acoustic wavelengths.

6) Composite Image

The image S.sub.total (r,.theta.) calculated in Section E.5 can be combined with the map .tau..sub.max (r,.theta.) to provide a composite map S.sub.comp (r,.theta.) which is a more accurate representation of the medium 12. The S.sub.comp (r,.theta.) image contains both the high spatial frequency components and low spatial frequency components of the medium 12.

Step (1). The computer 72 calculates .tau..sub.max (r,.theta.) as described by the procedures given in Section E.3 or E.4 and S.sub.total (r,.theta.) given by equation (20).

Step (2). Computer 72 then calculates: ##EQU16## In equation (22), In denotes the natural logarithm and .DELTA..omega. is the spacing between frequencies .omega..sub..alpha. (.DELTA..omega.=62.5 kHz in the preferred embodiment).

Step (3). Computer 72 then calculates S.sub.comp as follows:

S.sub.comp (r,.theta.)= .sup.-1 {.kappa. {e,cir S.sub.comp (r,.theta.)}}.(23)

Computer 72 then converts S.sub.comp (r,.theta.) to cartesian coordinates (x,y) and displays on computer monitor.

7) Further Refined Sound Speed Map Calculated From S.sup.(.alpha.) (r,.theta.)

This section discusses a method to further refine the accuracy of the sound speed maps c(r,.theta.) calculated in Sections E.3 and E.4. In Sections E.3 and E.4 we calculate at each point (r,.theta.) the transit time .tau..sub.max (r,.theta.) as the value .tau. in which the pulse P(.tau.) has a maximum value P.sub.max and derive sound speed maps from this information. In this Section, we calculate the transit time .tau..sub.av (r,.theta.) as the averaged arrival time of the power contained in the pulse P(.tau.). This map .tau..sub.av (r,.theta.) is then used to calculate refined maps of the sound speed c(r,.theta.) of the medium 12. Step (1). Computer 72 computes the following quantity: ##EQU17## where S.sub..DELTA.p,.DELTA.q.sup.(.alpha.) (r,.theta.) is calculated as in equation (18) and * denotes the complex conjugate. The preferred embodiment includes sums of 1).DELTA.p=-6 to +6; .DELTA.q=0; 2).DELTA.p=0, .DELTA.q=-6+6; or 3).DELTA.p=-6+6, .DELTA.q=-6 to +6. Other combinations of .DELTA.p and .DELTA.q could be used.

Step (2). The computer 72 calculates an average transit time .tau..sub.av (r,.theta.) image of medium 12 by the following calculation ##EQU18##

Step (3). The computer 72 then calculates a sound speed image c(r,.theta.) of medium 12 by the following calculation

c(r,.theta.)= .sup.-1 {.kappa. {2r.sub.o /.tau..sub.av (r,.theta.)}}+c.sub.o. (26)

Step (4). The computer 72 then creates an attenuation image .mu.(r,.theta.) of medium 12 by the following calculation ##EQU19##

F. Computer Simulations

We created computer simulations in order to check our algorithms. No hardware other than a computer was used in these reconstructions. FIGS. 4-7 display images reconstructed from simulated data at 1 MHz, using the same transducer configuration as our working prototype system. The simulated object shown in FIG. 3 is a 3-inch diameter cylinder 2 of a homogeneous material having a sound speed which is 5% greater than the surrounding coupling fluid. The two 1 cm diameter voids 1 in the object have the same sound speed as the coupling fluid 3 surrounding the object. The attenuation coefficients of the object and the coupling fluid are equal to zero (no attenuation). The images were produced by the following methods: FIG. 4 displays the real (right) and imaginary (left) components of a single frequency reconstruction as described in Section D. FIG. 5 displays the sound speed and attenuation images as described in Section E.4. FIG. 6 displays the real (right) and imaginary (left) components of a ten single frequency reconstructions which have been corrected for phase aberrations and summed as described in Section E.5. FIG. 6 displays a composite image as described in Section E.6.

G. Actual Images

FIGS. 8-12 display actual images obtained using the prototype device described in Sections A-D. FIG. 8 displays nine images of a 1.2 cm slice of a female breast. The acquired data is the same for all nine images, the reconstruction methods are different for each image. Beginning from the top left corner and progressing left to right the images are 1) real component of single frequency image (Section D), 2) imaginary component of a single frequency image (Section D), 3) magnitude of five single frequency images which have been summed (Section E.1), 4) real component of single frequency image which has been corrected for phase aberrations (Section E.5), 5) imaginary component of single frequency image which has been corrected for phase aberrations (Section E.5), 6) magnitude of five single frequency images each which has been corrected for phase aberrations (Section E.5), 7) sound speed image (Section E.7), 8) attenuation image (Section E.7), and 9) square root of sum of squares of sound speed and attenuation images (Section E.7). FIG. 9 displays nine images of an excised porcine kidney obtained using the prototype device described in Sections A-B and reconstructed with the methods described for FIG. 8.

FIG. 10 displays nine images of a female breast obtained with the prototype device described in Sections A-B anti with reference to FIG. 12. The female lays prone face down on a table and positioned her left breast through an eight inch hole into a water bath containing the transducer assembly 80. Beginning at the top left corner of FIG. 10, the nine images were obtained by sequentially moving the transducer assembly 80 from the area of the female breast 86 close to the chest wall to the nipple of the breast in 1 centimeter steps. The images are each magnitudes of the sum of five single frequency images each which has been corrected for phase aberrations as described in Section E.5. FIG. 11 displays similar images for another female breast in which ther is a large fluid filled cyst.

FIG. 12 displays an image of an excised porcine kidney obtained using the prototype device described in Sections A-B and reconstructed with the method described in Section E.5.

H. Biopsy Needle Guidance

The invention can be used for the monitoring of a biopsy needle as it is inserted in a volume of human tissue such as, for example, the female breast. FIG. 15 displays an illustration of the invention for this purpose. The coupling of the acoustic energy from the transducer ring 91 to the breast 90 may be accomplished with the use of a fluid filled bladder which conforms to the breast. Tomographic images of several slices of the breast are produced as in FIG. 10 to locate the biopsy location. The volume of the breast to be biopsied is determined from these tomographic images. A computer analyses the tomographic image and determines where tip of the biopsy needle 93 is to be placed. The computer then guides the biopsy needle tip to the predetermined location. A second set of tomographic images are obtained with the biopsy needle in place to verify the position of the needle tip.

I. Images of Temperature Distributions

This invention provides a device and a method for measuring variations in temperature of the medium 12. This is accomplished by virtue of the fact that temperature changes in human tissue result in changes in the sound speed of the tissue. The sound speed of human tissue which has been locally heated can be expressed as

c(r,.theta.){heated}=c(r,.theta.)+.DELTA.T(r,.theta.) [d/dT c(r,.theta.)](26) where .DELTA.T(r,.theta.) is the temperature distribution resulting from local heating of the tissue, c(r,.theta.) is the sound speed map of the tissue at the ambient body temperature, and d/dT c(r,.theta.) is the differential change in sound speed with respect to temperature. For example, in human liver tissue at 37.degree. C., c(r,.theta.)=1596 m/sec and d/dT c(r,.theta.)=0.96 m/(sec -.degree.C.). A .DELTA.T=5.degree. C. change in temperature due to local heating of the liver by laser or microwave heating would produce a change in sound speed of [c(r,.theta.).sup.-1 ].DELTA.T(r,.theta.) [d/dT c(r,.theta.)]=0.3%.

The invention can be used to produce a combination sound speed and temperature image by locally heating a volume of the tissue and then by directly imaging the sound speed of the tissue as prescribed in the procedures given in Sections A-F. Alternatively, one can produce a sound speed image of the tissue before and after locally heating a volume of the tissue and then subtracting or dividing the two images to produce an image of the temperature distribution.

J. Alternate Embodiments

An alternate embodiment involves the use of a periodic wideband insonification signal in which a periodic time-varying signal is transmitted into the medium sequentially from each transducer. The signal has a 1 MHz center frequency, a 600 kHz frequency bandwidth, and repeats every 16 microseconds. The Fourier transform of this signal is a comb of ten discrete frequencies evenly spaced in 62.5 kHz intervals from 687 KHz to 1.25 MHz. The relative acoustic signal strengths of the discrete frequency components is fairly uniform across the frequency range. The relative phases of the discrete frequency components are randomized in order to provide a fairly flat temporal response of the periodic time-varying signal. The wideband signal reaches a steady state condition after approximately 270 microseconds which is twice the acoustic travel time across the diameter of the transducer ring. After this time, each receiver sees a time varying signal that repeats every 16 microseconds. Data is acquired and stored from each receiver for 16 microseconds at a 4 MHz data rate. The acquisition method is repeated with each transducer acting as a transmitter of sound. A Fourier transform of the time-varying receiver signals provides multi-frequency data with a frequency resolution of 62.5 kHz (1/(16 microseconds)). The data set is comparable to the data taken by the acquisition method described in Sections C and D in which 10 separate data sets are acquired at ten discrete frequencies spaced 62.5 kHz apart from 650 kHz to 1.25 MHz. However, the total acquisition time is reduced to 3 seconds with the wideband acquisition method.

Alternate embodiments of the invention also include the geometry of the transducer array. These alternate embodiments may be appropriate to facilitate coupling of the acoustic waves to various parts of the human body including, for example, the abdomen, female breast, cranial cavity, neck, arm and thigh or to other non-human subjects. An embodiment designed to image the abdomen region of human patients is shown in FIG. 14. A device to perform mammography examination is shown in FIG. 13. In this embodiment the ring of 1024 transducers 80 is moveable vertically in an oil tank 82 surrounding water tank 84 so that a complete scan of female breast 86 can be obtained. These images could be viewed separately as shown in FIG. 10 or could be combined in the computer to display a three-dimensional image of the breast. An alternate embodiment includes a similar device to image arteries and veins in arms and legs. Alternate embodiments include a ring of acoustic transducers which do not lie on the locus of a circle, a ring of transducers which do not completely encircle the medium which is to be imaged, a set of transducers which is in segments of various lengths which partially circumscribe the medium, and a set of transducers arranged in parallel rows on opposite sides of the medium to be imaged. Another embodiment which concerns the use of a rubber fluid-filled bladder to couple the acoustic waves from the transducers to the medium which is shown in FIG. 13. Various numbers of transducers could be utilized but the number should be at least eight, and as indicated above, transducers should preferably be spaced no more farther apart than one half wavelength of the acoustic wave being transmitted. For non-circular arrays, the algorithm derived above would have to be recalculated to account for the appropriate boundary conditions.

Alternate embodiments concern the adaptation of the invention to probe medium with electromagnetic waves, with energy spectrums ranging from 1 Hz to x-rays energies, including radiowave, microwave, infrared light, visible light, ultraviolet light imaging devices.

We have tested our prototype at frequency in the range of 500 kHz to 1.25 MHz with good results and we expect that these range could be extended with our technique from 100 kHz to 1.5 MHz.

Although the preferred embodiments described above use the same transducers to broadcast and receive, separate means for broadcasting and receiving could be rotated around the medium to simulate a ring of broadcasting transducers. A similar arrangement could be provided to rotate a receiver to simulate a ring of receiving transducers.

The reader should construe the above described embodiments of this invention as examples and the scope of this invention shall be determined by the appended claim and their legal equivalents.

Claims

1. An acoustic imaging device for providing an image of at least a portion of a medium comprising:

A) an acoustic broadcasting means for broadcasting into said medium acoustic signals from at least seven broadcast locations on a circle at least partially surrounding said medium, one location at a time said acoustic signals being at least one discrete frequency and continuous for a period of time equal to at least the travel time of said acoustic signal across said medium,

B) an acoustic signal detection means for detecting, at a plurality of detection locations on said circle at least partially surrounding said medium, acoustic signals broadcast by said broadcasting means,

C) a data collection and comparison means for comparing signals received to said signals broadcast in order to provide with respect to each of said plurality of locations a set of data representing the phases and amplitudes of acoustic signals received at said plurality of locations relative to the phase and amplitude of said signals broadcast from said at least seven locations, and

D) a computing means for mathematically constructing at least one image of at least a portion of said medium utilizing said phase and amplitude data said computing means being programmed to accomplish the following steps:

1) transform said phase and amplitude data into a set of data defining an azimuthal data set in an azimuthal mode space said set of data representing azimuthal modes at the location of said broadcasting means and said detection means,

2) modify said azimuthal data set to account for transducer antenna patterns to produce an antenna independent azimuthal data set, and

3) multiply said antenna independent azimuthal data set by weighting function calculated at a large number of locations in said medium and sum over said azimuthal modes to produce a convolved image of said medium, and

4) deconvolved said convolved image of said medium to produce said image of said medium.

2. An acoustic imaging device as in claim 1 wherein said computing means is programmed to accomplish the additional step of sequentially accomplishing steps A-D at a plurality of discrete frequencies so as to construct a plurality of single frequency complex images.

3. An acoustic imaging device as in claim 2 wherein said computing means is programmed to accomplish the additional step of summing said plurality of single frequency complex images to produce a combined complex image.

4. An acoustic imaging device as in claim 2 wherein said computing means is programmed to accomplish the additional step of summing the magnitudes of said plurality of single frequency complex images to produce a combined image.

5. An acoustic imaging system as in claim 2 wherein said computing means is programmed to accomplish the following additional steps:

A) utilize said plurality of single frequency complex images to construct a set of time data representing a synthesized acoustic pulse which sequentially focuses at said large number of locations in said medium, and

B) utilizing said set of time data, calculate an image of said medium.

6. An acoustic imaging system as in claim 5 wherein said computing means is programmed to utilize said set of time data to determine the time-of-arrival data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said time-of-arrival data forming a time-of-arrival image.

7. An acoustic imaging system as in claim 6 wherein said computer means is programmed to accomplish the additional step of deconvolving a function of said time-of-arrival image to produce a sound speed image.

8. An acoustic imaging device as in claim 7 wherein said function is the diameter of said circle divided by said time-of-arrival image.

9. An acoustic imaging system as in claim 6 wherein said computing means is programmed to accomplish the following additional steps:

A) convolve each said single frequency complex image to produce a plurality of convolved single frequency complex images,

B) multiply each said convolved single frequency complex image by an exponential function of said time-of-arrival image to produce a plurality of phase aberration corrected convolved single frequency complex images,

C) sum said plurality of phase aberration corrected convolved single frequency complex images to produce a phase aberration corrected convolved complex image, and

D) deconvolve said phase aberration corrected convolved complex image to produce a phase aberration corrected complex image.

10. An acoustic imaging device as in claim 9 wherein said computing means is programmed to accomplish the additional step of calculating the magnitude of said phase aberration corrected complex image to produce a phase aberration corrected image.

11. An acoustic imaging system as in claim 9 wherein said computing means is programmed to accomplish the following additional steps:

A) sum a function of the logarithm of said convolved phase aberration corrected complex image with said time-of-arrival image to produce a convolved composite image, and

B) deconvolve said convolved composite image to produce a composite image.

12. An acoustic imaging system as in claim 5 wherein said computing means is programmed to utilize said set of time data to determine the maximum pulse power data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said maximum pulse power data forming a peak pulse power image.

13. An acoustic imaging system as in claim 12 wherein said computing means is programmed to accomplish the additional steps of a function of said maximum pulse power image to produce an attenuation image.

14. An acoustic imaging device as in claim 13 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

15. An acoustic imaging system as in claim 2 wherein said computing means is programmed to accomplish the following additional steps:

A) utilize said plurality of single frequency complex images to construct a set of time data representing a synthesized acoustic pulse which sequentially focuses at said large number of locations in said medium,

B) calculate refined set of time data from said set of time data by summing contributions to said synthesized acoustic pulse from said azimuthal modes numerically close to said azimuthal mode data set, and

C) utilizing said refined set of time data, calculate an image.

16. An acoustic imaging system as in claim 15 wherein said computing means is programmed to utilize said refined set of time data to determine the refined time-of-arrival data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said refined time-of-arrival data forming a refined time-of-arrival image.

17. An acoustic imaging system as in claim 16 wherein said computing means is programmed to accomplish the additional step of deconvolving a function of said refined time-of-arrival image to produce a refined sound speed image.

18. An acoustic imaging device as in claim 17 wherein said function is the diameter of said circle divided by said refined time-of-arrival image.

19. An acoustic imaging system as in claim 16 wherein said computing means is programmed to accomplish the following additional steps:

A) convolve each said single frequency complex image to produce a plurality of convolved single frequency complex images,

B) multiply each said convolved single frequency complex image by an exponential function of said refined time-of-arrival image to produce a plurality of refined phase aberration corrected convolved single frequency complex images,

C) sum said plurality of refined phase aberration corrected convolved single frequency complex images to produce a refined phase aberration corrected convolved complex image, and

D) deconvolve said refined phase aberration corrected convolved complex image to produce a refined phase aberration corrected complex image.

20. An acoustic imaging device as in claim 19 wherein said computing means is programmed to accomplish the additional step of calculating the magnitude of said refined phase aberration corrected complex image to produce a refined phase aberration corrected image.

21. An acoustic imaging system as in claim 19 wherein said computing means is programmed to accomplish the following additional steps:

A) sum a function of the logarithm of said refined convolved phase aberration corrected complex image with said refined time-of-arrival image to produce a refined convolved composite image, and

B) deconvolve said refined convolved composite image to produce a refined composite image.

22. An acoustic imaging system as in claim 15 wherein said computing means is programmed to utilize said refined set of time data to determine the refined maximum pulse power data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said refined maximum pulse power data forming a refined peak pulse power image.

23. An acoustic imaging system as in claim 22 wherein said computing means is programmed to accomplish the additional steps of deconvolving a function of said refined maximum pulse power image to produce a refined attenuation image.

24. An acoustic imaging device as in claim 23 wherein said function is the logarithm of said refined maximum pulse power image divided by the diameter of said circle.

25. An acoustic imaging device as in claim 1 wherein said at least discrete frequency is a plurality of discrete frequencies.

26. An acoustic imaging device as in claim 25 wherein said computing means is programmed to accomplish steps 1-4 with respect to said plurality of identifiable discrete frequencies so as to construct a plurality of single frequency complex images.

27. An acoustic imaging device as in claim 26 wherein said computing means is programmed to accomplish the additional step of summing said plurality of single frequency complex images to produce a combined complex image.

28. An acoustic imaging device as in claim 26 wherein said computing means is programmed to accomplish the additional step of summing the magnitudes of said plurality of single frequency complex images to produce a combined image.

29. An acoustic imaging system as in claim 26 wherein said computing means is programmed to accomplish the following additional steps:

A) utilize said plurality of single frequency complex images to construct a set of time data representing a synthesized acoustic pulse which sequentially focuses at said large number of locations in said medium, and

B) utilizing said set of time data, calculate an image.

30. An acoustic imaging system as in claim 29 wherein said computing means is programmed to utilize said set of time data to determine the time-of-arrival data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said time-of-arrival data forming a time-of-arrival image.

31. An acoustic imaging system as in claim 30 wherein said computer means is programmed to accomplish the additional step of deconvolving a function of said time-of-arrival image to produce a sound speed image.

32. An acoustic imaging device as in claim 31 wherein said function is the diameter of said circle divided by said time-of-arrival image.

33. An acoustic imaging system as in claim 29 wherein said computing means is programmed to utilize said set of time data to determine the maximum pulse power data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said maximum pulse power data forming a peak pulse power image.

34. An acoustic imaging system as in claim 33 wherein said computing means is programmed to accomplish the additional steps of deconvolving a function of said maximum pulse power image to produce an attenuation image.

35. An acoustic imaging device as in claim 34 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

36. An acoustic imaging system as in claim 26 wherein said computing means is programmed to accomplish the following additional steps:

A) utilize said plurality of single frequency complex images to construct a set of time data representing a synthesized acoustic pulse which sequentially focuses at said large number of locations in said medium,

B) calculate refined set of time data from said set of time data by summing contributions to said synthesized acoustic pulse from said azimuthal modes numerically close to said azimuthal mode data set, and

C) utilizing said refined set of time data, calculate an image of said medium.

37. An acoustic imaging system as in claim 36 wherein said computing means is programmed to utilize said refined set of time data to determine the refined time-of-arrival data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said refined time-of-arrival data forming a refined time-of-arrival image.

38. An acoustic imaging system as in claim 37 wherein said computing means is programmed to accomplish the additional step of deconvolving a function of said refined time-of-arrival image to produce a refined sound speed image.

39. An acoustic imaging device as in claim 38 wherein said function is the diameter of said circle divided by said time-of-arrival image.

40. An acoustic imaging system as in claim 37 wherein said computing means is programmed to accomplish the following additional steps:

A) convolve each said single frequency complex image to produce a plurality of convolved single frequency complex images,

B) multiply each said convolved single frequency complex image by an exponential function of said refined time-of-arrival image to produce a plurality of refined phase aberration corrected convolved single frequency complex images,

C) sum said plurality of refined phase aberration corrected convolved single frequency complex images to produce a refined phase aberration corrected convolved complex image, and

D) deconvolve said refined phase aberration corrected convolved complex image to produce a refined phase aberration corrected complex image.

41. An acoustic imaging system as in claim 45 wherein said computing means is programmed to accomplish the following additional steps:

A) sum a function of the logarithm of said refined convolved phase aberration corrected complex image with said refined time-of-arrival image to produce a refined convolved composite image, and

B) deconvolve said refined convolved composite image to produce a refined composite image.

42. An acoustic imaging system as in claim 36 wherein said computing means is programmed to utilize said refined set of time data to determine the refined maximum pulse power data of said synthesized acoustic pulse from said broadcast locations to said large number of locations within said medium and then to said detection locations said refined maximum pulse power data forming a refined peak pulse power image.

43. An acoustic imaging system as in claim 42 wherein said computing means is programmed to accomplish the additional steps of deconvolving a function of said refined maximum pulse power image to produce a refined attenuation image.

44. An acoustic imaging device as in claim 43 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

45. An acoustic imaging system as in claim 36 wherein said computing means is programmed to accomplish the following additional steps:

A) convolve each said single frequency complex image to produce a plurality of convolved single frequency complex images,

B) multiply each said convolved single frequency complex image by an exponential function of said time-of-arrival image to produce a plurality of phase aberration corrected convolved single frequency complex images,

C) sum said plurality of phase aberration corrected convolved single frequency complex images to produce a phase aberration corrected convolved complex image, and

D) deconvolve said phase aberration corrected convolved complex image to produce a phase aberration corrected complex image.

46. An acoustic imaging device as in claim 45 wherein said computing means is programmed to accomplish the additional step of calculating the magnitude of said phase aberration corrected complex image to produce a phase aberration corrected image.

47. An acoustic imaging device as in claim 45 wherein said computing means is programmed to accomplish the additional step of calculating the magnitude of the refined phase aberration corrected complex image to produce a refined phase aberration corrected image.

48. An acoustic imaging system as in claim 45 wherein said computing means is programmed to accomplish the following additional steps:

A) sum a function of the logarithm of said convolved phase aberration corrected complex image with said time-of-arrival image to produce a convolved composite image, and

B) deconvolve said convolved composite image to produce a composite image.