WO1995028883A1

WO1995028883A1 - Acoustic imaging device

Info

Publication number: WO1995028883A1
Application number: PCT/US1995/005078
Authority: WO
Inventors: Brett A. Spivey; Peter J. Martin; Douglas A. Palmer; Gregory Otto; Robert M. Cram
Original assignee: Thermotrex Corporation
Priority date: 1994-04-25
Filing date: 1995-04-24
Publication date: 1995-11-02
Also published as: JP3133764B2; JPH08508925A; EP0705073A1

Abstract

This invention provides an acoustic imaging device in which a large number of transducers (10) are spaced less than half acoustic wavelength apart on a circle surrounding an object (9) to be imaged. A signal generator (26) generates discrete acoustic frequencies in the range of 100 kHz to 1.5 MHz. Multiplexer systems (1-8) are provided to permit each transducer, one at a time, to broadcast a signal while the broadcast signal is detected by the other transducers. Electronic equipment records the detected signal and from the recorded information phase and amplitude data is calculated with respect to each transducer location. A computer (38) programmed with an algorithm computes images of slices through the object using the phase and amplitude data. We utilize single frequency, steady-state acoustic signals which enables the accurate determination of the phase and amplitude of the acoustic signals at each receiver transducer. Finally, in order to reduce the data acquisition time and reduce motion artifacts, we utilize acoustic signals comprised of a plurality of discrete identifiable frequencies each frequency providing an image of the medium.

Description

ACOUSTIC IMAGING DEVICE

BACKGROUND OF THE INVENTION

This application is a continuation in-part of Serial No. 07/891,851 filed June 1, 1992 which is a continuation in part of Serial No. 07/708,354 filed May 31, 1991. This invention relates to imaging systems and methods and in particular to systems and methods for creating internal images of objects using mathematical procedures and acoustic wave data.

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 (= 1 A). The mathematical reconstruction techniques assume a geometrical ray approximation and are quite simplified.

In ultrasonic tomography, the acoustic wavelengths (== 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 Sine 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 (US Patent 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.

The present invention provides an acoustic imaging device in which a large number of transducers are spaced less than a half acoustic wavelength apart on a circle surrounding an object to be imaged. A signal generator generates discrete acoustic frequencies in the range of 100 kHz to 1.5 MHz. Multiplexer systems are provided to permit each transducer, one at a time, to broadcast a signal while the broadcast signal is detected by the other transducers. Electronic equipment records the detected signal and from the recorded information phase and amplitude data is calculated with respect to each transducer location. A computer programmed with an algorithm computes images of slices through the object using the phase and amplitude data. We utilize single frequency, steady-state acoustic signals which enables the accurate determination of the phase and amplitude of the acoustic signals at each receiver transducer. Finally, in order to reduce the data acquisition time and reduce motion artifacts, we utilize acoustic signals comprised of a plurality of discrete indentifiable frequencies each frequency providing an image of the medium.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram for use in illustrating the mathematical reconstruction algorithm. FIG. 3 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,θ) (left) and attenuation μ(r,θ) (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₀ = 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,θ) (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 p₀ and sound speed c₀ 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 λ₀ = 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 J2 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,θ) 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,θ) 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-l 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 λo = 1.5 mm in water and approximately this wavelength in human tissue. The spacing between centerlines of adjacent transducers is somewhat less than J2 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 λjl, each transducer has a dipole attenna pattern in the (r,θ) plane. This means that the relative amplitude of the transmitted acoustic energy follows a cos ψ function of the angle ψ 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

devices.

3) Address Decoder

,.c FTV.

Our address decoder 79 is a model EB910 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 -JOT 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 λ/2 apart on a circle surrounding the object where λ 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 lO.j (j = 1,...,N) sequentially perform as transmitter of continuous wave, single frequency sound (frequency α>ι) with subsequent reception of the scattered waves by the remaining transducers lO.k (k = 1,...,N; k N j). The acquired information is in the form of an N x N complex matrix with (N-l)² coefficients which we call m_jk (ω_α) (j,k=l,..,N) for data acquired at insonification frequency ω_α. The (N-l)² = 1023² coefficients mj (ω_α) (j,k = 1,...,N; jNk) 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 (ωj = 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-l 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_D 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 x 10^-6 second) samples are taken at .25 x 10^-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-l (i.e. 1023) complex coefficients m_jk (ω (j =1; k = 2,...,1024). The data collection procedure continues by routing the α>ι = 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 mjk 0 = ²; k = 1,3,-, 1024). The measurement procedure continues until each transducer 10.1,...,10.1024 has acted as a transmitter of the acoustic wave. At this point, a complete set of data representing complex coefficients m_j (ωj) (j,k =

1 N; j N k) has been recorded for the frequency ωi = 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_jk (ω_α) at each of ten discrete frequencies ω_α (α = 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.l. The object 9 is removed and the entire measurement procedure described in the previous paragraph is repeated to acquire m^m _; ^cιdent (ω_α) ( ,k = l,..,N;j ≠ k) for the medium 12 with no object 9 (only the coupling fluid 11) at each of ten discrete frequencies ω (α = 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_jk (ω_α) (j,k = l,..,N;α = 1,2 10) acquired for medium 12 is uniquely determined by the sound speed c(r,θ) and attenuation coefficient μ(r,θ) of the medium. Our objective is to use this set of data m_jk (ω_α) 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^(α)(r,θ), at all locations (r,θ) 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 n-_jk(a>a) at each frequency ω_α. This image, defined as S^(α>(r,θ), 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^(α)(r,θ) is derived from the well known acoustic wave equation known as the Helmholtz equation:

V² f (r,θ) + S⁽«⁾(r,θ) f (r,θ) = 0 (1)

where: f(r,θ) = ' _1/2 represents a single frequency acoustic pressure wave p(r,θ) at p(r_»θ) frequency ω_α. p (r,θ) is the density distribution in the medium.

S(r,θ) = k²(r,θ) - k² - p(r,θ)^1/2 V²p(r,θ)^"1/2 is the complex scattering potential of the medium. ω_α = angular frequency of acoustic wave. ko = ω_α/c₀ c₀ = sound speed in coupling fluid 11 k(r,θ) = [ω_α ²/c²(r,θ) - iω N(r,θ)] ^{1 2} is the complex wavevector which accounts for the spatially varying sound speed c (r,θ) and includes an acoustic energy loss term N (r,θ) 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₀,θk) at transducer element lO.k (at r₀,θ_k) due to single frequency insonification (ω_α) of the medium 12 by transducer lO.j (at r₀,θ_j) can be expressed in terms of a two dimensional integral equation involving:

(i) the complex scattering potential S^(α)(r,θ),

(ii) the position of the transducers (r₀,θ_j) and (r₀,θ_k), and

(iii) tabulated Hankel (H) and Bessel (J) functions.

That solution is

m,n

The term tf^α)(r₀,θ_k) represents the pressure wave (frequency ω_α) at the position of transducer lO.k expressed in terms of amplitude and phase, relative to a pressure wave originating at the position of transducer lO.j. This is precisely equivalent to our acquired data m_jk(ω_α) which we described in Section B. Thus:

Therefore, we can write equation (2) as:

m_jk(ω_α) = ∑jH'_m(k₀r H (k₀r₀ -^imθ^^inθk J_m(k₀r)J_n(k₀r)e^iθ(m+n^^α)(r,θ)rdrdθ (4)

where:

S^(α)(r, θ) is a complex image of medium 12 determined from data obtained at frequency ω_α- It is a two dimensional matrix representing the acoustic scattering potential at (1024)² (r,θ) points within a 12 mm thick slice through medium 12 in the (r,θ) plane. Each of the (1024)² complex values of S^{( )}(r,θ) includes a real and an imaginary component. H'rn(k₀r₀) ^{l& me} first derivative of a Hankel function representing the antenna pattern of transmitter transducer lO.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, 1. Wiley and Sons, New York, 1962.

H'_n(k_or₀) is the first derivative of a Hankel function representing the antenna pattern of receiver transducer lO.k. For the dipole transducer used in the preferred embodiment the functions are the same as for the transmitter transducer.

Jm(ko ) and J_n(kor) are Bessel functions determined for all values r within the (r,θ) plane based on tabulated Bessel function values.

3 represents the Fourier transform of the indicated functions given by

g(κ,φ) -≡ 3{f(r,θ)} ≡ Jf(r,θ)e-ⁱ'^crcos(θ-^<,,)rdrdθ (5)

3"¹ represents the inverse Fourier transform of the indicated functions given by

f(r,θ) ≡ 3^"1{g(κ,φ)} ≡ Jg(κ,φ)e^iκrcos(θ-^,'^,)κdκdθ (6)

Equation (4) can be transformed by appropriate mathematical manipulations into the following form which gives S^(α)(r,θ) in terms of the measured acoustic data mj (cθα).

S^(α)(r,θ) = 3 -1

5 51122 (7)

-3 ∑∑JJp 'ς s{[J (k.r) } l p,q=-512

D . Creating the Images

The full process of producing the S^(α)(r,θ) 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_jk ^cιdent (ω_α) (j = l,..,N;j ≠ k) as described in Section B at each of 10 discrete frequencies ω_α (α = 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 ω_α. These calibration factors are stored as a separate 1024 x 1024 calibration matrix for each of the 10 discrete frequencies ω_α (α = 1,2,..., 10). These calibration matrices are stored in the memory of computer 72.

Step (2). Data representing mj (ω_α) 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 ωα, πjk ( α) is an 1024 x 1024 complex matrix of phase and amplitude data.

Step (3). For each discrete frequency ω , computer 72 multiplies the mjk (ω_α) data by the calibration matrix determined in Step (1) of this Section to provide calibrated mjk (Cύα) data.

Step (4). Using the calibrated mjk (ω_α) data, computer 72 calculates the quantity

for each frequency ω_α. This manipulation constitutes a fast Fourier transform of the calibrated complex matrix mjk (ω_α). The result of this transformation is that the calibrated mjk (ωα) data for each frequency ω_α is now in the form of a 1024 x 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_p(k₀r₀) and H_q(k₀r₀) 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 ω_α and each mode combination (p,q) to calculate the following quantity:

Step (7). We now calculate a weighting function J_p(k₀r)J_q(k₀r)e^"lθ^^p+q^ at each point r and θ which is a 1024 x 1024 dimensional matrix. J_p and J_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:

This step is done at each point (r,θ) in the medium 12. The result of this step is S^_v(r,θ). This is an image of S^(α)(r,θ) but it is convolved with [J₀(k₀r)]² which presents a smoothed image of

S⁽° (r,θ). For a more accurate image we need to deconvolve the quantity [J₀(kor)]² from S^_v(r,θ).

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₀(k₀r)]² is -p-. Thus we deconvolve π κ (4k₀ - κ)^1/2

[J_o(k_or)]² from sjjj_v(r,θ) to calculate an image S⁽«⁾(r,θ) as follows:

S^(α)(r,θ) = 3-^,{πκ(4k₀ - κ)^1/23{Sl _v(r,θ)}}. (11)

Using equation (11) we can produce an image S^(α)(r,θ) at each frequency ω_α of medium 12: To show the images on a digital computer monitor, computer 72 is programmed to convert from polar coordinates (r,θ) to cartesian coordinates (x,y) to produce S^(α)(x,y). Each of the S^(α)(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^(α)(r,θ) at discrete frequencies ω_α 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^(α)(r,θ) acquired at different frequencies ω_α ideally would present the same images of the medium 12, independent of frequency ω_α. However, the images can be complicated by frequency dependent artifacts resulting from reflections of acoustic energy from the faces of the receiver transducers lO.k and from multiple scattering events within the medium 12. We have found that combining images of S^{( )}(r,θ) reconstructed at different frequencies ω_α reduces the image artifacts and improves the image quality. The following discussion centers on seven preferred ways of combining single frequency images S^(α)(r,θ) reconstructed at different frequencies (Dα- For the purposes of combining images S^(α)(r,θ) formed at different frequencies ω_α, we have found it useful to define a term S^(r,θ) which we calculate as:

S^(α)(r,θ) = 3-¹{κ-^I3{s^(α)(r,θ)}} . (12)

In equation (12), S^(α)(r,θ) is calculated by first calculating an FFT of S^(α>(r,θ), then dividing by K and calculating the inverse FFT of the result.

1) Sum of the complex images S^(α)(r,θ).

From the acquired data mj_k (ω_α), we reconstruct separate images of S^{( )}(r,θ) at each frequency ω_α (α = 1,2,..., 10) as given in Steps (1) - (9) of Section D. We then sum these complex images to produce:

ιo , . S_totaiM) = ∑s^(α)(r,θ). (13) α=l

We convert to cartesian coordinates (x,y) to produce an image S_tota_l(r,θ) which has less artifact and more clarity than each single frequency image S^(α)(r,θ).

2) Sum of magnitudes of S^(α)(r,θ) Images

From the acquired data mj_k (ω_α), we reconstruct separate images of S^(α)(r,θ) at each frequency ω (α = 1,2,...,10) as given in Steps (1) - (9) of Section D. We then calculate the magnitudes of these complex images as IS⁽«⁾(r,θ)l = [S⁽«⁾(r,θ)S⁽«^)*(r,θ)]^{1 2} where S⁽«^)*(r,θ) is the complex conjugate of S^(α)(r,θ). We then sum the magnitudes of these images to produce:

¹⁰ , \ IS_total(r,θ)l= ∑IS^(α)(r,θ)l. (14) α=l

Again we convert to cartesian coordinates to produce an image IS_total(r,θ)l which has less artifacts and more clarity than each image IS^(α)(r,θ)l. 3) Sound Speed and Attenuation Maps Calculated from S^(α)(r,θ) Images

Sound speed maps and attentuation maps of medium 12 can be calculated from S^(α)(r,θ) as follows:

Step (1). From the acquired data m_jk (ω_α), we reconstruct separate images of S⁽ \r,θ) at each frequency ω_α (α = 1,2,..., 10) as given in Steps (1) - (9) of Section D. We then calculate the mapsS^(r,θ) at each frequency ω_α as given in equation (12).

Step (2). Computer 72 then calculates the following function

P(r,θ,τ) =

For each point (r,θ) in the medium 12, P(r,θ,τ) represents a synthesized acoustic pulse P(τ) focused at the point (r,θ).

Step (3). For each point (r,θ) in the medium 12, P(τ) has a maximum value P_max at a value of τ which we call τ_max. The computer 72 is programmed to find, for each point (r,θ), P_max(r,θ) and to determine the value of τ_max(r,θ)-

Step (4). The computer 72 then creates a sound speed map c(r,θ) of medium 12 by the following calculation c(r,θ) = 3-¹{κ3{2r₀ / τ_max(r,θ)}} + c₀. ( 16)

In equation (16), r₀ is the radius of the transducer ring and c₀ is the sound speed of the coupling fluid 11.

Step (5). The computer then creates an attenuation map N(r,θ) of medium 12 by the following calculation

In equation (17), μo is the attenuation coefficient of the coupling fluid 11. These values of c(r,θ) and N(r,θ) are 1024 x 1024 matrixes in r and θ. They are converted to cartesian coordinates (x,y) and displayed on a computer screen with the values of c and N depicted as gray scale values.

4) Refined Sound Speed and Attenuation Maps Calculated from S^(α)(r,θ)

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

Step (1). From the acquired data m_jk (ω_α), we reconstruct separate images of S^(α)(r,θ) at each frequency ω_α (α = 1,2,..., 10) as given in Steps (1) - (9) of Section D. We then calculate the maps S^^α^(r,θ) at each frequency ω_αas given in equation (12).

Step(2). We then calculate the maps S :(_Δα_p)_Δq(r,θ) at each frequency ω_α (α = 1,2,..., 10) as such:

The maps S_Δ_ _Δq(r,θ) represent coupling of acoustic energy from azimuthal modes p and q to modes p + Δp and q + !q.

Step (3). Computer 72 then calculates the following function:

P_corr(r,θ,τ) = (19)

We have had excellent results using sums of !p = -6 to +6, !q = 0; !p = 0, !q = -6 to +6; and !p = -6 to +6, !q = -6 to +6. Other combinations of !p and !q could be used. For each point (r,θ) in the medium 12, P_Corr(^r,θ,τ) represents a synthesized acoustic pulse P_corr(τ) focused at the point (r,θ). P_corr(r,θ,τ) is a refined version of the function P(r,θ,τ) given in equation (15).

Step (4). For each point (r,θ) in the medium 12, P_COrr(τ) has a maximum value P_max at a value of τ which we call τ_max- The computer 72 is programmed to find, for each point (r,θ), P_max(r,θ) and to determine the value of τ_max(^r,θ) as was done in Step (3) of Section E.3 except we use Pcorr (r,θ,τ) instead of P(r,θ,τ)).

Step (5). Using the refined values of τ_max(r,θ) and P_maχ(r,θ), we now calculate the sound speed c(r,θ)) and attenuation N(r,θ) maps as we did in equations (16) and (17).

5) Sum of Phase Aberration Corrected S^(α)(r,θ) Images

The complex images S^(α)(r,θ) 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,θ) 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^(α)(r,θ) and the combination of these images to reduce image artifacts.

Step (1). From the acquired data mj_k (ω_α), we reconstruct separate images of S^(α)(r,θ) at each frequency ω_α (α = 1,2,...,10) as given in steps (1) - (9) of Section D. We then calculate the maps s' ' (r, θ) at each frequency 0 as given equation (12).

Step (2). We calculate the map τm_ax(r,θ) as described by the procedures given in Section E.3 or section E.4.

Step (3). Computer 72 calculates the map S_total (r,θ) as such:

'total (r,θ)

Step (4). The final image S_tota_l(r,θ) is determined by the following calculation: S_total(r,θ) = 3-'{κ3{s_total(r,θ)}} (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 _ota_l(^r,θ) calculated in Section E.5 can be combined with the map τ_maχ(r,θ) to provide a composite map S_comp(r,θ) which is a more accurate representation of the medium 12. The S_COmp(r,θ) image contains both the high spatial frequency components and low spatial frequency components of the medium 12.

Step (1). The computer 72 calculates τ_max(r,θ) as described by the procedures given in Section E.3 or E.4 and S_total (r,θ) given by equation (20).

Step (2). Computer 72 then calculates:

S_c__mP.(r.θ) = + τ_max(r,θ) (22)

In equation (22), In denotes the natural logarithm and Δω is the spacing between frequencies ω_α (Δω = 62.5 kHz in the preferred embodiment).

Step (3). Computer 72 then calculates S_comp as follows:

S_comp ) = 3-¹{ 3{s_comp(r,θ)}}. (23)

Computer 72 then converts S_COmp(^r,θ) to cartesian coordinates (x,y) and displays on computer monitor.

7) Further Refined Sound Speed Map Calculated From S^(α)(r,θ) This section discusses a method to further refine the accuracy of the sound speed maps c(r,θ) calculated in Sections E.3 and E.4. In Sections E.3 and E.4 we calculate at each point (r,θ) the transit time τ_max (r,θ) as the value τ in which the pulse P(τ) has a maximum value P_maχ and derive sound speed maps from this information. In this Section, we calculate the transit time τ_av (r,θ) as the averaged arrival time of the power contained in the pulse P(τ). This map τ_av (r,θ) is then used to calculate refined maps of the sound speed c(r,θ) of the medium 12.

Step (1). Computer 72 computes the following quantity:

where S_Δ" _Δq(r,θ) is calculated as in equation (18) and * denotes the complex conjugate. The preferred embodiment includes sums of 1) !p = -6 to +6; !q = 0; 2) !p = 0, !q = - 6 + 6; or 3) !p = - 6 + 6, !q = - 6 to + 6. Other combinations of Δp and Δq could be used.

Step (2). The computer 72 calculates an average transit time τ_av (r,θ) image of medium 12 by the following calculation

τ_av(r,θ) = ^lnP_av(r,θ) (25)

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

c(r,θ) = 3^_1{κ3{2r₀ / τ_av(r,θ)}} + c₀. (26)

Step (4). The computer 72 then creates an attenuation image μ(r,θ) of medium 12 by the following calculation

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.l), 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 and 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,θ) { heated } = c(r,θ) + ΔT(r,θ) [d/dT c(r,θ)] (26)

where ΔT(r,θ) is the temperature distribution resulting from local heating of the tissue, c(r,θ) is the sound speed map of the tissue at the ambient body temperature, and d/dT c(r,θ) is the differential change in sound speed with respect to temperature. For example, in human liver tissue at 37 °C, c(r,θ) = 1596 m/sec and d/dT c(r,θ) = 0.96 m/(sec -°C). A ΔT = 5 °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,θ) -ⁱJΔTCr.θ) [d/dT c(r,θ)] = 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

We Claim:

1. An acoustic imaging device for producing a tomographic image along an image slice through a medium defining an average sound speed, said device comprising:

A) a large number of acoustic transducers located on a circle said large number being greater than N = 2πD/λ where D is the diameter of the circle and λ = 0.3 mm each of said large number of transducers being connected electronically so as to broadcast into said medium an acoustic signal when excited with a periodic signal in the range of 100 kHz to 1.5 MHz and also to produce an electrical signal when excited with an acoustic signal transmitted through said medium,

B) an electronic signal generation means for generating periodic electronic signals at at least one discrete frequency in the frequency range of 100 kHz to 1.5 MHz said signals being continuous for a time interval of at least D/c where c is the average sound speed of the medium,

C) a first multiplexing means for applying said periodic electronic signals to a each of said large number of acoustic transducers in order to cause a plurality of said acoustic transducers, one at a time, to broadcast acoustic energy into said medium at a discrete frequency in the frequency range of 100 kHz to 1.5 MHz,

D) a recording means for recording information carried by electronic transducer signals produced by a plurality of said large number of transducers excited by acoustic energy broadcast transmitted through said medium,

E) a second multiplexing means for connecting said large number of transducers to said recording means to permit said recording means to record said information carried by said electronic transducer signals, and

F) a computer means comprising an algorithm for calculating a tomographic image of said medium along said image slice utilizing said information recorded by said recording means.

2. An acoustic imaging device as in Claim 1 wherein said at least one identifiable discrete frequency is a plurality of identifiable discrete frequencies. 3. A device as in Claim 1 wherein said transducers are equally spaced on said circle.

4. A device as in Claim 1 wherein said large number is at least 256.

6. A device as in Claim 1 wherein said acoustic recording means comprises a plurality of analog-to- digital converter units each having a buffer memory unit for temporary storage of digital acoustic information.

7. A device as in Claim 6 wherein each of said buffer memory units comprise at least one million Bytes of storage capacity.

8. A process for producing a tomographic image of a slice through a medium defining an average speed of sound comprising the steps of:

A) positioning a large number of transducers on a circle surrounding said medium, said large number being greater than N = 2πD/λ where D is the diameter of the circle and λ = 0.3 mm, said transducers having the capability of broadcasting an acoustic signal into said medium upon being excited with an electrical signal and producing an electrical signal upon being excited with an acoustic signal from said medium,

B) exciting a plurality of said large number of transducers, one at a time, with an electronic signal defining at least one discrete frequency in order to cause each of said transducers to broadcast, one at a time, into said medium an acoustic signal containing said at least one discrete frequency, each of said acoustic signals being broadcast for a time interval of at least D/c where c is equal to the average speed of sound in said medium said signal to reach an equilibrium in said medium and excite other transducers on said circle to cause them to produce electrical transducer signals carrying information relating to acoustic properties of said medium,

C) record information carried by said electrical transducer signals from said large number of transducers,

D) using said recorded information calculate with a computer said tomographic image of said image slice of said medium. 9. A process as in Claim 8 wherein said at least one identifiable discrete frequency is a plurality of identifiable discrete frequencies.

10. A process as in Claim 7 comprising the additional step of producing a plurality of said tomographic images of said image slices of said medium.

1 1. A process as in Claim 10 comprising the additional step of producing a plurality of said tomographic images of said image slices of said medium said image slices being parallel to each other and equally spaced along said medium.

12. A process as in Claim 11 comprising the additional step of producing a plurality of said tomographic images of said image slices of said medium said image slices being parallel to each other and equally spaced along said medium and said image slices overlapping each other.

13. A process as in Claim 12 comprising the additional step of forming a three-dimensional image from said tomographic images of said image slices.

14. A process as in Claim 8 wherein said medium includes human tissue and comprising the following additional steps for monitoring the placement of a biopsy needle in said tissue:

A) calculate a sufficient number of tomographic images of said medium in order to identify a biopsy location,

B) insert a biopsy needle to said biopsy location,

D) calculate at least one tomographic image of said medium with said needle at said biopsy location,

E) verify that said biopsy needle is at said location, and

F) remove biopsy sample with said biopsy needle.

15. A process as in Claim 8 comprising the following additional steps for measuring temperature distributions of said medium comprising:

A) locally heat at least a portion of said image slice of said medium, B) calculate second said tomographic image of said image slice of said medium, and

C) subtract first said tomographic image from said second tomographic image to produce a final image representing the effect of heating the medium.

16. A process as in Claim 8 comprising the additional steps of measuring temperature distributions of said medium comprising:

A) locally heat at least a portion of said image slice of said medium,

B) calculate second said tomographic image of said image slice of said medium, and

C) divide first said tomographic image from said second tomographic image to produce a final image representing the effect of heating the medium.

17. A process as in Claim 8 wherein said medium comprises at least some portion of the human body.

18. A process as in Claim 8 wherein said medium comprises the human female breast.

19. A process as in Claim 8 wherein said medium comprises the human abdomen.

20. 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,

Q 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 beings 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 weighing 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) deconvolve said convolved image of said medium to produce said image of said medium. 21. An acoustic image device as in Claim 20 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.

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

23. An acoustic image device as in Claim 21 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.

24. An acoustic imaging system as in Claim 21 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 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.

25. An acoustic imaging system as in Claim 24 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.

26. An acoustic imaging system as in Claim 25 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.

27. An acoustic imaging device as in Claim 26 wherein said function is the diameter of said circle divided by said time-of-arrival image.

28. An acoustic imaging system as in Claim 24 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 medium and then to said detection locations said maximum pulse power data forming a peak pulse power image.

29. An acoustic imaging system as in Claim 28 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.

30. An acoustic imaging device as in Claim 29 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

31. An acoustic imaging system as in Claim 21 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

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

32. An acoustic imaging system as in Claim 31 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.

33. An acoustic imaging system as in Claim 32 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.

34. An acoustic imaging device as in Claim 33 wherein said function is the diameter of said circle divided by said refined time-of-arrival image. 35. An acoustic imaging system as in Claim 31 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.

36. An acoustic imaging device as in Claim 35 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.

37. An acoustic imaging device as in Claim 36 wherein said function is the logarithm of said refined maximum pulse power image divided by the diameter of said circle.

38. An acoustic imaging system as in Claim 25 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,

Q 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.

39. An acoustic imaging device as in Claim 38 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.

40. An acoustic imaging system as in Claim 32 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,

Q sum said plurality of refined phase aberration corrected convolve 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 device as in Claim 40 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.

42. An acoustic imaging system as in Claim 38 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.

43. An acoustic imaging system as in Claim 40 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. 44. An acoustic imaging device as in Claim 20 wherein said at least discrete frequency is a plurality of discrete frequencies.

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

46. An acoustic imaging device as in Claim 45 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.

47. An acoustic imaging device as in Claim 45 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.

48. An acoustic imaging system as in Claim 45 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.

49. An acoustic imaging system as in Claim 48 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.

50. An acoustic imaging system as in Claim 49 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.

51. An acoustic imaging device as in Claim 50 wherein said function is the diameter of said circle divided by said time-of-arrival image. 52. An acoustic imaging system as in Claim 48 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.

53. An acoustic imaging system as in Claim 52 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.

54. An acoustic imaging device as in Claim 53 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

55. An acoustic imaging system as in Claim 45 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 a said large number of locations in said medium,

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

56. An acoustic imaging system as in Claim 55 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.

57. An acoustic imaging system as in Claim 56 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. 58. An acoustic imaging device as in Claim 57 wherein said function is the diameter of said circle divided by said time-of-arrival image.

59. An acoustic imaging system as in Claim 55 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.

60. An acoustic imaging system as in Claim 59 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.

61. An acoustic imaging device as in Claim 60 wherein said function is the logarithm of said maximum pulse power image divided by the diameter of said circle.

62. An acoustic imaging system as in Claim 55 wherein said computing means is programmed to accomplish the following additional steps:

63. An acoustic imaging device as in Claim 62 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. 64. An acoustic imaging system as in Claim 63 wherein said computing means is programmed to accomplish the following additional steps:

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

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.

65. An acoustic imaging device as in claim 62 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.

66. An acoustic imaging system as in Claim 62 wherein said computing means is programmed to accomplish the following additional steps:

A) sum of 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.

67. An acoustic imaging system as in Claim 64 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 deconvolve said refined convolved composite image to produce a refined composite image.